# HERA Bound on x-ray luminosity when accounting for population III stars

Hovav Lazare<EMAIL_ADDRESS>Department of Physics, Ben-Gurion
University of the Negev, Be’er Sheva 84105, Israel  Debanjan Sarkar
<EMAIL_ADDRESS>Department of Physics, Ben-Gurion University of the
Negev, Be’er Sheva 84105, Israel  Ely D. Kovetz<EMAIL_ADDRESS>Department
of Physics, Ben-Gurion University of the Negev, Be’er Sheva 84105, Israel

###### Abstract

Recent upper bounds from the Hydrogen Epoch of Reionization Array (HERA) on
the cosmological 21-cm power spectrum at redshifts $z\approx 8,10$, have been
used to constrain $L_{\rm X<2\,keV}/{\rm SFR}$, the soft-band X-ray luminosity
measured per unit star formation rate (SFR), strongly disfavoring values lower
than $\approx 10^{39.5}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm
yr}$. This conclusion is derived from semi-numerical models of the 21-cm
signal, specifically focusing on contributions from atomic cooling galaxies
that host PopII stars. In this work, we first reproduce the bounds on $L_{\rm
X<2\,keV}/{\rm SFR}$ and other parameters using a pipeline that combines
machine learning emulators for the power spectra and the intergalactic medium
characteristics, together with a standard Markov chain Monte Carlo parameter
fit. We then use this approach when including molecular cooling galaxies that
host PopIII stars in the cosmic dawn 21-cm signal, and show that lower values
of $L_{\rm X<2\,keV}/{\rm SFR}$ are hence no longer strongly disfavored. The
revised HERA bound does not require high-redshift X-ray sources to be
significantly more luminous than high-mass X-ray binaries observed at low
redshift.

##  I Introduction

The 21-cm signal carries immense potential to probe the cosmic dawn and epoch
of reionization, two pivotal phases in the early Universe [1, 2, 3, 4, 5]. It
arises from the hyperfine transition of neutral hydrogen atoms, and its
measurement enables us to peer into the Universe’s infancy [6]. By studying
the fluctuations in the 21-cm signal, we can unravel the processes that drove
the formation of the first stars and galaxies [7, 8, 9], as well as the
subsequent ionization of the intergalactic medium (IGM) [10, 11]. The 21-cm
signal is a powerful tool for understanding cosmic evolution and unlocking the
mysteries of our cosmic origins [2, 3, 12].

Several groundbreaking experiments have been designed to target the 21-cm
signal [13, 14, 15, 16, 17, 18, 19], including the Hydrogen Epoch of
Reionization Array (HERA)111https://reionization.org/ [20]. HERA has recently
released two data sets [21, 22], that were used to set an upper bound on the
21-cm fluctuation power spectrum from the epoch of reionization
($z\\!\sim\\!8$ and $10$). This upper limit [21, 22], together with other
observations, including the galaxy ultraviolet (UV) luminosity function [23,
24, 25, 26], the reionization optical depth [27] and the IGM neutral fraction
[28] was then used to place the first constraints on properties related to the
evolution of the first galaxies and the epoch of reionization [29].

A major finding of Ref. [29] has to do with the parameter $L_{\rm
X<2\,keV}/{\rm SFR}$ that quantifies the X-ray luminosity sourced by the
endpoints of the stars that drive cosmic dawn and is defined as the ratio of
the integrated soft-band ($<2$ keV) X-ray luminosity to the star formation
rate (SFR). The HERA upper limit [21, 22] disfavors $L_{\rm X<2\,keV}/{\rm
SFR}\lesssim 10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$
values with high significance. Taken at face value, this would imply that
galaxies at high redshifts were more X-ray luminous than metal-enriched
galaxies in the local Universe [30, 31].

Figure 1: Marginalized posteriors for $L_{\rm X<2\,keV}/{\rm SFR}$, with (blue
and without (orange) HERA and with (solid) and without (dashed) molecular
cooling halos (MCGs). With the inclusion of MCGs, the preference for strong
X-ray heating weakens.

Crucially, the analysis in Ref. [29] considered only Population II (or PopII)
stars [32] as the sources of heating and ionization. The first generation of
stars (PopIII) [34, 33] are believed to have formed inside metal-poor galaxies
(sometimes called molecular cooling galaxies or MCGs) hosted by mini halos
where gas is cooled by molecular line transitions [35, 36, 37, 38, 34, 39, 40,
41]. While we have very limited knowledge about the formation and evolution of
MCGs, a number of simulations that have modelled them have shown that they can
significantly prepone the onset of cosmic dawn [42, 43, 44, 45, 46, 47]. In
addition, MCGs contribute to the total X-ray and ionizing photon budget of the
Universe [48, 49]. Ignoring MCGs in any analysis could therefore lead to an
overestimation of the X-ray and ionizing efficiencies of the atomic cooling
galaxies (ACGs) that are believed to host mostly PopII stars [46]. This
concern has also been raised in Ref. [22].

In this paper, we reanalyze the HERA data including MCGs, using machine
learning to emulate the 21-cm signal [50, 51, 52, 53, 54, 55, 56] and other
global quantities. The implication for the $L_{\rm X<2\,keV}/{\rm SFR}$
constraint is summarized in Fig. 1.

Simulations play a crucial role in modelling the 21cm signal and interpreting
the observational data [57, 58, 59, 61, 62, 63, 64, 65, 66, 60, 67, 68, 69,
70, 71, 72, 73, 74]. Due to the complexity of the astrophysical processes
involved, simulations help us understand the underlying physics and generate
theoretical predictions for the 21-cm signal. However, simulations can be
computationally expensive, if we need to simulate these complex processes
accurately. Furthermore, in the traditional Bayesian parameter-inference
pipeline, where we use simulations in each step to compare the outcomes
against observations, we require a significant amount of time before the
process converges. Even with faster semi-numerical codes such as
21cmFAST222github.com/21cmfast/21cmFAST [75], the MCMC process with a minimal
number of parameters can take up to several weeks to converge. We, therefore,
need alternative methods.

Machine learning (ML) techniques have emerged as powerful tools for emulating
the 21-cm signal and circumventing the time-consuming nature of traditional
simulators [50, 51, 52, 53, 54, 55, 56]. The importance of machine learning in
this context lies in its ability to significantly accelerate the parameter
inference process. By training a machine learning model on a large set of pre-
computed simulations, it becomes possible to generate fast and accurate
emulators that can predict the 21-cm signal for a given set of cosmological
and astrophysical parameters. Moreover, machine learning emulators can be
trained on a diverse range of simulations, enabling the exploration of
different cosmological scenarios and astrophysical processes.

With the goal of revisiting the analysis in Ref. [29] to check the robustness
of the $L_{\rm X<2\,keV}/{\rm SFR}$ constraint to the inclusion of the PopIII
stars or MCGs contribution, we have built an ML-based emulator, that emulates
the 21-cm power spectrum $\Delta^{2}_{21}(k)$ as a function of wave vector $k$
at different redshifts ranging from cosmic dawn to the epoch of reionization.
We also built separate emulators for the CMB optical depth and IGM neutral
fraction.

Our paper is organized as follows. In Section II we briefly describe the
models used in the HERA analysis [29] and the results of their analysis. In
Section III we describe the HERA phase-I data and the likelihood we use for
our analysis, including the additional observables. In Section IV we describe
the artificial neural network architecture our emulator is based on, how we
build our training, validation and test sets, and demonstrate the performance
of emulator. In Section V we explain how we include the PopIII contribution in
the modeling of the 21-cm signal. We then present our results in Section VI.
After reproducing the HERA parameter constraints with our pipeline when
considering only PopII stars, we include the MCGs in our simulations and train
corresponding emulators. We then perform the MCMC analysis by considering both
the ACG and MCG parameters. In particular, we present the marginalized
posterior of $L_{\rm X<2\,keV}/{\rm SFR}$ from this analysis (summarized in
Fig. 1). We find that with the mini halos, we can no longer discard the
$L_{\rm X<2\,keV}/{\rm SFR}\lesssim 10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm
M}_{\odot}^{-1}\;{\rm yr}$ X-ray luminosity values easily. Although there is a
small decline in posterior probability at $L_{\rm X<2\,keV}/{\rm
SFR}<10^{39.5}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$, the
$L_{\rm X}/{\rm SFR}\sim 10^{39}\,{\rm erg}\;{\rm s}^{-1}\;{\rm
M}_{\odot}^{-1}\;{\rm yr}$ values are still allowed. This suggests that the
relationship between star formation and soft X-ray luminosity for high
redshift galaxies may not be very different from the local galaxies [31, 30].
We discuss and elaborate on this conclusion in Section VII.

##  II Summary of the HERA Analysis

Ref. [29] have considered four different models of reionization to interpret
the HERA observation. We briefly discuss the models below, and the inferences
based on those.

  * •

Density-driven linear bias model: In this model, the 21-cm fluctuations are
assumed to follow the density fluctuations and the 21-cm power spectrum is
proportional to the matter power spectrum. A bias parameter, which is a
function of quantities like $x_{\rm HI}$, $T_{\rm S}$ and $T_{\rm CMB}$,
relates the two power spectra.

  * •

Phenomenological reionization-driven model: In this model, IGM is considered
as a two phased system consisting of fully ionized bubbles and medium with
uniform temperature outside of the bubbles [76].

  * •

Semi-numerical model I (21cmFAST): This model assumes that the star-forming
galaxies reside inside dark matter haloes and these galaxies are responsible
for the heating and ionization of the IGM. This utilizes empirical scaling
relations to relate galaxy properties to their host dark matter haloes. The
galaxy properties include stellar to halo mass ratio, escape fraction of
ionizing UV radiation, the star formation rate and X-ray luminosity [75].

  * •

Semi-numerical model II: This is an alternate semi-numerical model [77], built
on similar assumptions as 21cmFAST. On top of that, this model allows for an
excess radio background. Excess radio background, on top of CMB, produces an
enhanced 21-cm absorption signal and can also increase the amplitude of the
21-cm power spectrum.

The first two models relate the 21-cm signal directly to the IGM properties,
without trying to model the source properties. On the other hand, the semi-
numerical models try to model the properties of the radiation sources and
their evolution in redshift.

The analysis of Ref. [29] based on all of the above models suggest that the
IGM at $z=10.4$ is heated over the adiabatic cooling limit at $>95\%$
confidence, and this rule out the “cold reionization” scenarios in which the
IGM continues to adiabatically cool until it reionizes.

Considering the Semi-numerical model I, Ref. [29] quotes a lower limit on the
parameter $L_{\rm X<2\,keV}/{\rm SFR}$ which describes the heating ability of
the EoR galaxies per unit of star formation. We discuss the significance of
$L_{\rm X<2\,keV}/{\rm SFR}$ later in this paper.

HERA observations, Based on the Semi-numerical model II, and combined with the
Chandra X-ray background constraints, rule out most of the models that explain
the radio background observed by LWA [78] and ARCADE-2 [79] as originating at
$z\gtrsim 8$.

Note that the inferences in Ref. [29] are very much dependent on the choice of
physical models and the interpretations of different quantities change with
models. In our current analysis, we mainly focus on Semi-numerical model I
(21cmFAST), and we discuss the findings based on this model in detail.

##  III Data and Likelihood

###  III.1 The 21-cm Observables

The emission or absorption of the 21-cm line is characterized by the spin
temperature, $T_{S}$, and is usually measured as the differential brightness
temperature, $\delta T_{21}$, with respect to the brightness temperature of
the low-frequency radio background, $T_{\rm{rad}}$. In the usual scenario,
$T_{\rm{rad}}$ is taken to be the CMB temperature
$T_{\rm{cmb}}=2.7254\times(1+z)$ K where $z$ refers to the cosmological
redshift. The differential brightness temperature, $\delta T_{21}$, can be
expressed as [80, 81, 82]

$\displaystyle\delta T_{21}(\nu)$ $\displaystyle=\frac{T_{\mathrm{S}}-T_{\rm
rad}}{1+z}\left(1-e^{-\tau_{\nu_{0}}}\right)\,.$ (1)

Here, the factor $\left(1-e^{-\tau_{\nu_{0}}}\right)$ exhibits the effect of
propagation through a medium, where $\tau_{\nu_{0}}$ is the optical depth.
Note that, in this equation, the astrophysical information is mostly encoded
in $T_{\mathrm{S}}$, whereas the state of the propagating medium is captured
by the optical depth $\tau_{\nu_{0}}$.

The measurement of $\delta T_{21}$ is done in two different ways:
(i) Global signal : Here, $\delta T_{21}(\mathbf{n},z)$ is first measured
along different directions $\mathbf{n}$ of the sky, keeping the redshift (or
frequency) fixed. The measurements are finally averaged over the directions to
convert into an average quantity $\langle\delta T_{21}(z)\rangle$, which is
known as the global signal [83, 84].
(ii) Fluctuations: Here, telescopes measure the spatial fluctuations of the
$\delta T_{21}(\mathbf{x},z)$ field. The fluctuations are then interpreted in
terms of statistics like the power spectrum [85, 86, 87, 76, 88, 89, 90],
bispectrum [91, 92, 93, 94, 95, 96, 97], etc. The 21-cm power spectrum is the
most relevant quantity for this work, and is defined as

$\langle\tilde{\delta T}_{21}({\bf k}_{1})\tilde{\delta T}_{21}({\bf
k}_{2})\rangle=(2\pi)^{3}\delta^{D}({\bf k}_{1}-{\bf k}_{2})P_{21}({\bf
k}_{1})\,,$ (2)

where $\langle...\rangle$ denotes an ensemble average, $\tilde{\delta
T}_{21}({\bf k})$ is the Fourier transform of $\delta T_{21}(\bf x)$ and
$\delta^{D}$ is the Dirac delta function. Specifically, the quantity of
interest for us is $\Delta_{21}^{2}({\bf k})\equiv k^{3}P_{21}({\bf
k})/(2\pi^{2})$, which will be the output of HERA measurements, and of our
simulations and emulators.

###  III.2 HERA Phase-I data

The upper limits on the 21-cm signal from HERA were first published in Ref.
[21], and improved limits were published in Ref. [22] These limits are based
on 94 nights of observation using the HERA Phase I experimental configuration.
For detailed information on the configuration, the reader is referred to Ref.
[22]. Out of all the different bands of observation, in this work we
concentrate on two bands: Band 2, centered at $z=7.9$; and Band 1, centered at
$z=10.4$. These bands are largely free of radio frequency interference (RFI).
Following Ref. [22], we use the observed power spectrum computed using all
fields for each redshift band. For a comprehensive reading about the data, the
reader is referred to Ref. [22, 21].

