A Global Semi-Analytic Model of the First Stars and Galaxies Including Dark Matter Halo Merger Histories

We present a new self-consistent semi-analytic model of the first stars and galaxies to explore the high-redshift ($z{>}15$) Population III (PopIII) and metal-enriched star formation histories. Our model includes the detailed merger history of dark matter halos generated with Monte Carlo merger trees. We calibrate the minimum halo mass for PopIII star formation from recent hydrodynamical cosmological simulations that simultaneously include the baryon-dark matter streaming velocity, Lyman-Werner (LW) feedback, and molecular hydrogen self-shielding. We find an overall increase in the resulting star formation rate density (SFRD) compared to calibrations based on previous simulations (e.g., the PopIII SFRD is over an order of magnitude higher at $z=35-15$). We evaluate the effect of the halo-to-halo scatter in this critical mass and find that it increases the PopIII stellar mass density by a factor of ${\sim}1.5$ at $z{>}15$. Additionally, we assess the impact of various semi-analytic/analytic prescriptions for halo assembly and star formation previously adopted in the literature. For example, we find that models assuming smooth halo growth computed via abundance matching predict SFRDs similar to the merger tree model for our fiducial model parameters, but that they may underestimate the PopIII SFRD in cases of strong LW feedback. Finally, we simulate sub-volumes of the Universe with our model both to quantify the reduction in total star formation in numerical simulations due to a lack of density fluctuations on spatial scales larger than the simulation box, and to determine spatial fluctuations in SFRD due to the diversity in halo abundances and merger histories.


