Semi-analytic modelling of Pop. III star formation and metallicity evolution - I. Impact on the UV luminosity functions at 𝑧 = 9 − 16

We implemented Population III (Pop. III) star formation in mini-halos within the meraxes semi-analytic galaxy formation and reionisation model, run on top of a N-body simulation with 𝐿 = 10 ℎ − 1 cMpc with 2048 3 particles resolving all dark matter halos down to the mini-halos ( ∼ 10 5 𝑀 ⊙ ). Our modelling includes the chemical evolution of the IGM, with metals released through supernova-driven bubbles that expand according to the Sedov-Taylor model. We found that SN-driven metal bubbles are generally small, with radii typically of 150 ckpc at 𝑧 = 6. Hence, the majority of the first galaxies are likely enriched by their own star formation. However, as reionization progresses, the feedback effects from the UV background become more pronounced, leading to a halt in star formation in low-mass galaxies, after which external chemical enrichment becomes more relevant. We explore the sensitivity of the star formation rate density and stellar mass functions on the unknown values of free parameters. We also discuss the observability of Pop. III dominated systems with JWST, finding that the inclusion of Pop. III galaxies can have a significant effect on the total UV luminosity function at 𝑧 = 12 − 16. Our results support the idea that the excess of bright galaxies detected with JWST might be explained by the presence of bright top-heavy Pop. III dominated galaxies without requiring an increased star formation efficiency.


INTRODUCTION
The first episodes of star formation likely occurred at redshift 30-40 inside low-mass (10 5−7 M ⊙ ) halos with a pristine chemical composition (metal-free or extremely metal-poor).These environments were responsible for the chemical enrichment of the intergalactic medium (IGM) leading to the Pop.III/Pop II transition at redshift 15-10.However, due to the low mass and the high-, the study of Pop.III star formation in mini-halos is challenging both in terms of observations and simulations (see Klessen & Glover 2023 for a recent and comprehensive review).Even with the launch of JWST, we do not have confirmed observations of Pop III stars and only a few potential candidates (Welch et al. 2022;Maiolino et al. 2023).As pointed out by Trussler et al. (2023), a direct detection of a Pop.III galaxies will be extremely challenging as, in the best-case scenario, hundreds of hours of integration time are needed in order to detect an unlensed Pop.III system with 5 confidence.There are a number of Pop.III star formation simulations at different scales.Hydrodynamical simulations that follow the cooling of gas in a single pristine or very metal-poor cloud suggest that the inefficient cooling due to the lack of metals might favour an initial mass function (IMF) that is shifted to larger masses (e.g.Bromm et al. 1999;Hirano et al. 2014;Stacy et al. 2016;Chon et al. ★ E-mail: eventura@student.unimelb.edu.au2021); however, there is no general consensus (e.g.Wollenberg et al. 2020;Jaura et al. 2022;Prole et al. 2022 predict IMF shifted to lower masses).Simulations at larger scales (≥1 cMpc) are instead used to follow the global evolution of Pop.III stars with redshift and their impact on the cosmic metal enrichment and reionization.In the last decade, many of these simulations have been performed both with hydrodynamical (Pallottini et al. 2014;Xu et al. 2016;Sarmento et al. 2018;Jaacks et al. 2018;Skinner & Wise 2020;Schauer et al. 2021;Venditti et al. 2023;Correa Magnus et al. 2023;Yajima et al. 2023) and semi-analytical codes (Trenti et al. 2009;Crosby et al. 2013;Mebane et al. 2018;Visbal et al. 2018Visbal et al. , 2020;;Trinca et al. 2022;Magg et al. 2022;Hegde & Furlanetto 2023;Feathers et al. 2023).In this work, we chose to use the semi-analytical model of galaxy formation Meraxes (Mutch et al. 2016, M16 hereafter) within a simulation that allows us to resolve all the mini-halos down to a few 10 5 M ⊙ in a simulation of L = 10 Mpc ℎ −1 .The size of this simulation is not large enough to be representative of the Universe (we will miss the most massive galaxies).However, the scales are large enough to investigate the impact of the main physical processes on the Pop.III star formation.We incorporated a number of new physical processes that are relevant to Pop.III star formation in mini-halos including: (i) molecular hydrogen cooling functions for mini-halos, (ii) baryon-dark matter streaming velocities, (iii) photo-dissociation of H 2 molecules from the Lyman-Werner background and (iv) the chemical evolution of the intergalactic medium.This latter effect has been implemented assuming that metals are released through supernova explosions and within a bubble that expands accordingly to the Sedov-Taylor model (the approach is very similar to the analytical calculation shown in Furlanetto & Loeb 2003).We found that these bubbles are generally small and that they roughly agree with the previous estimate by Trenti et al. (2009), with typical radii of 150 ckpc at  = 6.The implementation of both internal and external metal enrichment allows us to understand whether a galaxy will form Pop. III or Pop.II stars and thus quantify the impact of Pop.III galaxies on the total luminosity function at  ≥ 5.This paper is structured as follows: In Section 2, we present the new physics implemented in Meraxes.In Section 3, we show the main global properties of the Pop.III star formation across cosmic history discussing which parameters have a stronger impact on the global star formation history and on the metallicity of the IGM.In Section 4 we discuss the SED evolution of Pop.III galaxies and the UV luminosity functions at  = 9 − 16 for different Pop.III IMFs.Finally, we discuss our main results and conclusions in Section 5.

POP. III STAR FORMATION IN MERAXES
The semi-analytical model meraxes was first developed by Mutch et al. (2016) (hereafter M16), Qin et al. (2017); Qiu et al. (2019) in order to study galaxy formation and growth through the Epoch of Reionization.meraxes assumes that all the galaxies form in a previously chemically enriched Universe and inside atomic cooling halos.Such an approximation does not allow us to study the first episodes of star formation that mostly occurred in mini-halos when the Universe did not have any metals.The version of meraxes presented in this work allows us to compute the physics of the first episodes of star formation from the initial molecular cooling of the gas to the external metal enrichment from the supernova feedback.As in the previous work, meraxes is coupled to the reionization so that all the radiative backgrounds are computed in a self-consistent way from the galaxy properties.This has been done by implementing a modified version of 21cmfast (Mesinger et al. 2011) that includes the local ionizing UV background from Sobacchi & Mesinger (2014) and the X-ray heating (Balu et al. 2023).In this version, we also included the Lyman-Werner background.

High-Resolution N-body simulation
The updates to meraxes for Pop.III necessitate high mass and spatial resolution.For this purpose, we introduce L10_N2048 (hereafter L10) from the Genesis suite of dark matter only N-body simulations.L10 is a periodic cubical simulation of side 10 ℎ −1 Mpc and consists of 2048 3 dark matter particles of mass m p = 9.935 ×10 3 ℎ −1  ⊙ resulting in a halo mass resolution of ∼ 3.18 ×10 5 ℎ −1  ⊙ (based on a minimum of 32 particles per halo).The simulation, run using the SWIFT (Schaller et al. 2018) cosmological code, evolves these dark matter particles from  = 99 down to  = 5.The halos were identified using the friends-of-friends phase space halo-finder VELOCIraptor (Elahi et al. 2019a) and the merger trees were constructed using TREEFROG (Elahi et al. 2019b).We note that the mass resolution achieved in this simulation allows us to resolve all the mini-halos below  ∼ 30 and is the highest resolution on which meraxes has been run.The output trees of the N-body simulation used in this work are available on Zenodo at Balu et al. (2024).