###  III.3 The data Likelihood

The form of the likelihood function is important for any Bayesian analysis.
The 21-cm power spectrum measurements from HERA suffer from (largely unknown)
systematics which need to be taken into account in the likelihood analysis.
Below, we briefly explain how to derive the likelihood in the presence of
systematics.

First, we divide our data into two vectors, one for every redshift bin,

$\textbf{d}_{1}=\left(\begin{array}[]{c}\Delta^{2}_{21}(k_{1},z_{1})\\\
\Delta^{2}_{21}(k_{2},z_{1})\\\
\vdots\end{array}\right),\textbf{d}_{2}=\left(\begin{array}[]{c}\Delta^{2}_{21}(k_{1}^{\prime},z_{2})\\\
\Delta^{2}_{21}(k_{2}^{\prime},z_{2})\\\ \vdots\end{array}\right),$ (3)

In general $k_{i}\neq k_{i}^{\prime}$. Considering a single data vector, in
order to reduce the covariance between the neighbouring $k$-bins, we follow
the method adopted in Ref. [29]. Starting from a minimum $k$-bin (which is
decided on the basis of the signal-to-noise ratio), we include the alternate
$k$-bins in our analysis such that the bins do not have significant overlap.
In the following derivation, we show the form of the likelihood function for a
single data vector, and the composite likelihood for the two data vectors may
be calculated as
$\mathcal{L}\propto\mathcal{L}_{1}(\mathbf{d}_{1})\cdot\mathcal{L}_{2}(\mathbf{d}_{2})$.

Following Ref. [29], we assume a Gaussian likelihood for the data, and denote
the unknown systematics as u, so that one can express the likelihood function
for a data vector d as

$\displaystyle\mathcal{L}(\mathbf{d}|\theta,\mathcal{M},\mathbf{u})\propto\exp\left(-\frac{1}{2}r(\theta,\mathbf{u})^{T}\Gamma
r(\theta,\mathbf{u})\right)$ (4)

where
$r(\theta,\mathbf{u})\equiv\textbf{d}-\textbf{W}\cdot\textbf{m}(\theta)-\textbf{u}$,
$\textbf{m}(\theta)$ is the emulator prediction for the data vector d, W is
the Window function (which accounts for the point spreading in Fourier k space
[98, 99, 100, 101, 102]), and $\Gamma\equiv\Sigma^{-1}$ where $\Sigma$ in the
covariance matrix. In the parameter estimation problem, we are interested in
finding the posterior probability for obtaining a parameter set given the
observation data. This can be expressed in terms of the likelihood, based on
Bayes’ theorem, as

$p(\theta|\textbf{d},\mathcal{M},\textbf{u})\propto\mathcal{L}(\textbf{d}|\theta,\mathcal{M},\textbf{u})p(\theta|\mathcal{M})p(\textbf{u}),$
(5)

where $p(\theta|\mathcal{M})$ is the prior function of the parameter set
$\theta$ given the model $\mathcal{M}$, and $p(\textbf{u})$ denotes the priors
on the systematics. In general, priors on systematics should come from the
experiments. However, in case of HERA, the systematics are largely unknown.
Therefore, we follow Ref. [29] and marginalize over the systematics prior
range in order to obtain the likelihood. We assume a multivariate uniform
distribution for u, where the systematic noise in each $k$ and $z$ bins are
independent. The marginalization can be done using,

$\displaystyle\mathcal{L}(\mathbf{d}|\theta,\mathcal{M},\mathbf{u})\propto\int_{\textbf{u}_{\rm
min}}^{\textbf{u}_{\rm
max}}\exp\left(-\frac{1}{2}\textbf{r}(\theta,\textbf{u})^{T}\Gamma\textbf{r}(\theta,\textbf{u})\right)d\textbf{u},\,.$
(6)

Assuming $\Gamma$ to be diagonal, the integral takes a simple form

$\displaystyle\mathcal{L}(\mathbf{d}|\theta,\mathcal{M},\mathbf{u})\propto\prod_{i}^{N_{d}}\int_{u_{i,{\rm
min}}}^{u_{i,{\rm
max}}}\exp\left(-\frac{\left[y_{i}-u_{i}\right]^{2}}{2\sigma_{i}^{2}}\right)du_{i},$
(7)

where the index ‘i’ denotes the $k$ bin number, $y_{i}=d_{i}-[Wm(\theta)]_{i}$
is the difference between the model and the data at each bin, $N_{d}$ is the
number of k bins, and $\sigma_{i}$ are the diagonal elements of the covariance
matrix. In order to calculate the integral, we set $\textbf{u}\in[0,\infty]$,
assuming systematics to be positive (although negative systematics can occur
in some cases [103] and those are subject to interpretation). Marginalization
over the systematics using the assumed bounds yields

$\displaystyle\mathcal{L}(\mathbf{d}|\theta,\mathcal{M})\propto\prod_{i}^{N_{d}}\frac{1}{2}\left(1+{\rm
erf}\left[\frac{y_{i}}{\sqrt{2}\sigma_{i}}\right]\right),$ (8)

where ‘$\rm erf$’ is the error function. We note that the marginalized
likelihood function is sensitive to the upper bounds of the data points, which
is particularly useful for the HERA power spectrum data.

#### Additional Observables

Note again that HERA measurements place only an upper bound on the 21-cm power
spectrum. Therefore, the Likelihood in Eq. (8) does not impose any tight
bounds on the model parameters. Following Ref. [29], we use additional
observables from other experiments to complement the HERA power spectrum data.
These additional observables are: (i) the upper bound on the neutral hydrogen
fraction $x_{HI}<0.06+0.05(1\sigma)$ at $z\sim 5.9$, measured by the dark
fraction on high redshift quasar spectra [28], (ii) the faint galaxy UV
luminosity functions at $z=6,7,8,10$ [23, 24, 26]
333https://github.com/21cmfast/21CMMC/tree/master/src/py21cmmc/data, and (iii)
the Thomson scattering optical depth of CMB photons
$\tau_{e}=0.0569^{+0.0081}_{-0.0066}$, using Planck [27] data from Ref. [104].
We shall see later that these additional observables help to place tight
bounds on some of the model parameters, while HERA helps constrain some of the
remaining parameters.

##  IV Artificial Neural Network (ANN)

In order to perform a Markov chain Monte Carlo (MCMC) analysis, we need to be
able to produce realisations of the 21-cm power spectra and other global
quantities, given a set of input parameters. However, to simulate the 21-cm
signal for a large cosmological volume with the required resolution using the
latest version of 21cmFAST [47] takes $\mathcal{O}(1\,{\rm hour})$. The
traditional MCMC pipeline requires to generate a large number of simulations
in a sequential manner to be able to construct the final posterior, which
takes a lot of time and can be prohibitive. In order to address this, we have
developed an artificial neural network (ANN)-based 21-cm signal emulator,
which we describe below.

###  IV.1 ANN Architecture

####  IV.1.1 Simulations

