# A Framework to Calibrate a Semi-analytic Model of the First Stars and
Galaxies to the Renaissance Simulations

Ryan Hazlett Ritter Astrophysical Research Center, Department of Physics and
Astronomy, University of Toledo, 2801 W. Bancroft Street, Toledo, OH 43606,
USA<EMAIL_ADDRESS>Mihir Kulkarni Ritter Astrophysical
Research Center, Department of Physics and Astronomy, University of Toledo,
2801 W. Bancroft Street, Toledo, OH 43606, USA<EMAIL_ADDRESS>Eli
Visbal Ritter Astrophysical Research Center, Department of Physics and
Astronomy, University of Toledo, 2801 W. Bancroft Street, Toledo, OH 43606,
USA<EMAIL_ADDRESS>John H. Wise Center for Relativistic
Astrophysics, Georgia Institute of Technology, 837 State Street, Atlanta, GA
30332, USA<EMAIL_ADDRESS>
###### Abstract

We present a method that calibrates a semi-analytic model to the _Renaissance
Simulations_ , a suite of cosmological hydrodynamical simulations with high-
redshift galaxy formation. This approach combines the strengths of semi-
analytic techniques and hydrodynamical simulations, enabling the extension to
larger volumes and lower redshifts inaccessible to simulations due to
computational expense. Using a sample of _Renaissance_ star formation
histories (SFHs) from an average density region of the Universe, we construct
a four parameter prescription for metal-enriched star formation characterized
by an initial bursty stage followed by a steady stage where stars are formed
at constant efficiencies. Our model also includes a treatment of Population
III star formation where a minimum halo mass and log-normal distribution of
stellar mass are adopted to match the numerical simulations. Star formation is
generally well reproduced for halos with masses
$\lesssim$$10^{9}M_{\mathrm{\odot}}$. Between $11<z<25$ our model produces
metal-enriched star formation rate densities (SFRDs) that typically agree with
_Renaissance_ within a factor of $\sim$2 for the average density region.
Additionally, the total metal-enriched stellar mass only differs from
_Renaissance_ by about $10\%$ at $z\sim 11$. For regions that are either more
overdense or rarefied not included in the calibration, we produce metal-
enriched SFRDs that agree with _Renaissance_ within a factor of $\sim$2 at
high-$z$, but eventually differ by higher factors for later times. This is
likely due to environmental dependencies not included in the model. Our star
formation prescriptions can easily be adopted in other analytic or semi-
analytic works to match our calibration to _Renaissance_.

Population III stars (1285); Galaxy formation (595); High-redshift galaxies
(734); Cosmology (343)

††facilities: Ohio Supercomputer Center (OSC)††software: numpy (Harris et al.,
2020), scipy (Virtanen et al., 2020), matplotlib (Hunter, 2007), cython
(Behnel et al., 2011), jupyter notebook (Kluyver et al., 2016), rockstar
(Behroozi et al., 2013a), consistent trees (Behroozi et al., 2013b), yt (Turk
et al., 2011), ytree (Smith & Lang, 2019).

##  1 Introduction

Understanding the first stars and galaxies is one of the major goals of modern
astronomy. In the context of the $\Lambda$-Cold Dark Matter ($\Lambda$CDM)
model of cosmology, simulations and analytic calculations predict that the
first Population III (Pop III) stars formed from primordial gas within
${\sim}10^{8}$ years after the Big Bang via molecular hydrogen cooling in
${\sim}10^{5-6}~{}M_{\odot}$ dark matter (DM) “minihalos” (Haiman et al.,
1996b; Tegmark et al., 1997; Abel et al., 2002; Yoshida et al., 2003; O’Shea &
Norman, 2007; Wise & Abel, 2007; Greif, 2015; Kulkarni et al., 2021; Schauer
et al., 2021; Nebrin et al., 2023; Hegde & Furlanetto, 2023). Due to less
efficient gas cooling during star formation, Pop III stars are expected to
have higher masses than Pop II/I, perhaps with typical values of
${\sim}10-1000~{}M_{\odot}$ (Greif et al., 2011; Hosokawa et al., 2011; Hirano
et al., 2014; Susa et al., 2014; Ishigaki et al., 2018), but the precise
initial mass function (IMF) remains uncertain.

We note that because the very first stars are predicted to form in small
numbers at very high redshifts ($z{\gtrsim}20$), they are unlikely to be
directly observed even with the _James Webb Space Telescope_ (_JWST_).
However, most models predict that Pop III star formation should continue well
into or beyond the epoch of reionization (Riaz et al., 2022; Borrow et al.,
2023) where they might be detected either as purely Pop III galaxies or in
mixtures of Pop III and metal-enriched star formation (e.g., in the case of
inefficient metal-mixing). In fact, there are already some exciting but
tentative/unconfirmed Pop III signatures from _JWST_ (Wang et al., 2022;
Maiolino et al., 2023). Other promising observational probes of Pop III stars
include He ii line intensity mapping (Visbal et al., 2015; Parsons et al.,
2022), pair-instability supernovae detection (Whalen et al., 2013a; Hartwig et
al., 2018), 21cm cosmology (Visbal et al., 2012; Fialkov et al., 2013; Cohen
et al., 2017; Bowman et al., 2018), and stellar archaeology (Frebel & Norris,
2015).

Due to the hierarchical character of large-scale structure formation in
$\Lambda$CDM, the minihalos hosting Pop III stars serve as building blocks for
the first metal-enriched galaxies. Thus, the properties of high-redshift low-
mass galaxies observable with _JWST_ are expected to depend on the details of
Pop III star formation. For example, Abe et al. (2021) found that a top-heavy
Pop III initial mass function (IMF), which results in more highly energetic
pair-instability supernovae (PISNe), results in suppressed star formation in
galaxies hosted by ${\sim}10^{9}~{}M_{\odot}$ halos. Understanding the first
metal-enriched galaxies is also interesting because they may drive the process
of cosmic reionization, they represent the earliest stages of more evolved
galaxies, and the existence of low-mass dark matter host halos, puts
constraints on alternative dark matter models, such as “fuzzy” or warm dark
matter (Gao & Theuns, 2007; Mocz et al., 2019, 2020; Kulkarni et al., 2022),
that suppress small-scale structure. Using _JWST_ we are able to study these
galaxies individually and determine their global properties through quantities
such as the UV luminosity function (LFs) (Bouwens et al., 2023; Atek et al.,
2023).

Theoretical models are required to interpret and guide the observations of the
first stars and galaxies described above. A primary tool which has been used
is hydrodynamical cosmological simulations. State-of-the-art simulations can
include a variety of physical processes such as gas chemistry/cooling,
radiative transfer of UV photons, and feedback from supernovae explosion.
These effects can be modeled with high fidelity, but have the drawback of
being numerically expensive. For instance, the _Renaissance Simulations_ ,
described in detail below, used around 10-million CPU hours on the NCSA Blue
Waters supercomputers to run a box size of 135 Mpc3 to a redshift of $z\sim
15$ (Xu et al., 2014; Chen et al., 2014).

This presents a serious challenge because in order to make predictions for the
observations described above one needs much larger volumes than can be modeled
with hydrodynamical simulations. For example, in order to understand the
impact of Pop III stars on low mass galaxies, one needs a large-statistical
sample of galaxies forming in ${\sim}10^{9}~{}M_{\odot}$ DM halos, but also
needs to resolve the low-mass halos that host Pop III star formation. Despite
its computational expense, the _Renaissance_ volume that can be simulated only
contains a few thousand of these objects, much less than required for robust
statistical analysis.

A different technique is to use computationally efficient semi-analytic models
which combine analytic prescriptions with dark matter halo merger trees
generated from $N$-body simulations or Monte Carlo methods. Semi-analytic
simulations have been used extensively to study the first stars and galaxies
(e.g., Trenti & Stiavelli, 2009; Crosby et al., 2013; Visbal et al., 2018;
Mebane et al., 2018; Yung et al., 2019; Liu & Bromm, 2020). These simulations
are able to rapidly scan the relevant parameter space (e.g., the Pop III IMF
and critical metallicity for Pop III star formation). However, the main
drawback of this type of modelling is that a large number of free parameters
are often employed and it can be unclear to what extent analytic assumptions
are justified.

In this paper, we explore using the _Renaissance_ simulations to perform a
detailed calibration of a semi-analytic model of the first stars (Visbal et
al., 2020). The goal of the exercise is to obtain a model that is both
computationally efficient and well-matched to detailed numerical simulations.
This enables extending the simulation results to much larger volumes such that
observational predictions become possible. To achieve this goal, we develop an
empirically calibrated star formation model based the the _Renaissance_
results. We include both Pop III and metal-enriched star formation, but the
emphasis is on metal-enriched since the halos hosting early Pop III star
formation are not well resolved. Our empirical calibration includes a number
of parameters such as metal-enriched starburst mass, quiescent period
following a burst, etc., that are measured directly from the simulations. This
model and parameterization can easily be applied to other semi-analytic models
or analytic calculations in the future.

This study is related to that of McCaffrey et al. (2023) which calibrates a
semi-analytic model to _Renaissance_. We note that in this case, only the
largest galaxy in _Renaissance_ was considered and the simulation datasets
used were unable to resolve the formation of Pop III. As explained in detail
below, we find that while our prescription is calibrated to a region at the
mean density of the Universe, it can reproduce the star formation rate density
typically within a factor of $\sim$2 without any calibration on such a region.
Thus, we have shown that a calibration framework such as that developed in
this paper can extend the reach of numerical simulations both in volume and
range of cosmic time.