Molecular cooling
The first process in order to enable star formation is the cooling of the gas.For dark matter halos with  vir ≥ 10 4 K, the main coolant is the atomic hydrogen, while in mini-halos (10 3 K ≤  vir < 10 4 K), the cooling occurs via roto-vibrational transitions of molecular hydrogen (e.g.Tegmark et al. 1997).For the details on how the cooling of the gas is implemented in meraxes, we refer the reader to M16; in this section, we only highlight the main differences between the atomic and the molecular cooling.As in M16, we compute the ratio of the specific thermal energy to the cooling rate per unit volume: where  hot () is the gas density profile (we assume a singular isothermal sphere), k is the Boltzmann constant, and T is the temperature of the gas (which we set to be the virial temperature of the halo).M16 included only atomic cooling halos.Thus they set the mean molecular weight μ = 0.59 assuming a fully ionized gas and the cooling function Λ() in units of erg cm −3 s −1 from Sutherland & Dopita (1993).For mini-halos we instead set μ = 1.22 (fully neutral gas) and implement the molecular hydrogen cooling functions Λ( H , ) from Galli & Palla (1998).This choice is valid as long as we assume that the cooling inside mini-halos occurs only due to H 2 .This may not be valid as, if a mini-halo is chemically enriched with metals by a nearby halo, metals are much more effective in the cooling of the gas (e.g.Nebrin et al. 2023).However, as we will show in section 3, this enrichment is almost ineffective at  > 10.The molecular cooling function Λ depends on the gas density of the halo (which can be directly computed for each galaxy assuming an isothermal sphere) and on the molecular hydrogen fraction   2 .For the latter we assumed a fixed value of 0.1%, which is consistent with results from Nebrin et al. (2023) for halos of T vir ≃ 5 × 10 3 K at  = 15 − 20 1 .As in M16, the cooling time is used to find a cooling radius that determines the amount of gas that is cooled down.

Streaming velocities
As a first approximation, gas in a mini-halo can start to cool down the gas once it reaches the virial temperature of ≃ 10 3 K.From Barkana & Loeb (2001), adopting  = 1.22 for a fully neutral gas, this requirement would correspond to a minimum virial mass of (Visbal et al. 2015): However, there are a number of effects that can decrease the amount of molecular hydrogen present in the halo, reducing the cooling efficiency and ultimately increasing the minimum virial mass of a halo capable of cooling (see Nebrin et al. 2023 for an extensive discussion and comparison between results found in different simulations.)A non-radiative process that can delay the gas cooling in very lowmass halos is the streaming velocity between baryons and dark matter (Tseliakhovich & Hirata 2010).This effect is a consequence of the different decoupling of baryons and dark matter particles from photons that results in a velocity difference between the two species of  bc that can be expressed as a constant multiple of the root-meansquare (rms) streaming velocity  bc = 30 × (1 + )/1000 km s −1 .The presence of a relative motion between baryons and dark matter particles will make it harder for baryons to fall into the potential wells of dark matter halos, delaying the accretion and, hence, the cooling of the gas in mini-halos.The main outcome of this physical process is to delay the very first episodes of star formation around  ∼ 40.
Throughout this work, we implemented the effect of the streaming velocities as per other semi-analytical models based on a fitting function found by Fialkov et al. (2012), which is calibrated to reproduce the results of the hydrodynamical simulations of Greif et al. (2011) and Stacy et al. (2011): where a = 3.714 km s −1 and b = 4.015 km s −1 .Equation 3 provides the minimum circular velocity that a halo needs to have in order to have enough H 2 to cool down the gas.V cool,H 2 can be easily converted into a virial mass M cool,H 2 using the formula in Barkana & Loeb (2001).In this work we fixed  bc () =  bc () and we assumed this value for the entire box. 2

LW background
A self-consistent model of star formation in mini-halos must consider the photo-dissociation of H 2 from UV photons in the Lyman-Werner (LW) band ([11.2 -13.6] eV).LW photons destroy H 2 and thus prevent mini-halos gas from cooling (Haiman & Bryan 2006).We implement this effect by changing the minimum mass for molecular cooling using the fitting from Visbal et al. (2015): This critical mass only applies to minihalos which are below the atomic cooling halo mass threshold.In the absence of an LW background, this equation simply gives M crit,MC as defined in Section 2.3.J LW is the LW flux in units of 10 −21 erg s −1 cm −2 Hz −1 sr −2 that reaches the minihalo.LW photons have a mean free path of ∼ 100Mpc, so each minihalo will be affected even by distant galaxies that formed at higher redshift.Thus, we modelled the LW background by integrating contributions across the cosmic history (see a similar approach in Qin et al. 2020 which is briefly summarized below.)As for all the radiative backgrounds in meraxes, we follow an excursion-set formalism (Furlanetto et al. 2004) which counts the number of photons in a certain band in spheres of radius R centred at location and redshift (x,z).We decrease the radius down to the size of the cell R cell .At each of these locations, we compute the LW emissivity  LW (x,  ′ ) using the spectral energy distributions from Barkana & Loeb (2005).These give   : the number of photons per solar mass per unit frequency for both Pop III and Pop II stars.We assume that LW photons are absorbed only at resonant frequencies (thus we are neglecting any absorption from H 2 in the IGM) where  n =  LL (1 −  −2 ) and  LL = 3.29 × 10 6 GHz is the Lyman limit frequency.Under these approximations, we can compute the LW emissivity smoothed over R at redshift  and location  for sources 2 This assumption is not entirely accurate as the streaming velocities are roughly constant up to a scale of a few Mpc while our box is larger than 10 Mpc.However, the effect of streaming velocities is dominant only at  ≥ 30 before the first stars form, and the Lyman-Werner background is built up.
emitted at redshift  ′ as: ′′ ℎ ′′ (5) In the equation above, we directly link the star formation rate density (SFRD) for both Pop.III and Pop.II stars (which we can compute from meraxes) to the LW emissivity summed over the Lyman series.This sum accounts for the resonances in the Lyman series (Lymann photons will be absorbed and re-emitted as a Lyman− photon).
Given the large mean free path of LW photons, we must also account for distant galaxies at higher redshift  ′ >  with the redshifted spectrum  ′ =  1+ ′ 1+ .The sum of these two effects causes the peculiar shape of the LW spectrum (see Fig. 2 in Qin et al. 2020).Following Barkana & Loeb (2005) we account for all the Lyman resonances with  ≤ 23.Under these approximations we can find the maximum redshift  max from which a LW photon can reach : Finally, we need to convert the emissivity into a flux following (Qin et al. 2020): The approach described above differs from the one in 21cmfast only for the computation of the SFRD.In Qin et al. (2020), the SFRD is estimated from the density field and the collapsed fraction, while in this work, it is computed directly from meraxes, which tracks the formation of each galaxy and its entire star formation history. 3 .
Together with the LW background, the main radiative feedback that suppresses star formation is the UV photo-ionization (the ionizing UV background inhibits gas accretion in low-mass halos Sobacchi & Mesinger 2013).We used the same prescriptions as in M16, which consists of reducing the baryon content through a baryon fraction modifier f mod (that depends on the local ionizing rate and the time when the nearby IGM becomes ionized) that stops the gas infall.In difference from the LW feedback, the UV photo-ionization is relevant only during reionization ( ≤ 10), and it also affects atomic cooling halos as massive as ∼ 10 9.5  ⊙ .

