© The Authors 2024

1. Introduction

The James Webb Space Telescope (JWST) plays a pivotal role in advancing our understanding of the high-redshift Universe. Its observation of galaxies at z ≳ 9 in unprecedented numbers offers a unique opportunity to unravel the properties of galaxies emerging in the first few hundred million years, thereby enhancing our knowledge and constraints on galaxy evolution and reionisation. Early release observations have given us glimpses into the properties of these galaxies, painting an initial picture of compact, metal-poor, young, star-forming galaxies with stellar populations emitting harder ionising radiation (e.g. Bradley et al. 2023; Bunker et al. 2023; Curtis-Lake et al. 2023; Heintz et al. 2023b; Schaerer et al. 2022; Vanzella et al. 2023). However, preliminary inferences based on photometric data suggest a higher abundance of UV bright galaxies than predicted by most theoretical models (Labbé et al. 2023; Adams et al. 2024; Atek et al. 2023; Austin et al. 2023; Boylan-Kolchin 2023) and that the corresponding UV luminosity functions evolves only mildly at z > 9 (Adams et al. 2024; Castellano et al. 2023; Donnan et al. 2023; Finkelstein et al. 2023; Harikane et al. 2023a, 2024; Naidu et al. 2022).

These discoveries have sparked theoretical investigations into potential selection biases and modifications in physical processes that could account for the observed abundances. Some studies have proposed that the observed galaxies might be biased samples, exclusively tracing the densest regions containing the most massive galaxies (McCaffrey et al. 2023) or detecting galaxies undergoing periods of intense star formation in their bursty star formation histories (Mason et al. 2023; Sun et al. 2023). While different simulations and observations agree on the presence of bursty star formation in early galaxies (e.g. Legrand et al. 2022; Gelli et al. 2023; Ciesla et al. 2024, and references therein), it remains unclear whether bursty star formation alone can explain the high abundance of bright z > 9 galaxies (e.g. Sun et al. 2023; Pallottini & Ferrara 2023). Conversely, others have explored the altered physical conditions in these very high-redshift galaxies. One possibility to reproduce the observed UV luminosity function at z > 9 could be the ejection of dust from the star formation site through radiatively driven outflows during the initial phases of galaxy formation, such that the reduced dust attenuation compensates for the increasing shortage of bright galaxies predicted in standard theoretical models at higher redshifts (Ferrara et al. 2023; Fiore et al. 2023; Mauerhofer & Dayal 2023; Yung et al. 2024; Ziparo et al. 2023). Another explanation for these UV-bright objects could be that early galaxies exhibited higher star formation efficiencies. For example, Dekel et al. (2023) estimated that if the gas in the most massive early galaxies is sufficiently dense and metal-poor, the free-fall time becomes shorter than the time required for low-metallicity stars to develop winds and SNe, resulting in feedback-free starbursts. Similarly, weaker stellar winds and fewer SNe, typical for low-metallicity stars, could weaken and delay the onset of mechanical feedback to around ∼10 Myr (Jecmen & Oey 2023; Yung et al. 2024). An alternative explanation could involve the UV luminosity produced by massive black holes (≳ 10⁸ M_⊙) that reside in UV-bright galaxies and accrete at or slightly above the Eddington rate quasars (Pacucci et al. 2022).

Alternatively, early galaxies could have a lower mass-to-light ratio due to a higher abundance of massive stars, implying a more top-heavy IMF as suggested in Haslbauer et al. (2022), Trinca et al. (2024), Woodrum et al. (2023), Harikane et al. (2023a, 2024), Finkelstein et al. (2023), Inayoshi et al. (2022), Pacucci et al. (2022). A more top-heavy IMF would also lead to higher ionising emissivities that would not only increase the strength of emission lines but also impact the morphology of the reionisation process of the intergalactic medium (IGM).

Indeed, simulations of the first (metal-free) stars predominantly result in stellar masses of ≳60 M_⊙ (Abel et al. 2002; Bromm et al. 2002; Yoshida et al. 2006; Fukushima et al. 2020), hinting that the fraction of massive stars may rise with increasing gas temperature of star-forming clouds. The gas temperature not only determines whether a region within a cloud will experience gravitational collapse but also influences the presence of substructure, which, in turn, regulates the amount of material available for accretion during the collapse and subsequent formation of the star (cf. results in Schneider & Omukai 2010; Chon et al. 2022; Sharda & Krumholz 2022). This dependence of the IMF on the gas temperature aligns with the results obtained from hydrodynamical simulations, showing that the IMFs of star-forming clouds become more top-heavy with a decreasing metallicity of the gas and increase in the background radiation intensity, such as the cosmic microwave background (CMB) (Chon et al. 2022). Indeed, the spectra of two Lyman-α emitters at z ≥ 5.9 (Cameron et al. 2023), along with the putative detections of Population III stellar populations at z = 6.6 and z = 10.6 characterised by the absence of metal lines and either strong He II or Lyman-α, Hγ, Hβ, and Hα lines (Maiolino et al. 2024; Vanzella et al. 2023), all observed with JWST, collectively argue for top-heavy IMFs. Furthermore, applying gas temperature-dependent IMF models to galaxy observations, where the temperature is treated as an additional parameter in the photometric template fitting, reveals that most star-forming galaxies at a fixed redshift have similar, but top-heavier IMFs than the Milky Way, with inferred gas temperatures increasing towards higher redshifts (Sneppen et al. 2022; Steinhardt et al. 2022b,a). In addition, observations of local star-forming regions suggest that the IMF varies across correlated star formation events. For instance, studies of globular clusters, ultra-compact dwarf galaxies and young massive clusters indicate that the IMF becomes top-heavy in low metallicity and dense gas environments (e.g. Dabringhausen et al. 2009, 2012; Marks et al. 2012; Zonoozi et al. 2016; Haghi et al. 2017; Kalari et al. 2018; Schneider et al. 2018; Dib 2023). In contrast, it becomes more bottom-heavy in metal-rich environments (e.g. Marks et al. 2012; Chabrier et al. 2014), such as in the centres of nearby elliptical galaxies (e.g. van Dokkum & Conroy 2010; Conroy et al. 2017). Averaging the IMF over all star-forming regions in a galaxy, this galaxy-wide IMF (gwIMF; Kroupa & Weidner 2003; Weidner et al. 2013; Yan et al. 2017; Jeřábková et al. 2018) is often found to be top-heavy in galaxies with high star formation rates and top-light in galaxies with low star formation rates (e.g. Meurer et al. 2009; Watts et al. 2018; Zhang et al. 2018; Gunawardhana et al. 2011; Fontanot et al. 2017, 2018).

The above-quoted studies highlight that the IMF in early galaxies is very likely to differ from the canonical Milky Way and present-day IMF that is used in nearly all state-of-the-art large-scale reionisation radiation hydrodynamics and semi-numerical simulations following the evolution of galaxies explicitly (Ocvirk et al. 2020; Lewis et al. 2022; Kannan et al. 2022; Gnedin 2014; Mutch et al. 2016; Hutter et al. 2021). To address this discrepancy, some simulations and semi-analytic models of galaxy evolution have introduced a distinction between the top-heavy IMF associated with the first stars (Population III) and the present-day canonical IMF for metal-enriched Population II stars (e.g. Maio et al. 2010; Norman et al. 2018; Visbal et al. 2020). However, these do not consider the IMF of Population II stars to evolve or depend on the galaxies’ properties. Only at lower redshifts, a few hydrodynamics simulations and semi-analytic galaxy evolution models have incorporated and explored the integrated galaxy-wide IMF model (Weidner & Kroupa 2005; Weidner et al. 2013) (assuming that the most massive star in a cluster is linked to its cluster mass) to explore the impact of an empirically informed varying IMF on the evolution of galaxy properties (Ploeckinger et al. 2014; Fontanot et al. 2017). More recently and at higher redshifts, only Trinca et al. (2024) have explored the effect of an IMF varying with stellar metallicity and redshift on the galaxy UV luminosity functions, finding the resulting UV LFs to better fit the observations at z > 9. However, these authors considered this evolving IMF only when deriving the galaxies’ UV luminosities in their semi-analytic galaxy evolution model, not when evaluating stellar feedback and metal yields for determining galaxy properties or the ionising photon production for reionisation.

Omitting the IMF’s dependency on stellar feedback can lead to severely erroneous conclusions. Firstly, a more top-heavy IMF will not only reduce the mass-to-light ratio, but also enhance the fraction of stars exploding as SNe, thereby reducing subsequent star formation more immediately while increasing the metal enrichment of the interstellar gas per exploding SN. Given the theoretical and observational hints towards an evolving IMF, we introduce the first model that follows the mutual evolution of galaxies and reionisation, while assuming an IMF that evolves in each galaxy according to the metallicity of its star-forming gas and redshift. For this purpose, we employed a parameterisation of the evolving IMF that follows the results from the spherical turbulent gas cloud simulations in the presence of the CMB at various gas metallicities and redshifts (Chon et al. 2022). We incorporated this IMF parameterisation into our ASTRAEUS framework, a semi-numerical model that tracks the interdependent evolution of galaxies and reionisation (Hutter et al. 2021; Ucci et al. 2023; Hutter et al. 2023) by incorporating IMF-dependence into all relevant physical processes. With this updated model, we investigate the following question during the Epoch of Reionisation (EoR), which spans between z ≃ 5 − 15 and during which most of the intergalactic medium (IGM) was ionised at z ≲ 7.5 (Planck Collaboration VI 2020): How does an evolving IMF affect the properties of early galaxies and reionisation, and can it explain the evolution of the observed UV LFs during the EoR?

This paper is organised as follows. In Sect. 2, we briefly describe the ASTRAEUS model and outline how we parameterise and integrate an evolving IMF into ASTRAEUS. Section 3 describes how we tune the free ASTRAEUS model parameters to reproduce the observed UV luminosity functions and reionisation history constraints for both the constant Salpeter IMF and our new evolving IMF scenarios. We then discuss how an evolving IMF changes the mass-to-UV light ratio, stellar-to-halo mass ratio, star formation sequence, stellar mass to metallicity and stellar mass to dust mass relations, as well as the UV luminosity to metallicity and UV luminosity to dust mass relations, compared to a constant Salpeter IMF (Sect. 4). We present our conclusions in Sect. 5. In this paper, we assume the AB magnitude system (Oke & Gunn 1983) and a ΛCDM universe with the following Planck Collaboration Int. XLVI (2016) cosmological parameters: Ω_Λ = 0.692885, Ω_m = 0.307115, Ω_b = 0.048206, H₀ = 100h = 67.77 km s⁻¹ Mpc⁻¹, n_s = 0.96, and σ₈ = 0.8228.

2. The model

The ASTRAEUS framework couples an enhanced version of the semi-analytic galaxy evolution model DELPHI (Dayal et al. 2014, 2022) with the semi-numerical reionisation scheme CIFOG (Hutter 2018). The resulting model runs on the outputs of a dark matter (DM) only N-body simulation. In this section, we briefly revisit the physical processes tracked in ASTRAEUS and described in Hutter et al. (2021), Ucci et al. (2023) and Hutter et al. (2023) in detail. Here, we also describe our implementation of an evolving IMF.

2.1. The N-body simulation

We used the VSMDPL (very small multidark planck) DM-only N-body simulation, which is part of the Multidark simulation project¹ and has been run with the GADGET-2 TREE+PM code (Springel 2005). The simulation encompasses a cubic volume with a side length of 160 h⁻¹ comoving Mpc (cMpc) and tracks the trajectories of 3840³ DM particles. Each DM particle carries a mass of 6 × 10⁶ h⁻¹ M_⊙. For 150 snapshots spanning from z = 25 to z = 0, halos and subhalos down to 20 particles or a minimum mass of 1.24 × 10⁸ h⁻¹ M_⊙ have been identified with the phase space ROCKSTAR halo finder (Behroozi et al. 2013a). Since ASTRAEUS includes time-synchronised processes like reionisation, we have used the pipeline internal CUTNRESORT scheme to re-sort the vertical merger trees generated by CONSISTENT TREES (Behroozi et al. 2013b) to local horizontal merger trees for all galaxies at z = 4.5 (for details see Appendix A in Hutter et al. 2021). For the first 74 snapshots spanning from z = 25 to z = 4.5, we generate the DM density fields required as input for the ASTRAEUS pipeline by mapping the DM particles onto 2048³ grids and subsequently resampling to these 512³ grids.

2.2. Galaxy evolution

ASTRAEUS² follows the key physical processes of early galaxy formation and reionisation. By post-processing the DM merger trees and density fields from the VSMDPL simulation, it tracks for each galaxy and at each time step (i.e. snapshot of the N-body simulation) the amount of gas accreted, the gas and stellar mass merged, the formation of stars and their feedback through SNe and metal enrichment, as well as the large-scale reionisation process and its feedback on the gas content in galaxies. The modelling of these processes is described below.

2.2.1. Gas and stars

Each galaxy that starts forming stars in a halo with mass M_h, is assumed to have a gas mass of $M_g^i(z)=f_g({\mathrm{\Omega }}_b/{\mathrm{\Omega }}_m)M_h(z)$. Here, f_g describes the gas fraction that is not evaporated by reionisation, assuming values of f_g < 1 and f_g = 1 as the galaxy forms in an ionised and neutral region, respectively³. In subsequent time steps, the gas mass of the galaxy includes the gas inherited from its progenitor, $M_g^{\text{mer}}(z)$ and gained by smooth accretion, $M_g^{\text{acc}}(z)$, but never exceeds the limit given by reionisation feedback (Gnedin 2000; Sobacchi & Mesinger 2013):

$$\begin{array}{c}\hfill M_g^i(z)=min(M_g^{\mathrm{mer}}(z)+M_g^{\mathrm{acc}}(z),f_g({\mathrm{\Omega }}_b/{\mathrm{\Omega }}_m)M_h),\end{array}$$(1)

with

$$\begin{array}{cc}& M_g^{\mathrm{acc}}(z)=\frac{{\mathrm{\Omega }}_b}{{\mathrm{\Omega }}_m}[M_h(z)-{\displaystyle \sum _{p=1}^{{\mathrm{N}}_{\mathrm{p}}}}{M}_{h,p}(z+\mathrm{\Delta }z)],\hfill \end{array}$$(2)

