# A self-consistent semi-analytic model for Population III star formation in
minihalos

Sahil Hegde1 and Steven R. Furlanetto1
1 Department of Physics & Astronomy, University of California, Los Angeles,
475 Portola Plaza, Los Angeles, CA 90095, USA
E-mail<EMAIL_ADDRESS>(SH)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The formation of the first stars marks a watershed moment in the history of
our universe. As the first luminous structures, these stars (also known as
Population III, or Pop III stars) seed the first galaxies and begin the
process of reionization. We construct an analytic model to self-consistently
trace the formation of Pop III stars inside minihalos in the presence of the
fluctuating ultraviolet background, relic dark matter-baryon relative
velocities from the early universe, and an X-ray background, which largely
work to suppress cooling of gas and delay the formation of this first
generation of stars. We demonstrate the utility of this framework in a semi-
analytic model for early star formation that also follows the transition
between Pop III and Pop II star formation inside these halos. Using our new
prescription for the criteria allowing Pop III star formation, we follow a
population of dark matter halos from $z=50$ through $z=6$ and examine the
global star formation history, finding that each process defines its own key
epoch: _(i)_ the stream velocity dominates at the highest redshifts ($z\gtrsim
30$), _(ii)_ the UV background sets the tone at intermediate times ($30\gtrsim
z\gtrsim 15$), and _(iii)_ X-rays control the end of Pop III star formation at
the latest times ($z\lesssim 15$). In all of our models, Pop III stars
continue to form down to $z\sim 7-10$, when their supernovae will be
potentially observable with forthcoming instruments. Finally, we identify the
signatures of variations in the Pop III physics in the global 21-cm spin-flip
signal of atomic hydrogen.

###### keywords:

population III stars – cosmology – high-redshift galaxies

††pubyear: 2023††pagerange: A self-consistent semi-analytic model for
Population III star formation in minihalos–References

##  1 Introduction

Formed in the pristine hydrogen and helium gas of the early universe, the
first stars, also known as Population III (Pop III) stars, are thought to have
been very different from the stars we observe today (Pop I/II) (Bromm, 2013,
Loeb & Furlanetto, 2013). In these metal-free clouds, Pop III stars must rely
on the radiative transitions of molecular hydrogen ($\mathrm{H}_{2}$) to cool,
fragment, and collapse to high densities. However, this cooling is
comparatively inefficient (relative to atomic hydrogen cooling at high
temperatures or cooling in metal-enriched gas via transitions of, e.g., CO,
H2O, etc.) and is believed to result in a population of supermassive stars.
These stars are expected to form in so-called dark matter (DM) ‘minihalos’ of
$10^{5}-10^{6}\ M_{\odot}$ and likely formed in isolation or groups of a few
(Abel et al., 2002, Bromm, 2013). At such large masses, these stars will have
short ($\sim 5$ Myr) lives, some of which may have ended their lives in
superluminous pair-instability supernovae (SNe; Barkat et al., 1967, Fryer et
al., 2001, Heger & Woosley, 2002, Heger et al., 2003). In so doing, they
ejected metals into the interstellar and intergalactic media, paving the way
for future generations of metal-enriched star formation (Ferrara et al., 2000,
Madau et al., 2001, Furlanetto & Loeb, 2003).

Observationally probing the epoch of the first stars, which is colloquially
referred to as ‘Cosmic Dawn’, is a difficult endeavor. Given their short
lifetimes and their formation in small groups, Pop III stars and their host
minihalos were likely very faint (Mebane et al., 2018). In addition, these
halos likely only dominated in the very early universe, requiring observations
that can probe to $z\gtrsim 20-30$. To that end, observations of the indirect
signatures of these stars, such as in their SNe or through the highly
redshifted 21-cm line, offer a far more promising avenue to study Cosmic Dawn
(e.g., Mebane et al., 2020, Magg et al., 2022, Muñoz et al., 2022, Lazar &
Bromm, 2022). However, even the most promising of these avenues, observations
of the 21-cm line, will be challenging. Therefore, leveraging robust theory
models will be crucial to the design of instruments that are sensitive to the
relevant areas of parameter space and development of software pipelines to
robustly infer those parameters.

Pop III halos have been studied theoretically with a variety of methods (in
order of decreasing computational cost): numerical simulations (e.g., Machacek
et al., 2001, Abel et al., 2002, Wise & Abel, 2007, O’Shea & Norman, 2008,
Maio et al., 2010, Stacy et al., 2012, Hirano et al., 2015, Xu et al., 2016,
Sarmento et al., 2018, Park et al., 2021, Kulkarni et al., 2021, Schauer et
al., 2021), semi-analytic models (e.g., Trenti et al., 2009, Jaacks et al.,
2018, Visbal et al., 2018, Mebane et al., 2018, Visbal et al., 2020, Magg et
al., 2022), and analytic calculations (e.g., Tegmark et al., 1997, Haiman et
al., 2000, McKee & Tan, 2008, Kulkarni et al., 2013, Ricotti, 2016). Because
the models are essentially unconstrained in the absence of observations,
exploring the breadth of the available parameter space is challenging,
especially with detailed numerical calculations. Semi-analytic models offer a
compromise between the computational efficiency of an analytic calculation and
the rigor of a numerical simulation. In practice, applications of these models
to high redshift star formation, such as the model described in this work,
begin with a numerical model for dark matter (DM) halo growth (such as through
abundance matching to a simulated halo mass function or through simulated halo
merger trees) and fold in the physics of star formation (initial mass
function, supernovae, feedback, etc.) through analytic approximations. In such
models, these two pieces are connected with the minimum DM halo mass for Pop
III stars to form.

The minimum mass is the host halo mass scale above which a star forming cloud
is able to form, cool, and collapse, driven by the presence of
$\mathrm{H}_{2}$ (Haiman et al., 1996a, Haiman et al., 1996b, Tegmark et al.,
1997, Machacek et al., 2001). As a result, an understanding of processes that
affect the amount of $\mathrm{H}_{2}$ present in a halo is crucial for
understanding when and where early star formation can occur. There are several
key processes believed to significantly affect the $\mathrm{H}_{2}$ fraction
in a halo—the relic velocity between baryons and DM in the early universe
(environmental) and the buildup of X-ray and photodissociating UV backgrounds
(feedback) (e.g., Loeb & Furlanetto, 2013, Ricotti, 2016, Mebane et al., 2018,
Schauer et al., 2021, Kulkarni et al., 2021, Muñoz et al., 2022, Nebrin et
al., 2023).

The effect of a strong UV background, specifically in the Lyman-Werner (LW)
bands of $\mathrm{H}_{2}$ (11.2-13.6 eV), is relatively well-studied, both
analytically and in numerical simulations (e.g., Tegmark et al., 1997,
Machacek et al., 2001, Shang et al., 2010, Wolcott-Green et al., 2011, Visbal
et al., 2014b, Kulkarni et al., 2021). As stars form, they radiate and build
up a metagalactic LW background, which breaks down $\mathrm{H}_{2}$ and
suppresses later generations of star formation. The relative DM-baryon motion,
often referred to as the stream velocity, has only come to be thought of as
significant in recent years and, as such, is relatively less studied (e.g.,
Tseliakhovich & Hirata, 2010, Dalal et al., 2010, Naoz et al., 2013, Fialkov,
2014). The stream velocity suppresses accretion and limits the gas fraction in
the least massive halos. As such, these two processes both negatively impact
the $\mathrm{H}_{2}$ content of a halo. X-rays, on the other hand, are more
complicated, as they can both positively and negatively feed back into the
star formation process through photoionization and heating of the IGM,
respectively (e.g., Machacek et al., 2003, Ricotti & Ostriker, 2004, Hummel et
al., 2015, Ricotti, 2016, Park et al., 2021). Each of these mechanisms can be
parameterized in terms of the minimum DM halo mass for star formation, and
most work has focused on studying these processes independently. The recent
simulations of Kulkarni et al. (2021) and Schauer et al. (2021) are some of
the first attempts to study the joint effects of the two negative processes,
the LW photodissociation and stream velocity, but present discrepant estimates
for this minimum mass scale.

Motivated by these uncertainties and the lack of a single model that
encompasses all of these effects,111We note, however, that while this paper
was in the final stages of preparation, Nebrin et al. (2023) presented a
calculation for the minimum halo mass to host Pop III star formation that
included the LW background and streaming mechanisms, in addition to several
others with more modest effects, but not the effects of X-rays. We will
compare to their model throughout. in this work we present an analytic
calculation for the minimum star-forming halo mass that incorporates the
aforementioned three processes. We update the semi-analytic model presented in
Mebane et al. (2018) to explore the effects of variations in the underlying
Pop III physics on the global star formation rate density. Finally, we use
these results to make preliminary observational predictions, specifically in
calculating Pop III SN rates and the global 21-cm signal.

The paper is structured as follows. In § 2 we outline our analytic calculation
of the minimum mass and present a fitting formula that summarizes these
calculations. In § 3, we describe the semi-analytic model and our improvements
to Mebane et al. (2018). In § 4 we present our results—namely the fiducial
model and the effects of varying the physics incorporated into the minimum
mass. In § 5 we compare our results to previous works and in § 6 we present
potential observable signatures of Pop III star formation. Finally, in § 7, we
summarize our main conclusions.

In this work we use a flat $\Lambda$CDM cosmology with $\Omega_{\rm
m}=0.3111$, $\Omega_{\Lambda}=0.6889$, $\Omega_{\rm b}=0.0489$,
$\sigma_{8}=0.8102$, $n_{s}=0.9665$, and $h=0.6766$, consistent with the
results of Planck Collaboration et al. (2020).

##  2 The Minimum Mass for Pop III Star Formation

In this section, we will describe the components of our simple analytic model
of the aforementioned processes—the metagalactic LW background, stream
velocity, and X-ray background—to generate an estimate of the minimum star-
forming halo mass scale. Schematically, this is broken down as follows: a halo
must first accrete gas (filter mass; § 2.1) and then must be able to
efficiently cool (cooling threshold; § 2.2). These two thresholds are modified
in the presence of a DM-baryon relative velocity (Sections 2.1.2 and 2.4).
Once the first stars form, they produce a UV background (dissociation
threshold; § 2.3) and an X-ray background (§ 2.5) that can affect the minimum
mass for subsequent generations of stars.

###  2.1 Accretion mass

In the absence of both a LW background and a stream velocity, there are two
relevant processes — a halo must be able to accrete baryons and must then be
able to cool and form stars. The former of these is governed by the so-called
‘filtering’ scale, which is the scale below which baryonic perturbations are
suppressed relative to the DM fluctuations (Gnedin & Hui, 1998) by thermal
pressure. This scale, or the associated mass scale (the filter mass), is
sometimes referred to as the cosmological Jeans mass because it resembles a
time-averaged Jeans mass and effectively sets a minimum halo mass necessary
for gas accretion.

####  2.1.1 No stream velocity

Building from the calculation of the filtering scale outlined in Naoz &
Barkana (2007), Naoz et al. (2013) calculate the filter mass using numerical
simulations and show that the results can be reproduced with an analytic
calculation from linear theory. Inspired by the calculations of Naoz et al.
(2013), in this section we describe a simple model that produces remarkable
agreement with the analytic model therein.

Following the discussion in Gnedin (2000), we calculate the filtering
wavenumber, $k_{F}$, defined as the scale at which the baryonic perturbations
grow significantly compared to the DM fluctuations. Assuming that the
perturbations satisfy $\delta_{b}(t,k=0)=\delta_{c}(t,k=0)$, the ratio of
these quantities can be expanded as

$\frac{\delta_{b}}{\delta_{c}}\approx 1-\frac{k^{2}}{k_{F}^{2}}$ (1)

where we have dropped the higher order terms and assumed that the
perturbations are small. If we introduce the growth factor $D(t)$ and define
$A(t)\equiv D(t)/k_{F}^{2}$, this becomes

$\frac{\delta_{b}}{\delta_{c}}=1-\frac{A(t)}{D(t)}k^{2}$ (2)

The linearized Euler equation can be recast in terms of the density
perturbations as

$\frac{\partial^{2}\delta}{\partial t^{2}}+2H\frac{\partial\delta}{\partial
t}=4\pi G\bar{\rho}\delta-\frac{c_{s}^{2}k^{2}}{a^{2}}\delta$ (3)

Substituting eq. 3 for baryon and DM overdensities into eq. 2 yields

$\frac{d^{2}A}{dt^{2}}+2H\frac{dA}{dt}=\frac{c_{s}^{2}}{a^{2}}D(t)$ (4)

which, when solved, gives the evolution of $A(t)$, or, equivalently, the
filtering scale $k_{F}(a)$. Noting that the Jeans wavenumber can be written as
$k_{J}^{2}=(3/2)\Omega_{m}H^{2}a^{2}/c_{s}^{2}$, the full solution of eq. 4
can be cast in terms of $k_{F}^{2}$:

$\frac{1}{k_{F}^{2}}=\frac{1}{D(t)}\int_{0}^{t}dt^{\prime}a^{2}(t^{\prime})\frac{\ddot{D}(t^{\prime})+2H(t^{\prime})\dot{D}(t^{\prime})}{k_{J}^{2}(t^{\prime})}\int_{t^{\prime}}^{t^{\prime\prime}}\frac{dt^{\prime\prime}}{a^{2}(t^{\prime\prime})}$
(5)

For large redshifts ($z\gtrsim 5$; the relevant regime here), we know that
$\Omega_{m}\to 1$, so we can use the results for an Einstein-de Sitter
cosmology; i.e., that $a\propto t^{2/3}$ and $D\propto a$. With these, eq. 5
simplifies to

$\frac{1}{k_{F}^{2}(t)}=\frac{3}{a}\int_{0}^{a}\frac{da^{\prime}}{k_{J}^{2}(a^{\prime})}\bigg{(}1-\sqrt{\frac{a^{\prime}}{a}}\bigg{)}$
(6)

####  2.1.2 Including the DM-baryon relative velocity