Metal Evolution
Once stars start to form in the Universe, they also explode as supernovae, releasing metals.Most of these will stay inside the same galaxy, contributing to the chemical enrichment of the galaxy itself.This process is commonly known as "genetic" enrichment and is likely to be the main mechanism of chemical enrichment of the Universe.However, some of the metals escape their parent galaxy.In this case, they will pollute the nearby IGM, and if later a galaxy forms in a region where the IGM was enriched, the new galaxy will be pre-enriched with metals.This latter mechanism is referred to as "external" metal enrichment (e.g.Pallottini et al. 2014;Smith et al. 2015;Hartwig et al. 2018;Visbal et al. 2020;Yamaguchi et al. 2023).Keeping track of the metallicity evolution of the Universe is crucial in order to put constraints on when the Pop.III/II transition occurred.We account for both processes, and we are able to follow the metallicity evolution of the IGM.Firstly, we choose a critical metallicity  crit = 10 −4  ⊙ as the threshold value below which a galaxy will form Pop. III stars.All new galaxies, unless externally polluted, will accrete pristine gas (without any metals) onto the hot gas reservoir.This gas, once it cools, will provide the reservoir for the star formation.Hence, a galaxy that is not externally polluted will always form Pop. III stars for the first time.At each snapshot, we compute the metallicity of the cold gas reservoir from the amount of metals released by earlier supernovae and if this is higher than  crit , the galaxy will form Pop. II stars, otherwise it will form Pop. III.While internal enrichment is the main mechanism that drives the Pop.III/II transition, some galaxies can be externally enriched through supernova winds originating in a nearby galaxy.Once several supernovae in a galaxy go off, they will form a "super-bubble" that will expand outside the galaxy escaping the binding energy of the dark matter halo.We followed the expansion in time of this "metal bubble" using the analytic approximation in Dĳkstra et al. (2014): All quantities that appear in this equation are computed in meraxes.
Δ SN is the total supernova energy released at a certain snapshot, and it is computed using the stellar mass that goes supernova in that snapshot multiplied by the supernova explosion energy (see Eq. 11 in M16).n gas is the number density of the gas to which the bubble has expanded (while the bubble is smaller than the virial radius, this is the gas density of the galaxy; otherwise, the density of the IGM), and t is the time since the explosion occurred.We assume that all the supernova events will occur in the middle of the snapshot 4 .Note that since meraxes accounts for both contemporaneous and delayed supernova feedback (see more in M16), each galaxy has several bubbles associated with the same star formation episode.However, we consider only the largest of these bubbles so that each galaxy has only one associated bubble.Having calculated the bubble size, we can predict if a nearby galaxy will accrete pristine or enriched gas.In order to reduce the computational cost, we avoid computing the distance between all the pairs of galaxies and instead use the grid-based approach outlined below: • Build a high-resolution grid with  3 grid cells smoothing the matter density grid from the N-body simulation.The size of the cell is chosen so that the volume of the cell is similar to the volume of the typical bubble at the end of the simulation ( ave,bubble ≃ 100ckpc, see Fig. 2).For our box, this requirement is satisfied when  grid = 128.This ensures that the metal bubble will not overflow outside the pixel where it originates.
• For each cell i, we compute the average metallicity  IGM,i as the ratio between the sum of the metals ejected by all galaxies inside the cell and the total gas in the cell: Here  IGM,i =  crit Ω  (1+  )    3 is the total baryonic mass inside each cell which depends on the average matter overdensity of the cell   and the cell size . EjGas and  gas are respectively the total mass of the ejected gas and the bounded gas of each galaxy.We only let the galaxies with  bubble ≥ 3 ×  vir contribute to  MetalsEj,j and  EjGas,j .For the computation of the metals ejected from Pop. II galaxies, we refer the reader to M16 and Qiu et al. ( 2019) and for Pop.III galaxies to Section 2.6.
• For each cell i, we compute the volume fraction filled with metals (the metal filling fraction), summing the volume of the largest bubble surrounding each galaxy within the pixel and dividing by the cell volume.As above, we account only for those galaxies with  bubble ≥ 3 ×  vir .This choice reflects the fact that a central galaxy needs to be surrounded by a bubble at least as large as its virial radius in order to pollute a nearby galaxy.The factor of 3 has been chosen given that the distance between two central galaxies in our simulation is always larger than three times the virial radius of the biggest galaxy.
Hence, the metal filling factor for the cell i is  metal,i =    bubble,j   .
• At the beginning of each snapshot, we assign the probability p for external metal enrichment to each galaxy inside the cell i.This probability will be given by the metal filling factor   .This choice assumes that the galaxies are randomly distributed within each pixel.This assumption holds only if we use a high-resolution grid where the pixel size is small enough such that clustering effects are negligible.
• We assign a random number  between 0 and 1 to each newly formed galaxy, and at each snapshot, we compare this random number to .When the condition  ≤  is satisfied, we label that galaxy as externally enriched, and it accretes enriched gas with a metallicity equal to the average metallicity of the cell that belongs to ( IGM,i ).Furthermore, when a galaxy experiences a star formation episode, we enforce  = 1 (in this latter case, we know that this galaxy will be inside its own metal bubble and thus cannot accrete pristine gas).With this latter condition and internal enrichment, we effectively stop Pop.III star formation inside a galaxy after the first supernova episode.
We illustrate an example in Fig. 1, where we selected a 2D slice of our L10 box at z = 10 with a thickness of 0.078ℎ −1 cMpc.The red (blue) dots are the galaxies with  vir ≥ 10 6  ⊙ (not) externally metal-enriched, while the green dots are previously formed galaxies with r bubble ≥ 0.05ℎ −1 cMpc.For a better visualization, we zoomed into a slice with 0.24ℎ −1 cMpc on a side showing the actual size of the metal bubble (right plot).Inside this region, all the 7 galaxies that formed inside the bubble are actually marked as externally enriched.At the same time, the 5 galaxies that are marked as not enriched lie outside the bubble, showing that the adopted approximation works quite well (see Appendix A for a more detailed discussion).This demonstrates that this method can reproduce the statistical properties of inhomogeneous metal enrichment.The main limitation of this technique is that we are not accounting for the overlap of the bubbles and thus are overestimating the metal-filling factor and underestimating the metallicity of those galaxies that are polluted from more than one galaxy.However, given the small size of the bubbles (see discussion in Section 3 and Fig. 2), this is not a major factor, especially prior to reionization.We note that the methodology shown here is similar to Sassano et al. (2021) with the main differences being that we divided the volume of our simulation in cells allowing a more precise computation of the metallicity of the enriched regions.Fig. 2 shows the bubble size distribution function (r bubble in ckpc) at redshift 20, 15, 10 and 5. From this plot, we see the bubble growth over time from a typical radius of ∼ 30 ckpc at  = 20 to ∼ 150 ckpc at  = 5.These values are consistent with the results shown in Trenti et al. ( 2009) (r bubble ≤ 150ℎ −1 ckpc at  = 6) with only a few bubbles at  = 5 larger than this value.Fig. 3 shows the IGM metallicity (top row), metal filling factor (middle row) and density (bottom row) maps at four different snapshots.
The red (blue) cells in the top row have  IGM ≥ (<) crit indicating that galaxies that will be externally polluted in those cells will accrete enriched (pristine) gas and will form Pop. II (III) stars.Comparing the maps at the same snapshot, we see that the more overdense regions have a larger filling factor and a higher metallicity as these hosts more galaxies, hence more metal bubbles.Fig. 3 also shows the progressive enrichment of the IGM over the cosmic history since, as we move to lower redshift (from left to right panels) we have more enriched cells and with larger filling factors.