$$\begin{array}{cc}& M_g^{\mathrm{mer}}(z)={\displaystyle \sum _{p=1}^{{\mathrm{N}}_{\mathrm{p}}}}{M}_{g,p}(z+\mathrm{\Delta }z).\hfill \end{array}$$(3)

Here, N_p is the number of progenitors of a galaxy, while M_h, p and M_g, p are each of its progenitor’s halo and gas masses at z + Δz. At each time step, we assume that a fraction of the merged and accreted gas mass $M_g^i$ forms stars over the time step’s length, Δt, amounting to a mass of newly formed stars of $M_{\star }^{\text{new}}(z)=f_{\star }^{\text{eff}}{M}_g^i(z)$. $f_{\star }^{\text{eff}}$ depicts the fraction of gas that forms stars in Δt, namely, the star formation efficiency. We note that its value depends on the gravitational potential of the galaxy: In massive galaxies the fraction of gas converted into stars reaches the maximum value, given by $f_{\star }\frac{{\tau }_{\mathrm{dyn}}(z=9)}{{\tau }_{\mathrm{dyn}}(z)}\frac{20\mathrm{Myrs}}{\mathrm{\Delta }t}$, where ${\tau }_{\mathrm{dyn}}(z)=\sqrt{r_{\mathrm{vir}}^3/(G{M}_h)}\propto {(1+z)}^{-3/2}$ is the dynamic time of the system at redshift z⁴. However, in lower mass galaxies, SN and radiative feedback from reionisation limit this fraction further. Our model assumes that a galaxy can form only as many stars as are necessary to eject all gas from the galaxy due to SN explosions, with the corresponding fraction given as:

$$\begin{array}{c}\hfill f_{\star }^{\mathrm{ej}}(z)=\frac{v_{\mathrm{c}}^2}{v_{\mathrm{c}}^2+f_w^{\mathrm{eff}}(z)E_{51}{\nu }_z}[1-\frac{f_w^{\mathrm{eff}}(z)E_{51}{\sum }_j{\nu }_j{M}_{\star ,j}^{\mathrm{new}}(z_j)}{M_{\mathrm{g}}^i(z)v_{\mathrm{c}}^2}].\end{array}$$(4)

Here, v_c is the rotational velocity of the halo, E₅₁ the energy released by a Type II supernova (SNII), ν_z the IMF-dependent fraction of stellar mass forming and exploding in the current time step, and

$$\begin{array}{c}\hfill f_w^{\mathrm{eff}}(z)=\frac{f_w}{1+(\frac{20\mathrm{Myr}}{\mathrm{\Delta }t}-1)\frac{M_g^{\mathrm{mer}}(z)}{M_g^i(z)}}\end{array}$$(5)

is the fraction of SN energy injected into the winds driving gas outflows. The right term in brackets arises from our delayed SN feedback scheme, where at each time step the effective star formation efficiency accounts also for the SN energy released from stars formed in previous time steps, following the mass-dependent stellar lifetimes in Padovani & Matteucci (1993). Thus, $M_{\star ,j}^{\text{new}}(z_j)$ is the stellar mass formed during previous time step j, and ν_j the fraction of stellar mass formed in previous time step j that explodes in the current time step given the IMF of the respective progenitor. Thus, the fraction of gas converted into stars, effectively the star formation efficiency, is given by $f_{\star }^{\mathrm{eff}}=min(f_{\star }\frac{{\tau }_{\mathrm{dyn}}(z=9)}{{\tau }_{\mathrm{dyn}}(z)}\frac{20\mathrm{Myrs}}{\mathrm{\Delta }t},f_{\star }^{\mathrm{ej}})$ at each time step. Both, f_⋆ and $f_w^{\text{eff}}$ are free model parameters. The total stellar mass at redshift, z, is then:

$$\begin{array}{c}\hfill M_{\star }(z)={\displaystyle \sum _{p=1}^{N_p}{M}_{\star ,p}(z+\mathrm{\Delta }z)+M_{\star }^{\mathrm{new}}(z)-M_{\star }^{\mathrm{destr}}(z),}\end{array}$$(6)

where $M_{\star }^{\text{destr}}(z)$ is the stellar mass that is returned into gas and metals in SNe and AGB stars.

2.2.2. Metals and dust

ASTRAEUS tracks the metal enrichment by stellar winds, SN Type II (SNII) and SN Type Ia (SNIa) explosions and AGB stars (for a detailed description see Ucci et al. 2023). At each time step, we assume the smoothly accreted gas to have the average metallicity of the intergalactic medium (IGM), Z_IGM. The quantity of newly formed metals depends on the number of massive stars exploding as SNe during the current time step. We use the stellar lifetimes from Padovani & Matteucci (1993), metal enrichment rates from stellar winds, SNII and SNIa as described in Yates et al. (2013) with the SNIa progenitor fraction and delay-time distribution as defined in Arrigoni et al. (2010) and Maoz et al. (2012), and the most recent metal yields from Kobayashi et al. (2020). Furthermore, we assume the galaxy’s gas and metals to be perfectly mixed. Thus, the amount of ejected metals is directly proportional to the ejected gas mass and the metallicity of the gas in the galaxy and contributes to Z_IGM.

ASTRAEUS also follows the formation, growth, destruction, astration, and ejection of dust in each galaxy, where our model assumes the dust mass reservoir to be part of the metal mass reservoir. Specifically, our model assumes that dust is produced by SNII and AGB stars through the condensation of metals in stellar ejecta, and its grains grow through the accretion of heavy elements in dense molecular clouds in the ISM. The main mechanisms reducing the dust mass in a galaxy are the destruction by SN blastwaves, the formation of new stars (astration) and the ejection of metal-enriched gas. All these processes account for the IMF and the varying lifetimes of stars with different masses. We refer the reader to Hutter et al. (2023) and Dayal et al. (2022) for the exact formalism that derives the dust mass of each galaxy.

For each galaxy, we follow Hutter et al. (2023) to derive the attenuation of UV continuum photons (λ ≃ 1500 Å) by dust from its dust mass, virial radius and spin. The dust mass, M_d, its distribution radius, r_d, and the radius, a = 0.03 μm, and material density, s = 2.25 g cm⁻³, of the graphite/carbonaceous dust grains provide an estimate for the optical depth to the UV continuum photons,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{UV},\mathrm{c}}=\frac{3{\sigma }_d}{4as},\end{array}$$(7)

where ${\sigma }_d=M_d/(\pi r_d^2)$ is the dust surface mass density. We assume the gas and dust in each galaxy to be perfectly mixed (r_d = r_g) where the radius of the gas is given by r_g = 4.5λr_vir[(1+z)/6]^1.8. Here, λ and r_vir are the spin parameter and virial radius of the simulated halo, respectively, while the third factor reflects the redshift evolution of the compactness of galaxies and is motivated by the non-evolution of [CII] sizes from z ≃ 7 to z ≃ 4.5 (Fudamoto et al. 2022). We derive the escape fraction of the UV continuum photons by assuming a slab-like geometry,

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{c}}=\frac{1-exp(-{\tau }_{\mathrm{UV},\mathrm{c}})}{{\tau }_{\mathrm{UV},\mathrm{c}}}.\end{array}$$(8)

The observed UV luminosity is then given by the dust attenuated intrinsic UV luminosity

$$\begin{array}{c}\hfill L_{\mathrm{c}}^{\mathrm{obs}}=f_{\mathrm{esc}}^{\mathrm{c}}{L}_{\mathrm{c}},\end{array}$$(9)

with the intrinsic UV luminosity L_c being computed as described in Hutter et al. (2021).

2.3. Incorporating an evolving IMF

In this work, we extended ASTRAEUS to include an IMF that evolves with redshift and depends on the metallicity of the star-forming gas as found in the hydrodynamics simulations presented in Chon et al. (2022). These simulations follow the collapse of a spherical turbulent gas cloud in the presence of the CMB at different redshifts and gas metallicities. They suggest that the IMF can be composed of a component with a Salpeter-like slope at the low stellar mass end and a log-flat component at the high stellar mass end, with the fraction of stellar mass produced in the log-flat component increasing with rising redshift and decreasing gas metallicity – originating from the CMB heating the gas to higher temperatures at higher redshifts and gas cooling becoming less efficient for metal-poorer gas. By fitting the results from these simulations, Chon et al. (2022) derives the respective scaling of the fraction of stellar mass in the log-flat component as

$$\begin{array}{cc}\hfill f_{\mathrm{massive}}& =1.07\ast (1-2^x)+0.04\ast 2.{67}^x\ast z,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill x& =1+lg\left(\frac{Z}{Z_{\odot }}\right),\hfill \end{array}$$(11)

with Z_⊙ = 0.0134 (Asplund et al. 2009), and z and Z being the redshift and gas metallicity of the star-forming cloud, respectively. To efficiently implement such an evolving IMF into the ASTRAEUS simulation framework, we define the IMF, $\varphi (M)=\frac{\mathrm{d}N}{\mathrm{d}M}$, as

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}M}=\{\begin{array}{cc}(1-f_{\mathrm{massive}})\frac{M_{\mathrm{i}}^{-0.35}-M_{\mathrm{f}}^{-0.35}}{M_{\mathrm{i}}^{-0.35}-M_{\mathrm{c}}^{-0.35}}{M}^{-2.35},\hfill & M_{\mathrm{i}}\le M<M_{\mathrm{c}}\hfill \\ f_{\mathrm{massive}}\frac{M_{\mathrm{f}}-M_{\mathrm{i}}}{M_{\mathrm{f}}-M_{\mathrm{c}}}{M}^{-1},\hfill & M_{\mathrm{c}}\le M\le M_{\mathrm{f}},\hfill \end{array}\end{array}$$(12)

where M_i and M_f are 0.1 M_⊙ and 100 M_⊙ respectively, M_c the cutoff mass that describes the mass where the Salpeter component transitions to the log-flat component, and f_massive the stellar mass fraction above M_c⁵. We show how the IMF and f_massive depend on redshift and gas metallicity in Fig. 1. In practice, we derive f_massive from Eq. (11), while we infer M_c by solving ${\int }_{M_{\text{i}}}^{M_{\text{f}}}\varphi (M) M \text{d}M=1M_{\odot }$ numerically. To speed up computations, we generate a corresponding look-up table for M_c that covers f_massive values ranging from 0 to 1 in steps of 0.001. We note that we choose an upper mass limit of M_f = 100 M_⊙ to be consistent with the parameterised upper mass limit of the IMF in Chon et al. (2022) as well as previous ASTRAEUS simulation runs.

Fig. 1. Evolving IMF compared with the Salpeter IMF, and the evolution of the massive fraction with redshift and metallicity. Left: evolving IMF at different absolute metallicities Z = 10−2 (purple lines), 10−3 (orange lines), 10−4 (green lines) and redshifts z = 5 (dashed lines), 10 (dot-dashed lines), and 15 (solid lines). For clarity, the evolving IMFs at z = 10 and 15 have been displaced by lg(ϕ) = 2 and 4, respectively. Thin black lines show the Salpeter IMF at the depicted redshifts. Right: massive fraction fmassive depicted as a function of redshift z and absolute metallicity Z. The top-left corner shows the evolving IMF developing into the Salpeter IMF, where fmassive = 0.

Several physical processes implemented in ASTRAEUS, such as SN feedback, metal enrichment and the ionising emissivity, depend on the assumed IMF. A key challenge in incorporating an evolving IMF is that these processes depend on the evolution of the IMF throughout the star formation history (SFH) of a galaxy. While our implementation of these processes can easily account for the redshift dependence of the IMF, tracking the gas metallicities of all progenitors of each galaxy would require advanced book-keeping structures in the simulation code. Therefore, for a galaxy at redshift z, we assume that its progenitors at redshift z_j formed stars with a metallicity mass-averaged over all N_p progenitors at z_j,

$$\begin{array}{cc}\hfill Z(z_j)& =\frac{{\sum }_{p=1}^{N_p}{M}_{\mathrm{metal}}^i(z_j)}{{\sum }_{p=1}^{N_p}{M}_{\mathrm{g}}^i(z_j)}.\hfill \end{array}$$(13)

Here, $M_{\text{metal}}^i$ and $M_{\text{g}}^i$ are the initial metal and gas masses of the progenitor at z_j, respectively. While this assumption becomes less accurate as the number of progenitors increases, it allows us to assign a single metallicity history to each simulated galaxy, making the implementation of the evolving IMF straightforward. We plan to develop a more accurate representation in future work, such as storing the metallicities of the progenitors in the last 30 Myrs.

Thus, in our delayed SN feedback scheme, the fraction of stellar mass forming and exploding in the current time step, ν_z, in Eq. (4) is derived from the IMF that is determined by the galaxy’s current redshift z and gas metallicity. However, the fraction of stellar mass forming in previous time steps and exploding in the current time step, ν_j, is determined by the redshift z_j of the previous time step and the metallicity mass-averaged over all progenitors at z_j.

We also account for the evolving IMF when deriving the metal yields from SNII and SNIa explosions as well as AGB stars by accounting for the evolving IMF, ϕ(M), in Eq. (8) in Ucci et al. (2023). For SNIa, we adjusted the fraction of stellar systems in the mass range 3 − 16 M_⊙, f_3 − 16, representing 2.8% of SNIa progenitors, and the number of stars produced per 1 M_⊙, k, to the evolving IMF. We adjusted the second term in the summation of Eq. (8) in Ucci et al. (2023) to $Ak{\int }_{\tau (0.85M_{\odot })}^{\tau (8M_{\odot })}{M}^{\mathrm{SNIa}}{f}_{3-16}(t-\tau )\mathrm{SFR}(t-\tau )\mathrm{DTD}(t)\mathrm{d}\tau$ with the analytic power-law SNIa delay-time distribution DTD as defined in Maoz et al. (2012) and Yates et al. (2013).

To derive the ionising emissivity and UV luminosity of each galaxy, we have used the stellar population synthesis code STARBURST99 (Leitherer et al. 1999) to generate the spectra for a starburst assuming a Salpeter IMF between 0.1 M_⊙ and M_c and a log-flat IMF between M_c and 100 M_⊙, for M_c values ranging from 0 − 100 M_⊙ and metallicity values ranging from Z = 0.001 to 0.008. From the STARBURST99 simulation outputs we fit the time evolutions of the ionising emissivity and UV luminosity as a function of redshift and metallicity, with the latter two defining the IMF. For the ionising emissivity, we have:

$$\begin{array}{cc}& lg\frac{\dot{Q}(t,z,Z)}{{\mathrm{s}}^{-1}{M}_{\odot }^{-1}}=45.93-1.17(lgZ+2)+0.064z(lgZ+3.46)\hfill \\ \hfill & -\max [0,lgt-6.35]\hfill \\ \hfill & \times \left([3.91+z((lgZ+2.97)\ast 0.051)+0.0293]\right).\hfill \end{array}$$(14)

Then, for the UV luminosity, we have:

$$\begin{array}{cc}& lg\frac{L_{\nu }(t,z,Z)}{\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}{M}_{\odot }^{-1}}=31.94-0.92(lgZ+1.68)\hfill \\ & +0.06z(lgZ+3.50)\hfill \\ \hfill & +0.69\max [0,\min (6.45,lgt)-6.0]\hfill \\ \hfill & -(2.0-0.032z)\max [0,lgt-7.32]\hfill \\ \hfill & -(1.56+0.06z)\hfill \\ \hfill & \times \max [0,\min (7.32,lgt)-6.45].\hfill \end{array}$$(15)

