# A halo model approach for the 21-cm power spectrum at cosmic dawn

Aurel Schneider<EMAIL_ADDRESS>Center for Cosmology and Theoretical
Astrophysics, Institute for Computational Science, University of Zurich,
Switzerland.  Sambit K. Giri<EMAIL_ADDRESS>Center for Cosmology and
Theoretical Astrophysics, Institute for Computational Science, University of
Zurich, Switzerland.  Jordan Mirocha<EMAIL_ADDRESS>McGill
University, Department of Physics & McGill Space Institute, Montréal, Canada.

###### Abstract

Prior to the epoch of reionisation, the 21-cm signal of the cosmic dawn is
dominated by the Lyman-$\alpha$ coupling and gas temperature fluctuations
caused by the first sources of radiation. While early efforts to model this
epoch relied on analytical techniques, the community quickly transitioned to
more expensive semi-numerical models. Here, we re-assess the viability of
simpler approaches that allow for rapid explorations of the vast astrophysical
parameter space. We propose a new analytical method to calculate the 21-cm
power spectrum based on the framework of the halo model. Both the
Lyman-$\alpha$ coupling and temperature fluctuations are described by
overlapping radiation flux profiles that include spectral red-shifting and
source attenuation due to look-back (light-cone) effects. The 21-cm halo model
is compared to the semi-numerical code 21cmFAST exhibiting generally good
agreement, i.e., the power spectra differ by less than a factor of three over
a large range of $k$-modes and redshifts. We show that the remaining
differences between the two methods are comparable to the expected variations
from modelling uncertainties associated with the abundance, bias, and
accretion rates of haloes. While these current uncertainties must be reduced
in the future, our work suggests that inference at acceptable accuracy will
become feasible with very efficient halo models of the cosmic dawn.

##  I Introduction

Current and upcoming radio interferometer telescopes such as LOFAR (Haarlem
_et al._ , 2006; Mertens  _et al._ , 2020), MWA (Tingay  _et al._ , 2013;
Trott  _et al._ , 2020), HERA (DeBoer  _et al._ , 2017), or SKA (Koopmans  _et
al._ , 2015) are expected to detect for the first time the strongly redshifted
21-cm clustering signal during and prior to the epoch of reionisation. These
measurements will open up a new window on the universe (Madau  _et al._ ,
1997; Shaver  _et al._ , 1999), providing insights into both astrophysics
(Ewall-Wice  _et al._ , 2018; Mirocha  _et al._ , 2019; Park  _et al._ , 2019)
and cosmology (McQuinn  _et al._ , 2006; Liu  _et al._ , 2019). Next to the
very first stars and galaxies (Barkana and Loeb, 2005a; Greig  _et al._ ,
2020), the 21-cm signal has the potential to find evidence for exotic sources
of radiation (Liu and Slatyer, 2018; Mena  _et al._ , 2019), new signs from
the dark matter sector (Sitwell  _et al._ , 2014; Lopez-Honorez  _et al._ ,
2016; Schneider, 2018; Lidz and Hui, 2018; Nebrin  _et al._ , 2019) or other
deviations from the standard $\Lambda$CDM cosmological model (Hill and Baxter,
2018; Yang  _et al._ , 2019; Lopez-Honorez  _et al._ , 2020).

Predicting the 21-cm signal of the cosmic dawn, however, remains a challenging
task. On the one hand, the complicated physics of radiation-hydrodynamics
combined with the enormous ranges of relevant scales make brute-force
simulations extremely difficult. On the other hand, the poorly known
characteristics of early star-forming sources as well as the complicated
interplay between gas cooling and feedback add important uncertainties that
need to be either understood or properly parametrised.

First attempts to predict the 21-cm clustering before reionisation were based
on analytical techniques, using a combination of cosmological perturbation
theory and excursion-set modelling prescriptions for the sources (Furlanetto
_et al._ , 2004; Barkana and Loeb, 2005a; Pritchard and Furlanetto, 2007) (for
a more recent attempt see also Ref. (Raste and Sethi, 2018)). While these
calculations were able to predict many important features of the 21-cm power
spectrum – such as for example the characteristic double-peaked shape due to
the Lyman-$\alpha$ coupling and heating epochs – it remains unclear how well
they agree with more detailed calculations (Santos  _et al._ , 2008).

A further important step towards more realistic predictions of the 21-cm
signal at cosmic dawn was taken with the development of semi-numerical methods
such as 21cmFAST (Mesinger and Furlanetto, 2007; Mesinger  _et al._ , 2010) or
simfast21 (Santos  _et al._ , 2009) (see also Refs. (Fialkov  _et al._ , 2014;
Mutch  _et al._ , 2016; Hutter, 2018) for other semi-numerical methods). These
models numerically evolve the matter perturbations, Lyman-$\alpha$ coupling,
and temperature fluctuations on a grid, where the source distributions are
either obtained via a halo catalogue from N-body simulations or via an
excursion-set recipe. Semi-numerical techniques consist of a major improvement
with respect to analytical approaches, mainly because they follow the
evolution of the spin temperature in configuration space and are therefore
able to produce maps of the 21-cm signal.

In principle, more accurate predictions can be obtained by post-processing
numerical $N$-body simulations using radiative-transfer (RT) calculations.
This is the strategy followed e.g. by Grizzly (Thomas  _et al._ , 2009; Ghara
_et al._ , 2015), C2-Ray (Mellema  _et al._ , 2006; Friedrich  _et al._ ,
2012), and CRASH (Graziani  _et al._ , 2013). While the former code is based
on approximate but faster one-dimensional RT calculations, the latter two
follow a full three dimensional ray-tracing approach. First results for the
epoch of pre-reionisation can be found in Refs. (Ghara  _et al._ , 2015, 2018)
as well as Ref. (Ross  _et al._ , 2019) and (Eide  _et al._ , 2020).

Full radiation-hydrodynamical simulations, including self-consistent formation
of sources, are becoming an increasingly common tool to study the process of
reionisation (e.g. Refs. Ocvirk  _et al._ , 2016; Doussot  _et al._ , 2019;
Ocvirk  _et al._ , 2020; Rosdahl  _et al._ , 2018; Katz  _et al._ , 2020).
However, for the pre-reionisation epoch of cosmic dawn, such simulations
remain rare (Semelin  _et al._ , 2017). This is because they require narrowly
binned multi-frequency RT calculations in a simulation box that resolves the
mini-haloes hosting the first sources and that accounts for the large
distances travelled by Lyman-$\alpha$ and X-ray radiation.

In this paper we propose an analytical approach to predict the 21-cm global
signal and power spectrum at cosmic dawn. While being potentially less
accurate than the different numerical techniques mentioned above, such a model
has the advantage of providing fast predictions based on a well defined
framework. Different parametrisations of the source modelling or effects from
varying cosmology can be readily implemented and tested. An analytical method
is also particularly well suited for parameter inference, which is an
important aspect of 21-cm cosmology due to the large uncertainties related to
early-universe galaxy formation that have to be parametrised and marginalised
over. Comparing multi-parameter models with observations requires a very large
number of calculations that can only be performed with fast prediction
routines.

Our method is inspired by earlier work of Holzbauer and Furlanetto (2011) who
used the halo model to predict the clustering from the Lyman-$\alpha$ coupling
between the gas and the first star-light at cosmic dawn. We extend this
approach adding a description for the temperature fluctuations as well as an
improved modelling of the halo accretion and the star formation rates. With
this at hand, we obtain a complete prediction of the 21-cm power spectrum at
cosmic dawn.

The paper is structured as follows. In Sec. II we discuss the source modelling
including halo mass function, bias, gas accretion, and star-formation
efficiency. Sec. III summarises the derivation of the differential brightness
temperature. In Sec. IV we present the halo model of flux profiles and we show
results for the 21-cm power spectrum assuming different star-formation
efficiencies. We then go on and compare our model to other analytical and
semi-numerical methods in Sec. V before concluding in Sec. VI.

##  II Source modelling

Accurately quantifying the abundance, distribution, and emission of sources is
a crucial step for any method predicting the 21-cm signal. The source
distribution, for example, can be readily quantified using prescriptions for
the halo mass function and the halo bias. Regarding the source emission,
analytical and semi-numerical models usually rely on estimates for the halo
accretion rates combined with a parametrisation of the star-formation
efficiency. While the former describes the accretion of primordial gas onto
the parent halo, the latter generically parametrises all star formation and
feedback processes.

In the present section, we will first discuss the halo mass function and
biasing model, before we describe and compare different estimates for halo
mass accretion rates. At the end we present our parametrisation for the star
formation efficiency parameter.

###  II.1 Halo mass function and bias

Both the halo mass function and the halo bias are central ingredients of the
halo model which our method is based on. For the halo mass function, we assume
the standard form motivated by the extended Press-Schechter (PS) approach
(Press and Schechter, 1974; Bond  _et al._ , 1991), i.e.,

$\displaystyle\frac{dn}{d\ln
M}=\frac{1}{2}\frac{\bar{\rho}}{M}f(\nu)\frac{d\ln\nu}{d\ln M}$ (1)

where the peak-height of perturbations is defined as
$\nu=\delta_{c}(z)^{2}/\sigma^{2}$ with $\delta_{c}(z)=1.686/D(z)$, $D(z)$
being the cosmological growth factor. The variance of the density field
($\sigma^{2}$) is defined as

$\sigma^{2}(M)=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}P(k)W^{2}(k|M),$ (2)

where $P(k)$ is the linear power spectrum at redshift zero and $W(k|M)$ is the
Fourier transform of the top-hat window function. The first crossing
distribution ($f$) us given by

$\displaystyle
f(\nu)=A\sqrt{2\frac{q\nu}{\pi}}\left[1+(q\nu)^{-p}\right]\exp(-q\nu/2),$ (3)

which comes with two free model parameters $q$ and $p$ (see Ref. Sheth and
Tormen, 2002). The factor A is fixed via the relation
$A=(1+2^{-p}\Gamma(1/2-p)/\sqrt{\pi})^{-1}$.

Using the peak-background split model introduced in Ref. (Sheth and Tormen,
1999), the halo bias can be described within the same extended PS formalism,
leading to the relation (see e.g Ref. Cooray and Sheth, 2002)

$\displaystyle
b(M)=1+\frac{q\nu-1}{\delta_{c}(z)}+\frac{2p}{\delta_{c}(z)[1+(q\nu)^{p}]}.$
(4)

In the following, we use the values $q=0.85$ and $p=0.3$ unless stated
otherwise. We have verified that this setup provides a better fit to the halo
mass function of the high-redshift $N$-body simulations from Refs. (Iliev  _et
al._ , 2011; Schneider, 2018) than the original Press-Schechter (Press and
Schechter, 1974) or Sheth-Tormen (Sheth and Tormen, 2002) mass functions
(which are characterised by $q=1$, $p=0$ and $q=0.707$, $p=0.3$,
respectively).