This paper is structured as follows. In Section 2, we review the main features
of the _Renaissance_ simulation suite. In Section 3, we describe our
parameterization of metal-enriched and Pop III star formation in _Renaissance_
along with our sample. In Section 4, we present the calibration of a semi-
analytic model to _Renaissance_ along with our results for metal-enriched and
Pop III star formation. Finally, we summarize the results and discuss our main
conclusions in Section 5. Throughout we assume a cosmology consistent with
_Renaissance_ using parameters from the WMAP7 $\Lambda$CDM+SZ+LENS best fit
(Komatsu et al., 2011) with $\Omega_{\mathrm{M}}=0.266$,
$\Omega_{\mathrm{\Lambda}}=0.734$, $\Omega_{\mathrm{b}}=0.0449$, $h=0.71$,
$\sigma_{\mathrm{8}}=0.81$, $n=0.963$.

##  2 The Renaissance Simulations

The _Renaissance Simulations_ are cosmological hydrodynamics simulations that
were performed with the adaptive mesh refinement code enzo (Bryan et al.,
2014; Brummel-Smith et al., 2019) to simulate the formation of the first stars
and galaxies. Radiative transfer from ionizing radiation is calculated with
the moray package (Wise & Abel, 2011) along with chemistry and cooling for
nine species of hydrogen and helium (Abel et al., 1997) and metallicity
dependent cooling tables (Smith et al., 2009).

_Renaissance_ consists of three different regions taken from a comoving volume
of $(40\ \mathrm{Mpc})^{3}$. The three regions named “RarePeak”, “Normal”, and
“Void” are in overdense, average, or underdense regions and evolved to
redshifts of $z=15,\ 12.5,\ \mathrm{and}\ 8$. RarePeak has a comoving volume
of $(3.8\times 5.4\times 6.6)\ \mathrm{Mpc}^{3}$ and Normal/Void have comoving
volumes of $(6.0\times 6.0\times 6.125)\ \mathrm{Mpc}^{3}$, respectively
(O’Shea et al., 2015). Our semi-analytic model calibration is focused on the
Normal region and then applied to RarePeak and Void for comparison.

It is important to emphasize key aspects of the metal-enriched and Population
III star formation prescriptions in _Renaissance_ due to their importance in
our star formation calibration. The metal-enriched star clusters are formed
from molecular clouds with a typical radius of $\sim 6$ pc. The clusters have
a total stellar mass populated with stars assuming a Salpeter IMF and emit
6000 hydrogen ionizing photons per stellar baryon over 20 Myr (the maximum
lifetime of an OB star). After 4 Myr, the clusters generate $6.8\times
10^{48}\ \mathrm{erg}\ \mathrm{M_{\mathrm{\odot}}^{-1}}$ for SNe feedback. For
Pop III, individual stars are formed before the metal-enriched star clusters
with a mass randomly sampled from a power-law initial mass function similar to
a Salpeter IMF with an exponential decrease below a characteristic mass of
$M_{\mathrm{char}}=40M_{\mathrm{\odot}}$. Stellar lifetimes along with Lyman-
Werner and hydrogen ionizing photon luminosities are drawn from Schaerer
(2002). Stars in the mass ranges of $11-20\ M_{\mathrm{\odot}}$, $20-40\
M_{\mathrm{\odot}}$, and $140-260\ M_{\mathrm{\odot}}$ explode as Type II
supernovae and hypernovae (Nomoto et al., 2006), and pair-instability
supernovae (Heger et al., 2003), respectively with energies exceeding
$10^{51}\ \mathrm{erg}$. Stars with masses outside these ranges directly
collapse into a black hole without an explosion. Further details can be found
in Wise et al. (2012a).

High-redshift galaxy observations from the CEERS (Finkelstein et al., 2023)
and JADES surveys (Eisenstein et al., 2023) imply that star formation rates
are likely higher than predicted by previous simulations (Arrabal Haro et al.,
2023; Bunker et al., 2023). McCaffrey et al. (2023) recently demonstrated that
early _Renaissance_ galaxies typically experience rapid star formation, with a
specific star formation rate of $\sim 10^{-8}\ \mathrm{yr}^{-1}$ which is
consistent with the galaxies observed by CEERS and JADES. This agreement
suggests that _Renaissance_ is a good choice for calibration to high-$z$
galaxy formation.

##  3 Calibration to Renaissance

###  3.1 Overview

We examined star formation histories for a wide selection of _Renaissance_
halos, building up a representative sample detailed in Section 3.2. In Figure
1, we show the metal-enriched star formation history for the most massive halo
in the Normal region. Stellar evolution in this halo is representative of
trends typically found in other _Renaissance_ halos. We find after the first
generation of stars, star formation proceeds through two distinct stages, a
bursty stage characterized by rapid star formation followed by a steady stage
with stars forming at a more constant rate. The division of metal-enriched
star formation into bursty and steady stages is the foundation for our
calibration to _Renaissance_ and later implementation in our semi-analytic
model.

Figure 1: Metal-enriched star formation history for a single _Renaissance_
halo. $M_{\mathrm{\star}}$ is the total mass of metal-enriched stars formed in
the halo. This halo experiences two episodes of bursty star formation followed
by a stage of steady star formation once the halo reaches a virial temperature
threshold. The burst stellar mass, quiescent period, ionization feedback
threshold temperature, and steady star formation efficiency are labeled
$M_{\mathrm{burst}}$, $t_{\mathrm{quies}}$, $T_{\mathrm{vir,steady}}$, and
$f_{\mathrm{II}}$ respectively. Bursty and steady stages are found throughout
the sample of halos we analyze.

The bursty stage begins with the first instance of metal-enriched star
formation following the death of Pop III stars and dispersal of their metals
into the interstellar medium (ISM). We do not include external enrichment in
this analysis, discussed further in Section 4. Forming in low-mass halos at
high redshift, metal-enriched stars form in a rapid starburst at the onset of
the bursty stage. Supernovae feedback from stars formed in the burst disperses
any dense clumps and expels most of the gas from the shallow gravitational
potential well. Following the expulsion of gas, star formation enters an
extended quiescent period as gas falls back into the halo and dense clumps
begin to reform. A halo can experience multiple bursts and quiescent periods
dependent on its virial mass, see Section 3.3.1 for more detail.

After the initial bursty stage, a change in star formation occurs as the
steady stage begins. This transition corresponds to a fixed virial temperature
which is consistent with the halo mass below which reionization has been shown
to suppress star formation as a result of gas photoheating (e.g., Shapiro et
al., 1994; Thoul & Weinberg, 1996; Gnedin & Hui, 1998; Gnedin, 2000; Dijkstra
et al., 2004; Hoeft et al., 2006; Okamoto et al., 2008; Sobacchi & Mesinger,
2013; Noh & McQuinn, 2014). Once the steady stage begins, the halo is
sufficiently massive to retain enough gas to fuel continuous future star
formation. Total metal-enriched star formation produced during the steady
stage eventually completely dominates at lower redshifts. More detail about
each component of the two stage bursty/steady model shown in Figure 1 is
covered in Section 3.3.

###  3.2 Sample

We quantified trends in star formation history by focusing on the evolution of
individual dark matter halos in _Renaissance_. Halos grow through a
combination of accretion and mergers with other halos, tracked using the halo
finder rockstar (Behroozi et al., 2013a) and merger trees generated with
consistent trees (Behroozi et al., 2013b). The main progenitor halos in these
trees represent the most massive progenitor at each step going back in time.
We preferentially select halos with smoother merger histories that typically
have uncomplicated star formation histories. We leave modeling of more violent
mergers for future work. Halos that experience large virial mass decreases in
their growth histories due to complicated mergers or other factors were
excluded for our sample.

To calibrate the semi-analytic model, we gathered a sample of 27 dark matter
halo and stellar mass histories illustrative of typical metal-enriched star
formation in _Renaissance_ halos. For each halo in our sample, we visually
identify when every starburst and quiescent period occurs along with the start
of steady star formation. In Figure 2, we show metal-enriched star formation
histories for seven different halos. A halo typically progresses through a
bursty stage with intermittent star formation followed by a steady stage, but
can sometimes bypass the bursty stage entirely if star formation begins after
the halo exceeds $T_{\mathrm{vir,steady}}$.

Figure 2: Metal-enriched star formation histories for seven _Renaissance_
halos in our sample. Halos can undergo one or multiple bursts and can also
skip the bursty stage, instead exhibiting star formation consistent with the
steady stage (green and brown).

###  3.3 Star Formation Parameterization

Our two-stage metal-enriched star formation model has four key physical
quantities. In the bursty stage this includes the mass formed in a starburst
$M_{\mathrm{burst}}$, and a quiescent period following a burst
$t_{\mathrm{quies}}$. The steady stage begins once the virial temperature
threshold, $T_{\mathrm{vir}}$, is exceeded, signalling the start of the steady
stage. Steady star formation then proceeds at a constant efficiency
$f_{\mathrm{II}}$.

We focus our calibration on three different aspects of Pop III formation in
_Renaissance_. When Pop III stars form inside halos, the total Pop III stellar
mass produced once formation conditions are met, and how much time must pass
after the formation of Pop III before the formation of metal-enriched stars.

Overall, we focus on connecting star formation to halo properties such as the
virial mass and do not generally consider environmental properties with the
exception of feedback from cosmic reionization. However, as we discuss below
in Section 4.2, environmental properties such as the large-scale overdensity
or different physics in massive halos could impact star formation. A summary
of all model parameters along with best fit values discussed in the following
sections can be found in Table 1.