Pop. III star formation
In this work we modify Pop.II description in M16 to allow star formation in Pop.III galaxies (see Section 2.4 in M16).In particular, once a galaxy reaches a mass of cold gas larger than a critical value  crit , this galaxy will convert the cold gas into stars according to the star formation efficiency  SF .The value of  crit is computed as in M16 and is determined by a critical surface density Σ norm of the gas in the cold disk.Since we are now considering two different stellar populations, we adopted two different free parameters, both for the star formation efficiency  SF,III and for the critical surface density Σ crit,III .Recent full hydrodynamical simulations that follow the collapse of a pristine gas cloud until a Pop.III star is formed (e.g.Hirano et al. 2014;Stacy et al. 2016;Chon et al. 2021), suggest that the Pop.III IMF is shifted to larger masses as a result of less fragmentation due to inefficient cooling.Within this work, we adopted the IMFs from Raiter et al. (2010) (see Table 2), while for Pop.II stars, we assumed a Kroupa (2001) IMF.For the fiducial model, we chose a Salpeter IMF with a mass between 1 and 500 M ⊙ .The choice of the IMF is crucial as it determines many properties of the stellar population, including what fraction of stars will explode as supernovae and, hence, the amount of energy and metals injected into the IGM.In this work, we explore the impact of the 4 free parameters discussed above (see also Table 1) since these have the strongest impact on the Pop.III star formation history.meraxes includes more free parameters describing the energy coupling efficiency of the supernova explosion (see M16 and Qiu et al. 2019) that will not be explored in this work (we will take the same fiducial values as for Pop.II).Since the focus of this work is on the first galaxies formed during the Cosmic Dawn, we do not explore parameters describing the UV ionizing and X-ray radiation (e.g.escape fraction), and we will take the same fiducial values as in the previous works (see Balu et al. 2023).Those will be discussed in a future work focused on the EoR using a larger box.The fate of a zero-metallicity star is quite uncertain due to the many poorly constrained physical mechanisms.Here, we adopt a simplified picture of the final fate of a Pop.III star depends only on its initial mass (Heger & Woosley 2002).For masses below 8 M ⊙ , there will be no SN event.If M ★ ∈ [8, 40]M ⊙ it will explode as a corecollapse SN (CCSN) leaving a remnant, if M ★ ∈ [140, 260]M ⊙ there will be a pair-instability SN (PISN) leaving no remnant and if M ★ ∈ [40, 140]M ⊙ or M ★ >260M ⊙ stars collapse directly into a black hole (BH) with negligible feedback.Given the different endings of a Pop.III star life, the choice of the IMF is crucial to compute the timing of delayed supernova feedback.While we use the technique of Qiu et al. (2019) for Pop.II stars with precomputed tables assuming a Kroupa IMF (Kroupa 2001), for Pop.III supernova feedback we estimated the amount of SN energy and metal yields with an analytic calculation as in M16.The total energy provided by Pop.
III supernovae explosions at the snapshot j is: where Δ i ★ is the new stellar mass formed at an earlier snapshot i,  energy is the energy coupling efficiency (Qiu et al. 2019),  CC = 10 51 erg and  PISN = 10 53 erg are the typical energy for a core collapse and a pair-instability supernova respectively and are respectively the fraction of stars formed during snapshot i that go core collapse and pair-instability supernova in snapshot j.These are computed by integrating the chosen IMF over the correct mass limits.Assuming that all-star formation occurred in the middle of the snapshot,   min,j and   max,j for CCSN are computed from the lifetimes for Pop.III stars using Schaerer (2002) assuming no mass loss and zero metallicity.We highlight that a zero metallicity star with  ∈ [140, 260] has a lifetime smaller than the time separation between two consecutive snapshots in meraxes; hence, it will explode as a PISN in the same snapshot in which it forms.For this reason, in Eq. 12, the mass limit for  PISN is the entire mass range of the PISN and when computing the total energy from PISN injected at snapshot j in Eq. 10 we consider only the stars that are formed at snapshot j.For CCSN, instead, we keep track of the star formation history over the last 17 snapshots, which correspond to ≳ 40 Myr, after which all stars with  ≥ 8 ⊙ will already be exploded.To compute the amount of gas that is reheated and ejected from the halo, we adopted the same prescription as Qiu et al. (2019).We also update the amount of ejected gas and metals from Pop. III stars.These are taken from Heger & Woosley (2010) assuming nonmixing (no rotation) and a supernova explosion of 1.2 ×10 51 erg.As in M16, when a star goes core-collapse supernova at snapshot j, it will release an amount Δ CCSN,Z,j of metals and Δ CCSN,X,j of gas into the interstellar medium: where  CCSN is the total fraction of stars that ends as a CCSN: and Δ CCSN,Z;i,j is the amount of metals released by a star formed at snapshot  and going supernova at snapshot j: In the above equation, Z is the metal mass fraction released into the ISM from supernovae, and it is incorporated inside the integral because it is not fixed over the mass range (i.e in stars with  ★ ≥ 25 ⊙ most of the metals will fall onto the remnant due to the massive fallback).With an entirely analogous calculation, it is possible to find Δ CCSN,X,j (we just need to replace Z with the gas mass fraction X).Once again, we are doing delayed SN recycling only for CCSN because, in the case of pair instability, all the feedback is contemporaneous.For PISN, we considered the gas and metal yields from Heger & Woosley (2002).In all the previous meraxes works, all the stellar mass locked up in remnants (neutron stars and BHs) was neglected as it was recycled into the gas mass budget of the galaxy.We decided to drop this approximation for Pop III stars as they are likely to leave more massive remnants (some might be massive enough to be the first "light" seeds of supermassive BHs).Firstly, we need to consider the BHs that formed after a "failed SN scenario" typical of a star with an initial mass  ★ ∈ [40, 140]  ⊙ and  ★ > 260 ⊙ .The BH mass arising from this scenario will be: Finally, we need to account for the BH remnants that are left after a CCSN.In this case, we use equations 13 and 15 substituting the metal mass fraction with the BH mass fraction  = (1 −  − ) (since all the remaining material that has not been ejected is locked onto the stellar remnant).meraxes, does not evolve the remnants (via accretion), as the main focus of this work is Pop.III stars.However, the impact of the first accreting BHs on the formation of the supermassive black holes and radiative backgrounds might be important even at high-redshift (Ventura et al. 2023).

Pop. III luminosity
In order to investigate the observability of the first Pop.III galaxies we implement spectral energy distributions (SED) for Pop.III stars.We use SEDs from Raiter et al. (2010) that have been computed for the IMFs listed in Table 2 assuming that star formation occurs in an instantaneous burst (see instantaneous burst in Raiter et al. 2010).These SEDs also include the nebular continuum emission and the UV ionizing properties.We used the latter to compute the number of UV ionizing photons per baryon    for each IMF provided.To compute the luminosity of galaxies at a specific wavelength , we used the model introduced by Qiu et al. (2019) (see also Mutch et al. 2023), and we extended it to Pop.III galaxies.Having the Pop.III SEDs, we can compute the luminosity of a Pop.III galaxy at time t as: Here reason, we expanded the calculation of the luminosity of Pop.III galaxies assuming that new stars Δ ★ form instantaneously at the time  −  =  0 .Hence, we can write ( − ) = Δ ★ ( 0 ) and Eq. 17 simplifies to: As we will discuss in Section 4, instantaneous (instead of continuous) Pop.III star formation has an impact on the estimated luminosity function of galaxies hosting Pop.III stars.This is because Pop.III stars have lifetimes that can be shorter than the duration of the snapshot (at  = 10 the snapshot will last 10 Myr, see Fig. 4).Hence, when we consider continuous star formation, we average the star formation over the entire snapshot, leading to lower luminosity with a higher duty cycle.
In order to account for instantaneous star formation, we assumed that a Pop.III star formation episode in a galaxy can occur at random Δ ∈ [0,  snap ] prior to the end of the snapshot.For the detailed evaluation of the UV luminosity function accounting for the stochasticity in the time at which the burst of star formation occurs in different galaxies, we refer the reader to Appendix B.