Figure 1: Halo growth, mass accretion rate, and star-formation rate density
(SFRD) for the three different models discussed in Sec. II.2. _Left:_ Halo
mass as a function of redshift for four different final masses $M$ (green,
brown, red, and blue lines). The orange lines correspond to the mean halo
growth from the B20 (Behroozi  _et al._ , 2020) simulations. _Centre:_ Mass
accretion rates for the same models and halo masses. _Right:_ Resulting star-
formation rate densities (SFRD) for the three different mass-accretion rate
models (see Eq. 13)

###  II.2 Halo mass accretion

There exists different ways to analytically quantify the growth of haloes.
Here we discuss three models that are used in the literature, and we compare
them to results from $N$-body simulations.

####  II.2.1 Exponential (EXP) halo growth

The simplest way to obtain a halo accretion history is to assume that haloes
grow exponentially. This assumption has been shown to provide a good match to
simulations, especially at high redshifts (see e.g. Refs McBride  _et al._ ,
2009; Dekel  _et al._ , 2013; Trac  _et al._ , 2015). We consider the
following equation

$\displaystyle M_{\rm ac}(M,z)=M\exp\left[\alpha(z_{0}-z)\right]$ (5)

where $M$ is the halo mass at the final redshift $z_{0}$. The factor $\alpha$
can be determined using simulations and has been shown to be only weakly mass
dependent at high redshift. Following Ref. (Dekel  _et al._ , 2013) we
therefore assume $\alpha=0.79$ independent of the halo mass $M$.

In the left-hand and centre panels of Fig. 1 we show the redshift evolution of
the halo mass ($M_{\rm ac}$) and the accretion rate ($\dot{M}_{\rm ac}$). The
dash-dotted lines correspond to the exponential (EXP) growth model of Eq. (5).
The model is in very good agreement with the numerical $N$-body simulations
from Behroozi  _et al._ (2020, B20) (orange lines). Note, however, that the
simulations do not cover all redshifts and mass ranges of interest to the
present study.

####  II.2.2 Extended Press-Schechter (EPS) method

Halo accretion rates can also be obtained by means of the extended Press-
Schechter (EPS) formalism. For example, Neistein and van den Bosch (2006)
calculated an ensemble of EPS merger trees and determined the average growth
of their main branches. Based on this, they proposed the equation

$\displaystyle\frac{dM_{\rm ac}}{dz}=-\sqrt{\frac{2}{\pi}}\frac{M_{\rm
ac}}{\sqrt{S(QM_{\rm ac})-S(M_{\rm ac})}}\frac{d\delta_{c}(z)}{dz},$ (6)

where $S\equiv\sigma^{2}$ at redshift zero (see also Ref. Correa  _et al._ ,
2015). This differential equation can be solved assuming $M_{\rm
ac}(M,z_{0})\equiv M$, where $z_{0}$ designates the final redshift of
interest. The value of $Q$ has to be selected empirically, and we use $Q=0.6$.
This number is larger than the range $Q=0.43-0.5$ proposed by Neistein and van
den Bosch (2006) but it provides a better agreement with the simulations from
B20 (see orange lines in Fig. 1).

The predicted EPS halo growth and accretion rates are shown as solid lines in
the left-hand and central panels of Fig. 1. The evolution is similar to the
exponential model except at small masses, where the halo growth is less
pronounced. The model is in good agreement with the B20 simulations, which
should not come as a surprise since we have recalibrated the $Q$-parameter
according to these simulations.

The advantage of the EPS model compared to the much simpler exponential model
lies in the fact that Eq. (6) is, in principle, sensitive to changes in
cosmology. Whether the true cosmology dependence of the mass accretion rate
can be accurately reproduced by the EPS model, remains, however, to be tested.

####  II.2.3 Abundance matching (AM) method

The final method we are investigating here is inspired by the abundance
matching (AM) technique (Vale and Ostriker, 2004) and was first applied in
Ref. (Furlanetto  _et al._ , 2017) as a measure of high-redshift galaxy
growth. The method aims to connect haloes between different redshift bins
($z_{n}$) by matching

$\displaystyle\int_{M_{n}}^{\infty}dM\frac{dn}{dM}(M,z_{n})=\int_{M_{n-1}}^{\infty}dM\frac{dn}{d\ln
M}(M,z_{n-1}),$ (7)

where the halo mass function $dn/dM$ is obtained via Eq. (1). Connecting all
masses $M_{n}$ at different redshifts $z_{n}$ allows us to estimate the
accretion rate of haloes. Note that the AM method connects haloes in a
strictly hierarchical way, which means that the $i$-largest halo at the final
redshift $z_{0}$ will be assumed to stay the $i$ largest halo at all higher
redshifts. As a consequence, the model implicitly assumes haloes to
exclusively grow via smooth accretion, since halo mergers cannot be accounted
for.

In Fig. 1 the halo growth and accretion rates of the AM model are shown as
dashed lines. Compared to the other methods, the AM model predicts
significantly slower halo growth over all redshifts and mass ranges. The AM
results do not match the B20 simulations very well. Note that our findings are
in qualitative agreement with Ref. (Mebane  _et al._ , 2018), where the AM
model was compared to a fitting function of Ref. (Trac  _et al._ , 2015).

###  II.3 Star-formation efficiency and stellar-to-halo mass ratio

The halo growth rate discussed above provides an estimate of the amount of
total matter accreted onto a halo. However, we still need to parametrise the
star-formation efficiency, i.e., how much of that accreted matter will be
transformed into stars that emit radiation. Following Refs. (Mason  _et al._ ,
2015; Sun and Furlanetto, 2016; Mirocha  _et al._ , 2017), we define

$\displaystyle f_{*}(M)\equiv\dot{M}_{*}/\dot{M}_{\rm ac},$ (8)

where $\dot{M}_{*}$ and $\dot{M}_{\rm ac}$ are the stellar and halo accretion
rates. Note that Eq. (8) is different from the stellar-to-halo mass ratio

$\displaystyle\tilde{f}_{*}(M)\equiv M_{*}/M_{\rm ac}$ (9)

often used in the literature (see e.g. Ref. Pritchard and Loeb, 2012; Park
_et al._ , 2019). The star-formation efficiency and the stellar-to-halo mass
ratio can be connected by the integral

$\displaystyle\tilde{f}_{*}(M)=\frac{1}{M_{\rm ac}}\int f_{*}\dot{M}_{\rm
ac}dt.$ (10)

Note that for the simple case of a redshift-independent star-formation
efficiency, Eq. (10) leads to $\tilde{f}_{*}=f_{*}$. A more general study of
the relation between star-formation efficiency and stellar-to-halo mass ratio,
including explicit redshift dependences motivated by feedback processes, can
be found in Refs. (Furlanetto  _et al._ , 2017; Mirocha, 2020).

In this paper, we use the functional form

$\displaystyle
f_{*}(M)=\frac{2(\Omega_{b}/\Omega_{m})f_{*,0}}{\left(M/M_{p}\right)^{\gamma_{1}}+\left(M/M_{p}\right)^{\gamma_{2}}}\times
S(M)$ (11)

as a parametrisation for the star-formation efficiency. Eq. (11) consists of a
double-power law, multiplied with a small-scale function

$\displaystyle S(M)=\left[1+(M_{t}/M)^{\gamma_{3}}\right]^{\gamma_{4}},$ (12)

that may provide either a suppression or a boost at the truncation mass scale
$M_{t}$. A suppression at small scales could naturally occur at scales where
atomic cooling processes become inefficient. A small-scale boost, on the other
hand, could emerge due to the presence of population III stars in mini-haloes.

Note that Eq. (11) does not depend on redshift, and we can therefore set
$f_{*}=\tilde{f}_{*}$ in this paper. We will nonetheless formally distinguish
between $f_{*}$ and $\tilde{f}_{*}$ in order to avoid confusion and to
acknowledge the fact that a more realistic description of the star-formation
efficiency may well include an explicit redshift dependence.

###  II.4 Star-formation rate density and collapse fraction

With the halo mass function, accretion rate, and star formation efficiency at
hand, it is possible to calculate the mean star-formation-rate density (SFRD),
which is an important ingredient for calculating the global 21-cm signal (see
Sec. III). We define the SFRD as the integral over the halo mass function
weighted by the star-formation efficiency and the accretion rate, i.e.,

$\displaystyle\dot{\rho}_{*}(z)=\int dM\frac{dn}{dM}f_{*}(M)\dot{M}_{\rm
ac}(M,z).$ (13)

The definition above differs from the relation
$\dot{\tilde{\rho}}_{*}(z)=\bar{\rho}df_{\rm coll}/dt$, which is based on the
global collapse fraction ($f_{\rm coll}$) and is often used in the literature
(e.g. Refs. Barkana and Loeb, 2005a; Pritchard and Furlanetto, 2006). The
collapse fraction is given by the integral

$\displaystyle f_{\rm coll}(z)=\frac{1}{\bar{\rho}}\int
dM\tilde{f}_{*}(M)\frac{dn}{dM}M,$ (14)

providing the total ratio of stars to matter in the universe. Note that Eq
(14) includes the stellar-to-halo mass ratio inside of the integral and is
therefore slightly different from the standard definition in the literature
(see e.g. Refs. Barkana and Loeb, 2004; Mirocha, 2014).

The star-formation rate densities based on the three different mass accretion
models introduced above are shown as black lines in the right-hand panel of
Fig. 1. While they have similar general trends, the exponential growth (EXP)
model is about a factor of two larger than the method based on abundance
matching (AM). The extended Press-Schechter (EPS) model lies in between, being
closer to the AM method at very high and closer to the EXP model at lower
redshifts.

The differences between the star-formation rate densities shown in Fig. 1
directly affect the modelling of the Lyman-$\alpha$ coupling, gas heating, and
ionisation. This is one of the main reasons why comparisons between different
methods to calculate the 21-cm signal are hard to interpret when they are
using different implementations for the star-formation rate (as we will see in
Sec. V).

##  III 21-cm brightness temperature

The 21-cm differential brightness temperature ($T_{21}$) is a function of the
background radiation ($T_{\gamma}$), the spin temperature of the gas
($T_{s}$), the neutral hydrogen fraction ($x_{\rm HI}$) and the gas density
field ($\delta_{b}$) which all depend on redshift $z$ and position
$\mathbf{x}$. Following e.g. Ref. (Furlanetto, 2006), the brightness
temperature can be written as

$T_{21}(\mathbf{x},z)=27x_{\rm HI}(1+\delta_{b})\\\
\times\left(\frac{\Omega_{b}h^{2}}{0.023}\right)\left(\frac{0.15}{\Omega_{m}h^{2}}\frac{(1+z)}{10}\right)^{\frac{1}{2}}\left(1-\frac{T_{\gamma}}{T_{s}}\right)$
(15)

in milli-Kelvin [mK]. Assuming a standard $\Lambda$CDM model without exotic
radio sources, the background temperature is dominated by the cosmic microwave
background (CMB) radiation. The rightmost expression in Eq. (15) can be
written as

$\displaystyle\left(1-\frac{T_{\gamma}}{T_{s}}\right)\simeq\frac{x_{\rm
tot}}{1+x_{\rm tot}}\left(1-\frac{T_{\gamma}}{T_{k}}\right),$ (16)