####  3.3.1 Metal-Enriched Star Formation

For the bursty stage, we evaluated multi-parameter log scale power law fits
using least squares regression from the scipy curve-fit package for the burst
masses and quiescent periods in our sample. We also determine standard
deviations to account for the distributions of $M_{\mathrm{burst}}$ and
$t_{\mathrm{quies}}$.

In Figure 3, we find that the metal-enriched burst masses in our sample are
closely correlated with both the halo virial mass and the redshift of the
burst. We parameterize the stellar mass formed in a burst as

$M_{\mathrm{burst}}=\alpha\left(\frac{M_{\mathrm{vir}}}{10^{7}M_{\mathrm{\odot}}}\right)^{\gamma}\left(\frac{1+z}{20}\right)^{\delta}M_{\mathrm{\odot}}$
(1)

where $M_{\mathrm{vir}}$ and $z$ are respectively the halo virial mass and
redshift before a starburst. The best fit parameters are $\alpha=5.364\times
10^{3}$, $\gamma=1.092$, and $\delta=0.578$. The distribution of the sample
around $M_{\mathrm{burst}}$ is approximately log-normal with a standard
deviation of
log${}_{\mathrm{10}}\left(M_{\mathrm{burst}}/M_{\mathrm{\odot}}\right)=0.408$.
As halos reach higher masses at lower redshift, the stellar mass of metal-
enriched stars formed in a burst increases.

Figure 3: Metal-enriched starburst masses from the sample of _Renaissance_
halos during the bursty stage. Best fit two parameter power law (solid-blue)
along with 1$\sigma$ lines (dashed-cyan). Burst mass increases for larger mass
halos and at lower redshifts.

For the quiescent periods following a starburst, we find a strong negative
correlation with the halo virial mass, as shown in Figure 4. The best fit is

$t_{\mathrm{quies}}=\epsilon\left(\frac{M_{\mathrm{vir}}}{10^{7}M_{\mathrm{\odot}}}\right)^{\kappa}\text{Myr}$
(2)

where $M_{\mathrm{vir}}$ is the virial mass following a starburst and the best
fit parameters are $\epsilon=4.481$ and $\kappa=-0.394$. The distribution of
the sample around $t_{\mathrm{quies}}$ is approximately log-normal with a
standard deviation of
log${}_{\mathrm{10}}\left(t_{\mathrm{quies}}/\mathrm{Myr}\right)=0.255$. As
halos increase in mass, building up deeper gravitational potential wells
better able to hold onto gas, a typical quiescent period declines from
$\sim$100 Myr to $\sim$10 Myr for halo masses between $\sim 10^{7}$–$10^{8}\
M_{\mathrm{\odot}}$.

If shorter quiescent periods are related to the retention of gas reservoirs
that fuels subsequent star formation, then the rate the halo is accreting mass
should impact the length of the quiescent period. A fast growing halo could
potentially recover gas ejected from a starburst and reduce the quiescent
period before the next instance of star formation accordingly. We
parameterized this effect by considering the time for a halo to double in
virial mass following a starburst, $t_{\mathrm{double}}$ (see Figure 4).
Incorporating $t_{\mathrm{double}}$ in our power law did not significantly
change the overall fit, with the variance of the sample data from the fit
being reduced by less than 20 percent. Therefore, we do not include
$t_{\mathrm{double}}$ in the final fits.

Figure 4: Quiescent times from the sample of _Renaissance_ halos during the
bursty stage. Power law best fit (solid-orange) and 1$\sigma$ lines (dashed-
gold). The color bar values correspond to the time in millions of years it
takes the halo to double in virial mass. Overall the quiescent time between
star formation events decreases for higher mass halos.

With the start of the steady stage, we found that the transition occurs at
virial temperatures within a factor of 2 of
$T_{\mathrm{vir,steady}}=1.76\times 10^{4}$ K corresponding to an ionization
feedback mass for almost all star formation histories in our sample. We assume
a mean molecular weight of $\mu=0.6$ consistent with fully ionized primordial
gas for the temperature (Barkana & Loeb, 2001). Due to the small range in
temperature, we define the onset of the steady stage to occur once the virial
temperature exceeds $T_{\mathrm{vir,steady}}$. The virial temperature at the
beginning of the steady stage showed weak or no correlation with other halo
properties.

Once the steady stage is underway, we found metal-enriched star formation to
be well described using a system of differential equations presented in
Furlanetto & Mirocha (2022). These coupled equations derive the evolution of
both the gas and stellar mass of a halo using a simplified implementation of
the ‘bathtub’ model (e.g., Bouché et al., 2010; Davé et al., 2012; Dekel &
Mandelker, 2014):

$\dot{M}_{\mathrm{*}}=\frac{f_{\mathrm{II}}}{t_{\mathrm{ff}}}M_{\mathrm{gas}},$
(3)

$\dot{M}_{\mathrm{gas}}=\dot{M}_{\mathrm{acc}}-\dot{M}_{\mathrm{*}}-\eta_{\mathrm{SN}}\dot{M}_{\mathrm{*}},$
(4)

where $\dot{M}_{\mathrm{*}}$ is the stellar mass growth rate,
$\dot{M}_{\mathrm{gas}}$ is the rate of change in halo gas mass,
$\dot{M}_{\mathrm{acc}}$ is the halo gas accretion rate, $f_{\mathrm{II}}$ is
star formation efficiency, $t_{\mathrm{ff}}$ is the free fall time assumed to
be a tenth of the Hubble time (which corresponds to the dynamical time at the
viral radius of a dark matter halo described in Loeb & Furlanetto (2013)), and
$\eta_{\mathrm{SN}}$ is the gas ejection efficiency for supernovae feedback.

Two parameters that warrant further discussion are the SNe ejection
efficiency, $\eta_{\mathrm{SN}}$, and steady star formation efficiency,
$f_{\mathrm{II}}$. We determine $\eta_{\mathrm{SN}}$ using a similar
prescription to Sassano et al. (2021), with a SNe ejection efficiency of

$\eta_{\mathrm{SN}}=\beta\times\frac{2E_{\mathrm{SN}}\epsilon_{\mathrm{SN}}R_{\mathrm{SN}}}{v_{\mathrm{esc}}^{2}},$
(5)

where $E_{\mathrm{SN}}$ is the average SNe explosion energy,
$\epsilon_{\mathrm{SN}}$ is the SNe kinetic energy conversion fraction,
$R_{\mathrm{SN}}$ is the SNe rate per mass, $v_{\mathrm{esc}}$ is the escape
velocity from the halo, and $\beta$ is the SNe ejection scale factor.

We assume the average energy for $E_{\mathrm{SN}}$ is $10^{51}$ erg (Nomoto et
al., 2006). With most of the energy from a supernovae lost as thermal energy,
only a small fraction is converted into kinetic energy coupled to the gas
(Walch & Naab, 2015). We adopt a value of $\epsilon_{\mathrm{SN}}=1.6\times
10^{-3}$ for this efficiency that is calibrated to observations (Sassano et
al., 2021). For $R_{\mathrm{SN}}$, we use metal-enriched SNe feedback energy
per mass in agreement with _Renaissance_ of $6.8\times 10^{48}$ erg
$\mathrm{M}_{\mathrm{\odot}}^{-1}$ (Wise et al., 2012b). Dividing this
quantity by $E_{\mathrm{SN}}$ yields a SNe rate of $R_{\mathrm{SN}}=6.8\times
10^{-3}\ \mathrm{M}_{\mathrm{\odot}}^{-1}$. Finally, $v_{\mathrm{esc}}$, is
set to be the circular velocity of a halo excluding a factor of $\sqrt{2}$ in
the numerator of the escape velocity using the relation
$\sqrt{GM/r_{\mathrm{vir}}}$ from Barkana & Loeb (2001).

We lower the SNe efficiency using a SNe ejection scale factor $\beta$
calibrated to high-mass halos. Values between $\beta=0$–$0.1$ generate star
formation consistent with our sample and we find $\beta=0.01$ produces total
metal-enriched star formation that agrees well with _Renaissance_. For
$\beta=1.0$, $\eta_{\mathrm{SN}}$ significantly suppresses star formation in
the most massive halos, producing stellar masses around an order of magnitude
lower than _Renaissance_ due to the efficient ejection of a large fraction of
gas early in the steady stage.

To calibrate the steady star formation efficiency in Equation 3 to
_Renaissance_ , we determine values for $f_{\mathrm{II}}$ in agreement with
our sample of metal-enriched star formation histories using the following
method. First, we separate the steady stage from the bursty stage for each
star formation history in our sample. We set the initial stellar mass equal to
_Renaissance_ and assume the gas mass is consistent with the cosmic baryon
fraction of the halo virial mass. Next, we evolve Equations 3 and 4 over the
entire steady stage and perform a regression using the scipy curve-fit package
to determine the best fit for $f_{\mathrm{II}}$.

Figure 5 shows the distribution of best fit steady stage star formation
efficiencies. We find the distribution is roughly log-normal and centered
around a mean efficiency of $f_{\mathrm{II}}=0.013$ with a standard deviation
of log $\left(f_{\mathrm{II}}\right)=0.422$. The predicted metal-enriched
stellar mass using the best fit efficiency with Equations 3 and 4 at most
deviates from _Renaissance_ by a factor of $\sim$1.5. We note that outliers
such as the low efficiency of $f_{\mathrm{II}}=0.0006$ have a large impact on
the variance of the distribution due to the small sample size.