Parameters | Description
---|---
$f_{\star}$ | Stellar to halo mass ratio for haloes
$\alpha_{\star}$ | Stellar-to-halo mass power-law index
$f_{\rm esc}$ | Escape fraction of ionizing photons from the haloes
$\alpha_{{\rm esc}}$ | Escape fraction of ionizing photons to halo mass power-law indices
$M_{\rm turn}$ | Mass scale below which inefficient cooling and/or feedback suppresses efficient star formation
$t_{\star}$ | Star formation time-scale in units of $H^{-1}(z)$
$L_{\rm X<2\,keV}/{\rm SFR}$ | Soft-band X-ray luminosity per SFR in units of ${\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$ for ACGs
$E_{0}$ | Minimum X-ray energy (in eV) escaping the galaxies into the IGM
$\alpha_{{\rm X}}$ | Spectral index of X-ray sources
Table 1: The main astrophysics parameters.

It is relatively well accepted that the galaxies form inside the dark matter
halos, with different feedback processes shape their evolution [105, 106, 107,
108, 109, 110]. These galaxies are the main sources of various radiation
fields that affect the 21-cm signal from the beginning of cosmic dawn to the
epochs that follow. We use the semi-numerical code 21cmFAST that simulates
various galaxy properties, depending on the host halo mass, and produces these
radiation fields [45, 46, 47, 75, 111, 112, 113]. The code then uses the same
radiation fields to deduce the 21-cm fluctuations at different redshifts. Note
that the semi-numerical calculations of 21cmFAST are tuned to reproduce the
galaxy UV luminosity functions [114, 115, 116] observed at reionization as
well as the observed opacity of IGM.

Given a set of cosmological and astrophysical parameters, 21cmFAST generates
realizations of the 21-cm fluctuations in coeval and/or light-cone boxes.
Those boxes are then used to compute the global quantities like the mean
brightness temperature of the 21-cm signal $\langle\delta T_{21}\rangle$,
optical depth of reionization $\tau_{\rm reio}$, etc., and 21-cm fluctuation
power spectrum $\Delta^{2}_{21}(\mathbf{k})$. In addition to this, 21cmFAST is
also able to output the UV luminosity function of galaxies at various
redshifts. We create an ensemble of the observables (which we call data sets)
by changing the model parameters, and this ensemble will be used to train and
validate our emulator, as we describe later in this section. For the 21-cm
power spectrum, in this work we are mainly concerned about two redshift bins
$z=7.9$ and $10.4$ where HERA observations have the best detection [29].

The definitions of the important astrophysical parameters used in this work
are given in Table 1. The 21cmFAST code models the effects of both the PopIII
and PopII stars. PopIII [34, 33] refers to the first generation of stars
believed to reside in molecular cooling halos (MCGs) or mini-halos, while
PopII [32] refers to the second generation of stars which predominantly reside
in atomically cooled galaxies or ACGs. We shall discuss these in more detail
in Section V. The parameters related to the ACGs (MCGs) are represented by a
‘10’ (‘7’ or ‘mini’) in the subscript (referring to the typical halo mass). In
the next few sections, we shall assume that the cosmic dawn and reionization
processes are driven solely by the PopII stars. We build the data set in the
following section using 21cmFAST where we consider only PopII stars. We note
that the same modelling has been considered in Ref. [29].

####  IV.1.2 Building the Data set

In this section, we discuss the method of composing the data set which will
then be divided into training, validation and testing sets. We start by
determining the desired range over which we want to vary each parameter.
Having its range fixed, we sample the whole parameter space using a Latin
Hypercube (LH) sampler which seeks to produce uniform sampling densities when
all points are marginalized onto any one dimension [117]. After getting the
sampled parameter sets, we run the 21cmFAST code to obtain the 21-cm power
spectrum and global quantities like $\tau_{e}$, $x_{\rm HI}$ for each
combination. We then save these quantities along with the parameter values.

We now discuss the data set which consists of the 21-cm power spectrum as a
function of $k$. Note that, we have a limited number of $k$-bins from the
simulations to begin with. Before the training process, we use interpolation
to increase the number of power spectrum $k$-bins to 70, ranging from $k=0.1$
Mpc-1 to $k=1.73$ Mpc-1, so that the variation with the wave number is
smoother. The data set is then classified into 6 different classes (denoted by
index ‘n’) according to the amplitude difference between k bins. The first 4
$(n=0-3)$ classes are for signals where the difference between the maximal
amplitude and the minimal amplitude of the power spectrum, is smaller than
$10^{n}$ mK2 but greater than $10^{(n-1)}$ mK2. Class 5 is for signals where
the difference between the maximum and minimum power spectrum values in
different bins, i.e. [max(power spectrum) - min(power spectrum)] is $>10^{3}$
mK2, Finally, class 6 is reserved for signals where the minimal amplitude is
not in the lowest $k$ bin (these are exotic-shaped signals). This
classification is done so that when the emulator is trained, we can take
approximately the same number of samples from each class, so that the trained
emulator will not be biased towards any of the classes. Without doing so,
emulating exotic signals tends to be more difficult, since they are rare in
the un-manipulated data set. Note also that the ‘class 0’ signals are also
very hard to emulate even after this process because most of them are actually
‘dead’ or flat signals having no features that the emulator can learn. This
happens when for some parameter combination, the fluctuations at all scales
cease at a specific redshift. For example, when the IGM becomes sufficiently
ionized at low redshifts that the fluctuations in the electron density field
become more or less uniform, and the power spectrum is flat (and even becomes
zero in some cases). Emulating the flat signals is extremely difficult, and we
have decided not to include these at all in the emulating process. Including
flat signals in the emulating process leads to unreliable predictions. A
method for predicting these flat signals (for use in the MCMC pipeline) will
be explained later. The bottom line is that we have discarded the class 0
signals when training the emulator.

For the training process, we take $\sim$ 2500 randomly sampled signals from
each class for $z=7.9$, and $\sim$ 3000 signals from each class for the
$z=10.4$ band. We choose the numbers such that we have a similar number of
samples in each class, and such that no class is over/under-represented.
Finally, the data set is divided into training (85%), validation (10%) and
testing (5%) sets. Before the training process, the astrophysical parameters
are normalised to be in the range $[-1,1]$. We do not normalize the power
spectrum as our initial tests suggested a better performance with the original
power spectrum.

####  IV.1.3 NN Architecture

In this section, we discuss the architecture of our NN emulators. All of our
NN emulators and classifiers are built using the TensorFlow [118] and Keras
[119] libraries defined in Python3. Here we consider a 9-parameter emulator
model based on 21cmFAST. According to this model, only the halos which are
cooled by atomic line transitions host sources that start cosmic dawn and
cosmic reionization. The choice of the parameters is based on Ref. [29]. As
explained above, there is another possible cooling mechanism, which makes it
possible to form sources inside mini-halos. Using 21cmFAST simulations, it is
possible to model these mini halos, and we shall discuss this in more detail
in the later sections.

Therefore, the emulator has 9 input parameters, 70 output $k$ bins, and 6
fully connected hidden layers \- of 288, 512, 512, 288, 512, 1024 dimensions,
respectively. The number of layers and the number of neurons in each was
chosen after an optimization process. We have tested the emulator for a number
of conventional activation functions, namely ReLU [120] and $tanh$, etc. We
finally settled with the following function (taken from Ref. [121])

$f(\textbf{x},\bm{\alpha})=\textbf{x}\odot\left(\bm{\alpha}+(1-\bm{\alpha})\odot\left(1+e^{-\textbf{x}}\right)^{-1}\right)\,,$
(9)

where x are the layer inputs, $\bm{\alpha}$ is a vector of trainable
parameters, and $\odot$ denotes the element-wise product. The activation is
followed by a batch normalization layer [122], which standardises the input
for the next layer. After the forward propagation, the loss is calculated
using a relative percentage error

$\mathcal{L}(y_{true},y_{pred}(\theta))=\frac{|y_{true}-y_{pred}(\theta)|}{y_{true}}\times
100\,.$ (10)

The training process is performed using mini batches of size 128 samples, so
that the trainable parameters are updated using the gradients of the loss
function, after forward propagating each batch. The training process was done
using the Adam [123] optimizer, where the initial learning rate is 0.01, and
is reduced when the prediction on the validation set does not improve over
more than 10 epochs. The training phase ends when there is no improvement over
more than 30 epochs, and the best weights are stored at that moment. Finally,
we examine our emulator performance on the test set, and the results are
presented in Fig. 2.

Figure 2: Emulator loss statistics on the testing sets, of size $\sim$ 650
signals for each band. The box represents all the samples that fall between
the 25th and the 75th percentile, and the line in the middle of it is the
median. The black circles are outliers that extend beyond the 95th percentile.

####  IV.1.4 Retraining

After training the emulator, we review its performance very carefully on the
test set. In the context of the power spectrum, we find that on rare occasions
(less than 5%), the emulator prediction differs from the real value by almost
100%. This may bias the MCMC results, particularly the median value. We do the
following to make sure that this error does not affect our results.

After the first training, we determine a rough $\sim 1\;\sigma$ region in the
parameter space by running the MCMC analysis using the emulator. We then
generate $\sim 5000$ more signals by sampling the rough $1\;\sigma$ parameter
space using the LH method, where the $1\;\sigma$ is considered for all the
parameters except $L_{\rm X<2\,keV}/{\rm SFR}$, for which the emulator has to
be accurate over the full prior range. We then re-train our emulator on this
re-sampled data set.

Note that, for some astrophysical parameters such as $\alpha_{\ast}$, the
$1\;\sigma$ constraint range is very small as we shall show later in Fig. 7.
For those parameters, we choose a wider range when sampling. We show the
performance results after the re-training process in Figs. 3 and 4. We see
that the re-trained emulator gives a median error of 2-3% in the $1\;\sigma$
parameter space. Now considering the whole range of parameters that we used
for initial training, we see that our re-trained emulator exhibits $<20\%$
median error outside of the $1\;\sigma$ parameter space. This suggests that
the emulator will have better precision, at least within the $1\;\sigma$
parameter space, such that it will not bias the MCMC analysis by much and will
provide better $1\;\sigma$ constraints.

We shall see later (see also Ref. [29]) that most of the astrophysical
parameters, courtesy of the additional observables used in addition to the
HERA power spectrum data, are very well constrained. Other parameters, like
$L_{\rm X<2\,keV}/{\rm SFR}$ (which is meaningfully constrained only at the
lower side) or $E_{0}$, have almost no constraints from the observations. Now,
due to the property of the Markov chains, the MCMC algorithm will mostly
sample the region of the parameter space where our retrained emulator performs
extremely well. This means we shall have many accurate MCMC samples in the
vicinity of the $1\;\sigma$ range, and very few less-accurate samples outside
this region. This weighs down any possible bias and ensures our predictions
are robust.

Figure 3: Retrained emulator loss statistics on our testing sets, of size
$\sim$ 250 signals for each band. The box represents all the samples that are
between the 25th and the 75th percentiles, and the line in the middle of it is
the median. The black circles are outliers that extend beyond the 95th
percentile. Figure 4: Different view of the loss statistics for the retrained
emulator. The dashed lines represent the 16th and 84th percentiles
respectively, and their exact value is written along with the median value in
the figure legend.

####  IV.1.5 Classification

We have already mentioned in Section IV.1.2 that the class 0 signals or flat
power spectra are very hard to emulate. The power spectrum amplitude in this
class is close to zero, $\lesssim 3$ mK2 for all the k bins. When we run the
MCMC process, for each set of parameters we need to know whether the
parameters will produce flat signals even before we feed those to the ANN. We,
therefore, need machinery that predicts the occurrence of the flat signal for
a given parameter set.

We have built a neural network classifier for the above purpose. This
classifier is trained at each frequency band to predict the occurrence of the
flat 21-cm power spectrum. We discuss the architecture of the classifier in
Appendix C. We find that the classifier’s prediction accuracy is $\gtrsim
97\%$. This indicates that, on average, we shall have a miss-classification
less than once in every 30 MCMC iterations. This does not boast a serious
error in the MCMC analysis.

Considering the flat signal, we choose the following path. When the classifier
predicts a flat signal for some parameter combination, we generate a random
power spectrum at each $k$-bin where the amplitude is drawn from a Gaussian
distribution having a mean of 2 mK2 and a standard deviation of 1 mK2. This we
do irrespective of the parameter values. Our goal here is not to predict the
flat signal accurately as these signals are still somewhat rare. Also,
considering the upper bound on the 21-cm power spectrum from HERA (and our
chosen likelihood function in Eq. (8)) it is clear that any flat signal below
the upper limit is allowed. Therefore, HERA data does not really discard these
parameter values, rather the constraints come from the additional likelihoods.
Therefore, we surmise that our solution for tackling the flat power spectra
does not pose any serious issues.

###  IV.2 Emulating the other observables

We extend our NN for emulating $x_{HI}$ and $\tau_{e}$. These are easy to
emulate as for a fixed redshift both of them are just single numbers. We have
therefore used a relatively simple neural network architecture for emulating
both $x_{H_{I}}$ and $\tau_{e}$. As we discuss in Appendix A, we find that the
emulators performs extremely well and the emulation error is negligible. As
for the simulated UV luminosity functions, we do not have any emulators for
these. UV luminosity functions can be calculated analytically very quickly via
21cmFAST, and in the MCMC analysis we generate these directly by plugging in
the parameters.

##  V The role of Mini-Halos

The current understanding of the stellar composition in galaxies at low
redshifts is quite comprehensive. However, our knowledge about the initial
galaxies that sparked cosmic reionization is limited. According to the
hierarchical model of structure formation, the first generation of stars,
known as PopIII stars, likely formed inside small molecular cooling galaxies
or MCGs around $z\sim 20-30$ [35, 36, 37, 38, 34, 39, 40, 41]. Inside MCGs,
the gas cools mainly through the ${\rm H_{2}}$ rotational–vibrational
transitions efficient at $T_{\rm vir}\sim 10^{3}-10^{4}\,{\rm K}$. These early
galaxies were believed to be hosted by halos with virial temperatures
$\lesssim 10^{4}$ K, corresponding to total halo masses of $\lesssim
10^{8}\,{\rm M}_{\odot}$ during the epochs of reionization and cosmic dawn
[124, 44, 41, 45, 47]. These small halos are often called mini-halos. Over
time, feedback mechanisms suppressed star formation in these MCGs, leading to
the emergence of heavier, atomically cooled galaxies (ACGs) with halo masses
$\gtrsim 10^{8}\,{\rm M}_{\odot}$. Therefore, most of the ACGs at high
redshifts are “second-generation” galaxies, forming out of MCGs. It is
believed that the second generation of stars, known as PopII stars,
predominantly reside in ACGs [46].

PopII stars likely drove most of the heating and ionization process at $z\sim
5-10$ [125, 126, 127]. However, there should be a small but significant
contribution from the PopIII stars or MCGs as well. The star formation in MCGs
is subject to a number of feedback processes like the Lyman-Werner radiation
feedback [74, 128, 129, 130, 131, 36, 132], baryon supersonic streaming
velocities [133, 134, 135, 136, 137, 138, 139], etc. Depending on these, the
radiation contribution from the MCGs varies. MCGs spark the cosmic dawn epoch
earlier compared to the scenarios where only ACGs exist [74]. X-ray photons
from MCGs also prepone the epoch of heating [74], although MCGs are believed
to have a negligible contribution to reionization. As emphasized above, in
Ref. [29], the 21-cm signal was modelled without taking the contribution of
the PopIII stars or MCGs into account. One of the main conclusions of Ref.
[29] is that the early galaxies were more X-ray luminous than their present-
day counterparts. Our main goal in this paper is to check the validity of this
result in the presence of a contribution to the 21-cm signal from PopIII
stars.

One of the most salient features in the latest version of 21cmFAST code [47]
is that it now includes the contributions of both the MCGs and ACGs. The code
assumes MCGs host PopIII stars with a simple stellar-to-halo mass relation
(SHMR), distinct from that of ACGs hosting PopII stars. For a comprehensive
reading about how 21cmFAST incorporates PopII and PopIII stars in the
calculations, the reader is referred to Ref. [47, 45].

As mentioned earlier, star formation in mini haloes depends on feedback
processes which further decide the molecular-cooling turnover halo mass.
Haloes below this mass are not capable of forming stars. In 21cmFAST, this
turnover mass is a product of three factors [47],

$M_{\rm mol}=M_{0}(z)f_{v_{\rm cb}}(v_{\rm cb})f_{\rm LW}(J_{21})\,,$ (11)

where $M_{0}(z)$ is the molecular cooling threshold in the absence of feedback
with $M_{0}(z)=\tilde{M_{0}}\times(1+z)^{-3/2}$ where $\tilde{M_{0}}=3.3\times
10^{7}\,{\rm M}_{\odot}$ (corresponding to a viral temperature of $10^{3}$ K
[37]), $f_{v_{\rm cb}}(v_{\rm cb})$ and $f_{\rm LW}(J_{21})$ encodes the
effects of respectively the baryon streaming velocity and Lyman-Werner
background radiation. $f_{\rm LW}(J_{21})$ is parametrized as,

$f_{\rm LW}(J_{21})=1+A_{\rm LW}(J_{21})^{\beta_{\rm LW}}\,,$ (12)

with $A_{\rm LW}$ and $\beta_{\rm LW}$ are two free parameters that quantifies
the strength of the feedback, and $J_{21}$ being the Lyman-Werner radiation
intensity. Similarly, $f_{v_{\rm cb}}(v_{\rm cb})$ is parametrized as,

$f_{v_{\rm cb}}(v_{\rm cb})=\left(1+A_{v_{\rm cb}}\frac{v_{\rm cb}}{v_{\rm
rms}}\right)^{\beta_{v_{\rm cb}}}\,,$ (13)

where $v_{\rm cb}$ is the local value of the streaming velocity, $v_{\rm rms}$
is the rms of the streaming velocity, $A_{v_{\rm cb}}$ and $\beta_{v_{\rm
cb}}$ are two free parameters that determine the strength of the feedback. We
have kept the feedback parameter values fixed to $A_{\rm LW}=2.0$, $\beta_{\rm
LW}=0.6$, $A_{v_{\rm cb}}=1$ and $\beta_{v_{\rm cb}}=1.8$ throughout all of
the simulations that include MCGs. Note that the feedback, or in other words,
the parameters, are largely unknown at the redshifts that we are interested
in. The contribution of MCGs in cosmic dawn and reionization vary with the
variation in these parameters. We do not vary these in our analysis as these
would increase the number of parameters. However, we expect that our results
will still be valid even when these parameters are varied in the vicinity of
the fiducial values.

The 21cmFAST code has a number of astrophysical parameters to characterize the
galaxy properties. In Table 1, we mention some of the important astrophysical
parameters used in this work. We follow the same method to build the data set
and utilize the same pipeline developed in Section IV to emulate the 21-cm
power spectrum and other observables in the presence of MCGs. We discuss the
performance of this emulator in Appendix B. Overall, we find that we get
similar results even after introducing the additional parameters for MCGs.

In Fig. 5 we show the importance of accounting for the MCGs when interpreting
the HERA data. plotting the power spectrum for $L_{\rm X<2\,keV}/{\rm
SFR}=10^{39}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$ where
we use only ACGs in the simulation, we see that this value of $L_{\rm
X<2\,keV}/{\rm SFR}$ is clearly disfavoured by the HERA observations. This is
consistent with the findings in Ref. [29]. However, when adding MCGs to our
simulations, we find that the same value of $L_{\rm X<2\,keV}/{\rm SFR}$ is
consistent with the HERA data. A prediction for this behaviour, was already
mentioned in Ref. [22]. The bottom line is that with MCGs, even the smallest
values of $L_{\rm X<2\,keV}/{\rm SFR}$ in our prior range cannot be ruled out
by current HERA data. Later, we shall see the same in the posterior
distribution of $L_{\rm X<2\,keV}/{\rm SFR}$ from our full MCMC analysis.

Figure 5: Comparison between simulated power spectra with and without MCGs,
and the emulator performance for each signal. The parameters taken for the
simulation are the medians of the marginalized posterior distributions
depicted in Fig. 8, for all the parameters (MCG parameters were not included
in the “without MCGs” simulation) except $L_{\rm X<2\,keV}/{\rm SFR}$, for
which we have taken $10^{39}\rm erg\,s^{-1}\,M_{\odot}^{-1}\,yr$. We quote
some of the other parameters used in this figure for convenience:
$\log_{10}f_{\ast,10}=-1.25$, $\log_{10}f_{\ast}=-2.65$, $\alpha_{\ast}=0.52$,
$\alpha_{\ast,{\rm mini}}=-0.03$, $\log_{10}f_{\rm esc,10}=-1.45$,
$\log_{10}f_{\rm esc,7}=-2.31$, $\alpha_{\rm esc}=0.55$. Note that, for the
“without MCGs” scenario (red curves), we have considered only ACGs in the
simulations (same as the scenario assumed in Fig. 7) with the appropriate
parameter values as mentioned above, ignoring any contribution from the MCGs.

##  VI Results

###  VI.1 Bayesian inference

We now perform the MCMC analysis based on our pipeline that combines the
emulator discussed in Section IV and the likelihood mentioned in Section III.
In the MCMC analysis, we have two main goals: (I) To reproduce the results of
Ref. [29] which includes only ACGs in the analysis, and (II) To study the
impact of MCGs on the X-ray heating constraints. Considering the two goals, we
have two scenarios for the likelihood: (A) Use a likelihood for only the
additional observables mentioned in Section III.3, which we denote as
‘External likelihood’, and (B) Consider a ‘Full likelihood’ that consists of
the external likelihood and the likelihood for the HERA power spectrum upper
bound. We, therefore, have four different MCMC inferences and these are:
(I-A): ACGs + External likelihood, (I-B): ACGs + Full likelihood, (II-A): MCGs
+ External likelihood, and (II-B): MCGs + Full likelihood. We discuss the
results below. For the MCMC runs, we use the Emcee [140] sampler. In table 2,
we show the prior ranges used for this work and we use uniform flat priors for
all the variables.

For inferences of type I-A and I-B, we have used the same prior range as in
Ref. [29]. For inferences with mini halos, we drop some parameters in order to
reduce the parameter space. We find that $t_{\ast}$ (0.5) and $\alpha_{X}$
(1.0) parameters are not very well constrained. Therefore, even if we drop
them, there should not be any major effect on the inference. We keep them
fixed at the values mentioned within the brackets. Furthermore, in principle,
PopIII stars should have different X-ray luminosity (denoted by
$L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$) than PopII stars (denoted by
$L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$). However, we have no prior knowledge
of $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$ of the PopIII stars, which reside
mainly in MCGs. Considering this uncertainty, and, in order to reduce the
number of parameters, we assume $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$ =
$L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$ and call this the effective $L_{\rm
X<2\,keV}/{\rm SFR}$. The prior range of the MCG parameters is taken to be the
same as the ACG parameters.

Note that the above choice may affect the prior volume non-trivially. Making
$L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$ = $L^{\rm(II)}_{\rm X<2\,keV}/{\rm
SFR}$ automatically reduces the effective value of $L_{\rm X<2\,keV}/{\rm
SFR}$. This would not be the case if $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$
is made to vary independently. However, as the PopIII stars are metal-poor,
our naive expectation is that $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}\gtrsim
L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$. Therefore, we expect our results to be
consistent and even somewhat conservative under the assumption we made.