where $x_{\rm tot}\equiv x_{\alpha}+x_{c}$ ($x_{\alpha}$ and $x_{c}$ denoting
the radiative and collisional coupling coefficients). Throughout this paper,
we set the collisional coupling to zero, as it is only important at very high
redshifts beyond $z\sim 30$.

The gas temperature ($T_{k}$) is obtained via the differential equation

$\displaystyle\frac{3}{2}\frac{dT_{k}(\mathbf{x},z)}{dz}=\frac{T_{k}(\mathbf{x},z)}{\rho(\mathbf{x},z)}\frac{d\rho(\mathbf{x},z)}{dz}-\frac{\Gamma_{h}(\mathbf{x},z)}{k_{B}(1+z)H},$
(17)

where $\rho$ is the matter density and $\Gamma_{h}$ the heating source term.
The latter is given by the sum

$\Gamma_{h}(\mathbf{x},z)=4\pi\sum_{i}f_{i}f_{X,h}\\\ \times\int_{\nu_{\rm
th}^{i}}^{\infty}d\nu\left(\nu-\nu_{\rm
th}^{i}\right)h_{P}\sigma_{i}(\nu)J_{X,\nu}(\mathbf{x},z)$ (18)

with $i=\\{\rm H,He\\}$, i.e. the hydrogen and helium components with
fractions $f_{i}=n_{i}^{0}/n_{b}^{0}$, threshold energies $\nu_{\rm
th}^{i}h_{P}$={13.6, 26.5} eV, and cross sections $\sigma_{i}(\nu)$. For the
fraction of X-ray energy deposited as heat, we assume
$f_{X,h}=\bar{x}_{e}^{0.225}$, where $\bar{x}_{e}$ is the free electron
fraction (see Fig. 4 in Ref. Furlanetto and Stoever, 2010) that is calculated
according to Eqs. 12 and 13 in Ref. (Mirocha, 2014).

The radiation coupling coefficient, induced by the Wouthuysen-Field effect
(Wouthuysen, 1952; Field, 1958), can be written as

$\displaystyle x_{\alpha}(\mathbf{x},z)=\frac{1.81\times
10^{11}}{(1+z)}S_{\alpha}J_{\alpha}(\mathbf{x},z),$ (19)

where $S_{\alpha}$ is given by Eq. (55) in Ref. (Furlanetto  _et al._ , 2006).

The gas temperature and Lyman-$\alpha$ coupling depend on the flux terms
$J_{X}$ and $J_{\alpha}$. While the former is dominated by X-ray radiation
between a few hundred to a few thousands eV, the latter stems from the narrow
spectral range between the Lyman-$\alpha$ and the Lyman-limit frequencies.

We parametrise the spectral energy distributions of the Lyman-$\alpha$ and
X-ray flux as simple power laws

$\displaystyle I_{s}(\nu)=A_{s}\nu^{-\alpha_{s}},$ (20)

with $s=\\{\alpha,X\\}$. The normalisation $A_{s}$ is defined so that
integrating $I_{s}(\nu)$ over the corresponding energy range becomes unity
(see Ref. Mirocha  _et al._ , 2015). The number emissivity of UV photons
between the Lyman-$\alpha$ and Lyman-limit range is given by

$\displaystyle\varepsilon_{\alpha}(\nu)=\frac{N_{\alpha}}{m_{p}}I_{\alpha}(\nu),$
(21)

where $N_{\alpha}$ is the number of photons per baryon emitted in the range
between the Lyman-$\alpha$ ($\nu_{\alpha}$) and the Lyman-limit ($\nu_{LL}$)
frequencies. The energy emissivity of X-ray photons is

$\displaystyle\varepsilon_{X}(\nu)=f_{X}c_{X}\frac{I_{X}(\nu)}{\nu h_{P}}$
(22)

where $f_{X}$ is a free parameter of order unity and $c_{X}$ is a
normalisation factor constrained by observations. It is set to
$c_{X}=3.4\times 10^{40}$ erg yr s-1M${}_{\odot}^{-1}$ based on the findings
of Ref. (Gilfanov  _et al._ , 2004).

###  III.1 Global Signal

The global differential brightness temperature is directly obtained by
averaging Eqs. (15) and (16). We thereby set $\delta_{b}$ to zero and assume
spatially averaged values for the gas temperature and coupling coefficient.
The Lyman-$\alpha$ coupling is obtained via Eq. (19), where the mean
Lyman-$\alpha$ flux is given by

$\displaystyle{\bar{J}}_{\alpha}(z)=\frac{(1+z)^{2}}{4\pi}\sum_{n=2}^{n_{m}}f_{n}\int_{z}^{z_{\rm
max}^{(n)}}dz^{\prime}\frac{c\varepsilon_{\alpha}(\nu^{\prime})}{H(z^{\prime})}{\dot{\rho}}_{*}(z^{\prime}),$
(23)

with $\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$. The recycling fractions $f_{n}$
are taken from Ref. (Pritchard and Furlanetto, 2006) with the sum truncated at
$n_{m}=23$. The integration limit is given by

$\displaystyle z_{\rm max}^{(n)}=(1+z)[1-(n+1)^{-2}]/(1-n^{-2})-1,$ (24)

designating the maximum redshift from which photons can Doppler-shift into the
Lyman resonances.

The global X-ray number flux per frequency is given by the relation (Pritchard
and Loeb, 2012)

$\displaystyle{\bar{J}}_{X,\nu}(z)=\frac{(1+z)^{2}}{4\pi}\int_{z}^{\infty}dz^{\prime}\frac{c\varepsilon_{X}(\nu^{\prime})}{H(z^{\prime})}{\rm
e}^{-\tau_{\nu}}{\dot{\rho}}_{*}(z^{\prime}).$ (25)

Here we have introduced the optical depth parameter $\tau$ defined as

$\displaystyle\tau_{\nu}(z,z^{\prime})=\int_{z}^{z^{\prime}}dz^{\prime\prime}\frac{dl}{dz^{\prime\prime}}\sum_{i}n_{i}\sigma_{i}(\nu^{\prime\prime}),$
(26)

where $i=\\{{\rm HI,HeI}\\}$ and where $\nu^{\prime\prime}$ is the frequency
redshifted from the source at $z^{\prime}$ to $z^{\prime\prime}$. Plugging Eq.
(25) into Eq. (18) leads to the mean heating rate $\bar{\Gamma}_{h}(z)$. The
global temperature evolution is finally obtained by solving

$\displaystyle\frac{3}{2}\frac{d\bar{T}_{k}}{dz}=\frac{3\bar{T}_{k}}{(1+z)}-\frac{\bar{\Gamma}_{h}(z)}{k_{B}(1+z)H},$
(27)

which corresponds to Eq. (17) at order zero in density perturbations.

The solution of Eq. (27) together with $\bar{x}_{\alpha}$ from Eq. (23) allows
us to obtain the global differential brightness temperature ($\bar{T}_{21}$).
Examples of the $\bar{T}_{21}$ signal are shown in the bottom-left panel of
Fig. 4.

###  III.2 Power spectrum

The 21-cm brightness temperature defined in Eq. (15) is a function of both
redshift and position. At linear order, the 21-cm perturbations are given by

$\displaystyle\delta_{21}(\mathbf{x},z)=\beta_{b}\delta_{b}+\beta_{\alpha}\delta_{\alpha}+\beta_{h}\delta_{h}+\beta_{p}\delta_{p}-\delta_{dv},$
(28)

This expansion is identical to the one proposed in Ref. (Barkana and Loeb,
2004), except that, for reasons that will become evident later on (see Sec.
IV.3), we furthermore separate the temperature fluctuations into a heating
($\delta_{h}$) and a primordial ($\delta_{p}$) term:

$\displaystyle\delta_{T}=f_{\rm T}\delta_{h}+(1-f_{\rm T})\delta_{p},$ (29)

where $f_{\rm T}\equiv(\bar{T}_{k}-\bar{T}_{p})/\bar{T}_{k}$. In this context
_primordial_ means prior to any heating from sources, i.e., the regime where
the gas cools adiabatically and the temperature fluctuations are seeded by the
matter perturbations. Note furthermore that Eq. (29) is a direct consequence
of the assumption $T_{k}=T_{h}+T_{p}$.

The other terms of Eq. (28) designate the baryon perturbations ($\delta_{b}$),
the Lyman-$\alpha$ coupling perturbations ($\delta_{\alpha}$), and the
perturbations due to the line-of-sight velocity gradient ($\delta_{dv}$)
caused by redshift-space distortion effects. At linear order and in Fourier
space, $\delta_{dv}$ is simply given by $\delta_{dv}=\mu\delta_{m}$, where
$\mu$ is the cosine of the angle between the wave vector $\mathbf{k}$ and the
line of sight (Bharadwaj and Ali, 2004).

The pre-factors of the individual perturbations in Eq. (28) are given by

$\displaystyle\beta_{b}$ $\displaystyle\simeq$ $\displaystyle 1,$ (30)
$\displaystyle\beta_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\bar{x}_{\alpha}}{\bar{x}_{\rm tot}(1+\bar{x}_{\rm
tot})},$ (31) $\displaystyle\beta_{h}$ $\displaystyle\simeq$ $\displaystyle
f_{\rm T}\frac{\bar{T}_{\gamma}}{(\bar{T}_{k}-\bar{T}_{\gamma})},$ (32)
$\displaystyle\beta_{p}$ $\displaystyle\simeq$ $\displaystyle(1-f_{\rm
T})\frac{\bar{T}_{\gamma}}{(\bar{T}_{k}-\bar{T}_{\gamma})},$ (33)

and only depend on redshift $z$ but not on the position vector $\mathbf{x}$.

Based on Eqs. (28) and (29), it is straight-forward to calculate the power
spectrum and to sort all components with respect to their power of $\mu$.
Taking the average over the angle then leads to (Barkana and Loeb, 2005b)

$P_{\rm 21}=P_{\rm\alpha\alpha}+P_{hh}+P_{pp}+P_{bb}\\\
+2\left(\frac{}{}P_{\alpha h}+P_{\alpha p}+P_{\alpha
b}+P_{hp}+P_{hb}+P_{pb}\right)\\\ +\frac{2}{3}\left(\frac{}{}P_{\alpha
m}+P_{hm}+P_{pm}+P_{bm}\right)+\frac{1}{5}P_{mm}.$ (34)

In the next section, we propose a method to calculate all individual auto and
cross power spectra of Eq. (34) in order to obtain a fast, analytical estimate
for the power spectrum of the 21-cm brightness fluctuations. Note furthermore,
that when talking about the 21-cm power spectrum, we will mean either $P_{21}$
or the expression $\bar{T}^{2}\Delta_{21}^{2}$ depending on the context. The
dimnesionless quantity $\Delta_{21}^{2}$ is defined as $\Delta_{21}^{2}\equiv
k^{3}P_{21}/(2\pi^{2}$).

##  IV Halo model