Figure 5: Metal-enriched star formation efficiencies for the steady stage in
_Renaissance_ halos. Each $f_{\mathrm{II}}$ were determined by finding the
best fit for Equations 3 and 4 during the steady stage of a halo in the
sample.

We found that a halo mass dependent $f_{\mathrm{II}}$ of the functional form
$f_{\mathrm{II}}M_{\mathrm{vir}}^{\tau}$ also produced good fits for the
_Renaissance_ star formation histories of most halos during the steady stage.
However for the largest halos at low redshift, this produces unphysically high
star formation efficiencies and therefore we do not include it in our model.

####  3.3.2 Population III Star Formation

Population III stars form in pristine metal-free halos that exceed the minimum
mass where dense gas clouds can condense due to cooling from roto-vibrational
molecular hydrogen transitions (Haiman et al., 1996a, b; Tegmark et al., 1997;
Machacek et al., 2001). The $\mathrm{H_{2}}$ fraction is fundamental to
determine if a halo forms Pop III stars and is included as a formation
condition in _Renaissance_ accordingly (see Section 2). Therefore, physical
processes that change the $\mathrm{H_{2}}$ fraction are important to consider.

Lyman-Werner photons with energies between 11.2-13.6 eV emitted by the first
stars can efficiently photodissociate $\mathrm{H_{2}}$ through the two-step
Solomon process (Stecher & Williams, 1967). The LW photons travel over
cosmological scales to neighboring metal-free halos, suppressing
$\mathrm{H_{2}}$ cooling and the formation of Pop III stars (Dekel & Rees,
1987; Haiman et al., 1997, 2000; Machacek et al., 2001; Wise & Abel, 2007;
O’Shea & Norman, 2008; Ahn et al., 2009; Visbal et al., 2014). We follow
Fialkov et al. (2013) to compute the minimum mass for molecular cooling that
can form stars in minihalos using the relation,

$M_{\mathrm{H_{2}}}=2.5\times
10^{5}\left(\frac{1+z}{26}\right)^{-3/2}\times\left[1+6.96(4\pi
J_{\mathrm{LW,21}})^{0.47}\right]M_{\mathrm{\odot}}$ (6)

where the first term is the minimum halo mass in the absence of a LW
background to suppress $\mathrm{H_{2}}$ determined using the same assumptions
as Visbal et al. (2020) and $J_{\mathrm{LW,21}}$ is the local LW background in
units of $10^{-21}\ \mathrm{erg}\ \mathrm{s}^{-1}\ \mathrm{cm}^{-2}\
\mathrm{Hz}^{-1}\ \mathrm{sr}^{-1}$. We use the background presented in Xu et
al. (2016) that self-consistently calculates a LW background from star
formation in the Normal region.

Once a halo reaches virial temperatures of $T_{\mathrm{vir,steady}}\gtrsim
10^{4}\ \mathrm{K}$, collisionally excited line emission from atomic hydrogen
becomes the dominant cooling process. For the atomic cooling mass threshold we
adopt the relation,

$M_{\mathrm{a}}=5.4\times
10^{7}\left(\frac{1+z}{11}\right)^{-3/2}M_{\mathrm{\odot}}$ (7)

which is determined using the cosmological hydrodynamical simulations of
Fernandez et al. (2014).

The _Renaissance_ simulations have a dark matter mass resolution of $2.9\times
10^{4}M_{\mathrm{\odot}}$ and the halo mass function is well-resolved for
halos with masses $\gtrsim 2\times 10^{6}M_{\mathrm{\odot}}$ (O’Shea et al.,
2015). This artificially sets a lower-limit on the minimum halo mass scale for
Pop III star formation. Simulations by Kulkarni et al. (2021) and Schauer et
al. (2021) find minimum masses of order $10^{5}$-$10^{6}M_{\mathrm{\odot}}$ at
high redshift for similar LW backgrounds which are unable to be resolved by
_Renaissance_. Therefore, we adopt a mass resolution limit of $7\times
10^{6}M_{\mathrm{\odot}}$ below which Pop III is assumed to be unable to form.
This allows our model to closely reproduce Pop III star formation by delaying
star formation in unresolved low mass halos similarly to _Renaissance_.

The left panel of Figure 6 shows the dark matter halo masses of _Renaissance_
halos that host the first Population III star formation,
$M_{\mathrm{Halo,PopIII}}$, as a function of redshift. We include
$M_{\mathrm{a}}$ and $M_{\mathrm{H_{2}}}$ using a $J_{\mathrm{LW,21}}$
equivalent to the global LW background (a larger local value for
$J_{\mathrm{LW,21}}$ than the background would result in a higher
$M_{\mathrm{H_{2}}}$).

To produce Pop III stellar masses consistent with _Renaissance_ , we
determined the total mass of Pop III stars formed during the initial instance
of star formation in _Renaissance_ halos. We found that the total Pop III
mass, $M_{\mathrm{Total,PopIII}}$, closely follows a log-normal distribution
centered on a mass of $\approx 150M_{\mathrm{\odot}}$ with a standard
deviation of
log${}_{\mathrm{10}}\left(M_{\mathrm{Total,PopIII}}/M_{\mathrm{\odot}}\right)=0.371$.
Both the virial mass of the halo and the redshift were found to have minimal
correlation with total Pop III mass formed.

Finally, we calculated the typical delay between the first Pop III supernovae
in a halo, occurring after a mass-dependent stellar lifetime from Schaerer
(2002), with the beginning of metal-enriched star formation. Figure 7 shows
the distribution of the delay times as a function of halo mass. We find three
distinct populations of delay times. One population has metal-enriched star
formation starting rapidly between 0.01–0.1 Myr after the first Pop III
supernova. Another population has delays between 0.1–3 Myr and a final smaller
group has longer delay times of 10–100 Myr. We investigated deeper potential
wells in high mass halos retaining metal-enriched gas, ionizing feedback
suppressing star formation, and powerful Pop III supernovae as physical
origins for these three populations. However, the delay times show little to
no correlation with halo mass, virial temperature, and total supernovae energy
emitted by Pop III stars before the onset of metal-enriched formation.

The halos that begin to form metal-enriched stars after less than a million
years are of particular interest due to the short timescales for metals to be
dispersed into the nearby ISM. In this case, star formation is triggered in
nearby clumps that have survived the radiative phase of the Pop III stars. A
SN blastwave will compress the clumps above the density threshold for star
formation, however, some of these cases could be artificial. It is possible a
number of these clumps meet the criteria for star formation in the simulation,
but are not gravitationally bound and should instead dissipate sometime later.
The typical _Renaissance_ output cadence of $\sim 4$ Myr hinders further
investigation into these clumps and the underlying cause of the quick
transition to metal-enriched star formation. To be as consistent with
_Renaissance_ as possible, we consider a $t_{\mathrm{delay}}$ ranging between
0.01–100 Myr with the understanding that numerical issues could be responsible
for the shortest delays.

Figure 6: Dark matter masses of halos that host the first Population III star
formation for Renaissance (left panel) and one realization of the semi-
analytic model (right panel). We show the median halo mass hosting Pop III for
Renaissance (magenta) and the semi-analytic model (green). The dashed and dot-
dashed curves show $M_{\mathrm{a}}$ and $M_{\mathrm{H_{2}}}$, respectively.
The sharp cut-off in halo mass for the semi-analytic model in the right panel
is due to the mass resolution limit of $7\times 10^{6}M_{\mathrm{\odot}}$.
Figure 7: The delay period before the formation of metal-enriched stars
following the death of the first Pop III star compared with the dark mass halo
mass at formation of Pop III.  Table 1: Parameters for Star Formation in the
Semi-analytic Model

Parameter | Description | Fit Values | Standard Deviation
---|---|---|---
$M_{\mathrm{burst}}$ | Metal-enriched starburst mass | $\alpha=5.364\times 10^{3}$ | $0.408$
|  | $\gamma=1.092$ |
|  | $\delta=0.578$ |
$t_{\mathrm{quies}}$ | Quiescent period following metal-enriched starburst | $\epsilon=4.481$ | $0.255$
|  | $\kappa=-0.394$ |
$f_{\mathrm{II}}$ | Star formation efficiency during steady stage | $0.013$ | $0.422$
$M_{\mathrm{Total,PopIII}}$ | Total Pop III Stellar Mass formed after halo exceeds $M_{\mathrm{min}}$ |  $150$ $M_{\mathrm{\odot}}$ | $0.371$
$t_{\mathrm{delay}}$ | Delay before metal-enriched star formation |  $0.01-100$ Myr | …

Note. — $M_{\mathrm{burst}}$ is given by Equation 1. $t_{\mathrm{quies}}$ is
given by Equation 2. $f_{\mathrm{II}}$ is given by the best fit metal-enriched
star formation efficiencies in Figure 5. $M_{\mathrm{Total,PopIII}}$ is
described in Section 3.3.2. $M_{\mathrm{burst}}$, $t_{\mathrm{quies}}$,
$f_{\mathrm{II}}$, and $M_{\mathrm{Total,PopIII}}$ each follow a log-normal
distribution. The last column shows the standard deviation of the log10 of the
corresponding quantity as mentioned in the text. $t_{\mathrm{delay}}$ is
described in Section 3.3.2 and follows a log-uniform distribution.

##  4 Model Modifications and Results

We modify the semi-analytic model introduced in Visbal et al. (2020) which is
based on dark matter halo merger trees from cosmological N-body simulations,
and incorporates star formation at a constant efficiency, feedback from Lyman-
Werner radiation and growth of H ii regions in the intergalactic medium. These
processes are calculated on a grid using fast-fourier transforms (FFTs) using
3D spatial positions and clustering of halos.