We note that since the ionising emissivity of a galaxy is computed on the fly, we used the galaxy’s mass-averaged metallicities as defined in Eq. (13). As we derive a galaxy’s UV luminosity in post-processing, we account for the exact gas metallicities of its progenitors. We also assume that the newly formed stellar mass, $M_{\star }^{\text{new}}$, is formed at a constant rate throughout the time step. This requires convolving Eqs. (14) and (15) with the respective top-heavy star formation history; we provide the resulting expressions in Appendix A.

2.4. Reionisation

ASTRAEUS follows the reionisation of the IGM on the fly. At each time step, it derives the number of ionising photons generated in each galaxy by convolving the star formation and metallicity histories with the evolving ionising emissivities of the respective stellar populations (see Eq. (14) for the ionising emissivities’ redshift and metallicity dependencies). The number of ionising photons that escape from each galaxy to contribute to the ionisation of the IGM is given by:

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{ion}}=f_{\mathrm{esc}}\dot{Q}.\end{array}$$(16)

Here, f_esc represents the fraction of ionising photons produced in the galaxy that escape into the IGM. The MHDEC model in Hutter et al. (2023) and the results in Ocvirk et al. (2021) suggest that an f_esc value effectively decreasing with halo mass provides a better fit to the global H I fraction at z < 6. For this reason, we adopt the physically motivated f_esc model from Hutter et al. (2021) where f_esc scales with the ejected gas fraction, i.e. $f_{\text{esc}}=f_{\text{esc}}^0\min (1,f_{\star }^{\text{eff}}/f_{\star }^{\text{ej}})$ with $f_{\text{esc}}^0$ being a free parameter. Using the resulting ionising emissivities Ṅ_ion of each galaxy at the time step and the distribution of the intergalactic gas density (assumed to follow the DM perfectly), ASTRAEUS computes the spatial distribution of the ionised regions in the simulation box. To determine whether a region within the simulation box is ionised, it compares the region’s cumulative number of ionising photons with its number of absorption events at different smoothing scales (see Hutter 2018, for the details of the algorithm), accounting in this way for the ionising radiation from neighbouring sources. Within ionised regions, ASTRAEUS also computes the photoionisation rate and the residual neutral hydrogen (H I) fraction, outputting photoionisation and ionisation fields. From these outputs, it determines (i) when the environment of a galaxy was reionised and (ii) the corresponding photoionisation rate. Both are then used to account for the radiative feedback from reionisation by computing the gas mass the galaxy can retain ($f_g{M}_g^i$). For the latter, we adopt the Photoionisation model described in Hutter et al. (2021), representing a weak to intermediate, time-delayed radiative feedback. In this model, a galaxy can only hold on to its gas if it exceeds a characteristic mass M_char. In turn, M_char increases as the photoionisation rate Γ_HI incident at the galaxy’s location at z_reion and/or the difference between z_reion and the galaxy’s current redshift z rises (see also Sobacchi & Mesinger 2013), where z_reion marks the redshift at which the environment of the galaxy became reionised.

3. Baselining the model against observed data sets

In the following, we analyse and compare the simulation results for two IMF models: (i) a constant Salpeter IMF ranging from 0.1 M_⊙ to 100 M_⊙, and (ii) an evolving IMF ranging from 0.1 M_⊙ to 100 M_⊙ and becoming increasingly top-heavy towards higher redshift and lower metallicities, characterised (as described in Sect. 2.3) by shifting the transition from a Salpeter to a log-flat slope to higher star masses (≤100 M_⊙). For both IMF models, we adjusted the three free model parameters in ASTRAEUS. The normalisation of the ionising escape fraction, $f_{\text{esc}}^0$, was adjusted to fit the IGM neutral hydrogen fraction constraints from Lyman-α emitters, quasar absorption spectra and gamma ray bursts (GRBs), and the optical depth measured by Planck (Planck Collaboration VI 2020), while the maximum star formation efficiency, f_⋆ and the SN wind coupling efficiency, f_w, were adjusted to reproduce the observed dust-attenuated ultraviolet (UV) luminosity functions (LFs) at z = 5 − 12.5. We chose these eight UV LFs to fit the galaxy model parameters f_⋆ and f_w. However, the more abundant observational constraints at lower redshifts (z ≲ 9) effectively determine these two redshift-independent parameter values; therefore, the UV LF at z = 12.5, which is of most interest, is not constraining the model. We go on to detail how we found the best-fit parameters and discuss their implications.

3.1. UV luminosity functions

For both IMF models, we derive the intrinsic and observed UV LFs from our simulations by computing the intrinsic (L_c) and dust-attenuated UV luminosities ($L_{\text{c}}^{\text{obs}}$) for each galaxy in our simulation box, as detailed in Sect. 4.3. The simulated observed UV LFs are then compared to the observational results at z = 5 − 12.5 to constrain the maximum star formation efficiency f_⋆, and SN feedback wind coupling efficiency f_w. The best-fit values for each of our two IMF model were found according to the following steps.

At each redshift, we derive a ${\chi }_z^2$ characteristic, which we calculate from the linearly interpolated logarithmic number density values of the simulated and observed UV LFs and assume Poisson uncertainties on the simulated values. We then minimise the total χ² value across all redshifts, ${\chi }_{\text{total}}^2=\sum _z{\chi }_z^2$ to identify the f_⋆ and f_w best-fit values listed in Table 1⁶. Although our χ² fitting procedure considers constraints from all redshifts, we stress that our results are biased to lower redshifts (z ≤ 9) where observational constraints are tighter. While increasing f_w suppresses the star formation in lower mass galaxies (M_h ≲ 10^10.3 M_⊙ and 10^10.4 M_⊙ for the Salpeter IMF and evolving IMF model at z = 6), thus repressing and flattening the faint end of the UV LF, increasing f_⋆ raises effectively the star formation in more massive galaxies shifting the bright end of the UV LF to higher UV luminosities. Since the evolving IMF model’s top-heavier IMF leads to a higher abundance of massive stars and, thus, a lower mass-to-UV light ratio, we need to reduce its star formation efficiency compared to the Salpeter IMF to fit the observed UV LFs. We identify a 2.5× lower star formation efficiency of f_⋆ = 0.01 and a 1.5× higher SN wind coupling efficiency of f_w = 0.3 than the Salpeter IMF model’s values of f_⋆ = 0.025 and f_w = 0.2.

We show the corresponding intrinsic (dashed lines) and observed (solid lines) UV LFs at z = 6 − 11, 12.5 and 15 for the evolving (blue lines) and Salpeter (red lines) IMF models along with observational data (grey and black points) in Fig. 2. Overall, both IMF models agree within the observational uncertainties with the observations at z = 5 − 11 but somewhat underestimate the abundance of M_UV < −19 galaxies at z = 12.5. Interestingly, despite the redshift dependency of the evolving IMF model, we find the observed UV LFs of both IMF models to be very similar across all redshifts; they consistently differ by no more than approximately 0.2 − 0.3 dex within the range where we possess sufficient statistics (lg(N)≳ − 5.8 with 10 galaxies per bin and a bin width of 0.5 mag) and star formation rate histories (SFHs) have converged⁷. One would intuitively expect that the UV LF of the evolving IMF evolves less significantly at higher redshifts than the UV LF of the Salpeter IMF, since its IMF becomes increasingly top heavy, reducing the mass-to-UV light ratio. However, the rising abundance of massive stars towards higher redshifts results in both a more efficient suppression of star formation through stronger SN feedback and a lower mass-to-light ratio, explaining why the redshift evolution of the UV LFs of the two IMF models are so similar.

Fig. 2. UV luminosity functions (LFs) at z = 6 − 11, 12.5 and 15 for the evolving IMF (blue lines) and the Salpeter IMF (red lines). Dashed lines show the intrinsic UV LFs, while solid lines depict the dust attenuated UV LFs. Vertical black solid and dashed lines denote the UV luminosity thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the star formation rates have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Black points show the observational constraints from JWST data: Adams et al. (2024) (hexagons), Bouwens et al. (2023a) (diamonds), Bouwens et al. (2023b) (squares), Casey et al. (2024) (triangles), McLeod et al. (2024) (x’s), Donnan et al. (2023) (circles), Finkelstein et al. (2023) (stars), Willott et al. (2024) (pentagons). Light and dark grey points depict the observational constraints before JWST: Atek et al. (2015) (light stars), Atek et al. (2018) (dark stars), Bouwens et al. (2015) (light squares), Bouwens et al. (2016) (dark squares), Bouwens et al. (2017) (light circles), Bouwens et al. (2021) (dark circles), Bowler et al. (2020) (light x’s), Calvi et al. (2016) (dark x’s), Finkelstein et al. (2015) (light diamonds), Ishigaki et al. (2018) (dark diamonds), Livermore (2016) (light ‘downward’ triangles), Livermore et al. (2017) (dark ‘downward’ triangles), McLure et al. (2009) (light ‘upward’ triangles), McLure et al. (2013) (dark ‘upward’ triangles), Schenker et al. (2013) (light ‘left’ triangles), Schmidt et al. (2014) (dark ‘left’ triangles), van der Burg et al. (2010) (light ‘right’ triangles), Willott et al. (2013) (dark ‘right’ triangles).

From Fig. 2, we can see that the UV LF reflects the hierarchical growth of galaxies. Firstly, as cosmic time progresses, galaxies grow in mass through mergers and gas accretion, moving the UV LF to higher luminosities. For instance, while the most luminous galaxies at z = 15 have UV luminosities of M_UV ≃ −20, the brightest at z = 6 are about 40× brighter reaching M_UV ≃ −24. Secondly, more low-mass galaxies form with cosmic time. These galaxies (i) exhibit bursty star formation caused by SN feedback and radiative feedback from reionisation (most effective in low-mass halos) and (ii) are likely to be consumed by mergers in the vicinity of more massive galaxies. These processes make the redshift evolution of the UV LF’s faint end more complex, effectively comprising a combination of positive and negative luminosity and number density evolution, with low-mass galaxies brightening and fading as well as forming and being consumed by merging. However, overall the number density of faint galaxies with M_UV = −16 increases by about a factor 30 from about 10⁻³ Mpc⁻³ at z = 15 to 3 × 10⁻² Mpc⁻³ at z = 6 for both IMF models. Moreover, as both SN and radiative feedback affect more massive galaxies with decreasing redshift and gas-poor low-mass galaxies can attain fainter UV luminosities, the slope of the faint end of the UV LF flattens. At the same time, the bright end steepens with decreasing redshift, as the UV radiation of more massive and thus more luminous galaxies is increasingly attenuated by dust⁸.

Finally, we briefly discuss our models’ predictions on the evolution of the UV LFs at z ≳ 10. We find the slope of the bright end (M_UV ≲ −19) to remain essentially constant, changing by no more than a few per cent, and the number density at M_UV ≳ −20 to decrease by less than 1 dex from z = 10 to z = 12.5. This minor redshift evolution is in rough agreement with the observational findings (e.g. Harikane et al. 2023a,b; Ferrara et al. 2023; Castellano et al. 2022) and continues up to z = 15. However, at z = 12.5, both our IMF models slightly underpredict the abundance of M_UV < −19 galaxies. Should spectroscopic follow-up observations validate these abundances and the prevalence of bright galaxies show only a minor decline towards even higher redshifts, as photometric detections of z ≳ 13 galaxies imply (Donnan et al. 2023; Finkelstein et al. 2023; Harikane et al. 2023a), we would need to increase our model’s star formation efficiency towards higher redshifts (e.g. Dekel et al. 2023) or invoke contributions from efficiently accreting massive black holes (e.g. Pacucci et al. 2022).

In summary, despite the evolving IMF becoming more top-heavy towards higher redshifts, we do not observe any alteration in the evolution of its UV LF compared to the Salpeter IMF. Any reduction in the mass-to-UV light ratio towards higher redshifts is compensated by stronger SN feedback, which leads to an increased suppression of star formation in addition to the overall assumed lower star formation efficiency.

3.2. Reionisation history

For both the Salpeter IMF and the evolving IMF models, we constrained the third free model parameter, $f_{\text{esc}}^0$, within our f_esc parameterisation, with the observational constraints on the volume-averaged ionisation history from quasar absorption lines, GRBs and Lyman-α emitters, and the electron optical depth derived from the CMB measurements with Planck (Planck Collaboration VI 2020). Here, we have prioritised obtaining the most similar evolutions of the volume-averaged ionisation fraction at z ≃ 7 for both IMF models. The corresponding values for $f_{\text{esc}}^0$ and evolutions of the global H I fraction ⟨χ_HI⟩(z) are shown in Table 1 and Fig. 3, respectively.

Fig. 3. Reionisation history. Top: ratio of the mass- and volume-averaged neutral hydrogen fraction. Bottom: volume averaged neutral hydrogen fraction as a function of redshift. In each panel, we show results for the Salpeter IMF (dash-dotted red line) and the evolving IMF (solid dark blue line) models. Grey points indicate observational constraints from: GRB optical afterglow spectrum analyses (light pentagons; Totani et al. 2006, 2014), quasar sightlines (faint squares; Fan et al. 2006) (medium bright circles; Bosman et al. 2022), Lyman-α LFs (dark upward triangles; Konno et al. 2018; Kashikawa et al. 2011; Ouchi et al. 2010; Ota et al. 2010; Malhotra & Rhoads 2004), Lyman-α emitter clustering (dark downwards triangles; Ouchi et al. 2010), the Lyman-α emitting galaxy fraction (dark leftwards triangles; Pentericci et al. 2011; Schenker et al. 2012; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014; Pentericci et al. 2014), Lyman-α equivalent widths (dark rightwards triangles; Mason et al. 2018, 2019; Bolan et al. 2022), dark pixels (light squares; McGreer et al. 2015), and damping wings (light diamonds; Davies et al. 2018; Greig et al. 2019).