In the framework of the original halo model haloes are considered as the
building blocks of the universe (see e.g. Ref. Cooray and Sheth, 2002). The
halo abundance, distribution, and internal profiles are used to calculate the
matter power spectrum at both linear and nonlinear scales. In the context of
21-cm clustering, however, the focus is not so much on haloes per se, but
rather on the sources inhabiting haloes, which emit radiation, thereby
affecting the gas cells around them. Instead of halo profiles, we therefore
consider radiation profiles around sources that extend into the intergalactic
space far beyond the halo limits. This means that there will typically be
several overlapping radiation profiles from different sources affecting any
single gas volume.

Although the picture of overlapping radiation profiles is rather different
from the original halo model, it turns out that the 21-cm power spectrum can
be described in a very similar way. In this section, we first introduce the
formalism, before going into the details of the source profiles and the
description of temperature fluctuations. At the end we show the resulting
21-cm power spectrum assuming three benchmark models with different
astrophysical parameters.

Figure 2: Normalised Lyman-$\alpha$ coupling (left), X-ray energy deposition
(centre), and heating (right) profiles. In the top-row we vary the redshift at
fixed halo mass and in the bottom-row we vary the halo mass at fixed redshift.

###  IV.1 Power spectrum description

In the context of the 21-cm halo model, the power spectra of different
components can be calculated in the following way:

$\displaystyle P_{XY}^{\rm 1h}(k,z)=$
$\displaystyle\frac{\beta_{X}\beta_{Y}}{({\bar{\rho}}f_{\rm coll})^{2}}\int
dM\frac{dn}{dM}\tilde{f}_{*}^{2}M^{2}|u_{X}||u_{Y}|,$ $\displaystyle
P_{XY}^{\rm 2h}(k,z)=$ $\displaystyle\frac{\beta_{X}}{{(\bar{\rho}}f_{\rm
coll})}\int dM\frac{dn}{dM}\tilde{f}_{*}M|u_{X}|b_{X}$ (35)
$\displaystyle\times$ $\displaystyle\frac{\beta_{Y}}{({\bar{\rho}}f_{\rm
coll})}\int dM\frac{dn}{dM}\tilde{f}_{*}M|u_{\rm Y}|b_{Y}\times P_{\rm lin},$
$\displaystyle P_{XY}(k,z)=$ $\displaystyle P_{XY}^{\rm 1h}(k,z)+P_{XY}^{\rm
2h}(k,z),$

where $P_{\rm lin}(k)$ the linear matter power spectrum, $u_{X,Y}(k,M,z)$ the
Fourier transformed flux profile, and $b_{X,Y}(M,z)$ the halo bias. The
subscripts $X$ and $Y$ refer to the mass ($m$), baryon ($b$), Lyman-$\alpha$
($\alpha$), and heating ($h$) components. Note furthermore that Eq. (IV.1)
describes both auto and cross spectra depending on whether $X=Y$ or $X\neq Y$.
This means that the halo model provides all components of Eq. (34) except the
ones that include the primordial gas temperature ($p$). We will derive these
in Sec. IV.3.

In the halo model framework, the one-halo term ($P_{XY}^{\rm 1h}$) describes
the clustering within one single source profile, while the two-halo term
($P_{XY}^{\rm 2h}$) accounts for the signal induced by different sources. As a
consequence $P_{XY}^{\rm 1h}$ and $P_{XY}^{\rm 2h}$ dominate at small and
large scales, respectively, with a transition region corresponding to the
typical size of the source profiles. Not surprisingly, only the two-halo term
carries information about the spatial distribution of sources via the
components $b_{X,Y}$ and $P_{\rm lin}$. The one-halo term, on the other hand,
carries information about the shot-noise of sources, and its shape is only
controlled by the Fourier transformed radiation profiles.

###  IV.2 Flux profiles

The key components of the halo model are the radiation profiles around
sources. Following Holzbauer and Furlanetto (2011), the profile of the
Ly-$\alpha$ radiation can be written as

$\displaystyle\rho_{\alpha}(r|M,z)=\frac{1}{4\pi
r^{2}}\sum_{n=2}^{n_{m}}f_{n}\varepsilon_{\alpha}(\nu^{\prime})f_{*}\dot{M}_{\rm
ac}(z^{\prime}|M,z)$ (36)

with $\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$. The profile is proportional to
the halo mass accretion $\dot{M}_{\rm ac}$ and exhibits the characteristic
$r^{2}$ decrease with radius $r$. It furthermore depends on the look-back
redshift $z^{\prime}=z^{\prime}(r)$ that corresponds to the redshift when a
photon at radius $r$ has been emitted at the source. At the emission redshift
$z^{\prime}$, the source is in an earlier stage of evolution (compared to the
redshift of the signal $z$), which means that its accretion rate is smaller as
well. The look-back redshift is obtained by inverting the co-moving distance

$\displaystyle
r(z^{\prime}|z)=\int_{z}^{z^{\prime}}\frac{c}{H(z^{\prime\prime})}dz^{\prime\prime}$
(37)

which has to be done numerically.

Figure 3: Normalised Fourier transforms of the Lyman-$\alpha$ (solid) and the
heating profiles (dashed). In the left- and right-hand panels we vary redshift
and halo mass while keeping the other one constant.

The shape of the Lyman-$\alpha$ flux profile is plotted in the left-hand
panels of Fig. 2, where the top and bottom panels specifically highlight the
dependencies on redshift ($z$) and halo mass ($M$). Here we have used the EPS
mass accretion model, but the plots look very similar if another model is
applied instead. All Lyman-$\alpha$ profiles are characterised by a $1/r^{2}$
decrease close to the source, which becomes gradually steeper towards the
outer parts. The steepening is a direct result of the radiation originating
from the source at a higher redshift ($z^{\prime}$) when the accretion rate
onto the source was smaller. At around 200 Mpc/h, the profiles exhibit a steep
drop that is due to the photons having redshifted out of the Lyman-$\alpha$
series. At very high redshift beyond z=25, the drop-off shifts towards much
smaller radii below 100 Mpc/h because the source is so young that the
radiation did not have time to expand further. Note that this effect is not
limited to the Lyman-$\alpha$ radiation, but is also visible in the X-ray and
heating profiles discussed below.

A closer inspection of the Lyman-$\alpha$ flux profiles in Fig. 2 reveals
small discontinuities in the form of step-like features between $r=0.1-100$
Mpc/h. These small steps are a consequence of the sum in Eq. 36 and have been
predicted in earlier work (see Ref. Holzbauer and Furlanetto, 2011).

The heating of the gas depends on the flux of X-ray radiation and can be
defined similarly to the Lyman-$\alpha$ profile, i.e.,

$\rho_{\rm xray}(r|M,z)=\frac{1}{r^{2}}\sum_{i}f_{i}f_{X,h}\\\
\times\int_{\nu_{\rm th}^{i}}^{\infty}d\nu\left(\nu-\nu_{\rm
th}^{i}\right)h_{P}\sigma_{i}(\nu)\varepsilon_{\rm
X}(\nu^{\prime})e^{-\tau_{\nu}}f_{*}\dot{M}_{\rm ac}(z^{\prime}|M,z),$ (38)

where $i$=$\\{\rm H,He\\}$ and $\nu_{\rm th}^{i}h_{P}$={13.6, 26.5} eV.
Compared to Eq. (36) there is an additional attenuation term due the optical
depth defined in Eq. (26).

The shape of the energy deposition profile ($\rho_{\rm xray}$) is shown in the
middle panels of Fig. 2, the top and bottom panels again highlighting its
redshift and mass dependencies. As a general trend, the profile extends
further in radius (to around 500 Mpc/h) and is more gradually suppressed than
the Lyman-$\alpha$ coupling profile. This is due to the fact that hard X-ray
radiation can travel large distances until it is deposited as heat, and X-ray
photons do not redshift out of a well defined spectral range. Since the large
distances traveled by X-ray photons also lead to increased look-back redshifts
(the difference between $z$ and $z^{\prime}$), the energy deposition profile
shows a stronger attenuation towards large radii.

The energy deposition from X-ray emission leads to a strong increase of the
temperature fluctuations around sources. The corresponding profile can be
obtained via the differential equation

$\frac{3}{2}\frac{d\rho_{h}(r|M,z)}{dz}=\frac{3\rho_{h}(r|M,z)}{(1+z)}-\frac{\rho_{\rm
xray}(r|M,z)}{k_{B}(1+z)H}.$ (39)

The first term on the right-hand-side of this equation describes the cooling
due to the expansion of space, the second term corresponds to the energy
deposition from X-ray radiation. The solution of Eq. (39) is called the
heating profile ($\rho_{h}$) in agreement with the notation of Sec. III.2,
where we separated the temperature fluctuations into a primordial and a
heating term (see Eq. 29).

The heating profile is illustrated on the right-hand-side of Fig. 2, where the
top and bottom panels again show the dependency with redshift and halo mass.
The shape of the heating profile ($\rho_{h}$) is similar to the energy
deposition profile ($\rho_{\rm xray}$), which is not a surprise, since X-ray
radiation is assumed to be the only source of heating in Eq. (39).

The source profiles of Eqs. (36) and (39) are key ingredients of the halo
model. However, in order to substitute them into Eq. (IV.1), we need to first
calculate their normalised Fourier transforms. They are given by

$\displaystyle u_{i}(k|M,z)=\frac{\int drr^{2}\rho_{i}(r|M,z)j_{0}(kr)}{\int
drr^{2}\rho_{i}(r|M,z)}$ (40)

where $j_{0}(x)=\sin(x)/x$ is the spherical Bessel function of order 0. The
subscript $i$ either stands for $\alpha$ or $h$, denoting the Lyman-$\alpha$
emission and heating profiles.

In Fig. 3 we show the Fourier transformed profiles with their redshift and
halo mass dependencies. All profiles go to unity at low $k$-modes, a key
characteristic that is guaranteed by the normalisation of Eq. (40). Towards
higher values of $k$, the profiles become strongly suppressed. In general, the
Lyman-$\alpha$ profiles are more suppressed than the heating profiles, which
is a consequence of the fact that soft X-ray photons do not travel far before
being absorbed by the gas. This early absorption leads to an excess of small-
scale clustering compared to the more freely emitted Lyman-$\alpha$ flux. One
exception to this behaviour becomes visible at $z=30$ (see red lines in the
left-hand panel), where both profiles have a very similar shape. This is not
surprising because at very high redshifts, the shape of the profiles is driven
by the emission of a very young source with photons that did not have time to
travel far.

Fig. 3 shows a clear redshift dependence of the normalised Fourier profiles
with more small scale clustering towards higher redshifts. The mass
dependence, on the other hand, is very weak. This means that, in principle,
the halo model of Eq. (IV.1) could be simplified considerably by assuming the
profiles $u_{i}$ not to depend on halo mass. Note, however, that for the sake
of completeness we keep the full mass dependence in our model.

###  IV.3 Temperature fluctuations

In Sec. III.2, we have separated the temperature fluctuations into a heating
($\delta_{h}$) and a primordial ($\delta_{p}$) term. While the heating term is
sourced by the X-ray flux emission and can therefore be readily described by
the halo model, the primordial fluctuations are driven by the matter
fluctuations and can be directly solved via Eq. (17) at the linear level.
Since we have separated out the contribution from the sources, we can set
$\Gamma=0$. After linearising $T=\bar{T}(1+\delta_{p})$ as well as
$\rho=\bar{\rho}(1+\delta)$ we obtain the solution