Star formation in the updated semi-analytic model consists of Pop III
formation followed by two stages of metal-enriched star formation with a
bursty and steady stage. The calibration to _Renaissance_ is implemented by
including a minimum mass for Pop III star formation, $M_{\mathrm{min}}$, along
with total Pop III stellar masses and delay periods before the beginning of
metal-enriched star formation that are consistent with _Renaissance_.
Additionally, for metal-enriched star formation we introduce our findings for
$M_{\mathrm{burst}}$, $t_{\mathrm{quies}}$, $T_{\mathrm{vir,steady}}$, and
$f_{\mathrm{II}}$ from the bursty and steady stages previously presented in
Section 3.3.1. See Table 1 for a concise overview of the implementation into
the model. Finally, we do not include external enrichment in our calibrations.
External enrichment of neighboring star-forming halos was found to marginally
increase the metal-enriched SFRD and slightly reduce the Population III SFRD
in Visbal et al. (2020). While this does not preclude external enrichment
being more important in dense, highly-clustered environments, for now we leave
calibration to external enrichment in _Renaissance_ for future work.

Now we will cover the modifications to the Visbal et al. (2020) model in more
detail. We assume the minimum mass for Pop III star formation to be the
minimum of $M_{\mathrm{H_{2}}}$ and $M_{\mathrm{a}}$ with a mass resolution
lower-limit of $\gtrsim 7\times 10^{6}M_{\mathrm{\odot}}$ from Section 3.3.2.
For the local LW background $J_{\mathrm{LW,21}}$ term in $M_{\mathrm{H_{2}}}$,
an external background is included in addition to the local component from Pop
III and metal-enriched stellar mass. We implement a background that self-
consistently calculates a LW background from star formation in the Normal
region discussed in Section 3.3.2. If a halo with no prior star formation is
located in an ionized region, we increase the minimum mass to
$M_{\mathrm{min}}=1.5\times
10^{8}\left(\frac{1+z}{11}\right)^{-3/2}M_{\mathrm{\odot}}$ (in agreement with
Dijkstra et al. (2004)).

The right panel of Figure 6 shows $M_{\mathrm{Halo,PopIII}}$ predicted by our
model. There are two differences between our model and _Renaissance_ that are
important to highlight. First, values of $M_{\mathrm{Halo,PopIII}}$ in our
model are much more tightly clustered at lower halo masses than _Renaissance_.
A potential consequence of this behavior is the same halo in our model
sometimes can form Pop III at a lower halo mass corresponding to an earlier
time. Second, $M_{\mathrm{Halo,PopIII}}$ in _Renaissance_ is typically smaller
for $z>17$ and higher for $z<17$ compared to the semi-analytic model.
Additionally, $M_{\mathrm{Halo,PopIII}}$ displays greater variability for
_Renaissance_. The consequences of these differences will be further discussed
in Section 4.2.

Once the mass of a halo exceeds $M_{\mathrm{min}}$, we use the findings for
Pop III formation discussed in Section 3.3.2. First, we randomly sample a
total Pop III stellar mass from the log-normal _Renaissance_ distribution
centered on $\approx 150\ M_{\mathrm{\odot}}$. Then, we assign a delay period
before metal-enriched star formation can begin, $t_{\mathrm{delay}}$, that
ranges between 0.01–100 Myr.

After $t_{\mathrm{delay}}$ following Pop III star formation has elapsed, we
determine if a halo is the bursty or steady stage by comparing its virial
temperature to the $T_{\mathrm{vir,steady}}=1.76\times 10^{4}$ K threshold. If
in the bursty stage, we determine the total stellar mass formed in a starburst
by randomly sampling from a log-normal distribution centered on
$M_{\mathrm{burst}}$ from Equation 1. A quiescent period sampled from another
log-normal distribution centered on $t_{\mathrm{quies}}$ from Equation 2 is
assigned to the halo and its descendants. The halo will then either experience
a subsequent starburst after the quiescent period has passed or it will exceed
the $T_{\mathrm{vir,steady}}$ threshold and enter the steady star formation
stage.

At the beginning of the steady stage, we initially assume the halo contains
gas consistent with the cosmic baryon fraction of its virial mass. A star
formation efficiency, $f_{\mathrm{II}}$, is then randomly sampled from the
roughly log-normal distribution in Figure 5. This efficiency will be passed to
all descendants of this halo unless a merger with other halos in the steady
stage occurs. Then, the average of progenitor halo star formation efficiencies
is assigned to the descendant instead. We calculate metal-enriched star
formation and follow the subsequent evolution of gas mass using Equations 3
and 4 respectively.

###  4.1 Individual Star Formation Histories

To determine the effectiveness of our calibration, we compare different
realizations of our model with _Renaissance_ for metal-enriched halo star
formation histories. We begin by running 10 different random seeds of our
semi-analytic model on the _Renaissance_ dark matter halo merger trees. Figure
8 shows a selection of metal-enriched SFHs from our simulations. The top-left
panel shows the metal-enriched SFH for the most massive halo in the
_Renaissance_ Normal region. The large second starburst and high
$f_{\mathrm{II}}\approx 0.1$ are extreme values in our distributions for
$M_{\mathrm{burst}}$ and $f_{\mathrm{II}}$, making them correspondingly
difficult to reproduce in our model. The SFHs in the top-right and bottom-
right panels transition from the bursty to the steady stage differently from
_Renaissance_. Our fixed threshold for the transition at
$T_{\mathrm{vir,steady}}=1.76\times 10^{4}$ K can lead to more starbursts or
an earlier transition to the steady stage, by the final simulation redshift,
the total star formation produced during the steady stage lies within the
scatter produced by the model. The bottom-left panel shows a SFH with bursty
and steady stages that fall well within the scatter of our model.

Generally, there is good agreement between star formation histories produced
by the semi-analytic model and _Renaissance_. By the final simulation redshift
when metal-enriched stars dominate, star formation for _Renaissance_ halos
typically lie within the range predicted by our model and agreement for the
best model is within a factor of 2. Additionally, the chaotic bursty stage is
also well reproduced if metal-enriched star formation begins at a similar
time.

While the overall agreement with _Renaissance_ is considerable, the model does
not perfectly capture the variance in individual star formation histories. Our
choice to calibrate to a sample of main progenitor halos without violent major
mergers in _Renaissance_ could potentially not be fully representative of all
halos. Additionally, the differential equations used during the steady stage
produce star formation that is smoother than _Renaissance_ , lacking
variations in star formation that often occur. Regardless, the simplified
prescription for metal-enriched star formation is often able to reproduce SFHs
in the mass range of $10^{6}$-$10^{8}$ $M_{\mathrm{\odot}}$ remarkably well.

Figure 8: Total metal-enriched stellar masses and redshifts of four different
main progenitor halos. The _Renaissance_ metal-enriched star formation history
(black) is compared with ten different realizations of the semi-analytic model
(dashed).

###  4.2 Total Star Formation Rate

Extending from the individual star formation histories discussed in the
previous section, we also compare star formation rate densities averaged
throughout the simulation box to evaluate how well 10 different realizations
of our model reproduce star formation across entire _Renaissance_ volumes.

First, we will focus on metal-enriched star formation shown in the top row of
Figure 9. In the top-left panel of Figure 9, we find that our model produces a
metal-enriched SFRD consistent with the _Renaissance_ Normal region within a
factor of $\sim$2\. The good agreement with the _Renaissance_ SFRD indicates
our two-stage metal-enriched star formation model, when aggregated for many
halo SFHs, is able to closely reproduce total star formation in a large
volume.

In addition to our simulations using the calibrations to the Normal region, we
also consider the Void and RarePeak regions not included in our sample of
halos (mentioned in Section 2). To ensure consistency with _Renaissance_ , we
introduce a different LW background produced by Pop III stars used in Wise et
al. (2012b) for the RarePeak region. For the Void region, we use the same LW
background as the Normal region (see Section 3.3.2).

Results for the Void region are shown in the middle column of Figure 9. The
metal-enriched SFRD of the model is within a factor of $\sim$2 of
_Renaissance_ for early times and at most differs by a factor of $\sim$3\. For
$z<16$, the steady stage begins in an increasing number of halos and will
begin to dominate the metal-enriched SFRD. The scatter in the metal-enriched
SFRD produced by our model begins to narrow as a consequence of the smooth
SFHs produced during the steady stage mentioned in Section 4.1. This results
in an SFRD that exhibits much less variation than _Renaissance_ and could
explain the difference with the semi-analytic model at $z<16$.

In the right column of Figure 9 we show our model results for the RarePeak
region. Metal-enriched star formation in _Renaissance_ agrees within a factor
of $\sim$2 for $z>19$ but eventually is more efficient than our model by a
factor of $\sim$3 by $z\sim 15$. The source of this discrepancy is currently
unknown, but we speculate a potential explanation is some physics missing from
our prescription. For example, it has been shown that accretion driven mergers
can boost star formation by increasing torques and driving gas towards the
center of a halo (Hopkins et al., 2013). This could explain the discrepancies
we find in both Void and RarePeak.