Firstly, we find the $f_{\text{esc}}^0$ value for the evolving IMF model (0.038) to be about 8× lower than for the Salpeter IMF (0.31), decreasing the f_esc values for all halo masses. The evolving IMF model’s decrease in f_esc values arises mainly from its higher ionising emissivities per stellar mass stemming from the higher abundance of massive stars generated by its top-heavier IMF.

Secondly, although the derived $f_{\text{esc}}^0$ values lead to very similar ionisation histories for both IMF models, we find the evolving IMF model to exhibit a higher electron optical depth (τ = 0.0580) than the Salpeter IMF model (τ = 0.0543). This difference is due to the electron optical depth tracing the mass-averaged and not volume-averaged ionisation fraction. In the evolving IMF model, the spatial distribution of the ionising emissivity emerging from the underlying galaxy population follows the underlying gas density distribution more closely (see the ratio of the global mass- and volume-averaged H I fractions in the top panel of Fig. 3), shifting the ionisation history to lower redshifts at fixed electron optical depth.

Finally, we note that both our IMF models underpredict the observational constraints on ⟨χ_HI⟩ at z < 5.4 by a factor ∼2 − 3. As highlighted in Ocvirk et al. (2021), reproducing these observational constraints requires the total ionising emissivity to drop towards lower redshifts at 4 ≲ z ≲ 6, driven by the strong radiative suppression of star formation in low-mass halos and the escape fractions steeply declining with halo mass. While our ASTRAEUS model includes radiative feedback from reionisation that significantly suppresses star formation in low-mass halos at z ≃ 6 (cf. Fig. 5 in Hutter et al. 2021), our models’ too small ⟨χ_HI⟩ values indicate that (i) our parameterisation of f_esc needs to reflect an even steeper decline towards higher halo masses, or (ii) our model of radiative feedback from reionisation needs to be stronger or less delayed, or (iii) we miss the dense small-scale gas distribution. We will explore this in detail in future work.

4. Effects of an evolving IMF on galaxy evolution and reionisation

We go on to a discussion of how going from a constant Salpeter IMF to an evolving IMF changes the relations between galaxy properties, including their star formation rates, stellar masses, gas metallicities, and dust masses.

4.1. Main characteristics of an evolving IMF in galaxy properties

First, to lay the foundation for understanding how all fundamental relations between galaxy properties change for an evolving IMF, we describe how such a change in the IMF alters the build-up of stellar mass (Sect. 4.1.1) and the light-to-stellar mass ratio (Sect. 4.1.2).

4.1.1. Shift to lower stellar masses

As mentioned in Sect. 3.1, we find that the evolving IMF model with its lower mass-to-light ratio requires overall lower star formation efficiencies than the Salpeter IMF model to reproduce the observed UV LFs. An important consequence of these lower star formation efficiencies is the slower buildup of stellar mass. As we can see from the relation between the stellar and halo masses in Fig. 4, this results in a lower stellar-to-halo mass ratio, by ∼0.4 − 1 dex, for the evolving IMF model across all redshifts.

Fig. 4. Redshift evolution of the stellar to halo mass relation. Left: relation between the stellar and halo mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of stellar masses (coloured contours) at a given halo mass. Horizontal black solid and dashed lines denote the stellar mass thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the stellar-to-halo mass ratio for different halo mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each halo mass bin.

The extent to which the median stellar mass is lower in the evolving IMF compared to the Salpeter IMF model for halos of the same mass exhibits two trends. Firstly, it increases towards lower halo masses, from ∼0.5 dex for M_h ≳ 10¹¹ M_⊙ to ∼1 dex for M_h ≃ 10^9.5 M_⊙ at z = 6 − 12 (cf. left panels in Fig. 4). Secondly, this difference decreases towards lower redshifts, particularly in more massive galaxies (M_h ≳ 10¹⁰ M_⊙) with the stellar-to-halo mass ratio (cf. right panel in Fig. 4) dropping by at least 0.1 dex from z = 12 to 5. The first trend stems from SN feedback in the evolving IMF model becoming more efficient in suppressing star formation in galaxies located in shallower gravitational potentials (i.e. lower halo masses). This increased efficiency is due to our evolving IMF becoming more top-heavy, characterised by a higher abundance of massive stars undergoing SN events, towards lower gas metallicities and, consequently, lower halo masses. The second trend reflects that our evolving IMF becomes less top-heavy with cosmic time, resulting in less efficient SN feedback and increased star formation (despite the gravitational potentials of halos becoming shallower).

We also note that the decrease in the stellar-to-halo mass ratio with decreasing redshift for lower halo masses (M_h ≲ 10^10.5 M_⊙) is due to gravitational potentials becoming shallower as the Universe expands, causing star formation to be SN feedback suppressed in galaxies with increasingly higher halo masses. In contrast, more massive galaxies where star formation is not SN feedback limited show a roughly constant stellar-to-halo mass ratio across all redshifts. We find the transition where a galaxy’s star formation is not affected by SN feedback to occur around M_h ∼ 10^10.3 M_⊙ (10^9.9 M_⊙) and 10^10.4 M_⊙ (10^10.1 M_⊙) for the Salpeter IMF and evolving IMF models at z = 6 (9), highlighting that the larger fraction of stars exploding as SNe is counteracted by the overall lower star formation efficiency in the evolving IMF model.

The shift of stellar masses to lower values in the stellar-to-halo mass relation is one of the main features of the evolving IMF. This characteristic explains to a first-order how the relations between the main galaxy properties change when transitioning from a Salpeter IMF to an evolving IMF, which we will discuss in the following sections in detail.

4.1.2. Reduced mass-to-light ratio

The second main feature of our evolving IMF is its lower galaxy-wide mass-to-light ratio compared to the constant Salpeter IMF. In Fig. 5, we show the observed UV luminosity-to-stellar mass relations for both IMF models at redshifts 6, 8, and 12 (left panels), as well as the redshift evolution of the mass-to-light ratios of galaxies at fixed stellar masses (right panel).

Fig. 5. Redshift evolution of the observed UV luminosity to stellar mass relation. Left: relation between the observed UV luminosity and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of observed UV luminosities (coloured contours) at a given stellar mass. Black solid and dashed lines denote the stellar mass (vertical) and UV luminosity (horizontal) thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses and observed UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Median values outside the corresponding contours reflect low number statistics. Right: redshift evolution of the observed UV luminosity-to-stellar mass ratio for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each stellar mass bin.

Overall, galaxies with the same stellar mass are about 0 − 2 mag brighter in the evolving IMF than in the Salpeter IMF model across all stellar masses and redshifts. The reason for this lies in the evolving IMF being top-heavier and thus forming a higher fraction of massive stars. As massive stars dominate the UV luminosity, galaxies of the same stellar mass contain more massive stars in the evolving IMF model and thus appear brighter.

However, we find the enhancement of the median UV luminosity at a given stellar mass to vary across stellar masses and redshifts. While the difference in the median UV luminosity at a fixed stellar mass remains around 2 − 3 mag for lower mass galaxies (M_⋆ ≲ 10⁷ M_⊙), it decreases for more massive galaxies (M_⋆ ≳ 10⁷ M_⊙) towards higher stellar masses with this decrease becoming more pronounced towards lower redshifts. For instance, at z = 12, this difference drops from ∼2 mag for M_⋆ = 10⁷ M_⊙ galaxies to ∼1 mag for M_⋆ = 10^8.5 M_⊙ galaxies, while it drops from ∼3 mag for M_⋆ = 10⁷ M_⊙ galaxies to ∼0 mag for M_⋆ = 10⁹ M_⊙ galaxies at z = 6. Similarly, we find the enhancement of the median mass-to-light ratio of M_⋆ = 10⁹ M_⊙ galaxies to decrease from ∼0.7 dex at z = 12 to 0 dex at z = 6 (cf. right panel).

These trends reflect that the evolving IMF becomes more top-heavy towards lower gas metallicities and higher redshifts (see Fig. 1). Firstly, due to having their star formation suppressed and their gas expelled by SN feedback, lower mass galaxies exhibit lower gas metallicities across all redshifts, which causes their top-heavy IMF and thus the low mass-to-light ratio of their newly forming stars to hardly evolve with cosmic time (cf. Fig. 1)⁹. This explains the evolving IMF’s nearly constant negative offset in median UV luminosities to the Salpeter IMF for lower mass galaxies (M_⋆ ≲ 10⁷ M_⊙). Secondly, as galaxies become more massive, they can retain their gas and metals, increasing their gas metallicities as they accumulate stellar mass. As a consequence, their IMFs become increasingly Salpeter-like. The resulting increasing mass-to-light ratio with rising stellar mass explains why the difference in the median UV luminosities between the evolving IMF and Salpeter IMF models decreases with increasing stellar masses. This decrease becomes more pronounced towards lower redshifts because the trend of the evolving IMF approaching the Salpeter IMF becomes more distinct. The latter is due to the redshift evolution of the evolving IMF, which is the strongest for high metallicities.

These metallicity and redshift dependencies of the evolving IMF also affect the redshift evolution of the absolute mass-to-light ratio values, as shown in the right panel of Fig. 5. While the mass-to-light ratio increases with decreasing redshifts across all stellar masses for both models, the evolving IMF model exhibits a steeper increase with redshift than the Salpeter model. Generally, an increase in the mass-to-light ratio towards lower redshifts is expected: galaxies accumulate stellar mass as they evolve, but their UV luminosities are dominated by the short-lived (≲30 Myr) massive stars and thus only trace the relatively newly forming stars. Given that our evolving IMF approaches the Salpeter IMF towards lower redshifts and higher stellar masses, the abundance of UV-bright short-lived massive stars drops. This drop in the massive star abundance reduces the mass-to-light ratio of the newly forming stars and thus the overall mass-to-light ratios of galaxies as they grow with cosmic time.

Finally, it is important to note that the evolving IMF model shows a larger scatter in the UV luminosity values for a given stellar mass than the Salpeter IMF model. For instance, while the UV luminosities for galaxies with M_⋆ ∼ 10⁸ M_⊙ range only from M_UV ∼ −12 to −19 at z = 6 in the Salpeter IMF model, they extend up to M_UV ∼ −21 in the evolving IMF model. There are two reasons for this enhanced scatter. Firstly, since the more top-heavy IMF in the evolving IMF model produces a higher abundance of massive stars, fewer stars are left over after the short-lived massive stars die, giving rise to UV-fainter galaxies between consecutive starburst events. Secondly, as the evolving IMF is sensitive to the metallicity, with metal-poor galaxies having a more top-heavy IMF than metal-rich galaxies, any scatter in the galaxies’ metallicity values due to their different accretion and merger histories enhances the scatter in the UV.

In summary, the main characteristics of an IMF that becomes increasingly top-heavy towards higher redshifts and lower metallicities are: (1) the slower build-up of stellar mass resulting in a shift to lower stellar masses in the stellar-to-halo mass relations and (2) the reduced mass-to-light ratio due to the higher abundance of massive stars.

4.2. Star formation main sequence

We next discuss how the star formation main sequence depends on the assumed IMF model. For this purpose, we show the relation between the star formation rate (SFR) and stellar mass at z = 6, 8 and 12 and the redshift evolution of the specific star formation rate sSFR = SFR/M_⋆ of galaxies with fixed stellar masses in Fig. 6. We derive the SFR of each galaxy in our simulations by evaluating how much stellar mass formed in the current time step, $M_{\star }^{\text{new}}/\mathrm{\Delta }t$ with Δt being the current time step’s length.

Fig. 6. Redshift evolution of the SFR to stellar mass relation. Left: relation between SFR and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of SFRs (coloured contours) at a given stellar mass. Vertical black solid and dashed lines denote the stellar mass threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the specific SFR for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The value in the panel shows the central value of each stellar mass bin.

Firstly and most notably, we find that the median star formation main sequence remains unchanged for galaxies where star formation is driven by gas accretion and not SN feedback limited (M_⋆ ≳ 10⁷ M_⊙), when transitioning from the Salpeter IMF to the evolving IMF model. This independence of the SFR – stellar mass relation from the chosen IMF model can be explained as follows: The build-up of stellar mass depends on (i) the SFR and (ii) the dark matter assembly history. For each IMF, the SFR is essentially constrained by the UV luminosities of the observed galaxy populations; that is, for a more top-heavy IMF the lower stellar mass-to-observed UV luminosity ratio is compensated by a lower star formation efficiency to reproduce the observed UV LFs as described in Sect. 3.1. Since the build-up of dark matter and thus gas available for star formation in galaxies is self-similar (i.e. Ṁ_h/M_h ≃ const.) the fractional increase in halo or gas mass is essentially independent of the halo mass. This self-similarity has the consequence, that for different assumed star formation efficiencies (f_⋆), galaxies of similar stellar masses will have different halo masses, but similar SFRs. In scenarios where the star formation efficiency is reduced, galaxies of a given stellar mass are located in more massive halos where more gas is available for star formation, such that the stellar mass gained (SFR ∝ f_⋆M_g) remains the same as in a scenario with a higher star formation efficiency where galaxies of the same stellar mass reside in less massive halos. While this picture holds for the median star formation main sequence relation, it ignores the halo-mass-dependent scatter of the DM assembly histories that cause the scatter in SFR values. As, on average, DM assembly histories become less diverse towards higher halo masses, we find that the scatter in SFR values in scenarios with a lower star formation efficiency is smaller. This dependence explains why the scatter in SFR values is smaller in the evolving IMF than in the Salpeter IMF model.