To this point, we have ignored the effects of a relative velocity between DM
and baryons. However, Tseliakhovich & Hirata (2010) showed that the relative
velocity between DM and baryon density fluctuations following recombination
can have important consequences for structure formation in the early universe.
This relative velocity, often referred to as the stream velocity, is a product
of the differing growth histories of DM and baryon perturbations prior to
recombination and decoupling. That is, because they only interact
gravitationally, DM particles began collapsing into their potential wells well
before the end of recombination. Baryons, on the other hand, were prevented
from such collapse by the thermal pressure that resulted from their coupling
to the photon field. Once recombination had concluded, the photons and baryons
decoupled and the baryons were free to fall into DM potential wells. However,
they retained memory of their motion prior to that point and a spatially-
varying distribution of relative velocities was produced as a result. It has
since been shown, analytically and with numerical simulations, that this
relative stream velocity has a variety of effects on structure formation in
the early universe, ranging from reducing the number density of DM halos to
suppressing the gas content of halos (e.g., Dalal et al., 2010, Tseliakhovich
& Hirata, 2010, O’Leary & McQuinn, 2012, McQuinn & O’Leary, 2012, Naoz et al.,
2012, Naoz et al., 2013, Fialkov, 2014, Williams et al., 2022, Lake et al.,
2023). It is the latter phenomenon that is especially relevant for Pop III
star formation, as a suppression of gas accretion can delay $\mathrm{H}_{2}$
cooling and collapse.222The suppression of halo number density has important
implications for structure formation, especially at the highest redshifts,
when we expect the very first stars to form. However, we defer analysis of
this contribution to future work and here limit our focus to the effect of the
stream velocity on the gas being used to form stars.

Naoz et al. (2013) carry out the full self-consistent calculation of this
effect, which requires modification of the definition of the filtering scale
to account for large scale structure and inclusion of the temperature and
density fluctuations in the coupled differential equations describing the
growth of perturbations (i.e., eq. 3). As a simple approximation to this
calculation, we propose taking the stream velocity as a modification to the
sound speed (which appears only in calculation of the Jeans wavenumber) and
adding it in quadrature:

$c_{s}\mapsto\sqrt{Ac_{s}^{2}+v_{\mathrm{bc}}^{2}}$ (7)

where we include the constant $A$ to account for differences between our
calculation of the sound speed and that of Naoz et al. (2013). Calibrating our
calculation to the results of Naoz et al. (2013), we set $A=0.64$.

In Figure 1, we present the results of this calculation for various values of
the stream velocity.333We conventionally refer to the stream velocity in
multiples of the root-mean-square value of its Maxwellian distribution, which
is $\sigma_{\rm vbc}=30\ {\rm km\ s^{-1}}$ at $z_{\rm rec}=1100$. While our
results do not precisely agree with the analytic calculation of Naoz et al.
(2013), the approximation we employ performs remarkably well, as our estimates
of the filter mass are discrepant by less than a factor of 2 and the agreement
improves with increasing $v_{\rm bc}$. In Figure 2, we compare our analytic
approximation to the numerical simulations of Naoz et al. (2013),
demonstrating similar agreement to that seen in Figure 1. Note that to compare
the results of this three-dimensional analytic calculation to the simulated
results, we scale the velocity by a factor of $1/\sqrt{3}$ because the stream
velocity acts in a particular direction, whereas the sound speed is isotropic.

We find that the following fitting formula reproduces our treatment of the
filter mass quite well:

$M_{F}(v_{\rm bc},z)\simeq 1.66\times 10^{4}\bigg{(}1+\frac{v_{\rm
bc}}{\sigma_{\rm vbc}}\bigg{)}^{5.02}\bigg{(}\frac{1+z}{21}\bigg{)}^{0.85}$
(8)

Figure 1: The filter mass calculated using our approximation (solid) compared
to the results of Naoz et al. (2013) (dashed) and the fit to our calculation
(eq. 8; dotted) for various values of the stream velocity (different colors).
Figure 2: The filter mass calculated using our approximation (solid curves)
compared to the numerically simulated results of Naoz et al. (2013) (points)
for various values of the stream velocity (different colors).

###  2.2 Molecular hydrogen cooling threshold

Once a halo is able to accrete gas, we must then find the threshold mass
necessary to build up sufficient molecular hydrogen so that the gas can cool
and collapse to high densities. Molecular hydrogen is primarily formed through
a two step process, as follows:

$\displaystyle{\rm H}+{\rm e}^{-}$ $\displaystyle\to{\rm H}^{-}+h\nu\ (k_{9})$
(9) $\displaystyle{\rm H}^{-}+{\rm H}$ $\displaystyle\to{\rm H}_{2}+{\rm
e}^{-}\ (k_{10})$ (10)

where the associated reaction rate coefficients $k_{i}$ are given in Shang et
al. (2010). Note that for this calculation we use the T4 spectrum rate
coefficients given therein. Accounting for the formation of $\mathrm{H}_{2}$
via these two processes and the destruction of ${\rm H}^{-}$ by CMB photons
(with associated rate coefficient $k_{25}$), the net production rate
coefficient of $\mathrm{H}_{2}$ is

$k_{\mathrm{form}}=k_{9}\bigg{[}\frac{k_{10}n_{\rm H}}{k_{10}n_{\rm
H}+k_{25}}\bigg{]},$ (11)

where $n_{H}$ is the hydrogen number density and we have limited our focus to
$\mathrm{H}_{2}$ formation in the gas-phase.444Once a halo has been enriched
by early generations of star formation, $\mathrm{H}_{2}$ formation could, in
principle, be catalyzed by the presence of dust grains (Nakatani et al.,
2020). However, for $T_{\rm vir}\leq 10^{4}\ {\rm K}$, where $\mathrm{H}_{2}$
cooling will be dominant, Nebrin et al. (2023) find that the dust-catalyzed
formation rate will be smaller than the gas-phase formation rate except in the
cases of a highly-enriched halo forming $\mathrm{H}_{2}$ with maximal
efficiency. Yamaguchi et al. (2023) find that the volume-filling fraction of
metal-enriched wind bubbles will be $\lesssim 10$% and the average metallicity
of the universe $\langle Z\rangle/Z_{\odot}\lesssim 10^{-3}$ for $z\gtrsim 6$
in both their optimistic Pop III and Pop II models. Therefore, even if a
significant number of halos are forming Pop II stars, it is unlikely that
these halos will be able to enrich nearby Pop III halos to a degree that dust
catalysis will be the dominant $\mathrm{H}_{2}$ formation channel.

The halo density is therefore a key input to the cooling process. Motivated by
their numerical simulations, Visbal et al. (2014a) argued that high-$z$
minihaloes fall into two regimes. At very small masses, the maximum gas
density in a halo follows $n_{H}\propto T_{\rm vir}^{3/2}$, which reflects the
maximum temperature allowed by adiabatic compression from the ambient IGM
temperature to the halo virial temperature. For these small haloes, the IGM
entropy is sufficiently large, compared to the entropy generated during halo
collapse, that this adiabatic limit is a good approximation. However, at
larger halo masses the entropy is dominated by halo formation, and in this
regime the gas settles into a “universal” profile in which the central core
density is independent of halo mass (at a fixed redshift). To normalize the
scalings in these two regimes, we fit to the gas densities found by O’Leary &
McQuinn (2012), which yields555We have chosen to normalize to these results,
which fall a factor of a few below that of Visbal et al. (2014a), in order to
estimate the typical halo gas density rather than the maximum. Variations in
this normalization are effectively folded into the choice of our free
parameter, $\zeta$ (introduced later). We discuss the effect of making a
different density choice (e.g., Nebrin et al., 2023) in more detail in § 5.1.

$n_{\rm H}\simeq\begin{dcases}6.19\ \bigg{(}\frac{T_{\rm vir}}{10^{3}\ {\rm
K}}\bigg{)}\ \mathrm{cm}^{-3},&{\rm for}\ T_{\rm vir}<2\times 10^{3}\ {\rm
K}\\\ 12.38\ \bigg{(}\frac{1+z}{21}\bigg{)}^{3}\ \mathrm{cm}^{-3},&{\rm for}\
T_{\rm vir}\geq 2\times 10^{3}\ {\rm K}.\end{dcases}$ (12)

We note that the transition between these regimes occurs at a redshift-
dependent halo mass,

$M_{\rm turn}=9.64\times 10^{5}\bigg{(}\frac{1+z}{21}\bigg{)}^{-3/2}\
M_{\odot}.$ (13)

Defining the ionized fraction $x_{\rm HII}\equiv n_{\rm HII}/n_{\rm H}$ and
molecular fraction $f_{\mathrm{H}_{2}}\equiv n_{\mathrm{H}_{2}}/n_{\rm H}$,
the above pair of reactions (eq. 9) yield the evolution equations

$\displaystyle\dot{x}_{\rm HII}$ $\displaystyle=-\alpha_{B}n_{\rm H}x_{\rm
HII}^{2}$ (14) $\displaystyle\dot{f}_{\mathrm{H}_{2}}$ $\displaystyle=k_{\rm
form}n_{\rm H}x_{\rm HII}(1-x_{\rm HII}-2f_{\mathrm{H}_{2}})$ (15)

where we have introduced the Case B recombination coefficient $\alpha_{B}$ to
account for electron depletion due to recombination with H atoms.

Following the discussion outlined in Tegmark et al. (1997), in the limit of
inefficient cooling (where $n$ and $T$ are constant), these can be solved to
find that

$f_{\mathrm{H}_{2}}(t)-f_{\mathrm{H}_{2}}^{i}\approx\frac{k_{\rm
form}}{\alpha_{B}}\ln(1+t/t_{\rm rec}^{i})$ (16)

where $t_{\rm rec}^{i}\equiv(\alpha_{B}n_{\rm H}x_{\rm HII})^{-1}$ is the
initial recombination timescale. For $t\ll t_{\rm rec}^{i}$, the density of
electrons is sufficiently large that $\mathrm{H}_{2}$ is produced at a
constant rate $k_{\rm form}/\alpha_{B}$. For $t\gg t_{\rm rec}^{i}$, electrons
will have been sufficiently depleted and the $\mathrm{H}_{2}$ fraction will
grow slowly. Therefore, when $t\sim t_{\rm rec}^{i}$, the $\mathrm{H}_{2}$
fraction will reach a saturation level

$f_{\mathrm{H}_{2},\rm sat}\approx\frac{k_{\rm form}}{\alpha_{B}}\approx
4.97\times 10^{-4}\bigg{(}\frac{T}{10^{3}\ {\rm K}}\bigg{)}^{1.52}$ (17)

which is a factor of $\sim 1.5$ larger than the value reported in Tegmark et
al. (1997), where the difference results from the incorrect use of the Case A
recombination coefficient in that calculation (as noted by Nebrin et al.,
2023).666To calculate the second equality, we have approximated $k_{\rm
form}\approx k_{9}$ (because $k_{25}\ll nk_{10}$ for the gas densities of
relevance here) and made a power-law fit to the Case B recombination
coefficient given in Draine (2011) between $30\ {\rm K}\leq T\leq 10^{4}\ {\rm
K}$.

With this in hand, we need to find the critical level of $\mathrm{H}_{2}$ that
needs to build up for cooling to become efficient. This can be found by
comparing the cooling time of a cloud of gas to (a fraction $\zeta$ of) the
then-current Hubble time; i.e., $t_{\rm cool}<\zeta t_{H}$. The cooling time
is given by

$t_{\rm cool}=\frac{3k_{B}T_{\rm vir}}{2\Lambda(n_{\rm H},T_{\rm
vir})f_{\mathrm{H}_{2}}}$ (18)

where for the molecular hydrogen cooling function we use the approximation to
that found by Galli & Palla (1998), valid between the temperatures $120\ {\rm
K}$ and $6400\ {\rm K}$ (Trenti & Stiavelli, 2009):

$\Lambda(n_{\rm H},T)\simeq 10^{-31.6}\bigg{(}\frac{T}{100\ {\rm
K}}\bigg{)}^{3.4}\bigg{(}\frac{n_{\rm H}}{10^{-4}{\rm cm^{-3}}}\bigg{)}\ {\rm
erg\ s^{-1}}$ (19)

For the era of the first stars, when $z\gg 1$, we can again approximately use
the results for an Einstein-de Sitter cosmology, for which the Hubble time is

$t_{H}=\frac{2}{3H_{0}}\Big{[}\Omega_{m}(1+z)^{3}\Big{]}^{-1/2}\approx
6.52\times 10^{9}\ \Big{[}\Omega_{m}h^{2}(1+z)^{3}\Big{]}^{-1/2}\ \mathrm{yr}$
(20)

Enforcing the condition that cooling must occur on timescales shorter than a
fraction $\zeta$ of the Hubble time (meant to represent the rate at which a
halo accumulates thermal energy through accretion), we find that the critical
$\mathrm{H}_{2}$ density needed for efficient cooling at a virial temperature
$T$ is

$n_{\mathrm{H}_{2},\rm crit}=1.534\times
10^{-4}\zeta^{-1}\big{(}\Omega_{m}h^{2}\big{)}^{1/2}\bigg{(}\frac{T}{10^{3}\
{\rm K}}\bigg{)}^{-2.4}\bigg{(}\frac{1+z}{21}\bigg{)}^{3/2}\ {\rm cm^{-3}}$
(21)

From this, we can compute the critical $\mathrm{H}_{2}$ fraction as
$f_{\mathrm{H}_{2},\rm crit}=n_{\mathrm{H}_{2},\rm crit}/n_{\rm H}$. With the
virial temperature given in Barkana & Loeb (2001) and plugging in the values
for our chosen cosmology, this can be written as

$f_{\mathrm{H}_{2},\rm crit}\simeq 8.17\times
10^{-7}\zeta^{-1}\bigg{(}\frac{M}{10^{6}M_{\odot}}\bigg{)}^{-2.27}\bigg{(}\frac{1+z}{21}\bigg{)}^{-1.9}$
(22)

It turns out that the critical $\mathrm{H}_{2}$ fraction will be achieved in
the low-mass regime (in eq. 12) for the redshifts of interest, so we have
omitted the solution associated with the high-mass halo density. We will see
that once we introduce a photodissociating Lyman-Werner background, as is
discussed in § 2.3, the high-mass regime will become important.

Comparing this threshold to the molecular hydrogen fraction in a halo (eq.
17), we can solve for the critical virial temperature (or mass) necessary for
efficient cooling (setting $\zeta=0.25$ for our fiducial calculations; see
Visbal et al., 2014b). This yields the cooling threshold

$M_{\rm cool}(z)\simeq 1.55\times
10^{5}\bigg{(}\frac{\zeta}{0.25}\bigg{)}^{-0.3}\bigg{(}\frac{1+z}{21}\bigg{)}^{-1}\
M_{\odot}$ (23)

This threshold (along with the filter mass; see § 2.1) sets a baseline value
for the minimum halo mass for efficient cooling in the absence of any
radiation backgrounds or other external effects and is displayed in the flat
part of the curves in Figure 3. That is, for each $z$, there is a critical
value of the Lyman-Werner background intensity below which the cooling
threshold sets the minimum mass.

###  2.3 The photodissociating Lyman-Werner background