GLOBAL EVOLUTION OF POP. III STARS
To explore the Pop.III contributions to the cosmic star formation history, we ran two simulations on the L10 box, one with all the updates described in Section 2 adopting the fiducial parameters (see Table 1 and one without the new physics that we labelled as "NoMini").Given that there are no observational constraints on Pop.III, the choice of a fiducial model is arbitrary.The one adopted in this work has both a low star formation efficiency and a Salpeter-like IMF, which will result in a relatively small global impact of Pop.III star formation compared to Pop.II.In Figure 5 we show the stellar mass function (SMF) at different redshifts accounting for all (black), only Pop.III dominated (cyan), and only Pop.II dominated (red) galaxies.We classified a galaxy as Pop.III or Pop.II dominated based on which population is brighter in the UV band.Thin (thick) lines are computed from the fiducial (NoMini) model.The total SMF has two components, separated by their stellar types.In this fiducial model, as Pop.III SFE is 10 times lower than Pop.II, we see a resulting bimodality (which we will further expand discussion in Sec 3.2.1).As shown in the upper panels, the new updates on meraxes mostly affect the lower end of the SMF with a larger impact at high-.This reflects the star formation in mini-halos (we are now considering molecular cooling), which is dominant at  ≥ 15 before the Lyman-Werner background becomes strong enough to photodissociate all the molecular hydrogen.The low mass of Pop.III systems is mainly a consequence of both the shorter lifetimes of Pop.III stars (given their larger mass) and the lower SF efficiency, but also suggests that most of the Pop.III star formation must occur in mini-halos.In order to investigate this, we checked which halos form Pop. III and Pop.II stars.In Fig. 6, we show the halo mass function for Pop.III (cyan line) and Pop.II (red line) star-forming halos at  = 20, 15, 10 and 5.The dashed line corresponds to the limit of atomic cooling ( vir = 10 4 K).Given the small size of the box, we have only a few atomic cooling halos at  ∼ 20, so all star-forming systems are mini-halos (and are mostly Pop.III).At lower redshift, the peak of the distribution shifts to higher masses (by  = 10 most of the stars form in atomic cooling halos except for some low mass halos that are undergoing a merger event) and the impact of Pop.II increases.Finally, Fig. 6 shows that Pop.III stars mostly form in mini-halos (except at  = 5 when the radiative feedback has suppressed both the gas infall and cooling in mini-halos).This is because the more massive halos are more likely to have already experienced star formation and so will be internally chemically enriched.
In Figure 7 we show the total star formation rate density (SFRD) as a function of redshift for both the fiducial and the NoMini model.In the case of the fiducial model, we also show the Pop.III (Pop.II) contribution in the upper (lower) panel.Mini-halos have an appreciable impact on the SFRD only up to  ≥ 20, while at lower redshift, most star formation occurs in atomic cooling halos  The early flattening occurs due to the build-up of the Lyman-Werner background that affects Pop.III star formation in minihalos.The sharp drop at  ≤ 10 is mostly caused by the photoionizing feedback from reionization that also affects the atomic cooling halos.The feedback from reionization also mildly affects Pop.II as can be seen from the flattening of the red line at  ∼ 8 in the upper panel 5

Chemical enrichment of the IGM
Compared to reionization, the chemical enrichment of the IGM is a much slower process due to the lower velocity of the expansion of the metal bubbles.As already found in previous works (e.g.Visbal et al. 2020;Yamaguchi et al. 2023), the impact of external metal enrichment is subdominant compared to internal metal enrichment.However, external enrichment might still be important for low-mass halos that do not form stars until  ∼ 10.The top left panel of Fig. 8 shows the redshift evolution of the average metallicity of the box (in solar units) computed using Eq. 9 and averaging through all 128 3 pixels.As expected, the average metallicity increases monotonically.It crosses the critical value at  ∼ 11 and at the end of the simulation is ∼ 5 × 10 −3  ⊙ .The redshift evolution of  IGM is in very good agreement with Yamaguchi et al. (2023), especially for their "bursty" models, which is the one that 5 Given the small size of the box, the feedback from reionization it is slightly overestimated as we (i) are missing the most massive systems that would not be affected by the photo-ionizing feedback and (ii) the box gets completely ionized at  ∼ 6.9 which is much earlier than what suggested by the latest Lyman- forest measurements (e.g.Qin et al. 2021;Bosman et al. 2022).
most closely resembles the star formation in meraxes.However, the average metallicity of the IGM does not completely inform us of the average metallicity of the galaxies since, when averaging through the entire box, we are considering voids that have no galaxies and hence zero IGM metallicity.For this reason, in the top right panel, we computed the black solid line, averaging only through cells that have at least one metal bubble (meaning that there must be at least one galaxy that formed stars).The shaded region shows the scatter in the metallicity of these cells from the highest to the lowest metallicity.The average metallicity, in this case, is much higher, and it is always above the critical value.This suggests that once the first galaxies form, the ejection of metals from Pop. III is quite effective in the nearby regions, allowing a fast Pop.III/II transition in those cells where star formation already occurred (and the galaxies are internally enriched).The lower panels instead show the fraction of the IGM that reached the critical metallicity (or equivalently, the metal filling factor ( ≥  crit )).Considering the entire box (left panel), less than 1% of the volume gets enriched above  crit by  ∼ 5.This result is fairly consistent with Visbal et al. (2020) and Yamaguchi et al. (2023) (∼ 1% by  = 6).These small differences are likely due to differences in the star formation models and in the choice of parameters, such as the star formation efficiency.If we focus only on the regions with at least one galaxy that formed stars, we find larger filling factors, and by the end of the simulation, those pixels are all completely enriched.This final result reflects our choice of pixel size that is designed to be similar to the average volume of the metal bubble at  ∼ 5.
Finally, in Fig. 9, we show the halo mass function for the externally (black line) and internally (grey line) metal enriched halos at  = 20, 15, 10 and 56 .While most halos get internally enriched by their own star formation, at  = 5 halos with virial masses  vir ≤ 10 7.5  ⊙ get their metals mostly from a nearby supernova bubble.This picture reflects the fact that low-mass halos during reionization are not able to form stars because of the LW and photoionizing feedback.Hence, those objects can be chemically enriched only from an external source.In conclusion, when looking at the global evolution of star formation, the effect of the external metal enrichment is quite negligible as it is important only for low-mass halos at low redshift that will hardly form stars due to the radiative feedback effects.We verified this by running a simulation without external metal enrichment (a halo can only get enriched by its own star formation).We plot in Fig. 10 the total, Pop.III and Pop.II SFRD (dashed lines) compared to the fiducial model (solid lines).There are no appreciable differences between the two models except at  ≤ 8 when the dashed line is slightly larger (about 0.1 dex difference).

Pop. III star formation parameters
The results in previous sections were obtained adopting quite conservative assumptions for Pop.III stars given the very low SF efficiency and the Salpeter-like IMF (although with  min = 1 ⊙ ).In the following sections, we explore the four main free parameters that regulate Pop.III star formation in meraxes listed in Table 1: the Pop.III star formation efficient  SF,III , the critical metallicity for Pop.III/II transition  crit , the critical surface density of cold gas for Pop.III star formation to occur Σ norm,III and the IMF.