Moreover, the SFR in galaxies with lower stellar masses (M_⋆ ≲ 10⁷ M_⊙) is limited by SN feedback and therefore no longer scales with the gas accretion rate, which causes the median SFRs of our IMF models having different SN feedback efficiencies to diverge. The transition at which SN feedback becomes significant in suppressing star formation is visible for the Salpeter IMF; its median SFR at z = 6 (12) starts to drop more below M_⋆ ≲ 10^7.5 M_⊙ (10^6.5 M_⊙), corresponding to M_h ≲ 10¹⁰ M_⊙ (10^9.5 M_⊙). However, for the evolving IMF, we find this drop to be shifted to stellar masses below the convergence limit (see Appendix B in Hutter et al. 2021)¹⁰. Instead of dropping abruptly, the median SFR-stellar mass relation in the evolving IMF model becomes steeper towards lower redshifts at M_⋆ ≲ 10⁷ M_⊙. At first sight, this median SFR, which exceeds that of the Salpeter IMF model, seems to be at odds with the higher SN wind coupling efficiency f_w and the higher fraction of stars turning into SN due to the more top-heavy IMF. However, a higher fraction of massive stars also implies that a higher fraction of a stellar population’s SN energy is injected sooner (i.e. in the same time step) rather than being delayed (to subsequent time steps). As a result, the SFR is less driven or affected by SN feedback from earlier stellar populations (at previous time steps) compared to the Salpeter IMF model, resulting in the evolving IMF model’s SFR histories being less extreme in its minimum and maximum SFR values and contributing to the smaller scatter in SFR values.

Finally, we find the median sSFR values in all stellar mass bins to decrease with decreasing redshift for both IMF models (see the right panel in Fig. 6), which is in agreement with previous observational and theoretical works (e.g. Labbé et al. 2013; Stark et al. 2013; Hutter et al. 2014). This decrease in the sSFR values is expected because the gravitational potentials of galaxies become shallower towards lower redshifts, which again allows SN feedback to limit star formation in more massive galaxies.

4.3. Evolution of metallicity and dust mass

During reionisation, SNe play a crucial role as primary sources of metals and thus contribute to the formation of dust. While the abundance and types of elements produced depend on the IMF of the forming stars, in this work, we focus on how the evolution of the total metal and dust content in galaxies changes as we transition from the constant Salpeter IMF to the more top-heavy evolving IMF. For this purpose, we explored how the galaxies’ gas-phase metallicity and dust masses correlate with their stellar masses and UV luminosities. While the existing literature primarily discusses the stellar mass-metallicity relation, it is important to note that inferring the stellar mass from observations requires assuming a specific IMF. Therefore, we have also included how the gas-phase metallicity depends on the observed UV luminosity – with the latter being a directly observable quantity and, thus, independent from the IMF in observational inferences. We note that such a relation will not directly probe the underlying 3D relation between stellar mass, metallicity, and SFR (e.g. Mannucci et al. 2010; Lara-López et al. 2010), as the observed UV luminosity is a proxy for dust-attenuated and not intrinsic SFR.

4.3.1. Fundamental metallicity relation

We start by discussing how the relation between the gas-phase oxygen abundance, 12 + lg(O/H), (referred to as metallicity in the following) and stellar mass, M_⋆, the mass-metallicity relation (MZR), evolves for both our IMF models. We show the respective relations at z = 6, 8 and 12, and the redshift evolution of the median metallicity for different stellar mass bins in Fig. 7. We note that since ASTRAEUS follows the oxygen abundance within a galaxy’s metal reservoir explicitly, we derive the metallicity directly from the oxygen, (M_O), and gas masses, (M_g), $lg(\mathrm{O}/\mathrm{H})=lg\left(\frac{M_{\mathrm{O}}}{M_{\mathrm{g}}}\frac{m_{\mathrm{H}}}{m_{\mathrm{O}}(1-Y)}\right),$ with m_H and m_O being the masses of hydrogen and oxygen atoms and Y as the helium mass fraction. We assume the solar metallicity from Asplund et al. (2009), yielding 12 + lg(O/H) = 8.76 as the solar metallicity normalisation. From the left panels of Fig. 7, we can see that the gas in a galaxy becomes increasingly metal-enriched as the galaxy builds up its stellar mass through gas accretion-driven star formation and mergers. For instance, at z = 6, the median oxygen-based metallicity increases from 7.8 and 7.6 at M_⋆ = 10⁸ M_⊙ to 8.1 and 8.0 at M_⋆ = 10¹⁰ M_⊙ for the Salpeter IMF and evolving IMF models, respectively. While each episode of star formation leads to an injection of metals into the ISM, only galaxies not limited by SN feedback can retain their metal-enriched gas. In our model, low-mass galaxies (M_⋆ ≲ 10⁷ M_⊙ at z = 6) have most of their gas ejected and cannot hold onto their metal-enriched gas, causing the drop of the median metallicity at M_⋆ ≲ 10⁷ M_⊙ in both IMF models. However, as galaxies become more massive, the fraction of gas and metals ejected by SN feedback decreases (see Fig. 3 in Ucci et al. 2023), causing the median gas-phase metallicity to rise with increasing stellar mass despite accreting more pristine gas from the IGM.

Fig. 7. Redshift evolution of the oxygen-based gas metallicity to stellar mass relation. Left: relation between the oxygen-based gas metallicity and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of gas-phase metallicities (coloured contours) at a given stellar mass. Black squares and circles show the observational results from Heintz et al. 2023a (z ≃ 7.1 − 9.5) and Curti et al. 2023 (z ≃ 5.8 − 9.4) converted to a Salpeter IMF stellar masses according to Speagle et al. (2014), respectively. Observations at z < 7 are shown in the z = 6 panel and observations at 7 ≤ z < 10 in the z = 8 panel. Vertical black dashed lines denote the stellar mass threshold in the Salpeter IMF simulations, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the oxygen-based gas-phase metallicity for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each stellar mass bin.

Comparing the results of the evolving IMF and Salpeter IMF models, we find that galaxies with the same stellar masses contain more metal-enriched gas in the evolving IMF model. For galaxies not limited by SN feedback (M_⋆ ≳ 10⁸ M_⊙), the enhancement in the median metallicity is ∼0.2 dex and hardly varies with redshift and stellar mass and can be explained as follows: The more top-heavy IMF in the evolving IMF model leads to a higher metal enrichment rate per stellar mass formed and to a higher abundance of oxygen than in the Salpeter IMF model, both increasing the median oxygen-based metallicity. We also note that the halo mass-metallicity relations are very similar for both IMFs, possibly because we adjust the star formation efficiencies to reproduce the observed UV LFs (leading to very similar halo mass-to-UV luminosity ratios) and the enhancement in the UV luminosity per stellar mass for a more top-heavy IMF is similar to the enhancement of the corresponding metal enrichment rate per stellar mass.

Figure 7 also shows that galaxies with the same stellar mass can have a range of metallicities, with the width of the simulated galaxies’ metallicity values decreasing towards larger stellar masses, i.e. from Δ(12 + lg(O/H)) ≃ 1.8 (1.2) at M_⋆ = 10⁸ M_⊙ to 0.8 (0.6) at M_⋆ = 10¹⁰ M_⊙ for the evolving (Salpeter) IMF model at z = 6. This scatter in gas-phase metallicity values is driven by these galaxies different SFR values: Galaxies with higher SFR values exhibit lower metallicity values (see Fig. B.1), as the associated enhanced accretion of metal-poor gas from the IGM that drives the star formation dilutes the gas and decreases the metallicity (see also Ucci et al. 2023). Moreover, as galaxies become more massive, the scatter in their SFR values decreases as their DM assembly and thus star formation histories become less diverse (i.e. a massive halo will always sit in a highly overdense region while a less massive halo can be located in an overdense or underdense region where its accretion of matter is reduced or enhanced, respectively; see also discussion for the scatter in the SFR values in Sect. 4.2).

This dependence of the gas-phase metallicity on the SFR highlights – as has also been noted in previous studies (e.g. Mannucci et al. 2010; Lara-López et al. 2010; Ucci et al. 2023) – that the MZR is a projection of the 3D fundamental metallicity relation (FMR), which connects the stellar mass, SFR and gas-phase metallicity. As demonstrated in previous works (Hunt et al. 2012, 2016; Tortora et al. 2022), the fundamental metallicity relation can be distilled to a fundamental plane of metallicity. We perform a linear regression analysis and find our data to lie on a plane in the four-dimensional parameter space of gas-phase metallicity, stellar mass, SFR, and redshift. We derive this high-z fundamental plane of metallicity (HFPZ) for both IMF models. While the one for the Salpeter IMF yields

$$\begin{array}{cc}& 12+lg(\mathrm{O}/\mathrm{H})=(0.382\pm 0.064)\ast lg\left(\frac{M_{\star }}{M_{\odot }}\right)\hfill \\ & -(0.067\pm 0.053)\ast lg\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)\hfill \\ \hfill & +(0.003\pm 0.012)\ast z+(4.421\pm 0.57),\hfill \end{array}$$(17)

the HPFZ for the evolving IMF shifts to higher metallicities for any point on the star formation main sequence:

$$\begin{array}{cc}& 12+lg(\mathrm{O}/\mathrm{H})=(0.683\pm 0.186)\ast lg\left(\frac{M_{\star }}{M_{\odot }}\right)\hfill \\ & -(0.331\pm 0.138)\ast lg\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)\hfill \\ \hfill & +(0.032\pm 0.042)\ast z+(2.113\pm 1.652).\hfill \end{array}$$(18)

We note that for both IMF models, the HFPZ depends only weakly on redshift and would also be consistent with no evolution with redshift. The reason for this weak redshift evolution lies in the star formation main sequence evolving only slightly towards higher stellar masses with decreasing redshift; the corresponding slow drop of the sSFR values and thus gas fraction M_gas/(M_gas + M_⋆) towards lower redshifts causes the slightly rising metallicities. Finally, we briefly comment on the drop in the median metallicity for galaxies with lower stellar masses (M_⋆ ≲ 10^7.5 M_⊙) across all redshifts. In the Salpeter IMF model, this drop occurs at slightly higher stellar masses compared to the evolving IMF model. For instance, at z = 6 (z = 12), the drop arises at M_⋆ ≃ 10⁷ M_⊙ (10^6.6 M_⊙) in the Salpeter model versus M_⋆ ≃ 10^7.3 M_⊙ (10⁷ M_⊙) in the evolving IMF model. This drop results from our treatment of gas, star formation, and SN feedback in low-mass galaxies, as detailed in Sect. 2. In these galaxies, either SN events have expelled all their gas, leaving them gas-depleted, or they retain gas when their star formation remains low enough to prevent SN events from ejecting all gas. As the gravitational potentials of these galaxies become shallower over cosmic time, the halo mass at which SN feedback can expel the entire gas content rises. However, despite the evolving IMF’s stronger SN feedback, its drop occurs at lower stellar masses. As explained in Sect. 4.2, its more top-heavy IMF results in a larger fraction of SN events occurring in the current time step, maintaining a low SFR without allowing it to oscillate as much as the Salpeter IMF model’s SFR.

Comparison against observational data. We compare our stellar mass – metallicity relations at z ≃ 6 − 10 with results from JWST observations, including the JWST Advanced Extragalactic Survey (JADES; Curti et al. 2024) as well as observations from the gravitational lensing clusters Abell 2744, RXJ-2129 and the Cosmic Evolution Early Release Science (CEERS) survey (Heintz et al. 2023a). The stellar masses inferred from these JWST observations depend on the IMF assumed when fitting the galaxy spectra, where Heintz et al. (2023a) used the Kroupa IMF and Curti et al. (2024) adopted the Chabrier IMF. To compare this observational data with our simulated galaxy population, we convert the “observed” stellar masses to a Salpeter IMF according to Speagle et al. (2014). At this stage, an analogous conversion to our evolving IMF is not possible, as it would require a redshift- and metallicity-dependent conversion factor obtained from corresponding spectral fitting analyses. Therefore, we only compare our Salpeter IMF model to the available data, shown as black points in the left panels in Fig. 7. Overall, we find that the stellar mass-metallicity combinations of nearly all observed galaxies¹¹ occur also in our simulated galaxy population.

Comparison against theoretical models. We briefly discuss how our stellar mass-metallicity relations at z > 5 compare to those found in other theoretical models, including the first generation of ASTRAEUS simulations (Ucci et al. 2023), ILLUSTRIS TNG (Torrey et al. 2019), FIRSTLIGHT (Langan et al. 2020), FIRE (Ma et al. 2016), and OBELISK (Trebitsch et al. 2021) (see Fig. 8). Comparing our Salpeter IMF model results with those of the first generation of ASTRAEUS simulations (Ucci et al. 2023), we found that the MZR has shifted by about 0.1 dex towards higher metallicities for M_⋆ ≳ 10⁸ M_⊙. This shift is due to the 2.5× higher star formation efficiency, which was necessary to compensate for the dust attenuation of UV light at the UV LF’s bright end that was not accounted for in Ucci et al. (2023).

Fig. 8. Comparison of the mass-metallicity relations at z = 6 derived from different simulations. Solid blue and red lines show the median values for our evolving IMF and Salpeter IMF models. Black lines show the relations from ASTRAEUS without dust (solid with 1 − σ dispersion; Ucci et al. 2023), ILLUSTRIS-TNG (dot-dot-dashed with 1 − σ dispersion; Torrey et al. 2019), FIRSTLIGHT (dotted; Langan et al. 2020), FIRE (dot-dashed; Ma et al. 2016) and OBELISK (dashed with 1 − σ; Trebitsch et al., in prep.). We have rescaled the metallicities from all these works to our solar value of 8.76.

ILLUSTRIS TNG and OBELISK serve as the most comparable simulations in terms of metal enrichment modelling (metal production from SNe and AGB stars) and modelled galaxy masses (M_⋆ ≃ 10⁶⁻¹⁰ M_⊙ and M_⋆ ≃ 10⁸⁻⁹ M_⊙), with solar normalisations corresponding to 12 + lg(O/H) = 8.6 and 8.69. Interestingly, the ILLUSTRIS TNG’s MZR aligns more with our evolving IMF’s MZR at z = 6, possibly a consequence of the ILLUSTRIS TNG simulations assuming: (i) a Chabrier IMF, resulting in a higher abundance of massive stars per solar mass formed than our Salpeter IMF (and thus a shift of ∼0.2 dex to lower stellar mass) and (ii) a lower metal loading of their winds (0.4 Z_ISM). The smaller scatter in metallicity values at fixed stellar masses in ILLUSTRIS TNG is likely a consequence of its smaller simulation box of 100h⁻¹cMpc, which does not cover more extreme mass assembly histories in over-dense or under-dense regions. Extrapolations of ILLUSTRIS TNG results suggest a more pronounced evolution of the MZR (increase by 0.25 dex from z = 10 to 6) due to the galaxies’ gas fraction M_gas/M_⋆ rising toward higher redshifts. However, similarly to the findings in Ucci et al. (2023), our simulations show no significant redshift evolution of the gas fraction, and thus, our MZR evolves very weakly with redshift.