$\displaystyle\delta_{p}=(1+\delta)^{2/3}-1.$ (41)

for the primordial (prior to sources) temperature fluctuations. Note that in
Fourier space, the matter perturbations are readily obtained by setting
$\delta=\sqrt{P_{\rm mm}}$. With this at hand, we can now derive the reminding
auto and cross power spectra of Eq. (34) which include the primordial heating
term and are thus not covered by Eq. (IV.1). They are given by

$\displaystyle P_{pp}(k,z)$ $\displaystyle=$
$\displaystyle\beta_{p}(z)^{2}\delta_{p}(k,z)^{2},$ (42) $\displaystyle
P_{Xp}(k,z)$ $\displaystyle=$
$\displaystyle\beta_{p}(z)\sqrt{P_{XX}(k,z)}\delta_{p}(k,z),$ (43)

where $X$ again stands for $\\{\alpha,h,b,m\\}$.

Figure 4: _Leftmost panels:_ Star-formation efficiency (top) and resulting
global signal (bottom) for the three benchmark models with suppressed,
unchanged and boosted small-scale behaviour (blue, cyan, and magenta).
_Remaining panels:_ Power spectra as a function of $k$-modes (top) and
redshift (bottom) for the same models.

###  IV.4 Results

Based on the formalism derived above, it is now possible to calculate the
21-cm power spectrum for a given set of model parameters. In this section, we
will first present our choices of parametrisation before showing the power
spectrum for a selected set of redshift and $k$-modes.

For the Lyman-$\alpha$ sources, we assume $N_{\alpha}=10000$ evenly
distributed ($\alpha_{\alpha}=0$) over the energy range between the
Lyman-$\alpha$ and Lyman-limit frequencies (see Eq. 21). The X-ray energy
emission is defined by $f_{X}=1$ plus a power-law spectral energy distribution
with $\alpha_{X}=1.5$ over the range $E=0.5-5$ keV (see Eq. 22). Regarding the
halo mass function, we assume the Sheth and Tormen (2002) prescription (Eq. 1)
with modified parameter $q=0.85$. The halo bias is modelled based the peak-
background split approach (Eq. 4) with the same parameters than what is used
for the halo mass function.

In order to highlight the sensitivity of the results to the source
parametrisation, we focus on three benchmark models (A, B, and C) that are
characterised by different star formation efficiencies ($f_{*}$). All models
have the same large-scale behaviour given by the double power law of Eq. (11)
with $f_{*,0}=0.3$, $\gamma_{1}=0.49$, $\gamma_{2}=-0.61$, and $M_{p}=2\times
10^{11}$ M⊙/h (see Ref. Mirocha  _et al._ , 2017). At small mass scales,
however, Model A assumes a strong additional suppression (with $M_{t}=5\times
10^{7}$ M⊙/h, $\gamma_{3}=1$, $\gamma_{4}=-4$) mimicking the effects of
inefficient cooling processes. Model C, on the other hand, is characterised by
a boost of $f_{*}$ towards very small masses (with $M_{t}=10^{7}$ M⊙/h,
$\gamma_{3}=1$, $\gamma_{4}=1$). Such a behaviour can be motivated by the
presence of Population-III stars in mini-haloes. Model B finally shows neither
additional suppression nor boost of $f_{*}$ but a continuation of the power-
law decrease down to the smallest masses (i.e. $\gamma_{4}=0$). All models are
truncated at $M_{\rm min}=5\times 10^{5}$ M⊙/h (which roughly corresponds to
the smallest halo mass where stars can form via molecular gas cooling).

The star-formation efficiencies of the three benchmark models are plotted in
the top-left panel of Fig. 4. The typical shape from the double-power law
prescription is visible at large masses above $\sim 10^{9}$ M⊙/h. At smaller
halo masses the models diverge showing the characteristic suppression, power-
law continuation, and boost of the benchmark models A, B, and C described
above.

The effect of the different star-formation efficiencies on the global
differential brightness signal is illustrated in the bottom-left panel of Fig.
4. While model A leads to a narrow absorption signal at a rather low redshift
of $z\sim 14$, model B and C show wider troughs shifted towards $z\sim 17$ and
$z\sim 21$, respectively 111Note that although our parametrisation is similar
to that of Mirocha  _et al._ (2017), the global signal is shifted to slightly
higher redshift due to differing assumptions about the star-formation
efficiency normalisation and the minimum halo mass..

The 21-cm power spectra as a function of $k$-values are plotted in the
remaining top-panels of Fig. 4 with increasing redshift from left to right.
Many of the lines have a wave like feature with a local flattening or minimum
at $k\sim 0.1-1$ h/Mpc. This feature is located at the transition where the
1-halo starts to dominate over the 2-halo term.

The remaining panels at the bottom of Fig. 4 show the power spectra as a
function of redshift for the three specific $k$-modes. At $k=$ 0.01 h/Mpc
(left) and 0.1 h/Mpc (middle) the power spectra show the characterised double-
peak feature, which indicate the two characteristic epochs where the
Lyman-$\alpha$ coupling and the X-ray heating dominate. At $k=1$ h/Mpc the two
peaks have merged into one single, broader peak in qualitative agreement with
other work from the literature (see e.g. Refs. Santos  _et al._ , 2008; Ross
_et al._ , 2019).

The results presented in Fig. 4 are based on the EPS mass accretion model plus
the halo mass function and bias parameters described above. However, it is
important to notice that these modelling choices introduce significant
uncertainties regarding the 21-cm power spectrum. In the Appendices A and B,
we quantify the effects due to the choice of the mass accretion model, the
halo mass function and halo bias. The bottom line of this analysis is that
modelling choices, such as switching from Sheth-Tormen to Press-Schechter halo
prescription, or using the AM instead of the EPS accretion rate modelling, has
an effect on the 21-cm power spectrum that can in some cases be as large as a
factor of $\sim 10$. The difference can be significantly reduced if the flux
parameters $N_{\alpha}$ are re-adapted so that the global signals are forced
to match, but a remaining difference of the power spectrum of about a factor
of $\sim 2$ remains.

##  V Comparison with other approaches

In this section we compare the 21 cm halo model with the analytical model from
Refs. (Barkana and Loeb, 2005a; Pritchard and Furlanetto, 2007) (abbreviated
to BLPF model) and with the semi-numerical code 21cmFAST (Mesinger and
Furlanetto, 2007; Mesinger  _et al._ , 2010). While the former is more of a
general consistency check, the latter consist of a true test for our model.
Note, however, that the comparison with 21cmFAST is not straight-forward,
mainly because there are subtle differences in the parametrisation of source
properties which cannot be fully accounted for without changing the code
itself, something we postpone to future work.

Figure 5: Comparison between the 21-cm halo model developed in this paper
(solid lines) and the BLPF method from Refs. (Barkana and Loeb, 2005a;
Pritchard and Furlanetto, 2007) (coloured bands) assuming astrophysical
parameters from the benchmark model B (see Fig 4). _Top panels:_ global signal
(left) and the total 21-cm power spectrum with respect to $k$-modes (centre)
and redshift (right). _Middle panels:_ Auto power spectra of the
Lyman-$\alpha$ coupling (left), temperature (centre), and mass (right)
components. _Bottom panels:_ Corresponding cross power spectra of the
individual components.

###  V.1 Comparing to the analytical approach of BLPF

A first analytical calculation of the 21-cm power spectrum at cosmic dawn has
been performed by Barkana and Loeb (2005a), focusing on the perturbations
induced by the Lyman-$\alpha$ coupling. Their model has been extended to
include temperature fluctuations by Pritchard and Furlanetto (2007) which is
why we abbreviate it as BLPF. In contrast to the halo model approach, which is
centred on the radiation sources as building blocks, the BLPF approach focuses
on the gas, calculating the light-cone effects from the surrounding sources. A
good summary of the model can be found in Ref. (Pritchard and Loeb, 2012).

For this comparison we have implemented the BLPF model, mostly following the
descriptions in Ref. (Pritchard and Furlanetto, 2007). The source
parametrisation, however, which includes the halo mass function, the bias, and
the spectral energy distribution, has been adapted to the description
presented in this paper. We furthermore omit any shot-noise corrections, since
there is no fully worked-out model that includes both shot-noise of the
Lyman-$\alpha$ and temperature fluctuations (Pritchard and Furlanetto, 2007).
Finally, the star-formation rate density is based on the EPS accretion rate
modelling (see Eq. 13) and not on the time derivative of the collapse
fraction. These changes with respect to the original work of Refs. (Barkana
and Loeb, 2005a; Pritchard and Furlanetto, 2007) allow us to carry out a fair
comparison, where any resulting discrepancies can be fully attributed to
differences between the methods and are not the result of different source
descriptions.

The comparison between the BLPF and the halo model is performed using the
source parametrisation of benchmark model B. This means we assume the
radiation flux parameters $N_{\alpha}=10000$ and $f_{X}=1$, a non-truncated
double-power law for the star formation efficiency, as well as power-law
spectra with indices $\alpha_{\alpha}=0$ and $\alpha_{X}=1.5$ for the UV and
X-ray radiation. See Sec. IV.4 for more details about the parametrisation.

The global signal is shown in the top-left panel of Fig. 5. Since we use the
exact same source modelling and the same calculation for the SFRD, it is not
surprising that both approaches yield the exact same result. This perfect
agreement is very convenient because it guarantees that any changes at the
level of the power spectrum are not induced by different amplitudes of the
global differential brightness temperature.

The central-top panel of Fig. 5 shows the 21-cm power spectrum as a function
of $k$-modes for a selection of four different redshifts. The redshift values
are chosen to lie in the heating dominated regime ($z=13$), in the transition
regime of maximum absorption ($z=16$), in the Lyman-$\alpha$ regime ($z=19$),
and at the epoch of the very first stars ($z=22$). The solid lines correspond
to the 21-cm halo model, while the coloured bands show the results from the
BLPF model. In general, the BLPF model predicts less power, especially during
the heating (black) and, to a lesser extend, during the Lyman-$\alpha$ epochs
(dark green). The differences between the models are typically of the order of
a few, but they can grow to about an order of magnitude for specific redshifts
and $k$-ranges.

The right-hand panel in the top row of Fig. 5 illustrates the 21-cm power
spectrum, this time as a function of redshift for the selected modes $k\sim
0.01$, $0.06$, $0.2$, and $0.95$ h/Mpc. As before, we observe a good
qualitative agreement between the models, both of them showing a
characteristic double-peak feature for low $k$-modes merging into one single
peak at $k\sim 1$ h/Mpc. A closer look reveals, however, that the double-peak
feature survives to higher $k$-values in the halo model compared to the BLPF
model. In general, we conclude that the two models differ by no more than a
factor of a few in their redshift evolution, with some exceptions where the
difference can grow to about an order of magnitude at most.