Pop. III SF efficiency
The star formation efficiency determines the conversion of the cold gas into stars and is the free parameter with the largest impact on the Pop.III global SFRD and SMF.This parameter is largely unconstrained, with some simulations supporting very low values (10 −4 − 10 −3 Skinner & Wise 2020) and others suggesting higher values (Fukushima et al. 2020).We ran two simulations keeping all the Pop.III free parameters unchanged except for  SF,III , which we boosted by one order of magnitude (to the same value we use for Pop.II  SF = 0.08) and to a very high value ( SF = 0.5).
As shown in Fig. 11, Pop.III galaxies become more massive as  SF increases, erasing the "double peak" feature in the total SMF.For the model with  SF,III =  SF,II , the peak of the Pop.III SMF is still shifted to the left by half dex compared to the Pop.II.This is because, despite having the same star formation efficiency, Pop.III galaxies mostly form in mini-halos.For  SF,III = 0.5, the peak of the Pop.III and Pop.II SMF is at the same mass value (10 5 − 10 6  ⊙  at  ∼ 15 − 5) with the only difference being the larger variance for Pop.II galaxies.Increasing the Pop.III SF efficiency does not affect the Pop.II SMF (except for small variations at  = 5) meaning that the mechanical feedback does not change significantly.
The star formation efficiency also regulates the Pop.III/II transition as it directly correlates with the SFRD.Looking at the SFRD in the three different models (Fig. 12) we see that for  SF,III = 0.5, the SFRD is dominated by Pop.III up to  ∼ 15 compared to the fiducial model where Pop.II SF becomes larger than Pop.III at  > 20.The Pop. II SFRD does not significantly change between the three models, and so the change in the Pop.III SFRD also affects the total SFRD (the dashed line is almost one order of magnitude higher than the solid one).All the models converge at  ∼ 13 when even with the highest  SF,III the total SFRD is entirely dominated by the Pop.II contributions.

Critical metallicity
The critical metallicity defines the metallicity at which there is a change in the IMF (metal-free gas clouds have inefficient cooling that leads to less fragmentation and higher masses Chiaki & Yoshida 2022).There is currently no consensus on the value of  crit , with two competing models that assume the fragmentation driven by either the carbon and oxygen line (e.g.Chon et al. 2021) or the dust cooling (e.g.Schneider et al. 2012;Chiaki & Yoshida 2022).The first class of models determine a  crit ∼ 10 −2 −10 −3  ⊙ while the seconds give a lower value (∼ 10 −6  ⊙ ).Some studies also argue that the value of  crit evolves with redshift due to the effect of the CMB (Chon et al. 2022).Following simulations in the literature (Schneider et al. 2006;Visbal et al. 2020), in this work we adopted an intermediate value  crit = 10 −4  ⊙ as the fiducial value, hereafter we will also show results for the two extreme values of 10 −2  ⊙ and 10 −6  ⊙ .The SMF between the fiducial model and Z crit = 10 −6 Z ⊙ does not change.This is because, as we saw in the top right panel of Fig. 8, in the overdense regions where there is star formation, the metallicity becomes larger than 10 −4  ⊙ very quickly so that changing  crit to 10 −6  ⊙ does not further accelerate the Pop.III/II transition.However, considering a higher value changes the total SMF as we end up with many more Pop.III stars.As shown in Fig. 13, in the Pop.III SMF in this case, many of the previously identified as Pop.II galaxies are now Pop.III galaxies across the entire mass range.The high-mass tail of Pop.III SMF extends up to 10 7.5  ⊙ at  = 5 when  crit = 10 −2  ⊙ , while the Pop.II systems are less massive as they start to form later and so have less time to build up their mass.The combined effect of Pop.III and Pop.II makes the total SMF computed from the simulation with  crit = 10 −2  ⊙ single peaked as the high-mass peak coming from Pop. II stars are washed out even at  = 5 so that the total SMF peaks at  ★ ∼ 10 4  ⊙ at  = 10 − 5 with a high-mass tail that extends up to 10 9  ⊙ .
The critical metallicity also has a strong impact on the evolution of the SFRD (see Fig. 14) as it affects both the internal and external enrichment.While decreasing the value of Z crit from the fiducial value only mildly decreases the Pop.III SFRD without altering the total result, choosing Z crit = 10 −2  ⊙ strongly changes the SFRD history at  ≤ 18.The Pop. III SFRD increases by 1-2 orders of magnitude while the Pop.II decreases by the same amount.It also starts to increase steadily only from  ∼ 18. Increasing the Pop.III star formation and simultaneously decreasing the Pop.II impacts on the total SFRD.Since Pop.III stars form less efficiently ( SF,III = 0.1 ×  SF,II ) we predict a lower total SFRD compared to the fiducial model (and in this case, the two models do not converge until the very end of the simulation).Overall the critical metallicity heavily impacts both the SFRD and SMF during the Cosmic Dawn and the Epoch of Reionization and affects when the universe transitions from being Pop.III dominated Pop.II.

Critical surface density
In order to trigger the star formation in a galaxy, our model requires the gas density in the disk to be above a certain threshold.This result has been obtained in the observational results of Kennicutt (1998); Kennicutt & De Los Reyes (2021) and theoretically motivated by Kauffmann (1996); however, the value of this parameter is highly uncertain, especially at high- as it depends on the structure of the disk (e.g.thickness, turbulence etc.).We decided to duplicate this free parameter so that we have one for Pop.III and one for Pop.II.While keeping the Pop.II parameter fixed, we explore the extreme case of Σ norm,III = 0, which is equivalent to assuming that Pop.III stars start to form as soon as the cooling of the gas begins, and all cold gas is available to form Pop. III stars (in proportion to the star formation efficiency).This choice results in a larger abundance of low-mass Pop.III systems at high- (the largest difference between the cyan dashed and solid line in Fig. 15 are at  = 20 and 15).However, since the main changes in the Pop.III SF are in low-mass halos, setting Σ norm,III = 0 does not significantly change the SFRD history (see Fig. 16 where the cyan dashed line is only slightly lower than the fiducial model).
Overall, this free parameter is the one with the smallest impact on the star formation history, as it only impacts the low-end of the stellar mass function at high- for which there are no observational constraints.

IMF
The last parameter we explore is the IMF of Pop.III stars.Our fiducial model is quite similar to the Pop.II IMF as we assume a Salpeter IMF that favours the formation of low-mass stars.The only difference is that we take larger lower and upper limits (1 M ⊙ and 500 M ⊙ ).We consider two log-normal IMFs labelled logA and logE as in Tumlinson (2006) with the parameters summarized in Table 2.For these IMFs we have spectral energy distributions (from Raiter et al. 2010) which will be used in Section 4).The IMF logA is centred at  = 10 ⊙ , which enhances the probability of a Pop.III stars ending its life as a core-collapse SN. logE is centred at even a higher mass ( = 60 ⊙ ), so we expect most of the stars with very short lifetimes to end their lives by collapsing into a black hole without any supernova explosion.Compared to the Salpeter IMF, both log-normal IMFs do not have many low mass stars ( < 8 ⊙ ).
Decreasing the lifetime of Pop.III stars make the Pop.III SMF in Fig. 17 shift to lower masses, increasing the separation between the two peaks of the total SMF (the effect of changing the IMF towards higher masses is the opposite of increasing the SF efficiency).The difference between the fiducial model and one adopting the logA IMF is less than one dex.However, in the logE IMF model, most of the stars die in the same snapshot in which they form (the lifetime of a 60 ⊙ star is ∼ 4.5 Myr, which at  < 15 is less than the time separation between two consecutive snapshots in meraxes, see also Fig. 4).This effect strongly prevents the build-up of Pop.III systems, and in Fig. 17, we can see that the most massive Pop.III galaxies at  = 5 have  ∼ 10 5  ⊙ .There are no differences in the Pop.II SMF with different IMFs.
In difference from the SF efficiency, the IMF mostly affects the stellar mass function.This is because the main change is with respect to the lifetime of Pop.III stars.In Fig. 18 we see that there are no significant differences in the SFRD from adopting different IMFs, except for the first few snapshots of our simulation where the star formation is mildly reduced in the logE model.This result comes from the stronger feedback typical of the more top-heavy IMF (a larger number of massive stars implies more supernovae), and it has a stronger impact at  ≥ 24 when Pop. III star formation is dominating.
In conclusion, the global evolution of Pop.III stars are strongly dependent on the star formation efficiency because high  SF,III increases both the SMF and the SFRD.The evolution is also dependent on the IMF adopted (although the IMF affects only the SMF).The critical metallicity has an impact only if we consider a high value, and it mostly affects the self-enrichment as most of the galaxies in our simulation get their metals from their own star formation rather than from an external metal bubble.The critical surface density instead has only a small impact on the Pop.III star formation history and it is restricted to small mass halos and high-.