OBELISK’s MZR, derived from the halo-averaged gas-phase metallicities (excluding dust) and stellar masses, aligns the closest with our Salpeter IMF model’s MZR compared to other non-ASTRAEUS simulations. It exceeds our Salpeter IMF model’s MZR only by 0.1 − 0.2 dex, likely due to assuming a Kroupa IMF, and exhibits a comparable slope at lower stellar masses (M_⋆ ≲ 10^8.5 M_⊙). Notably, its MZR also shows a change in slope, although at slightly larger stellar masses.

The FIRE and FIRSTLIGHT simulations are both zoom-in simulations that track galaxies with stellar masses ranging from M_⋆ = 10³⁻⁹ M_⊙ for FIRE and M_⋆ = 10^6−9.5 M_⊙ for FIRSTLIGHT. Both simulations adopt a Kroupa IMF, SNII yields from Woosley & Weaver (1995), and a solar metallicity normalisation of 12 + lg(O/H) = 8.9. Despite their different implementations of physics and simulation methods, we find the FIRSTLIGHT’s MZR to surpass our Salpeter IMF model’s MZR by only ≲0.2 dex at z = 6 − 8 for M_⋆ ∼ 10^7.5−8.5 M_⊙. Meanwhile, the MZR from the FIRE simulation lies approximately 0.3 dex below our Salpeter IMF results. Notably, FIRSTLIGHT shows a comparable redshift evolution and slope of the MZR, while the predictions from FIRE suggest a more pronounced redshift evolution of the MZR, but only based on extrapolations beyond z = 6.

4.3.2. The stellar mass to dust mass relation

We now discuss how the relation between the dust and stellar mass alters for an evolving IMF. For this purpose, we show the respective relations for the Salpeter and evolving IMF models at z = 5, 7 and 12 and the redshift evolution of the dust-to-stellar mass ratio for galaxies with fixed stellar masses in the left and right panels of Fig. 9, respectively.

Fig. 9. Redshift evolution of the dust to stellar mass relation. Left: relation between the dust and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of dust masses (coloured contours) at a given stellar mass. Black circles show the observational results at z ≃ 7 from the REBELS survey shown in Dayal et al. (2022) assuming a Salpeter IMF. Vertical black dashed lines denote the stellar mass threshold in and Salpeter IMF simulations, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the dust-to-stellar mass ratio for different stellar mass bins with a width of 1 dex. The values in the panel shows the central value of each stellar mass bin. Circles show median values of the REBELS galaxies.

Firstly, since the majority of dust is also produced in SN events, the dust mass evolves with the metal mass. Therefore, for both IMF models, the median dust mass rises with increasing stellar mass and drops at the same stellar masses towards lower stellar masses, analogous to the increasing trend of the median gas-phase metallicity with stellar mass.

Secondly, however, we can see that the dust masses in the evolving IMF model are ∼0.3 − 0.5 dex higher than in the Salpeter IMF model for M_⋆ ≳ 10^7−7.5 M_⊙ at z = 12 − 6. This shift to higher dust masses is not solely due to the ∼0.5 − 1 dex lower stellar-to-halo mass ratio discussed in Sect. 4.1.1. Unlike metals, dust – which we consider to be a subreservoir of metal mass – is also destroyed by SN shock waves. Therefore, due to the more top-heavy IMF, the evolving IMF model yields about 0.3 dex lower dust-to-metal mass ratio and, thus, lower dust masses for a given halo mass than in the Salpeter IMF model. This explains (i) why the evolving IMF model’s offset of the dust masses from those in the Salpeter IMF model is smaller than the shift towards lower stellar masses and (ii) why the evolving IMF model’s difference between the intrinsic and dust attenuated UV LFs shown in Fig. 2 are smaller.

Thirdly, for galaxies with M_⋆ ≳ 10⁸ M_⊙, the dust-to-stellar mass ratio hardly evolves with redshift, for the same reason that the median metallicity barely evolves. For less massive galaxies (M_⋆ ≲ 10⁸ M_⊙), the dust-to-stellar mass ratio decreases towards lower redshifts, as the gravitational potentials of galaxies of similar halo masses become shallower towards lower redshifts and, together with the gas, dust is more easily ejected by SN-driven winds.

Finally, we also compare our dust and stellar masses to the observational results from the ALMA REBELS survey (see gray points in Fig. 9 from Dayal et al. 2022). Both our IMF models agree with the REBELS data within the observational uncertainties. However, it should be noted that the stellar and dust masses derived from the REBELS galaxies assume a Salpeter IMF.

4.3.3. Metallicity and dust mass relations with UV luminosity

Finally, we explore how the gas-phase metallicities and dust masses of galaxies are related to the galaxies’ observed UV luminosities. Our goal is to assess to which degree these relations differ between our two IMF models, as these relations are mostly independent of the IMF assumed in observational analyses, thereby providing a potential means to distinguish between different IMFs. To this end, we show the corresponding relations at z = 6, 8, 12 and the redshift evolutions of the gas-phase metallicity (in M_UV bins with a width of 1 dex) and dust mass-to-UV luminosity ratio (in dust mass bins of 1 dex) in Figs. 10 and 11, respectively.

Fig. 10. Redshift evolution of the oxygen-based metallicity to observed UV luminosity relation. Left: relation between the oxygen-based gas metallicity and the observed UV luminosity for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of gas metallicities (coloured contours) at a given UV luminosity. Black squares show the observational results from Heintz et al. (2023a) at z ≃ 7.1 − 9.5. Vertical black solid and dashed lines denote the UV luminosity threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the oxygen-based gas metallicity for different UV luminosity bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each UV luminosity bin.

Fig. 11. Redshift evolution of the dust mass to observed UV luminosity relation. Left: relation between the dust mass and the UV luminosity for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of dust masses (coloured contours) at a given UV luminosity. Vertical black solid and dashed lines denote the UV luminosity threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the dust mass-to-UV luminosity ratio for different dust mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each dust mass bin.

Firstly, in contrast to the metallicity and dust mass relations with stellar mass, we find that the metallicities and dust masses of galaxies with a given observed UV luminosity are lower in the evolving IMF than in the Salpeter IMF model. This switch is due to the evolving IMF model’s about 1 dex lower stellar mass-to-light ratio.

Secondly, brighter galaxies (M_UV ≲ −17) exhibit similar median metallicities in both IMF models, with the median value in the evolving IMF model being slightly lower, by less than 0.1 − 0.2 dex, than in the Salpeter IMF model. This similarity results from the comparable halo mass-to-light ratio and halo mass – metallicity relations discussed in Sect. 4.3.1. However, at the same time, similar bright galaxies have lower dust masses, by about 0.4 − 1 dex, in the evolving IMF than in the Salpeter IMF model, a consequence of the lower dust-to-metal mass ratio in the evolving IMF model (see Sect. 4.3.2).

Thirdly, towards lower UV luminosities (M_UV ≳ −17), both the median metallicity and dust mass sharply decrease in both IMF models, albeit at higher UV luminosities in the evolving IMF model (M_UV ≃ −17) than in the Salpeter IMF model (M_UV ≃ −15). The shift of the location of the sharp decrease to higher UV luminosities in the evolving IMF model is the combined result of the model’s higher SN rate and more efficient SN-wind coupling leading to the ejection of gas and dust in more massive halos, as well as the similar halo-mass – observed UV luminosity relationships of the IMF models.

Then, from the right panel in Fig. 10, we see that the median metallicities of galaxies with a given observed UV luminosity increase slightly towards lower redshifts, by about 0.3 dex from z = 12 to 5 for all luminosities. This increase is primarily due to galaxies becoming more dust attenuated with decreasing redshift (cf. intrinsic and observed UV LFs in Fig. 2). It is slightly steeper for the bright and thus more metal-enriched galaxies (M_UV ≲ −20) in the evolving IMF model, since the evolving IMF approaches the Salpeter IMF towards lower redshifts and higher metallicities.

Next, we find similar trends for the median dust mass-to-light ratios of galaxies with fixed dust masses (see right panel in Fig. 11), albeit the evolving IMF model’s increase with decreasing redshift is more prominent. The dust mass-to-light ratio follows essentially the stellar mass-to-light ratio. However, due to the higher rate of dust mass production per stellar mass formed in the evolving IMF model, the relations for the two IMF models are more similar, within about 0.3 dex for galaxies with M_dust ≳ 10⁵ M_⊙.

Finally, we compare our model predictions for the observed UV luminosity – metallicity relation at z = 8 to the observed galaxy sample analysed in Heintz et al. (2023a) (black points in Fig. 10). We note that the inferred UV luminosities of these galaxies are not corrected for magnification; we disregard galaxies with magnifications of > 2.5 (one galaxy) and where the limited wavelength coverage prohibits the inferral of the UV. These observations agree with both IMF models, and given the similarity of our IMF models’ medians and scatter, the UV luminosity – metallicity relation does not provide constraints on the IMF.

4.4. Ionisation topology

Large radio interferometers, such as the Square Kilometre Array (SKA) or Hydrogen Epoch of Reionisation Array (HERA), will measure the 21cm signal from the intergalactic H I gas, providing maps of the evolving spatial distribution of the H I densities and thus of the ionised regions around galaxies during reionisation. A key statistic quantifying the size distribution of ionised regions and thus constraining the ionising emissivity distribution of the galaxy population is the 21cm power spectrum. In this section, we discuss how assuming an evolving IMF affects the topology of the ionised regions around the galaxies and the 21cm power spectrum during reionisation. For this purpose, we show the ionisation maps of both IMF models at z ≃ 9 and z ≃ 7.4 in Fig. 12 and the power spectra of the simulated 21cm signal (see Hutter et al. 2020b, for the derivation of the 21cm signal from the simulated density and ionisation fields) in Fig. 13.

Fig. 12. Neutral hydrogen fraction fields at z ≃ 9.0 (top) and z ≃ 7.0 (bottom) for the evolving IMF (left) and Salpeter IMF (right) models. The volume-averaged value of the neutral fraction in each cell is marked in the panels. For each redshift and model, we show a slice through the centre of the simulation box.

To understand how and why the ionisation topologies differ between our two IMF models, we first analyse the dependence of the number of ionising photons escaping from a galaxy, the escaping ionising emissivity, on its halo mass. As discussed in Sects. 3.2 and 4.1, a more top-heavy IMF leads to a higher abundance of massive stars and lower mass-to-light ratio in the UV; thus, the average intrinsic ionising emissivity per stellar mass ($\dot{Q}/M_{\star }$) is by about a factor ∼7 higher in the evolving IMF than in the Salpeter model. While compensating the higher intrinsic ionising emissivities in the evolving IMF model with a lower $f_{\text{esc}}^0$ value leads to a similar ionisation history compared to the Salpeter IMF model (see Fig. 3), the median relation between the escaping ionising emissivities (Ṅ_ion) and halo mass is flatter in the evolving IMF model, i.e. more low-mass (massive) galaxies have larger (lower) intrinsic ionising emissivity values (see Fig. C.1). This altered distribution of the escaping ionising emissivities stems from (i) the flatter distribution of star formation rates and thus intrinsic ionising emissivities ($\dot{Q}$) for low-mass sources (M_h ≲ 10¹⁰ M_⊙, see Fig. 6) and (ii) the lower f_esc values for massive galaxies (M_h ≳ 10¹⁰ M_⊙). However, the maximum escaping ionising emissivity values for low-mass galaxies (M_h ≲ 10^9.5 M_⊙) are higher for the Salpeter than the evolving IMF model, as the SFR reaches higher values at similar halo masses (cf. Fig. 6 and accounting for the shift of ∼1 dex of the evolving IMF model to lower stellar masses). Together, these trends lead to the following differences in the ionisation topology.

Firstly, the higher Ṅ_ion values around 10¹⁰ M_⊙ in the evolving IMF model lead to a higher abundance of medium-sized ionised regions compared to those in the Salpeter model. While this size difference is barely visible in the ionisation maps (Fig. 12), we can see an excess power around k ≃ 0.3 − 3 h Mpc⁻¹ corresponding to spatial scales of ∼2 − 20 h⁻¹ Mpc in the 21cm power spectra at z ≃ 9 − 7 (Fig. 13). At lower redshifts, the differences in the 21cm power spectra between our two IMF models arise from the models’ diverging ionisation histories, with the evolving IMF showing a more accelerated ionisation of the IGM and, thus, lower 21cm power spectra amplitudes.

Fig. 13. 21cm power spectra at fixed redshifts z using the best-fit parameters noted in Table 1. In each panel, we show results for our two IMF models studied in this work: evolving IMF (dark blue solid line) and Salpeter IMF (red solid line). Each panel shows the average volume-averaged H I fraction in each model at that redshift.

Secondly, the higher maximum Ṅ_ion values for low-mass galaxies in the Salpeter model give rise to more small-sized (≲1 h⁻¹ Mpc) ionised regions than in the evolving IMF model. These regions can be seen in the predominantly neutral regions in the ionisation maps. They also cause the 21cm power spectra of the Salpeter IMF model to slightly exceed that of the evolving IMF model at small scales of k ≲ 6 h Mpc⁻¹ at ⟨χ_HI⟩≳0.5.

In summary, once we adjust for the f_esc values in the evolving IMF model to reproduce similar ionisation histories, we find the ionisation topologies to hardly differ with minor differences occurring on medium scales of 3 − 30 Mpc. Given that the median relation between the escaping ionising emissivity, Ṅ_ion, and halo mass is flatter in the evolving IMF than in the Salpeter IMF model, an evolving IMF results in an ionising emissivity distribution that is stronger correlated to the underlying density distribution, leading to a higher electron optical depth for similar ionisation histories.

5. Conclusions