In the middle and bottom rows of Fig. 5 we show the auto and cross power
spectrum for the individual components $\alpha$, $T$, and $m$. We observe as a
general rule that the halo model and BLPF power spectra are well converged at
the largest scales (lowest $k$ values) before they start to slowly diverge
towards higher $k$-values. Beyond $k\sim 1$ h/Mpc, the divergence is
accentuated which is due to the dominance of the one-halo term. The largest
differences are visible in the $P_{\alpha\alpha}$ and $P_{TT}$ auto spectra as
well as the $P_{\alpha T}$ cross spectrum.

We conclude that the 21-cm halo model predictions are in qualitative agreement
with the results from the BLPF model. In general, the halo model power
spectrum is larger by a factor of a few, and, for exceptional $k$ and
$z$-values the differences can grow to about an order of magnitude. We want to
emphasise, however, that at this point we do not know which of the two models
is more accurate. Although former findings by Ref. (Santos  _et al._ , 2008)
suggest that the BLPF model lacks power with respect to semi-numerical
calculations as well, we want to remind that these findings were based on a
different source parametrisation and are therefore not directly comparable.

Figure 6: Comparison between the halo model and 21cmFAST. _Top-left:_ Star
formation efficiency of Eq. 11 fitted to the one from 21cmFAST (blue line and
band, respectively). _Top-right:_ Global differential brightness temperature
obtained by assuming the EXP and EPS accretion rates (dash-dotted and solid
lines) and matched to the 21cmFast results (band) by refitting the total
Lyman-$\alpha$ and X-ray photon numbers. _Centre and Bottom:_ Corresponding
power spectra as a function of $k$-modes and redshift. Note that redshifts
below $z\sim 10$ should not be trusted since effects related to reionisation
become important.

###  V.2 Comparing to 21cmFAST

The code 21cmFAST (Mesinger  _et al._ , 2010; Murray  _et al._ , 2020) is
based on a semi-numerical approach to predict 21-cm maps by solving the spin
temperature evolution and the reionization process on a three-dimensional
grid. The matter field is evolved according to a first-order Zel’Dovich
displacement (Zel’Dovich, 1970). The sources are not resolved individually,
but their distribution is calculated using an excursion-set method directly
applied to the matter field.

A one-to-one comparison between our model and 21cmFAST is not straight-forward
because of subtle differences in the parametrisation and implementation that
may significantly affect the results. For example, we do not exactly know
which halo bias and halo mass function agrees best with the excursion set
implementation used in 21cmFAST. Furthermore, there are small differences in
the parametrisations of the spectral energy distributions and the star-
formation efficiency, that might be of relevance.

For the 21cmFAST run, we assume the fiducial model from Ref. (Park  _et al._ ,
2019) whose astrophysical parameters guarantee an agreement with the observed
luminosity function at $z=6-10$. At the same time, it leads to a pre-
reionisation signal that is shifted to somewhat lower redshifts compared to
the benchmark models introduced above. The stellar-to-halo mass ratio assumed
by Ref. (Park  _et al._ , 2019) corresponds to a power law, followed by an
exponential cutoff towards small scales (with parameters $f_{*,10}=0.05$ for
the amplitude, $\alpha_{*}=0.5$ for the power-law index, and $M_{\rm
turn}=5\times 10^{8}$ M⊙ for the mass scale of the exponential downturn, see
Eq. (2) in Ref. (Park  _et al._ , 2019)). We attempt to reproduce the same
functional form without changing the parametrisation described in Eq. (11).
The best agreement is found with the parameters $f_{*,0}=0.05$,
$M_{p}=10^{10}$ M⊙, $\gamma_{1}=-0.5$, $\gamma_{2}=-0.5$, $M_{t}=2\times
10^{8}$ M⊙, $\gamma_{3}=1.4$, and $\gamma_{4}=-4$. A comparison of the two
functions is shown in the top-left panel of Fig. 6.

Regarding the halo mass function, we use the Sheth-Tormen model with $q=0.707$
and $p=0.3$. Although this mass function is not such a good fit to high-
redshift $N$-body simulations (see Ref. (Iliev  _et al._ , 2011; Schneider,
2018) and discussions in Sec. II.1), it is used as the reference for the
collapse fraction calculated in 21cmFAST (see e.g. Eq. 14 in Mesinger  _et
al._ , 2010). For the halo bias, on the other hand, we apply Eq. (4) with the
Press-Schechter (PS) parameters $q=1$ and $p=0$. This is because in 21cmFAST,
sources are populated with respect to the PS conditional mass function.
Finally, we do not know which prescription for the halo-accretion rate is
supposed to match best with the algorithm of 21cmFAST. For this reason we
decide to show the EPS as well as the EXP approach, since both are in good
agreement with high-redshift $N$-body simulations (see Sec. II.2).

In order to allow for a meaningful comparison, we select the flux parameters
$N_{\alpha}=1400$ and $f_{X}=0.26$ so that we obtain a good match to the
21cmFAST global signal. The spectral energy range and power-law index of the
X-ray radiation, on the other hand, is kept fix at the default values of
21cmFAST (which are $E_{\rm min}=0.5$ keV, $E_{\rm max}=2$ keV, and
$\alpha_{X}=1$).

In the top-right panel of Fig. 6, we plot the global differential brightness
temperature. The results from the 21-cm halo model with EPS and EXP accretion
are shown as solid and dash-dotted lines, while the global signal from
21cmFAST is plotted as broad blue band. Although the agreement between the
lines and the band is very good (which is not so surprising, since we have
selected the flux parameters $N_{\alpha}$ and $f_{\rm X}$ to obtain the best
fit to the 21cmFAST global signal), the curves are not identical. The
differences are of order 10 percent or less, and could either stem from
unaccounted deviations in the source parametrisation or from the different
ways the global signal is calculated.

The central panel of Fig. 6 shows the power spectrum as a function of
$k$-modes for three different redshifts, representing the regime dominated by
heating ($z\sim 10.8$, purple), the transition regime where the global
absorption signal is at its maximum ($z\sim 12.9$, blue), and the regime
dominated by the Lyman-$\alpha$ coupling ($z\sim 15.0$, cyan). The results
from 21cmFAST are again shown as broad coloured bands. Note that we plot two
runs with box-size $L=300$ and $450$ Mpc, respectively (where the number of
low- and high-resolution cells are kept constant at $300^{3}$ and $1200^{3}$).
The resulting power spectra from the 21-cm halo model are plotted as solid and
dashed lines.

The bottom panel of Fig. 6 illustrates the power spectrum, this time as a
function of redshift. At large scales (low $k$-values) the two characteristic
peaks due to the Lyman-$\alpha$ coupling and the heating epochs are clearly
visible at $z\sim 11$ and 14. Arounf $k\sim 0.3$ h/Mpc, a continuous
transition from two to one prominent peak is visible in both models. The
agreement between 21cmFAST and the 21-cm halo model is best at very large and
very small scales. In between, at the transition scale (see $k=0.2$ and $0.44$
h/Mpc), there is some visible discrepancies between the two models, with
21cmFAST predicting more power by a factor of a few between $z\sim 11$ and
$13$. Furthermore, there are increasing differences towards very low and very
high redshifts. Note, however, that at redshifts below $z\sim 10$ the halo
model results cannot be trusted because effects from reionisation are ignored.
At very high redshifts, on the other hand, both predictions from 21cmFAST and
the halo model become increasingly uncertain due to the poor knowledge of
source numbers and distributions.

As a summary, let us emphasise that we find good qualitative agreement between
21cmFAST and the 21-cm halo model. At the quantitative level, the differences
do not exceed a factor of $\sim 3$, except at very low and very high
redshifts, where the models cannot be fully trusted. Note that this level of
disagreement is of the same order than the expected systematic effects due to
modelling choices regarding the halo accretion, the halo abundance, and the
halo bias (see Appendix A and B for more details). We therefore conclude that,
at the current stage, it is impossible to know if the observed differences
between the 21-cm halo model and 21cmFAST are a result of different modelling
choices or if they hint towards more fundamental issues regarding the halo
model approach.

##  VI Conclusions

Fast and sufficiently accurate predictions for the 21-cm clustering signal are
important in order to develop a better understanding of the cosmic dawn, the
epoch of the high-redshift universe right before the phase transition from
neutral to ionised hydrogen. Many important aspects of the prediction
pipeline, especially related to the source modelling, remain unknown, and fast
models may help to explore the vast parameter space of possibilities.

In this paper we present a new analytical method based on the framework of the
halo model, where Lyman-$\alpha$ coupling and temperature fluctuations are
described with the help of overlying flux profiles that properly include red-
shifting and source attenuation due to the expansion of the universe and the
finite speed of light. The model provides a natural framework to predict all
auto and cross power spectra of the Lyman-$\alpha$ coupling, the temperature
fluctuations, and the matter perturbations. The temperature fluctuations are
separated in a primordial component sourced by the matter perturbations and a
heating component induced by the first sources. The effects form the process
of reionisation could, in principle, be added to the framework, but this is
left as future work.

A distinctive advantage compared to other analytical methods is that the halo
model approach naturally includes shot-noise effects induced by the low number
of sources in the early universe. This is important for high redshifts and
small scales, especially when investigating models with very luminous sources.
The influence of shot-noise on the power spectrum is model-dependent and can
only be quantified in combination with the assumed astrophysical parameters.

Next to presenting the framework of a 21-cm halo model, we investigate the
effects of the stellar accretion, the source abundance, and the halo bias. We
show that different choices regarding these model ingredients can strongly
affect both the global signal absorption trough and the 21-cm power spectrum.
If the flux parameters are re-fitted to correct for the shift in the global
signal, the effect on the power spectrum becomes smaller, but the difference
can still be as large as a factor of a few. This means that analytical and
semi-numerical methods require a precise understanding of structure formation
at very early times in order to avoid significant systematic errors.

In order to check the general validity of the 21-cm halo model, we compare it
to the earlier analytical model by Barkana, Loeb, Pritchard, and Furlanetto
(BLPF) introduced in Refs. (Barkana and Loeb, 2005a; Pritchard and Furlanetto,
2007). In general, we find good qualitative agreement between the two
approaches regarding both the $k$-mode and redshift evolution. At the
quantitative level, the halo model predicts a somewhat larger clustering
signal compared to the BLPF model. However, the difference is no more than a
factor of a few, except for specific redshifts and $k$-values where it can
grow to about an order of magnitude.

Furthermore, we compare the 21-cm halo model to the semi-numerical code
21cmFAST (Mesinger  _et al._ , 2010; Murray  _et al._ , 2020). Although an
exact comparison between the two methods is currently unfeasible due to small
differences in the source parametrisations, we nevertheless find a very
encouraging agreement. At redshifts before the onset of reionisation, the
difference between the halo model and 21cmFAST stays below a factor of $\sim
3$ (being considerably better than this at most $k$-modes and redshifts).
Whether the remaining differences between the models are due to the specifics
of the source modelling, or whether they are a sign for a more fundamental
failure of the halo model remains to be investigated in the future.

In general, a detailed apple-to-apple comparison including analytical, semi-
numerical, and full simulation-based calculations of the 21-cm signal would be
extremely beneficial for the community. Only a combined use of accurate but
very expensive simulations together with much faster approximate methods will
lead to a more complete understanding of the complex and rich signal from the
epoch of the cosmic dawn.