POP. III OBSERVABILITY WITH JWST
Unconstrained Pop.III free parameters in our model have an impact on the redshift evolution of the SMF and SFRD.However, our model predicts a consistent number density of Pop.III dominated systems with stellar masses between 10 3 − 10 5 M ⊙ .Here we will explore the UV luminosity of these systems using the model described in Section 2.7 focusing on the luminosity functions at very high-.
Given the same total stellar mass, Pop.III systems have larger luminosities compared to Pop.II galaxies, but shine for shorter times.This is demonstrated in figure 19 where we show the evolution in the first 10Myr of the SED of a Pop.III galaxy that forms 10 5  ⊙ Pop.III stars at z = 8 (top panels) and the absolute UV magnitude (bottom panels) at different times Δ since the star formation burst.Left, mid and right panels assume a Pop.III Salpeter, logA and logE IMF respectively.The black horizontal line corresponds to the magnitude computed assuming a continuous star formation throughout the snapshot.For reference, we also show the magnitude of a Pop.II galaxy forming the same amount of Pop.II stars (red dashed line).Firstly, we notice that Pop.III galaxies with star formation and a log-normal IMF are one magnitude brighter than Pop.II galaxies with the same SFR (while there is no significant difference if Pop.III stars follow a Salpeter IMF).Once we focus on the instantaneous star formation model for Pop.III stars, the brightness of the galaxy can be boosted (or reduced) by several magnitudes depending on the time Δ following the burst at which we are observing the system (see also Trussler et al. 2023).Even when we consider a Salpeter IMF, the brightness of a Pop.III galaxy changes ∼ 2 mags in ∼ 10Myr.This change is more dramatic for the log-normal IMFs (especially for the logE IMF where in 10 Myr Δ AB ≃ 7 mags).The large difference between the instantaneous and continuous star formation models shown in Fig. 19 for a single Pop.III galaxy will reflect in the entire population and potentially lead to a higher number density of extremely bright Pop.III galaxies.Fig. 20 shows the luminosity function (computed as described in the Appendix B).As a result of instantaneous Pop.III star formation, there are chances to observe these galaxies at a time when their luminosities are higher than what is given by the continuous scenario.This is reflected in Fig. 20 where we see the bright end of the solid cyan line shifting to higher luminosities.When looking at the total luminosity function we see a much significant Pop.III contribution to the number density of galaxies around M AB ∼ −16 when assuming a LogE IMF for their SF.
Given their larger brightness, a population of Pop.III galaxies with a top-heavy IMF has been suggested as a possible explanation for the abundance of bright galaxies at  > 10 (see, e.g.Trinca et al. 2023;Yung et al. 2023 andHarikane et al. 2023).Without invoking any exotic physics or a revision of the standard ΛCDM model of cosmology, other possible explanations are a combination of increased star formation efficiencies at high- and reduced feedback (Qin et al. 2023), bursty star formation (Sun et al. 2023), cosmic variance (Shen et al. 2023) and a modified ΛCDM power spectrum (Parashari & Laha 2023;Padmanabhan & Loeb 2023).Here, we show the predicted UV luminosity functions at  = 16, 12, 11, 10, 9 (we chose these values as the ones for which there are more JWST observations available).Results for different IMFs and different star formation efficiencies are summarized in Fig. 21 where all the other Pop.III parameters are taken as in the fiducial model.For all the models considered below, our total luminosity function agrees with early JWST observations at  ≤ 12 (except for the points at  AB < −19 where we do not have galaxies in meraxes).
With our fiducial parameters (dashed lines), both with the Salpeter and with the logA IMF (upper and middle row), we find that Pop.III systems are not the brightest galaxies at the redshift considered (M AB ≥ −14).However, when we consider a log-normal IMF with   21), and found a good agreement at  ≤ 12 and M AB ≲ −16.5 for all the models considered.The main differences are that at  = 9, we do not have any galaxy with M AB ∼ −20, and their model has a steeper evolution predicting many more faint galaxies.At  = 16 only models with  SF,III = 0.08 and a log-normal IMF reproduce their results at −14 ≤ M AB ≤ −18.
Even though the total UVLF are similar, their Pop.III contribution is significantly lower, especially when compared to the models with  SF,III = 0.08 (their results predict that less than 10% of bright galaxies at  = 15 − 16 host an active Pop.III stellar population).
Given that the choice of Pop.III free parameters are similar (they also take  SF,III ∼ 0.1 and account for a critical surface density of the cold gas before forming Pop.III stars).The main differences in the model that can explain the different results are second point, our model accounts for the fact that a Pop.III galaxies might be observed immediately after they formed stars.As shown in Fig. 20, this enables us to boost the luminosity of some Pop.III dominated galaxies by several magnitudes.Without accounting for this effect, our Pop.III galaxies would have similar luminosity to those predicted by Trinca et al. (2023).The main limitation of our results is that, given the small size of the simulation, we intrinsically miss the brightest galaxies (at  ≥ 9 we do not get galaxies brighter than M AB ≃ −19).Larger simulation volumes will allow us to be conclusive, but this result suggests that if Pop.III stars are still forming at  ≥ 10 − 12, and they have a top-heavy IMF, they do not require increased star formation efficiencies to outshine the Pop.II systems.As Pop.III dominated galaxies become much more rare at  ≲ 10, this scenario allows us to increase the abundance of bright galaxies only at  > 10 without boosting the bright end of the UV luminosity function at lower redshift hence we achieve a good agreement with early JWST observations from Donnan et al. ( 2023

CONCLUSIONS
In this paper, we studied Pop.III star formation in mini-halos with an updated version of the meraxes semi-analytic model of galaxy formation that includes Lyman-Werner background and the streaming velocities.We also implemented external metal enrichment following the growth of the supernova bubbles according to the Sedov-Taylor model.We computed the bubble size distribution function and found that our results agree with Trenti et al. (2009) (most bubbles are smaller than 150ℎ −1 ckpc at  = 6 and have a typical size of 100 − 200 ckpc at  = 5) hence most halos get self-enriched rather than externally polluted.Only low-mass halos (M vir ≤ 107.5  ⊙ ) at  < 10 are more likely to get their metals from the IGM.This is a consequence of the small size of the supernova bubbles that results in a small filling factor (0.1% at  = 12 and less than 1% at  = 5.)The free parameters allow us to explore the global properties of Pop.III dominated galaxies.We ran this model on top of a dark matter-only N-body simulation able to resolve all the mini-halos down to ∼ 3 × 10 5  ⊙ at  ≤ 30.We explored the impact on the SMF and SFRD of the main free parameters of our model.All models converge at  ≲ 10.However, Pop.III star formation efficiency and IMF lead to differences at high-.Finally, we investigated the SED evolution of a Pop.III dominated galaxy for different IMFs.The shorter lifetime of a Pop.III galaxy motivated us to use an instantaneous star formation model when computing the luminosity function.We computed the total and the Pop.III galaxy UV luminosity functions and compared these to the early JWST results in order to study whether the excess of bright galaxies at high- could be explained by a population of Pop.III dominated galaxies.Having explored different IMFs and star formation efficiencies, our model predicts, for a log-normal IMF with a characteristic mass of 60 M ⊙ and a Pop.III star formation efficiency comparable to the Pop.II, most of the brightest galaxies at  ≥ 12 and all at  = 16 (M AB = −18) are Pop.III dominated.Even with a smaller mass or a Salpeter IMF, Pop.III dominated systems still have an impact on the bright end of the UVLF at  = 16.In conclusion, this work supports the scenario for which top-heavy Pop.III dominated galaxies might explain the abundance of bright JWST galaxies at  ≥ 12 without requiring very high star formation efficiencies or extremely weak feedback at high-.

DATA AVAILABILITY
The main data presented in this work have been analysed with the dragons 7 and sector8 python packages publicly available.The latest version of meraxes is now publicly available on GitHub here9 .The output of the high-resolution N-body simulation used in this work can be downloaded from Zenodo at Balu et al. (2024).Further data and codes will be shared on reasonable request to the corresponding author. .left: mass fraction of false pristine (solid line) and false enriched (dashed line) galaxies (see text for details).Blue, red, green, black and grey lines refer to different grid resolutions (N = 16, 32, 64, 128, 256 respectively).center: absolute value of the difference between the mass fraction of false pristine and false enriched galaxies.The same colour coding is adopted.right: sum of the mass fraction of false pristine and false enriched galaxies.The same colour coding is adopted.

Figure 1 .Figure 2 .
Figure1.left: 2D slice of the L10 box at  = 10 with the position of all the galaxies with  vir ≥ 10 6  ⊙ (blue dots) highlighting the ones that are externally metal enriched (red dots).The green dots are galaxies with an associated metal bubble larger than 0.05 ℎ −1 cMpc.The slice is 10 ℎ −1 Mpc on a side and ∼ 0.078ℎ −1 Mpc thick.right: zoom-in of a single pixel (0.24 ℎ −1 Mpc on a side).The circle represents the actual size of the bubble around the green dot.

Figure 3 .
Figure 3. 2D projections of the IGM metallicity log( IGM ) in solar units (top row), the metal filling factor log() (middle row) and the overdensity  (bottom row).From left to right, the columns correspond to the same redshift:  = 20, 16, 12, 8.Each slice is 10 ℎ −1 Mpc on a side and considers the average contribution over the entire thickness of the box.In the top row, red (blue) represents enriched (pristine) pixels:  IGM ≥ (<)  crit .Green pixels have a larger filling factor , indicating a larger probability of a galaxy accreting enriched gas.Dark cells in the bottom row have larger mean overdensity .

Figure 6 .
Figure 6.Halo mass function at z = 20, 15, 10 and 5 of Pop.III and Pop.II hosting systems.The dashed line corresponds to the  vir = 10 4 K that defines the atomic cooling limit.

Figure 7 .
Figure 7. SFRD vs redshift for all stars (black line), Pop.III (cyan) and Pop.II (red) using the fiducial model (thin line) and the NoMini model (thick line).

Figure 8 .
Figure 8. top: Average metallicity of the IGM (in solar units) vs redshift highlighting  crit for the Pop.III/II transition.bottom: fraction of the volume with a metallicity  ≥  crit vs redshift.The left panels show these quantities averaged through the entire box, while the right panels consider only those cells with at least one metal bubble in them.Shaded regions in the right panels highlight the minimum and the maximum value (see text).

Figure 9 .
Figure 9. Halo mass function at z = 20, 15, 10 and 5 of Pop external (black) and internal (grey) metal enriched halos.The dashed line corresponds to the  vir = 10 4 K that defines the atomic cooling limit.

Figure 10 .
Figure 10.As Fig. 7 for the fiducial model (solid lines) and turning off the external metal enrichment (dashed lines).

Figure 19 .
Figure 19.top: time evolution of the SED of the brightest Pop.III galaxy at  = 8 assuming a Salpeter (left), logA (mid), logE (right) IMF.Lighter lines assume a smaller Δ since the star formation burst.bottom: Absolute UV magnitude of the same galaxy as a function of Δ in Myr.The black (red) horizontal line shows the magnitude value of the galaxy assuming a Pop.III (Pop.II) continuous star formation.

Figure 20 .
Figure 20.Total (black) and Pop.III (cyan) UV luminosity function at  = 16 assuming the Salpeter (top) and the logE (bottom) IMF.Solid lines are computed assuming that Pop.III stars form in a single burst, dashed lines take a continuous star formation model.
Figure A1.left: mass fraction of false pristine (solid line) and false enriched (dashed line) galaxies (see text for details).Blue, red, green, black and grey lines refer to different grid resolutions(N = 16, 32, 64, 128, 256 respectively).center: absolute value of the difference between the mass fraction of false pristine and false enriched galaxies.The same colour coding is adopted.right: sum of the mass fraction of false pristine and false enriched galaxies.The same colour coding is adopted.

Table 1 .
Qiu et al. (2019)10)( − ) the stellar mass formed at  −  with age between  and  + ,   () is the luminosity of a stellar population per unit mass and   () is the transmission function of the interstellar medium.For the latter term, we refer the reader toQiu et al. (2019)while   for each Pop.III IMF has been taken fromRaiter et al. (2010)as discussed above.This calculation is nearly identical toQiu et al. (2019), with the difference being that we do not have the metallicity dependence because Pop.III stars SEDs are Free parameters of Pop.III star formation.
defined with a zero metallicity.In this framework, we assume that the star formation occurs continuously throughout the snapshot.To compute ( − ), we take the stellar mass formed in that snapshot and average it over the entire duration of the snapshot.While this is a good approximation for Pop.II stars, Pop.III star formation is expected to occur in a single burst because of the feedback from Pop. III stars are likely to prevent continuous star formation.For this ★ see