INTRODUCTION
Theoretical calculations based on the standard model of cosmology predict that the first stars began illuminating the Universe within the first ∼100 Myr that followed the Big Bang.Simulations indicate these stars, categorized as Population III (PopIII) stars, formed from primordial metal-free gas within dark matter "minihalos" (M vir = 10 5 − 10 6 M ⊙ ), with much higher masses (M * = 10 − 1000 M ⊙ ) than metal-enriched stars (for a recent review, see Klessen & Glover 2023).Thus, the first PopIII stars likely had stellar lifetimes of only a few colton.feathers@rockets.utoledo.edumihir.kulkarni@utoledo.eduelijah.visbal@utoledo.eduryan.hazlett@rockets.utoledo.eduMyr (Schaerer 2002).Following their short lives, they injected metal-enriched material into their surroundings via supernovae (SNe) winds, resulting in the formation of Population II (PopII) stars (see e.g., Smith et al. 2015) and ultimately the first galaxies.
While there are currently no confirmed detections of PopIII stars, a variety of upcoming observations have the potential to constrain their abundance and properties (e.g., the PopIII initial mass function (IMF)).For example, the recently launched James Webb Space Telescope (JWST ) may observe pair-instability SNe (Hartwig et al. 2018a;Whalen et al. 2013) from PopIII stars with initial masses of ∼140 -250 M ⊙ .If large clusters of PopIII stars form in so-called PopIII galaxies (e.g., due to inefficient mixing of metals), they may be directly observable with JWST (Visbal et al. 2017;Kulkarni et al. 2019;Sarmento et al. 2019a;Vikaeus et al. 2022).We note that large PopIII starbursts detectable with JWST have recently been predicted to oc-cur in alternative dark matter scenarios that suppress small scale structure, such as "fuzzy dark matter" (e.g.Kulkarni et al. 2022;Ferreira 2021;Hui et al. 2017;Hu et al. 2000).While PopIII stars currently remain elusive, JWST has already observed metal-enriched galaxies out to very high redshifts (e.g.Naidu et al. 2022;Labbé et al. 2023;Finkelstein et al. 2023).The properties of lowmass galaxies may even depend on the characteristics of PopIII stars due to the hierarchical nature of structure formation in ΛCDM (e.g., Abe et al. 2021).
Another promising route to constrain properties of the first stars and galaxies is to measure their impact on the thermal and ionization properties of the intergalactic medium (IGM).For example, the optical depth due to electron scattering of the cosmic microwave background (CMB) has been used to put upper limits on the efficiency of PopIII star formation (e.g.Haiman & Holder 2003;Shull & Venkatesan 2008;Ahn et al. 2012;Visbal et al. 2015a;Miranda et al. 2017).Additionally, the cosmological 21cm signal is sensitive to the properties of early star formation.This is true both for global experiments such as EDGES (Bowman et al. 2018) and LEDA (Price et al. 2018), as well as those aiming to measure 3-dimensional (3D) spatial fluctuations, such as HERA (DeBoer et al. 2017) and SKA (Mellema et al. 2013).Additional observational probes of PopIII stars include gamma-ray bursts (Bromm & Loeb 2007;Burlon et al. 2016;Kinugawa et al. 2019), line-intensity mapping of the 1640 Å He II recombination line (Visbal et al. 2015b;Parsons et al. 2022), and studying extremely metal-poor stars in the local universe that could have been formed from the SNe ejecta of PopIII stars ("stellar archaeology" e.g., Magg et al. 2018;Hartwig et al. 2018bHartwig et al. , 2015;;Frebel & Norris 2015).
The measurements described above require accurate theoretical models of early star and galaxy formation to both suggest optimal observing strategies and interpret the results.Previous modelling efforts have taken many forms, from large-scale cosmological hydrodynamic simulations (e.g.Xu et al. 2013;Jaacks et al. 2018) to analytic calculations for quick predictions of global properties (e.g Haiman & Bryan 2006;Wyithe & Cen 2007;Trenti & Stiavelli 2009;Visbal et al. 2015a;Furlanetto & Mirocha 2022).While hydrodynamical cosmological simulations, such as the Renaissance Simulations (Xu et al. 2013(Xu et al. , 2014;;O'Shea et al. 2015;Smith et al. 2018;Xu et al. 2016;Ahn et al. 2015), include most of the physical processes relevant to early star formation, such as radiative, mechanical, and chemical feedback, they are very numerically expensive (e.g., ∼10 7 CPU hours per realization in the case of Renaissance).While these simulations follow the physics with high fidelity, they still typically involve free parameters in sub-grid prescriptions (e.g., the minimum metallicity for PopII star formation or the PopIII IMF), and are too costly to explore the parameter space in detail.
Semi-analytic models of the first stars and galaxies represent an alternate approach that can be used to rapidly survey the uncertain parameter space, while accurately following the hierarchical assembly of dark matter structure (e.g., Trenti & Stiavelli 2009;Yung et al. 2019a;Crosby et al. 2013;Valiante et al. 2016).The dark matter "backbone" can be computed either through cosmological N-body simulations or analytic Monte Carlo (MC) methods based on the Extended Press-Schechter (EPS) formalism (Press & Schechter 1974;Bond et al. 1991).Star formation and astrophysical feedback is then implemented on top of this structure using analytic techniques, allowing for rapid simulation across cosmological volumes and large redshift ranges while still including most of the important physics (see Magg et al. 2018;Mebane et al. 2018;Visbal et al. 2018Visbal et al. , 2020;;Liu & Bromm 2020, for examples).
In this paper, we present a new global high-redshift (z ≥ 15) semi-analytic model of the first stars and galaxies.The main new features of this model are that it includes complete dark matter (DM) halo merger histories (from MC methods) and is calibrated to state-ofthe-art hydrodynamical cosmological simulations of the critical halo mass for PopIII star formation (Kulkarni et al. 2021).We utilize our model to make updated predictions of the high-redshift PopIII and metal-enriched star formation rate densities (SFRDs) and to test the impact of various physical effects/parameters (e.g., the baryon-dark matter streaming velocity).Given that our model incorporates fully detailed DM halo merger history, we are also able to test the effect of approximations that have been adopted previously in literature to estimate the SFRD, such as smooth merger histories determined by halo abundance matching (e.g.Furlanetto et al. 2017;Mebane et al. 2018) and simple integrals of the halo mass function (e.g., Visbal et al. 2015a;Muñoz et al. 2022;Muñoz 2023).Finally, we use our model to simulate sub-regions of the Universe with different volumes, accounting for the Poisson fluctuations in the number of halos and their diverse merger histories.This allows us to estimate the impact of merger history on 3D spatial fluctuations of the PopIII and metal-enriched SFRDs relevant to future predictions of the cosmological 21cm signal.It also allows us to show the impact on the SFRDs predicted by numerical simulations of finite box size that effectively ignore density fluctuations on spatial scales larger than the simulation box.
The remainder of this paper is structured as follows.
In Section 2, we contextualize our new model by reviewing previous work in the literature.In Section 3, we discuss the details of our fiducial model and the physical processes included.The results of this work are divided into Sections 4, 5, and 6.Section 4 presents the main results of our fiducial model, and details the impact of various modelling fits and parameters on the global SFRD.Section 5 focuses on the impact of modelling choices for the dark matter halo evolution and the prescription for determining PopII star formation.Here we also check the accuracy of previous assumptions made in the literature regarding the effects of halo assembly on predicted SFRDs.In Section 6 we simulate finite, physically representative volumes of the Universe using our fiducial model.We compare our fiducial SFRDs with those previously published in the literature in Section 7. Finally, in Section 8 we summarize and discuss our results and explore future research.Unless otherwise stated, all distances are in comoving units.Throughout this work, we use a ΛCDM cosmology, consistent with Planck Collaboration et al. ( 2020) and with the following parameters: Ω m = 0.32, Ω Λ = 0.68, Ω b = 0.049, h = 0.67.

PREVIOUS WORKS
In this section, we briefly review semi-analytic work in the literature to provide context for our new model.The general approach in semi-analytic modelling of galaxy evolution is to generate dark matter halo merger trees, either with MC methods or cosmological N-body simulations, and then apply analytic prescriptions for star formation and astrophysical feedback processes.This has been widely used to model galaxies at lower redshifts than those explored here (for a review see Somerville & Davé 2015).It has also been used to study high-redshift galaxies and reionization, but without including PopIII stars (e.g., Poole et al. 2016;Mutch et al. 2016;Yung et al. 2019bYung et al. , 2020bYung et al. ,a, 2021Yung et al. , 2022)).
Described in detail in Section 3, we have adopted this framework to develop a global model of the first PopIII stars and early galaxies that includes full dark matter halo merger history, incorporates self-consistent LW feedback (Haiman et al. 1997(Haiman et al. , 2000;;Machacek et al. 2001;O'Shea & Norman 2008;Ahn et al. 2009), and is calibrated to recent hydrodynamical simulations.Most of the previous merger tree-based semi-analytic modelling of the first stars has been utilized to study the assembly of individual halos, rather than the global quantities which we focus on.This includes the formation of high-redshift supermassive black holes in rare high-mass halos (e.g., Valiante et al. 2016;Lupi et al. 2021;Trinca et al. 2022) as well as Milky Way-like halos to make pre-dictions for stellar archaeology (e.g., Magg et al. 2018;Griffen et al. 2018;Ishiyama et al. 2016;de Bennassuti et al. 2017;Graziani et al. 2017).We note that there have been previous global merger-tree-based models related to ours.In Magg et al. (2016), the PopIII SFRD was the focus; this is updated here with newly simulated critical halo masses for PopIII star formation (Kulkarni et al. 2021).More recently, Trinca et al. (2022) utilized EPS merger trees to study black hole and PopIII star formation, but this work focused on lower redshifts than what are explored here.We also note that towards the completion of this work, Ventura et al. (2022) and Trinca et al. (2023) presented similar global models of the first stars and galaxies.Their respective focuses were on the predicted 21cm signal and UV luminosity functions, whereas we utilize our model to explore the impact of various physical effects and modelling techniques on the PopIII and metal-enriched SFRDs.An additional novel aspect of our work is that we use our model to understand the behavior of various sub-volumes of the Universe as a step towards accurate modelling of the 3D spatial fluctuations of high-redshift star formation.
Several recent semi-analytic works have included 3D spatial information from N-body simulations (in addition to merger trees) to model the effects of ionization feedback and "external" metal-enrichment between neighboring halos (e.g., Sarmento et al. 2019b;Liu & Bromm 2020;Visbal et al. 2020;Hartwig et al. 2022).This is complementary to the work presented here.Although our fiducial model cannot properly include 3D feedback (since MC EPS methods are used to generate the merger trees), it has been found that 3D feedback does not have a significant impact of the SFRDs at z ≳ 15 (Visbal et al. 2020), which is when we have decided to stop our models presented below.With our semi-analytic framework, we are able to rapidly explore much larger volumes than what is possible when N-body simulations are used.
In addition to these semi-analytic methods, analytic work has been used to model the global abundance of the first stars and galaxies.These include what we refer to as "smooth accretion models" where halo growth histories are approximated using abundance matching of the halo mass function (Furlanetto et al. 2017;Mebane et al. 2018Mebane et al. , 2020)).A number of previous works have also estimated the SFRDs of PopIII and metal-enriched stars via analytic integration of the halo mass function (e.g., Visbal et al. 2015a;Muñoz et al. 2022;Muñoz 2023).We use our model below to test the impact that these approximations related to halo assembly have on the global star formation history.
We also note that our new model relates to largescale semi-numerical models which have been used to predict the high-redshift 21cm signal (Mesinger et al. 2011;McQuinn & O'Leary 2012;Visbal et al. 2012;Fialkov et al. 2013;Fialkov & Barkana 2014;Kaur et al. 2022;Magg et al. 2022).These models simulate vast cosmological volumes (e.g., boxes ∼ 1 Gpc across) with individual resolution elements that are a few Mpc across.Within these volume elements, the local SFRD is typically computed with integrals of the halo mass function (e.g.Jaacks et al. 2018), though recent work has been calibrated with merger tree-based semi-analytic models (Magg et al. 2022).Utilizing our model to simulate variations within individual sub-volumes (as we present in Section 6) can lay the groundwork to improve future sub-grid prescriptions within large-scale semi-numerical models.

FIDUCIAL SEMI-ANALYTIC MODEL
In this section, we introduce our fiducial semi-analytic model whose purpose is to self-consistently determine the evolution of the global SFRD (along with other aspects of early star formation such as the scatter in SFRD within sub-volumes of the Universe) given some input DM halo model, star formation prescription, and parameter values.We note that the efficiency of this framework allows each realization to be completed in ∼75 minutes on a single CPU.Our fiducial modelling parameters are summarized in Table 1 and described in more detail below.

Dark Matter Halo Merger History
Our model uses MC merger trees based on EPS formalism (Press & Schechter 1974;Bond et al. 1991) to model DM halo growth.We generate merger trees following the prescription in Lacey & Cole (1993), which assumes binary mergers of DM halos at each redshift step, separated by ∆z.For a given final halo mass and redshift, the merger history is modelled by iteratively stepping back through time, determining the first progenitor halo mass via the EPS mass-weighted progenitor function, and setting the second progenitor halo mass such that the sum of the two progenitors is equal to the merged descendent halo mass.This is done for each halo by solving equation 5 of Visbal et al. (2014) for σ M1 , the root-mean-square (rms) density fluctuation on a scale corresponding to mass M 1 , (1) Here, σ M0 is the rms density fluctuation for the scale corresponding to the present halo mass, δ crit = 1.686 is the critical overdensity in linear perturbation theory, D(z) is the linear growth factor, and x is a randomly drawn value between 0-1.From σ M1 we determine the progenitor halo mass, M 1 , and the merging halo progenitor mass, M 2 = M 0 − M 1 .This process is repeated at each step back through cosmic time until a halo falls below the resolution halo mass, M res , after which no further progenitors are modelled.To account for any statistically rare merger histories, we generate multiple merger trees for each halo mass bin whose masses are set at the end of the simulation, z = 15.
In our fiducial model we use 36 logarithmically spaced mass bins from 10 5.6 −10 9.1 M ⊙ , finding this mass range to be sufficiently converged for estimating the global SFRD between z = 15 − 35.These precise mass values have little physical motivation, we simply find that adding halo mass bins above 10 9.1 M ⊙ only affects the results at the earliest redshifts (z ≳ 42), and that the critical mass for PopIII star formation is higher than 10 5.6 M ⊙ at z = 15 in all reported runs, indicating lower mass bins would not contribute to the resulting SFRDs.We also varied the number of halo mass bins and found that our fiducial number of bins yields SFRD values that are converged to within 10% on average between z = 15 − 35.Within each mass bin we generate 10 merger trees each initialized at z = 15, finding this number of trees to be the most reasonable trade-off between the required computation time and the statistical noise in our resulting SFRDs.We adopt a resolution halo mass of M res = 5 × 10 4 M ⊙ for our merger trees, a threshold well below typical values of M crit so that no early star formation is missed.We run our models from z = 60-15, with redshift steps of ∆z = 0.05.Overall, we find the results of our work to be sufficiently converged for predicting the globally-averaged SFRD over the redshift range z = 15 − 35.

Star Formation
Here we discuss the methods used to determine the formation of both PopIII and PopII stellar mass within merger tree halos.Note that we do not track individual stars in our model, and so any PopII and PopIII stellar masses discussed are total mass values.We focus mainly on the high-redshift transition from PopIII to PopII and leave any analysis of the transition from PopII to PopI for future work.

Critical Mass & PopIII Star Formation
We use a simple "instantaneous" prescription for PopIII star formation in our fiducial model.A given halo without previous star formation will undergo a single episode of pristine star formation, creating a constant stellar mass of M III,new = 200 M ⊙ immediately after it surpasses the critical mass for star formation, M crit .This is roughly consistent with the average PopIII stellar mass found within halos by Skinner & Wise (2020) (see their Figure 8), and approximately corresponds to a constant PopIII star formation efficiency of f III = 0.0001 for a halo of mass 10 6.1 M ⊙ .
The critical mass is a global threshold for PopIII star formation, determined at each time step by finding the smaller of the atomic and molecular hydrogen cooling halo masses, i.e. (2) Since the gas within early DM halos is pristine, molecular hydrogen (H 2 ) is initially the most efficient gas coolant.As star formation proceeds, the resulting radiative feedback dissociates H 2 molecules, effectively raising the halo virial temperature and therefore M H2 .For halos with virial temperatures T vir ≳ 10 4 K, atomic hydrogen cooling becomes the more efficient mode of halo cooling, and so we adopt our atomic cooling mass from Visbal et al. (2020), This mass is based on the threshold found in hydrodynamic cosmological simulations, with precise values from Fernandez et al. (2014).The H 2 cooling mass used in our fiducial model is adopted from the hydrodynamical simulations of Kulkarni et al. (2021), which are the first (along with Schauer et al. 2021) to include a dependency on both the relative baryon-DM streaming velocity, v bc , and the LW radiation background intensity, J LW (also see Nebrin et al. 2023).They also consider an improved treatment of H 2 self-shielding, ultimately leading to an analytic fitting function that is well-calibrated to the hydrodyamical simulations: Here, the critical mass at z = 20, M z20 , and the power index, α, are both functions of the LW background intensity and baryon-DM streaming velocity (see Kulkarni et al. (2021) for the best-fit values of these parameters to the simulations).At each time step, our semi-analytic model determines which halos have masses above M crit and forms M III,new within each.In Section 4.1, we compare the results of our fiducial model, with M H2 determined by equation 4, to a previous model of M H2 based on the simulations of Machacek et al. (2001) (also see Stacy et al. 2011;Greif et al. 2011;Fialkov et al. 2012).This molecular mass fit, or variations thereof, has commonly been used in simulations of early star formation throughout the literature.Here we implement the form used in Visbal et al. (2014) which includes a dependence on the redshift: M H2 = 2.5 × 10 5 1 + z 26 −1.5 1 + 6.96(4πJ LW ) 0.47 .
(5) The key difference between this and our fiducial equation for M H2 is that the fit presented in Kulkarni et al. (2021) includes the effects of H 2 self-shielding and a stronger relative baryon-DM streaming velocity dependence.Fitting functions based on previous simulations, including those of Machacek et al. (2001), either do not include v bc , or assumed that the LW intensity and v bc were independent modelling processes and their effects on M crit were multiplicative (Fialkov et al. 2012).As shown by Kulkarni et al. (2021), their effects on M crit are not multiplicative, hence our choice of this M H2 fit for our fiducial model.

Delay Time for Enriched Star Formation
A key parameter of our semi-analytic model is the delay time between the formation of PopIII and metalenriched PopII stars, t delay .This delay period is a result of the radiative and supernovae feedback from the first stars which heats and ejects gas from within the hosting minihalos.Once the halo gas reservoir recollects and cools once more, the enriched material begins forming PopII stars (Chiaki et al. 2017(Chiaki et al. , 2013;;Jeon et al. 2014).We choose t delay = 10 Myr for our fiducial model, and explore the impact of different t delay values on the resulting SFRDs in Section 4.4.

PopII Star Formation
We assume that once PopIII stars form in a halo and the delay time has elapsed, it begins forming metalenriched PopII stars.For PopII star formation, we adopt the system of differential equations described in Furlanetto & Mirocha (2022) that simultaneously determines the gas and stellar masses of a given halo through a simplified "bathtub" model (e.g.Bouché et al. 2010;Davé et al. 2012;Dekel & Mandelker 2014): Here, Ṁ * ,II is the halo PopII star formation rate, f II is the PopII star formation efficiency per halo freefall time, t ff , which we assume to be 10% of the Hubble time (i.e.t ff = 0.1t H ), Ṁgas represents the time derivative of the halo gas mass, Ṁacc is the mass accretion rate of the halo, and η SN is the supernova ejection efficiency (described in more detail below).By analytically solving these differential equations at each time step, we obtain the gas and stellar masses of each halo, and determine the gas mass ejected through SNe by subtracting these two masses from the initial total baryonic mass of the halo.This prescription forms stars such that a fraction of the total halo gas mass, f II , is converted into stellar mass over one halo freefall time, allowing isolated DM halos to continue forming PopII stars as long as there is gas within them.While we note that there is large uncertainty in the star formation efficiency (SFE) of unobserved low-mass galaxies at very high redshifts (z ≳ 15), we have aimed to select a realistic value for our fiducial choice of f II = 0.0025.When combined with the SNe feedback prescription described below, we find that this choice of f II leads to a SFRD consistent with Visbal et al. (2018), which was calibrated to observations of z∼6 UV galaxy luminosity functions (Bouwens et al. 2015).In a related upcoming study (Hazlett et al., in prep.), we have created a model of PopII star formation directly calibrated to the Renaissance simulations.We find that this calibration yields PopII SFRDs that agree to within a factor of two at z ≲ 32 with those predicted in our fiducial model presented here.
At each time step, we determine the PopII and PopIII star formation rates (SFRs) for each merger tree by dividing any new stellar mass formed by the length of the time step.We average the SFRs of all merger trees in each mass bin, and weight each average by the corresponding halo mass bin number density, n i , as determined from the Sheth-Tormen halo mass function (Sheth & Tormen 1999) at bin center, i.e.
Here, dn/dM is the Sheth-Tormen halo mass function and n i is the number density of halos at z = 15 in mass bin i, bounded by upper and lower halo masses m u and m l .For the remainder of this paper, all DM halo mass bin number densities are the z = 15 values unless otherwise stated.The weighted averages of all mass bins are then summed to give the global PopII and PopIII SFRD values for that redshift step, i.e.

SF RD(z)
Here, Ṁ * ,tot is the total SFR of every halo at redshift z within a merger tree in mass bin i. Averaging the total SFRs of all merger trees in mass bin i gives ⟨ Ṁ * ,tot ⟩ i which we then weight by n i and sum across all mass bins.

Feedback Processes
Our semi-analytical framework also includes various crucial feedback processes that impact global SFRD evolution.In this subsection, we will discuss the auxiliary features and physical feedback processes of our fiducial model and how they affect early star formation.

Lyman-Werner Background Intensity
LW feedback must be considered in any semi-analytic model of early star formation.As the first stars begin to shine and emit radiation, photons in the energy range E LW = 11.2-13.6eV freely stream out up to a horizon of ∼100 Mpc and into other DM halos, dissociating H 2 molecules within them (Haiman et al. 1997(Haiman et al. , 2000;;Machacek et al. 2001;O'Shea & Norman 2008;Ahn et al. 2009).In our fiducial model, we self-consistently determine the LW background intensity at each time step in units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 using the J LW calculation from Visbal et al. (2020): (10) Here, ϵ LW (z ′ ) is the mean LW emissivity, t H is the Hubble time, and f LW (z ′ , z) is the attenuation of the LW flux from z ′ to z as LW photons redshift into Lyman series lines and are absorbed (Haiman et al. 1997).We approximate this LW attenuation using equation 22 in Ahn et al. (2009), and determine ϵ LW (z ′ ) by summing the contributions of both stellar populations.
We assume the number of LW photons produced per stellar baryon to be η III = 65000 for PopIII stars (Schaerer 2002), and η II = 4000 for metal-enriched stars (Samui et al. 2007;Incatasciato et al. 2023).We set these to be equal to the number of ionizing photons produced per stellar baryon for simplicity, and also note that our assumed value of η III may be an overestimate, but we have found that the PopII SFRD quickly dominates the LW background in all runs, and varying η III only changes the PopIII SFRD by a few percent.An increase in J LW results in an increase in both M z20 and α in equation 4, giving a larger M H2 for equation 2. Increased LW radiation intensity therefore results in an overall increase to the critical mass for PopIII star formation, suppressing further star formation.
In our model, we do not consider the effects of Xray feedback on PopIII star formation.We note that Hegde & Furlanetto (2023) recently found that the Xray background is not a dominant effect for PopIII star formation in minihalos at the redshifts considered in this work (z > 15, also see Ricotti 2016;Ricotti & Ostriker 2004).

Relative Baryon-DM Streaming Velocity
Another important parameter in modelling early star formation is the relative streaming velocity between baryonic and dark matter (Tseliakhovich & Hirata 2010;Tseliakhovich et al. 2011;Fialkov et al. 2013;Greif et al. 2011;Stacy et al. 2012).At cosmic recombination, photons decoupled from the baryonic matter which was then free to interact gravitationally.Dark matter density perturbations, however, could grow under the influence of gravity before this.As a result, we are left with a Maxwell-Boltzmann distribution of relative baryon-DM streaming velocities with a rms value of σ vbc = 30 km s −1 at recombination, which decreases with redshift as v bc (z) = v bc,0 ((1 + z)/1100).This relative motion is roughly constant over scales of ∼3 Mpc, and is crucial for modelling early star formation because the DM halos are unable to efficiently capture and cool the high-speed gas streaming by them.An increase in streaming velocity leads to a decrease in the amount of gas bound to halos and delays halo gas cooling, thus leading to an overall increase in M crit .The value of this velocity is commonly referred to in terms of σ vbc , and while we explore the impact of various σ vbc values below, we adopt a fiducial global streaming velocity of v bc = 1σ vbc = 30 km s −1 .

Feedback from Supernovae
A fraction of the gas mass within a given halo will be ejected by the star formation and SNe occurring at a given time step.The gas mass ejected by PopII SNe is determined by the SN ejection efficiency of the halo, η SN .In this work, we follow a similar prescription for determining η SN to the one presented in Sassano et al. (2021): where E SN is the average energy per PopII supernova, ϵ SN is an efficiency parameter, R SN is the rate of SNe per solar mass, and v esc is the escape velocity of the halo, which we assume to be the halo circular velocity determined by equation 25 of Barkana & Loeb (2001).We note that the circular and escape velocities at a halo's virial radius differ by a factor of √ 2, however the effect that SNe ejection has on the global SFRD is comparatively low (i.e.removing all SNe feedback results in a ∼ 30% increase to the z = 15 PopII SFRD), and since we do not assume a preferred radius for SNe ejection from our DM halos, we adopt the circular velocity for simplicity.We adopt values of E SN = 10 51 ergs and ϵ SN = 0.0016 from Wise et al. (2012).R SN is determined by dividing the SN energy per solar mass, E SN,M⊙ , by E SN .Adopting a value of E SN,M⊙ = 6.8×10 48ergs M −1 ⊙ (Wise et al. 2012) gives us R SN = 0.0068 M −1 ⊙ .The numerator of equation 11 can then be described in terms of velocity, and so we normalize this expression for a halo at z = 20 with a mass equal to the atomic cooling mass (equation 3).This gives us η SN = 5.23 × (v esc /14.46 km s −1 ) −2 .
The ejected gas mass, M ej , is determined at each time step by subtracting the final stellar and gas masses within the halo from its initial total baryonic mass.M ej is then removed from the available halo gas mass until it is reintroduced for future star formation after one freefall time, which we have assumed to be 10% of the Hubble time at the moment of ejection.

RESULTS I: FIDUCIAL MODEL
In this section we present the global SFRDs of our fiducial model.We illustrate the effects of a new critical mass model (equation 4) that more accurately considers the effects of LW feedback, v bc , and H 2 self-shielding, as well as the impact of a halo-to-halo scatter on M crit .We also assess the impact that various physical processes and parameter value choices have on the results.Unless otherwise stated, all SFRD values plotted throughout this work have been smoothed post-simulation with a running average over ∆z = 1.We also provide analytic fits for a sample of the LW backgrounds that result from this research in Appendix A.

Critical Mass Model
Figure 1 shows the SFRDs and corresponding critical masses resulting from both our fiducial model with M crit based on the new simulations of Kulkarni et al. (2021), and those resulting from the Visbal et al. (2020) adaptation of the commonly used critical mass model consistent with the simulations of Machacek et al. (2001), Greif et al. (2011), andStacy et al. (2012) in which the effects of H 2 self-shielding are not considered.This "no self-shielding" M crit model stems from the work of Fialkov et al. (2012), which updated the critical mass threshold from Machacek et al. (2001) to include more sophisticated treatments for v bc -dependence (also see equation 3 of Visbal et al. (2020) and accompanying text for more information on including LW and v bcdependence).However, since recent works have revised the critical mass threshold for PopIII star formation to include the effects of H 2 self-shielding (e.g.Kulkarni et al. 2021;Schauer et al. 2021;Nebrin et al. 2023), we refer to this earlier M crit as the "no self-shielding" critical mass model for the remainder of this paper.
The left set of plots is for a global streaming velocity of 1σ vbc .Here we see that our fiducial critical mass is a factor of a few lower than the no self-shielding model throughout, falling to a factor of ∼6 times lower at z ∼ 15.This lower M crit leads to earlier star formation and an overall larger number of star-forming halos in the fiducial model.We thus find a PopIII SFRD that is over an order of magnitude higher than the no self-shielding PopIII SFRD throughout most of the redshifts shown.The fiducial PopII SFRD is also consistently higher than the one given by the no self-shielding M crit , but this difference falls from an order of magnitude at z = 35 to a factor of two times the no self-shielding PopII SFRD at z = 15.
The right set of plots in Figure 1 show the results for no baryon-DM streaming.Here we see a somewhat smaller difference in the SFRDs.At z = 35, our fiducial model gives SFRD values that are factors of 30 and 8 times higher than the no self-shielding values for PopIII and PopII, respectively.The final PopII SFRD of the fiducial model is only 30% higher than the no self-shielding case, whereas the final PopIII SFRD is a factor of ∼2 higher.Here, the effects that cause the differences in the 1σ vbc case are reduced because the critical masses fall within a factor of five to one another throughout, owing to the stronger v bc dependence in the no self-shielding model.

Critical Mass Scatter
In our fiducial model we assume M crit to be a global constant at each time step; we now consider the impact of a scatter on M crit (z) which effectively gives each DM halo an individual critical mass threshold for PopIII star formation, M crit,halo .To more accurately account for variations in halo mass assembly and geometry, Kulkarni et al. (2021) quantified this scatter in terms of J LW and v bc which we may implement via log 10 (M crit,halo ) = log 10 (M crit ) + R. (12) Here, the value of R is determined for each DM halo after every halo freefall time by randomly drawing from a Gaussian distribution centered on zero and with a standard deviation of σ = 0.15.The M crit scatter values reported in Figure 7 of Kulkarni et al. (2021) represent the difference between the 25 th and 75 th percentile of critical masses, equivalent to 1.34σ.For a streaming velocity of 1σ vbc , this value ranges from 0.15-0.3,and we adopt the value 0.2, roughly corresponding a LW background intensity of J LW = 1 × 10 −21 erg s −1 cm −2 Hz −1 sr −1 .Although we self-consistently determine J LW at each time step, its values at z ≲ 30 are comparable to unity, and so we conservatively adopt σ = 0.15 (≃ 0.2/1.34)for the standard deviation of our distribution.(for more information, see Figure 7 of Kulkarni et al. 2021).
We illustrate the impact of this relative scatter on M crit in Figure 2. The solid lines in the top panel show the stellar mass densities (SMDs) for both PopII and PopIII as determined by our fiducial model.The SMDs of ten simulations including M crit scatter are shown in light gray for each stellar population, and the average of these ten SMDs are represented by the corresponding dashed lines.We choose to show the SMDs in Figure 2 rather than the SFRDs to more clearly show the impact that M crit scatter has on the star formation resulting from our semi-analytic model.We see from the bottom panel of Figure 2 that introducing M crit scatter increases the PopIII SMD by a factor of ∼2 at early times (z ≳ 30).Below z ∼ 30, the average PopIII SMD with scatter maintains a roughly 50% increase over the no scatter case.This increase is due to halos drawing low random R values for equation 12, resulting in halo critical masses lower than the global mean value and earlier star formation than without scatter.The average PopII SMD with scatter is also increased by a factor of ∼2 at the highest redshifts shown, before tapering off and coming into agreement with the no scatter SMD at the percent level by z = 15.PopIII star formation relies solely on a halo overcoming M crit , and halos are much more likely to randomly draw a low enough R value to cause earlier star formation than they are to consistently draw high R values that keep M crit,halo > M crit .We therefore see a consistent increase in the PopIII SMD that is not strongly reflected by the PopII as the PopII SFRD does not depend as directly on M crit .

Baryon-Dark Matter Streaming
In Figure 3, we present PopII and PopIII SFRDs and their corresponding M crit values for global streaming velocities of 0-3σ vbc (0-90 km s −1 ).Recall that the velocities stated are those at recombination, which fall off with decreasing redshift.We also show the average SFRD determined by using many global streaming velocities ranging from 0-3σ vbc and weighting each by the probability of the v bc value used.The probabilities are determined from the Maxwell-Boltzmann distribution of v bc values found at recombination (Fialkov et al. 2012;Tseliakhovich et al. 2011;Tseliakhovich & Hirata 2010), and these weighted averages are shown by the dashed purple curves in the top and middle panels of Figure 3 for PopII and PopIII SFRDs, respectively.
We see from the top two panels of Figure 3 that the SFRDs of both stellar populations are increasingly suppressed with faster streaming velocities.This is because higher v bc values increase M crit , which delays any PopIII and subsequent PopII star formation.This effect is greater for PopIII since it is directly governed by M crit , whereas PopII star formation begins after t delay and depends on the gas mass of the host halo thereafter.The weighted average SFRDs each agree with those of the fiducial 1σ vbc case to within 50% at z < 35, indicating that our fiducial treatment is a reasonable representation of the global streaming velocity.We note that the true baryon-DM streaming velocity varies over spatial distances of a few Mpc, a detail that we plan to incorporate into future work.

Additional Parameter Variations
In this subsection, we explore the effect of various parameter value choices on the resulting SFRDs.To show the impact of our fiducial choices for the PopII SFE and PopIII formation mass, we vary the values of f II and M III,new and plot the resulting SFRDs in the top row of Top: PopII SFRDs for four models with global streaming velocities of 0σ vbc = 0 km s −1 (blue), 1σ vbc = 30 km s −1 (orange), 2σ vbc = 60 km s −1 (green), 3σ vbc = 90 km s −1 (red).The average PopII SFRD is also shown (purple dashed), with each SFRD weighted by the probabilities of the streaming velocity used, which is set by a Maxwell-Boltzmann distribution.Middle: Same as the top panel, but for the PopIII SFRDs.Bottom: The value of Mcrit over time for each of the constant velocity runs.The close agreement between the weighted average SFRDs and those of the 1σ vbc case justifies our fiducial streaming velocity of 30 km s −1 .
PopIII by the corresponding vertical black lines.We see that higher f II values cause the PopII SFRD to first dominate over the PopIII at earlier times, ranging from z ≈ 27 for 0.5f II to z ≈ 32 for 4f II .
Looking to the top right panel of Figure 4, we see the SFRDs for three realizations with varying PopIII formation mass values ranging from M III,new = 100-800 M ⊙ , as well as a realization with a constant PopIII SFE of f III = 0.001 which instantly converts a fraction f III   We show both the MC EPS approach implemented in our fiducial model (gray curves) as well as the smooth accretion method (blue curve) from Furlanetto et al. ( 2017) (each described in detail in Sections 3 and 5, respectively).
To show the effect of the supernova ejection efficiency, we vary the fiducial value of η SN in our model and present the global SFRDs in the bottom right panel of Figure 4.Here we see that η SN mostly affects the PopII SFRD at late times, with lower efficiencies yielding higher PopII SFRDs overall.At redshift 15, the 0.5η SN case is a factor of ∼5 higher than the 10η SN case.The effect of supernovae feedback increases with time, partially due to the redshift dependence of v esc in equation 11 causing more gas to escape over cosmic time.

RESULTS II: IMPACT OF ENRICHED STAR FORMATION & DM HALO MODELS
In this section, we address the SFRD changes that result from alternative methods for determining the DM halo mass evolution and enriched star formation.We present a direct comparison of the SFRDs resulting from MC merger trees and smooth accretion models of DM halo evolution (as in Furlanetto et al. 2017).We also compare the results of our fiducial enriched star formation prescription (equations 6 & 7) to those of a simpler instantaneous model of star formation (e.g.Sun & Furlanetto 2016;Park et al. 2019;Magg et al. 2018Magg et al. , 2022)).For completeness, we also include an analytic approach reliant on halo mass function integration to estimate the SFRD (e.g.Muñoz 2023;Muñoz et al. 2022;Mashian et al. 2016;Visbal et al. 2015b).

Descriptions of Alternative Models
Starting with our "smooth accretion model" of DM halo evolution, we adopt the abundance matching model described in Furlanetto et al. (2017), where n h (m|z i ) is the number density n h , of halos with masses ≥ m, at redshift z i , as determined by the Sheth-Tormen halo mass function (Sheth & Tormen 1999).For a given starting mass m 1 at redshift z 1 , we determine the halo mass growth history for a constant n h by solving for m 2 at each subsequent time step (see Yamaguchi et al. 2022, for a recent example).Figure 5 shows an example smooth accretion halo mass growth track along with a sampling of MC merger trees with the same final halo mass.While the smooth accretion model follows a constant cosmic number density through time for each halo mass, MC merger trees initialize n h at z = 15 and determine the merger history for a given halo through EPS formalism.The two models produce halo mass evolution tracks that differ on average by a factor of ∼2 at z ≲ 35 due to this discrepancy between Sheth-Tormen and EPS formalism.We note, however, that this discrepancy in the halo mass growth between models is outweighed by the uncertainties in other high redshift astrophysical parameters implemented here, and thus represents a viable comparison to our fiducial model.We intend to explore alternative DM halo evolution models in future work.
For the "instantaneous" star formation prescription, we compare the SFRDs of our fiducial model with those determining enriched star formation via Here, M * ,II is the newly-formed PopII stellar mass, f II,i is the instantaneous PopII SFE, and M gas,new is the gas mass that has accreted into the halo over the last time step.This method uses a fraction, f II,i , of the gas mass accreting onto a halo at a given time step to instantly form PopII stellar mass within it.This means that if a halo already has PopII stars, it will not form any more PopII stellar mass until a halo with no prior enriched star formation merges with it and introduces new gas.
In the case of two halos with previous PopII star formation merging, no new stellar mass is formed.We include this method into our comparison as it is implemented by semi-analytic models throughout the literature in various forms that resemble equation 14; using a SFE and gas mass to instantly introduce new stellar mass (e.g.Sun & Furlanetto 2016;Park et al. 2019;Magg et al. 2018Magg et al. , 2022)).For the rest of this paper, we refer to this PopII star formation prescription as the "instantaneous" method.for determining halo merger history and for determining PopII star formation.Models utilizing the EPS merger tree method are denoted by "Trees" in the legend, and those utilizing the smooth accretion method for halo mass growth (equation 13) are labeled "Smooth".Models that use the fiducial set of ODEs for PopII star formation (equations 6 & 7) are labeled as "ODE", and those using the Instantaneous method (equation 14) are denoted by "Inst" in the legend.Also shown is the fully analytic Integral method (purple), as described in Section 5.1 with a PopII SFE of fII = 0.005 and a new PopIII stellar mass of MIII,new = 400 M⊙.Finally, the black "Trees -S-T" curve uses the fiducial ODE PopII star formation method and EPS merger trees for the DM halo evolution, but the halo mass bin number densities, ni, are weighted to match the Sheth-Tormen halo mass function at each redshift step, as is the case for the smooth accretion models.As discussed in the main text, we find that the smooth accretion gives similar results to the full merger tree models.Additionally, we find that instantaneous PopII SFRDs change less rapidly with redshift than our fiducial ODE-based prescription.
For completeness, we also consider the analytic model detailed in Visbal et al. (2015a) to compare with the results of our semi-analytic framework.This approach determines the SFRD by integrating the Sheth-Tormen halo mass function, dn/dM , at each time step (also see Mashian et al. 2016;Muñoz et al. 2022;Muñoz 2023, for another example of this method).Halo mass ranges are defined for both PopII and PopIII such that the cosmic mass fraction collapsed into DM halos, F coll , can be determined for stellar population i via Here, Ω m is the cosmological density parameter for matter and ρ c is the critical density.= 400 M ⊙ .We also note that we smooth the value of F coll,i over the previous ten time steps (∆z = 0.5) to avoid numerical feedback effects when determining its derivative.

Comparison of Methods
Figure 6 shows the global SFRDs by stellar population for various combinations of DM halo model and star formation prescription, as well as the SFRDs of the fully analytic approach.A key takeaway of Figure 6 is that the merger trees and smooth accretion models for DM halo evolution give similar PopII SFRDs.For a given PopII star formation prescription, instantaneous or ODE, the PopII SFRDs of both DM halo models agree to within a factor of two over the redshifts shown.
We also find that the fiducial ODEs for enriched star formation yield steeper PopII SFRD growth than the instantaneous treatment.Calculating the total change in PopII SFRD over the change in time shown in Figure 6 (i.e.δ ρ * ,II /δt), the slopes given by the ODE prescription are steeper than those of the instantaneous prescription by over a factor of two, given the same DM halo model.In our fiducial prescription, an isolated halo will form stars at a continually decreasing rate, whereas the in-stantaneous model requires fresh infalling gas, hence the shallower growth of the instantaneous PopII SFRDs.
Looking to the right panel of Figure 6, we find that the PopIII SFRDs of the merger trees have similar qualitative behavior to the smooth accretion SFRDs, but are roughly double their value over the entire redshift range shown.This is a consequence of the discrepancy in halo number densities between Sheth-Tormen and EPS formalism as described in Section 5.1, causing the smooth accretion models to underestimate the PopIII SFRDs by a factor of ∼2.We verified this by determining the halo mass function given by our EPS merger trees at each redshift then calculating its ratio with respect to the Sheth-Tormen halo mass function.We then used these ratio values to weight the SFRs of each individual halo before calculating the SFRD, effectively constraining the halo number densities of each mass bin to match the Sheth-Tormen halo mass function throughout the simulation.The resulting PopIII SFRD in Figure 6 (black) agrees with those of the smooth accretion model to within 15% between z = 20 − 35.At lower redshifts the Smooth -ODE curve decreases due to its higher PopII SFRD and LW background (discussed further in section 5.3), and the Sheth-Tormen curve approaches the fiducial SFRD values as the weighting ratios approach unity at z = 15.
Looking to the SFRDs determined by analytic integration in Figure 6, we see overall reasonable agreement with the SFRDs of the semi-analytic models.In fact, when compared to the fiducial SFRDs, the analytic values agree to within a factor of four and to within ∼ 65% over the redshifts shown for PopII and PopIII, respectively.The integral PopII SFRD shows better agreement with the instantaneous star formation models, falling within a factor of two and to within 40% over the plotted redshifts for the merger trees and smooth accretion models, respectively.This agreement, however, is a result of tuning the star formation in the analytic integration model to match the results of our semi-analytic model.We used a PopII SFE and new PopIII stellar mass that were each twice their fiducial value, and set the boundary between PopIII and PopII star formation to be 2.5M crit (z).Variations of these parameter values can alter the resulting SFRDs by factors of a few.So while it is able to rapidly give order of magnitude agreement with the SFRDs of our semi-analytic framework, purely analytic models miss out on key astrophysics for early star formation, and so tuning is required to match semi-analytic values.

High J LW Limiting Smooth Accretion Models
In Figure 6, we see a clear downturn in the "Smooth -ODE" PopIII SFRD at z ≲ 20 (green, right panel).Although this effect is also seen in the fiducial PopIII SFRD (blue, same panel), it is much more pronounced for the smooth accretion model than it is for the merger trees.This behavior results from a limitation of the smooth accretion model that we have discovered in cases of high LW background intensity.To illustrate the breakdown that smooth accretion models experience in high J LW environments, we artificially increase the LW background by a factor of ten and compare the resulting SFRDs for merger trees and smooth accretion halos in Figure 7.
The top panel of Figure 7 shows the PopII and PopIII SFRDs for both MC merger trees and smooth accretion halo models, given a LW background intensity an order of magnitude higher than in our fiducial framework.The bottom panel shows the critical mass of the smooth accretion realization on top of a sample of the smooth accretion halo mass growth tracks used.The green mass tracks depict halos that cross M crit at any single point in the simulation, whereas the dashed purple mass tracks do not ever cross M crit .
With a boosted LW background intensity, the critical mass increases more rapidly than in our fiducial model.At z ∼ 19, the last smooth accretion track crosses M crit , after which no further PopIII star formation occurs.This suggests that smooth accretion models for DM halo mass evolution break down in high-radiation environments, and are unable to form stars as the growth of M crit outpaces the growth of the halos.This is particularly relevant for M crit models with a z-dependence that scales with J LW , as in the Kulkarni et al. (2021) model.MC merger trees, with their stochastic halo mass growth histories, cover a wider range of halo masses at all time steps, making the rapid growth of M crit a nonissue for MC merger tree models.Thus, the complete termination of PopIII star formation in high LW environments represents a limitation to semi-analytic models with smooth halo mass growth determined via cosmic number density abundance matching.

RESULTS III: REALISTIC BOX VOLUME SIMULATIONS
In addition to studying the global star formation history, we also use our model to represent sub-volumes of the Universe.We accomplish this by building "boxes" that simulate an ensemble of halos with realistic merger histories for a specified volume.While we cannot model spatial information within these boxes (e.g., the 3D positions of halos), we can compute the variation from region-to-region due to differences in halo number and assembly history.As mentioned above, this is relevant to the large-scale modelling of spatial fluctuations in the high-redshift 21cm signal which are sourced by fluctuations in star formation.This approach also allows us to determine the extent to which finite simulation boxes (e.g., for hydrodynamical cosmological simulations) systematically underestimate the SFRD due to their lack of density fluctuations on spatial scales larger than the simulation box.
To simulate a finite-volume box, we begin by determining the mean number of halos in each mass bin (which are the same as described in Section 3) at the final redshift of our simulation, z = 15.This number is computed using the Sheth-Tormen halo mass function with a power spectrum truncated (i.e.set to zero) at k min = 2π/( √ 3L box ), corresponding to the longest length scale in the box.This truncation is applied to approximate the impact of boxes with density equal to the cosmic mean as is typically assumed for numerical cosmological simulations.We note that in future work we intend to explore simulating similar boxes with different mean overdensities, which could then inform the sub-grid modelling in large-scale semi-numerical models of the 21cm signal (Visbal et al. 2012;Fialkov et al. 2013;Magg et al. 2022).We sample the number of halos for each mass bin from a Poisson distribution and generate the MC merger trees for each as described in Section 3.1.Our fiducial semi-analytic model is then applied to these merger trees assuming an external LW background, J LW (z), computed from our global fiducial model.
We present the results of our finite-volume boxes in Figure 8.Here we show the mean PopII and PopIII SFRDs for 10 boxes of various volumes, and compare them to the SFRDs of our fiducial global model presented above.For boxes smaller than 3 Mpc across, we find that removing power associated with scales larger than the boxes results in a systematic reduction of the PopII SFRD.The PopIII SFRD however, is less impacted due to it being sourced by smaller halos (with abundance set by the power on smaller scales).From this, we conclude that a 3 Mpc box is sufficiently large to provide an unbiased estimate of the global star formation at z = 15 − 35.One caveat is that this assumes an accurate model of J LW (i.e. from our global fiducial semi-analytic simulation) and would not hold if the LW background were computed self-consistently from the finite-volume box.We also note that if the SFRD were to be dominated by more massive halos than in our model, larger-volume boxes would be required for an unbiased SFRD (as likely occurs at lower redshift).
In Figure 8, we also show the scatter of the PopII and PopIII SFRDs for different box volumes.In order to prevent the scatter from depending on our redshift step size, we assign each PopIII star formation burst to a random time between z and z + ∆z, then smooth the results over a 3 Myr period (roughly corresponding to a typical PopIII stellar lifetime).This process effectively allows for star formation to proceed continuously over time instead of at discrete redshift values which bin star formation events into finite time spans, thereby introducing a dependence on the bin size (i.e.∆z) into the error calculations.As expected, the SFRDs increase with higher redshift as halos become increasingly rare.We also note that the fluctuations in PopIII are significantly higher, in large part due to PopIII star formation occurring in instantaneous bursts rather than smoothly across time as is the case for PopII in our model.The scatter is the error one would expect in a single simulation of the given box size.Additionally, these are the levels of scatter one expects in sub-regions of large-scale semi-numerical simulations, as these works typically assume (3 Mpc) 3 volume resolution elements corresponding to the scales where v bc does not vary spatially.We find order unity scatter at z ∼ 35 and 10% at z ∼ 22.5 for the PopII SFRD.For the PopIII, we predict order unity at z ∼ 30 and a value above 20% for all our redshifts (z > 15).
In the future, we plan to extend this analysis to overdense/underdense regions to improve the modelling of large-scale 3D spatial fluctuations in star formation history.We note that the size fluctuations shown here represent lower limits because we assume fixed star formation efficiencies (i.e.constant f II and M III,new ).In reality these will vary from halo to halo, and vary over time for a specific halo.In future work, we intend to calibrate the SFEs with hydrodynamical cosmological simulations including star formation and radiative feedback (Hazlett et al., in prep).

STAR FORMATION HISTORY COMPARISON WITH OTHER WORKS
We now look to previous works throughout the literature and compare their SFRD histories to those found with our semi-analytic framework.Figure 9 shows the PopII and PopIII SFRDs determined by our fiducial model, as well as those using the no self-shielding critical mass (see Figure 1), compared to several recently reported SFRDs from the literature.These SFRDs are the results of a variety of modelling prescriptions including semi-analytic simulations (e.g.Magg et al. 2018;Mebane et al. 2018;Visbal et al. 2018Visbal et al. , 2020;;Liu & Bromm 2020), cosmological hydrodynamic frameworks (Jaacks et al. 2018;Sarmento et al. 2019b), and even a fully analytic model (Muñoz et al. 2022).We note that the majority of the curves shown in Figure 9 Hartwig et al. 2022H&F 2023Jaacks et al. 2018L&B 2020Magg et al. 2018Mebane et al. 2018Munoz et al. 2021Sarmento et al. 2019Visbal et al. 2018Visbal et al. 2020 Figure 9.A comparison of the PopII (left) and PopIII (right) SFRDs found in this work with those from the literature.We show our fiducial SFRDs (thick black) and those found using the no self-shielding critical mass model (thick gray, see Figure 1) alongside recently published SFRD histories.These SFRDs alternate between dashed and solid curves for clarity, and are labelled alphabetically by first-author surname.Note, for conciseness, Hegde & Furlanetto (2023) and Liu & Bromm (2020) are labelled as "H&F 2023" and "L&B 2020", respectively.We also note that not all referred works include PopII SFRD histories.
ing that only broad trends are captured.Exceptions include Jaacks et al. (2018) PopII SFRD which has an assumed functional form, Liu & Bromm (2020) provide a fit for their PopIII SFRD, and the authors of Hegde & Furlanetto (2023) kindly provided us with raw SFRD values resulting from their fiducial framework.
In the right panel of Figure 9, we see that the majority of the PopIII SFRDs shown fall below our fiducial SFRD; however, many of them come into agreement when compared with our no self-shielding M crit realization.The SFRDs reported by Mebane et al. (2018), Visbal et al. (2018), Liu &Bromm (2020), andVisbal et al. (2020) all agree with the no self-shielding M crit PopIII SFRD to within a factor of a few throughout the redshift range shown.In particular, Mebane et al. (2018) and Visbal et al. (2018) are consistent to within a factor of ≤2 across the entire redshift range, despite their respective use of smooth accretion and N-body prescriptions for DM halo evolution.We therefore attribute the general agreement of these semi-analytic models with our no self-shielding M crit PopIII SFRD to their use of the same prescription for the critical mass (equation 5).
Focusing on a few particular cases in Figure 9, the PopIII SFRD found by Hartwig et al. (2022) shows qualitative agreement with our fiducial model at z ≲ 27, albeit a factor of ∼3 higher.Although they adopt a critical mass model with H 2 self-shielding (Schauer et al. 2021), their J LW is predetermined as a function of redshift which sets their M crit to five times the threshold found in our self-consistent model at z = 35.Once PopIII star formation begins, however, we infer that their high SFE of f III = 0.38 per freefall time (calibrated from observations), paired with the limited negative feedback of their analytic J LW (z) prescription, gave this overall higher PopIII SFRD.
Most recently, the semi-analytic model presented in Hegde & Furlanetto (2023) has one of the steepest PopIII SFRD evolutions shown in Figure 9.While their critical mass agrees with our fiducial M crit to within a factor of a few throughout, their star formation prescription allows for multiple generations of PopIII stars per halo, a feature that most likely sources such rapid growth.
The PopIII SFRDs resulting from cosmological hydrodynamic simulations in Figure 9 (i.e.Jaacks et al. 2018;Sarmento et al. 2019b) are also significantly steeper than those of our semi-analytic framework.Since these models mainly cover lower redshifts than what is explored here (each ending at z ∼7), their SFRDs at the redshifts shown in Figure 9 are likely products of the largest halos found in the simulation.Without resolving minihalos at higher redshifts, star formation is confined to the few resolved DM halos that persist at such early times, causing the SFRD to fall off more rapidly.If unresolved minihalos were included, they likely would supplement early star formation and carry the SFRDs at higher redshifts to give more qualitative agreement with our model.Now looking to the left panel of Figure 9, we see the available corresponding PopII SFRDs which result from a diverse range of star formation prescriptions.On top of this diversity, PopII star formation is frequently affected by other modelling aspects as well (such as the M crit model, f II , t delay , SNe feedback, etc.), which makes diagnosing variations between the resulting SFRDs more challenging than with PopIII.For example, Hegde & Furlanetto (2023) determine PopII star formation using our fiducial ODEs adopted from Furlanetto & Mirocha ( 2022), yet their PopII SFRD grows more consistently than the SFRDs resulting from our framework.This is likely a result of their PopII SFE which is determined at each time step using a star formation duty cycle, unlike the constant f II value used here.Conversely, the SFRD of Muñoz et al. (2022) is determined via analytic integration of the halo mass function which, as discussed in Section 5.2, heavily depends on the halo mass ranges over which one integrates and the SFEs used.The PopII SFRDs given by Mebane et al. (2018), Visbal et al. (2018, 2020), and Hartwig et al. (2022) each arise from prescriptions similar to our instantaneous PopII star formation model, but most likely differ in evolution as a result of varying PopII SFEs, varying the halo gas mass used, the inclusion of bursty star formation, and the DM halo model itself.
The SFRDs resulting from this research have extensively shown that one's modelling choices can significantly alter the resulting star formation history.In particular, the right panel of Figure 9 tells us that one must use the improved M crit from Kulkarni et al. (2021) or similar works (e.g.Schauer et al. 2021;Nebrin et al. 2023) that accurately accounts for the effects of H 2 selfshielding and baryon-DM streaming.In general, varying the SFEs of either stellar population can shift their SFRD values up or down with time, and the longer star formation is delayed, the steeper the SFRD growth once it commences.

CONCLUSIONS & FUTURE WORK
In this paper, we presented a new global semi-analytic model of the first stars and galaxies at redshifts z ≥ 15.Our model includes complete DM halo merger histories (typically not included in global star formation models at such high redshifts) and is calibrated to the latest numerical simulations (Kulkarni et al. 2021).Here we conclude by summarizing key results from our study and mention directions of future work.
We found that, compared to previous calibrations in the literature, the updated critical halo mass significantly changes the high-redshift SFRDs (mainly due to lower M H2 values as a consequence of the molecular hydrogen self-shielding).For instance in our fiducial model, the PopIII SFRD is increased by more than an order of magnitude compared to the no self-shielding critical masses based on the simulations of Machacek et al. (2001); Greif et al. (2011);Stacy et al. (2012).The PopII SFRD is increased by roughly one order of magnitude at z = 35, but falls to a factor of two disagreement by z = 15 (see Figure 1).
We also assessed the impact of including the scatter in M H2 due to individual differences in halo merger history and geometry between halos from the simulations of Kulkarni et al. (2021).Figure 2 shows that including this scatter modestly increases the abundance of stars formed.It increased the final PopIII stellar mass density by a factor of ∼1.5, and the final PopII SMD only increased by ∼10% relative to the no scatter case.
We also studied the effect that DM halo mergers have on star formation history by comparing the SFRDs resulting from MC merger trees to those resulting from smooth accretion based on abundance matching, as described in Furlanetto et al. (2017) (Figure 5).For our fiducial astrophysical parameters, we found that the resulting SFRDs (both PopIII and PopII) each agreed to within a factor of two from z = 18-31 for these DM halo models (see Figure 6).This suggests that the merger history is likely not required to model global SFRDs, provided one does not need an accuracy much greater than a factor of two, but it is required for spatial fluctuations such as those we show in Figure 8.We note, however, that our model does not incorporate changes in astrophysics due to mergers (e.g., an enhanced PopII star formation efficiency following mergers), which could occur in more sophisticated star formation models.We also point out that in models with very high J LW (for instance if the star formation efficiency was significantly higher than we assume), smooth accretion models have an unphysical shutoff in primordial star formation and a detailed treatment of the merger history may be necessary (see Figure 7).Additionally, we tested how the global star formation history is affected by our PopII star formation prescription.We compared our fiducial model, based on equations 6 and 7 (Furlanetto & Mirocha 2022) to a simpler "instantaneous" method used previously in the literature and found that the instantaneous model leads to a shallower metal-enriched SFRD slope throughout the simulation (see Figure 6).Our fiducial prescription yields steeper PopII SFRDs because it allows halos to continue forming stars in isolation, unlike the instantaneous model which requires infalling gas at each time step (resulting in relatively more late-time formation).The effects of this steeper enriched SFRD manifest in the PopIII SFRD at late times (z ≲ 25) when the increasing LW background intensity sufficiently increases M crit to alter the PopIII star formation history as well.
We also compared our fiducial merger-tree model to a simple integration of the halo mass function based on Visbal et al. (2015b).We found that this analytic calculation closely reproduces the PopII SFRD in the instantaneous prescription, and can reproduce our fiducial enriched SFRD to within a factor of ∼2.It does a more reasonable job for PopIII, agreeing with our fiducial PopIII SFRD to within a ∼ 65% at the redshifts we explored, but we again note that these SFRDs were the result of tuning the analytic prescription to more closely align with the semi-analytic SFRD values.
We showed a comparison of the SFRDs resulting from our semi-analytic framework to those from the literature in Figure 9, and found that one's modelling choices can significantly alter the resulting star formation history.In particular, we found that the reduction of M crit due to the effects of H 2 self-shielding and baryon-DM streaming allows for a PopIII SFRD that is about an order of magnitude higher than without the inclusion of such effects.We therefore recommend the use of the critical mass threshold from Kulkarni et al. (2021) or similar works (e.g.Schauer et al. 2021;Nebrin et al. 2023) that accurately accounts for the effects of H 2 self-shielding and baryon-DM streaming for predicting global SFRD evolution.
However, since these works investigating the effects of H 2 self-shielding are relatively recent, many uncertainties remain.Although Kulkarni et al. (2021) and Schauer et al. (2021) both utilize hydrodynamical cosmological simulations (Enzo and Arepo, respectively) to probe various redshifts, LW background intensities, and streaming velocities to determine the critical mass for star formation, the thresholds reported by the two works differ by a factor a few.While differences in the simulations used, the parameter spaces probed, and the assumed thresholds for halo cooling/collapse may have contributed to this, exact sources of this discrepancy remain unclear.More recently, Nebrin et al. (2023) published their critical mass model that builds up from halo gas physics to include dependencies on J LW and v bc using (semi-)analytic methods.However, despite hopes that a more analytic approach could shed light on the discrepancy, their resulting M crit seemingly falls between the two models of (Kulkarni et al. 2021) and Schauer et al. (2021).We therefore caution the reader to take care when selecting a particular critical mass threshold for future star formation models.
Finally, we used our model to simulate finite-volume boxes each with their own unique set of halos and merger histories (though not including 3D spatial positions).For boxes which are forced to have a density equal to the cosmic mean (as is typically done in non-zoom cosmological numerical simulations), we tested how the lack of power on scales larger the simulation box can impact estimates of the star formation history.We found that a (3 Mpc) 3 box is sufficiently large to provide an unbiased estimate of the global star formation at z = 15−35, provided that a reliable external LW background is known (for analytic fits of various LW backgrounds found by our semi-analytic model, see Appendix A).We also determined the typical scatter in the PopIII/PopII SFRDs as a function of redshift in regions of various volumes that are set to the mean cosmic density.
In future work, we will apply our model to predict large-scale spatial fluctuations in high-redshift star formation, which can then be used to model the cosmological 21cm signal.We will extend our finite-volume box analysis presented in Section 6 to include variations in star formation as a function of local overdensity and baryon-dark matter streaming velocity in (3 Mpc) 3 volume elements (as is typically done in other 21cm predictions, e.g.Visbal et al. 2012;Fialkov et al. 2013).We will then explore ways to rapidly emulate our model, such that it is possible to self-consistently determine the LW background and associated feedback across ∼Gpc distance scales while also accounting for variations in halo abundance and merger history down to the critical halo mass for forming PopIII stars.
We acknowledge support from NSF grant AST-2009309 and NASA ATP grant 80NSSC22K0629.All numerical calculations were carried out at the Ohio Supercomputer Center (OSC).We thank Sahil Hegde for kindly providing us with SFRD data from Hegde & Furlanetto (2023) for Figure 9. Table 2. Polynomial coefficients used in equation A1 for the JLW(z) fits shown in Figure 10.

FiducialFigure 1 .
Figure 1.The impact of different critical mass calculation methods for realizations including (left column) and not including (right column) relative baryon-DM streaming.Top left: Comparison of the PopII (blue) and PopIII (orange) SFRDs with different methods for determining Mcrit including a global streaming velocity of 1σ = 30 km s −1 .Our fiducial framework using the Kulkarni et al. (2021) Mcrit model (solid) is compared to a realization using the critical mass model from Visbal et al. (2020), based on the simulations of Machacek et al. (2001) (dashed).Top right: Same as top left, but for realizations with no global baryon-DM streaming (i.e.v bc = 0 km s −1 ).Bottom left: The value of Mcrit(z) for each method with a streaming velocity of 1σ vbc .Bottom right: Same as bottom left but with no baryon-DM streaming.

Figure 2 .
Figure 2. The impact of a critical mass scatter on the results of our fiducial model.Top: The evolution of the stellar mass density over time.The fiducial PopII (blue) and PopIII (orange) SMDs without scatter are shown by the solid lines.The densities for the sample of ten runs including Mcrit scatter are plotted in gray, and the average of these ten runs for both stellar populations are shown by the corresponding dashed lines.Bottom: Ratio of the average SMD including Mcrit scatter to the fiducial model SMD for each stellar population.

Figure 4 .Figure 3 .
Figure 3.The impact of the relative baryon-dark matter streaming velocity on the global SFRD and critical mass.Top: PopII SFRDs for four models with global streaming velocities of 0σ vbc = 0 km s −1 (blue), 1σ vbc = 30 km s −1 (orange), 2σ vbc = 60 km s −1 (green), 3σ vbc = 90 km s −1 (red).The average PopII SFRD is also shown (purple dashed), with each SFRD weighted by the probabilities of the streaming velocity used, which is set by a Maxwell-Boltzmann distribution.Middle: Same as the top panel, but for the PopIII SFRDs.Bottom: The value of Mcrit over time for each of the constant velocity runs.The close agreement between the weighted average SFRDs and those of the 1σ vbc case justifies our fiducial streaming velocity of 30 km s −1 .

Figure 5 .
Figure5.Comparison of mass growth histories for DM halos with a final masses of log 10 (M halo /M⊙) = 7.5 at z = 15.We show both the MC EPS approach implemented in our fiducial model (gray curves) as well as the smooth accretion method (blue curve) fromFurlanetto et al. (2017) (each described in detail in Sections 3 and 5, respectively).

Figure 6 .
Figure6.Time-evolution of the global PopII (left) and PopIII (right) SFRDs for different combinations of the methods used for determining halo merger history and for determining PopII star formation.Models utilizing the EPS merger tree method are denoted by "Trees" in the legend, and those utilizing the smooth accretion method for halo mass growth (equation 13) are labeled "Smooth".Models that use the fiducial set of ODEs for PopII star formation (equations 6 & 7) are labeled as "ODE", and those using the Instantaneous method (equation 14) are denoted by "Inst" in the legend.Also shown is the fully analytic Integral method (purple), as described in Section 5.1 with a PopII SFE of fII = 0.005 and a new PopIII stellar mass of MIII,new = 400 M⊙.Finally, the black "Trees -S-T" curve uses the fiducial ODE PopII star formation method and EPS merger trees for the DM halo evolution, but the halo mass bin number densities, ni, are weighted to match the Sheth-Tormen halo mass function at each redshift step, as is the case for the smooth accretion models.As discussed in the main text, we find that the smooth accretion gives similar results to the full merger tree models.Additionally, we find that instantaneous PopII SFRDs change less rapidly with redshift than our fiducial ODE-based prescription.

Figure 7 .
Figure 7.The breakdown of smooth accretion DM halo models in high LW background intensity cases.Top: The PopII (blue) and PopIII (orange) SFRDs for models utilizing MC merger trees (solid) and smooth accretion (dashed) for DM halo mass evolution, and an artificially boosted LW background of 10JLW.Bottom: The critical mass for the smooth accretion model shown in the top panel (thick black) alongside a sample of the halo mass growth tracks used in said model.The mass growth tracks that cross Mcrit at any point (and therefore form PopIII stellar mass) are shown in green, and those that do not cross Mcrit at all are shown by the purple dashed lines.

Figure 8 .
Figure8.Top panels: The average PopII (left) and PopIII (right) SFRDs for a range of box sizes.Each curve is a 10realization average and is denoted by the length of one side of the box volume.The shaded regions corresponding to each curve represent the standard deviation of the ten constituent runs.Also shown are the SFRDs of our fiducial model (black, both panels), which is a globally-averaged model and hence, a supposed upper limit for an ever-increasing box size.Bottom panels:The error on the mean for each of the PopII (left) and PopIII (right) SFRDs, i.e. the ratio of the standard deviation to the average for each corresponding curve in the top panels.

Table 1 .
Parameters and fiducial values adopted in our framework for modelling the global SFRD.
governed by the onset of PopIII star formation and t delay .The t delay parameter represents the delay between PopIII and metal-enriched star formation resulting from feedback associated with star formation and SNe winds.To show the impact of this, we vary our fiducial delay time and present the resulting global SFRDs in the bottom left panel of Figure4.Here we see the SFRDs of five different time delays ranging from 5-100 Myr.Unsurprisingly, changes in t delay mostly affect the timing of the initial PopII star formation event.Regardless of the length of the delay, however, enriched star formation quickly overtakes the PopIII SFRD afterwards and gives similar final PopII SFRD values, all agreeing to within 75% at z = 15.
(dashed dark orange, for clarity).Bottom left: SFRDs for five different values of t delay : 5 Myr (dotted), 10 Myr (solid, fiducial), 20 Myr (dashed), 50 Myr (dot-dashed), and 100 Myr (dot-dot-dashed).Bottom right: SFRDs for different supernova gas ejection efficiencies: 0.5ηSN (dotted), ηSN (solid, fiducial), 2ηSN (dashed), 5ηSN (dot-dashed), and 10ηSN (dot-dot-dashed).crit and suppresses further star formation.The PopII SFRDs, however, are virtually unchanged by the varying PopIII star formation prescrip-tions since enriched star formation is mainly In this work, we determine F coll at each time step by integrating the halo mass function from M crit -2.5M crit for PopIII, and from 2.5M crit -10 9 M ⊙ for PopII.The time derivatives for both stellar populations, (dF coll /dt) i , are then calculated analytically at each step and are multiplied by ρ b f II and Ω m ρ c M III,new to give the PopII and PopIII SFRDs, respectively.To more closely match the semi-analytic SFRD values, we double the fiducial value of both the PopII SFE and the new PopIII stellar mass in this analytic integration model, i.e. f II = 0.005 and M III,new