###### Acknowledgements.

We thank Romain Teyssier and Kai Hofmann for very helpful suggestions related
to the model. This work is supported by the Swiss National Science Foundation
via the grant PCEFP2_181157.

Figure 7: Mass accretion modelling and how it affects the 21-cm signal
assuming the source parameters of fiducial model B. _Top-left:_ Halo growth
for the three accretion models EPS, EXP, and AM (see Sec. II.2). _Top-centre:_
Star-formation rate density for the same models. _Top-right:_ Resulting global
21-cm signal, where the flux parameters of the EXP and AM models
($N_{\alpha}$, $f_{X}$) are either kept the same (solid lines) or where they
are modified so that the minimum of the absorption signal lie at exactly the
same redshift the the one from the EPS model (dashed lines). _Centre and
Bottom:_ Resulting power spectra as a function of $k$-modes and redshift.

##  Appendix A Effects from the mass accretion modelling

In Sec II.2 we have investigated different modelling choices for the halo
growth rate based on an extended Press-Schechter (EPS) prescription, an
exponential accretion rate (EXP), and an abundance-matching technique (AM). We
have compared these cases to simulations of B20 (Behroozi  _et al._ , 2020)
before selecting the EPS model as our standard method for the paper.

In this appendix we have a closer look at the EPS, EXP, and AM halo accretion
rates, focusing on how they affect the 21-cm global signal and power spectrum.
We thereby assume the fiducial model B for the source parameters. More details
about model B can be found in Sec. IV.4.

The top panels of Fig. 7 show the halo growth (left), the star-formation rate
density (centre), and the global 21-cm signal (right) for the EPS, EXP, and AM
models. Note that the absorption trough of the global signal is shifted
towards higher redshifts when going from the AM to the EXP prescription.

In order to discriminate between effects originating from the global signal
and the ones affecting the power spectrum, we allow the Lyman-$\alpha$ and
X-ray flux parameters $N_{\alpha}$ and $f_{X}$ to be renormalised so that the
global signal absorption troughs of the AM and EXP models are aligned with the
ones from the EPS model. The resulting signal from such a renormalised flux
analysis are given by the dashed lines in Fig. 7.

The power spectrum as a function of $k$-modes and redshift obtained from the
different accretion rate models are plotted in the middle and bottom panels of
Fig. 7. There are large differences of up to an order of magnitude between the
models, clearly exceeding the differences observed in the global signal alone.
The differences are smaller but still significant when looking at the models
with matching global signals due to renormalised flux parameters (see dashed
and dark-green solid lines). The power spectrum from these models can still
differ by about a factor of 3 at most. We therefore conclude that the mass
accretion rate has to be modelled at high accuracy in order to obtain reliable
predictions of the 21-cm power spectrum.

##  Appendix B Effects from halo bias and mass function prescriptions

In Appendix A we have shown that modelling choices regarding the mass
accretion rate may significant affect the 21-cm signal. Here we focus on
another central modelling component: the halo mass function and corresponding
bias prescription introduced in Eqs. (1-4).

Figure 8: Effects from the mass function and bias modelling on the 21-cm
signal, assuming source parameters from the fiducial model B. _Top-left and
centre:_ Halo bias and mass functions for different values of $q$ and $p$. The
purple and brown lines correspond to the Press-Schechter ($q=1$, $p=0$) and
Sheth-Tormen ($q=0.7$, $p=0.3$) mass functions, while the red line shows an
case in-between ($q=0.85$, $p=0.3$) which is in better agreement with high
redshift simulations. _Top-right:_ Global 21-cm signal resulting from these
mass function and bias prescriptions, where the flux parameters ($N_{\alpha}$,
$f_{X}$) are either kept the same (solid lines) or where they are changed so
that the absorption troughs lie at the same redshift (dashed lines). _Centre
and Bottom:_ Resulting power spectra as a function of $k$-modes and
redshsifts.

In order to test the sensitivity of the 21-cm global signal and power
spectrum, we vary the free model parameters $q$ and $p$ of the first crossing
distribution $f(\nu)$ of Eq. (3). We thereby investigate three main cases, the
original Press and Schechter (1974) model ($q=1$, $p=0$), the Sheth and Tormen
(2002) model ($q=0.707$, $p=0.3$), and an intermediate model ($q=0.85$,
$p=0.3$) that is in best agreement with results from high-redshift simulations
(see discussion in Sec. II.1).

The corresponding halo bias and mass functions of the three models are shown
in Fig. 8 (top-left and top-centre panels). Note that the difference between
the models become quite large in both cases, especially towards large halo
masses. As a result, the global 21-cm signal is shifted substantially, the
maximum of the absorption trough moving from $z\sim 16$ to $z\sim 18$ for the
case of the ST instead of the PS prescription. In order to counteract this
effect at the level of the global signal, we show two more models where the
flux parameters $N_{\alpha}$ and $f_{X}$ are modified in order to bring the
absorption trough back in line with the intermediate model (dashed lines).
Note, however, that there are remaining differences in the shape of the
absorption signal that cannot be absorbed by such a simple recalibration.

In the remaining panels of Fig. 8, we show the 21-cm power spectrum as a
function of $k$-modes and redshift. There are significant differences,
especially between the ST model (brown line) and the other two cases (solid
red and purple lines). The recalibrated models, on the other hand, are much
closer together, but differences of about a factor of two remain.

We conclude that next to the halo accretion rate, the halo bias and mass
functions need to be known to good accuracy in order to avoid large modelling
errors with the halo model approach. Note, that similar uncertainties also
exist for semi-numerical methods or even simulation results that may only be
trusted if they have realistic halo abundance and distributions.