Table 2 Table 2 .
IMF model parameters.Duration of the snapshot dt in Myr as a function of redshift.Dashed lines correspond to the lifetime of a zero metallicity star of 8M ⊙ , 10M ⊙ , 40M ⊙ and 60M ⊙ .

Table 1
(at  < 15, the total SFRD in the fiducial and NoMini model are superimposed).We see that accounting for Pop.III star formation in mini-halos is crucial during the Cosmic Dawn because Pop.III stars dominate the global star formation history at  ≥ 20 and are still relevant up to  = 18.As we can see in the lower panel, the Pop.III SFRD flattens at  ∼ 18 and starts to decrease at  ∼ 10.
Trinca et al. (2023)s of 60M ⊙ (lower row), Pop.III galaxies become significantly brighter (the cyan line approaches the black line) and at  = 12 and 16 when some of the brightest systems (M AB ∼ −16) are Pop.III dominated.Even though most of the Pop.III galaxies are still very faint (well below the sensitivity of JWST).This result may indicate that Pop.III dominated galaxies are present at  > 12.This result becomes more robust when considering  SF,III =  SF,II (solid lines).The bright end of the total and Pop.III dominated galaxy UVLF shifts of 1 − 2 magnitudes (depending on the IMF) at  = 16 and 12.At  ≤ 11, there is no or little difference in the total UVLF with different star formation efficiency even though Pop.III dominated galaxies are still ∼ 1 mag brighter.Models with a Salpeter or a log-normal IMF centered at 10 ⊙ predict that the bright end of the UVLF at  = 16 is impacted by Pop.III dominated systems while at  ≤ 12 none or very few of the brightest galaxies are Pop.III dominated.Overall, the abundance of bright galaxies at  ≥ 12 is better explained by the model with the logE Pop.III IMF and  This may suggest lower star formation efficiencies than what considered in this work making Pop.III galaxies much fainter and very unlikely to be observed.We compared our results withTrinca et al. (2023)(black pentagons in Fig. Xu et al. 2016;Sarmento & Scannapieco 2022) 12 we find Pop.III dominated galaxies with M AB ≃ −18 and at  = 16 all the brightest galaxies (M AB ≤ −16) are Pop.III dominated.The results discussed above considered a star formation efficiency  SF,III ≥ 0.008 which allows the formation of Pop.III stellar systems in meraxes with masses up to ∼ 10 5  ⊙ .These values are substantially higher than what predicted by some hydro-dynamical simulations (e.g.Xu et al. 2016;Sarmento & Scannapieco 2022).