Once the first stars form, they produce radiation backgrounds that make
subsequent generations of star formation more complex by affecting the amount
of $\mathrm{H}_{2}$ available to cool. The primary process that suppresses the
$\mathrm{H}_{2}$ content of a halo is negative feedback associated with the
Lyman-Werner background; i.e. radiation of UV photons with energies of
11.2-13.6 eV that can dissociate molecular hydrogen. Therefore, as more stars
form, a growing LW background builds up (e.g., Visbal et al., 2014b) and
future generations of star formation are delayed and suppressed. In practice,
the buildup of the LW background drives an increase in the minimum mass for
star formation.

In order to get a minimum mass scale for a halo to be able to efficiently cool
in the presence of a photo-dissociating LW background, we compare two
$\mathrm{H}_{2}$ number density thresholds: the collapse threshold, which is
set by the cooling time, $n_{H_{2}}^{\mathrm{crit}}$ (eq. 21), and the
dissociation equilibrium threshold, $n_{H_{2}}^{\mathrm{eq}}$ (i.e., the
equilibrium $\mathrm{H}_{2}$ density that results from a fixed LW background
intensity). The former of these can be parameterized by the virial temperature
of the halo, $T_{\mathrm{vir}}$, and the latter by the intensity in the LW
bands, which is conventionally expressed in units of
$J_{21}=10^{-21}\mathrm{erg}\ \mathrm{s}^{-1}\ \mathrm{cm}^{-2}\
\mathrm{Hz}^{-1}\ \mathrm{sr}^{-1}$. Equating the two thresholds yields a
maximum value of the LW intensity that allows a cloud with a fixed virial
temperature to collapse or, equivalently, yields the critical virial
temperature needed for collapse under a fixed value of the LW background
intensity.

The dissociation equilibrium threshold is set by balancing the rate of
$\mathrm{H}_{2}$ formation and destruction, comparing the rates of
photodissociation (pd) and collisional dissociation (cd):

$f_{H_{2}}^{\mathrm{eq}}=\min\bigg{(}\frac{k_{\mathrm{form}}}{k_{\mathrm{pd}}}f_{e},\frac{k_{\mathrm{form}}}{k_{\mathrm{cd}}}f_{e}\bigg{)}$
(24)

where we use the recombination code CosmoREC (Chluba & Thomas, 2013) to
compute the electron fraction $f_{e}$ and the reaction rate coefficients are
the same as defined above (Shang et al., 2010).

This equilibrium is modified by the effects of self-shielding—in halos where a
sufficiently high column density of $\mathrm{H}_{2}$ builds up, the gas can
become optically thick to LW radiation and ‘shield’ itself against
photodissociation. This is usually parameterized in the form of a ‘shielding
factor’, $f_{\rm shield}$, which reduces the photodissociation rate:

$k_{\mathrm{pd}}(N_{H_{2}},J_{\rm
LW})=f_{\mathrm{shield}}(N_{H_{2}},T)k_{\mathrm{pd}}(N_{H_{2}}=0,J_{\rm LW})$
(25)

so a value of $f_{\rm shield}=1$ (0) would correspond to no (complete)
shielding. We approximately calculate the column density using the peak halo
gas density given in eq. 12—i.e.,
$N_{\mathrm{H}_{2}}=f_{\mathrm{H}_{2}}\lambda_{\mathrm{J}}n_{\mathrm{\rm
H}}/2$, where $\lambda_{J}$ is the Jeans length in the halo gas cloud. We
compute $f_{\mathrm{shield}}$ following Wolcott-Green et al. (2011) (WG11) and
Wolcott-Green & Haiman (2019) (WG19). Noting that the Wolcott-Green & Haiman
(2019) expression is a correction to Wolcott-Green et al. (2011) at higher
densities, our shield factor is given by