We have developed the first framework that simulates the evolution of galaxies and the ionisation of the IGM and accounts for an evolving IMF in all relevant physical processes. The IMF considered here parameterises the results found in Chon et al. (2022). It is composed of a Salpeter IMF and a log-flat IMF at the massive end, with the latter increasing towards lower gas metallicities and higher redshifts, reflecting the increasing gas temperature in galaxies due to less efficient cooling and a higher CMB temperature respectively. We have integrated this IMF into our ASTRAEUS framework (Hutter et al. 2021, 2023; Ucci et al. 2023; Trebitsch et al. 2023), where it shapes all IMF-sensitive processes: SN feedback, metal enrichment, and the emission of H I ionising and UV photons. With this updated simulation framework, we ran two simulations: one with a non-evolving, constant Salpeter IMF and the other one assuming our evolving IMF model. The free model parameters were calibrated to match the observed UV luminosity functions at z = 5 − 12.5 and the observational constraints on the global H I fraction, ⟨χ_HI⟩ of the IGM from quasar sightlines, and Lyman-α emission from galaxies and GRBs. We analysed these two IMF scenarios to explore how an evolving IMF alters the properties of galaxies, their evolution through cosmic time, and the topology of the evolving spatial distribution of the ionised regions during reionisation. Our main results are as follows.

1. The evolving IMF scenario hardly changes the redshift evolution of the UV LFs compared to the Salpeter IMF scenario. Any reduction in the mass-to-light ratio due to a more top heavy IMF towards higher redshifts is compensated by stronger SN feedback, which leads to an increased suppression of star formation in addition to the overall assumed lower star formation efficiency.

2. Due to the more top-heavy IMF and thus stronger SN feedback, stellar mass accumulates slower in the evolving IMF than in the Salpeter IMF model; the star formation efficiency is about 2.5× lower. This results in a drop of the stellar-to-halo-mass ratio of about 0.5 − 1 dex throughout reionisation, shifting the stellar mass functions by ∼0.5 − 1 dex to lower stellar masses at z = 6 − 12.

3. As a consequence of the altered stellar-to-halo-mass relation, the evolving IMF model produces (i) UV brighter galaxies when comparing galaxies with the same stellar masses and (ii) together with its stronger and more immediate SN feedback, a broader range of UV luminosities in galaxies with similar stellar masses. Moreover, the decrease of the IMF’s top-heaviness in sufficiently metal-enriched and thus massive galaxies (M_h ≳ 10⁸ M_⊙) with increasing cosmic time causes the evolving IMF’s mass-to-light ratio to approach that of the Salpeter IMF.

4. The median star formation sequence for galaxies where star formation is driven by gas accretion (M_h ≳ 10¹⁰ M_⊙, M_⋆ ≳ 10⁷ M_⊙) does not change as we go from the constant Salpeter IMF to the evolving IMF. This robustness results from the self-similar growth of structures where the fractional mass growth of a DM halo is, on average, independent of its mass. For lower-mass galaxies (M_⋆ ≲ 10⁷ M_⊙), however, SFR rate values depend on the interplay between the efficiency of SN feedback and gravitational potential. Despite having more efficient SN feedback, our evolving IMF model shows fewer low SFR values for these galaxies due to them living in more massive halos with less diverse merger histories.

5. The halo mass to gas-phase metallicity relation remains unchanged as we transition from the constant Salpeter to the evolving IMF. The independence of this relation on the IMF is attributed to calibrating star formation and SN feedback efficiencies to match the observed UV luminosity functions in both models. In contrast, the stellar mass to metallicity relation in the evolving IMF shifts towards higher metallicities compared to the Salpeter IMF, primarily due to the slower accumulation of stellar mass (partially counteracted by the increased oxygen production due to the explosion of more massive stars). Nevertheless, the lower mass-to-light ratio in the evolving IMF model causes its UV luminosity-metallicity relation to shift to metallicity values comparable to or below those of the Salpeter IMF.

6. The dust mass assembly closely follows that of the metal mass. Yet, it is also subject to destruction by SN shock waves, leading to a lower dust-to-metal mass ratio in the evolving IMF model and galaxies where SN feedback does not suppress star formation (M_h ≳ 10¹⁰ M_⊙). This lower ratio leads to a more distinct difference in the UV luminosity-dust mass relation between the two IMF models, with the evolving IMF being shifted towards lower dust masses.

7. The lower mass-to-light ratio in the evolving IMF model results in higher production of ionising photons within galaxies than the Salpeter IMF model. Therefore, when reproducing similar reionisation histories, the average escape fraction f_esc in the evolving IMF is ∼8× lower.

8. For similar reionisation histories, ⟨χ_HI⟩(z), the topology of the ionised regions during reionisation changes only slightly when assuming the evolving IMF compared to the constant Salpeter IMF model. During the first half of reionisation, the evolving IMF model leads to slightly more medium-sized ionised regions (3 − 30 Mpc) due to an enhancement of the escaping ionising emissivity of medium-massive galaxies (M_h ≃ 10¹⁰ M_⊙) when the escape fraction f_esc follows the fraction of gas ejected by SN explosions.

The results from our modelling of an evolving IMF within galaxy evolution have important implications.

Firstly, our findings indicate that an IMF that becomes increasingly top-heavy towards higher redshifts does not increase the abundance of bright galaxies at z ≳ 10 compared to a constant IMF. This lack of discrepancy arises because the lower mass-to-light ratio driven by the greater abundance of massive stars in a top-heavier IMF is counteracted by stronger stellar feedback, which more efficiently suppresses star formation. While we acknowledge the necessity of exploring other plausible IMF parameterisations, such as the one assumed in Steinhardt et al. (2022b), to confirm this self-regulation and ensure the robustness of the star formation main sequence and the metallicity to halo mass relation against different IMF models, we also emphasise that our chosen parameterisation for an evolving IMF is derived from hydrodynamic star-forming cloud simulations. Because of this self-regulation, the current observational constraints on the UV LF at z ≃ 12.5 also remain slightly underpredicted for our evolving IMF model, even when assuming no attenuation of the UV light by dust. This implies that neither a top-heavier IMF nor the ejection of dust by radiative processes (Ferrara et al. 2023) or a combination of the two can explain the early luminous galaxies alone, suggesting that these galaxies need to exhibit higher star formation efficiencies, which could be feedback-free starbursts Dekel et al. (2023), or that our model underpredicts the stochasticity of star formation. Assuming the IMF becomes top-heavier towards higher redshifts, the higher abundance of massive stars could increase stochasticity in star formation due to their shorter lifetimes and stronger stellar feedback, affecting the surrounding ISM more dynamically. Hence, assuming that stars form constantly across a time step – as we do in our model – may not reflect this burstiness sufficiently.

Secondly, the lower average f_esc values in the evolving IMF model, paired with the increased gas-phase metallicity and possibly increased star formation efficiency, suggest that z ≳ 10 galaxies may show strong Lyman-α and [O III] emissions. While these lower f_esc values may seem to be in contradiction with first expectations of a top heavier IMF leading to higher f_esc values due to the increased abundance of massive stars, they may reflect that the formation of fewer stars in the evolving IMF results in higher gas-to-stellar mass ratios in galaxies and fewer low-density tunnels through which ionising photons can escape. Potentially, the increased covering fraction of dense gas clouds and velocities of outflowing gas, resulting from a higher number of SNe exploding at a time, could cause a redshift in the Lyman-α line escaping from these early galaxies that is sufficient for a significant fraction of the Lyman-α radiation to be transmitted through the predominantly neutral IGM. However, to fully understand how the formation of more massive stars at a slower rate changes f_esc and the Lyman-α line emerging from high-redshift galaxies, radiation hydrodynamic simulations and semi-analytic models are needed that follow how stellar feedback of stellar populations with different IMFs affect the structure and gas densities of the interstellar and circumgalactic media.

We note that, in theory, instead of reducing the star formation efficiency f_⋆ in the evolving IMF model to fit the observed UV LFs, we could have assumed the same f_⋆ value as in the Salpeter IMF model, essentially baselining our galaxy model parameters to the z > 10 observed UV LFs, and adjust our dust model (i.e. assumptions on the dust radius and/or grain sizes) to reproduce the UV LFs at lower redshifts. However, with our current model, this option would result in a poorer fit at lower redshifts due to the f_w and f_⋆ values determining the position of the change in the UV LF slope that develops towards lower redshift (e.g. M_UV ≃ −18 at z = 6). The evolving IMF model’s lower mass-to-light ratio would shift this change in slope to higher UV luminosities (M_UV ≃ −20). Yet, galaxies fainter than the UV luminosity where the slope changes are predominantly SN feedback-limited exhibit minimal dust contents. Hence, aligning our model with observations in this scenario would require a dust model that suppresses the UV very efficiently in galaxies with a low dust content but only moderately in galaxies with higher dust contents, which would imply that the composition of dust and/or its concentration in the halo changes rapidly as galaxies grow in mass from M_h ≃ 10^10.5 M_⊙ to 10¹² M_⊙.

Furthermore, our upper mass limit of the IMF of 100 M_⊙ could be adjusted, given the observational evidence revealing the presence of very massive stars with masses exceeding 100 M_⊙ and sub-solar metallicities (see e.g., Martins et al. 2023; Meštrić et al. 2023; Wofford et al. 2023; Upadhyaya et al. 2024). Martins & Palacios (2022) demonstrated that adding the spectra of very massive stars ranging from 100 M_⊙ to 300 M_⊙ and with a metallicity of Z = 0.4 Z_⊙ enhances the ultraviolet continuum of a BPASS 0.1 − 100 M_⊙ spectra by approximately a factor two. Extending the mass range of our IMF models to include very massive stars would likely lead to a further decrease of the mass-to-light ratio, resulting in a lower star formation efficiency, f_⋆, and lower ionising escape fraction normalisation, $f_{\text{esc}}^0$. If we were to extend our current assumption – where each star more massive than 8 M_⊙ explodes and injects an energy of 10⁵¹ erg s⁻¹ – to stars with masses exceeding 100 M_⊙, an IMF including such very massive stars would yield a higher energy injected per stellar mass, compensating for at least some of the decrease in the mass-to-light ratio and driving the aforementioned self-regulation of galactic processes in the presence of an evolving IMF. However, this assumption might oversimplify the stellar evolution process. Stars with masses ≳40 M_⊙ and lower metallicities Z ≲ 10⁻² − 10⁻¹ Z_⊙ are expected to directly collapse into black holes (Heger et al. 2003). While most intermediate and massive galaxies in our simulation attain high enough metallicity values, this assumption may not hold for the lower-mass galaxies, where the energy injected by supernovae would be lower than currently assumed (allowing for a lower supernovae coupling efficiency, f_w). Yet, the evolution of metal-poor massive stars remains highly uncertain owing to uncertainties surrounding their mass loss, rotation, and binary interactions.

Observationally distinguishing between our IMF models poses challenges, given that common galaxy scaling relations often rely on stellar mass, whose inference depends on the assumed IMF. To address this, we need to explore IMF-independent observables or quantities derived from such observables. While the galaxy population’s gas-phase metallicities show minor differences in their UV luminosity dependency between our IMF models, the UV luminosity dependencies of their dust masses differ increasingly towards higher redshifts. This suggests that combining dust mass estimates from ALMA observations with JWST UV measurements of the same galaxies could provide valuable constraints on the IMF during reionisation.

Our implementation of an evolving IMF into the ASTRAEUS framework is a step towards a more-self-consistent model for galaxy evolution and reionisation and enables us now to further explore which observable galactic properties, such as various emission lines or UV slopes, and the relations among them, are best suited to constrain the nature of the IMF of the galaxies in the first billion years of evolution.

Acknowledgments

We thank the anonymous referee for their comments which improved the quality of the paper. We also thank Koki Kakiichi and Viola Geli for useful discussions and comments. AH and CM acknowledge support by the VILLUM FONDEN under grant 37459. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. PD and MT 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. CM also acknowledges support from the Carlsberg Foundation under grant CF22-1322. GY acknowledges Ministerio de Ciencia e Innovación (Spain) for partial financial support under research grant PID2021-122603NB-C21. The VSMDPL simulation has been performed at LRZ Munich within the project pr87yi. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (https://www.lrz.de). The CosmoSim data base (https://www.cosmosim.org) provides access to the simulation and the Rockstar data. The data base is a service by the Leibniz Institute for Astrophysics Potsdam (AIP). This research made use of matplotlib, a Python library for publication quality graphics (Hunter 2007); and the Python library numpy (Oliphant 2006).

References

Appendix A: Luminosities for continuous star formation

Equations. 15 and 14 describe the time evolution of the ionising emissivity and UV luminosity of a starburst for the evolving IMF. However, in ASTRAEUS we assume that stars form continuously with a constant SFR over a time step. Here, we provide the expression we use in ASTRAEUS. The ionising emissivity, $\dot{Q(t)}$, and UV luminosity, L_ν(t), of stars forming with a constant SFR s₀ from gas with metallicity Z between times t_i and t_f at a given time t ≥ t_f > t_i and z being the redshift at time t_i are

$$\begin{array}{cc}\hfill \frac{\dot{Q}(t,t_i,t_f,s_0,z,Z}{{\mathrm{s}}^{-1}{\mathrm{M}}_{\odot }^{-1}}& =\frac{s_0}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}{10}^{e_0}([t_f-\max (t_i,\min [t-{10}^{6.35}\mathrm{yr},t_f])]\hfill \\ \hfill & -\frac{{10}^{-6.35\ast e_1}}{e_1+1}[{(t-\min (t_f,\max [t-{10}^{6.35}\mathrm{yr},t_i]))}^{e_1+1}\hfill \\ \hfill & -{(t-t_i)}^{e_1+1}]),\hfill \end{array}$$(A.1)

with

$$\begin{array}{cc}\hfill e_0& =45.93-1.17(lgZ+2)+0.064z(lgZ+3.46),\hfill \\ \hfill e_1& =-3.91-z[0.051(lgZ+2.97)+0.0239],\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \frac{L_{\nu }(t,z,Z)}{\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}{\mathrm{M}}_{\odot }^{-1}}& =\frac{s_0}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\ast {10}^{e_0}([t_f-\max (t_i,\min [t-{10}^6\mathrm{yr},t_f])]\hfill \\ \hfill & -\frac{{10}^{-6\ast e_1}}{e_1+1}[{(t-\min (t_f,\max [t-{10}^6\mathrm{yr},t_i]))}^{e_1+1}\hfill \\ \hfill & -{(t-\max (t_i,\min [t-{10}^{6.45}\mathrm{yr},t_f]))}^{e_1+1}])\hfill \\ \hfill & -\frac{s_0}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\ast {10}^{e_0+0.45e_1}(\frac{{10}^{-6.45\ast e_2}}{e_2+1}\hfill \\ \hfill & \times [{(t-\min (t_f,\max [t-{10}^{6.45}\mathrm{yr},t_i]))}^{e_2+1}\hfill \\ \hfill & -{(t-\max (t_i,\min [t-{10}^{8.5}\mathrm{yr},t_f]))}^{e_2+1}])\hfill \\ \hfill & -\frac{s_0}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\ast {10}^{e_0+0.45e_1+2.05e_2}(\frac{{10}^{-8.5\ast e_3}}{e_3+1}\hfill \\ \hfill & \times [{(t-\min (t_f,\max [t-{10}^{8.5}\mathrm{yr},t_i]))}^{e_3+1}-{(t-t_i)}^{e_3+1}]),\hfill \end{array}$$(A.2)