Next, we will discuss Pop III star formation. Again, it is important to note
that _Renaissance_ has difficulties in resolving halos with masses $\lesssim
7\times 10^{6}M_{\mathrm{\odot}}$. Therefore, we focused mainly on calibration
without going into the same level of detail as we do for metal-enriched star
formation. As shown in the bottom row of Figure 9, the Pop III SFRDs produced
by the semi-analytic model for Normal, Void, and RarePeak demonstrate similar
behavior. The Pop III SFRD first increases, then, rests at the same value
before gradually decreasing. This is due to small discrepancies in the
$M_{\mathrm{min}}$ prescription in the semi-analytic model compared to
_Renaissance_.

We conclude that our model is able to reproduce the _Renaissance_ SFRDs within
a factor of $\sim$2 for the Normal region. For the Void and RarePeak regions
not used during the calibration of our model, the agreement is similar to the
Normal at high-$z$ and varies at low-$z$ to a factor of $\sim$2 and $\sim$3
for Void and RarePeak, respectively. In future work, incorporating additional
physics dependent on environmental effects could improve the agreement between
our model and _Renaissance_.

Figure 9: Cosmic SFRDs for both metal-enriched (orange top-row) Population III
(blue bottom-row) compared to _Renaissance_ (black). Results for the Normal,
Void and RarePeak volumes are shown in the left-, middle-, and right-columns
respectively. Overall, our model produces metal-enriched SFRDs typically
within a factor of $\sim$2 of _Renaissance_ for a variety of environments. The
larger discrepancy found at lower redshift in Void and RarePeak will be
investigated in future work.

###  4.3 Comparison with Updated Prescription for Population III Stars

Another question interesting to consider is the impact of incorporating recent
findings updating the minimum halo mass of Population III star formation,
$M_{\mathrm{min}}$, with new physics. In Kulkarni et al. (2021),
$M_{\mathrm{min}}$ is obtained as a function of LW background (including
$H_{\mathrm{2}}$ self-shielding), redshift, and baryon-dark matter streaming.
While we were able to calibrate to the _Renaissance_ Pop III SFRD within a
factor of $\sim$2, it is important to emphasize again that _Renaissance_ is
limited in resolving minihalos that form Pop III.

When updating our prescription for $M_{\mathrm{min}}$ to be consistent with
Kulkarni et al. (2021) and removing the mass resolution limit, we find in
Figure 10 a Pop III SFRD approximately an order of magnitude higher than the
prescription calibrated to _Renaissance_. This is due to the inclusion of
$H_{\mathrm{2}}$ self-shielding which lowers the impact of
$J_{\mathrm{LW,21}}$ on $M_{\mathrm{min}}$. The impact on metal-enriched star
formation due to the increased amount of Pop III was marginal. The metal-
enriched SFRD from our semi-analytic prescription both with and without the
updated $M_{\mathrm{min}}$ along _Renaissance_ agree within a factor of
$\sim$2 for $z<18$. This result can be heavily influenced by any changes to LW
feedback and reionization along with localized metal enrichment of halos not
included in our model. Determining the relative impact of these processes
after further calibration to _Renaissance_ will be the goal of future work.
Our method will allow this work to be done with minimal computational expense
and is a key benefit of semi-analytic models.

Figure 10: Cosmic Metal-enriched and Population III SFRDs compared with
_Renaissance_ for different $M_{\mathrm{min}}$ prescriptions in our model. The
$M_{\mathrm{min}}$ with the mass resolution limit used in this paper (dashed)
produces a Pop III SFRD that is roughly an order of magnitude lower than the
Pop III SFRD using $M_{\mathrm{min}}$ from Kulkarni et al. (2021). The metal-
enriched SFRDs agree within a factor of $\sim$2 except for $z>18$.

###  4.4 H ii Regions

To determine the size and growth of ionized H ii regions we apply the same FFT
method as Visbal et al. (2020). Using a $256^{3}$ cubic grid, we find the
total number density of ionizing photons produced by dark matter halos that
have escaped in the surrounding IGM in each cell. This is calculated as
$P_{\mathrm{cell}}=f_{\mathrm{esc,II}}\eta_{\mathrm{II}}M_{\mathrm{\star,II,cell}}/(m_{\mathrm{proton}}V_{\mathrm{cell}})+f_{\mathrm{esc,III}}\eta_{\mathrm{III}}M_{\mathrm{\star,III,cell}}/(m_{\mathrm{proton}}V_{\mathrm{cell}})$,
where $f_{\mathrm{esc,II}}$ and $f_{\mathrm{esc,III}}$ are the escape
fractions of hydrogen ionizing photons from halos in a cell containing metal-
enriched and Pop III stars. The number of hydrogen ionizing photons produced
per baryon for metal-enriched and Population III stars is assumed to be
$\eta_{\mathrm{II}}=4,000$ and $\eta_{\mathrm{III}}=65,000$ respectively. The
metal-enriched value corresponds to a Salpeter IMF from $0.1$ to
$100M_{\mathrm{\odot}}$ and metallicity $Z=0.0004$ (see Table 1 in Samui et
al. (2007)). Note, if values were used for $\eta$ that correspond directly
with _Renaissance_ , we would expect $f_{\mathrm{esc}}$ to change accordingly
but the overall ionization fraction and structure would remain the same. The
Pop III value is expected to emitted from a $\sim 40M_{\mathrm{\odot}}$ star
over its lifetime (Schaerer, 2002). Both $M_{\mathrm{\star,II,cell}}$ and
$M_{\mathrm{\star,III,cell}}$ are the total stellar masses ever formed by
halos in a cell. Finally, $V_{\mathrm{cell}}$ is the comoving volume of a
single cell. For each cell, $P_{\mathrm{cell}}$ is then smoothed over a range
of spherical bubble sizes and is selected as the center of an H ii bubble if
it exceeds a threshold of $i_{\mathrm{thres}}=4\bar{n}_{\mathrm{H}}$ (Visbal
et al., 2020).

The extent of ionized H ii regions determines when halos first form stars
through the minimum halo mass for Pop III star formation and if a halo’s
descendants ever form metal-enriched stars in our model with no external
enrichment. The top panel of Figure 11 compares the total ionization fraction
in our model to _Renaissance_. We calibrated the escape fraction of ionizing
radiation to $f_{\mathrm{esc},\mathrm{II}}=0.02$ for metal-enriched stars and
$f_{\mathrm{esc},\mathrm{III}}=0.1$ for Population III stars to produce a
total fraction of the IGM ionized by our model within a factor of 2 of
_Renaissance_ by $z\sim 11$.

Figure 11: Top panel: mean fraction of the volume of the IGM that is ionized
in _Renaissance_ (black) compared with the the semi-analytic model (green).
Projections of the H ii fraction for a region with a width and depth of
$\approx$2.7 and $\approx$0.33 comoving Mpc, respectively. Middle panel:
_Renaissance_. Bottom panel: semi-analytic model.

In the middle and bottom panels of Figure 11, we compare ionized bubble sizes
in our model and _Renaissance_ for halos in the same thin volume at $z=11.6$.
The difference in the H ii regions produced by our model in this region
demonstrates the limitations of our current approach in reproducing ionization
in simulations like _Renaissance_. We find that our model typically produces
larger ionized bubbles than _Renaissance_. These H ii regions extend from
massive halos with substantial star formation and are generally found at the
center of clustered environments. To address this discrepancy, we are
currently improving the prescription to produce more realistic bubbles and
will be the focus of future work (Behling et al., in prep).

##  5 Summary and Discussion

In this paper we explored calibrating a semi-analytic model to the
_Renaissance Simulations_ , radiative hydrodynamical simulations of the first
stars and galaxies. An accurate semi-analytic model can make computationally
efficient predictions in a fraction of the time, finishing in less than a day
on a laptop compared to the $\sim$10 million core hours on the Blue Waters
supercomputer used for a single _Renaissance_ simulation. The efficiency of
semi-analytic models allows them to complement other simulations by extending
over larger volumes, reaching lower redshifts, and exploring parameter space
(e.g., Pop III IMF). While this work is focused on the _Renaissance
Simulations_ , the same methodology can be straightforwardly implemented in
other semi-analytic calculations and readily applied to other simulations.

To construct a semi-analytic model consistent with _Renaissance_ , we
assembled a sample of star formation histories and calibrate our model to
them. For metal-enriched star formation, we introduced a two-stage model (that
does not include external enrichment) consisting of an early bursty formation
stage that transitions to a steady stage with star forming at a constant
efficiency. Notably, we implemented stochastic elements during star formation
not typically included in other models by sampling parameters from the
distributions in Table 1. We find that our stochastic two-stage model
typically produced metal-enriched star formation histories consistent with
_Renaissance_ , but completely capturing the variance in individual SFHs was
not the intent of our model.

To calibrate Population III star formation, we implemented a minimum halo mass
consistent with when the first stars form in _Renaissance_. Once halos exceed
this mass, a total Pop III stellar mass was assigned to halos along with a
delay period before the onset of metal-enriched star formation. These values
are sampled from the distributions described in Table 1. Our approach can be
used to easily calibrate to different Pop III prescriptions and incorporate
them into the model.

We determined how closely our model is able to reproduce the _Renaissance_
cosmic star formation history by comparing evolution of SFRDs and found
agreement within a factor of $\sim$2 for redshifts $11<z<25$ for halo masses
$\lesssim$$10^{9}M_{\mathrm{\odot}}$. We also tested our model on both the
rarefied Void and overdense RarePeak regions in _Renaissance_. The metal-
enriched SFRDs produced by our model typically agrees by a factor of $\sim$2
in Void and $\sim$3 in RarePeak, respectively. The difference in SFRDs based
on environment suggests there is physics not currently included in our model
that either leads to higher star formation efficiencies or increased external
enrichment. Therefore, we expect semi-analytic models that do not consider
these effects will fail to reproduce star formation by similar factors.