$f_{\mathrm{shield}}=\left\\{\begin{array}[]{ll}f(\alpha=1.1)\
(\mathrm{WG11})&n\leq 10^{3}\ \mathrm{cm}^{-3}\\\ f(\alpha=\mathrm{eq\ 8})\
(\mathrm{WG19})&\begin{aligned} 10^{3}\leq n/\mathrm{cm}^{-3}\leq 10^{7},\\\
N_{H_{2}}\leq 10^{17}\ \mathrm{cm}^{-2},\\\ T\leq 8000\
\mathrm{K}\end{aligned}\\\ 1&\mathrm{otherwise}\\\ \end{array}\right.$ (26)

where the limits reflect those set by WG11/19. Equating eq. 21 and
$n_{\mathrm{eq}}=n_{\rm H}f_{\mathrm{H}_{2}}^{\mathrm{eq}}$ (given by eqs. 12
and 24), we numerically solve for the virial temperature necessary for a gas
cloud to collapse in the presence of a fixed LW background intensity. In the
presence of a strong LW background, halos may be in the high mass limit of eq.
12, so we fit those independently and summarize the associated minimum masses
as follows. For the low-mass limit, we find

$\begin{split}M_{\rm LW,low}\simeq\ &3.35\times
10^{4}\bigg{(}\frac{\zeta}{0.25}\bigg{)}^{-0.19}\\\
&\times\Bigg{[}1+13\bigg{(}\frac{J_{\rm
LW}}{J_{21}}\bigg{)}^{0.38}\Bigg{]}\bigg{(}\frac{1+z}{21}\bigg{)}^{-1}M_{\odot}\end{split}\\\
$ (27)

At high masses, we find

$\begin{split}M_{\rm LW,high}\simeq\ &1.05\times
10^{4}\bigg{(}\frac{\zeta}{0.25}\bigg{)}^{-0.33}\\\
&\times\Bigg{[}1+26\bigg{(}\frac{J_{\rm
LW}}{J_{21}}\bigg{)}^{0.62}\Bigg{]}\bigg{(}\frac{1+z}{21}\bigg{)}^{-3.91}M_{\odot}\end{split}\\\
$ (28)

Noting that the transition for these two regimes in the halo density occurs at
$M_{\rm turn}$ (as defined in eq. 13), the photodissociation mass is set by

$\begin{split}M_{\rm LW}=\max\Big{[}\min\big{(}&M_{\rm LW,low},M_{\rm
turn}\big{)},\\\ &M_{\rm LW,high}\mathcal{H}\big{(}M_{\rm LW,high}-M_{\rm
turn}\big{)}\Big{]}\end{split}$ (29)

where $\mathcal{H}(t-t_{0})$ is the Heaviside step function (which returns 0
for $t<t_{0}$ and 1 for $t\geq t_{0}$). At a fixed redshift as a function of
increasing $J_{\rm LW}$, this means we select the low-density result until
that crosses the turnover mass, at which point we select the larger of the
turnover mass or the high-density result.

The results of this calculation for various redshifts are given in Figure 3.
Because this portion of the calculation does not account for the buildup of
$\mathrm{H}_{2}$ (eq. 17), we must still also compare this to the cooling
threshold presented in § 2.2 and check that a halo has enough $\mathrm{H}_{2}$
to efficiently cool (i.e., we choose the larger of eqs. 23 and 29). The
cooling threshold sets the baseline for the minimum mass for small values of
the LW intensity, with the turnover decreasing with redshift. Once the LW
background takes over, however, the critical mass increases with increasing LW
background intensity at a fixed $z$.

We note that the amplitude of the LW background will be computed self-
consistently by our semi-analytic model, as described in § 3.

Figure 3: The minimum mass for a halo to form $\mathrm{H}_{2}$ in the presence
of a Lyman-Werner background, calculated for a range of LW backgrounds
following the steps described in Sections 2.2 and 2.3 at various redshifts
(colored lines). The fitting formulae for the cooling and dissociation
thresholds (eqs. 23 and 29, respectively) are overlaid with dashed lines.

###  2.4 Thermalization threshold

The cooling argument of § 2.2 essentially requires that the halo produce
enough molecular hydrogen in order for its radiative cooling to shed the
thermal energy generated with gravitational infall. However, in the presence
of streaming, it must also shed the excess kinetic energy carried by the
baryons (which will also be thermalized during collapse). This effect can be
described as follows (Fialkov et al., 2012):

$T_{\rm vir}\mapsto T_{\rm vir}+\frac{\mu m_{p}(\alpha_{\rm vbc}v_{\rm
bc})^{2}}{2k_{B}}$ (30)

where $\alpha_{\rm vbc}$ parameterizes the magnification of the stream
velocity as gas collapses into the halo potential well and has been shown to
take on values of $\sim 4-6$ in simulations (Fialkov et al., 2012, McQuinn &
O’Leary, 2012).

We can understand the numerical value of $\alpha_{\rm vbc}$ by the following
two arguments (McQuinn & O’Leary, 2012). First, we consider the case of
accreting shock heated gas as it collapses into the halo. In adiabatic
collapse, $\rho\propto T^{3/2}$, so the overdensity of the virialized gas with
no stream velocity is

$\delta_{b}\sim\bigg{(}\frac{T_{\rm vir}}{T_{\rm IGM}^{\rm ad}}\bigg{)}^{3/2}$
(31)

The thermal energy associated with the stream velocity provides an extra
contribution to the IGM temperature, corresponding to $T_{\rm IGM}\mapsto
T_{\rm IGM}(1+5\mathcal{M}_{\rm bc}^{2}/9)$, where $\mathcal{M}_{\rm
bc}=v_{\rm bc}/c_{\rm s,IGM}$ is the mach number of the streaming gas.777This
can be seen if we assume that the energy in the stream velocity is converted
into thermal energy in the IGM. That is, the kinetic energy associated with
streaming gives $\mu mv_{\rm bc}^{2}/2=\mu mc_{s}^{2}\mathcal{M}_{\rm
bc}^{2}/2$. Using the expression for the sound speed of a monatomic ideal gas
and assuming some fraction $\beta$ of this energy is converted into thermal
energy (as seen in the simulations of O’Leary & McQuinn, 2012 and McQuinn &
O’Leary, 2012), then the temperature is increased from $T\to T^{\prime}$ by a
factor of $f\equiv T^{\prime}/T=1+5\beta\mathcal{M}_{\rm bc}^{2}/9$. We have
taken $\beta=1$ to compute the ballpark value above. With eq. 31, this means
that the gas density in the halo scales as $n_{\rm
H}\propto\delta_{b}\propto(1+5\mathcal{M}_{\rm bc}^{2}/9)^{-3/2}$. We can fold
this additional contribution to the gas density into equation 21 and calculate
the new virial temperature necessary for efficient cooling with the same
argument as that outlined in § 2.2. Combining this result with equation 30
suggests $\alpha_{\rm vbc}\sim 6$ at $z=20$. However, as the critical
$\mathrm{H}_{2}$ fraction, IGM temperature, and stream velocity all evolve
with redshift, the value of $\alpha_{\rm vbc}$ given by this calculation will
as well.

Next, we can consider that velocities will be boosted by a factor of
$(1+\delta)^{1/3}$ in an adiabatically collapsing region, so the effective
stream velocity will be magnified to $\approx 6v_{\rm bc}$ for overdensities
of $\delta\sim 200$. If the circular velocity of the halo is to be larger than
this local stream velocity, the circular velocity necessary for the high
velocity gas to be captured by the halo potential well is

$v_{\rm c}(v_{\rm bc})=\sqrt{v_{\rm c}^{2}+(6v_{\rm bc})^{2}/3}$ (32)

where we divide by a factor of 3 because the local stream velocity will have a
particular directionality, whereas the halo gas moves isotropically after
thermalization. Converting this velocity to the associated virial temperature
corresponds to a redshift-independent value of $\alpha\sim 4$.

These analytic arguments give us a ballpark value of $\alpha_{\rm
vbc}\sim\mathcal{O}(1-10)$ on average. For our fiducial calculations, however,
we defer to the more detailed numerical simulations of McQuinn & O’Leary
(2012) (who find $\alpha_{\rm vbc}\sim 4-6$) and therefore set $\alpha_{\rm
vbc}=5$.

In practice, we incorporate this into our model by adding this extra
contribution to the virial temperature (eq. 30) to the result of the
calculation in Sections 2.3 and 2.2. This contribution to the mass scale can
be parameterized as a modification to the virial velocity of the halo. That
is, for our chosen cosmology, the virial velocity of a halo without the
effects of streaming is given by (Barkana & Loeb, 2001)

$v_{0}=5.28\times 10^{-2}M_{0}^{1/3}\bigg{(}\frac{1+z}{21}\bigg{)}^{1/2}\ {\rm
km\ s^{-1}}$ (33)

where $M_{0}$ is the larger of eqs. 23 and 29. Then, the thermalization
threshold is given by

$M_{\rm bc,LW,cool}\simeq 6.79\times 10^{3}\Big{[}v_{0}^{2}+(\alpha_{\rm
vbc}v_{\rm
bc})^{2}\Big{]}^{3/2}\bigg{(}\frac{1+z}{21}\bigg{)}^{-3/2}M_{\odot}$ (34)

###  2.5 X-ray background

Unlike the Lyman-Werner dissociation and stream velocity, the buildup of a
metagalactic X-ray background both negatively and positively feeds back into
the minimum mass (and thus the star formation rate). A background of X-ray
radiation can photoionize the IGM, boosting the electron fraction and heating
the gas. The increased density of free electrons can then promote the
formation of $\mathrm{H}_{2}$, supporting star formation, while the heating
can make accretion of gas more difficult. Some energy from this radiation goes
into ionizing the hydrogen and helium atoms and the rest goes to the electron,
heating it. Following the discussion in Furlanetto (2006), the hot electron
distributes its energy by collisionally ionizing other atoms, producing
secondary electrons, collisionally exciting H and He, and undergoing Coulomb
collisions with thermal electrons. The cross sections of these processes set
the fraction of the energy that goes into heating and the fraction that goes
into ionization, for which we use the expressions of Shull & van Steenberg
(1985):

$\displaystyle f_{\mathrm{X,h}}$
$\displaystyle=C_{1}\Big{[}1-(1-x_{e}^{a_{1}})^{b_{1}}\Big{]}$ (35)
$\displaystyle f_{\rm X,ion}$ $\displaystyle=C_{2}(1-x_{e}^{a_{2}})^{b_{2}}$
(36)

Note that we test the effects of more accurately including secondary
ionizations (see e.g., Furlanetto & Stoever, 2010, Ricotti, 2016) and find
that they do not meaningfully change the results.

Therefore, calculating the effects of an X-ray background requires a
calculation of the IGM temperature and electron fraction. As described in
Furlanetto (2006), the X-ray contribution to the temperature and electron
fraction of the IGM can be written in terms of the X-ray emissivity
$\epsilon_{X}$ as follows:

$\displaystyle\frac{dT_{X}}{dz}$ $\displaystyle=\frac{2}{3}f_{\rm
X,h}\frac{\epsilon_{X}}{k_{B}n}\frac{dT_{X}}{dt}\frac{dt}{dz}=-\frac{2}{3}\left(\frac{f_{\rm
X,h}}{0.2}\right)\frac{\epsilon_{X}}{k_{B}nH(z)(1+z)}$ (37)
$\displaystyle\frac{dx_{e}}{dz}$ $\displaystyle=-f_{\rm
X,ion}\frac{\epsilon_{X}}{E_{H}nH(z)(1+z)}$ (38)

where $E_{H}=13.6\ {\rm eV}$ is the ionization energy of hydrogen. We
calculate a spectrum-averaged X-ray emissivity from the $L_{X}-{\rm SFR}$
relation, which we scale by a factor $f_{X}$ to account for deviations from
the local calibration (e.g., Mineo et al., 2012, Lehmer et al.,
2016)888Theoretical models (e.g., Mesinger et al., 2013, Fragos et al., 2013)
predict that the $L_{X}-{\rm SFR}$ relation will increase with redshift, with
$f_{X}\sim 10-50$ at $z\sim 10$ suggested by the results of Greig et al.
(2021) and The HERA Collaboration et al. (2022). An upper limit on the value
of $f_{X}$ is placed by the Chandra Deep Field-South observations of the
unresolved soft X-ray background (SXB; Hickox & Markevitch, 2007);
$f_{X}\gtrsim 100$ would saturate the SXB (McQuinn, 2012), so we bound our
models by testing values between $f_{X}\sim 1-100$.:

$\epsilon_{X}(z)=2.6\times 10^{39}f_{X}\left(\frac{\rho_{\rm
SFR}(z)}{M_{\odot}\mathrm{yr}^{-1}\mathrm{cMpc}^{-3}}\right)\mathrm{erg}\
\mathrm{s}^{-1}\mathrm{cMpc}^{-3}$ (39)

where $\rho_{\rm SFR}$ is the star formation rate density computed from our
model.

Initializing the calculation with the IGM temperature and electron fraction
calculated with CosmoREC, we integrate the full evolution equations (given in
Furlanetto, 2006) simultaneously, yielding the histories displayed in Figure
4.

Figure 4: The evolution of the IGM temperature (top) and electron fraction
(bottom) with redshift for various values of $f_{X}$ using the star formation
history computed with the fiducial parameters of Mebane et al. (2018). We
compare our calculation to that of CosmoREC (i.e., no X-ray heating) with a
black dot-dashed curve and show the CMB temperature with a black dotted curve.

With these in hand, we compute the X-ray contribution to the minimum mass
through its effects on the previously discussed quantities: the filter mass
and LW mass. To incorporate the effects of heating, we modify the temperature
evolution of the IGM used in the calculation of the filter mass (see § 2.1).
This will have the effect of raising the filter mass at late times as the IGM
temperature rises steeply for $z\lesssim 15$ (Figure 4). Ionization raises the
free electron fraction of the IGM. The free electrons catalyze the formation
of $\mathrm{H}_{2}$, hastening collapse and lowering the minimum mass scale
derived from LW photodissociation (see § 2.3).

However, because of the higher density and thus recombination rate within a
halo, the relevant electron fraction within a DM halo will be smaller than the
IGM value that we calculate from integrating equation 38. We can estimate the
electron fraction in a halo by appealing to a similar argument as that
outlined in § 2.2. That is, if we assume that $n_{\rm HII}=n_{\rm e}$, then
eq. 14 will describe the evolution of $x_{\rm e}$ as well. This then admits
the same solution as before; that is:

$x_{\rm e,rec}=\frac{x_{\rm e}^{\rm IGM}}{1+t_{H}/t_{\rm rec}^{i}}$ (40)

where $t_{\rm rec}$ is defined as before and we have taken $t=t_{H}$ to be the
relevant timescale for comparison here.

Therefore, there are two relevant limits for this calculation. If the
recombination time is long compared to the Hubble time (so $t_{\rm rec}\gg
t_{H}$), then the cloud will not reach equilibrium and thus the electron
fraction in the halo will be given by eq. 40. If the recombination time is
sufficiently short compared to the Hubble time (i.e., $t_{\rm rec}\ll t_{H}$),
then the solution (eq. 40) no longer applies—because the ionizing background
could raise the electron fraction during that time—and we must instead
calculate the equilibrium level of ionization.

For the equilibrium calculation, we balance ionization and recombination,
solving

$\alpha_{B}n_{e}n_{p}=\Gamma_{X}n_{\rm HI}$ (41)

where $\Gamma_{X}$ is the ionization rate. Because we do not directly model
the X-ray source spectrum (i.e., $\epsilon_{X}$ is taken to be frequency-
independent), the ionization rate can be found by solving the cosmological
radiative transfer equation (Haardt & Madau, 2012):

$\Gamma_{X}(z)=c\int_{\nu_{\rm thresh}}^{\infty}d\nu\frac{\sigma_{\rm
HI}(\nu)}{h\nu}\int_{z}^{\infty}dz^{\prime}\frac{(1+z)^{3}}{H(z^{\prime})(1+z^{\prime})^{4}}\epsilon_{X}(z^{\prime})e^{-\tau}$
(42)

where the threshold frequency $h\nu_{\rm thresh}=13.6\ {\rm eV}$,
$\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$, and the optical depth is
$\tau=\int_{z}^{z^{\prime}}(d\ell/dz^{\prime\prime})dz^{\prime\prime}[n_{\rm
HI}\sigma_{\rm HI}(\nu^{\prime\prime})+n_{\rm HeI}\sigma_{\rm
HeI}(\nu^{\prime\prime})]$. For the cross sections $\sigma_{\rm HI}$ and
$\sigma_{\rm HeI}$, we use the fits of Verner et al. (1996).

As we evolve the X-ray background throughout our calculation, we fold the
newfound halo electron fraction into our calculation of the LW mass (in eq.
24), and expect that it will suppress the mass scale as $x_{\rm e,halo}$
increases.

Including the effects of X-rays in a simple fitting formula is comparatively
more difficult than the preceding effects. In particular, to incorporate the
effects of heating, one must integrate over the full temperature history to
accurately estimate the filter mass. Therefore, we only present a simple
modification to the fitting formulae for the photodissociation threshold to
account for the effect of X-ray ionization:

$\displaystyle M_{\rm LW,Xray,low}$ $\displaystyle\simeq M_{\rm
LW,low}\bigg{(}\frac{x_{\rm e,halo}}{x_{\rm e,no\ Xray}}\bigg{)}^{-0.19}$ (43)
$\displaystyle M_{\rm LW,Xray,high}$ $\displaystyle\simeq M_{\rm
LW,high}\bigg{(}\frac{x_{\rm e,halo}}{x_{\rm e,no\ Xray}}\bigg{)}^{-0.33}$
(44)

where $M_{\rm LW,low}$ and $M_{\rm LW,high}$ are given by equations 27 and 28,
and $x_{\rm e,no\ Xray}$ is the electron fraction of the IGM without the
effects of an X-ray background.999For this calculation, one can refer to a
recombination code such as RECFAST or CosmoREC (Seager et al., 1999, Chluba &
Thomas, 2013). For the relevant redshifts for Pop III star formation (i.e.,
between $z=5$ and $50$), however, the electron fraction without X-rays is well
fit by a power-law of the form $x_{\rm e,no\ Xray}\approx 2.19\times
10^{-4}[(1+z)/21]^{0.12}$. From these formulae, we can see that an increase in
the electron fraction from the IGM electron fraction without the effects of an
X-ray background will indeed suppress the dissociation mass, though the
strength of this contribution will depend on the density regime we are in.

###  2.6 Summary

Our minimum mass model can be summarized as follows. A halo must first exceed
the filter mass threshold (§ 2.1; eq. 8) in order to be able to accrete
baryons. Once this gas is accreted, a sufficiently large fraction of
$\mathrm{H}_{2}$ must build up such that the gas can efficiently cool and
collapse to high densities (§ 2.2; eq. 23). If the halo is in a region of the
universe with a large DM-baryon relative velocity, the halo must also overcome
the additional thermal energy associated with the relative motion (§ 2.4; eq.
34). As the first generations of stars form, they produce a photodissociating
LW background (parameterized in terms of the specific intensity $J_{\rm LW}$)
and a photo-heating and -ionizing X-ray background (which increases the
temperature of the IGM and the electron fraction of the halo $x_{\rm
e,halo}$). Strong LW intensities and efficient X-ray heating can significantly
boost the minimum mass (§ 2.3 and 2.5; eq. 29), while ionization can promote
the formation of $\mathrm{H}_{2}$ and lower the mass (eq. 43).

In practice, this means the minimum mass is given by

$M_{\rm min}=\max[M_{F},M_{\rm bc,LW,cool}]$ (45)

where $M_{F}$ is given by eq. 8 and $M_{\rm bc,LW,cool}$ is given by eq. 34
(which depends on eqs. 23, 29, and 43). These fitting formulae presented in
the preceding several sections are accurate to within 40% of the full
calculation, and we find that using them does not change the results of the
semi-analytic calculation (if we do not include the effects of IGM heating).
IGM heating appears to be important at the latest times and in the case of the
strongest X-rays; that is, when $f_{X}\gtrsim 50-100$. This contribution to
the minimum mass is discussed in more detail in § 4.2.

Calculating the minimum mass and understanding how and when individual
components of the physics are important requires robust characterization of
the star formation histories and associated radiation backgrounds. As such, we
now explore the effects of this minimum mass model in the context of a semi-
analytic model for high-redshift star formation.

##  3 Semi-analytic Model

We base our calculations on an updated version of the semi-analytic model
presented in Mebane et al. (2018) (hereafter, M18). This model employs a
feedback-limited star-formation prescription, wherein star formation is
tracked in individual halos until feedback shuts off the formation process. We
follow a sample of 100 halos from $z=50$ to 6 with masses ranging from
$10^{6}-10^{13}M_{\odot}$ at $z=6$ with timesteps of 1 Myr. The growth history
of these halos is calculated using abundance matching, where we assume that
halos maintain their comoving number density through cosmic time, with the
halo mass function given in Trac et al. (2015). Star formation is followed in
individual halos by comparing the halo mass to the evolving minimum mass
outlined in § 2. Once a halo passes this star formation threshold, we compare
its gas mass (which is given by $M_{\rm gas}=\Omega_{b}/\Omega_{m}M_{\rm h}$
until supernovae evacuate some fraction of this gas) to the local Jeans mass.
If the gas mass exceeds the local Jeans mass, it will randomly form one or two
stars (determined by a fixed binary probability $f_{\rm bin}=0.5$) with masses
drawn from a Chabrier-like IMF (Chabrier, 2003):

$\frac{dN}{dm}\propto M^{-\alpha}\exp\Bigg{[}-\bigg{(}\frac{M_{\rm
char}}{M}\bigg{)}^{\beta}\Bigg{]}$ (46)

where $\alpha$ is the Salpeter-like power-law slope that characterizes the
high-mass end of the IMF and $M_{\rm char}$ is the characteristic mass. These
massive stars live for 5 Myr and their end-of-life behavior is determined by
their stellar mass; i.e., stars with $M/M_{\odot}\in[8,40]\cup[140,260]$ end
their lives in SNe and the rest collapse directly to BHs.

Because Pop III stars form in minihalos at high redshifts, supernova feedback
can significantly affect their environments through metal enrichment and gas
ejection. When a supernova occurs in a halo, the circumstellar gas is enriched
according to the metal yields given by Heger & Woosley (2002) and Heger &
Woosley (2010) for pair-instability and core-collapse SNe, respectively. We
then assume that 10% of the released kinetic energy couples to gas in the halo
and ejects it. The mass and metallicity of the ejected gas is then followed
and allowed to reaccrete after a free-fall time (roughly 50 Myr for a
$10^{5}M_{\odot}$ minihalo at $z=15$; see §6.4 in Mebane et al., 2018 for a
detailed discussion of this choice). In the meantime, the halo continues to
grow via accretion from the pristine IGM, enabling the potential formation of
subsequent generations of Pop III stars. We note that the ability for halos to
form multiple generations of Pop III stars is a unique feature of our model
that we discuss in more detail in § 5.2. The simulations of Abe et al.
(2021)—which modeled the effects of SN feedback and metal enrichment on Pop
III star formation—found that the delay in reaccretion of enriched gas allowed
for multiple generations of Pop III star formation before a halo crossed the
critical metallicity and transitioned to the Pop II regime, which motivates
our choice here.

The transition to Pop II star formation is an important piece of our model, as
Pop II stars will dominate in the buildup of the metagalactic LW and X-ray
backgrounds. A halo will transition to Pop II star formation once it has
reached the atomic cooling threshold ($T_{\rm vir}>10^{4}$ K) or has been
sufficiently enriched (i.e., the mean metallicity of the halo exceeds the
critical CII and OI concentrations given in Bromm & Loeb, 2003). We note that
we assume that atomic-cooling haloes form Pop II stars without attempting to
model metal mixing in those haloes, but even if metal mixing is inefficient,
their rapid radiative cooling will modify the process of star formation so
that $\mathrm{H}_{2}$ is no longer required. Any Pop III stars forming in
massive haloes are thus likely to be more similar to Pop II stars than those
forming inside minihaloes.

Once a halo transitions, we use the bursty star formation prescription
outlined in Furlanetto & Mirocha (2022). This modifies the feedback-regulated
models of Furlanetto et al. (2017) to include a feedback delay that accounts
for the nonzero lifetimes of SN progenitors. That is, rather than
instantaneously injecting SN feedback into the system after stars form, we
include a 5-30 Myr delay (accounting for stellar lifetimes) that allows for
continued star formation in the meantime. As a result, the star formation rate
can ‘overshoot’ the expectation of a simple feedback-limited model and halos
will ultimately go through cycles of star formation in characteristic
‘bursts.’ This effect is most pronounced in the least massive halos, which are
the most susceptible to the effects of feedback. As a halo grows, these star
formation cycles damp out and the star formation efficiency approaches the
equilibrium result used in M18 (originally presented in Furlanetto et al.,
2017).

We use the aforementioned framework to evolve the system of halos between
redshifts 50 and 6. We extrapolate these results (which are calculated for a
sample of 100 halos) to calculate global quantities, such as the star-
formation rate density (SFRD), by averaging over the halo mass function. To
minimize noise in our calculations (because we track only 100 halos for
computational efficiency), at every timestep, we inject 10,000 ‘fake’ halos
into our calculation. These halos form stars with a probability given by the
star formation duty cycle $f_{\rm duty}(z)=N_{\rm III,on}(z)/N_{\rm III}(z)$,
where $N_{\rm III,on}(z)$ is the number of halos that are actively forming
stars and $N_{\rm III}(z)$ is the number of halos that are able to form Pop
III stars at a given $z$ (i.e., have yet to transition to Pop II star
formation). We compare this procedure to calculating the SFRD with 1,000
‘real’ halos (compared to the 100 in our fiducial runs) and find that the
‘fake’ halo injection robustly reproduces the expected SFRD.

Figure 5: The minimum DM halo mass for Pop III star formation as a function of
redshift in our fiducial model, which has $v_{\rm bc}=\sigma_{\rm vbc}$ and
$f_{X}=10$ (solid black curve) and the associated fitting formula (dashed
black curve). We present the individual components of the calculation
following § 2: the cooling threshold (red), the accretion threshold (green),
and the photodissociation threshold (blue). For the cooling and dissociation
thresholds, the equivalent calculation with $v_{\rm bc}=0$ is given as a
dashed curve of the same color. The dissociation threshold without the
calculation of the X-ray background (§ 2.5) is given in orange.

##  4 Results

We now use the updated prescription for the minimum mass (§ 2) and the
improvements to the semi-analytic framework (§ 3) to follow star formation
during the Pop III era. For the fiducial model:

  1. (i)

In calculating the minimum mass, we assume a stream velocity magnitude of
$v_{\rm bc}=1\sigma_{\rm vbc}$ and an ${\rm L}_{X}-{\rm SFR}$ scaling of
$f_{X}=10$.

  2. (ii)

We assume a Chabrier IMF (equation 46) with a maximum mass of $500M_{\odot}$,
$\alpha=2.35$ (Salpeter, 1955), $\beta=1.6$, and $M_{\rm char}=20M_{\odot}$.

  3. (iii)

We include bursty star formation for Pop II halos.

We note that our fiducial value for the magnitude of the stream velocity
($1\sigma_{\rm vbc}$) was chosen for ease of comparison with existing work.
The stream velocity follows a Maxwell-Boltzmann distribution, for which the
mean is $\sqrt{8/3\pi}\sigma_{\rm vbc}\sim 0.92\sigma_{\rm vbc}$ and the most
probable value is $\sqrt{2/3}\sigma_{\rm vbc}\sim 0.82\sigma_{\rm vbc}$. We
find that choosing one of these values would lower the minimum mass at early
times by a factor of a few but would be qualitatively similar to the
$1\sigma_{\rm vbc}$ case. Our fiducial choice of $f_{X}=10$ is motivated by
the findings of The HERA Collaboration et al. (2022) (and the predictions of
Fragos et al., 2013). We also present variations of many of these parameters
around these fiducial values.

Figure 6: (top) The minimum DM halo mass for Pop III star formation as a
function of redshift in our fiducial model ($v_{\rm bc}=1\sigma_{\rm
vbc},f_{X}=10$; solid black curve) and the previous incarnation of the semi-
analytic model (M18; dashed black curve). (bottom) The star formation rate
density (SFRD) as a function of redshift comparing our fiducial model (solid)
to M18 (dashed). Note that we are using a different fiducial parameter set
from M18 for this comparison. The Pop III (II) contribution to the SFRD is
given in red (blue) and the total SFRD is given in black.

In Figure 5, we present the minimum mass calculated using the radiation
backgrounds computed with our fiducial parameter set (see lower panel of
Figure 6 for the associated SFRD). We can understand the evolution of the
minimum mass with time by inspecting the evolution of the individual
contributions to the minimum mass (colored curves) in various epochs.

  1. (I)

Streaming dominance, $z\gtrsim 30$: At the earliest times, the stream velocity
is the dominant limiting factor and the shape of the minimum mass traces the
evolution of the cooling threshold for a region of 1$\sigma_{\rm vbc}$
streaming (red). The comparison between the cooling threshold with and without
streaming (solid and dashed red curves, respectively) demonstrates the
significant effect of the stream velocity on the minimum mass (nearly two
orders of magnitude magnification). During this time, the minimum mass decays
with the stream velocity as $v\propto(1+z)$ and the calculations with and
without streaming (dashed vs solid curves) rapidly converge.

  2. (II)

LW dominance, $10\lesssim z\lesssim 30$: By $z=30$, the SFRD has grown
(largely dominated by Pop III stars; the solid red curve in the lower panel of
Figure 6) and the LW background has built up, so the LW mass (blue) takes over
and the minimum mass begins to grow. Coincident with this is an almost plateau
in the Pop III SFRD, as the minimum mass is now growing at a rate comparable
to the growth rate of low mass halos. This redshift evolution of the LW mass
is tied to the efficiency of self-shielding in massive halos. In other words,
because of the steep decrease in halo gas density as the universe expands
(i.e., $n_{\rm gas}\propto(1+z)^{3}$), self-shielding is less efficient at
lower $z$ and the LW mass can climb steadily as a result.

  3. (III)

Bursty star formation, $7\lesssim z\lesssim 15$: At later times, Pop II star
formation dominates and the effect of the bursty SF cycles is apparent in both
the SFRD and the minimum mass. As a result, the Pop II SFRD overshoots the
equilibrium case (black dashed curve in lower panel of Figure 6) and induces a
slight increase in the minimum mass, as we expect from the results of
Furlanetto & Mirocha (2022). We discuss this effect in more detail in § 4.3.

  4. (IV)

X-ray era, $z\lesssim 15$: At the latest times, we can compare the positive
and negative effects of the X-ray background. The dominant effect of the X-ray
background is the positive feedback from photoionization, which suppresses the
minimum mass more efficiently at lower $z$ (orange vs. black curve) as the
universe becomes ionized (Figure 4). With $f_{X}=10$ in our fiducical model,
the negative effect of X-ray heating is subdominant and is manifested as an
upturn in the filter mass at the latest times. Because of this minimum mass
suppression (from ionization), Pop III star formation persists through the end
of our calculation to $z\sim 6$.

While the redshifts of transition between the different effects will vary with
$v_{\rm bc}$ and $f_{X}$, this evolution of the minimum mass is qualitatively
generic to all models, as is discussed below.

###  4.1 Comparison to M18

Comparing our fiducial minimum mass to that of M18 during this epoch, we can
isolate the key differences that result from the updated calculations
described in § 2. The M18 model used a prescription for the minimum mass fit
to the results of simulations (O’Shea & Norman, 2008) that did not include the
effects of self-shielding. As the Lyman-Werner background builds up, we see a
nearly order of magnitude suppression in the mass scale. The inclusion of
self-shielding, following the updated criterion of WG11/WG19 (eq. 26),
introduces more efficient shielding in the densest halos, suppressing the
effects of the photodissociating background. Coupled with the positive effect
of X-ray ionization,101010The minimum mass prescription used in M18 does not
include the effects of an X-ray background — without this positive feedback,
the minimum mass would be increased by a factor of a few at late times; see §
4.2. this has the effect of allowing a sustained and slowly growing level of
Pop III star formation down to the latest times, when the SFRD begins to
decline. Moreover, the peak Pop III SFRD is two orders of magnitude larger in
the new model, largely because of the self-shielding in massive halos. As a
result, Pop III star formation will contribute significantly for a much longer
time. The difference in the shapes of the minimum mass curves at $z\gtrsim 25$
is due to the inclusion of the stream velocity in the calculation of the
cooling threshold. However, because of the additional requirement that the
halo gas mass must exceed the local Jeans mass in order for star formation to
begin, the nearly order of magnitude increase in the minimum mass only
slightly filters into the global SFRD.

In Figure 7, we compare the calculation of the Pop III SFRD using our fiducial
cosmology to that of M18 and consider some individual parameter variations to
diagnose the differences. Disentangling the sources of discrepancy induced by
changes in the overall cosmology is challenging, but it is clear that the
early phases of Pop III star formation are sensitive to this choice. Once the
radiation backgrounds build up and the SFRD reaches its peak, however, the
differences are washed out and the various models converge to the same result.

For a more detailed discussion of the star formation process in individual
minihalos we refer the reader to the analysis in M18, for which our model is
analogous. Here, we instead now turn to an exploration of the variations in
the global SFRD that result from variations in the individual star formation
physics.

Figure 7: The redshift evolution of the Pop III SFRD for different choices of
the underlying cosmology (different colors). Our fiducial choice is given in
black, the cosmology used in M18 in orange, and variations of one or a few
parameters from our fiducial choice in the other colors. Figure 8: (top) The
minimum mass for three different values of the stream velocity (different
colors), comparing our calculation (solid curves) to that based on Nebrin et
al. (2023) (dashed curves; i.e., using their density threshold, taking
$\zeta=0.16$, and setting $\alpha_{\rm vbc}=6$). (bottom) The associated Pop
III SFRD for the same cases as above. Figure 9: (top) The minimum mass as a
function of redshift for $f_{X}=0,1,10,100$ (green, black, blue, red,
respectively) and $v_{\rm bc}=0,1,3\sigma_{\rm vbc}$ moving from left to
right. The Jeans mass (effectively the minimum halo mass for the first star to
form; dashed black curve) and atomic cooling threshold ($T_{\rm vir}=10^{4}\
{\rm K}$; dotted black curve) are shown for reference. (bottom) The associated
Pop III SFRD for the same cases as the upper panels.

###  4.2 The effects of the minimum mass on Pop III star formation

The minimum mass model outlined in § 2 allows us to naturally test the effects
of varying the details of Pop III star formation physics. To this end, we
examine the evolution of the minimum mass and Pop III SFRD for representative
values of the stream velocity and X-ray background strength (which is
parameterized by $f_{X}$; see equation 39) in Figures 8 and 9. Note that in
these figures, we have omitted the Pop II SFRD because the timing of the
transition is effectively independent of the details of Pop III star formation
(and our model assumes these halos are not subject to the same feedback
mechanisms)—in other words, most halos transition to Pop II star formation
once they cross the atomic cooling threshold, rather than as a result of metal
enrichment. 111111For example, a $10^{13}M_{\odot}$ halo at $z=6$ will have
crossed the atomic cooling threshold at $z\sim 43$, while for
$10^{10}M_{\odot}$ and $10^{8}M_{\odot}$ halos, that crossing point is $z\sim
23$ and $z\sim 10$, respectively. We have also zoomed in on the times during
which the effects of X-ray feedback will be dominant ($z\leq 25$).

First, we can isolate the effects of introducing and varying the stream
velocity strength with no X-ray background (i.e., $f_{X}=0$; Figure 8). We
also demonstrate the effects of variations in our choice of halo gas density
by overlaying the evolution of the minimum mass and SFRD using the fiducial
parameters adopted in Nebrin et al. (2023) (i.e., their density prescription,
$\zeta=0.16$, and $\alpha_{\rm vbc}=6$; dashed curves). However, we defer
discussion of those differences to § 5.1 and focus on variations in $v_{\rm
bc}$ and $f_{X}$ here.

Broadly, the minimum mass is dominated by the cooling mass at early times and
the joint LW-X-ray mass once star formation has built up for all values of
$v_{\rm bc}$. As we increase the stream velocity, however, the filter and
cooling masses grow in strength at early times and the minimum mass
demonstrates a qualitatively different shape as a result (as it decays with
the stream velocity $v_{\rm bc}\propto(1+z)$). Increasing the stream velocity
from $0\to 1$ or $1\to 2\sigma_{\rm vbc}$ induces an increase in the minimum
mass by a factor of $\sim 10-20$. This drastic increase in the minimum mass is
clearly reflected in the Pop III SFRD for regions with a moderate or strong
relative velocity ($v_{\rm bc}\gtrsim 2\sigma_{\rm vbc}$). Note that the
effect of this in the SFRD is diminished somewhat in going from $0$ to
$1\sigma_{\rm vbc}$ because of the additional requirement that a halo must
build up enough gas to exceed the local Jeans mass (dashed black curve in
Figure 9), which takes a significant time, thereby raising the effective
minimum mass for the first generation of star formation. Nevertheless, the
first stars do not form until $z\sim 30-35$ in regions of the strongest
relative velocity (compared to $z\sim 40-45$ in weak streaming regions)—a
nearly 50 Myr delay in the onset of star formation. However, despite this
delay, once minihalos begin to form Pop III stars, the SFRD quickly approaches
the ‘equilibrium’ level seen in the no/weak streaming case, so the late time
behavior is indistinguishable from the weak streaming limit in most cases. In
the case of strong relative velocity (yellow curves), the minimum mass is
already growing steeply by the time the SFRD builds up to a high level, so it
is not able to reach the same peak level as in regions with $v_{\rm
bc}\lesssim 2\sigma_{\rm vbc}$.

In Figure 9, we isolate the effects of X-ray feedback at fixed $v_{\rm bc}$.
First, across all three panels, it is evident that the inclusion of X-ray
feedback significantly enhances Pop III star formation. That is, moving from
no X-ray feedback (green) to even a weak background (black) induces a factor
of a few change in the minimum mass and flattens the redshift evolution
somewhat. Indeed, the clearest difference between the runs with and without
X-ray feedback are the duration of high levels of Pop III star formation—when
X-ray suppression of the minimum mass is not present, Pop III star formation
is quenched by $z\sim 7-8$.

Across all three streaming runs, both the positive and negative effects of
X-ray feedback are apparent in the minimum mass. During the epoch where the
minimum mass is set by LW feedback, positive feedback is manifested as a
suppression to the minimum mass but the effect is relatively weak—the minimum
mass is lowered by a factor of $\sim 2-3$ over two orders of magnitude
increase in $f_{X}$. Despite this seemingly weak effect, the peak level of
star formation achieved is boosted as we increase $f_{X}$. At the latest
times, however, with $f_{X}=100$ (red), the X-ray background has built up
sufficiently for the IGM heating to translate into a steep increase in the
filter mass, which can be seen in the smooth upturn of the minimum mass curves
at $z\sim 5-10$, where it switches back from tracing the LW or cooling curves
to again following the filter mass. The timing of this transition is dependent
on the strength of the X-ray background; with $f_{X}\gtrsim 50$ the IGM is
efficiently heated (Figure 4) after $z\sim 15$ and the minimum mass rises
steeply in response. This negative feedback is reflected in the Pop III SFRD
for the strongest X-rays—as heating overtakes the suppression due to enhanced
positive X-ray feedback, $M_{\rm min}$ approaches the same late time value
($\sim 10^{8}M_{\odot}$) as the no X-ray model and Pop III star formation is
again shut off at the latest times, though the transition happens more sharply
than the no X-ray case.

Analyzing the effects of $v_{\rm bc}$ and $f_{X}$ together, it is apparent
that the independent variations in the minimum mass and SFRD persist in the
presence of both effects. That is, as we increase $v_{\rm bc}$, the peak SFRD
is increasingly suppressed and as we increase $f_{X}$, it is boosted. In the
case of the strongest X-rays and strong streaming (red curve in rightmost
panels), X-ray heating is significant at nearly the same time that the LW mass
takes over from the cooling threshold. The steeper growth of the minimum mass
in this case results in the SFRD plummeting just as quickly as it grew.

We note that, regardless of $f_{X}$ and the local stream velocity, all of
these models yield a peak Pop III SFRD $\sim 10^{-3}\ M_{\odot}$ yr-1 Mpc-3.
These parameter choices primarily determine the duration of that peak era.
However, the peak SFRD will be sensitive to other choices in the model; for
example, we assume that each halo forms only one or two stars at a time. The
SFRD will be directly proportional to the average mass of stars formed in each
event.

###  4.3 The effects of bursty Pop II star formation

In Figure 10 we explore the effects of bursty Pop II star formation on the
global Pop III and II SFRD. For clarity in isolating the effects of bursty
star formation, we present the SFRD calculated for the minimum mass
prescription used in M18. As described in § 3, cycles of bursty star formation
should cause Pop II halos to overshoot the equilibrium SFR, as is reflected in
the two red curves. The signature of these cycles is most clearly seen at low
redshifts (when our 1 Myr timesteps are much more finely spaced in redshift
space) but they affect the first Pop II halos as well. An increase in the Pop
II SFRD leads to a quicker buildup of the LW background and thus an increase
in the minimum mass. This in turn has the effect of suppressing Pop III star
formation—the Pop III SFRD plateau achieved in the bursty case is smaller than
that of the equilibrium case and star formation shuts off at a redshift of
$\sim 16\ (12)$ in the bursty (equilibrium) case. In this sense, bursty Pop II
star formation acts as another negative ‘feedback’ mechanism for Pop III star
formation.

However, we include bursty star formation in our fiducial model and in the
parameter variations shown in Figure 9 and find that the negative effects are
subdominant relative to the other effects setting the minimum mass. That is,
the inclusion of X-ray ionization and updated models for self-shielding, for
example, suppress the minimum mass to a level that bursty star formation in
Pop II halos does not meaningfully affect the ability of a minihalo to form
Pop III stars.

Figure 10: The Pop III (black) and Pop II (red) SFRD as a function of redshift
comparing the bursty (solid) and equilibrium (dotted) models for the minimum
mass prescription used in M18.

##  5 Comparison to other works

Several other groups have considered the question of when and where Pop III
stars can form, both from a numerical and analytic standpoint. Because this
epoch of the universe’s history has yet to be observed, the uncertainties are
large and it is useful to consider a range of models that explore the entire
parameter space. In this section we compare our results to those obtained by
other groups, first focusing on the minimum mass calculation and next on the
semi-analytic model.

###  5.1 Minimum mass model

The minimum halo mass for molecular hydrogen cooling has been studied in a
number of contexts. One of the first attempts to analytically describe early
star formation came in Tegmark et al. (1997), wherein the collapse threshold
argument was introduced to identify when $\mathrm{H}_{2}$ cooling will become
efficient. We base the broad structure of our photodissociation criterion—see
§ 2.3—on this argument. This argument was shown to produce minimum masses that
qualitatively agree with the results of simulations (e.g., Machacek et al.,
2001, O’Shea & Norman, 2008, Visbal et al., 2014b) though in detail the
results differ by a factor of a few. The effect of streaming on gas accretion
has been studied many times as well; we have already compared our results to
Naoz et al. (2013) in § 2.1.2, demonstrating good quantitative agreement.

The first simulations to study the _joint_ effects of LW photodissociation and
the stream velocity on the minimum mass were K21 and S21. In Figure 11, we
compare our calculations to the simulated results of K21 and S21, isolating
the effects of photodissociation and streaming. These two simulations tested
various values for the LW background and stream velocity with the goal of
estimating the minimum DM halo mass for star formation in the presence of
these two effects. However, their simulations probe different regions of the
LW parameter space—the K21 group used values of $J_{\rm LW}/J_{21}=0,1,10,30$
and S21 used $J_{\rm LW}/J_{21}=0,10^{-2},10^{-1}$—and are discrepant in their
limited overlap. Therefore, we remain agnostic about committing to either set
of results and defer to our analytic calculation of the minimum
mass.121212During the era of LW dominance in the minimum mass ($z\lesssim 25$;
see Figure 6 and § 4), the LW intensity is of order $J_{\rm
LW}/J_{21}\sim\mathcal{O}(1-10)$, so the K21 results likely probe the more
relevant parameter space for regulating Pop III star formation. We find that
we can more closely reproduce their results with a modification of $\zeta$ in
our cooling threshold, perhaps in part because of differences in the
definition of the critical (K21) or average (S21) mass thresholds. In
particular, for the following comparison, we take $t_{\rm cool}<t_{\rm
ff}=\zeta t_{H}\implies\zeta=0.16$.

Figure 11: The redshift evolution of the minimum mass (solid colored curves)
for various values of $J_{\rm LW}$ with no streaming (left) and for various
values of $v_{\rm bc}$ with no LW background (right). These panels are the
equivalent of Figures 2 and 3 in K21 and Figure 5 in S21. The simulated data
from K21 (pentagons) and S21 (triangles) are overlaid in each case. Note that
we have omitted the error bars on these data for readability. For $J_{\rm
LW}<0.1J_{21}$, the cooling mass is larger than the photodissociation
threshold for all $z\geq 10$, so the blue and green curves lie atop each other
in the left panel.

In Figure 11, we show that our analytic calculations agree reasonably well
with the K21 simulations and underpredict the S21 results in the case of no
streaming and no LW background (where the mass scale should be set by the
$v_{\rm bc}=0$ cooling threshold described in § 2.2). This discrepancy
continues as we introduce a nonzero LW intensity and stream velocity.

First, we fix $v_{\rm bc}=0$ and isolate the effects of varying $J_{\rm LW}$
(left panel). Comparing the simulated results, we see that those of S21
demonstrate a larger minimum mass for $J_{\rm LW}/J_{21}=0.1$ than those of
K21 for $J_{\rm LW}/J_{21}=1$. We find that the redshift evolution of our
results for $M_{\rm crit}$ agrees qualitatively well with those of K21,
especially for large $J_{\rm LW}$. For $J_{\rm LW}/J_{21}=1$, however, our
calculations differ from that of K21 by a factor of a few. One key difference
in the simulations of S21 and K21 is in their calculation of $f_{\rm shield}$.
In particular, K21 use the same self-shielding criterion as in this work
(WG11/19), which was introduced as a correction to that used in S21 (namely,
Draine & Bertoldi, 1996) for warm, dense gas. Previous calculations of $M_{\rm
min}$ that ignored the effects of self-shielding (e.g., Machacek et al., 2001,
O’Shea & Norman, 2008) predict minimum masses nearly an order of magnitude
larger than those found in this work. We find that the minimum mass is fairly
sensitive to the magnitude and redshift evolution of the density threshold
which in turn feeds into the strength of the self-shielding effect and
speculate that this is perhaps one source of the discrepancy between K21 and
S21.

Next, we fix $J_{\rm LW}=0$ and isolate the effects of streaming. We find that
the calculation of the thermalization threshold (§ 2.4) with our choice of
$\alpha_{\rm vbc}=5$ produces excellent agreement with the simulated results
of K21. For $v_{\rm bc}\neq 0$, our calculations underpredict those of S21 by
a factor of a few, with the discrepancy increasing with increasing $v_{\rm
bc}$. However, given that K21 and S21 disagree by similar amounts, it is not
clear how to resolve this difference.

Despite the small quantitative differences noted in the evolution of the
minimum mass displayed in Figure 11, we find that the resulting Pop III SFRD
changes only modestly if we replace our minimum mass model with the fitting
formula reported in K21.

On the other hand, we find that the large variations in $M_{\rm crit}$ with
$v_{\rm bc}$ reported by S21 do significantly change the Pop III SFRD, if we
use their fitting formula. For example, with $v_{\rm bc}=3\sigma_{\rm vbc}$
their procedure makes the minimum mass so large that no Pop III stars form
(cf. the star formation delay we see in Figure 9).

The impact of an X-ray background on the Pop III star formation process is
similarly uncertain. Ricotti (2016) studied the joint effects of LW
photodissociation and X-ray feedback analytically, using a cooling and Jeans
mass criterion to set the minimum mass. They found that this criterion defined
a global feedback loop, wherein Pop III stars can positively and negatively
feed back into the star formation process, with the dominant effect depending
on the strength of the X-ray background. Both of their conditions for defining
the minimum mass are similar in spirit to ours but differ in detail, due to
differences in the details of our cooling criterion, X-ray calculations, and
our use of the filter mass rather than the Jeans mass. Qualitatively, however,
we find a similar feedback loop wherein a critical value for the X-ray
strength ($f_{X}\gtrsim 40$ for our model) defines the transition between
positive and negative feedback dominance. The simulations of Park et al.
(2021) follow up on this analytic calculation to investigate the effects of
incorporating local feedback processes in more detail and find broad
qualitative agreement with the results of Ricotti (2016).

Figure 12: A comparison of the redshift evolution of the minimum mass for
three representative values of the stream velocity ($0,1,2\sigma_{\rm vbc}$ —
purple, black, and orange, respectively) between our model and that of Nebrin
et al. (2023). The three sets of curves compare our calculation (eq. 45; solid
curves), the Nebrin et al. (2023) model (dashed curves), and our model with
the Nebrin et al. (2023) density criterion (dotted curves). This is analogous
to their Figure 9. For ease of comparison, all take the Lyman-Werner radiation
background presented in Incatasciato et al. (2023).

While this project was being completed, a related study developing a semi-
analytic estimate for the minimum mass (akin to our § 2) was released (Nebrin
et al., 2023). In Figure 12, we compare our results to their minimum mass
model, adopting the same Lyman-Werner background evolution (Incatasciato et
al., 2023) and cooling threshold (i.e., $t_{\rm cool}<6t_{\rm ff}=\zeta
t_{H}\implies\zeta\simeq 0.16$) to simplify the comparison. While their model
is fairly similar in spirit to that presented in this work, in detail, there
are some key differences that result in distinct calculations of the minimum
mass. For example, they consider the joint effects of $\mathrm{H}_{2}$
cooling, LW feedback, and the stream velocity on the minimum mass, as we do.
They include the effects of reionization feedback—which we do not
include—though it is likely to be significant only at fairly late times
($z\lesssim 10$) and be patchy until reionization is complete (at $z\lesssim
6$). We do not include such inhomogeneous mechanisms here. They do not include
the effects of an X-ray background, which we show to be important in both
promoting and suppressing high redshift star formation between $z\sim 6-20$ (§
2.5). Therefore, we limit our comparison with their results to the overlap:
the calculations of the cooling, photodissociation, and thermalization
thresholds, for which there are some differences.

Based on the comparison in Figure 12, the similarities between our models are
clear. That is, the redshift evolution is comparable and the minimum mass
clearly transitions between the cooling mass and LW mass at roughly the same
time in the two models. Despite the broad similarity, there is a difference of
a factor of a few in the magnitude of the minimum mass at the highest
redshifts (most notably in the case with no streaming). This difference can be
attributed to their choice of central halo gas density (i.e., comparing our
equation 12 to their equation 7). As a result, our model has a critical
$\mathrm{H}_{2}$ cooling fraction that is two orders of magnitude smaller than
their equation 16.131313Evaluating our eq. 22 for $M=10^{6}M_{\odot}$ and
converting from $\zeta\to\eta$ at $z=9$, we find $f_{\mathrm{H}_{2},\rm
crit}\simeq 2.08\times 10^{-5}$. This also results in a somewhat different
redshift evolution of the cooling mass (we find a different temperature
scaling for the low-mass density threshold). This difference will also feed
into the calculation of the dissociation threshold. Indeed, if we modify our
density threshold to match theirs (dotted curves in Figure 12), much of the
difference between our models at small streaming velocities is erased.

Our treatment of the contribution to the minimum mass due to the stream
velocity is identical to that of Nebrin et al. (2023) (i.e., we both employ
the threshold identified in Fialkov et al., 2012), but we also include the
calculation of the filtering scale (accretion threshold) as a baseline for our
model and choose $\alpha_{\rm vbc}=5$ for our fiducial calculations (while
they set $\alpha_{\rm vbc}=6$). As such, we see similar evolution of the
minimum mass with increasing $v_{\rm bc}$, though the magnitude of this
increase differs for the same reasons as outlined above. Despite these
differences, at late times our calculations converge to the same minimum mass
value $M_{\rm min}\to 10^{8}M_{\odot}$ as $z\to 6$.

In Figure 8, we overlay the results of our semi-analytic calculation assuming
the same fiducial parameters as Nebrin et al. (2023) (dashed curves; i.e.,
with their density threshold, taking $\zeta=0.16$, and setting $\alpha_{\rm
vbc}=6$) to test how these assumptions feed into the predictions made with our
semi-analytic model. At high redshifts, the evolution of the minimum mass is
qualitatively similar between the two cases, with the increased $\alpha_{\rm
vbc}$ setting the discrepancy. As expected, this leads to a slightly earlier
onset of star formation in our models than we would see with a stronger effect
from the stream velocity. At late times, the differences between the models
are a result of the different density parameterizations. In particular, our
high-mass density threshold is normalized to a somewhat smaller value than
that chosen in Nebrin et al. (2023) (cf. we choose $n_{\rm H}\sim 12\ {\rm
cm}^{-3}$ while they choose $n_{\rm H}\sim 32\ {\rm cm}^{-3}$ at $z=20$). With
all else fixed, this means that, once we are in the high-mass regime, a halo
needs to be more massive to approach the same critical dissociation threshold.
Because the density evolves with redshift, this is equivalent to an earlier
rise in the minimum mass to reach the same fixed density, which we see.
Combining the two effects, the resulting SFRDs are effectively discrepant by a
constant offset of $\Delta z\sim 1$. While this would in turn produce
different observable signatures, we emphasize that such discrepancies are
significantly smaller than the uncertainties in other parameters outside of
the minimum mass model, such as the selection of $f_{X}$ and the Pop III IMF.

###  5.2 The Pop III star formation history

Next, we focus our comparison on two recent works studying semi-analytic
modeling of the formation of the first stars: Visbal et al. (2020) and Magg et
al. (2022). We note that both of these works use the same LW feedback
prescription used in M18, which does not include the effects of self-shielding
and to which we compare our calculation in Figure 6. We demonstrate that
inclusion of $\mathrm{H}_{2}$ self-shielding can suppress the minimum mass by
nearly an order of magnitude during this era and so expect that these results
will overpredict the suppression of Pop III star formation by the LW
background (and hence underpredict the resulting SFRD). We also note that
neither of these works includes X-ray feedback.

Visbal et al. (2020) apply analytic models for star formation to DM halo
merger trees from cosmological simulations to model the Pop III star formation
process. This approach, combined with their grid-based calculation of the
radiation backgrounds, allows them to account for local, inhomogeneous effects
on the SFRD, such as clustering and mergers. Rather than following star
formation in individual halos (as we do), they assume a fixed fraction (0.001)
of the gas in a halo is converted into Pop III stars in a single star
formation event, corresponding to roughly $100M_{\odot}$ of stars in a
$10^{6}M_{\odot}$ minihalo. After a fixed delay time, these halos are assumed
to have re-accreted metal-enriched gas and will transition to equilibrium Pop
II star formation with a similar fixed star formation efficiency of 0.05.
While qualitatively similar to our approach, their model differs in its
details: we have a longer delay for re-accretion of enriched gas ($\sim 50$
Myr, c.f. their Fig. 3), and we treat the subsequent enriched star formation
more carefully (without imposing a fixed star formation efficiency). But the
most important difference is in our use of a self-shielding prescription.
Ultimately, they find sustained Pop III star formation to $z=6$ (as we do),
but find a Pop III SFRD nearly an order of magnitude smaller than our results,
likely primarily because of differences in the minimum mass. In their model,
the transition to Pop II star formation also occurs at earlier times than ours
($z_{\rm trans}\sim 25-30$ for our model). Nevertheless, the overall histories
are qualitatively similar.

Magg et al. (2022) similarly develop their semi-analytic model using DM merger
trees from cosmological N-body simulations with the goal of analyzing the
effects of variations in the recovery time (time to transition to metal
enriched star formation) on the 21-cm signal (see § 6.2). In a procedure akin
to that of Visbal et al. (2020), they assume a fixed star formation efficiency
for Pop II stars and for their total mass of Pop III stars. However, they do
not explicitly include the LW feedback or stream velocity in their calculation
of the minimum mass, and they ignore X-ray feedback. Instead, they
parameterize the minimum mass in terms of a critical temperature and calculate
this threshold (which is based on the local circular velocity) using
individual pixels in their 21-cm simulation box. They also employ a different
IMF, choosing a power-law with $\phi\propto M^{0.5}$ (i.e., a very top-heavy
IMF), compared to our Chabrier IMF, which has a power-law slope of $-2.35$ at
large masses. Similarly to Visbal et al. (2020), they specify a single event
of Pop III star formation, followed by a fixed recovery time before a
transition to metal-enriched star formation (for which they consider several
possibilities; our fiducial model corresponds most closely to their
‘intermediate’ (30 Myr) transition case.

In this case, we find qualitative agreement with their model, although their
transition redshift is somewhat later than ours ($z_{\rm transition}\sim
20-25$ compared to $25-30$ in our fiducial model) and they find much higher
peak levels of Pop III star formation at a level about an order of magnitude
larger than the peak late-time value we see. These are likely a result of
their much higher star formation efficiencies ($\approx 0.02$, which is a
factor of 20 larger than the value chosen in Visbal et al. 2020,and
corresponds to a $M_{\rm III}\sim 10^{3}M_{\odot}$ starburst in the smallest
halos). Because many more Pop III stars form in their model than ours, it
takes longer for the Pop II star formation to catch up and they find a
somewhat later transition redshift. Despite the differences, Pop III star
formation persists to late times in both of our calculations, though they find
nearly an order of magnitude larger Pop III contribution to the SFRD than in
our model. This illustrates the sensitivity of the SFRD to the assumed Pop III
star formation efficiency. Otherwise, the qualitative similarity of our
histories suggest that it is the total Pop III star formation efficiency,
integrated over all star formation episodes, that sets the amplitude of the
SFRD, whether or not a single event is assumed or repeated cycles (as found in
the simulations of Abe et al. 2021).

##  6 Observational Implications

We extend our semi-analytic model to generate preliminary predictions of
observable signatures of Pop III star formation in minihalos in the early
universe. Unfortunately, direct observations of Pop III halos will be very
challenging with both current and forthcoming telescopes. Indeed, M18
calculate that these halos will have absolute magnitudes between $M_{\rm
AB}\approx-10$ and $-5$. Indirect signatures of Pop III stars—i.e., in their
transients or in their effect on their surroundings—offer a far more promising
avenue to constrain the physics of Pop III star formation with current
technology.

Figure 13: Redshift evolution of the Pop III SN rate for several models. We
show results for our fiducial model (solid curves), the fiducial model in a
region with no streaming (dotted curves), and a model with a much higher
characteristic Pop III mass (dashed curves). In each case, we show the PISN
and core-collapse rates separately (black and red curves, respectively).

###  6.1 Pop III supernovae

Stars with masses between $140-260M_{\odot}$ are expected to end their lives
in pair-instability supernovae (PISNe; Barkat et al., 1967, Fryer et al.,
2001, Heger & Woosley, 2002, Heger & Woosley, 2010). These superenergetic
supernovae can produce nearly two orders of magnitude more energy than
traditional core-collapse SNe and are so powerful that they completely tear
the progenitor star apart, leaving no compact remnant behind. If Pop III stars
have a top-heavy IMF—as we expect from numerical simulations—it is likely that
some of them will form in this mass range and could produce PISNe. We can
calculate the transient rate as:

$\frac{d^{2}N}{dt_{\rm obs}d\Omega_{\rm obs}}(z)=\frac{\eta_{\rm
IMF}}{1+z}\frac{d^{2}V}{dzd\Omega_{\rm obs}}\rho_{\rm SFR}(z)$ (47)

where $d^{2}V/dzd\Omega$ is the differential comoving volume element,
$\eta_{\rm IMF}$ is the number of progenitors per unit stellar mass (and is
determined by the IMF), the factor of $1/(1+z)$ accounts for cosmological time
dilation, and $\rho_{\rm SFR}(z)$ is the usual SFRD calculated from our semi-
analytic model. From this, we can identify the rates associated with different
types of SNe by integrating the IMF over our range of interest

$\eta_{\rm IMF}=\frac{\int_{X}m\phi(m)dm}{\int m\phi(m)dm}$ (48)

where the numerator is an integral over the range of progenitor masses for the
transient of interest (Lazar & Bromm, 2022).

We summarize the results of this basic calculation in Figure 13. By comparing
runs with two different characteristic masses in the IMF—moderately top-heavy
(solid) and extremely top-heavy (dashed)—we demonstrate the sensitivity of the
supernova rates to the details of the IMF. That is, if the Pop III IMF is
truly very top-heavy, as is shown in the dashed curve, the PISN rate will be
several orders of magnitude larger than the core-collapse SN rate.

Of course, these supernovae are only visible in very deep surveys. Kasen et
al. (2011) simulated the light curves of PISNe and found a wide range of
expected luminosities (see also Hummel et al. 2012). Their results suggest
that surveys reaching limiting magnitudes of $m\ga 30$ will be required to
identify these sources, which (depending on the PISN model) could reach events
at $z\sim 10$–25. Our model demonstrates that, even in optimistic models, such
surveys must span many square degrees to accumulate a significant number of
events.

Figure 14: The differential 21-cm brightness temperature (relative to the CMB)
for three representative Pop III star formation models—no stream velocity and
weak X-rays (blue), no streaming and strong X-rays (red), and average
streaming with moderate X-rays (the fiducial model; black)—and one with no Pop
II star formation (dashed black). For each of these models, the transition
from Pop III to Pop II dominance in the SFRD is indicated with a colored
vertical dotted line.

###  6.2 21-cm global signal

Another potential indirect signature of Pop III stars is in their impact on
their environment, the IGM. In particular, the sky-averaged 21-cm ‘spin-flip’
signal of atomic hydrogen provides a tracer of the IGM that is sensitive to
the electron fraction, temperature, and Ly$\alpha$ intensity. We have already
demonstrated that high redshift star formation has the potential to very
efficiently heat and ionize the IGM (Figure 4). Mirocha et al. (2018) and
Mebane et al. (2020) demonstrated that the inclusion of Pop III stars in
calculations in the 21-cm signal can introduce a characteristic asymmetry in
the absorption trough due to early star formation.

Here we extend our semi-analytic model to include calculation of the
Ly$\alpha$ background following Holzbauer & Furlanetto (2012):

$J_{\alpha}(z)=\sum_{n=2}^{n_{\rm max}}\int_{z}^{z_{\rm
max}(n)}dz^{\prime}f_{\rm
recycle}(n)\frac{(1+z)^{2}}{4\pi}\frac{c}{H(z^{\prime})}\epsilon(\nu_{n}^{\prime},z^{\prime})$
(49)

where $z_{\rm max}$ is the maximum redshift that a Ly$n$ photon can experience
before it redshifts into the next transition, and $\nu_{n}^{\prime}$ is
defined as above (equation 42). We sum over Ly$n$ levels because a fraction
$f_{\rm recycle}(n)$ of Ly$n$ photons will cascade down to become Ly$\alpha$
photons via a series of radiative transitions (Hirata, 2006, Pritchard &
Furlanetto, 2006). In practice, we truncate this sum at $n=23$ to exclude the
levels for which the photon horizon lies within the size of a ${\rm H\ II}$
region of a typical galaxy at this time. We calculate the emissivity with a
procedure akin to our calculation of the LW emissivity:

$\epsilon(z)=\int_{M_{\rm
min}}^{\infty}n(M)\frac{\Omega_{b}}{\Omega_{m}}\frac{\dot{M}_{*}}{m_{p}}\bigg{(}\frac{N_{{\rm
Ly}n}E_{{\rm Ly}n}}{\Delta\nu_{{\rm Ly}n}}\bigg{)}dM$ (50)

where $M_{\rm min}$ is the minimum mass for star formation, $\dot{M}_{*}$ is
the star formation rate of a halo of mass $M$, $N_{{\rm Ly}n}$ is the number
of Ly$n$ photons produced per baryon in stars (Barkana & Loeb, 2005), $E_{{\rm
Ly}n}$ is the average energy of a Ly$n$ photon, and $\Delta\nu_{{\rm Ly}n}$ is
the frequency spacing between the $n$ and $(n+1)$st Lyman line.

The spin temperature of the 21-cm line is set by a competition between
scattering of CMB photons, collisions, and scattering of Ly$\alpha$ photons
(Wouthuysen, 1952, Field, 1958, Furlanetto et al., 2006). Following Furlanetto
(2006) and Mebane et al. (2020) (and the references therein), we calculate the
coupling coefficients $x$ to compute the spin temperature

$1-\frac{T_{\gamma}}{T_{S}}=\frac{x_{c}+x_{\alpha}}{1+x_{c}+x_{\alpha}}\bigg{(}1-\frac{T_{\gamma}}{T_{K}}\bigg{)}$
(51)

and from this the 21-cm brightness temperature (relative to the backlight of
the CMB)

$\delta T_{b}=27x_{\rm
HI}\bigg{(}\frac{\Omega_{b}h^{2}}{0.023}\bigg{)}\bigg{(}\frac{0.15}{\Omega_{m}h^{2}}\frac{1+z}{10}\bigg{)}^{1/2}\bigg{(}\frac{T_{S}-T_{\gamma}}{T_{S}}\bigg{)}\
{\rm mK}$ (52)

where $x_{\rm HI}$ is the neutral fraction.

With this, we compute the differential 21-cm brightness temperature for three
representative models (Figure 14). Note that we have not included the ionizing
UV background that will drive the process of cosmic reionization, so these
calculations are meant to be diagnostics of the early universe signatures of
Pop III physics on the 21-cm global signal. We have also ignored the potential
effects of a cosmic radio background from accreting black holes.

Indeed, the shape and timing of the signal are clearly sensitive to the global
processes that govern the SFR. When star formation is strongly coupled to the
X-ray luminosity (red curve; strong X-rays), the IGM is efficiently heated and
the depth of the absorption trough is suppressed. As a result, the signal
turns over to emission more quickly, yielding a very sensitive probe of the
X-ray background strength. In kind, there is a clear signature of the stream
velocity in the width and timing of the absorption trough. That is, with the
stream velocity, star formation is delayed and the absorption trough in turn
is delayed as well. At late times, however, the star formation rate quickly
climbs to match the non-streaming levels and the signals are
indistinguishable.

For $z\gtrsim 20-25$, the contributions of Pop III and II stars to the SFRD
(and thus global signal) are comparable. Indeed, comparing our fiducial model
to one with no Pop II contribution to the radiation backgrounds (dashed
black), the depth and shape of the absorption trough are different at late
times (following the transition to Pop II star formation). At early times, Pop
III stars tightly couple the spin temperature to the gas kinetic temperature
driven by the Wouthuysen-Field effect. At later times, the spin temperature
remains tightly coupled to the temperature of the IGM, but the enhanced SFRD
(in response to the growth of Pop II halos) drives more rapid heating of the
IGM and thus the absorption signal is somewhat suppressed.

Comparing these calculations to the results of Mebane et al. (2020), we find
higher-redshift features than seen before. Namely, the updated calculation of
the LW mass introduces an earlier onset of the absorption trough than seen in
Mebane et al. (2020). We find that variations in the X-ray background strength
can suppress the depth of the absorption trough, as they do.

We can also compare our 21-cm global signal predictions to those made in Magg
et al. (2022). We note that they neglect to include the effects of X-rays
(which will dominate the heating of the IGM) in their calculation of the 21-cm
signal. As expected, we find a similar absorption depth and qualitative global
signal to their ‘intermediate’ transition model. However, we find that the
nadir of the absorption trough is achieved at higher redshifts in our model
($z_{\rm abs}\sim 20-25$ compared to 10-15 in their model) and find a
characteristic high redshift contribution to the absorption from Pop III stars
that is not reflected in their results. This discrepancy is puzzling because
we see comparable (or larger) contributions to the Pop III SFRD between their
model and ours (see § 5.2). We speculate that this is perhaps a result of the
lack of X-ray heating included in their model. For example, if we compare the
weak and strong X-ray curves (blue and red in Figure 14, respectively), we see
that the inclusion of a strong X-ray background moves the nadir of the
absorption trough earlier and suppresses its depth. In Magg et al. (2022), the
only source of heating is from Ly$\alpha$ photons, whereas we include both
contributions, so heating in their model proceeds more slowly.

Though this calculation makes several simplifying assumptions, this estimate
demonstrates that variations in the Pop III star formation physics are
reflected in the 21-cm global signal, especially at the highest redshifts.
Indeed, in all of the displayed cases, we find the characteristic asymmetry in
the absorption trough produced by early IGM heating from high-redshift Pop III
star formation identified by Mirocha et al. (2018). Forthcoming low-frequency
radio telescopes, such as the lunar FarView array, plan to study the early
universe through the highly redshifted 21-cm line. Measurements of the global
signal at frequencies below 70 MHz will offer a strong and promising probe of
the details of Pop III star formation physics.

##  7 Conclusions

We have presented a simple analytic model for the minimum DM halo mass for Pop
III star formation that incorporates the combined effects of a relic relative
DM-baryon ‘stream’ velocity from the early universe and feedback from the
buildup of UV/LW and X-ray backgrounds. Such a criterion is the crux of any
semi-analytic model for high-redshift star formation and allows us to self-
consistently model the formation of the first stars and the transition to
subsequent generations of metal-enriched star formation. We incorporate this
analytic calculation, a Chabrier-like IMF, and a criterion for ‘bursty’ Pop II
star formation as updates to the semi-analytic model of M18.

From this model, we identify three key epochs of Pop III star formation. At
the earliest times ($z\gtrsim 30$), the stream velocity is the dominant
environmental factor in setting the Pop III SFRD—so the SFRD is most sensitive
to choice of $v_{\rm bc}$ in this era. However, the stream velocity decays as
$v_{\rm bc}\propto(1+z)$, so for intermediate times ($30\gtrsim z\gtrsim 10$),
the LW background sets the Pop III SFRD, coincident with the transition to Pop
II star formation dominating the SFRD. Finally, at the latest times
($z\lesssim 10$), the persistence of Pop III star formation is dictated by the
strength of the X-ray background. While bursty star formation can suppress Pop
III star formation through the LW background, we find that this is a secondary
effect to those of the aforementioned three processes.

These distinct epochs in the minimum mass are translated into variations in
the Pop III SFRD. In particular, the onset of early star formation is directly
responsive to the strength of the DM-baryon relative velocity and the late-
time magnitude of the SFRD (i.e., the peak of the SFR) is sensitive to the
physics of the LW and X-ray backgrounds. Therefore, observations that even
indirectly probe the Pop III SFRD will be a powerful tool for understanding
the physics governing star formation in that era.

While Pop III halos will likely not be directly observable with current (or
near-future) technology, indirect signatures of these stars are a promising
avenue to probe this epoch. We find that superluminous Pop III SNe are
theoretically observable with forthcoming deep JWST and RST surveys and can
shed light into the IMF for these first stars. Similarly, global 21-cm
experiments that are able to reach the lowest frequencies will offer us a
window into Cosmic Dawn and will potentially be able to constrain the physics
of Pop III star formation. Recent work has shown that Pop III star formation
is not limited to Pop III halos alone; accretion of pristine gas from the IGM
can result in Pop III star formation in Pop II halos, though these stars will
likely be very rare and short-lived (Venditti et al., 2023).

Although our model includes most of the key physical processes regulating Pop
III star formation, we have ignored processes for which small-scale
inhomogeneities are essential (such as metal enrichment and photoheating from
UV photons during reionization). In the future, incorporating these processes
will be essential for modeling spatial fluctuations in the Pop III population.

## Data Availability

No new data were obtained as part of this work. Results used to generate the
figures are available from the authors upon reasonable request.

## Acknowledgements

The authors thank Matt McQuinn, Rick Mebane, Smadar Naoz, Claire Williams,
William Lake, Massimo Ricotti, and Roy Zhao for useful discussions. This work
was supported by the National Science Foundation through award AST-1812458. In
addition, this work was directly supported by the NASA Solar System
Exploration Research Virtual Institute cooperative agreement number
80ARC017M0006.

Software: numpy (van der Walt et al., 2011), astropy (Astropy Collaboration et
al., 2013), matplotlib (Hunter, 2007), scipy (Virtanen et al., 2020)

## References

  * Abe et al. (2021) Abe M., Yajima H., Khochfar S., Dalla Vecchia C., Omukai K., 2021, MNRAS, 508, 3226
  * Abel et al. (2002) Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Barkana & Loeb (2001) Barkana R., Loeb A., 2001, Phys. Rep., 349, 125
  * Barkana & Loeb (2005) Barkana R., Loeb A., 2005, ApJ, 626, 1
  * Barkat et al. (1967) Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379
  * Bromm (2013) Bromm V., 2013, Reports on Progress in Physics, 76, 112901
  * Bromm & Loeb (2003) Bromm V., Loeb A., 2003, Nature, 425, 812
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chluba & Thomas (2013) Chluba J., Thomas R. M., 2013, CosmoRec: Cosmological Recombination code, Astrophysics Source Code Library, record ascl:1304.017 (ascl:1304.017)
  * Dalal et al. (2010) Dalal N., Pen U.-L., Seljak U., 2010, J. Cosmology Astropart. Phys., 2010, 007
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton University Press
  * Draine & Bertoldi (1996) Draine B. T., Bertoldi F., 1996, ApJ, 468, 269
  * Ferrara et al. (2000) Ferrara A., Pettini M., Shchekinov Y., 2000, MNRAS, 319, 539
  * Fialkov (2014) Fialkov A., 2014, International Journal of Modern Physics D, 23, 1430017
  * Fialkov et al. (2012) Fialkov A., Barkana R., Tseliakhovich D., Hirata C. M., 2012, MNRAS, 424, 1335
  * Field (1958) Field G. B., 1958, Proceedings of the IRE, 46, 240
  * Fragos et al. (2013) Fragos T., et al., 2013, ApJ, 764, 41
  * Fryer et al. (2001) Fryer C. L., Woosley S. E., Heger A., 2001, ApJ, 550, 372
  * Furlanetto (2006) Furlanetto S. R., 2006, MNRAS, 371, 867
  * Furlanetto & Loeb (2003) Furlanetto S. R., Loeb A., 2003, ApJ, 588, 18
  * Furlanetto & Mirocha (2022) Furlanetto S. R., Mirocha J., 2022, MNRAS, 511, 3895
  * Furlanetto & Stoever (2010) Furlanetto S. R., Stoever S. J., 2010, MNRAS, 404, 1869
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Furlanetto et al. (2017) Furlanetto S. R., Mirocha J., Mebane R. H., Sun G., 2017, MNRAS, 472, 1576
  * Galli & Palla (1998) Galli D., Palla F., 1998, A&A, 335, 403
  * Gnedin (2000) Gnedin N. Y., 2000, ApJ, 542, 535
  * Gnedin & Hui (1998) Gnedin N. Y., Hui L., 1998, MNRAS, 296, 44
  * Greig et al. (2021) Greig B., et al., 2021, MNRAS, 501, 1
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  * Haiman et al. (1996a) Haiman Z., Thoul A. A., Loeb A., 1996a, ApJ, 464, 523
  * Haiman et al. (1996b) Haiman Z., Rees M. J., Loeb A., 1996b, ApJ, 467, 522
  * Haiman et al. (2000) Haiman Z., Abel T., Rees M. J., 2000, ApJ, 534, 11
  * Heger & Woosley (2002) Heger A., Woosley S. E., 2002, ApJ, 567, 532
  * Heger & Woosley (2010) Heger A., Woosley S. E., 2010, ApJ, 724, 341
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, The Astrophysical Journal, 591, 288
  * Hickox & Markevitch (2007) Hickox R. C., Markevitch M., 2007, ApJ, 661, L117
  * Hirano et al. (2015) Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS, 448, 568
  * Hirata (2006) Hirata C. M., 2006, MNRAS, 367, 259
  * Holzbauer & Furlanetto (2012) Holzbauer L. N., Furlanetto S. R., 2012, MNRAS, 419, 718
  * Hummel et al. (2012) Hummel J. A., Pawlik A. H., Milosavljević M., Bromm V., 2012, ApJ, 755, 72
  * Hummel et al. (2015) Hummel J. A., Stacy A., Jeon M., Oliveri A., Bromm V., 2015, MNRAS, 453, 4136
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Incatasciato et al. (2023) Incatasciato A., Khochfar S., Oñorbe J., 2023, arXiv e-prints, p. arXiv:2301.08242
  * Jaacks et al. (2018) Jaacks J., Thompson R., Finkelstein S. L., Bromm V., 2018, MNRAS, 475, 4396
  * Kasen et al. (2011) Kasen D., Woosley S. E., Heger A., 2011, ApJ, 734, 102
  * Kulkarni et al. (2013) Kulkarni G., Rollinde E., Hennawi J. F., Vangioni E., 2013, ApJ, 772, 93
  * Kulkarni et al. (2021) Kulkarni M., Visbal E., Bryan G. L., 2021, ApJ, 917, 40
  * Lake et al. (2023) Lake W., et al., 2023, ApJ, 943, 132
  * Lazar & Bromm (2022) Lazar A., Bromm V., 2022, MNRAS, 511, 2505
  * Lehmer et al. (2016) Lehmer B. D., et al., 2016, ApJ, 825, 7
  * 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
  * Machacek et al. (2003) Machacek M. E., Bryan G. L., Abel T., 2003, MNRAS, 338, 273
  * Madau et al. (2001) Madau P., Ferrara A., Rees M. J., 2001, ApJ, 555, 92
  * Magg et al. (2022) Magg M., et al., 2022, MNRAS, 514, 4433
  * Maio et al. (2010) Maio U., Ciardi B., Dolag K., Tornatore L., Khochfar S., 2010, MNRAS, 407, 1003
  * McKee & Tan (2008) McKee C. F., Tan J. C., 2008, ApJ, 681, 771
  * McQuinn (2012) McQuinn M., 2012, MNRAS, 426, 1349
  * McQuinn & O’Leary (2012) McQuinn M., O’Leary R. M., 2012, ApJ, 760, 3
  * Mebane et al. (2018) Mebane R. H., Mirocha J., Furlanetto S. R., 2018, MNRAS, 479, 4544
  * Mebane et al. (2020) Mebane R. H., Mirocha J., Furlanetto S. R., 2020, MNRAS, 493, 1217
  * Mesinger et al. (2013) Mesinger A., Ferrara A., Spiegel D. S., 2013, MNRAS, 431, 621
  * Mineo et al. (2012) Mineo S., Gilfanov M., Sunyaev R., 2012, MNRAS, 419, 2095
  * Mirocha et al. (2018) Mirocha J., Mebane R. H., Furlanetto S. R., Singal K., Trinh D., 2018, MNRAS, 478, 5591
  * Muñoz et al. (2022) Muñoz J. B., Qin Y., Mesinger A., Murray S. G., Greig B., Mason C., 2022, MNRAS, 511, 3657
  * Nakatani et al. (2020) Nakatani R., Fialkov A., Yoshida N., 2020, ApJ, 905, 151
  * Naoz & Barkana (2007) Naoz S., Barkana R., 2007, MNRAS, 377, 667
  * Naoz et al. (2012) Naoz S., Yoshida N., Gnedin N. Y., 2012, ApJ, 747, 128
  * Naoz et al. (2013) Naoz S., Yoshida N., Gnedin N. Y., 2013, ApJ, 763, 27
  * Nebrin et al. (2023) Nebrin O., Giri S. K., Mellema G., 2023, arXiv e-prints, p. arXiv:2303.08024
  * O’Leary & McQuinn (2012) O’Leary R. M., McQuinn M., 2012, ApJ, 760, 4
  * O’Shea & Norman (2008) O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14
  * Park et al. (2021) Park J., Ricotti M., Sugimura K., 2021, MNRAS, 508, 6176
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Pritchard & Furlanetto (2006) Pritchard J. R., Furlanetto S. R., 2006, MNRAS, 367, 1057
  * Ricotti (2016) Ricotti M., 2016, MNRAS, 462, 601
  * Ricotti & Ostriker (2004) Ricotti M., Ostriker J. P., 2004, MNRAS, 352, 547
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sarmento et al. (2018) Sarmento R., Scannapieco E., Cohen S., 2018, ApJ, 854, 75
  * Schauer et al. (2021) Schauer A. T. P., Glover S. C. O., Klessen R. S., Clark P., 2021, MNRAS, 507, 1775
  * Seager et al. (1999) Seager S., Sasselov D. D., Scott D., 1999, ApJ, 523, L1
  * Shang et al. (2010) Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249
  * Shull & van Steenberg (1985) Shull J. M., van Steenberg M. E., 1985, ApJ, 298, 268
  * Stacy et al. (2012) Stacy A., Greif T. H., Bromm V., 2012, MNRAS, 422, 290
  * Tegmark et al. (1997) Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F., 1997, ApJ, 474, 1
  * The HERA Collaboration et al. (2022) The HERA Collaboration et al., 2022, arXiv e-prints, p. arXiv:2210.04912
  * Trac et al. (2015) Trac H., Cen R., Mansfield P., 2015, ApJ, 813, 54
  * Trenti & Stiavelli (2009) Trenti M., Stiavelli M., 2009, ApJ, 694, 879
  * Trenti et al. (2009) Trenti M., Stiavelli M., Shull J. M., 2009, ApJ, 700, 1672
  * Tseliakhovich & Hirata (2010) Tseliakhovich D., Hirata C., 2010, Phys. Rev. D, 82, 083520
  * Venditti et al. (2023) Venditti A., Graziani L., Schneider R., Pentericci L., Di Cesare C., Maio U., Omukai K., 2023, arXiv e-prints, p. arXiv:2301.10259
  * Verner et al. (1996) Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Visbal et al. (2014a) Visbal E., Haiman Z., Bryan G. L., 2014a, MNRAS, 442, L100
  * Visbal et al. (2014b) Visbal E., Haiman Z., Terrazas B., Bryan G. L., Barkana R., 2014b, MNRAS, 445, 107
  * Visbal et al. (2018) Visbal E., Haiman Z., Bryan G. L., 2018, MNRAS, 475, 5246
  * Visbal et al. (2020) Visbal E., Bryan G. L., Haiman Z., 2020, ApJ, 897, 95
  * Williams et al. (2022) Williams C. E., et al., 2022, arXiv e-prints, p. arXiv:2211.02066
  * Wise & Abel (2007) Wise J. H., Abel T., 2007, ApJ, 671, 1559
  * Wolcott-Green & Haiman (2019) Wolcott-Green J., Haiman Z., 2019, MNRAS, 484, 2467
  * Wolcott-Green et al. (2011) Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Xu et al. (2016) Xu H., Norman M. L., O’Shea B. W., Wise J. H., 2016, ApJ, 823, 140
  * Yamaguchi et al. (2023) Yamaguchi N., Furlanetto S. R., Trapp A. C., 2023, MNRAS, 520, 2922
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22