with

$$\begin{array}{cc}\hfill e_0& =31.947-0.921(lgZ+1.676)+0.0569z(lgZ+3.498),\hfill \\ \hfill e_1& =0.6876,\hfill \\ \hfill e_2& =-1.56-0.0579z,\hfill \\ \hfill e_3& =-2+0.03z,\hfill \end{array}$$

respectively.

Appendix B: SFR dependence of the stellar mass to metallicity relation

In Fig.B.1, we present the stellar mass-metallicity relations at z = 6, 8, and 12 for both the evolving IMF (left) and the Salpeter IMF (right) models. In contrast to Fig.7, which illustrates the number distribution of galaxies in the stellar–metallicity plane along with the resulting median relations, this figure showcases (i) the median SFR value in each stellar mass–metallicity bin as coloured tiles and (ii) the high-redshift fundamental plane of metallicity (HFPZ) at different SFR values ranging from lg(SFR) = − 6 to 2. We note two differences in the HFPZs of the evolving IMF and Salpeter IMF models: Firstly, the evolving IMF exhibits less scatter in SFR values across the MZR than the Salpeter IMF at all redshifts. This reduced scatter can be attributed to (i) the 2.5× lower star formation efficiency for more massive galaxies and (ii) a higher fraction of the stars forming exploding as SNe in the same time step (and thus reducing the SFR in that time step), as the IMF is more top-heavy, leading to massive stars exploding earlier (see Sect. 4.2 for the detailed explanation). Secondly, the SFR for a fixed metallicity and stellar mass is lower in the evolving IMF than the Salpeter IMF model due to its lower mass-to-light ratio. This lower mass-to-light ratio also results in a weaker correlation between metallicity and SFR.

Fig. B.1. Best-fit high-redshift fundamental planes of metallicity (HFPZ) at z = 6, 8, 12 for the evolving IMF (left) and Salpeter IMF (right) models. Coloured lines show the HFPZ at fixed SFR values, while the tiles depict the median SFR values in the respective stellar mass – metallicity bin according with their colours encoding the SFR value according to the colour bar on the right side of each panel.

Appendix C: Intrinsic and escaping ionising emissivities

In Fig. C.1, we show the distribution of intrinsic ($\dot{\mathrm{Q}}$, top panels) and escaping ionising emissivities (${\dot{N}}_{\mathrm{ion}}=f_{\mathrm{esc}}\dot{Q}$, bottom panels) as a function of halo mass for the evolving IMF (left) and Salpeter IMF (right) models. We note the following differences between the two IMF models: Firstly, due to its lower mass-to-light ratio, the evolving IMF model’s intrinsic ionising emissivity is always higher than the one in the Salpeter IMF model (see Sect. 4.1.2). Secondly, the drop of the ionising escape fraction, f_esc, towards lower stellar masses occurs at higher stellar masses in the evolving IMF than in the Salpeter IMF model. This results from the evolving IMF model’s stronger SN feedback, expelling a larger fraction of gas in more massive galaxies. Thirdly, the evolving IMF model’s relation between the median escaping ionising emissivity and halo mass is flatter (i.e. low-mass galaxies (M_h ≲ 10¹⁰ M_⊙) have larger escaping ionising emissivity values, while more massive galaxies (M_h ≳ 10¹⁰ M_⊙) exhibit lower ones) and has less scatter for low-mass galaxies (M_h ≲ 10¹⁰ M_⊙) than in the Salpeter IMF model. These differences arise from the more immediate nature of SN feedback in the evolving IMF model. The more immediate feedback is attributed to the fact that more stars that explode as SNe are more massive and have shorter lifetimes. Consequently, SN feedback in the evolving IMF regulates star formation efficiency within the same time step. In contrast, the Salpeter IMF model experiences a broader range of star formation efficiencies because more of the energy released by SN explosions originates from stars in previous time steps.

Fig. C.1. Intrinsic (top) and escaping (bottom) ionising emissivities as a function of halo mass Mh for the evolving IMF (left) and Salpeter IMF (right) models at z = 6 (dark blue), 8 (violet), and 12 (orange). Solid lines show the median ionising emissivity and transparent shaded coloured regions the 1 − σ limits of the ionising emissivity values at a given halo mass.

All Tables

All Figures

Fig. 1. Evolving IMF compared with the Salpeter IMF, and the evolution of the massive fraction with redshift and metallicity. Left: evolving IMF at different absolute metallicities Z = 10−2 (purple lines), 10−3 (orange lines), 10−4 (green lines) and redshifts z = 5 (dashed lines), 10 (dot-dashed lines), and 15 (solid lines). For clarity, the evolving IMFs at z = 10 and 15 have been displaced by lg(ϕ) = 2 and 4, respectively. Thin black lines show the Salpeter IMF at the depicted redshifts. Right: massive fraction fmassive depicted as a function of redshift z and absolute metallicity Z. The top-left corner shows the evolving IMF developing into the Salpeter IMF, where fmassive = 0.

In the text

Fig. 2. UV luminosity functions (LFs) at z = 6 − 11, 12.5 and 15 for the evolving IMF (blue lines) and the Salpeter IMF (red lines). Dashed lines show the intrinsic UV LFs, while solid lines depict the dust attenuated UV LFs. Vertical black solid and dashed lines denote the UV luminosity thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the star formation rates have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Black points show the observational constraints from JWST data: Adams et al. (2024) (hexagons), Bouwens et al. (2023a) (diamonds), Bouwens et al. (2023b) (squares), Casey et al. (2024) (triangles), McLeod et al. (2024) (x’s), Donnan et al. (2023) (circles), Finkelstein et al. (2023) (stars), Willott et al. (2024) (pentagons). Light and dark grey points depict the observational constraints before JWST: Atek et al. (2015) (light stars), Atek et al. (2018) (dark stars), Bouwens et al. (2015) (light squares), Bouwens et al. (2016) (dark squares), Bouwens et al. (2017) (light circles), Bouwens et al. (2021) (dark circles), Bowler et al. (2020) (light x’s), Calvi et al. (2016) (dark x’s), Finkelstein et al. (2015) (light diamonds), Ishigaki et al. (2018) (dark diamonds), Livermore (2016) (light ‘downward’ triangles), Livermore et al. (2017) (dark ‘downward’ triangles), McLure et al. (2009) (light ‘upward’ triangles), McLure et al. (2013) (dark ‘upward’ triangles), Schenker et al. (2013) (light ‘left’ triangles), Schmidt et al. (2014) (dark ‘left’ triangles), van der Burg et al. (2010) (light ‘right’ triangles), Willott et al. (2013) (dark ‘right’ triangles).

In the text

Fig. 3. Reionisation history. Top: ratio of the mass- and volume-averaged neutral hydrogen fraction. Bottom: volume averaged neutral hydrogen fraction as a function of redshift. In each panel, we show results for the Salpeter IMF (dash-dotted red line) and the evolving IMF (solid dark blue line) models. Grey points indicate observational constraints from: GRB optical afterglow spectrum analyses (light pentagons; Totani et al. 2006, 2014), quasar sightlines (faint squares; Fan et al. 2006) (medium bright circles; Bosman et al. 2022), Lyman-α LFs (dark upward triangles; Konno et al. 2018; Kashikawa et al. 2011; Ouchi et al. 2010; Ota et al. 2010; Malhotra & Rhoads 2004), Lyman-α emitter clustering (dark downwards triangles; Ouchi et al. 2010), the Lyman-α emitting galaxy fraction (dark leftwards triangles; Pentericci et al. 2011; Schenker et al. 2012; Ono et al. 2012; Treu et al. 2012; Caruana et al. 2012, 2014; Pentericci et al. 2014), Lyman-α equivalent widths (dark rightwards triangles; Mason et al. 2018, 2019; Bolan et al. 2022), dark pixels (light squares; McGreer et al. 2015), and damping wings (light diamonds; Davies et al. 2018; Greig et al. 2019).

In the text

Fig. 4. Redshift evolution of the stellar to halo mass relation. Left: relation between the stellar and halo mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of stellar masses (coloured contours) at a given halo mass. Horizontal black solid and dashed lines denote the stellar mass thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the stellar-to-halo mass ratio for different halo mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each halo mass bin.

In the text

Fig. 5. Redshift evolution of the observed UV luminosity to stellar mass relation. Left: relation between the observed UV luminosity and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of observed UV luminosities (coloured contours) at a given stellar mass. Black solid and dashed lines denote the stellar mass (vertical) and UV luminosity (horizontal) thresholds in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses and observed UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Median values outside the corresponding contours reflect low number statistics. Right: redshift evolution of the observed UV luminosity-to-stellar mass ratio for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each stellar mass bin.

In the text

Fig. 6. Redshift evolution of the SFR to stellar mass relation. Left: relation between SFR and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of SFRs (coloured contours) at a given stellar mass. Vertical black solid and dashed lines denote the stellar mass threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the specific SFR for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The value in the panel shows the central value of each stellar mass bin.

In the text

Fig. 7. Redshift evolution of the oxygen-based gas metallicity to stellar mass relation. Left: relation between the oxygen-based gas metallicity and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of gas-phase metallicities (coloured contours) at a given stellar mass. Black squares and circles show the observational results from Heintz et al. 2023a (z ≃ 7.1 − 9.5) and Curti et al. 2023 (z ≃ 5.8 − 9.4) converted to a Salpeter IMF stellar masses according to Speagle et al. (2014), respectively. Observations at z < 7 are shown in the z = 6 panel and observations at 7 ≤ z < 10 in the z = 8 panel. Vertical black dashed lines denote the stellar mass threshold in the Salpeter IMF simulations, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the oxygen-based gas-phase metallicity for different stellar mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each stellar mass bin.

In the text

Fig. 8. Comparison of the mass-metallicity relations at z = 6 derived from different simulations. Solid blue and red lines show the median values for our evolving IMF and Salpeter IMF models. Black lines show the relations from ASTRAEUS without dust (solid with 1 − σ dispersion; Ucci et al. 2023), ILLUSTRIS-TNG (dot-dot-dashed with 1 − σ dispersion; Torrey et al. 2019), FIRSTLIGHT (dotted; Langan et al. 2020), FIRE (dot-dashed; Ma et al. 2016) and OBELISK (dashed with 1 − σ; Trebitsch et al., in prep.). We have rescaled the metallicities from all these works to our solar value of 8.76.

In the text

Fig. 9. Redshift evolution of the dust to stellar mass relation. Left: relation between the dust and stellar mass for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of dust masses (coloured contours) at a given stellar mass. Black circles show the observational results at z ≃ 7 from the REBELS survey shown in Dayal et al. (2022) assuming a Salpeter IMF. Vertical black dashed lines denote the stellar mass threshold in and Salpeter IMF simulations, below which the stellar masses have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the dust-to-stellar mass ratio for different stellar mass bins with a width of 1 dex. The values in the panel shows the central value of each stellar mass bin. Circles show median values of the REBELS galaxies.

In the text

Fig. 10. Redshift evolution of the oxygen-based metallicity to observed UV luminosity relation. Left: relation between the oxygen-based gas metallicity and the observed UV luminosity for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of gas metallicities (coloured contours) at a given UV luminosity. Black squares show the observational results from Heintz et al. (2023a) at z ≃ 7.1 − 9.5. Vertical black solid and dashed lines denote the UV luminosity threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the oxygen-based gas metallicity for different UV luminosity bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each UV luminosity bin.

In the text

Fig. 11. Redshift evolution of the dust mass to observed UV luminosity relation. Left: relation between the dust mass and the UV luminosity for the Salpeter IMF (red) and evolving IMF (blue). Coloured thick solid lines depict the medians of the distributions of dust masses (coloured contours) at a given UV luminosity. Vertical black solid and dashed lines denote the UV luminosity threshold in the evolving IMF and Salpeter IMF simulations, respectively, below which the UV luminosities have not converged for all galaxies due to the mass resolution limit of the VSMDPL simulation. Right: redshift evolution of the dust mass-to-UV luminosity ratio for different dust mass bins with a width of 1 dex for the Salpeter IMF (dashed lines) and evolving IMF (solid lines) models. The values in the panel show the central value of each dust mass bin.

In the text

	Fig. 12. Neutral hydrogen fraction fields at z ≃ 9.0 (top) and z ≃ 7.0 (bottom) for the evolving IMF (left) and Salpeter IMF (right) models. The volume-averaged value of the neutral fraction in each cell is marked in the panels. For each redshift and model, we show a slice through the centre of the simulation box.
In the text

	Fig. 13. 21cm power spectra at fixed redshifts z using the best-fit parameters noted in Table 1. In each panel, we show results for our two IMF models studied in this work: evolving IMF (dark blue solid line) and Salpeter IMF (red solid line). Each panel shows the average volume-averaged H I fraction in each model at that redshift.
In the text

Fig. B.1. Best-fit high-redshift fundamental planes of metallicity (HFPZ) at z = 6, 8, 12 for the evolving IMF (left) and Salpeter IMF (right) models. Coloured lines show the HFPZ at fixed SFR values, while the tiles depict the median SFR values in the respective stellar mass – metallicity bin according with their colours encoding the SFR value according to the colour bar on the right side of each panel.

In the text

	Fig. C.1. Intrinsic (top) and escaping (bottom) ionising emissivities as a function of halo mass Mh for the evolving IMF (left) and Salpeter IMF (right) models at z = 6 (dark blue), 8 (violet), and 12 (orange). Solid lines show the median ionising emissivity and transparent shaded coloured regions the 1 − σ limits of the ionising emissivity values at a given halo mass.
In the text