Additionally, while reproducing the _Renaissance_ star formation was a primary
goal, our Population III prescription could be updated to include more recent
findings. For example, it has been shown that baryon-dark matter streaming,
self-shielding from LW radiation, and X-ray feedback impact the
$\mathrm{H_{2}}$ fraction and therefore halo minimum mass for Pop III
formation. A combination of these effects can raise the minimum halo mass and
delay the formation of the first stars. We found using the minimum dark matter
halo mass for Pop III formation from Kulkarni et al. (2021) that the Pop III
SFRD increases by about an order of magnitude compared to _Renaissance_.
Meanwhile, for $z<18$ the metal-enriched SFRDs agree within a factor of
$\sim$2\. Both LW feedback and metal enrichment could impact star formation
and will be the focus of future work implementing additional calibrations to
_Renaissance_.

Using the grid-based approach of Visbal et al. (2020) to calculate feedback
from LW radiation using fast Fourier transforms, we found our model produces
total ionized fractions that agree within a factor of $\sim$3 to _Renaissance_
between redshifts $14\lesssim z\lesssim 11$ and even closer agreement between
$18\lesssim z\lesssim 14$. When comparing individual H ii regions, we found
that our model typically forms larger ionized bubbles than those in
_Renaissance_. H ii regions in highly clustered environments can suppress the
onset of star formation in neighboring halos and could reduce the overall
amount of metal-enriched and Population III stars.

In future work, we plan to investigate the impact of the environment on star
formation in greater detail. Including dependence on large scale overdensity,
external enrichment, and a more detailed treatment of ionized H ii bubble
formation in our calibration could potentially improve the agreement with
_Renaissance_ for different environments. Our two-stage prescription for
metal-enriched star formation would then be more accurate for a range of
cosmological environments.

We also intend to make observational predictions by applying our calibrated
model to dark matter merger trees that extend over a large statistically
representative volume. Some of these predictions include supernovae rates
(e.g., Scannapieco et al., 2005; Weinmann & Lilly, 2005; Wise & Abel, 2005;
Mesinger et al., 2006; Pan et al., 2012; Hummel et al., 2012; Johnson et al.,
2013; Tanaka et al., 2013; Whalen et al., 2013b; de Souza et al., 2014;
Hartwig et al., 2018), luminosity functions (e.g., Bouwens et al., 2023; Atek
et al., 2023), and abundance predictions in strong gravitational lenses.
Incorporating direct collapse black hole seeds in our model and tracking their
growth could also potentially lead to new insights into the supermassive black
holes found at the centers of most galaxies.

##  6 Acknowledgments

We thank Greg Bryan and Zoltan Haiman for useful discussions. R.H. and E.V.
acknowledge support from NASA ATP grant 80NNSSC22K0629 and NSF grant
AST-2009309. J.W. acknowledges support by NSF grant AST-2108020 and NASA
grants 80NSSC20K0520 and 80NSSC21K1053. The numerical simulations in this
paper were run on the Ohio Supercomputer Center (OSC).