Priors
---
Parameter name | Lower bound | Upper bound
$\log_{10}f_{\ast,10}$ | -3.0 | 0.0
$\log_{10}f_{\ast,7}$ | -3.5 | -1.0
$\alpha_{\ast}$ | -0.5 | 1.0
$\alpha_{\ast,\rm mini}$ | -0.5 | 0.5
$\log_{10}f_{{\rm esc},10}$ | -3.0 | 0.0
$\log_{10}f_{{\rm esc},7}$ | -3.0 | 0.0
$\alpha_{\rm esc}$ | -1.0 | 1.0
$\log_{10}[M_{\rm turn}/M_{\odot}]$ | 8.0 | 10.0
$t_{\ast}$ | 0.0 | 1.0
$\log_{10}\frac{L_{\rm X<2\,keV}/{\rm SFR}}{\rm erg\,s^{-1}\,M_{\odot}^{-1}\,yr}$ | 38.0 | 42.0
$E_{0}/{\rm keV}$ | 0.1 | 1.5
$\alpha_{X}$ | -1.0 | 3.0
Table 2: Prior range summary for all the astrophysical parameters used in the
different MCMC runs. Note that only a subset of them is used in each run, as
explained in the text.

###  VI.2 Weaker constraints on $L_{\rm X}$

In the top panel of Fig. 6, we show the marginalized posteriors for $L_{\rm
X<2\,keV}/{\rm SFR}$, the parameter most constrained by HERA, for inferences
I-A and I-B. This figure is directly comparable to Fig. 28 of Ref. [22] (see
also Fig. 7 in Ref. [29]) and our result is largely in agreement with theirs.
The external likelihoods do not impose any constraint on $L_{\rm
X<2\,keV}/{\rm SFR}$. The HERA data helps to place bounds at low $L_{\rm
X<2\,keV}/{\rm SFR}$, and we can see that values below $L_{\rm X<2\,keV}/{\rm
SFR}\sim 10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$ are
heavily disfavoured, while higher values are preferred.

Figure 6: $L_{\rm X<2\,keV}/{\rm SFR}$ marginalized posterior distribution.
The dark shaded area denotes 68% highest posterior density (HPD), and the
light shaded area denotes 95% HPD. The dashed black line represents the
average $L_{\rm X<2\,keV}/{\rm SFR}$ as inferred from observations of high-
mass X-ray binaries (HMXBs) in local, star-forming galaxies [30]. The dotted
black line denotes the theoretical prediction for metal-free HMXBs presented
in Ref. [148] (such galaxies are expected to be common in the early Universe).
The average observed value is outside the 2$\sigma$ region of the posterior
distribution in the top plot, suggesting a possible tension between the
observations. However, the bottom plot implies that this tension is relaxed
when including MCGs in the 21-cm simulations.