## References

  * Haarlem  _et al._ (2006) M. V. Haarlem _et al._ , Electronics 12, 54 (2006).
  * Mertens  _et al._ (2020) F. Mertens _et al._ , Mon. Not. Roy. Astron. Soc. 493, 1662 (2020), arXiv:2002.07196 [astro-ph.CO] .
  * Tingay  _et al._ (2013) S. Tingay _et al._ , Publications of the Astronomical Society of Australia 30, 1 (2013).
  * Trott  _et al._ (2020) C. M. Trott _et al._ , Mon. Not. Roy. Astron. Soc. 493, 4711 (2020), arXiv:2002.02575 [astro-ph.CO] .
  * DeBoer  _et al._ (2017) D. R. DeBoer _et al._ , Publ. Astron. Soc. Pac. 129, 045001 (2017), arXiv:1606.07473 [astro-ph.IM] .
  * Koopmans  _et al._ (2015) L. Koopmans _et al._ , PoS AASKA14, 001 (2015), arXiv:1505.07568 [astro-ph.CO] .
  * Madau  _et al._ (1997) P. Madau, A. Meiksin, and M. J. Rees, Astrophys. J. 475, 429 (1997), arXiv:astro-ph/9608010 .
  * Shaver  _et al._ (1999) P. Shaver, R. Windhorst, P. Madau, and A. de Bruyn, Astron. Astrophys. 345, 380 (1999), arXiv:astro-ph/9901320 .
  * Ewall-Wice  _et al._ (2018) A. Ewall-Wice, T.-C. Chang, J. Lazio, O. Dore, M. Seiffert, and R. Monsalve, Astrophys. J. 868, 63 (2018), arXiv:1803.01815 [astro-ph.CO] .
  * Mirocha  _et al._ (2019) J. Mirocha _et al._ , (2019), arXiv:1903.06218 [astro-ph.CO] .
  * Park  _et al._ (2019) J. Park, A. Mesinger, B. Greig, and N. Gillet, Mon. Not. Roy. Astron. Soc. 484, 933 (2019), arXiv:1809.08995 [astro-ph.GA] .
  * McQuinn  _et al._ (2006) M. McQuinn, O. Zahn, M. Zaldarriaga, L. Hernquist, and S. R. Furlanetto, Astrophys. J. 653, 815 (2006), arXiv:astro-ph/0512263 .
  * Liu  _et al._ (2019) A. Liu _et al._ , (2019), arXiv:1903.06240 [astro-ph.CO] .
  * Barkana and Loeb (2005a) R. Barkana and A. Loeb, Astrophys. J. 626, 1 (2005a), arXiv:astro-ph/0410129 [astro-ph] .
  * Greig  _et al._ (2020) B. Greig _et al._ , (2020), arXiv:2006.03203 [astro-ph.CO] .
  * Liu and Slatyer (2018) H. Liu and T. R. Slatyer, Phys. Rev. D 98, 023501 (2018), arXiv:1803.09739 [astro-ph.CO] .
  * Mena  _et al._ (2019) O. Mena, S. Palomares-Ruiz, P. Villanueva-Domingo, and S. J. Witte, Phys. Rev. D 100, 043540 (2019), arXiv:1906.07735 [astro-ph.CO] .
  * Sitwell  _et al._ (2014) M. Sitwell, A. Mesinger, Y.-Z. Ma, and K. Sigurdson, Mon. Not. Roy. Astron. Soc. 438, 2664 (2014), arXiv:1310.0029 [astro-ph.CO] .
  * Lopez-Honorez  _et al._ (2016) L. Lopez-Honorez, O. Mena, A. Moliné, S. Palomares-Ruiz, and A. C. Vincent, JCAP 08, 004 (2016), arXiv:1603.06795 [astro-ph.CO] .
  * Schneider (2018) A. Schneider, Phys. Rev. D 98, 063021 (2018), arXiv:1805.00021 [astro-ph.CO] .
  * Lidz and Hui (2018) A. Lidz and L. Hui, Phys. Rev. D 98, 023011 (2018), arXiv:1805.01253 [astro-ph.CO] .
  * Nebrin  _et al._ (2019) O. Nebrin, R. Ghara, and G. Mellema, JCAP 04, 051 (2019), arXiv:1812.09760 [astro-ph.CO] .
  * Hill and Baxter (2018) J. C. Hill and E. J. Baxter, JCAP 08, 037 (2018), arXiv:1803.07555 [astro-ph.CO] .
  * Yang  _et al._ (2019) W. Yang, S. Pan, S. Vagnozzi, E. Di Valentino, D. F. Mota, and S. Capozziello, JCAP 11, 044 (2019), arXiv:1907.05344 [astro-ph.CO] .
  * Lopez-Honorez  _et al._ (2020) L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz, P. Villanueva-Domingo, and S. J. Witte, JCAP 06, 026 (2020), arXiv:2004.00013 [astro-ph.CO] .
  * Furlanetto  _et al._ (2004) S. Furlanetto, M. Zaldarriaga, and L. Hernquist, Astrophys. J. 613, 1 (2004), arXiv:astro-ph/0403697 .
  * Pritchard and Furlanetto (2007) J. R. Pritchard and S. R. Furlanetto, Mon. Not. Roy. Astron. Soc. 376, 1680 (2007), arXiv:astro-ph/0607234 .
  * Raste and Sethi (2018) J. Raste and S. Sethi, Astrophys. J. 860, 55 (2018), arXiv:1711.03827 [astro-ph.CO] .
  * Santos  _et al._ (2008) M. G. Santos, A. Amblard, J. Pritchard, H. Trac, R. Cen, and A. Cooray, Astrophys. J. 689, 1 (2008), arXiv:0708.2424 [astro-ph] .
  * Mesinger and Furlanetto (2007) A. Mesinger and S. Furlanetto, Astrophys. J. 669, 663 (2007), arXiv:0704.0946 [astro-ph] .
  * Mesinger  _et al._ (2010) A. Mesinger, S. Furlanetto, and R. Cen, Monthly Notices of the Royal Astronomical Society 411, 955 (2010).
  * Santos  _et al._ (2009) M. Santos, L. Ferramacho, M. Silva, A. Amblard, A. C. C. . Ist, and U. Irvine, Monthly Notices of the Royal Astronomical Society 406, 2421 (2009).
  * Fialkov  _et al._ (2014) A. Fialkov, R. Barkana, A. Pinhas, and E. Visbal, Mon. Not. Roy. Astron. Soc. 437, 36 (2014), arXiv:1306.2354 [astro-ph.CO] .
  * Mutch  _et al._ (2016) S. Mutch, P. M. Geil, G. Poole, P. Angel, A. Duffy, A. Mesinger, and J. Wyithe, Monthly Notices of the Royal Astronomical Society 462, 250 (2016).
  * Hutter (2018) A. Hutter, Mon. Not. Roy. Astron. Soc. 477, 1549 (2018), arXiv:1803.00088 [astro-ph.CO] .
  * Thomas  _et al._ (2009) R. M. Thomas _et al._ , Mon. Not. Roy. Astron. Soc. 393, 32 (2009), arXiv:0809.1326 [astro-ph] .
  * Ghara  _et al._ (2015) R. Ghara, T. R. Choudhury, and K. K. Datta, Mon. Not. Roy. Astron. Soc. 447, 1806 (2015), arXiv:1406.4157 [astro-ph.CO] .
  * Mellema  _et al._ (2006) G. Mellema, I. Iliev, M. Alvarez, and P. Shapiro, New Astron. 11, 374 (2006), arXiv:astro-ph/0508416 .
  * Friedrich  _et al._ (2012) M. Friedrich, G. Mellema, I. Iliev, and P. Shapiro, Monthly Notices of the Royal Astronomical Society 421, 2232 (2012).
  * Graziani  _et al._ (2013) L. Graziani, A. Maselli, and B. Ciardi, Monthly Notices of the Royal Astronomical Society 431, 722 (2013).
  * Ghara  _et al._ (2018) R. Ghara, G. Mellema, S. K. Giri, T. R. Choudhury, K. K. Datta, and S. Majumdar, Mon. Not. Roy. Astron. Soc. 476, 1741 (2018), arXiv:1710.09397 [astro-ph.CO] .
  * Ross  _et al._ (2019) H. E. Ross, K. L. Dixon, R. Ghara, I. T. Iliev, and G. Mellema, Mon. Not. Roy. Astron. Soc. 487, 1101 (2019), arXiv:1808.03287 [astro-ph.CO] .
  * Eide  _et al._ (2020) M. B. Eide, B. Ciardi, L. Graziani, P. Busch, Y. Feng, and T. Di Matteo, Mon. Not. Roy. Astron. Soc. 498, 6083 (2020), arXiv:2009.06631 [astro-ph.CO] .
  * Ocvirk  _et al._ (2016) P. Ocvirk _et al._ , Mon. Not. Roy. Astron. Soc. 463, 1462 (2016), arXiv:1511.00011 [astro-ph.GA] .
  * Doussot  _et al._ (2019) A. Doussot, H. Trac, and R. Cen, Astrophys. J. 870, 18 (2019), arXiv:1712.04464 [astro-ph.CO] .
  * Ocvirk  _et al._ (2020) P. Ocvirk _et al._ , Mon. Not. Roy. Astron. Soc. 496, 4087 (2020), arXiv:1811.11192 [astro-ph.GA] .
  * Rosdahl  _et al._ (2018) J. Rosdahl, H. Katz, J. Blaizot, T. Kimm, L. Michel-Dansac, T. Garel, M. Haehnelt, P. Ocvirk, and R. Teyssier, Monthly Notices of the Royal Astronomical Society 479, 994 (2018).
  * Katz  _et al._ (2020) H. Katz _et al._ , Mon. Not. Roy. Astron. Soc. 494, 2200 (2020), arXiv:1905.11414 [astro-ph.GA] .
  * Semelin  _et al._ (2017) B. Semelin, E. Eames, F. Bolgar, and M. Caillat, Mon. Not. Roy. Astron. Soc. 472, 4508 (2017), arXiv:1707.02073 [astro-ph.CO] .
  * Holzbauer and Furlanetto (2011) L. N. Holzbauer and S. Furlanetto, Monthly Notices of the Royal Astronomical Society 419, 718 (2011).
  * Press and Schechter (1974) W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
  * Bond  _et al._ (1991) J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Astrophys. J. 379, 440 (1991).
  * Sheth and Tormen (2002) R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 329, 61 (2002), arXiv:astro-ph/0105113 [astro-ph] .
  * Sheth and Tormen (1999) R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 308, 119 (1999), arXiv:astro-ph/9901122 .
  * Cooray and Sheth (2002) A. Cooray and R. K. Sheth, Phys. Rept. 372, 1 (2002), arXiv:astro-ph/0206508 .
  * Iliev  _et al._ (2011) I. Iliev, G. Mellema, P. Shapiro, U.-L. Pen, Y. Mao, J. Koda, and K. Ahn, Monthly Notices of the Royal Astronomical Society 423, 2222 (2011).
  * Behroozi  _et al._ (2020) P. Behroozi _et al._ ,  (2020), 10.1093/mnras/staa3164, arXiv:2007.04988 [astro-ph.GA] .
  * McBride  _et al._ (2009) J. McBride, O. Fakhouri, and C. Ma, Monthly Notices of the Royal Astronomical Society 398, 1858 (2009).
  * Dekel  _et al._ (2013) A. Dekel, A. Zolotov, D. Tweed, M. Cacciato, D. Ceverino, and J. Primack, Mon. Not. Roy. Astron. Soc. 435, 999 (2013), arXiv:1303.3009 [astro-ph.CO] .
  * Trac  _et al._ (2015) H. Trac, R. Cen, and P. Mansfield, Astrophys. J. 813, 54 (2015), arXiv:1507.02685 [astro-ph.CO] .
  * Neistein and van den Bosch (2006) E. Neistein and F. C. van den Bosch, Mon. Not. Roy. Astron. Soc. 372, 933 (2006), arXiv:astro-ph/0605045 .
  * Correa  _et al._ (2015) C. A. Correa, J. S. B. Wyithe, J. Schaye, and A. R. Duffy, Mon. Not. Roy. Astron. Soc. 450, 1514 (2015), arXiv:1409.5228 [astro-ph.GA] .
  * Vale and Ostriker (2004) A. Vale and J. Ostriker, Mon. Not. Roy. Astron. Soc. 353, 189 (2004), arXiv:astro-ph/0402500 .
  * Furlanetto  _et al._ (2017) S. R. Furlanetto, J. Mirocha, R. H. Mebane, and G. Sun, Mon. Not. Roy. Astron. Soc. 472, 1576 (2017), arXiv:1611.01169 [astro-ph.GA] .
  * Mebane  _et al._ (2018) R. H. Mebane, J. Mirocha, and S. R. Furlanetto, Mon. Not. Roy. Astron. Soc. 479, 4544 (2018), arXiv:1710.02528 [astro-ph.GA] .
  * Mason  _et al._ (2015) C. Mason, M. Trenti, and T. Treu, Astrophys. J. 813, 21 (2015), arXiv:1508.01204 [astro-ph.GA] .
  * Sun and Furlanetto (2016) G. Sun and S. R. Furlanetto, Monthly Notices of the Royal Astronomical Society 460, 417 (2016).
  * Mirocha  _et al._ (2017) J. Mirocha, S. R. Furlanetto, and G. Sun, Mon. Not. Roy. Astron. Soc. 464, 1365 (2017), arXiv:1607.00386 .
  * Pritchard and Loeb (2012) J. R. Pritchard and A. Loeb, Rept. Prog. Phys. 75, 086901 (2012), arXiv:1109.6012 [astro-ph.CO] .
  * Mirocha (2020) J. Mirocha,  (2020), 10.1093/mnras/staa3150, arXiv:2008.04322 [astro-ph.GA] .
  * Pritchard and Furlanetto (2006) J. R. Pritchard and S. R. Furlanetto, Mon. Not. Roy. Astron. Soc. 367, 1057 (2006), arXiv:astro-ph/0508381 [astro-ph] .
  * Barkana and Loeb (2004) R. Barkana and A. Loeb, Astrophys. J. 609, 474 (2004), arXiv:astro-ph/0310338 .
  * Mirocha (2014) J. Mirocha, Mon. Not. Roy. Astron. Soc. 443, 1211 (2014), arXiv:1406.4120 [astro-ph.GA] .
  * Furlanetto (2006) S. Furlanetto, Mon. Not. Roy. Astron. Soc. 371, 867 (2006), arXiv:astro-ph/0604040 .
  * Furlanetto and Stoever (2010) S. Furlanetto and S. J. Stoever, Monthly Notices of the Royal Astronomical Society 404, 1869 (2010).
  * Wouthuysen (1952) S. A. Wouthuysen, AJ 57, 31 (1952).
  * Field (1958) G. B. Field, Proceedings of the IRE 46, 240 (1958).
  * Furlanetto  _et al._ (2006) S. Furlanetto, S. Oh, and F. Briggs, Phys. Rept. 433, 181 (2006), arXiv:astro-ph/0608032 .
  * Mirocha  _et al._ (2015) J. Mirocha, G. J. A. Harker, and J. O. Burns, Astrophys. J. 813, 11 (2015), arXiv:1509.07868 [astro-ph.CO] .
  * Gilfanov  _et al._ (2004) M. Gilfanov, H. Grimm, and R. Sunyaev, Mon. Not. Roy. Astron. Soc. 347, L57 (2004), arXiv:astro-ph/0301331 .
  * Bharadwaj and Ali (2004) S. Bharadwaj and S. S. Ali, Mon. Not. Roy. Astron. Soc. 352, 142 (2004), arXiv:astro-ph/0401206 .
  * Barkana and Loeb (2005b) R. Barkana and A. Loeb, Astrophys. J. Lett. 624, L65 (2005b), arXiv:astro-ph/0409572 .
  * Note (1) Note that although our parametrisation is similar to that of Mirocha  _et al._ (2017), the global signal is shifted to slightly higher redshift due to differing assumptions about the star-formation efficiency normalisation and the minimum halo mass.
  * Murray  _et al._ (2020) S. G. Murray, B. Greig, A. Mesinger, J. B. Muñoz, Y. Qin, J. Park, and C. A. Watkinson,  (2020), 10.21105/joss.02582, arXiv:2010.15121 [astro-ph.IM] .
  * Zel’Dovich (1970) Y. B. Zel’Dovich, Astron. Astrophys. 500, 13 (1970).