## References

  * Abe et al. (2021) Abe, M., Yajima, H., Khochfar, S., Dalla Vecchia, C., & Omukai, K. 2021, MNRAS, 508, 3226, doi: 10.1093/mnras/stab2637
  * Abel et al. (1997) Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New A, 2, 181, doi: 10.1016/S1384-1076(97)00010-9
  * Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93, doi: 10.1126/science.295.5552.93
  * Ahn et al. (2009) Ahn, K., Shapiro, P. R., Iliev, I. T., Mellema, G., & Pen, U.-L. 2009, ApJ, 695, 1430, doi: 10.1088/0004-637X/695/2/1430
  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, Nature, 622, 707, doi: 10.1038/s41586-023-06521-7
  * Atek et al. (2023) Atek, H., Labbé, I., Furtak, L. J., et al. 2023, arXiv e-prints, arXiv:2308.08540, doi: 10.48550/arXiv.2308.08540
  * Barkana & Loeb (2001) Barkana, R., & Loeb, A. 2001, Phys. Rep., 349, 125, doi: 10.1016/S0370-1573(01)00019-9
  * Behnel et al. (2011) Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science and Engineering, 13, 31, doi: 10.1109/MCSE.2010.118
  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18, doi: 10.1088/0004-637X/763/1/18
  * Borrow et al. (2023) Borrow, J., Kannan, R., Garaldi, E., et al. 2023, MNRAS, 525, 5932, doi: 10.1093/mnras/stad2523
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001, doi: 10.1088/0004-637X/718/2/1001
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145
  * Bowman et al. (2018) Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67, doi: 10.1038/nature25792
  * Brummel-Smith et al. (2019) Brummel-Smith, C., Bryan, G., Butsky, I., et al. 2019, The Journal of Open Source Software, 4, 1636, doi: 10.21105/joss.01636
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Bunker et al. (2023) Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, 677, A88, doi: 10.1051/0004-6361/202346159
  * Chen et al. (2014) Chen, P., Wise, J. H., Norman, M. L., Xu, H., & O’Shea, B. W. 2014, ApJ, 795, 144, doi: 10.1088/0004-637X/795/2/144
  * Cohen et al. (2017) Cohen, A., Fialkov, A., Barkana, R., & Lotem, M. 2017, MNRAS, 472, 1915, doi: 10.1093/mnras/stx2065
  * Crosby et al. (2013) Crosby, B. D., O’Shea, B. W., Smith, B. D., Turk, M. J., & Hahn, O. 2013, ApJ, 773, 108, doi: 10.1088/0004-637X/773/2/108
  * Davé et al. (2012) Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98, doi: 10.1111/j.1365-2966.2011.20148.x
  * de Souza et al. (2014) de Souza, R. S., Ishida, E. E. O., Whalen, D. J., Johnson, J. L., & Ferrara, A. 2014, MNRAS, 442, 1640, doi: 10.1093/mnras/stu984
  * Dekel & Mandelker (2014) Dekel, A., & Mandelker, N. 2014, MNRAS, 444, 2071, doi: 10.1093/mnras/stu1427
  * Dekel & Rees (1987) Dekel, A., & Rees, M. J. 1987, Nature, 326, 455, doi: 10.1038/326455a0
  * Dijkstra et al. (2004) Dijkstra, M., Haiman, Z., Rees, M. J., & Weinberg, D. H. 2004, ApJ, 601, 666, doi: 10.1086/380603
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465, doi: 10.48550/arXiv.2306.02465
  * Fernandez et al. (2014) Fernandez, R., Bryan, G. L., Haiman, Z., & Li, M. 2014, MNRAS, 439, 3798, doi: 10.1093/mnras/stu230
  * Fialkov et al. (2013) Fialkov, A., Barkana, R., Visbal, E., Tseliakhovich, D., & Hirata, C. M. 2013, MNRAS, 432, 2909, doi: 10.1093/mnras/stt650
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Frebel & Norris (2015) Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631, doi: 10.1146/annurev-astro-082214-122423
  * Furlanetto & Mirocha (2022) Furlanetto, S. R., & Mirocha, J. 2022, MNRAS, 511, 3895, doi: 10.1093/mnras/stac310
  * Gao & Theuns (2007) Gao, L., & Theuns, T. 2007, Science, 317, 1527, doi: 10.1126/science.1146676
  * Gnedin (2000) Gnedin, N. Y. 2000, ApJ, 542, 535, doi: 10.1086/317042
  * Gnedin & Hui (1998) Gnedin, N. Y., & Hui, L. 1998, MNRAS, 296, 44, doi: 10.1046/j.1365-8711.1998.01249.x
  * Greif (2015) Greif, T. H. 2015, Computational Astrophysics and Cosmology, 2, 3, doi: 10.1186/s40668-014-0006-2
  * Greif et al. (2011) Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75, doi: 10.1088/0004-637X/737/2/75
  * Haiman et al. (2000) Haiman, Z., Abel, T., & Rees, M. J. 2000, ApJ, 534, 11, doi: 10.1086/308723
  * Haiman et al. (1996a) Haiman, Z., Rees, M. J., & Loeb, A. 1996a, ApJ, 467, 522, doi: 10.1086/177628
  * Haiman et al. (1997) —. 1997, ApJ, 476, 458, doi: 10.1086/303647
  * Haiman et al. (1996b) Haiman, Z., Thoul, A. A., & Loeb, A. 1996b, ApJ, 464, 523, doi: 10.1086/177343
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hartwig et al. (2018) Hartwig, T., Bromm, V., & Loeb, A. 2018, MNRAS, 479, 2202, doi: 10.1093/mnras/sty1576
  * Hegde & Furlanetto (2023) Hegde, S., & Furlanetto, S. R. 2023, MNRAS, 525, 428, doi: 10.1093/mnras/stad2308
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288, doi: 10.1086/375341
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60, doi: 10.1088/0004-637X/781/2/60
  * Hoeft et al. (2006) Hoeft, M., Yepes, G., Gottlöber, S., & Springel, V. 2006, MNRAS, 371, 401, doi: 10.1111/j.1365-2966.2006.10678.x
  * Hopkins et al. (2013) Hopkins, P. F., Cox, T. J., Hernquist, L., et al. 2013, MNRAS, 430, 1901, doi: 10.1093/mnras/stt017
  * Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250, doi: 10.1126/science.1207433
  * Hummel et al. (2012) Hummel, J. A., Pawlik, A. H., Milosavljević, M., & Bromm, V. 2012, ApJ, 755, 72, doi: 10.1088/0004-637X/755/1/72
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ishigaki et al. (2018) Ishigaki, M. N., Tominaga, N., Kobayashi, C., & Nomoto, K. 2018, ApJ, 857, 46, doi: 10.3847/1538-4357/aab3de
  * Johnson et al. (2013) Johnson, J. L., Dalla Vecchia, C., & Khochfar, S. 2013, MNRAS, 428, 1857, doi: 10.1093/mnras/sts011
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in IOS Press, 87–90, doi: 10.3233/978-1-61499-649-1-87
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18
  * Kulkarni et al. (2021) Kulkarni, M., Visbal, E., & Bryan, G. L. 2021, ApJ, 917, 40, doi: 10.3847/1538-4357/ac08a3
  * Kulkarni et al. (2022) Kulkarni, M., Visbal, E., Bryan, G. L., & Li, X. 2022, ApJ, 941, L18, doi: 10.3847/2041-8213/aca47c
  * Liu & Bromm (2020) Liu, B., & Bromm, V. 2020, MNRAS, 497, 2839, doi: 10.1093/mnras/staa2143
  * Loeb & Furlanetto (2013) Loeb, A., & Furlanetto, S. R. 2013, The First Galaxies in the Universe (Princeton University Press)
  * Machacek et al. (2001) Machacek, M. E., Bryan, G. L., & Abel, T. 2001, ApJ, 548, 509, doi: 10.1086/319014
  * Maiolino et al. (2023) Maiolino, R., Uebler, H., Perna, M., et al. 2023, arXiv e-prints, arXiv:2306.00953, doi: 10.48550/arXiv.2306.00953
  * McCaffrey et al. (2023) McCaffrey, J., Hardin, S., Wise, J., & Regan, J. 2023, arXiv e-prints, arXiv:2304.13755, doi: 10.48550/arXiv.2304.13755
  * Mebane et al. (2018) Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2018, MNRAS, 479, 4544, doi: 10.1093/mnras/sty1833
  * Mesinger et al. (2006) Mesinger, A., Johnson, B. D., & Haiman, Z. 2006, ApJ, 637, 80, doi: 10.1086/498294
  * Mocz et al. (2019) Mocz, P., Fialkov, A., Vogelsberger, M., et al. 2019, Phys. Rev. Lett., 123, 141301, doi: 10.1103/PhysRevLett.123.141301
  * Mocz et al. (2020) —. 2020, MNRAS, 494, 2027, doi: 10.1093/mnras/staa738
  * Nebrin et al. (2023) Nebrin, O., Giri, S. K., & Mellema, G. 2023, MNRAS, 524, 2290, doi: 10.1093/mnras/stad1852
  * Noh & McQuinn (2014) Noh, Y., & McQuinn, M. 2014, MNRAS, 444, 503, doi: 10.1093/mnras/stu1412
  * Nomoto et al. (2006) Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424, doi: 10.1016/j.nuclphysa.2006.05.008
  * Okamoto et al. (2008) Okamoto, T., Gao, L., & Theuns, T. 2008, MNRAS, 390, 920, doi: 10.1111/j.1365-2966.2008.13830.x
  * O’Shea & Norman (2007) O’Shea, B. W., & Norman, M. L. 2007, ApJ, 654, 66, doi: 10.1086/509250
  * O’Shea & Norman (2008) —. 2008, ApJ, 673, 14, doi: 10.1086/524006
  * O’Shea et al. (2015) O’Shea, B. W., Wise, J. H., Xu, H., & Norman, M. L. 2015, ApJ, 807, L12, doi: 10.1088/2041-8205/807/1/L12
  * Pan et al. (2012) Pan, T., Kasen, D., & Loeb, A. 2012, MNRAS, 422, 2701, doi: 10.1111/j.1365-2966.2012.20837.x
  * Parsons et al. (2022) Parsons, J., Mas-Ribas, L., Sun, G., et al. 2022, ApJ, 933, 141, doi: 10.3847/1538-4357/ac746b
  * Riaz et al. (2022) Riaz, S., Hartwig, T., & Latif, M. A. 2022, ApJ, 937, L6, doi: 10.3847/2041-8213/ac8ea6
  * Samui et al. (2007) Samui, S., Srianand, R., & Subramanian, K. 2007, MNRAS, 377, 285, doi: 10.1111/j.1365-2966.2007.11603.x
  * Sassano et al. (2021) Sassano, F., Schneider, R., Valiante, R., et al. 2021, MNRAS, 506, 613, doi: 10.1093/mnras/stab1737
  * Scannapieco et al. (2005) Scannapieco, E., Madau, P., Woosley, S., Heger, A., & Ferrara, A. 2005, ApJ, 633, 1031, doi: 10.1086/444450
  * Schaerer (2002) Schaerer, D. 2002, A&A, 382, 28, doi: 10.1051/0004-6361:20011619
  * Schauer et al. (2021) Schauer, A. T. P., Glover, S. C. O., Klessen, R. S., & Clark, P. 2021, MNRAS, 507, 1775, doi: 10.1093/mnras/stab1953
  * Shapiro et al. (1994) Shapiro, P. R., Giroux, M. L., & Babul, A. 1994, ApJ, 427, 25, doi: 10.1086/174120
  * Smith & Lang (2019) Smith, B., & Lang, M. 2019, The Journal of Open Source Software, 4, 1881, doi: 10.21105/joss.01881
  * Smith et al. (2009) Smith, B. D., Turk, M. J., Sigurdsson, S., O’Shea, B. W., & Norman, M. L. 2009, ApJ, 691, 441, doi: 10.1088/0004-637X/691/1/441
  * Sobacchi & Mesinger (2013) Sobacchi, E., & Mesinger, A. 2013, MNRAS, 432, L51, doi: 10.1093/mnrasl/slt035
  * Stecher & Williams (1967) Stecher, T. P., & Williams, D. A. 1967, ApJ, 149, L29, doi: 10.1086/180047
  * Susa et al. (2014) Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32, doi: 10.1088/0004-637X/792/1/32
  * Tanaka et al. (2013) Tanaka, M., Moriya, T. J., & Yoshida, N. 2013, MNRAS, 435, 2483, doi: 10.1093/mnras/stt1469
  * Tegmark et al. (1997) Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1, doi: 10.1086/303434
  * Thoul & Weinberg (1996) Thoul, A. A., & Weinberg, D. H. 1996, ApJ, 465, 608, doi: 10.1086/177446
  * Trenti & Stiavelli (2009) Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879, doi: 10.1088/0004-637X/694/2/879
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Visbal et al. (2012) Visbal, E., Barkana, R., Fialkov, A., Tseliakhovich, D., & Hirata, C. M. 2012, Nature, 487, 70, doi: 10.1038/nature11177
  * Visbal et al. (2020) Visbal, E., Bryan, G. L., & Haiman, Z. 2020, ApJ, 897, 95, doi: 10.3847/1538-4357/ab994e
  * Visbal et al. (2015) Visbal, E., Haiman, Z., & Bryan, G. L. 2015, MNRAS, 450, 2506, doi: 10.1093/mnras/stv785
  * Visbal et al. (2018) —. 2018, MNRAS, 475, 5246, doi: 10.1093/mnras/sty142
  * Visbal et al. (2014) Visbal, E., Haiman, Z., Terrazas, B., Bryan, G. L., & Barkana, R. 2014, MNRAS, 445, 107, doi: 10.1093/mnras/stu1710
  * Walch & Naab (2015) Walch, S., & Naab, T. 2015, MNRAS, 451, 2757, doi: 10.1093/mnras/stv1155
  * Wang et al. (2022) Wang, X., Cheng, C., Ge, J., et al. 2022, arXiv e-prints, arXiv:2212.04476, doi: 10.48550/arXiv.2212.04476
  * Weinmann & Lilly (2005) Weinmann, S. M., & Lilly, S. J. 2005, ApJ, 624, 526, doi: 10.1086/428106
  * Whalen et al. (2013a) Whalen, D. J., Fryer, C. L., Holz, D. E., et al. 2013a, ApJ, 762, L6, doi: 10.1088/2041-8205/762/1/L6
  * Whalen et al. (2013b) Whalen, D. J., Even, W., Frey, L. H., et al. 2013b, ApJ, 777, 110, doi: 10.1088/0004-637X/777/2/110
  * Wise & Abel (2005) Wise, J. H., & Abel, T. 2005, ApJ, 629, 615, doi: 10.1086/430434
  * Wise & Abel (2007) —. 2007, ApJ, 671, 1559, doi: 10.1086/522876
  * Wise & Abel (2011) —. 2011, MNRAS, 414, 3458, doi: 10.1111/j.1365-2966.2011.18646.x
  * Wise et al. (2012a) Wise, J. H., Abel, T., Turk, M. J., Norman, M. L., & Smith, B. D. 2012a, MNRAS, 427, 311, doi: 10.1111/j.1365-2966.2012.21809.x
  * Wise et al. (2012b) Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012b, ApJ, 745, 50, doi: 10.1088/0004-637X/745/1/50
  * Xu et al. (2014) Xu, H., Ahn, K., Wise, J. H., Norman, M. L., & O’Shea, B. W. 2014, ApJ, 791, 110, doi: 10.1088/0004-637X/791/2/110
  * Xu et al. (2016) Xu, H., Wise, J. H., Norman, M. L., Ahn, K., & O’Shea, B. W. 2016, ApJ, 833, 84, doi: 10.3847/1538-4357/833/1/84
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645, doi: 10.1086/375810
  * Yung et al. (2019) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241