The amplitude of the power spectrum increases as $L_{\rm X<2\,keV}/{\rm SFR}$
is decreased, and below $L_{\rm X<2\,keV}/{\rm SFR}\sim 10^{39}\,{\rm
erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$ (as can be seen from Fig.
5, where we change the $L_{\rm X<2\,keV}/{\rm SFR}$ values while keeping the
other parameters at the median values as mentioned in Fig. 7), the power
spectrum amplitude exceeds the HERA $1\;\sigma$ upper bound. On the other
hand, we find an almost flat posterior towards the large $L_{\rm
X<2\,keV}/{\rm SFR}$. This result somewhat disagrees with Ref. [22] where the
$L_{\rm X<2\,keV}/{\rm SFR}$ posterior shows a steady decline at $L_{\rm
X<2\,keV}/{\rm SFR}>10^{40.8}\,{\rm erg}\;{\rm s}^{-1}\;{\rm
M}_{\odot}^{-1}\;{\rm yr}$. This is possibly due to the unequally-weighted
sample points of the MultiNest sampler [142] used by 21cmMC [143]. The key
point in Fig. 6 is that for models where only ACGs (or PopII stars) are
present, the HERA upper limit suggests more X-ray luminosity for the high-
redshift ($z>6$) sources.

The bottom panel of Fig. 6 shows the marginalized posteriors for $L_{\rm
X<2\,keV}/{\rm SFR}$, for inferences II-A and II-B. As discussed earlier, both
these inferences are based on models where we include MCGs (or PopIII stars)
and ACGs (or PopII stars), and $L_{\rm X<2\,keV}/{\rm SFR}$ here represents an
effective value for the X-ray luminosity of these galaxies. The External
likelihoods do not place any bounds on $L_{\rm X<2\,keV}/{\rm SFR}$ within the
prior range we used, which is expected. However, with HERA data, we find a
striking difference when compared with the top panel. We see that the strong
preference for the high $L_{\rm X<2\,keV}/{\rm SFR}$ ($L_{\rm X<2\,keV}/{\rm
SFR}>10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$) values
is severely relaxed as soon as we include MCGs. This is also evident from Fig.
5. Although there is a decline in posterior probability below $L_{\rm
X<2\,keV}/{\rm SFR}\sim 10^{39.5}\,{\rm erg}\;{\rm s}^{-1}\;{\rm
M}_{\odot}^{-1}\;{\rm yr}$, this is not as significant as in the top panel
where the decline is very sharp at $L_{\rm X<2\,keV}/{\rm SFR}<10^{40}\,{\rm
erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$. Overall, Fig. 6 clearly
indicates that the lower $L_{\rm X<2\,keV}/{\rm SFR}$ values are still
possible when we add MCGs to the analysis.

###  VI.3 Posteriors for other astrophysical parameters

In Fig. 7, we show 1D and 2D marginalized posterior distributions for
inferences I-A and I-B. This figure is directly comparable to Fig. 6 of Ref.
[29]. We find that our results resemble the results of Ref. [29] except for
one detail. In Ref. [29], we see that there is a sudden decline in the 1D
posterior at the two ends for parameters $\alpha_{\rm X}$ and $E_{0}$, and at
the highest end for $t_{\ast}$ and $L_{\rm X<2\,keV}/{\rm SFR}$. The reason
for this in the case of $L_{\rm X<2\,keV}/{\rm SFR}$ was discussed previously,
and it is likely due to the unequally-weighted sample points in 21cmMC, and we
suspect the same is happening here for the other parameters.

Figure 7: Posterior distribution with and without HERA data. Posteriors
without HERA data are produces using upper bound on the neutral fraction
$\bar{x}_{HI}$ at $z=5.9$, $\tau_{e}$ \- CMB optical depth, and galaxy UV
luminosity function at redshift $z$ = 6 - 10 ($M_{\rm UV}>-20$). The most
significant effect of adding the HERA data, is ruling out the possibility of
low values of $L_{\rm X<2\,keV}/{\rm SFR}$ and setting the 68% highest
posterior density (HPD) on $L_{\rm X<2\,keV}/{\rm SFR}=[10^{40.6},10^{42}]$
$\rm erg\,s^{-1}\,M_{\odot}^{-1}\,yr$, which is in good agreement with the
results presented at [22], $L_{\rm X<2\,keV}/{\rm SFR}=[10^{40.4},10^{41.7}]$
$\rm erg\,s^{-1}\,M_{\odot}^{-1}\,yr$. Figure 8: Posterior distribution with
and without HERA data, including MCGs. The most significant difference from
what presented at figure 7, are the relaxed constraints on $L_{\rm
X<2\,keV}/{\rm SFR}$, setting the 68% highest posterior density (HPD) on
$L_{\rm X<2\,keV}/{\rm SFR}=[10^{40},10^{42}]$ $\rm
erg\,s^{-1}\,M_{\odot}^{-1}\,yr$, and the 95% HPD on $[10^{38.6},10^{42}]$
$\rm erg\,s^{-1}\,M_{\odot}^{-1}\,yr$.

Fig. 8 is the same as Fig. 7 but for inferences II-A and II-B. Here also most
of the parameters are largely constrained by the external data, and HERA data
only helps to put a lower bound on $L_{\rm X<2\,keV}/{\rm SFR}$. Median values
of most of the ACG parameters are almost similar to Fig. 7 except for $f_{{\rm
esc},10}$. With MCGs in the scenario, we find that a smaller median value of
$f_{{\rm esc},10}$ is now preferred. This is expected as less radiation is
required to escape from the ACGs since MCGs also contribute to the ionizing
photon budget. We also see that $\log_{10}f_{\ast,10}$ here is more tightly
constrained than in Fig. 7, although the median values are similar in both
figures. Considering the MCG parameters, we find that the total likelihood
places only upper bounds on $\log_{10}f_{\ast,7}$ and $\log_{10}f_{{\rm
esc},7}$, and almost no statistically significant bound on $\alpha_{\ast,\rm
mini}$.

###  VI.4 Considering different priors on MCGs parameters

As discussed in [45], it is important to note that the constraints over
$L_{\rm X<2\,keV}/{\rm SFR}$, achieved when introducing MCGs, are highly
dependent on the allowed values of $\log_{10}f_{\ast,7}$ and $\log_{10}f_{{\rm
esc},7}$, i.e. on the choice of priors. In the above analysis, we did not
allow $\log_{10}f_{\ast,7}$ to be smaller than -3.5. When lower values are
permitted, the external likelihoods push $\log_{10}f_{\ast,7}$ towards them.
The star-formation-rate density (SFRD) of MCGs depends linearly on
$f_{\ast,7}$ [47], thus it decreases as well. In turn, that decreases the
number of X-ray producing remnants in MCGs, which increases the overall
$L_{\rm X<2\,keV}/{\rm SFR}$ required so that the power spectrum amplitude
remains below the HERA upper bounds. The resulting constraints are very
similar to the ones produced when taking only ACGs into account, as shown in
Fig. 9 (left panel). The opposite effect occurs when increasing the lower
bound on $\log_{10}f_{\ast,7}$. Higher values of $\log_{10}f_{\ast,7}$, i.e.
higher MCGs SFRD, yield more X-ray producing sources, thus decreasing the
required $L_{\rm X<2\,keV}/{\rm SFR}$ for the signal to be consistent with the
current HERA data. This behaviour is illustrated as well in Fig. 9 (left
panel).

Figure 9: Marginalized posterior distribution of $L_{\rm X<2\,keV}/{\rm SFR}$
obtained when including MCGs, while varying the lower bound of
$\log_{10}f_{\ast,7}$ (Left) or the lower bound of $\log_{10}f_{{\rm esc},7}$
(Right). When lower values of $\log_{10}f_{\ast,7}$ are enabled, the posterior
distribution tends towards the results of I-B (ACGs only). In case we restrict
$\log_{10}f_{\ast,7}$ to a higher range, the $L_{\rm X<2\,keV}/{\rm SFR}$
constraints tend to vanish, as less X-ray luminosity is now required to lower
the amplitude of the signal at the redshifts measured by HERA. Similarly, when
higher values of $\log_{10}f_{{\rm esc},7}$ are forced, the posterior
distribution tends towards the results of I-B (ACGs only). If we allow
$\log_{10}f_{{\rm esc},7}$ to be smaller, the $L_{\rm X<2\,keV}/{\rm SFR}$
constraints tend to relax.

Similar behaviour is seen when the lower bounds on $\log_{10}f_{{\rm esc},7}$
are varied. Our fiducial value is $\log_{10}f_{{\rm esc},7}$ $>-3$. When
$\log_{10}f_{{\rm esc},7}$ is forced to be larger, ionization occurs earlier,
thus increasing the optical depth ($\tau$ is calculated as an integral, where
the integrand is proportional to $1-x_{HI}$). In order to fit for the optical
depth constraint from Ref. [27], $\log_{10}f_{\ast,7}$ must be smaller, and so
the $L_{\rm X<2\,keV}/{\rm SFR}$ posterior tends towards the one achieved in
analysis I-B, as discussed above. In the opposite scenario, we enable lower
values of $\log_{10}f_{{\rm esc},7}$, and two related processes take place.
Since $\log_{10}f_{{\rm esc},7}$ is smaller, ionization occurs later, such
that the optical depth is smaller as well. In addition, the neutral fraction
at $z=5.9$, also increases when $\log_{10}f_{{\rm esc},7}$ is lowered. The
constraints on these two observables force a higher value of
$\log_{10}f_{\ast,7}$, and as discussed above this enhances the effect of
including MCGs is the analysis. The $L_{\rm X<2\,keV}/{\rm SFR}$ posteriors in
the different cases are shown in Fig. 9 (right panel). The bottom line is that
the above analysis is prior dependent, and caution is warranted when reaching
conclusions for models with MCGs.

##  VII Conclusions

In this work, we have introduced an ANN-based emulator of the 21-cm signal
which was trained on data produced by the 21cmFAST semi-numerical code. We
found this emulator to be accurate enough to predict the 21-cm signal over a
range of redshifts $z$ and wave numbers $k$.

We then used this emulator in an MCMC pipeline to reproduce the parameter
constraints presented in Refs. [29, 22], based on the HERA phase-I upper
limits on the 21-cm power spectrum and the external likelihoods containing
information on the high-redshift UV luminosity function, the IGM neutral
fraction and the reionization optical depth. Here, the underlying
astrophysical model assumption is that the atomically cooled galaxies (ACGs)
(which host the PopII stars) sparked the cosmic dawn and are responsible for
the subsequent heating and reionization of the IGM. Our emulator-based MCMC
pipeline seems to produce similar results as presented in Ref. [29], which
validates our pipeline. One of the most important results of this analysis,
based on the posterior of the parameter $L_{\rm X<2\,keV}/{\rm SFR}$, can be
stated as follows. If PopII stars (or ACGs) dominate the X-ray heating and
reionization of the IGM, then we expect the high redshift ($z>6$) galaxies to
be more X-ray luminous (with $L_{\rm X<2\,keV}/{\rm SFR}\gtrsim 10^{40}\,{\rm
erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$) than their present-day
counterparts (for which $L_{\rm X<2\,keV}/{\rm SFR}\approx 10^{39}\,{\rm
erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$) [30, 31].

Our primary goal in this work was to check the validity of the above result
when we include molecular-cooled galaxies (MCGs) in the analysis. A number of
recent simulations show that MCGs, which predominantly host PopIII stars, also
contribute to the total photon budget during cosmic dawn and reionization.
Ignoring MCGs in any analysis could therefore lead to an overestimation of the
X-ray and ionizing efficiencies of the atomic cooling galaxies (ACGs). The
same concern has also been mentioned in Ref. [22].

We therefore trained our emulator with the results obtained after including
PopIII stars in the simulations, and finally ran an MCMC analysis that fits
for additional parameters corresponding to the PopIII stars. The most
interesting result that emerges out of the final posterior distributions we
obtained is that including both ACGs and MCGs in the simulations relaxes the
preference for the high X-ray luminosity of the high redshift sources as we
discussed in the previous paragraph, and now values as low as $L_{\rm
X<2\,keV}/{\rm SFR}\lesssim 10^{39}\,{\rm erg}\;{\rm s}^{-1}\;{\rm
M}_{\odot}^{-1}\;{\rm yr}$ are still allowed by the HERA power spectrum data.
This is due to the fact that MCGs contribute to the X-ray heating and
therefore the ACGs do not need to be as X-ray efficient. This indicates that
the X-ray luminosity of the high redshift sources may not be very different
from their low redshift counterparts. It is important to note, that although
we do see a small decline in the posterior distribution at the lower $L_{\rm
X<2\,keV}/{\rm SFR}$ regime when including MCGs (as in Figs. 1 and 6) we
cannot determine confidently that higher values are preferred. This is because
the chosen likelihood function (see Eq. (8)) gives higher probability for
power spectra well below the HERA upper bounds. In fact, as we see in Fig. 5,
the spectra of this low X-ray signal is just below the HERA bound, making it a
valid scenario. Contrary to that, the same figure shows that when MCGs are not
included, the spectra obtained with large $L_{\rm X<2\,keV}/{\rm SFR}$ is well
above the upper bound, As mentioned above, this result depends on the choice
of priors, and specifically on the allowed ranges of $\log_{10}f_{\ast,7}$ and
$\log_{10}f_{\ast,7}$, and caution has to be taken when applying this
conclusions.

Before closing, we check the validity of our results by comparing them with
the 21-cm global signal $\langle\delta T_{21}\rangle$ measured by the Shaped
Antenna Measurement of the Background Radio Spectrum 3 (SARAS 3) [144]
experiment in the frequency band $55\rm MHz$ to $85\rm MHz$. Note that SARAS 3
did not yield a statistically significant detection of $\langle\delta
T_{21}\rangle$, rather they provide RMS residuals of amplitude $213\rm mK$
which can be considered as upper and lower bounds of the signal. Fig. 10 shows
that our results do not contradict the SARAS 3 results. Note that existing
global signal [145, 83, 146, 147, 144] measurements do not have much
constraining power and hence are not included in our analysis. However, with
more precise experiments in the future, global signal measurements will
achieve more constraining power. In that case, our pipeline can be easily
extended to emulate the global signals [53, 52, 54]. We shall explore this
possibility in future work.

Figure 10: Comparison of the 21-cm global signal obtained with median values
from the posterior distribution from the inference which accounts for MCGs,
and different $L_{\rm X<2\,keV}/{\rm SFR}$ values. These values are compared
to the SARAS 3 brightness temperature bounds, demonstrating that introducing
MCGs while enabling lower values of $L_{\rm X<2\,keV}/{\rm SFR}$ does not
create a conflict with earlier cosmological measurements.

Note that in our analysis, we have assumed that the X-ray luminosities of both
the PopII and PopIII stars are same, i.e. $L^{\rm(II)}_{\rm X<2\,keV}/{\rm
SFR}$ = $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$ = $L_{\rm X<2\,keV}/{\rm
SFR}$. We call $L_{\rm X<2\,keV}/{\rm SFR}$ as the effective X-ray luminosity.
This choice is made just to reduce the total number of parameters. However,
this choice affects the prior volume in a non-trivial manner. Making
$L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$ = $L^{\rm(III)}_{\rm X<2\,keV}/{\rm
SFR}$ automatically reduces the effective value of $L_{\rm X<2\,keV}/{\rm
SFR}$. This would not be the case if $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}$
is made to vary independently. In this case, if the likelihood had preferred a
lower value $L^{\rm(III)}_{\rm X<2\,keV}/{\rm SFR}<<L^{\rm(II)}_{\rm
X<2\,keV}/{\rm SFR}$, we would have obtained the same conclusion as in Ref
[29]. However, this is unlikely as the theoretical modelling indicates that
the low-metallicity stars, like PopIII, are expected to exhibit high $L_{\rm
X<2\,keV}/{\rm SFR}$ [148], and we anticipate $L^{\rm(III)}_{\rm
X<2\,keV}/{\rm SFR}\gtrsim L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$. Therefore,
we expect that our conclusions will still be valid under the assumption
$L^{\rm(II)}_{\rm X<2\,keV}/{\rm SFR}$ = $L^{\rm(III)}_{\rm X<2\,keV}/{\rm
SFR}$ = $L_{\rm X<2\,keV}/{\rm SFR}$.

Future study will also try to address the limitations of the emulator used in
this work, most significant of which is the inability to predict the flat or
zero signals reliably. We have proposed a workaround considering the present
scenario. However, with a more precise measurement of the power spectrum, we
will need to predict signals of all kinds with better precision. Below, we
discuss a few possible approaches to improving our emulation pipeline: (i)
Using Generative Adversarial Networks (GANs) [149], which implement some
unsupervised learning methods to train models. (ii) Training the ANN on the
likelihood function instead of the power spectrum. In that case, the ANN does
not have to learn the complicated features of the signal. Rather, it learns to
make decisions as to how acceptable the signal is based on the likelihood.
(iii) Since in future experiments we expect to have more redshift bands, we
can treat the power-spectra observations as two dimensional data, where one
dimension is the redshift, and the other is the wave number. We therefore can
exploit the developing research in the field of Convolutional Neural Networks
(CNNs) [150, 151], which uses for machine learning algorithms where the data
has more than one dimension, for our goal. A different approach to the
computational problem of producing samples of our model in a reasonable time,
is using analytic codes like Zeus21 [71]. This code provides good
approximations for all of 21cmFAST (ACG-only) outputs, in
$\mathcal{O}(1-10)\,{\rm sec}$.

The conclusions of our work may have a number of consequences which we note
below: (i) MCGs spark cosmic dawn and consequently, the Ly-$\alpha$ coupling
starts early. This allows the spin temperature to catch up with the kinetic
temperature at a higher redshift compared to the scenario where only ACGs
exist. Now, a combination of less X-ray luminous ($L_{\rm X<2\,keV}/{\rm
SFR}\lesssim 10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm
yr}$) MCG and ACGs would make the absorption signal deeper and extended over
several redshifts (compared to more luminous $L_{\rm X<2\,keV}/{\rm
SFR}>10^{40}\,{\rm erg}\;{\rm s}^{-1}\;{\rm M}_{\odot}^{-1}\;{\rm yr}$
sources). This bodes well for global signal experiments as deeper and more
extended absorption signals are easier to detect. (ii) Since the amplitude of
the 21-cm power spectrum increases as we decrease $L_{\rm X<2\,keV}/{\rm
SFR}$, this could also be fortunate for experiments like HERA and SKA that
target 21-cm fluctuations. (iii) Low X-ray efficiency of the early sources
provides room for other heating mechanisms to be effective, such as
Ly-$\alpha$ heating due to resonant scattering of the Ly-$\alpha$ photons with
the atoms in IGM [77, 152].

It is crucial to take caution when interpreting results related to the
modelling of the 21cm signal at high redshifts, as these outcomes heavily
hinge on the intricate modelling of various processes. Notably, different
codes employing similar source populations but employing distinct algorithms
may yield disparate results. 21cmFAST is one such model that uses a specific
algorithm to model the cosmic dawn and reionization. It would be worth
checking the consistency of the results presented here with other similar
codes.

Finally, we emphasize that the emulator-based analysis pipeline developed in
this work is flexible, fast and does not require a prohibitive number of
simulations. We therefore plan to employ it to place 21-cm bounds on models of
dark matter [153, 154, 155], primordial magnetic fields [156] and various
other extensions to $\Lambda$CDM.

###### Acknowledgements.

We thank Rennan Barkana, Daniela Breitman, Anastasia Fialkov, Jordan Flitter
and Caner Unal for useful discussions. We are grateful to Andrei Mesinger and
Julian Muñoz for insightful comments on the manuscript. We also acknowledge
the efforts of the entire 21cmFAST team in producing a state-of-the-art 21-cm
public simulation code. HL is supported by a Dkalim excellence M.Sc.
scholarship from Ben-Gurion University. EDK is supported by an Azrieli
Foundation faculty fellowship.

##  Appendix A Emulators for $\tau_{e}$ and $x_{HI}$

The emulators for $\tau_{e}$ and $x_{\rm HI}$ are built on the same structure
described in Section IV. The data set for each observable contain $\sim$ 20000
simulations where the parameters are sampled using the LH [117] sampler. The
whole data is randomly divided into training (85%), validation (10%) and
testing (5%) sets. The astrophysical parameters are normalized to fit in the
range $[-1,1]$ for convenience. The NN here is composed of 3 layers of sizes
256, 512, 512 respectively. We have used the ReLU [120] activation function,
and we have checked that it performs really well. As done for the emulators
described in section IV, we added a batch Normalization layer [122] after each
activation. For the loss function, we rely on the standard L2 norm, defined as

$\mathcal{L}(y_{true},y_{pred}(\theta))=\frac{1}{m}\sum_{i=1}^{m}(y_{true}-y_{pred}(\theta))^{2}\,,$
(14)

where $m$ is the number of samples in a batch which is taken to be 128. Here,
the Adam [123] optimizer is used in the training with an initial learning rate
of 0.01. Depending on the performance of the ANN on the validation set, we
even reduce the learning rate by a factor of two. In the end, we restore the
best weights.

Figure 11: Loss statistics for the $\tau$ emulator based on testing sets
having $\sim$ 1000 samples each. The dashed lines represent the 16th and 84th
percentiles respectively, and their exact value is written along with the
median value in the figure legend. Note that the X-axis here is not normalized
with the true values.  Figure 12:  Loss statistics for the $x_{\rm HI}$
emulator based on testing sets having $\sim$ 1000 samples each. The dashed
lines represent the 16th and 84th percentiles respectively, and their exact
value is written along with the median value in the figure legend. Note that
the X-axis here is not normalized with the true values.

As can be seen from Figures 11 and 12, the emulators perform really well. Even
the rare errors are not significant. This convinces us that the constraints
obtained from these observables will not suffer from emulation errors.

##  Appendix B Emulator with MCGs

Here we present the performance of our emulator where the simulations used to
train it include both ACGs and MCGs. The pipeline and procedures used here are
the same as described in Section IV, and we do not repeat these. We present
the emulation results in Figure 13. We find that the results are almost
similar to what we found for the emulator where only ACGs are present.

Figure 13:  Loss statistics for the power spectrum emulator, where MCGs are
included in the simulations. Results are based on testing sets having $\sim$
250 samples in each. The dashed lines represent the 16th and 84th percentiles
respectively, and their exact value is written along with the median value in
the figure legend.

##  Appendix C Classifiers

In this work, we are using a total of four different signal classifiers for
the power spectrum. One for each redshift $(7.9,10.4)$, and for each redshift,
there are two sub-classifiers for emulation with (i) only ACGs and (ii)
ACGs+MCGs. For each classification, we have used $\sim 10000$ signals and
divided those into training (85%) validation (10%) and testing (5%) sets. We
have also normalized the astrophysical parameters to fit the range $[-1,1]$
for faster and more efficient learning.

The NN for all the classifiers is composed of three layers of sizes 128, 512,
and 256 respectively. For activation, we have used ReLU [120] followed by a
batch Normalization layer [122] which standardizes the input for the next
layer, and a Dropout [hinton2012improving] layer, which randomly turns off
$10\%$ of the neurons in the previous layer in order to prevent over-fitting.
The end layer contains one neuron that uses a ‘sigmoid’ activation function.
This is chosen in order to have the classifier output in the range $\in[0,1]$.
Here, an output smaller than $0.5$ is classified as class 0 and $\geq 0.5$ is
classified as class 1 (as defined in Section IV.1.5). We use the Binary cross
Entropy loss function, which is the standard choice for classifying a data set
into two classes, defined as,

$\mathcal{L}=-\frac{1}{m}\sum_{i=1}^{m}y_{i}\log(P(y_{i}))+(1-y_{i})\log(1-P(y_{i}))\,,$
(15)

where $y_{i}$ are the labels, $P(y_{i})$ is the classifier probability
prediction for this label, and $m$ is the number of samples in a batch which
we have chosen to be $512$. We have used the Adam [123] optimizer with an
initial learning rate of 0.01, which is reduced by a factor of two based on
the performance. We store the best weights at the end. We found that the
accuracy is $\gtrsim 97\%$ for all the classifiers, and we can confidently use
this in the MCMC process.

## References

  * [1] P. Madau, A. Meiksin and M. J. Rees, “21-CM tomography of the intergalactic medium at high redshift,” Astrophys. J. 475, 429 (1997) [arXiv:astro-ph/9608010 [astro-ph]].
  * [2] A. Loeb and M. Zaldarriaga, “Measuring the small - scale power spectrum of cosmic density fluctuations through 21 cm tomography prior to the epoch of structure formation,” Phys. Rev. Lett. 92, 211301 (2004) [arXiv:astro-ph/0312134 [astro-ph]].
  * [3] J. R. Pritchard and A. Loeb, “21-cm cosmology,” Rept. Prog. Phys. 75, 086901 (2012) [arXiv:1109.6012 [astro-ph.CO]].
  * [4] J. C. Pober, A. Liu, J. S. Dillon, J. E. Aguirre, J. D. Bowman, R. F. Bradley, C. L. Carilli, D. R. DeBoer, J. N. Hewitt and D. C. Jacobs, et al. “What Next-Generation 21 cm Power Spectrum Measurements Can Teach Us About the Epoch of Reionization,” Astrophys. J. 782, 66 (2014) [arXiv:1310.7031 [astro-ph.CO]].
  * [5] A. Liu et al. [Snowmass 2021 Cosmic Frontier 5 Topical Group], “Snowmass2021 Cosmic Frontier White Paper: 21cm Radiation as a Probe of Physics Across Cosmic Ages,” [arXiv:2203.07864 [astro-ph.CO]].
  * [6] S. Furlanetto, S. P. Oh and F. Briggs, “Cosmology at Low Frequencies: The 21 cm Transition and the High-Redshift Universe,” Phys. Rept. 433, 181-301 (2006) [arXiv:astro-ph/0608032 [astro-ph]].
  * [7] R. Barkana and A. Loeb, “In the beginning: The First sources of light and the reionization of the Universe,” Phys. Rept. 349, 125-238 (2001) [arXiv:astro-ph/0010468 [astro-ph]].
  * [8] R. Barkana and A. Loeb, “Detecting the earliest galaxies through two new sources of 21cm fluctuations,” Astrophys. J. 626, 1-11 (2005) [arXiv:astro-ph/0410129 [astro-ph]].
  * [9] X. L. Chen and J. Miralda-Escude, “The 21cm Signature of the First Stars,” Astrophys. J. 684, 18-33 (2008) [arXiv:astro-ph/0605439 [astro-ph]].
  * [10] J. Miralda-Escude, M. Haehnelt and M. J. Rees, “Reionization of the inhomogeneous universe,” Astrophys. J. 530, 1-16 (2000) [arXiv:astro-ph/9812306 [astro-ph]].
  * [11] S. Furlanetto, M. Zaldarriaga and L. Hernquist, “The Growth of HII regions during reionization,” Astrophys. J. 613, 1-15 (2004) [arXiv:astro-ph/0403697 [astro-ph]].
  * [12] X. Chen, P. D. Meerburg and M. Münchmeyer, “The Future of Primordial Features with 21 cm Tomography,” JCAP 09, 023 (2016) [arXiv:1605.09364 [astro-ph.CO]].
  * [13] G. Mellema, L. V. E. Koopmans, F. A. Abdalla, G. Bernardi, B. Ciardi, S. Daiboo, A. G. de Bruyn, K. K. Datta, H. Falcke and A. Ferrara, et al. “Reionization and the Cosmic Dawn with the Square Kilometre Array,” Exper. Astron. 36, 235-318 (2013) [arXiv:1210.0197 [astro-ph.CO]].
  * [14] A. P. Beardsley, B. J. Hazelton, I. S. Sullivan, P. Carroll, N. Barry, M. Rahimi, B. Pindor, C. M. Trott, J. Line and D. C. Jacobs, et al. “First Season MWA EoR Power Spectrum Results at Redshift 7,” Astrophys. J. 833, no.1, 102 (2016) [arXiv:1608.06281 [astro-ph.IM]].
  * [15] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen and N. Mahesh, “An absorption profile centred at 78 megahertz in the sky-averaged spectrum,” Nature 555, no.7694, 67-70 (2018) [arXiv:1810.05912 [astro-ph.CO]].
  * [16] C. M. Trott and J. C. Pober, “The status of 21cm interferometric experiments,” [arXiv:1909.12491 [astro-ph.IM]].
  * [17] B. Greig, C. M. Trott, N. Barry, S. J. Mutch, B. Pindor, R. L. Webster and J. S. B. Wyithe, “Exploring reionization and high-z galaxy observables with recent multiredshift MWA upper limits on the 21-cm signal,” Mon. Not. Roy. Astron. Soc. 500, no.4, 5322-5335 (2020) [arXiv:2008.02639 [astro-ph.CO]].
  * [18] F. G. Mertens, M. Mevius, L. V. E. Koopmans, A. R. Offringa, G. Mellema, S. Zaroubi, M. A. Brentjens, H. Gan, B. K. Gehlot and V. N. Pandey, et al. “Improved upper limits on the 21-cm signal power spectrum of neutral hydrogen at $\bm{z\approx 9.1}$ from LOFAR,” Mon. Not. Roy. Astron. Soc. 493, no.2, 1662-1685 (2020) [arXiv:2002.07196 [astro-ph.CO]].
  * [19] R. Ansari, “Current status and future of cosmology with 21cm Intensity Mapping,” [arXiv:2203.17039 [astro-ph.CO]].
  * [20] D. R. DeBoer, A. R. Parsons, J. E. Aguirre, P. Alexander, Z. S. Ali, A. P. Beardsley, G. Bernardi, J. D. Bowman, R. F. Bradley and C. L. Carilli, et al. “Hydrogen Epoch of Reionization Array (HERA),” Publ. Astron. Soc. Pac. 129, no.974, 045001 (2017) [arXiv:1606.07473 [astro-ph.IM]].
  * [21] Z. Abdurashidova et al. [HERA], “First Results from HERA Phase I: Upper Limits on the Epoch of Reionization 21 cm Power Spectrum,” Astrophys. J. 925, no.2, 221 (2022) [arXiv:2108.02263 [astro-ph.CO]].
  * [22] Z. Abdurashidova et al. [HERA], “Improved Constraints on the 21 cm EoR Power Spectrum and the X-Ray Heating of the IGM with HERA Phase I Observations,” Astrophys. J. 945, no.2, 124 (2023) [arXiv:2210.04912 [astro-ph.CO]].
  * [23] R. J. Bouwens, G. D. Illingworth, P. A. Oesch, M. Trenti, I. Labbé, L. Bradley, M. Carollo, P. G. van Dokkum, V. Gonzalez and B. Holwerda, et al. “UV Luminosity Functions at redshifts $z\sim$4 to $z\sim$10: 10000 Galaxies from HST Legacy Fields,” Astrophys. J. 803, no.1, 34 (2015) [arXiv:1403.4295 [astro-ph.CO]].
  * [24] R. J. Bouwens, P. A. Oesch, I. Labbé,G. D. Illingworth, G. G. Fazio, D. Coe, B. Holwerda, R. Smit, M. Stefanon, P. G. van Dokkum, M. Trenti, M. L. N. Ashby, J. S. Huang, L. Spitler, C. Straatman, L. Bradley, D. Magee, “The Bright End of the z 9 and z 10 UV Luminosity Functions Using All Five CANDELS Fields,” Astrophys. J. 830, no.2, 67 (2016) [arXiv:1506.01035 [astro-ph.GA]]
  * [25] P. A. Oesch, R. J. Bouwens, G. D. Illingworth, I. Labbe, M. Stefanon “The Dearth of z 10 Galaxies in All HST Legacy Fields—The Rapid Evolution of the Galaxy Population in the First 500 Myr”, Astrophys. J. 855, 105 (2017) [arXiv:1710.11131 [astro-ph.CO]].
  * [26] N. J. F. Gillet, A. Mesinger and J. Park, “Combining high-z galaxy luminosity functions with Bayesian evidence,” Mon. Not. Roy. Astron. Soc. 491, no.2, 1980-1997 (2020) [arXiv:1906.06296 [astro-ph.GA]].
  * [27] N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * [28] I. McGreer, A. Mesinger and V. D’Odorico, “Model-independent evidence in favour of an end to reionization by $z\approx$ 6,” Mon. Not. Roy. Astron. Soc. 447, no.1, 499-505 (2015) [arXiv:1411.5375 [astro-ph.CO]].
  * [29] Z. Abdurashidova et al. [HERA], “HERA Phase I Limits on the Cosmic 21 cm Signal: Constraints on Astrophysics and Cosmology during the Epoch of Reionization,” Astrophys. J. 924, no.2, 51 (2022) [arXiv:2108.07282 [astro-ph.CO]].
  * [30] S. Mineo, M. Gilfanov and R. Sunyaev, “X-ray emission from star-forming galaxies - I. High-mass X-ray binaries,” Mon. Not. Roy. Astron. Soc. 419, 2095 (2012) [arXiv:1105.4610 [astro-ph.HE]].
  * [31] B. D. Lehmer, D. M. Alexander, F. E. Bauer, W. N. Brandt, A. D. Goulding, L. P. Jenkins, A. Ptak and T. P. Roberts, “A Chandra Perspective On Galaxy-Wide X-ray Binary Emission And Its Correlation With Star Formation Rate And Stellar Mass: New Results From Luminous Infrared Galaxies,” Astrophys. J. 724, 559-571 (2010) [arXiv:1009.3943 [astro-ph.CO]].
  * [32] B. Smith, J. Wise, B. O’Shea, M. Norman and S. Khochfar, “The first Population II stars formed in externally enriched mini-haloes,” Mon. Not. Roy. Astron. Soc. 452, no.3, 2822-2836 (2015) [arXiv:1504.07639 [astro-ph.GA]].
  * [33] V. Bromm, “Formation of the first stars,” IAU Symp. 228, 121 (2005) [arXiv:astro-ph/0509354 [astro-ph]].
  * [34] V. Bromm and R. B. Larson, “The First stars,” Ann. Rev. Astron. Astrophys. 42, 79-118 (2004) [arXiv:astro-ph/0311019 [astro-ph]].
  * [35] Z. Haiman, M. J. Rees and A. Loeb, “H(2) cooling of primordial gas triggered by UV irradiation,” Astrophys. J. 467, 522 (1996) [arXiv:astro-ph/9511126 [astro-ph]].
  * [36] Z. Haiman, M. J. Rees and A. Loeb, “Destruction of molecular hydrogen during cosmological reionization,” Astrophys. J. 476, 458 (1997) [arXiv:astro-ph/9608130 [astro-ph]].
  * [37] M. Tegmark, J. Silk, M. J. Rees, A. Blanchard, T. Abel and F. Palla, “How small were the first cosmological objects?,” Astrophys. J. 474, 1-12 (1997) [arXiv:astro-ph/9603007 [astro-ph]].
  * [38] T. Abel, G. L. Bryan and M. L. Norman, “The formation of the first star in the Universe,” Science 295, 93 (2002) [arXiv:astro-ph/0112088 [astro-ph]].
  * [39] N. Yoshida, K. Omukai, L. Hernquist and T. Abel, “Formation of Primordial Stars in a lambda-CDM Universe,” Astrophys. J. 652, 6-25 (2006) [arXiv:astro-ph/0606106 [astro-ph]].
  * [40] Z. Haiman and G. L. Bryan, Astrophys. J. 650, 7-11 (2006) doi:10.1086/506580 [arXiv:astro-ph/0603541 [astro-ph]].
  * [41] M. Trenti, “Population III Star Formation During and After the Reionization Epoch,” AIP Conf. Proc. 1294, no.1, 134-137 (2010) [arXiv:1006.4434 [astro-ph.CO]].
  * [42] M. Kulkarni, E. Visbal and G. L. Bryan, “The Critical Dark Matter Halo Mass for Population III Star Formation: Dependence on Lyman–Werner Radiation, Baryon-dark Matter Streaming Velocity, and Redshift,” Astrophys. J. 917, no.1, 40 (2021) [arXiv:2010.04169 [astro-ph.GA]].
  * [43] R. H. Mebane, J. Mirocha and S. R. Furlanetto, “The effects of population III radiation backgrounds on the cosmological 21-cm signal,” Mon. Not. Roy. Astron. Soc. 493, no.1, 1217-1226 (2020)
  * [44] J. Mirocha, R. H. Mebane, S. R. Furlanetto, K. Singal and D. Trinh, “Unique signatures of Population III stars in the global 21-cm signal,” Mon. Not. Roy. Astron. Soc. 478, no.4, 5591-5606 (2018) [arXiv:1710.02530 [astro-ph.GA]].
  * [45] Y. Qin, A. Mesinger, J. Park, B. Greig and J. B. Muñoz, “A tale of two sites – I. Inferring the properties of minihalo-hosted galaxies from current observations,” Mon. Not. Roy. Astron. Soc. 495, no.1, 123-140 (2020) [arXiv:2003.04442 [astro-ph.CO]].
  * [46] Y. Qin, A. Mesinger, B. Greig and J. Park, “A tale of two sites – II. Inferring the properties of minihalo-hosted galaxies with upcoming 21-cm interferometers,” Mon. Not. Roy. Astron. Soc. 501, no.4, 4748-4758 (2021) [arXiv:2009.11493 [astro-ph.CO]].
  * [47] J. B. Muñoz, Y. Qin, A. Mesinger, S. G. Murray, B. Greig and C. Mason, “The impact of the first galaxies on cosmic dawn and reionization,” Mon. Not. Roy. Astron. Soc. 511, no.3, 3657-3681 (2022) [arXiv:2110.13919 [astro-ph.CO]].
  * [48] H. Xu, K. Ahn, J. H. Wise, M. L. Norman and B. W. O’Shea, “Heating the Intergalactic medium by X-rays from Population III Binaries in High Redshift Galaxies,” Astrophys. J. 791, no.2, 110 (2014) [arXiv:1404.6555 [astro-ph.CO]].
  * [49] H. Xu, K. Ahn, M. L. Norman, J. H. Wise and B. W. O’Shea, “X-ray Background at High Redshifts from Pop III Remnants: Results from Pop III star formation rates in the Renaissance Simulations,” Astrophys. J. Lett. 832, no.1, L5 (2016) [arXiv:1607.02664 [astro-ph.GA]].
  * [50] N. S. Kern, A. Liu, A. R. Parsons, A. Mesinger and B. Greig, “Emulating Simulations of Cosmic Dawn for 21 cm Power Spectrum Constraints on Cosmology, Reionization, and X-Ray Heating,” Astrophys. J. 848, no.1, 23 (2017) [arXiv:1705.04688 [astro-ph.CO]].
  * [51] W. D. Jennings, C. A. Watkinson, F. B. Abdalla and J. D. McEwen, “Evaluating machine learning techniques for predicting power spectra from reionization simulations,” Mon. Not. Roy. Astron. Soc. 483, no.3, 2907-2922 (2019) [arXiv:1811.09141 [astro-ph.CO]].
  * [52] A. Cohen, A. Fialkov, R. Barkana and R. Monsalve, “Emulating the Global 21-cm Signal from Cosmic Dawn and Reionization,” [arXiv:1910.06274 [astro-ph.CO]].
  * [53] C. H. Bye, S. K. N. Portillo and A. Fialkov, “21cmVAE: A Very Accurate Emulator of the 21 cm Global Signal,” Astrophys. J. 930, no.1, 79 (2022) [arXiv:2107.05581 [astro-ph.CO]].
  * [54] H. T. J. Bevins, W. J. Handley, A. Fialkov, E. d. Acedo and K. Javid, “globalemu: a novel and robust approach for emulating the sky-averaged 21-cm signal from the cosmic dawn and epoch of reionization,” Mon. Not. Roy. Astron. Soc. 508, no.2, 2923-2936 (2021) [arXiv:2104.04336 [astro-ph.CO]].
  * [55] M. Choudhury, A. Datta and S. Majumdar, “Extracting the 21-cm power spectrum and the reionization parameters from mock data sets using artificial neural networks,” Mon. Not. Roy. Astron. Soc. 512, no.4, 5010-5022 (2022) [arXiv:2112.13866 [astro-ph.CO]].
  * [56] S. Sikder, R. Barkana, I. Reis and A. Fialkov, “Machine Learning to Decipher the Astrophysical Processes at Cosmic Dawn,” [arXiv:2201.08205 [astro-ph.CO]].
  * [57] B. Ciardi, A. Ferrara, S. Marri and G. Raimondo, Mon. Not. Roy. Astron. Soc. 324, 381 (2001) doi:10.1046/j.1365-8711.2001.04316.x [arXiv:astro-ph/0005181 [astro-ph]].
  * [58] M. G. Santos, L. Ferramacho, M. B. Silva, A. Amblard and A. Cooray, Mon. Not. Roy. Astron. Soc. 406, 2421-2432 (2010) doi:10.1111/j.1365-2966.2010.16898.x [arXiv:0911.2219 [astro-ph.CO]].
  * [59] P. R. Shapiro, I. T. Iliev, G. Mellema, K. Ahn, Y. Mao, M. Friedrich, K. Datta, H. Park, E. Komatsu and E. Fernandez, et al. AIP Conf. Proc. 1480, no.1, 248-260 (2012) doi:10.1063/1.4754363 [arXiv:1211.0583 [astro-ph.CO]].
  * [60] S. Majumdar, G. Mellema, K. K. Datta, H. Jensen, T. R. Choudhury, S. Bharadwaj and M. M. Friedrich, Mon. Not. Roy. Astron. Soc. 443, no.4, 2843-2861 (2014) doi:10.1093/mnras/stu1342 [arXiv:1403.0941 [astro-ph.CO]].
  * [61] S. Hassan, R. Davé, K. Finlator and M. G. Santos, Mon. Not. Roy. Astron. Soc. 457, no.2, 1550-1567 (2016) doi:10.1093/mnras/stv3001 [arXiv:1510.04280 [astro-ph.CO]].
  * [62] D. Sarkar, S. Bharadwaj and S. Anathpindika, Mon. Not. Roy. Astron. Soc. 460, no.4, 4310-4319 (2016) doi:10.1093/mnras/stw1111 [arXiv:1605.02963 [astro-ph.CO]].
  * [63] R. Ghara, G. Mellema, S. K. Giri, T. R. Choudhury, K. K. Datta and S. Majumdar, Mon. Not. Roy. Astron. Soc. 476, no.2, 1741-1755 (2018) doi:10.1093/mnras/sty314 [arXiv:1710.09397 [astro-ph.CO]].
  * [64] A. Mesinger, IAU Symp. 333, 3-11 (2017) doi:10.1017/S1743921317011139 [arXiv:1801.02649 [astro-ph.CO]].
  * [65] D. Sarkar and S. Bharadwaj, Mon. Not. Roy. Astron. Soc. 476, no.1, 96-108 (2018) doi:10.1093/mnras/sty206 [arXiv:1801.07868 [astro-ph.CO]].
  * [66] P. Ocvirk, D. Aubert, J. G. Sorce, P. R. Shapiro, N. Deparis, T. Dawoodbhoy, J. Lewis, R. Teyssier, G. Yepes and S. Gottlöber, et al. Mon. Not. Roy. Astron. Soc. 496, no.4, 4087-4107 (2020) doi:10.1093/mnras/staa1266 [arXiv:1811.11192 [astro-ph.GA]].
  * [67] D. Sarkar and S. Bharadwaj, Mon. Not. Roy. Astron. Soc. 487, no.4, 5666-5678 (2019) doi:10.1093/mnras/stz1691 [arXiv:1906.07032 [astro-ph.CO]].
  * [68] M. Molaro, R. Davé, S. Hassan, M. G. Santos and K. Finlator, Mon. Not. Roy. Astron. Soc. 489, no.4, 5594-5611 (2019) doi:10.1093/mnras/stz2171 [arXiv:1901.03340 [astro-ph.CO]].
  * [69] R. Kannan, E. Garaldi, A. Smith, R. Pakmor, V. Springel, M. Vogelsberger and L. Hernquist, Mon. Not. Roy. Astron. Soc. 511, no.3, 4005-4030 (2022) doi:10.1093/mnras/stab3710 [arXiv:2110.00584 [astro-ph.GA]].
  * [70] A. K. Shaw, A. Chakraborty, M. Kamran, R. Ghara, S. Choudhuri, S. S. Ali, S. Pal, A. Ghosh, J. Kumar and P. Dutta, et al. J. Astrophys. Astron. 44, no.1, 4 (2023) doi:10.1007/s12036-022-09889-6 [arXiv:2211.05512 [astro-ph.CO]].
  * [71] J. B. Muñoz, doi:10.1093/mnras/stad1512 [arXiv:2302.08506 [astro-ph.CO]].
  * [72] B. Maity and T. R. Choudhury, Mon. Not. Roy. Astron. Soc. 521, no.3, 4140-4155 (2023) doi:10.1093/mnras/stad791 [arXiv:2211.12909 [astro-ph.CO]].
  * [73] J. B. Muñoz, Phys. Rev. Lett. 123, no.13, 131301 (2019) doi:10.1103/PhysRevLett.123.131301 [arXiv:1904.07868 [astro-ph.CO]].
  * [74] J. B. Muñoz, “Robust Velocity-induced Acoustic Oscillations at Cosmic Dawn,” Phys. Rev. D 100, no.6, 063538 (2019) doi:10.1103/PhysRevD.100.063538 [arXiv:1904.07881 [astro-ph.CO]].
  * [75] A. Mesinger, S. Furlanetto and R. Cen, “21cmFAST: A Fast, Semi-Numerical Simulation of the High-Redshift 21-cm Signal,” Mon. Not. Roy. Astron. Soc. 411, 955 (2011) [arXiv:1003.3878 [astro-ph.CO]].
  * [76] J. Mirocha, J. B. Muñoz, S. R. Furlanetto, A. Liu and A. Mesinger, Mon. Not. Roy. Astron. Soc. 514, no.2, 2010-2030 (2022) doi:10.1093/mnras/stac1479 [arXiv:2201.07249 [astro-ph.CO]].
  * [77] I. Reis, A. Fialkov and R. Barkana, “The subtlety of Ly $\alpha$ photons: changing the expected range of the 21-cm signal,” Mon. Not. Roy. Astron. Soc. 506, no.4, 5479-5493 (2021) [arXiv:2101.01777 [astro-ph.CO]].
  * [78] J. Dowell and G. B. Taylor, Astrophys. J. Lett. 858, no.1, L9 (2018) doi:10.3847/2041-8213/aabf86 [arXiv:1804.08581 [astro-ph.CO]].
  * [79] D. J. Fixsen, A. Kogut, S. Levin, M. Limon, P. Lubin, P. Mirel, M. Seiffert, J. Singal, E. Wollack and T. Villela, et al. Astrophys. J. 734, 5 (2011) doi:10.1088/0004-637X/734/1/5 [arXiv:0901.0555 [astro-ph.CO]].
  * [80] S. Bharadwaj and S. S. Ali, “The CMBR fluctuations from HI perturbations prior to reionization,” Mon. Not. Roy. Astron. Soc. 352, 142 (2004) [arXiv:astro-ph/0401206 [astro-ph]].
  * [81] S. Furlanetto, A. Lidz, A. Loeb, M. McQuinn, J. Pritchard, J. Aguirre, M. Alvarez, D. Backer, J. Bowman and J. Burns, et al. “Astrophysics from the Highly-Redshifted 21 cm Line,” [arXiv:0902.3011 [astro-ph.CO]].
  * [82] P. Madau and F. Haardt, “Cosmic Reionization after Planck: Could Quasars Do It All?,” Astrophys. J. Lett. 813, no.1, L8 (2015) [arXiv:1507.07678 [astro-ph.CO]].
  * [83] R. A. Monsalve, A. E. E. Rogers, J. D. Bowman and T. J. Mozdzen, “Results from EDGES High-Band: I. Constraints on Phenomenological Models for the Global $21$ cm Signal,” Astrophys. J. 847, no.1, 64 (2017) [arXiv:1708.05817 [astro-ph.CO]].
  * [84] H. T. J. Bevins, A. Fialkov, E. d. Acedo, W. J. Handley, S. Singh, R. Subrahmanyan and R. Barkana, “Astrophysical constraints from the SARAS 3 non-detection of the cosmic dawn sky-averaged 21-cm signal,” Nature Astron. 6, no.12, 1473-1483 (2022) [arXiv:2212.00464 [astro-ph.CO]].
  * [85] D. Vrbanec, B. Ciardi, V. Jelic, H. Jensen, I. T. Iliev, G. Mellema and S. Zaroubi, Mon. Not. Roy. Astron. Soc. 492, no.4, 4952-4958 (2020) doi:10.1093/mnras/staa183 [arXiv:2001.08814 [astro-ph.CO]].
  * [86] D. Karagiannis, J. Fonseca, R. Maartens and S. Camera, Phys. Dark Univ. 32, 100821 (2021) doi:10.1016/j.dark.2021.100821 [arXiv:2010.07034 [astro-ph.CO]].
  * [87] M. Rahimi, B. Pindor, J. L. B. Line, N. Barry, C. M. Trott, R. L. Webster, C. H. Jordan, M. Wilensky, S. Yoshiura and A. Beardsley, et al. Mon. Not. Roy. Astron. Soc. 508, no.4, 5954-5971 (2021) doi:10.1093/mnras/stab2918 [arXiv:2110.03190 [astro-ph.CO]].
  * [88] J. H. Cook, C. M. Trott and J. L. B. Line, Mon. Not. Roy. Astron. Soc. 514, no.1, 790-805 (2022) doi:10.1093/mnras/stac1330 [arXiv:2205.04644 [astro-ph.CO]].
  * [89] M. Berti, M. Spinelli and M. Viel, Mon. Not. Roy. Astron. Soc. 521, no.3, 3221-3236 (2023) doi:10.1093/mnras/stad685 [arXiv:2209.07595 [astro-ph.CO]].
  * [90] J. S. W. Lewis, A. Pillepich, D. Nelson, R. S. Klessen and S. C. O. Glover, [arXiv:2305.09721 [astro-ph.GA]].
  * [91] S. Yoshiura, H. Shimabukuro, K. Takahashi, R. Momose, H. Nakanishi and H. Imai, Mon. Not. Roy. Astron. Soc. 451, no.1, 266-274 (2015) doi:10.1093/mnras/stv855 [arXiv:1412.5279 [astro-ph.CO]].
  * [92] H. Shimabukuro, S. Yoshiura, K. Takahashi, S. Yokoyama and K. Ichiki, Mon. Not. Roy. Astron. Soc. 458, no.3, 3003-3011 (2016) doi:10.1093/mnras/stw482 [arXiv:1507.01335 [astro-ph.CO]].
  * [93] D. Sarkar, S. Majumdar and S. Bharadwaj, Mon. Not. Roy. Astron. Soc. 490, no.2, 2880-2889 (2019) doi:10.1093/mnras/stz2799 [arXiv:1907.01819 [astro-ph.CO]].
  * [94] R. Durrer, M. Jalilvand, R. Kothari, R. Maartens and F. Montanari, JCAP 12, 003 (2020) doi:10.1088/1475-7516/2020/12/003 [arXiv:2008.02266 [astro-ph.CO]].
  * [95] S. Cunnington, C. Watkinson and A. Pourtsidou, Mon. Not. Roy. Astron. Soc. 507, no.2, 1623-1639 (2021) doi:10.1093/mnras/stab2200 [arXiv:2102.11153 [astro-ph.CO]].
  * [96] M. Kamran, S. Majumdar, R. Ghara, G. Mellema, S. Bharadwaj, J. R. Pritchard, R. Mondal and I. T. Iliev, [arXiv:2108.08201 [astro-ph.CO]].
  * [97] D. Karagiannis, R. Maartens and L. F. Randrianjanahary, JCAP 11, 003 (2022) doi:10.1088/1475-7516/2022/11/003 [arXiv:2206.07747 [astro-ph.CO]].
  * [98] M. Tegmark, “How to measure CMB power spectra without losing information,” Phys. Rev. D 55, 5895-5907 (1997) [arXiv:astro-ph/9611174 [astro-ph]].
  * [99] A. Liu and M. Tegmark, “A Method for 21cm Power Spectrum Estimation in the Presence of Foregrounds,” Phys. Rev. D 83, 103006 (2011) [arXiv:1103.0281 [astro-ph.CO]].
  * [100] J. S. Dillon, A. Liu, C. L. Williams, J. N. Hewitt, M. Tegmark, E. H. Morgan, A. M. Levine, M. F. Morales, S. J. Tingay and G. Bernardi, et al. “Overcoming real-world obstacles in 21 cm power spectrum estimation: A method demonstration and results from early Murchison Widefield Array data,” Phys. Rev. D 89, no.2, 023002 (2014) [arXiv:1304.4229 [astro-ph.CO]].
  * [101] A. Liu and J. R. Shaw, “Data Analysis for Precision 21 cm Cosmology,” Publ. Astron. Soc. Pac. 132, no.1012, 062001 (2020) [arXiv:1907.08211 [astro-ph.IM]].
  * [102] S. Cunnington, L. Wolz, A. Pourtsidou and D. Bacon, “Impact of foregrounds on HI intensity mapping cross-correlations with optical surveys,” Mon. Not. Roy. Astron. Soc. 488, no.4, 5452-5472 (2019) [arXiv:1904.01479 [astro-ph.CO]].
  * [103] M. Kolopanis, D. C. Jacobs, C. Cheng, A. R. Parsons, S. A. Kohn, J. C. Pober, J. E. Aguirre, Z. S. Ali, G. Bernardi and R. F. Bradley, et al. “A simplified, lossless re-analysis of PAPER-64,” [arXiv:1909.02085 [astro-ph.CO]].
  * [104] Y. Qin, V. Poulin, A. Mesinger, B. Greig, S. Murray and J. Park, “Reionization inference from the CMB optical depth and E-mode polarization power spectra,” Mon. Not. Roy. Astron. Soc. 499, no.1, 550-558 (2020) [arXiv:2006.16828 [astro-ph.CO]].
  * [105] J. Silk, A. Di Cintio and I. Dvorkin, Proc. Int. Sch. Phys. Fermi 186, 137-187 (2014) doi:10.3254/978-1-61499-476-3-137 [arXiv:1312.0107 [astro-ph.CO]].
  * [106] M. Vogelsberger, S. Genel, D. Sijacki, P. Torrey, V. Springel and L. Hernquist, Mon. Not. Roy. Astron. Soc. 436, 3031 (2013) doi:10.1093/mnras/stt1789 [arXiv:1305.2913 [astro-ph.CO]].
  * [107] R. S. Somerville and R. Davé, Ann. Rev. Astron. Astrophys. 53, 51-113 (2015) doi:10.1146/annurev-astro-082812-140951 [arXiv:1412.2712 [astro-ph.GA]].
  * [108] A. Pillepich, V. Springel, D. Nelson, S. Genel, J. Naiman, R. Pakmor, L. Hernquist, P. Torrey, M. Vogelsberger and R. Weinberger, et al. Mon. Not. Roy. Astron. Soc. 473, no.3, 4077-4106 (2018) doi:10.1093/mnras/stx2656 [arXiv:1703.02970 [astro-ph.GA]].
  * [109] M. Vogelsberger, F. Marinacci, P. Torrey and E. Puchwein, Nature Rev. Phys. 2, no.1, 42-66 (2020) doi:10.1038/s42254-019-0127-2 [arXiv:1909.07976 [astro-ph.GA]].
  * [110] A. Wetzel, C. C. Hayward, R. E. Sanderson, X. Ma, D. Angles-Alcazar, R. Feldmann, T. K. Chan, K. El-Badry, C. Wheeler and S. Garrison-Kimmel, et al. Astrophys. J. Suppl. 265, no.2, 44 (2023) doi:10.3847/1538-4365/acb99a [arXiv:2202.06969 [astro-ph.GA]].
  * [111] S. G. Murray, B. Greig, A. Mesinger, J. B. Muñoz, Y. Qin, J. Park and C. A. Watkinson, J. Open Source Softw. 5, no.54, 2582 (2020) doi:10.21105/joss.02582 [arXiv:2010.15121 [astro-ph.IM]].
  * [112] D. Sarkar and E. D. Kovetz, Phys. Rev. D 107, no.2, 023524 (2023) doi:10.1103/PhysRevD.107.023524 [arXiv:2210.16853 [astro-ph.CO]].
  * [113] G. Shmueli, D. Sarkar and E. D. Kovetz, [arXiv:2305.07056 [astro-ph.CO]].
  * [114] C. Mason, M. Trenti and T. Treu, Astrophys. J. 813, no.1, 21 (2015) doi:10.1088/0004-637X/813/1/21 [arXiv:1508.01204 [astro-ph.GA]].
  * [115] J. Park, A. Mesinger, B. Greig and N. Gillet, Mon. Not. Roy. Astron. Soc. 484, no.1, 933-949 (2019) doi:10.1093/mnras/stz032 [arXiv:1809.08995 [astro-ph.GA]].
  * [116] A. Rudakovskyi, A. Mesinger, D. Savchenko and N. Gillet, Mon. Not. Roy. Astron. Soc. 507, no.2, 3046-3056 (2021) doi:10.1093/mnras/stab2333 [arXiv:2104.04481 [astro-ph.CO]].
  * McKay et al. [1979] McKay, M. D., Beckman, R. J., & Conover, W. J. 1979, Technometrics, 21, 239
  * [118] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow et al., “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.” (2015), Software available from: https://www.tensorflow.org/.
  * [119] F. Chollet et al., “Keras.” (2015), Software available from: https://keras.io.
  * [120] K. Fukushima. “Visual Feature Extraction by a Multilayered Network of Analog Threshold Elements”, IEEE Transactions on Systems Science and Cybernetics, 5, no. 4, pp. 322-333 (1969).
  * [121] A. Spurio Mancini, D. Piras, J. Alsing, B. Joachimi and M. P. Hobson, “CosmoPower: emulating cosmological power spectra for accelerated Bayesian inference from next-generation surveys,” Mon. Not. Roy. Astron. Soc. 511, no.2, 1771-1788 (2022) [arXiv:2106.03846 [astro-ph.CO]].
  * [122] S. Ioffe and C. Szegedy, “Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift,” [arXiv:1502.03167 [cs.LG]].
  * [123] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” [arXiv:1412.6980 [cs.LG]].
  * [124] B. W. O’Shea and M. L. Norman, Astrophys. J. 673, 14 (2008) doi:10.1086/524006 [arXiv:0706.4416 [astro-ph]].
  * [125] C. A. Mason, A. Fontana, T. Treu, K. B. Schmidt, A. Hoag, L. Abramson, R. Amorin, M. Bradac, L. Guaita and T. Jones, et al. Mon. Not. Roy. Astron. Soc. 485, no.3, 3947-3969 (2019) doi:10.1093/mnras/stz632 [arXiv:1901.11045 [astro-ph.CO]].
  * [126] T. R. Choudhury, A. Paranjape and S. E. I. Bosman, Mon. Not. Roy. Astron. Soc. 501, no.4, 5782-5796 (2021) doi:10.1093/mnras/stab045 [arXiv:2003.08958 [astro-ph.CO]].
  * [127] E. Visbal, Z. Haiman and G. L. Bryan, Mon. Not. Roy. Astron. Soc. 453, no.4, 4456-4466 (2015) doi:10.1093/mnras/stv1941 [arXiv:1505.06359 [astro-ph.CO]].
  * [128] A. Fialkov, R. Barkana, E. Visbal, D. Tseliakhovich and C. M. Hirata, “The 21-cm signature of the first stars during the Lyman-Werner feedback era,” Mon. Not. Roy. Astron. Soc. 432, 2909 (2013) [arXiv:1212.0513].
  * [129] E. Visbal, Z. Haiman, B. Terrazas, G. L. Bryan and R. Barkana, “High-redshift star formation in a time-dependent Lyman–Werner background,” Mon. Not. Roy. Astron. Soc. 445, no.1, 107-114 (2014) [arXiv:1402.0882].
  * [130] C. Safranek-Shrader, M. Agarwal, C. Federrath, A. Dubey, M. Milosavljevic and V. Bromm, “Star Formation in the First Galaxies I: Collapse Delayed by Lyman-Werner Radiation,” Mon. Not. Roy. Astron. Soc. 426, 1159 (2012) [arXiv:1205.3835].
  * [131] M. Ricotti, N. Y. Gnedin and J. M. Shull, “Feedback from galaxy formation: production and photodissociation of primordial molecular hydrogen,” Astrophys. J. 560, 580 (2001) [arXiv:astro-ph/0012335].
  * [132] K. Ahn, P. R. Shapiro, I. T. Iliev, G. Mellema and U. L. Pen, “The Inhomogeneous Background of Hydrogen-Molecule Dissociating Radiation during Cosmic Reionization,” Astrophys. J. 695, no.2, 1430-1445 (2009) [arXiv:0807.2254].
  * [133] A. Fialkov, “Supersonic Relative Velocity between Dark Matter and Baryons: A Review,” Int. J. Mod. Phys. D 23, no.08, 1430017 (2014) [arXiv:1407.2274].
  * [134] R. Barkana, “The Rise of the First Stars: Supersonic Streaming, Radiative Feedback, and 21-cm Cosmology,” Phys. Rept. 645, 1-59 (2016) [arXiv:1605.04357].
  * [135] D. Tseliakhovich and C. Hirata, “Relative velocity of dark matter and baryonic fluids and the formation of the first structures,” Phys. Rev. D 82, 083520 (2010) [arXiv:1005.2416].
  * [136] J. Bovy and C. Dvorkin, “Low-mass suppression of the satellite luminosity function due to the supersonic baryon–cold-dark-matter relative velocity,” Astrophys. J. 768, 70 (2013) [arXiv:1205.2083].
  * [137] A. Stacy, V. Bromm and A. Loeb, “Effect of Streaming Motion of Baryons Relative to Dark Matter on the Formation of the First Stars,” Astrophys. J. Lett. 730, no.1, L1 (2011) [arXiv:1011.4512].
  * [138] A. Fialkov, R. Barkana, D. Tseliakhovich and C. M. Hirata, “Impact of the Relative Motion between the Dark Matter and Baryons on the First Stars,” Mon. Not. Roy. Astron. Soc. 424, 1335-1345 (2012) [arXiv:1110.2111].
  * [139] F. Schmidt, “Effect of relative velocity and density perturbations between baryons and dark matter on the clustering of galaxies,” Phys. Rev. D 94, no.6, 063508 (2016) [arXiv:1602.09059].
  * [140] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, Publ. Astron. Soc. Pac. 125, 306-312 (2013) doi:10.1086/670067 [arXiv:1202.3665 [astro-ph.IM]].
  * [141] T. Fragos, B. Lehmer, M. Tremmel, P. Tzanavaris, A. Basu-Zych, K. Belczynski, A. Hornschemeier, L. Jenkins, V. Kalogera and A. Ptak, et al. “X-ray Binary Evolution Across Cosmic Time,” Astrophys. J. 764, 41 (2013) [arXiv:1206.2395 [astro-ph.HE]].
  * [142] F. Feroz, M. P. Hobson and M. Bridges, Mon. Not. Roy. Astron. Soc. 398, 1601-1614 (2009) doi:10.1111/j.1365-2966.2009.14548.x [arXiv:0809.3437 [astro-ph]].
  * [143] Y. Qin, private communication.
  * [144] S. Singh, J. Nambissan T., R. Subrahmanyan, N. Udaya Shankar, B. S. Girish, A. Raghunathan, R. Somashekar, K. S. Srivani and M. Sathyanarayana Rao, “On the detection of a cosmic dawn signal in the radio background,” Nature Astron. 6, no.5, 607-617 (2022) [arXiv:2112.06778 [astro-ph.CO]].
  * [145] J. D. Bowman and A. E. E. Rogers, “A lower limit of dz $>$ 0.06 for the duration of the reionization epoch,” Nature 468, 796-798 (2010) [arXiv:1209.1117 [astro-ph.CO]].
  * [146] R. A. Monsalve, B. Greig, J. D. Bowman, A. Mesinger, A. E. E. Rogers, T. J. Mozdzen, N. S. Kern and N. Mahesh, “Results from EDGES High-Band: II. Constraints on Parameters of Early Galaxies,” Astrophys. J. 863, no.1, 11 (2018) [arXiv:1806.07774 [astro-ph.CO]].
  * [147] R. A. Monsalve, A. Fialkov, J. D. Bowman, A. E. E. Rogers, T. J. Mozdzen, A. Cohen, R. Barkana and N. Mahesh, “Results from EDGES High-Band: III. New Constraints on Parameters of the Early Universe,” Astrophys. J. 875, no.1, 67 (2019) [arXiv:1901.10943 [astro-ph.CO]].
  * [148] T. Fragos, B. Lehmer, M. Tremmel, P. Tzanavaris, A. Basu-Zych, K. Belczynski, A. Hornschemeier, L. Jenkins, V. Kalogera and A. Ptak, et al. Astrophys. J. 764, 41 (2013) doi:10.1088/0004-637X/764/1/41 [arXiv:1206.2395 [astro-ph.HE]].
  * [149] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville and Y. Bengio, “Generative Adversarial Networks,” [arXiv:1406.2661 [stat.ML]].
  * [150] D. C. Ciresan, U. Meier, J. Masci, L. M. Gambardella, J. Schmidhuber, “Proceedings of the 22nd International Joint Conference on Artificial Intelligence,” Vol. 22 (2011), pp. 1237–1242
  * [151] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke and A. Rabinovich, “Going Deeper with Convolutions,” [arXiv:1409.4842 [cs.CV]].
  * [152] D. Sarkar, J. Flitter and E. D. Kovetz, “Exploring delaying and heating effects on the 21-cm signature of fuzzy dark matter,” Phys. Rev. D 105, no.10, 103529 (2022) [arXiv:2201.03355 [astro-ph.CO]].
  * [153] E. D. Kovetz, V. Poulin, V. Gluscevic, K. K. Boddy, R. Barkana and M. Kamionkowski, “Tighter limits on dark matter explanations of the anomalous EDGES 21 cm signal,” Phys. Rev. D 98, no.10, 103529 (2018) [arXiv:1807.11482 [astro-ph.CO]].
  * [154] E. D. Kovetz, I. Cholis and D. E. Kaplan, “Bounds on ultralight hidden-photon dark matter from observation of the 21 cm signal at cosmic dawn,” Phys. Rev. D 99, no.12, 123511 (2019) [arXiv:1809.01139 [astro-ph.CO]].
  * [155] J. Flitter and E. D. Kovetz, “Closing the window on fuzzy dark matter with the 21-cm signal,” Phys. Rev. D 106, no.6, 063504 (2022) [arXiv:2207.05083 [astro-ph.CO]].
  * [156] T. Adi, S. Libanore, H. A. G. Cruz and E. D. Kovetz, “Constraining Primordial Magnetic Fields with Line-Intensity Mapping,” [arXiv:2305.06440 [astro-ph.CO]].