spectrum. At higher temperatures, the modification of the continuum spectrum
gets milder, while the falloff in injected photons is sharper. In both cases,
scattering effects become weaker at higher temperatures, approaching the no-
scattering limit (constant and step function spectrum, respectively for
continuum and injected photons).

Figure 4.3: Spectrum for continuum (continuous lines) and injected (dashed
lines) photons around the Ly$\alpha$ line, for three different temperatures at
redshift $z=10$.

#####  4.3.2 Ly$\alpha$ scattering rate

The deviation of the Ly$\alpha$ intensity due to line scattering can be
quantified defining the quantity

$S_{\alpha,a}=\int dx\phi(x)\frac{J_{a,\nu}(x)}{J_{a,\infty}},$ (4.44)

for $a=c,i$. Note that $S_{\alpha,a}$ depends only on the shape of the
spectrum and not on its normalization. Then, the Ly$\alpha$ scattering rate
can be written from Eq. (4.33) as

$P_{\alpha}=4\pi\sigma_{0}\int
dx\phi(x)J(x)=4\pi\sigma_{0}S_{\alpha}J_{\alpha,\infty},$ (4.45)

where it has been implicitly summed over the continuum and injected
contributions. Hence, from Eqs. (4.27), (4.32) and (4.45), the Ly$\alpha$
coupling coefficient can be written as

$x_{\alpha}=\frac{16\pi\sigma_{0}T_{*}}{27A_{10}T_{\gamma}}S_{\alpha}J_{\alpha,\infty},$
(4.46)

which shall be used to compute the spin temperature. The flux far from the
line, $J_{\alpha,\infty}$, is determined by the astrophysical model and
cosmological evolution, and it will be computed in the next chapter. However,
the factor $S_{\alpha}$ is only determined by the integrated spectrum, and
therefore it does not depend on the specific astrophysical parameters, being a
function of the temperature and redshift (or, more specifically, of
temperature and the Sobolev parameter $\gamma_{S}$). This quantity accounts
for scattering effects in the flux, and is usually of order $\sim 1$. The
integrals of Eqs. (4.39) and (4.42) cannot be solved analytically, although
there are numerical fits for $S_{\alpha,a}$ in the literature, such as the one
from Ref. [371], employed in the code 21cmFAST [378] (widely used in Part II).
It is possible however to approximate the shape of the line profile to obtain
analytic results, as done in the so-called wing approximation, with a
remarkable agreement with numerical computations [376, 372]. Although that
solution is given in terms of hypergeometric functions, an accurate fit
following this approximation can be derived [376], writing

$S_{\alpha}\simeq\exp{[-1.79\,\alpha]},$ (4.47)

with
$\alpha=\eta\left(3A_{\alpha}/(8\pi\Delta\nu_{D}\gamma_{S})\right)^{1/3}=0.717T_{K}^{-2/3}(\tau_{GP}/10^{6})^{1/3}$
101010Other approximate fits slightly more accurate can be found in Ref.
[372].. From this it can be seen that $S_{\alpha}$ is close to unity at large
temperatures and low optical depths (which correspond to low redshifts in a
mostly neutral medium, see Eq. (4.36)), being damped in the opposite case.
This behavior can be verified in Fig. 4.4, which depicts the variation of
$S_{\alpha}$ with the temperature at different redshifts, for the exact
calculation (i.e., integrating Eq. (4.41)) and for the fit in Eq. (4.47). In
the shown cases, spin exchange effects have been neglected, which is a good
enough approximation as long as $T_{K},T_{S}\gg T_{se}=0.40$ K. Otherwise,
$S_{\alpha}$ depends on $T_{S}$, via the modifications of $\eta$ and
$\gamma_{S}$ described above.

Figure 4.4: $S_{\alpha}$ (for continuum photons) as a function of the kinetic
temperature $T_{K}$, at several redshifts. Continuous lines stand for the
exact case, while dotted lines correspond to the simple fit from Eq. (4.47)
[376].

#####  4.3.3 Color temperature

Although the Ly$\alpha$ spectrum is not an equilibrium spectrum, and therefore
it does not have a definite temperature, one can define an effective color
temperature as [372]

$\frac{h}{k_{B}T_{c}(\nu)}\equiv-\frac{1}{f(\nu)}\frac{df(\nu)}{d\nu},$ (4.48)

where $f(\nu)=J/(2\nu^{2}/c^{2})$ is the photon occupation number. Note that,
with the above definition, if $T_{c}$ is a constant of the frequency, one
recovers a thermal spectrum with that temperature (a Maxwell-Boltzmann
distribution for high energies, since stimulated emission is neglected). It
must be emphasized that the color temperature from Eq. (4.48) would not have
to be necessarily the same as the one implicitly defined in Eq. (4.25).
However, it can be shown that both temperatures coincide if the Ly$\alpha$
line profile is narrow enough [379]. One can employ the calculations from the
previous section to predict the color temperature for Eq. (4.48) (and
equivalently for Eq. (4.25) if the previous argument holds). Employing Eq.
(4.39), this temperature can be written as

$\frac{h}{k_{B}T_{c}(x)}=\frac{2\eta}{\Delta\nu_{D}}+\frac{2\gamma_{S}}{\Delta\nu_{D}\phi(x)}\left(1-\frac{J_{\infty}}{J(x)}\right).$
(4.49)

Including the effect of spin exchange (i.e., substituting
$\gamma_{S}\rightarrow\gamma_{S}(1+T_{se}/T_{K})^{-1}$ and
$\eta\rightarrow\eta(1+T_{se}/T_{S})/(1+T_{se}/T_{K})$, as commented above),
and evaluating the color temperature at the center of the line, one gets

$T_{c}=T_{K}\left[\left(\frac{1+T_{se}/T_{S}}{1+T_{se}/T_{K}}\right)+\frac{\gamma_{S}\left(1+T_{se}/T_{K}\right)^{-1}}{\eta\phi(0)}\left(1-\frac{J_{\infty}}{J(0)}\right)\right]^{-1},$
(4.50)

where it can be shown that the second term becomes negligible at both low and
high temperatures [372]. The above one is an implicit equation: one needs
$T_{c}$ to compute $T_{S}$ with Eq. (4.26), but $T_{c}$ also depends on
$T_{S}$. Therefore, it has to be solved by an iterative procedure, as done in
21cmFAST [378]. From the previous equation, we can see that when the coupling
to the kinetic temperature is efficient enough, $T_{c}\simeq T_{K}$ is a good
approximation, as stated in Sec. 4.2. In any case, since $T_{K},T_{S}\gg
T_{se}$ in most cases of interest, $T_{K}$ is a very good approximation for
the color temperature even without a strong Ly$\alpha$ coupling.

#####  4.3.4 Ly$\alpha$ heating

Scattering around the Ly$\alpha$ line, besides modifying the spectral shape of
the radiation field, can also produce an additional heating contribution to
the gas temperature. Part of the energy of the Ly$\alpha$ photon can be left
in the atom as thermal energy, although depending on the gas temperature, the
atoms can also lose energy in this process, becoming a cooling mechanism. This
contribution was firstly studied in Ref. [380], under the assumption of
neglecting the thermal motion of the atoms. This simplification overestimated
the real heating rate, as was shown by the authors of Ref. [370], who
calculated it properly from the kinetic equation.

In order to derive the heating contribution from scatterings around the
Ly$\alpha$ line, we have to multiply the kinetic equation Eq. (4.37) by the
energy $h\nu$ and integrate over frequency. In the case of injected photons,
we have to take into account the energy carried by the photons injected in the
line, $h\nu_{\alpha}\dot{n}_{\alpha}$ [381], where $\dot{n}_{\alpha}=4\pi
n_{\rm HI}\sigma_{0}\mathcal{C}$ is the rate of change in the number of
photons, obtained from the integration of Eq. (4.37) over frequency. This
contribution vanishes in the case of considering only continuum photons, since
scattering processes conserve photon number. Thus, for either the continuum,
$k=c$, or injected case, $k=i$, the heating rate of the gas by Ly$\alpha$
scattering reads

$\displaystyle\mathcal{Q}|_{\rm Ly\alpha,k}$
$\displaystyle=-\frac{4\pi}{c}\int d\nu\;h\nu\;\left(\frac{\partial
J}{\partial t}-H\nu\frac{\partial
J}{\partial\nu}\right)+h\nu_{\alpha}\dot{n}_{\alpha}$ (4.51)
$\displaystyle=-4\pi n_{\rm HI}\sigma_{0}\int
d\nu\;h\nu\;\frac{1}{2}\frac{\partial}{\partial
x}\left[\phi\left(\frac{\partial J}{\partial x}+2\eta J\right)\right]$ (4.52)
$\displaystyle=4\pi n_{\rm HI}\sigma_{0}h\Delta\nu_{D}\int
d\nu\;\frac{1}{2}\phi\left(\frac{\partial J}{\partial x}+2\eta J\right),$
(4.53)

where in the last equality we have integrated by parts. We have also taken
into account the fact that the spectral profile vanishes far away from the
resonance. Employing Eq. (4.39) for continuum photons, or Eq. (4.40) for
injected ones, we can write the heating rate for each channel as [370, 372,
382]

$\mathcal{Q}|_{\rm Ly\alpha,c}=\frac{4\pi
Hh\nu_{\alpha}\Delta\nu_{D}}{c}\;\mathcal{I}_{k}\;J^{(k)}_{\infty},$ (4.54)

where $k=c,i$, and we have defined the continuum and injected heating
integrals, $\mathcal{I}_{c}$ and $\mathcal{I}_{i}$ respectively, as

$\mathcal{I}_{c}=\int_{-\infty}^{\infty}dx\;\frac{\left(J_{c,\infty}-J_{c}(x)\right)}{J_{c,\infty}}$
(4.55)

and

$\mathcal{I}_{i}=\int_{-\infty}^{0}dx\;\frac{\left(J_{i,\infty}-J_{i}(x)\right)}{J_{i,\infty}}-\int^{\infty}_{0}dx\frac{J_{i}(x)}{J_{i,\infty}}.$
(4.56)

The integrals $I_{c}$ and $I_{i}$ encode the details on the scattering effects
and, as happened with the $S_{\alpha}$ factor, do not depend on the
normalization, but only on the redshift and the temperature, or equivalently,
on the temperature and the Sobolev parameter $\gamma_{S}$. Numerical fits for
these quantities, based on the wing approximation, can be found in Ref. [372].
The importance of this heating mechanism in the evolution of the IGM is
however mild, as shall be discussed in the next chapter.

####  4.4 21 cm power spectrum

The most obvious measurement of the 21 cm signal comes from averaging the
incoming radiation from all sky directions, which is known as the global
signal $\overline{\delta T_{b}}$, corresponding to take the mean value of Eq.
(4.21). While this spatially averaged signature could provide a superb source
to understand the IGM evolution, it lacks some information about
inhomogeneities. The signal is enhanced in overdense regions, and thus
sensitive to the pattern of the cosmic web and distribution of galaxies and
halos. Processes such as X-ray heating and reionization are expected to be
highly inhomogeneous, since radiation fields would expand around galaxies, and
their complex behavior cannot be completely encoded in a spatially averaged
quantity. Moreover, the main experimental caveats to measure the global signal
are the strong foregrounds from the galactic synchrotron emission and the
effect of the ionosphere, which are more easily removed if different spatial
scales are considered, as shall be discussed in Sec. 4.5. It is therefore
mandatory to study not only the global signal but also its spatial
fluctuations, via higher order statistics. The two-point correlation function
is the most obvious choice. Given that radiointerferometers extract
information about the Fourier modes of the signal, it is natural to work in
Fourier space, rather than in real space, employing thus the power spectrum,
as done in Sec. 2.2.1 for the matter density field. The two-point statistics
are of remarkable importance, since for gaussian fields, they present all the
information of the distribution.

Defining the normalized fluctuation of the 21 cm signal from Eq. (4.21) as

$\delta_{21}(z,\mathbf{r})=\frac{\delta T_{b}(z,\mathbf{r})-\overline{\delta
T_{b}}(z)}{\overline{\delta T_{b}}(z)},$ (4.57)

one can take its Fourier transform, as defined in Chapter 2, denoted as
$\widetilde{\delta}_{21}(z,\mathbf{k})$. Therefore, equivalently to Eq.
(2.10), the 21cm differential brightness temperature power spectrum is defined
as

$\langle\widetilde{\delta}_{21}(z,\mathbf{k})\widetilde{\delta}_{21}^{*}(z,\mathbf{k}^{\prime})\rangle\equiv(2\pi)^{3}\delta^{\rm
D}(\mathbf{k}-\mathbf{k}^{\prime})P_{21}(z,k)\leavevmode\nobreak\ .$ (4.58)

In the previous equation, the delta function accounts for invariance under
translations. On the other hand, invariance under rotations imply that the
power spectrum only depends on the modulus of the wavevector $\mathbf{k}$.
These facts are the consequences of statistical homogeneity and isotropy,
respectively. Instead of $P_{21}$, it is customary to work with the
dimensionless power spectrum, defined as
$\Delta^{2}_{21}(k,z)=(k^{3}/2\pi^{2})P_{21}(k,z)$. Since this quantity may
present singularities when $\overline{\delta T_{b}}(z)$ vanishes, we rather
employ the product $\overline{\delta T_{b}}^{2}(z)\Delta^{2}_{21}(k,z)$, which
has units of temperature squared.

One can gain insight into the different contributions to the brightness
temperature field by linearizing Eq. (4.21). Employing the spin temperature
equation, Eq. (4.26), the fluctuation of the 21 cm signal can be written as
[354, 353]

$\delta_{21}=\beta\delta+\beta_{x}\delta_{x}+\beta_{\alpha}\delta_{\alpha}+\beta_{T}\delta_{T}-\delta_{\partial
v},$ (4.59)

with each $\delta_{i}$ accounts for the fluctuations of the different fields:
$\delta$ for the matter density field, $\delta_{x}$ for the neutral fraction,
$\delta_{\alpha}$ for the Ly$\alpha$ flux, $\delta_{T}$ for the kinetic
temperature of the gas and $\delta_{\partial v}$ for the line-of-sight
gradient of the peculiar velocities. The functions $\beta_{i}$ weight the
importance of each fluctuation in the 21 cm one, and a straightforward
computation shows

$\beta=1+\frac{x_{c}}{x_{\rm tot}(1+x_{\rm tot})},$ (4.60)
$\beta_{x}=1+\frac{x_{c}^{\rm HH}-x_{c}^{\rm eH}}{x_{\rm tot}(1+x_{\rm
tot})},$ (4.61) $\beta_{\alpha}=\frac{x_{\alpha}}{x_{\rm tot}(1+x_{\rm
tot})},$ (4.62) $\beta_{T}=\frac{T_{\gamma}}{T_{K}-T_{\gamma}}+\frac{1}{x_{\rm
tot}(1+x_{\rm tot})}\left(x_{c}^{\rm eH}\frac{d\,\ln\kappa^{\rm
eH}_{10}}{d\,\ln T_{K}}+x_{c}^{\rm HH}\frac{d\,\ln\kappa^{\rm HH}_{10}}{d\,\ln
T_{K}}\right),$ (4.63)

with $x_{\rm tot}=x_{\alpha}+x_{c}$, and $x_{c}^{\rm eH}$ and $x_{c}^{\rm HH}$
accounting only for the eH and HH contributions to the collisional coupling,
respectively. Therefore, fluctuations in the brightness temperature can be
regarded as a weighted sum over the fluctuations of the different fields on
which it depends. To derive the previous equations, $T_{c}=T_{K}$ has been
assumed for the sake of simplicity, which is a reasonable assumption, as
stated in the previous section.

The unity terms of $\beta$ and $\beta_{x}$ come from the proportionality of
the brightness temperature to the density and the ionization fraction,
respectively, while the second term includes the contribution from collisional
coupling, only relevant at high redshifts or overdense regions. The first term
of the temperature coefficient $\beta_{T}$ accounts for how fast the spin
temperature changes under fluctuations in $T_{K}$, while the second one comes
from the temperature dependence on the collisional coupling coefficients. The
first one presents a singularity when $T_{K}=T_{\gamma}$, but this is
unphysical and disappears in the quantity $\overline{\delta
T_{b}}^{2}(z)\Delta^{2}_{21}(k,z)$. The $\beta_{\alpha}$ is proportional to
the Ly$\alpha$ coupling coefficient. Understanding the behavior of these
coefficients can shed light on the evolution of the power spectrum, which
shall be discussed in the next chapter.

In principle, it is possible to study the evolution of first order
perturbations, by linearizing the evolution equations of the temperature and
ionization fraction (which shall be explained in the next chapter), as done
with the matter density field. This approach has been followed in the
literature [383, 384], but since IGM processes become non-linear at some
moment, its range of validity is limited. Although we shall compute the 21 cm
power spectrum from the non-linear evolved fields in our numerical
simulations, it is useful to overview these linear quantities to qualitatively
understand the evolution of the power spectrum.

It can be argued that the different $\delta_{i}$ may be proportional to the
overdensity $\delta$ in the linear limit. This is easily seen in the case of
the temperature perturbation $\delta_{T}$, before the onset of heating, as
long as the sound speed is spatially constant. Writing the pressure as $P=\rho
k_{B}T_{K}/\mu$, with $\mu$ the molecular weight and $\rho$ the mass density,
$c_{s}^{2}=dP/d\rho=k_{B}T_{K}/\mu+(\rho k_{B}/\mu)dT_{K}/d\rho$ and thus
$\delta_{T}=(\rho/T_{K})dT_{K}/d\rho\delta=\left[c_{s}^{2}/(k_{B}T_{K}/\mu)-1\right]\delta$
[383]. From the linearized ionization equation, before the galactic sources
are switched on, the evolution is dictated by recombinations, whose rate is
proportional to $\delta$, and thus $\delta_{x}\propto\delta$ (the details of
ionization evolution shall be properly discussed in the next chapter). When
X-rays are starting to be efficiently emitted, the previous arguments would
not be valid. Nevertheless, X-ray photons are radiated from the galaxies
placed at the overdensities, and thus would also be related to $\delta$. An
analogous reasoning can be applied to the Ly$\alpha$ flux, $\delta_{\alpha}$.
Finally, in linear theory, the redshift-space distortions term can be easily
computed, getting $\delta_{\partial v}\simeq-{\rm cos}^{2}\theta f\delta$,
where $\theta$ is the angle between the direction of the line of sight and the
wavevector, and $f=d\,\ln D/d\,\ln a$ ($=1$ in the Einstein-de Sitter
universe) [30], with $D$ the growth factor defined in 2.2.1. Therefore, the 21
cm power spectrum approximately traces the matter power spectrum, weighted by
combinations of the $\beta_{i}$ functions. However, the actual behavior of the
different fluctuations implies extra dependences with the scale for each
$\delta_{i}$, since astrophysical processes are highly inhomogeneous [384].
For this reason, the matter power spectrum cannot be trivially extracted from
the 21 cm one.

Finally, it is noteworthy that, unlike the CMB, the 21 cm signal is expected
to be non-gaussian. This is due to the growth of perturbations into the non-
linear regime, where deviations from gaussianity are developed (as happens in
standard perturbation theory at higher order, see, e.g., [385]). Therefore,
the power spectrum could not provide the full statistical information about
its fluctuations. Although not considered in this thesis, other summary
statistics have been exploited in the literature, and could be extremely
useful to understand the signal. For instance, the bispectrum (three-point
correlation function in Fourier space) encodes information not only about the
magnitude of the field but also about its complex phases (which is related to
the filamentarity) [386, 387, 388]. On the other hand, Minkowski functionals
supply information regarding the morphology and shapes of the structures,
especially useful to gain insight into the growth of HII bubbles during the
EoR [389, 390, 391, 392].

####  4.5 Observational tests of the 21 cm cosmological signature

The observation of the redshifted HI cosmological signal is one of the
greatest challenges in current observational astronomy. Thus, the requisites
to measure this signature need to be critically evaluated. Given the strong
foregrounds and the faintness of the signal, achieving a great sensitivity is
a crucial point, and large collecting areas and long integration times are
needed. Unlike the CMB, 21 cm photons may provide snapshots for different
redshifts, which allows attempting a 3D tomographic view of the IGM. In order
to make a 3D map of the universe and to obtain information from different
epochs, broadband antennas are required, spanning over a large range of
frequencies. As already stated, photons observed at a frequency $\nu$ are
redshifted from its absorption or emission by the hyperfine transition with a
frequency $\nu_{0}=1420$ MHz at redshift $z$, such that $\nu=\nu_{0}/(1+z)$.
Therefore, the range of frequencies observed provides information from a given
epoch (although higher redshifts can also be reached employing Fourier modes
along the line of sight in scale dependent measurements [393]).

Depending on the task, either single stations or interferometric arrays of
antennas could be employed to detect radio frequencies. One of the key
differences between both approaches is the spatial resolution. In a single
antenna, usually designated as dish (although often they do not have that
shape), the angular resolution $\theta$ is determined by the wavelength
$\lambda$ of the observed radiation and the characteristic size of the antenna
$L$, as [393]

$\theta\simeq\frac{\lambda}{L}\,.$ (4.64)

On the other hand, instead of employing one unique dish, it is possible to
employ interferometric techniques to correlate many of them, acting as a huge
effective antenna which enlarges the collecting area. These telescopes are
known as interferometers, and are characterized, among other parameters, by
the number of antennas and the maximum baseline distance, which is roughly the
maximum length covered by the array. The observation of a given Fourier mode
can be performed by the correlation of a pair of antennas separated by a
baseline distance $b$. In this case, the angular resolution reads

$\theta\simeq\frac{\lambda}{b}\,.$ (4.65)

While single dish experiments usually have sizes of tens of meters (being the
Green Bank Telescope with 100 m the largest of the world [394]),
interferometers may spread over tens of kilometers, which allow them to reach
much larger resolutions. Selecting pairs of antennas with different baseline
distances, distinct Fourier modes could be probed. For this reason, in order
to observe the 21 cm fluctuations, interferometers are the proper choice.
However, these installations may be extremely expensive. Experiments committed
to measure the global average signal do not need such a high resolution, and
could work with single dishes. In the case of the interferometers, the maximum
baseline is of great relevance since it sets the maximum angular resolution
which can be achieved at a given wavelength. And equivalently, the minimum
baseline between two antennas determines the largest spatial scales which
interferometers are sensitive to. More explicitly, the wavenumber $k$ is
proportional to the baseline distance $b$ [355], and thus, inversely
proportional to the angular resolution $\theta$. In interferometers, the size
of antennas has also a great importance since it determines, together with the
number of elements, the total collecting area. The larger it is, the better
the sensitivity gets, and the fainter the signals that can be detected are.
Finally, from the equations above it can be noticed that the wavelength of
observation also determines the observable scale. As a consequence, low
frequency radiation, redshifted from older epochs such as the Dark Ages,
cannot be as well resolved as photons with larger frequencies.

#####  4.5.1 Foregrounds

Unlike the CMB, whose temperature is relatively free from foregrounds at the
observed frequencies, in the case of the 21 cm line they become the major
issue to overcome. These contaminants can be 3 or 4 orders of magnitude larger
than the cosmological brightness temperature, which makes understanding its
behavior mandatory in order to extract cosmological information.

There are two main sources of foregrounds from cosmic origins [355]: i)
galactic foregrounds, given by diffuse synchrotron and free-free emission from
the Milky Way (such as thermal bremsstrahlung of free electrons scattering off
free ions without being captured), and ii) extragalactic foregrounds, as radio
emission from Active Galactic Nuclei (AGN) and star-forming galaxies. While
the former are mostly important at large scales, being important for angular
scales larger than a degree, the latter dominates at smaller scales. The
synchrotron emission in the Milky Way is produced by relativistic cosmic-ray
electrons accelerated in supernovae remnants and other sources, which
typically present a power law energy spectrum. This implies that their
synchrotron emission has also a power-law spectral shape [238]. Since the
galactic plane contains most of stars, supernovae are more common there, where
cosmic rays originate. This leads to a stronger synchrotron emission along
this plane. On the other hand, AGN are powered by accretion of matter onto a
central supermassive BH, which can present strong radio emission. Star forming
galaxies can also present synchrotron and free-free emission in a similar way
to the Milky Way. In each case, power laws can describe their spectra.
Therefore, in both galactic and extragalactic foregrounds, given the
characteristic spectra of AGN, synchrotron and free-free emission, both
contributions are expected to be smoothly varying in frequency, which is a key
point for their removal. The characteristic decaying power-law shape provides
larger foregrounds at lower frequencies, which makes harder to detect signals
from distant epochs, such as the Dark Ages, whose photons have redshifted
more.

Employing scale-dependent statistics, such as the power spectrum, allows
removing the foregrounds much more easily than in the global signal case. This
is because, at some scales, gradients in the inhomogeneous signal could
present a sharp evolution. This is, for instance, what happens at the edges of
the ionized bubbles, where the signal varies abruptly along short distances.
These sudden changes contrast with the astrophysical foregrounds, which vary
much more smoothly, due to the synchrotron and free-free spectra, which can be
easily separated from the cosmological signal.

Besides the cosmic foregrounds, the presence of the ionosphere can also impact
the observed signal. It is specially important at the lowest radio
frequencies, becoming an even harder task to observe signals from distant
epochs, such as the Dark Ages. The most effective way to mitigate their effect
would be to observe from outside the Earth, as proposed in several telescopes
discussed in the following sections.

In interferometric searches, an alternative to foreground removal is
foreground avoidance, based on observing modes where foregrounds are
suppressed. The Fourier space can be decomposed in wavenumbers parallel and
perpendicular to the line of sight, $(k_{\|},k_{\bot})$. The frequency axis
corresponds to the line of sight direction. Thus, spectrally smooth
foregrounds would be mostly compacted in the lowest $k_{\|}$ modes. This
defines a region in the Fourier space, where foregrounds are restricted, and
then can be avoided by looking away. However, due to chromatic instrumental
effects which imprint unsmooth spectral features, the dominance of foregrounds
extends towards large $k_{\bot}$ at higher $k_{\|}$, shaping a wedge in
Fourier space [395]. Luckily, this region is relatively well delimited,
allowing to evade the contaminated modes. Foreground avoidance presents an
advantage since allows to elude possible biases from foregrounds fitting. The
drawback of this method is that the number of accessible scales becomes
significantly reduced. The subtleties of this procedure are discussed in,
e.g., Refs. [393, 396, 397].

#####  4.5.2 Global signal experiments

Attempting to measure the global averaged signal is one of the primary goals
of observational 21 cm cosmology. These measurements correspond to the
monopole in a spherical harmonics decomposition, typical in CMB studies.
Experiments designed to detect the global signature do not need to resolve
spatial details on the sky, and are thus commonly carried out by single
antennas that have wide beams, rather than large interferometric arrays, as
happens with experiments focused on measuring fluctuations (although there are
exceptions, such as LEDA, as discussed below). For the very same reason, these
experiments can also employ coarser frequency channels with larger bandwidths,
with no need for a spectral resolution as large as the experiments searching
for fluctuations. The mitigation of cosmic and ionosphere foregrounds, which
requires a precise calibration, becomes the greatest challenge. However, an
epoch of rapid reionization, with an abrupt drop in the signal, may be
different enough from the smoothly varying foregrounds, fact which has been
used to rule out too fast and sharp reionization processes, placing
constraints on the timing and duration of the EoR. The most relevant
experiments pursuing the global signature are listed in what follows (see,
e.g., Refs. [398, 399] for more details).

  * •

The Experiment to Detect the Global EoR Signature (EDGES) [400] is located at
Murchison Radio-astronomy Observatory in Western Australia, and composed by a
single dipole antenna. Early measurements were able to rule out very rapid EoR
[401]. The most remarkable result was announced in 2018, when its team
reported the evidence of an absorption profile in the low-band spectrum
(50-100 MHz), centered around 78 MHz (corresponding to $z\sim 17$), with an
amplitude of $\sim 0.5$ K [36]. This deep trough is roughly twice the maximum
predicted within the context of the $\Lambda$CDM model. The implications and
caveats of this measurement are discussed in Sec. 4.5.3. Regardless of this
absorption feature, high-band measurements spanning 90-190 MHz [402] have been
employed to constrain the astrophysical parameter space from Cosmic Dawn to
the EoR [403, 404, 405].

  * •

The Shaped Antenna to measure the background RAdio Spectrum (SARAS) [406] was
sited at the Raman Research Institute in India, although it has been its
upgraded version, SARAS 2 [407] at Timbaktu Collective in Southern India, the
one which have offered more competitive results. Using its observations at
frequencies 110-200 MHz, rapid periods of reionization have been constrained
[408]. Furthermore, SARAS 2 has ruled out several astrophysical models which
present a very cold IGM at the epoch of the Cosmic Dawn, which would lead to
large amplitudes in the global signal [408, 409]. Note that this may be in
tension with the EDGES results.

  * •

The Large-Aperture Experiment to Detect the Dark Ages (LEDA) [410, 411] is
placed on the Long Wavelength Array stations at Owens Valley Radio Observatory
(OVRO-LWA), in California, USA. Unlike other global signal single-antenna
experiments, it is formed by $\sim 250$ dipole antennas spread over an area of
$\sim 200$ m of diameter, constituting an interferometric array capable of
observing frequencies between 30 and 88 MHz, which may probe the Cosmic Dawn
epoch (redshifts $16<z<34$). LEDA observations allowed to place limits on the
amplitude and width of the 21 cm absorption profile at $z\sim 20$ [412].

  * •

The Dark Ages Polarimeter Pathfinder (DAPPER) [413, 414] is an ambitious space
based project to put into orbit a spacecraft in the dark side of the Moon,
avoiding therefore the ionosphere foregrounds. Unlike the previous
experiments, it is at the planning stage and has not been built yet. It is
designed to operate at very low frequencies, between 10-110 MHz, being able in
principle to reach the Dark Ages, at redshifts $z\sim 36-83$. It was preceded
by the proposal of the Dark Ages Radio Explorer (DARE) [415, 416].

#####  4.5.3 The EDGES affair

Figure 4.5: Global 21 cm brightness temperature for the best fit to explain
the EDGES absorption feature shown in red, with its 99% CL interval [36].
EDGES low-band data lie in the range covered by the green region. For
comparison, an extreme case with adiabatic cooling, $T_{S}$ completely coupled
to $T_{K}$ and instantaneous reionization at $z\sim 10$ is also included in
blue, as well as a fiducial model shown in purple computed from 21cmFAST [378]
such that it presents the absorption dip at roughly the same redshifts (see
Chapter 5 for a detailed discussion on the evolution of the brightness
temperature and the astrophysical processes involved).

The reported absorption profile by the EDGES collaboration in the low-band
spectrum (50-100 MHz) has supposed a major result in the last years, given its
unexpected shape. In order to explain such a deep absorption trough, there are
two main possibilities, as can be seen from equation Eq. (4.21): either i) the
spin temperature (and thus the kinetic one) is very cold, or ii) the CMB
background $T_{\gamma}$ has extra contributions at radio frequencies.

Standard evolution histories within the $\Lambda$CDM scenario limit the
coldest kinetic temperature attainable, since the only cooling process
efficient at this regime is due to adiabatic expansion, $T_{K}\propto a^{-2}$.
Figure 4.5 shows a sample of such extremely optimistic scenario in blue, with
the largest possible amplitude within the standard picture, compared to the
EDGES result in red. Even with the spin temperature completely coupled to the
adiabatic one at these epochs, the largest possible signals could not exceed
$\sim-250$ mK at $z\sim 17$. The issue gets worse for more realistic
scenarios, as the one shown in purple, where heating of the IGM and the non-
perfect coupling of the spin temperature to the kinetic one have been taken
into account properly (as shall be reviewed in Chapter 5). The EDGES
measurement may thus imply that extra cooling processes are present during the
Cosmic Dawn and Dark Ages to drive the spin temperature to such a low value.
Following this idea, many authors have proposed non-standard scenarios to
achieve a cooler IGM. This is the case of DM-baryon interactions, which would
transfer thermal energy from the gas to the DM particles [417, 418, 419, 420,
421, 422]. An early Dark Energy period could provide an earlier decoupling of
the gas temperature from the Compton interactions, giving more time to cool
adiabatically [423].

The excess of radio radiation at the epoch of the Cosmic Dawn has also been
suggested as a possible way to increase $T_{\gamma}$ and thus explain the
EDGES signal [424]. However, it has been argued that the diffuse radio
background from astrophysical origins, by either synchrotron or inverse-
Compton scattering due to relativistic electrons, seems unlikely to explain
the signal. This is because the radio emissivity required to enhance
$T_{\gamma}$ should be enlarged by a factor from $10^{3}$ to $10^{6}$ respect
to the local $z=0$ estimates [425, 426]. However, the existence of soft photon
emission at low radio frequencies from light DM, such as axion-like particles,
may account for the required background [427, 428, 429]. Accretion onto
intermediate-mass black holes has also been proposed to produce the required
extra radio emission [430, 431].

Besides trying to explain the signal, the EDGES result has also been used to
constrain different DM models, such as DM annihilation [432, 433], DM decay
[434], WDM [435, 436], PBHs [350, 351, 352] or FDM [437, 436]. In Part II of
this thesis, an example of this is shown, constraining IDM and WDM models
making use of the timing of the EDGES signal [5]. The consistency of the EDGES
result with other experimental constraints from EoR and CMB have also been
explored, even though that work is not included in this thesis [8].

Thus, although some (exotic) solutions may be able to explain the EDGES
result, it also presents some controversy within the scientific community for
other reasons. It has not been confirmed yet by other experiments.
Observations of SARAS 2 at similar range of frequencies have ruled out models
which should show such a deep absorption trough, being in tension with the
EDGES finding [408, 409]. The shape of the signal extracted from the data is
also highly non-trivial and not well understood. The absorption model to fit
the data employed in [36] was chosen to be a flattened gaussian, with no clear
interpretation from physical grounds. The broad dip suggests a relatively
extended epoch around $z\sim 17$ with constant temperature, which can not be
explained within current IGM evolutionary models. In fact, when one attempts
to accommodate the astrophysical modeling of the signal to the EDGES data, the
width of the absorption dip can only be fitted in models with high heating,
opposite to what is required to obtain large amplitudes, suggesting an
additional inconsistency between the measurement and canonical astrophysical
scenarios [8].

In addition, several caveats regarding the treatment of foregrounds and
calibration have been claimed, which could invalidate the measurement.
Although the EDGES calibration has been carefully scrutinized by the team
during 10 years, some concerns about possible systematics have been raised. A
possible ground plane artifact may produce such a broad absorption signal, as
was shown by Ref. [438]. This would induce one or several resonances in the
data, which could be confused with a cosmological signal. The authors found
that three of these resonant features with no 21 cm signal provided a fit as
good as the one obtained by EDGES. Remarkably, the best-fit resonant
frequencies were in good agreement with those estimated for a ground plane
artifact, considering the dimensions of the square patch and properties of the
soil.

Furthermore, Ref. [439] found that different foreground and absorption models
could fit the EDGES signal. It was shown that the best fit parameters for the
ionospheric foregrounds employed by EDGES were not physical, since they imply
optical depths and temperatures of the ionosphere with negative values. In
addition, it was found that two gaussians or a sinusoidal feature with no
absorption signal can provide a fit to the spectrum as good as the one
presented by the EDGES collaboration. These non-cosmological features could be
a consequence of some sort of overlooked systematics, such as the one proposed
by [438].

The authors of Ref. [440] further analyzed these questions, by employing a
Bayesian evidence-based comparison over a broad amount of models. They
discussed different models for the signal, including physical motivated ones
from simulations varying astrophysical parameters. Cosmic and ionospheric
foregrounds, as well as possible calibration errors with sinusoidal shapes,
were also considered. It was found that models including calibration
systematics were decisively preferred. Moreover, models with amplitudes
consistent with the standard cosmological scenario were present among the best
fitting models. Furthermore, they did not find strong evidence to favor models
exhibiting 21 cm signal over those without one. All these facts place
reasonable doubts about whether it is possible to infer an actual 21 cm
cosmological measurement from the EDGES data.

All the above references show the great impact that the EDGES result have left
on the astrophysics and cosmology communities. However, the aforementioned
caveats place important concerns regarding the validity of this measurement,
which can only be confirmed by upcoming observations from future
radiotelescopes.

#####  4.5.4 21 cm power spectrum experiments

Facility | Location |  $\nu$ [MHz] | $z$ | $N_{\rm ant}$ |  $b_{\rm max}$ [km]
---|---|---|---|---|---
GMRT | India | 50-1420 | 0-27 | 30 | 30
MWA | Australia | 70-90, 135-195 | 15-19, 6-10 | 128 | 5
LOFAR | Netherlands | 30-80, 120-190 | 17-46, 6-11 | 50–60 | 50
PAPER | South Africa | 110-180 | 7–12 | 64 | 0.210
OVRO-LWA | USA | 27-88 | 15-52 | 288 |  $<$ 1.5
LEDA | USA | 45-88 | 15-30 | 256 |  $<$ 10
HERA* | South Africa | 50-250 | 5-27 | 350 |  $<$ 1
SKA-low* | Australia | 50-350 | 4-27 | $\sim 10^{5}$ | 65
SKA-mid* | South Africa | 350-15.3$\times 10^{3}$ | 0-4 | 197 | 150

Table 4.1: List of current and future interferometers attempting to measure
the spatial fluctuations of the 21 cm signal, including its location, range of
frequencies of operation $\nu$, accessible redshifts $z$, number of antennas
$N_{\rm ant}$, and maximum baseline attainable $b_{\rm max}$. Experiments
marked with an asterisk (*) are not completed yet. Note that LEDA, despite
being an interferometer, is devoted to detect the global signal rather than 21
cm fluctuations.

Figure 4.6: Summary of upper limits on the 21 cm power spectrum. The dots
denote the bounds found by the experiments GMRT [441, 442], MWA [397, 443,
444, 445, 446, 447, 448], LOFAR [449, 450, 451], OVRO-LWA [452, 453], PAPER
[454] and AARTFAAC (an extension to LOFAR) [455]. Most of them lie between the
wavenumbers $k=0.1$ and $0.5$ $h/$Mpc, except the ones from Refs. [442] and
[450] (note also that Ref. [449] presents stronger limits at $k=0.053$
$h/$Mpc). The continuous black line corresponds to a fiducial theoretical
model computed with 21cmFAST [378], at $k=0.2$ $h/$Mpc (details regarding its
evolution and the astrophysical modeling are discussed in Chapter 5). Gray
bands depict the forecasted sensitivities for the HERA and SKA experiments,
computed via the code 21cmSense1313footnotemark: 13[456]. The script to
generate this figure has been made publicly available1414footnotemark:
14[457].

Although detecting the global signature is of significative importance, the
information which can be extracted from it is limited, and the experimental
challenges and caveats are numerous. For this reason, most observational
efforts are devoted to detect the spatial fluctuations of the 21 cm signal,
which would enable a multi-scale view of the Hydrogen map, and an easier
removal of the foregrounds. Interferometry based on large arrays of antennas
conforms the proper tool to access small scales. Interferometers offer a
natural way to directly estimate the power spectrum of an astronomical image.
This is because the voltage induced in a pair of antennas can be correlated,
providing the visibility, which is nothing but the Fourier transform of the
perceived intensity of the radiation [458, 355]. Moreover, as already
mentioned, foregrounds can be more easily treated than when dealing with
global measurements.

The achieved sensitivity is subject to different aspects of the experimental
design and systematics. It is also limited by the spatial scales and redshifts
to be observed. Among the possible sources of noise, the most important one is
thermal noise. Thermal fluctuations in the apparatus produce a gaussian white-
noise signal (i.e., scale independent) characterized by its root-mean-square
$T_{N,rms}$ [458]. Thus, the thermal noise at a single baseline goes as
$\overline{\delta T_{b}}^{2}(z)\Delta^{2}_{21}(k,z)\sim k^{3}T_{N,rms}^{2}$.
This implies that the sensitivity worsens at very small scales. Furthermore,
due to the factor to convert from angles and bandwidth to cosmological
distances, this noise also increases with redshift, roughly as $(1+z)^{5/2}$
[459]. Therefore, distant sources (i.e., remote epochs) become more difficult
to observe by a limitation in the sensitivity, and not only by the effect of
foregrounds. In order to calculate the sensitivity of a full interferometric
array, one must sum up the contributions by each baseline over all the Fourier
modes. Combining independent $k$ modes, the sensitivity can be improved, with
a noise dependence of $\sim k^{5/2}$ rather than $k^{3}$. Moreover, if
different baselines measure the same Fourier modes, having redundant
observations, the total noise decreases with the number of antennas as $N_{\rm
ant}^{-1}$. This dependence is easy to interpret, given that the power
spectrum is a sort of variance, which typically decreases with $N_{\rm
ant}^{-1}$ when independent measurements are taken into account. It is thus
auspicious to have interferometers with a large number of dishes, placed at
regular grids to probe the same Fourier modes, in order to improve the
sensitivity, see Refs. [458, 459] for more details on these regards.

Thus far, none of the existing experiments has achieved a detection of the 21
cm fluctuations from the pre-EoR universe (although some cross-correlation
with other probes have been detected from the post-EoR low redshift universe,
as discussed below). However, current radiointerferometers have placed
remarkable upper limits on the amplitude of the fluctuations, increasingly
tighter over the years, which have constrained the possible thermal evolution
histories. The current upper bounds are shown in Fig. 14, together with a
fiducial model of the 21 cm power spectrum. In the following, we comment on
the interferometer experiments which have been able to place constraints on
the 21 cm fluctuations amplitude from epochs between the Cosmic Dawn and the
EoR, which are also listed in Tab. 4.1. For more information, see, e.g., Refs.
[393, 355, 398, 460].

  * •

The Giant Metrewave Radio Telescope (GMRT) [461] is an array located in
Maharashtra, India, for low-frequency general purposes. It consists of 30
dishes of 45 m of diameter, spread over 25 km. It allows observing in 6
frequency bands, covering from 50 to 1420 MHz. Reionization and post-EoR
epochs have been studied with this instrument, being pioneer in the treatment
of foregrounds and their removal employing spatial correlations. Although not
optimized for the cosmological 21 cm signal, it has been the first one to
place upper limits at redshifts relevant for the EoR, at $z\sim 9$ [462],
being subsequently revised (and significantly weakened) [441]. A stringent
constraint has also been recently proposed at $z\simeq 8$ [442].

  * •

The Murchison Widefield Array (MWA) [37] is located in the Murchison Radio-
astronomy Observatory, at the Western Australian desert. It consisted of 128
elements in its Phase I (2013-2016), being increased to 256 with longer
baselines since 2016. Each of these elements is formed by a tile of $4\times
4$ grid of crossed dipole antennas, ranging from 70 to 300 MHz. The MWA has
published several upper limits covering a broad range of epochs, from the EoR,
at redshifts $z\sim 6$ to $11$ [397, 443, 444, 446, 447, 448], to Cosmic Dawn
[445], constraining redshifts $z\sim 12,15$ and $17$.

  * •

The Donald C. Backer Precision Array for Probing the Epoch of Reionization
(PAPER) [463, 464] is an interferometer designed and optimized for
observations of the EoR, operating between 100 and 200 MHz. It consisted of
two stations, a preliminary one with 32 antennas sited at Green Bank, West
Virginia, USA, for testing, and the main installations for scientific research
with 64 antennas placed at the South African Karoo desert. Its team pioneered
ways to improve the sensitivity, proposing methods such as maximizing
redundant observations of the same Fourier modes to reduce the noise [458],
and the so-called delay spectrum technique to avoid foregrounds [465]. It
provided several upper limits at the EoR [466, 467, 468]. However, posterior
re-analysis of their data showed errors in their analysis related to the
covariance matrices and to an underestimation of the level of signal loss,
which invalidated these limits [469, 470]. Nevertheless, they were later
updated properly accounting for these errors, significantly weakening the
bounds [454] (see [355] for further details on this issue). The PAPER
experiment is already deccomissioned, but its infrastructure in South Africa
currently serves for the HERA interferometer.

  * •

The LOw Frequency Array (LOFAR) [38] is a multi-purpose interferometer focused
on low frequencies sited at the Netherlands151515The Netherlands has a long
tradition related to the 21 cm signal. Among other achievements, it was Dutch
astronomers, van de Hulst and Wouthuysen, who predicted the hyperfine
transition and the Ly$\alpha$ coupling, respectively, becoming aware of the
astronomical importance of the 21 cm signal in radioastronomy. Furthermore,
Muller, Oort and Van de Hulst observed the first rotation curve and HI map of
the Milky Way, showing for the first time its spiral arms. See Ref. [471] for
more information.. It includes two types of antennas, covering high-band
(110-250 MHz) and low-band (10-90 MHz) observations. While the low-band
elements are standard dipole antennas, the high-band ones consist of a tile of
16 individual antennas. Upper limits on the 21 cm power spectrum have been
placed using its observations, both at the EoR, covering redshifts $z\sim 8$
to $10$ [449, 451], and at the Cosmic Dawn, reaching redshifts from $z\sim 19$
to $25$ [450]. An add-on real time transient detector facility appended to
LOFAR called Amsterdam-ASTRON Radio Transients Facility And Analysis Center
(AARTFAAC) [472] has recently found a stringent upper limit at $z\simeq 18$
[455]. An extension of LOFAR coined NenUFAR, located in Nançay, France, is
ongoing, which may serve as a pathfinder for future experiments such as SKA
[473].

  * •

The Owens Valley Radio Observatory Long Wavelength Array (OVRO-LWA) is sited
in Owens Valley, California, USA. It is formed by 288 dipole antennas spanning
from 27 to 85 MHz, arranged in an pseudo-random disposition in a compact core
of 200 m of diameter, being thus well suited for transient searches such as
gamma ray bursts or exoplanets. Besides that, it has also been employed for
constraining the 21 cm fluctuations, finding an upper limit at $z\simeq 18$
[452], and recently in the range $25<z<31$, being this one the bound located
at the farthest redshift so far [453]. The LEDA facilities are also located at
the OVRO-LWA.

All of the aforementioned upper limits on the amplitude of fluctuations are
placed at a range of scales spanning $0.1\,h{\rm Mpc}^{-1}<k<1\,h{\rm
Mpc}^{-1}$ (although the LOFAR ones include also modes down to $k\sim 0.01$
$h{\rm Mpc}^{-1}$), and are about 2 orders of magnitude larger than the
expected fiducial models for the 21 cm signal [393]. Besides the
interferometers mentioned above, some single dish experiments have been
employed to constrain the 21 cm fluctuations at low redshift, from the post-
reionization universe. Although these epochs lie beyond the scope of this
thesis, it is worth mentioning some results due to their significance. The
Green Bank Telescope (GBT) [394], located at Green Bank, West Virginia, USA,
is the world’s largest fully steerable radio telescope, with 100-meter
diameter. Operating at frequencies between 100 and 115 GHz, it has been the
first to detect 21 cm fluctuations in the post-EoR universe, via cross-
correlations with other surveys: with the optical galaxy survey DEEP2 at
$z\simeq 0.5-1.1$ [474], and WiggleZ at $z\sim 0.8$ [475]. The Parkes Radio
Telescope [476] in New South Wales, Australia, has 64-meter diameter, and
detects frequencies from 1230 to 1530 MHz. Similarly to the GBT, the cross-
correlation of its observations with those from the 2dF galaxy survey provided
detections of the cosmological HI line at redshifts $0.057<z<0.098$ [477].

On the other hand, several radiointerferometers are expected to be built
within this decade, which may have enough sensitivity to attain a positive
detection of the signal in the next few years. Some of them are devoted to
perform HI line intensity mapping, i.e., measuring the integrated 21 cm
emission from unresolved gas clouds, accessing clumpy neutral regions in the
post-reionization era, up to $z\sim 2.5$. Examples of these are the Canadian
CHIME [478], HIRAX in South Africa [479], the Tianlai experiment in China
[480] or the BINGO observatory in South America [481]. Among the proposals for
studying the Cosmic Dawn and EoR periods, here we highlight two
interferometers which by their sensitivity and international effort, may
become a great step forward towards the detection of the HI signature.

  * •

The Hydrogen Epoch of Reionization Array (HERA) [39] is placed at the same
facilities than PAPER in the Karoo desert, South Africa, specifically designed
to observe the EoR. For that, it will observe in a range between 50 and 250
MHz, and is expected to have 350 dishes with 14 m diameter each, when
completed. It has been deployed following the ideas developed for PAPER, with
a highly redundant array to enhance the sensitivity. Although not finished
yet, it is already running and collecting data.

  * •

The Square Kilometer Array (SKA) [40], is the most ambitious
radiointerferometer designed so far, with the world’s largest area ($\sim$ 1
km2) when completed. It will consist of two different facilities operating at
different frequencies. On the one hand, SKA-low will be comprised of $\sim
10^{5}$ antennas grouped into stations in the Murchison Radio-astronomy
Observatory at the Western Australian desert, spanning frequencies from 50 to
350 MHz. The SKA-mid, on the other hand, will be formed by 197 dishes in the
South African Karoo desert, observing frequencies from 350 MHz to 15.3 GHz.
The pseudo-random arrangement of antennas will allow reaching baselines up to
150 and 65 km, respectively. It will employ the facilities previously occupied
by MWA and HERA. Although far from being finished yet, pathfinder telescopes
are already constructed and gathering data, such as MeerKAT in South Africa
[482] and ASKAP in Australia [483].

  * •

The Farside Array for Radio Science Investigations of the Dark ages and
Exoplanets (FARSIDE) is a proposal for building an interferometric array based
on the Moon [484], in addition to the orbiting satellite DAPPER. It would be
composed by 128 dipoles across a 10 km2 area. Covering from $0.1$ to 40 MHz,
it is aimed to detect the faint 21 cm power spectrum during the Dark Ages.
Contrary to previous interferometers, it is still an idea under design and its
construction has not started yet.

The above lists of current and future interferometers are not intended to be
comprehensive, but to provide an outlook of experimental status, as well as
about the increasing efforts of the scientific community for pursuing the
cosmological 21 cm signal. A lot of work remains to be done, but these are
exciting times for seeking the hyperfine line of the IGM.

###  Chapter 5 Evolution of the Intergalactic Medium

The physics of the IGM is highly rich and complex. Its study is fundamental in
order to understand how and when galaxies formed, evolving from a mostly
homogeneous fluid towards the cosmic web of galaxies and quasars immersed in
an ionized medium we can see today. At $z\sim 1100$, the CMB decouples from
the thermal plasma, marking the beginning of the so-called Dark Ages. After
Recombination, most atoms become neutral, only a small fraction of free
electrons remaining. The gas temperature keeps cooling, coupled to the photon
temperature due to Compton scattering until $z\sim 150$, when these
interactions become inefficient. Thereafter, the gas cools adiabatically.
During the Dark Ages, overdensities around primordial seeds keep growing,
leading eventually to the formation of galaxies. The newly born stars emit
highly energetic radiation capable of altering the IGM evolution.

In order to understand the IGM evolution in detail, full numerical simulations
are required, properly computing the hydrodynamics and radiative transfer,
which are very expensive computationally wise. However, some semi-analytical
approximations can be assumed, which simplify significantly the picture, and
allow us to gain insight into the physical processes involved. Simplifying the
full picture, we can treat the IGM as a two-phase medium [353]. One of the
phases is the mostly neutral IGM, characterized by the gas temperature
$T_{K}({\bf x},z)$, and the free electron fraction $x_{e}({\bf x},z)$, which
accounts for the free electrons remaining from Recombination. After the onset
of star formation, energetic X-ray radiation can reach these neutral areas,
heating and ionizing the gas, although still remaining mostly neutral. On the
other hand, the immediately surrounding regions of galaxies and ionizing
sources are expected to be almost completely ionized, forming HII bubbles.
Once the star formation processes are efficient enough, these bubbles can grow
and eventually merge among them, ionizing the full IGM during the so called
Epoch of Reionization (EoR). We characterize the fraction of volume occupied
by the HII bubbles with the ionized filling factor $Q_{i}$, where the local
ionization fraction is $\sim 1$. In these completely ionized regions,
temperatures reach $\sim 10^{4}$ K, where atomic cooling peaks. Therefore, the
global ionization fraction including the (mostly) neutral IGM plus the HII
regions can be written as

$\bar{x}_{i}=Q_{i}+(1-Q_{i})\bar{x}_{e},$

where $\bar{x}_{e}$ is the spatial average over $x_{e}({\bf x},z)$ (in the
neutral IGM). Hence, during the EoR, the mostly neutral IGM becomes completely
ionized. The two-phase IGM is a simplification of the real IGM, but accurate
enough and widely employed in semi-analytical computations, as in the code
21cmFAST [378, 485], which is widely employed throughout this thesis. In this
chapter, we review the evolution of the ionized fraction and gas temperature
at each phase, discussing the different sources of energy injection, as well
as the current constraints on the reionization epoch. For further details,
there are many reviews in the literature treating the physics of reionization
and the IGM, see Refs. [52, 486, 487, 488, 57, 58].

####  5.1 Ionized IGM

In this section, we discuss the details of the reionization processes,
starting by considering only the globally averaged ionized fraction, and
commenting latter on the effects of inhomogeneities. Current observational
constraints on the EoR are summarized at the end of the section.

#####  5.1.1 Global Reionization

The evolution of the HII number density $n_{\rm HII}({\bf x},t)$ is given by
its continuity equation (e.g., [487]):

$\frac{dn_{\rm HII}}{dt}+3Hn_{\rm HII}=\Gamma_{\rm HI}n_{\rm
HI}-\alpha_{A}n_{\rm HII}n_{e},$ (5.1)

where the left-hand side stands for the adiabatic evolution and the right-hand
side accounts for the interplay between ionizations and recombinations, being
$\Gamma_{\rm HI}$ the ionization rate and $\alpha_{A}$ the so-called case-A
recombination coefficient. Peculiar velocities may also be taken into account,
but we neglect them in this derivation for the sake of simplicity. The
recombination coefficient considered depends upon the medium. In neutral IGM,
ionizing photons are rapidly absorbed, and transitions to the ground state do
not effectively contribute to the ionization. This is the so-called case-B
recombination, where recombinations to the ground state are excluded.
Contrarily, the HII phase is a thin medium where these UV photons can freely
travel through. For this reason, employing case-A recombination there is more
adequate, such that recombinations to the lowest state are taken into account
[354]. In this case, the coefficient takes the value $\alpha_{A}=4.2\times
10^{-13}$cm3 s-1, evaluated at $T_{K}=10^{4}$ K, a factor of $\sim 2$ larger
than case-B.

We characterize the state of ionization of the IGM by means of the fraction of
ionized volume, defining the filling factor as the volume-average of the HII
fraction in the ionized medium: $Q_{i}=\langle x_{\rm HII}\rangle$, where
$x_{\rm HII}=n_{\rm HII}/n_{\rm H}$. This quantity provides the fraction of
volume of the IGM which is completely ionized by UV photons. To obtain an
evolution equation for $Q_{i}$, we must spatially average the above equation.
Since in this phase, both Hydrogen and Helium are mostly ionized,
$n_{e}\simeq(1+\chi)n_{\rm HII}$, with $\chi=\mathfrak{f}_{\rm
He}/\mathfrak{f}_{\rm H}$, with $\mathfrak{f}_{\rm H}=n_{\rm H}/n_{b}=0.92$
and $\mathfrak{f}_{\rm He}=n_{\rm He}/n_{b}=0.08$ the fraction of Hydrogen and
Helium atoms per baryon, respectively.111Doubly-ionized Helium, HeIII, can be
usually neglected for the redshifts of interest, due to its high ionization
threshold. Decomposing the HII number density as $n_{\rm HII}=x_{\rm
HII}\,\bar{n}_{\rm H}\,(1+\delta)$, with $\delta$ the density contrast respect
to the mean value, as defined in Sec. 2.2.1, the squared averaged HII density
can be written as

$\langle n_{\rm HII}^{2}\rangle\simeq C\,Q_{i}\,\bar{n}_{\rm H}^{2},$ (5.2)

where we have approximated $\langle x_{\rm HII}^{2}\rangle\simeq\langle x_{\rm
HII}\rangle=Q_{i}$, given that $x_{\rm HII}$ takes values very close to 1 in
these mostly ionized regions. The clumping factor $C=\langle n_{\rm
H}^{2}\rangle/\langle n_{\rm H}\rangle^{2}=\langle(1+\delta)^{2}\rangle$
accounts for the fact that the IGM is not completely homogeneous and has a
density distribution. These inhomogeneities have an impact in the ionization
history through this clumping factor, since recombinations are enhanced in
overdense regions. There are several prescriptions to compute this
coefficient, depending on the specific definition and method of simulation,
and it depends upon the scale at which it is evaluated [489, 490]. While the
modeling of the clumping factor is highly challenging, most estimations
predict values of order unity at the redshifts of interest [491, 492].

We focus now on the ionizing term of Eq. (5.1). As long as the mean free path
is much smaller than the size of the universe, or equivalently, the opacity
$\kappa$ is larger than the expansion rate, $1/\kappa\ll c/H$, ionizing
photons are absorbed soon after being released. This is the so-called on-the-
spot approximation, which can be safely used for the ionizing field at the
EoR. Employing $J$ and $\varepsilon$ as the flux and comoving emissivity by
number, respectively,222Note that in Sec. 4.1, $I$ and $\epsilon$ were the
flux and emissivity by energy, rather than by number, i.e., with an additional
factor $h\nu$. this simplification implies that Eq. (4.1) can be replaced by

$0\simeq-c\kappa J+\frac{c(1+z)^{3}\varepsilon}{4\pi}.$ (5.3)

Hence, the ionizing UV flux, $J_{\rm UV}(\nu)$ (in units of s-1 keV-1 cm-2
sr-1), reads

$J_{\rm UV}(\nu)\simeq\frac{(1+z)^{3}\varepsilon(\nu)}{4\pi\kappa(\nu)}$ (5.4)

with $\kappa(\nu)=\sigma_{\rm HI}(\nu)n_{\rm HI}$ and $\sigma_{\rm HI}(\nu)$
the photoionization cross section for Hydrogen. Introducing the above flux in
the source term of Eq. (5.1), one obtains:

$n_{\rm HI}\Gamma_{\rm HI}=4\pi n_{\rm HI}\int^{\infty}_{\nu_{\rm
HI}}d\nu\,\sigma_{\rm HI}(\nu)J_{\rm
UV}(\nu)\simeq(1+z)^{3}\int^{\infty}_{\nu_{\rm
HI}}d\nu\,\varepsilon(\nu)=\dot{n}_{\gamma},$ (5.5)

with $\dot{n}_{\gamma}$ the proper rate density of emitted photons and
$h\nu_{\rm HI}=13.6$ eV the threshold ionization energy for Hydrogen. Thus,
spatially averaging Eq. (5.1), and making use of the results from Eqs. (5.2)
and (5.5), we obtain the evolution equation for the filling factor,

$\frac{dQ_{i}}{dt}=\frac{\dot{n}_{\gamma}}{\bar{n}_{\rm
H}}-\frac{Q_{i}}{\bar{t}_{rec}},$ (5.6)

being $\bar{t}_{rec}^{-1}=\alpha_{A}(1+\chi)C\bar{n}_{\rm H}$ the mean
recombination time. The source term $\dot{n}_{\gamma}$ is completely
determined by the emissivity, which can be modeled as proportional to the
comoving star formation rate $\dot{\rho}_{*}$,

$\varepsilon(\nu)=\mathcal{N}(\nu)f_{esc}\frac{\dot{\rho}_{*}}{\mu
m_{p}}=\mathcal{N}(\nu)f_{*}f_{esc}\bar{n}_{b,0}\frac{df_{coll}(z,M_{*})}{dt},$
(5.7)

with $f_{*}$ the fraction of baryon into stars, $f_{esc}$ the fraction of
photons which escape from the galaxies to the IGM, $\mu$ the mean baryonic
weight and $\mathcal{N}(\nu)$ the the number of emitted photons per baryon per
frequency unit, which accounts for the energy spectrum and it is usually
assumed to be a power-law function. In the second equality we have written the
star formation rate $\dot{\rho}_{*}$ in terms of the derivative of the the
fraction of mass collapsed in halos above a mass $M_{*}$, $f_{coll}$, defined
in Eq. (2.23). Here, since we assume that ionizing photons come from stars
within galaxies, $M_{*}(z)$ should be taken as the minimum mass required to
host star formation in a halo. Therefore, only halos with masses large enough
to host star formation can contribute to the UV flux. Using the Press-
Schechter prescription, Eq. (2.21), the fraction of collapsed mass can be
computed analytically, obtaining Eq. (2.24), while for the Sheth-Tormen
function, Eq. (2.22), it must be computed numerically. This threshold can be
related to the minimum virial mass $M_{\rm vir}^{\rm min}$ corresponding to a
virial temperature $T_{\rm vir}^{\rm min}$, which can be defined from Eq.
(2.17) as [52]

$M_{\rm vir}^{\rm min}(z)=10^{8}\left(\frac{T_{\rm vir}^{\rm min}}{1.98\times
10^{4}\,{\rm
K}}\frac{0.6}{\mu}\right)^{3/2}\left(\frac{1+z}{10}\right)^{-3/2}M_{\odot}/h\leavevmode\nobreak\
,$ (5.8)

where $\mu$ is the mean molecular weight which is equal to $1.2$ ($0.6$) for a
neutral (fully ionized) primordial gas. The minimum virial temperature is
usually chosen at the atomic cooling threshold, $T_{\rm vir}^{\rm min}\simeq
10^{4}$ K, minimum value required to have an efficient atomic cooling,
allowing the gas to fragment, condense and form stars. However, molecular
cooling would be relevant in metal-free clouds, cradle of the so-called
Population III stars, a hypothetical first generation of very massive stars
lacking metals. The formation of H2 would present the most important mechanism
to lose thermal energy, leading to lower values of the temperature, down to
$\sim 10^{3}$ K [52, 58]. Given the uncertainty on the stellar population and
cooling mechanism in halos, in the following we consider either $T_{\rm
vir}^{\rm min}$ or $M_{\rm vir}^{\rm min}$ as a free parameter. The density
rate of ionizing photons is then

$\dot{n}_{\gamma}(z)=N_{\gamma/b}f_{*}f_{esc}\bar{n}_{b}(z)\frac{df_{coll}(z,M_{\rm
vir}^{\rm min})}{dt},$ (5.9)

being $N_{\gamma/b}=\int_{\nu_{\rm HI}}d\nu\mathcal{N}(\nu)$ the number of
ionizing photons per baryon. Note that this formalism is specially appropriate
for high redshifts, since at later times, $z<10$, mergers of halos can also
contribute to star formation [353]. Eq. (5.6) can thus be rewritten as

$\frac{dQ_{i}}{dt}=\xi\frac{df_{coll}}{dt}-\frac{Q_{i}}{\bar{t}_{rec}},$
(5.10)

being $\xi$ the _ionization efficiency_ , defined as the product of the
astrophysical parameters

$\xi=N_{\gamma/b}f_{*}f_{esc}\frac{1}{\mathfrak{f}_{\rm H}}.$ (5.11)

This is the final equation that must be solved in order to predict the
evolution of the ionized medium during the reionization era.333Eq. (5.10) can
also apply to the post-reionization era by accounting for the residual neutral
clouds in the ionized IGM as an extra term in the opacity, correctly
reproducing the redshift evolution [493]. A further simplification can be
taken when recombinations are slow enough. Assuming that the product
$N_{\gamma/b}f_{*}f_{esc}$ is redshift independent, we can formally solve Eq.
(5.10) as [494]

$Q_{i}(z)=\frac{\xi\;f_{coll}(z,M_{\rm vir}^{\rm min})}{1+N_{rec}},$ (5.12)

where
$N_{rec}=Q_{i}(z)^{-1}\int_{z}(dt/dz^{\prime})dz^{\prime}Q_{i}(z^{\prime})\bar{t}_{rec}^{-1}(z^{\prime})$
is the mean cumulative number of recombinations per Hydrogen atom. The filling
factor is therefore proportional to the fraction of matter collapsed in star-
forming halos, with an efficiency factor given by a product of astrophysical
parameters. In principle, these parameters may be redshift and mass dependent,
but due to the uncertainties in their values and for the sake of simplicity,
we shall consider them as constant both in mass and in redshift. The number of
ionizing photons per baryon $N_{\gamma/b}$ can be estimated from the spectrum
and the initial mass function of stars, which depend upon the details of star
formation (burst formation or continuous) and on their metallicity [495, 496].
Making use of the present-day initial mass function of low-metallicity stars,
known as Population II stars, $N_{\gamma/b}\sim 4\times 10^{3}$. Nonetheless,
it can reach $\sim 10^{5}$ if metal-free and very massive stars of
$>10^{2}M_{\odot}$, the so-called Population III, are considered as the
initial composition. Regarding the fraction of baryons into stars, previous
works on the comparison of the star formation rate with the one derived from
UV luminosity function measurements (see, e.g., Refs. [497, 498] and
references therein), as well as radiation-hydrodynamic simulations of high-
redshift galaxies (see, e.g., Ref. [499]), have found values of this quantity
around $f_{*}\sim 0.01$. On the other hand, most photons produced in a galaxy
would remain within, being consumed in ionizing its neutral matter.
Nonetheless, a small fraction of photons, encoded by $f_{esc}$, could escape
from the galactic boundaries into the IGM and contribute to the reionization
mechanism. The estimation of this parameter is challenging since it depends on
the distribution of column densities and galactic morphologies. Although it is
neither well observed nor well modeled, it takes typical values of
$f_{esc}\sim 0.1$ (see, e.g., Ref. [487] and references therein). Finally, in
Eq. (5.12), the mean number of recombinations $N_{rec}$ can also be taken as a
free parameter, provided that it does not vary much with time. It is expected
to take low values, $N_{rec}<1$ (e.g., Ref. [500]). Note that, at least as far
as ionization is concerned, the relevant quantity is the product of these
parameters, and not the specific value of each one. It also implies
degeneracies between these parameters, which could only be broken if other
observables besides the ionization fraction are considered. Due to the lack of
observations of galaxies at high redshifts, these parameters are still poorly
constrained and remain highly uncertain, and thus $\xi$ is usually taken as a
free parameter. For the default values commented above with Population II
stars, the fiducial value for the ionizing efficiency is $\xi\simeq 40$.
Therefore, from Eqs. (5.6) or (5.12), one can conclude that the parameters
$\xi$ and $M_{\rm vir}^{\rm min}$ (or $T_{\rm vir}^{\rm min}$) rule the
ionization history, and completely determine the filling factor in this simple
picture of global ionization.

To roughly estimate the conditions at which the IGM becomes completely
ionized, one can see in Fig. 2.1 that $f_{coll}(z,M_{\rm vir}^{\rm min})\sim
0.1$ at $z\sim 6$ for $T_{\rm vir}^{\rm min}=10^{4}$ K. Thus, from Eq. (5.12),
to ensure that the IGM becomes completely ionized at that redshift, as
expected from data, it is required that $\xi\gtrsim 10$, which roughly agrees
with the above estimation with the fiducial values. To properly evaluate the
evolution of the filling factor, one must numerically solve Eqs. (5.6) or
(5.12). Examples of ionization histories, computed with 21cmFAST [378, 485],
can be found in Fig. 5.1, for several values of the relevant astrophysical
parameters, $\xi$ and $M^{\rm min}_{\rm vir}$. The plot shows also the current
observational constraints on the ionized fraction, which will be explained in
Sec. 5.1.3. The minimum virial masses of $M^{\rm min}_{\rm vir}=10^{8}$ and
$10^{9}\,M_{\odot}/h$ have been chosen to correspond to minimum virial
temperatures of $T_{\rm min}^{\rm vir}=10^{4}$ (atomic cooling threshold) and
$10^{5}$ K respectively at $1+z\simeq 10$. On the other hand, the considered
values of the ionizing factor $\xi=4$ and $10$ may be related to Population II
($N_{\gamma/b}=4\times 10^{3}$) and Population III ($N_{\gamma/b}=10^{5}$)
stars, respectively, in case $f_{*}\simeq 0.01$ and $f_{esc}\simeq 0.1$.
Notice that both parameters are somewhat degenerate, since increasing
(decreasing) the ionizing efficiency (the minimum virial mass) leads to an
earlier reionization period. Efficiencies of $\xi\gtrsim 4$ with threshold
masses of $\lesssim 10^{9}\,M_{\odot}/h$ are needed in order to have the
universe ionized at $z\sim 6$, as indicated by data (see Sec. 5.1.3).
Molecular cooling, with $T_{\rm min}^{\rm vir}\sim 10^{3}$ K, may lead to a
very early reionization era inconsistent with observations, unless extremely
low ionizing efficiencies were considered.

Many candidates have been considered as the sources of reionization. QSOs
(Quasi-Stellar Objects, also known as quasars) have been largely considered,
since they are among the brightest sources in the universe, being able to emit
in the UV range. However, their number density drops at high redshifts, and
seems to be not enough to contribute significantly to the global ionization of
the IGM, accounting only for 1-5% of the required ionizing background [501].
Standard galaxies become then the preferred option, since we have
observational confirmation of bright galaxies for $z\gtrsim 6$ up to $z\sim
11$ thanks to the Hubble Space Telescope (HST) Ultra Deep Field [502, 503].
Future telescopes such as the JWST [504] could be able to detect farther and
fainter galaxies. Given the slope of the typical luminosity function, dwarf
faint galaxies may be the most relevant ones for producing an ionizing
background due to their high abundance [58]. Regarding the stellar components,
early redshift stars responsible for reionization may belong to the Population
II group, whose stars present low metallicity, about one tenth of the solar
metallicity. However, the first stars of the universe could have had even
lower metallicities, being mostly composed only by Hydrogen and Helium, the
most abundant elements after Primordial Nucleosynthesis. The so-called
Population III stars would constitute a metal-free first stellar generation,
rather different to the other known stellar types [58]. The lack of heavy
elements reduces the cooling capability of the initial collapsing cloud,
leading to much more massive stars than the usually observed, with hotter
central temperatures. As stated above, the number of emitted photons could be
larger than in the Population II case. Moreover, the cooling of the gas
happens by $H_{2}$ formation, leading to lower virial temperatures. This kind
of stars has not been directly observed, yet some indirect evidences are
already present [505].

Figure 5.1: Evolution of the filling factor $Q_{i}$ as a function of redshift,
for different values of $M_{\rm vir}^{\rm min}$ and $\xi$. The included
observational data is listed in Tab. 5.1 and explained in Sec. 5.1.3, arising
from different probes: Gunn-Peterson effect, Ly$\alpha$ emission in galaxies,
damping wings in Ly$\alpha$ spectra from QSOs, dark gaps in quasar spectra and
estimates from the Thomson optical depth of the CMB. The upwards (downwards)
arrows denote lower (upper) limits on $Q_{i}$.

#####  5.1.2 Inhomogeneous Reionization

The discussion above stands for a global reionization scenario, and it is
useful to predict the average ionized volume. However, numerical simulations
agree on the fact that reionization is highly inhomogeneous [506, 507]. In
principle, only hydrodynamic simulations with detailed radiative transfer can
account for the complex physics of the inhomogeneous reionization process.
However, several semi-analytical procedures have been proposed in the
literature in order to account for inhomogeneities with good agreement with
numerical simulations (see, e.g., Ref. [508] for a comparison). In Ref. [509],
the _FZH model_ 444After its authors Furlanetto, Zaldarriaga and Hernquist.
was proposed in order to understand the growth of ionized regions, making use
of the excursion set formalism. The basic idea is to promote Eq. (5.12) to be
scale dependent, writing the filling factor $Q_{i}(z;\delta,M)$ of a region
with a mass $M$ and an overdensity $\delta$ (smoothed within the scale given
by $M$) as

$Q_{i}(z;\delta,M)=\frac{\xi\;f_{coll,cond}(z,M_{\rm vir}^{\rm
min};\delta,M)}{1+N_{rec}},$ (5.13)

where the global collapsed fraction has been replaced by the _conditional_
one, $f_{coll,cond}(z,M_{\rm vir}^{\rm min};\delta,M)$, the fraction of mass
in collapsed objects above $M_{*}$ at a redshift $z$ within a region of mass
scale $M$ and overdensity $\delta$. For the Press-Schechter function, this
reads

$f_{coll,cond}(z,M_{\rm vir}^{\rm
min};\delta,M)=erfc\left(\frac{\delta_{c}(z)-\delta}{\sqrt{2(\sigma(M_{\rm
vir}^{\rm min})^{2}-\sigma(M)^{2})}}\right).$ (5.14)

This is equivalent to assume that the ionized mass is proportional to the mass
of collapsed objects, with the proportionality constant $\xi/(1+N_{rec})$.
Therefore, as done in the standard excursion-set formalism [175], a point of
space ${\bf x}$ with a density contrast $\delta({\bf x})$ smoothed over a
scale $M$ is asigned to an ionized region of mass $M$ if $M$ is the largest
scale at which the condition $Q_{i}(z;\delta({\bf x}),M)\geq 1$ is fulfilled.
This is equivalent to the criterion $\xi f_{coll,cond}(z,M_{\rm vir}^{\rm
min};\delta({\bf x}),M)/(1+N_{rec})\geq 1$, where this last inequality imposes
a condition over $\delta$ in the ionized bubbles, implying that only regions
dense enough are capable of self-ionize, i.e the condition:

$\delta\geq\delta_{c}(z)-erfc^{-1}\left((1+N_{rec})/\xi)\sqrt{2(\sigma(M_{*})^{2}-\sigma(M)^{2})}\right)$
(5.15)

must hold within the region in order to be ionized. Note that the expression
above is valid only for the PS prescription. For these reasons, this is
considerered as an inside-out scenario of reionization, where, on average,
high-density regions are firstly ionized. The threshold overdensity at a scale
$M$ can be considered a moving barrier to construct a HII mass function in the
excursion set formalism [509]. This prescription is well suited for numerical
implementation and it is computationallt very efficient, since it only needs
to filter the density field, marking as ionized the cells above the threshold
at a given filter scale. For this reason, it is the method used in many codes
such as 21cmFAST [378]. It has been shown to be in good agreement with more
computationally expensive radiative transfer codes at scales $\gtrsim 1$ Mpc
[510].

There are other analytical models to treat inhomogenous reionization, such as
the known as _MHR model_ 555After Miralda-Escudé, Haehnelt and Rees. proposed
in [511]. In this prescription, it is assumed that only regions with
overdensities above a certain treshold remain neutral, as a consequence of an
enhanced recombination rate. Contrarily to the FZH model, this is an outside-
in scenario, where the most dense regions ionize later. Although it may seem
in contradiction with FZH, the application of MHR is appropriate in the later
stages of reionization, when only small overdense spots of the IGM remain
neutral. Moreover, it is evaluated in a single, small scale, instead of the
multiscale approach of FZH, which is able to describe the morphology of the
reionization process on large scales. Both methods can be combined and provide
an accurate description of the reionization from different perspectives [508].

Finally, although our focus in this thesis is on semi-analytic methods to
solve the ionization history, a proper determination of the IGM evolution
requires more accurate though computationally expensive methods. Some of them
involve N-body simulations to compute the density fields and halo population,
combined with a detailed radiative transfer treatment, via, e.g., adaptive ray
tracing algorithms [512, 513]. However, hydrodynamic radiative transfer
simulations are the most reliable at the smallest scales, able to reach
spatial resolutions of $\sim 100$ pc - $\sim 10$ kpc, although are limited to
boxes of tens of Mpc. They account for different baryonic effects absent in
other prescriptions, which may include stellar and supernovae feedback, star
formation, fully coupled three-dimensional radiative transfer and detailed
cooling and heating for the metal enriched gas [514, 515, 516]. Radiative
transfer can be treated by three main approaches, namely employing moments
methods (e.g., [517]), Monte Carlo sampling (e.g., [518]) or ray tracing
(e.g., [519]). See Refs. [520, 521] for detailed comparisons between codes and
radiative transfer methods, and also, e.g., Refs. [522, 523] for overviews and
more details on this topic.

#####  5.1.3 Constraints on the Epoch of Reionization

Data |  Redshift |  $Q_{i}$ |  Reference
---|---|---|---
Gunn-Peterson effect | 5.03 |  0.9999451${}_{-0.0000165}^{+0.0000142}$ | Fan et al., 2006 [524]
5.25 | 0.9999330${}_{-0.0000244}^{+0.0000207}$
5.45 | 0.9999333${}_{-0.0000301}^{+0.0000247}$
5.65 | 0.9999140${}_{-0.0000460}^{+0.0000365}$
5.85 | 0.9998800${}_{-0.0000490}^{+0.0000408}$
6.10 | 0.99957$\pm$0.00030
Dark gaps in quasar spectra | 5.6 |  $>$0.91 | McGreer et al., 2015 [525]
5.9 |  $>$0.89
6.1 |  $>$0.42
Ly$\alpha$ Emission in Galaxies | 7 | 0.66 ${}^{+0.12}_{-0.09}$ | Schenker et al., 2014 [526]
8 | $<0.35$
6.6 | $0.7\pm 0.2$ | Konno et al., 2018 [527]
7 | $0.5^{+0.3}_{-0.1}$ | Inoue et al., 2018 [528]
$\sim 7$ | $0.41^{+0.15}_{-0.11}$ | Mason et al., 2018 [529]
7.3 | $0.75\pm 0.25$ | Itoh et al., 2018 [530]
$7.6\pm 0.6$ | $0.12^{+0.10}_{-0.05}$ | Hoag et al., 2019 [531]
$7.9\pm 0.6$ | $<0.26$ | Mason et al., 2019 [532]
Clustering of Ly$\alpha$ Emitters | 6.6 | $\gtrsim 0.5$ | McQuinn et al., 2007 [533]
6.6 | $\gtrsim 0.5$ | Ouchi et al., 2010 [534]
7 | $\gtrsim 0.5$ | Sobacchi and Mesinger, 2015 [535]
QSO Ly$\alpha$ damping wings | 6.24-6.42 | <0.9 (2$\sigma$) | Schroeder et al., 2013 [536]
7.00 | $0.30^{+0.23}_{-0.20}$ | Wang et al., 2020 [537]
7.09 | $0.60^{+0.19}_{-0.21}$ | Greig et al., 2017 [538]
7.52 | $0.61^{+0.13}_{-0.22}$ | Yang et al., 2020 [539]
7.54 | $0.79^{+0.19}_{-0.17}$ | Greig et al., 2018 [540]
7.54 |  $0.44^{+0.18}_{-0.21}$ (95% CL) | Bañados et al., 2018 [541]
7.09 | $0.52\pm 0.26$ | Davies et al., 2018 [542]
7.54 | $0.40^{+0.23}_{-0.20}$
Thomson Optical Depth | $8.8\pm 0.9$ | 0.5 | Adam et al., 2016 [543]
$\tau_{T}=0.0519^{+0.0030}_{-0.0079}$ | Aghanim et al., 2020 [24]

Table 5.1: List of representative constraints on $Q_{i}$ at $1\sigma$ (unless
stated otherwise). Data points are plotted in Fig. 5.1, see, e.g., Refs. [544,
545] for comprehensive (though slightly outdated) compilations of data.

Currently, several estimations of the global ionized fraction of the IGM at
different redshifts have been performed. Although the end of reionization is
well determined, the precise evolution at lower redshifts becomes much more
uncertain, being some of the estimations strongly dependent on the underlying
assumed models. In this section we review some of the most important methods
to bound the onset and timing of the EoR. Table 5.1 shows a summary of the
most important constraints, which are also plotted in Fig. 5.1.

######  5.1.3.1 Gunn-Peterson effect

The first, and one of the most important probes of the EoR, is the Gunn-
Peterson effect, named after the physicists who predicted it in 1965 [546].
Due to the fact that most of the Hydrogen atoms in the IGM are in their
fundamental state, the relevant transitions are those which relate its ground
state with the excited ones, as the Ly$\alpha$ transition. Many astrophysical
objects like galaxies and quasars can emit high energy radiation which, due to
the expansion of the universe, redshift to lower frequencies along its
journey. In their passing across the IGM, these photons could eventually
redshift to the wavelength corresponding to the Ly$\alpha$ transition. Thus,
any neutral Hydrogen atom around could absorb it, preventing us from seeing
it. Since the IGM is densely filled by Hydrogen atoms, these absorptions would
be likely even if a considerable fraction of Hydrogen atoms were ionized. The
efficiency of this absorption is given by the optical depth of the Ly$\alpha$
transition, the so-called Gunn-Peterson optical depth $\tau_{\rm GP}$, already
computed in Eq. (4.36), which can be approximated as $\tau_{\rm GP}\simeq
5.0\times 10^{5}\,x_{\rm HI}\left((1+z)/7\right)^{3/2}$. From that formula,
one can clearly see that the medium is optically thick ($\tau_{\rm GP}>1$)
unless it is highly ionized, e.g. $x_{\rm HI}<10^{-5}$ at $z\sim 6$. Hence,
observing QSO spectra embedded in a mostly neutral IGM, we would see an
absorption trough (known as Gunn-Peterson trough) just above the Ly$\alpha$
transition (since photons are redshifted from higher frequencies). However, we
actually can observe light from relatively distant quasars without this
expected suppression, up to $z\sim 6$, when a significant damping in the
spectrum appears. This means that, for $z\lesssim 6$, essentially there is no
neutral gas in the IGM, being most of the atoms in a highly ionized state. For
higher redshifts, however, the medium appears to be still predominantly
neutral. This transition between these two regimes is the so-called EoR.

The first evidence of this GP trough at $z\sim 6$ was reported in Ref. [31]
from spectroscopic observations of quasars in the Keck Observatory. The Gunn-
Peterson effect has been extensively used in order to constrain the end of
this epoch, observing the light emitted by distant quasars. The detection of
the Gunn-Peterson trough in the spectra allows us to determine the ionized
fraction with great precision. This method was used in Ref. [524] to estimate
the neutral fraction $1-Q_{i}$ between $\sim 10^{-5}$ and $\sim 10^{-4}$ at
six different redshifts $z=5.03,$ $5.25,$ $5.45,$ $4.65,$ $5.85$ and $6.10$,
by means of the observation of spectra of 19 bright quasars. These
measurements are not, however, model-independent, since assumptions such as
the temperature-density relation of the IGM or the PDF of overdensities were
assumed.

On the other hand, there is a related method to constrain the ionized fraction
less subject to the modeling of the IGM. In Ref. [525], lower limits on
$Q_{i}$ were placed from the distribution of dark gaps in quasar spectra at
redshifts $5.6$, $5.9$ and $6.1$ [525] (indicated in Tab. 5.1 as dark gaps in
quasar spectra). These results were achieved by determining the fraction of
binned spectra in Ly$\alpha$ and Ly$\beta$ forests which were dark, i.e.,
presented zero flux. Although less restrictive, these are robust limits
obtained in a model-independent way, unlike the estimations from Ref. [524].
All these data indicate that reionization has to be completed by $z\sim 6$. It
is worth mentioning that observations of the Ly$\alpha$ forest work under the
same idea than the GP effect. While the GP trough is a deep suppression of the
spectrum due to a not completely neutral IGM, the Ly$\alpha$ forest arises
from neutral clouds of gas embedded in the ionized IGM, which imprint narrow
absorption dips in the QSO spectrum determined by their size. This allows to
tightly constrain the residual neutral fraction in the post-EoR universe via
estimations of the optical depth [547, 548, 549].

######  5.1.3.2 Thomson Optical Depth

Once reionization has started, the increase in the number of free electrons
makes Thomson scattering to be relevant again, as it happened before the
Recombination era. This scattering isotropizes the radiation field, being the
angular power spectrum of the CMB suppressed by a damping factor
$e^{-2\tau_{T}}$, where $\tau_{T}$ is the Thomson optical depth, defined as

$\tau_{T}=\int\frac{cdz}{(1+z)H(z)}\,\sigma_{T}\,n_{e}(z)\leavevmode\nobreak\
,$ (5.16)

with $\sigma_{T}=0.6652\times 10^{-24}$ cm2 the Thomson cross section. This
damping enforces the radiation field to become isotropic for large enough
$\tau_{T}$, erasing anisotropies. As a consequence, the reionization epoch
leaves imprints on the CMB spectrum, being $\tau_{T}$ the quantity which
characterizes its impact. This is one of the six parameters of the baseline
$\Lambda$CDM cosmological model, determined with great precision by the Planck
mission [550]. Although in the temperature $TT$ spectrum there is a degeneracy
between the Thomson optical depth and the amplitude of the primordial
fluctuations, it is broken in the polarization of the CMB spectrum, since
Thomson scattering induces polarization at low multipoles [551].
Parameterizing the evolution of $Q_{i}$ with $z$ as a $\tanh$ function, the
latest Planck analysis offers a value for this quantity of [24]

$\tau_{T}=0.0519^{+0.0030}_{-0.0079}\leavevmode\nobreak\ ,$ (5.17)

making use of the low-$l$ $EE$ modes likelihood. Previous Planck analyses
[543] evaluated also the beginning, the final and the duration of reionization
based on a slightly larger, earlier estimation of the optical depth,
$\tau_{T,2016}=0.058\pm 0.012$. Within the $tanh$ parameterization, the
reionization process occurs suddenly enough around a characteristic redshift
of $z_{re}=8.8\pm 0.9$ (where the ionization fraction reaches 50%), with a
total duration of $\Delta z<4.6$ (at $95\%$ CL), with the additional
constraint that the IGM must be fully ionized at $z<6$, due to the Gunn-
Peterson effect. These results change only slightly if an asymmetric power-law
parameterization is assumed, rather than a $tanh$ function [543]. Besides the
functional parameterization, other approaches have been used in the literature
to exploit the $\tau_{T}$ constraint, such as the PCA (Principal Component
Analysis) decomposition [552, 553, 554, 3], the PCHIP (Piecewise Cubic Hermite
Interpolating Polynomial) [3] or CMB plus reionization modeling [555]. In Part
II of this thesis, CMB data is exploited in order to reconstruct the
ionization history by means of the aforementioned PCA and PCHIP analyses,
studying the plausibility of an early reionization period [3].

Measurements of the Thomson optical depth constitute a robust indication of
reionization and allow us to parameterize some general properties of the EoR.
However, in order to infer timing, duration or patchiness aspects of the EoR,
specific parameterizations are needed, since this does not provide us with
information about the evolution of the free electron fraction at each
redshift, but instead about its integral along time.

######  5.1.3.3 Other constraints

While the Gunn-Peterson trough and the CMB measurements provide the most
robust bounds on the EoR, other methods may also lead to constraints on the
ionization history. The most relevant ones are summarized in the following.

  * •

Ly$\alpha$ emission in galaxies

One of the most profitable approaches is based on the observation of
Ly$\alpha$ emitting galaxies, usually presenting large star formation rates,
which can constitute a successful approach for exploring how reionization has
proceeded at higher redshifts ($z\gtrsim 7$) [556]. The Subaru Hyper Suprime-
Cam has provided measurements of thousands of luminosities of Ly$\alpha$
emitters through several wavelength bands, which allows placing constraints on
the reionization history. Examples of this are the survey SILVERRUSH [557],
from which several estimates of the neutral fraction have been derived at
redshifts $z=6.6$ [527] and $z=7.3$ [528], and the CHORUS survey [558], which
has provided constraints at $z\simeq 7$ [530]. Other constraints rely on
Ly$\alpha$ emission from faint Lyman break galaxy candidates, those detected
by a drop in the spectrum at frequencies above the Lyman limit, caused because
ionizing photons are absorbed in their local environment. The detection and
non-detection of Ly$\alpha$ emission in such galaxies allows placing bounds
around $z\sim 7$ [526, 529], $z\sim 7.6$ [531] and $z\sim 8$ [526, 532]. All
the above results indicate that reionization is not yet complete at those
epochs, $z\sim 7$. One of the caveats of employing Ly$\alpha$ emitters to
infer the ionization fraction is that radiative transfer modeling is needed,
these results being, therefore, quite sensitive to the assumptions considered
and to the details of the numerical implementations.

  * •

Clustering of Ly$\alpha$ emitters

Constraints closely related to the above ones can be derived from the
clustering properties of Ly$\alpha$ emitters, rather than their luminosity
functions. An enhancement in the clustering of Ly$\alpha$ emitters could be
regarded as a consequence of reionization, since radiation from clumped
galaxies would suffer less absorption than from isolated ones. This effect can
be hardly mimicked via other processes, thus indicating the progress of the
EoR. It has allowed to place lower limits on the ionized fraction at $z=6.6$
[533, 534] and $z\simeq 7$ [535]. The main difficulty of this method is that
these estimates require a great number of Ly$\alpha$ emitters, in order to
properly infer the clustering properties.

  * •

Ly$\alpha$ damping wing

Ly$\alpha$ photons absorbed in the tails (_damping wings_) of the cross
section instead of the resonance could also provide a powerful probe of the
EoR. Ly$\alpha$ cross section can be modeled by a Voigt profile, Eq. (4.43),
which near the resonance behaves similar to a Gaussian distribution, but far
enough from the center of the line, decays much slower as a Lorentzian
function. One can split the optical depth in two different contributions, one
from the center of the line, and the other from the tails of the distribution,
which is a smoother function of the frequency. The strength of the damping
wing absorption strongly depends on the neutral fraction of the environment,
fact which allows placing constraints on $Q_{i}$ [559, 560]. Taking advantage
of its spectral smoothness, it was possible to set upper limits on $Q_{i}$
between $z=6.24$ and $6.42$ [536]. Stronger bounds can be obtained by studying
the absorption on the red side of the Ly$\alpha$ line, where resonant
absorption is negligible, and thus must be produced at the tails. The
discovery of four of the farthest observed quasars has allowed placing
constraints using that method at redshifts $z=7.00$ [537], $z=7.09$ [538,
542], $z=7.52$ [539] and $z=7.54$ [540, 541, 542]. The derivation of all these
constraints based on the Ly$\alpha$ damping wing also rely on radiative
transfer and hydrodynamical simulations. For the same redshifts, different
methodologies and model assumptions may provide differing constraints (see
Tab. 5.1).

On the other hand, there is recent evidence of overlapping ionized bubbles
around three galaxies at $z\sim 7.7$ from spectroscopic observations [561].
Finally, as commented in the previous chapter, although the 21 cm signal have
not been measured at these redshifts, some constraints have been obtained on
the duration and completion redshift of the EoR. The EDGES high-band
observations are capable of ruling out too abrupt and rapid reionization
processes [403]. Other compilations of ionization data can be found in Refs.
[544, 545].

####  5.2 Neutral IGM

In this section, the relevant physics of the neutral IGM is outlined. Before
ionized regions around each galaxy grow and merge, the surrounding medium is
still in a neutral state, with a remnant ionized fraction from Recombination.
The astrophysical processes taking part in galaxies and other sources can
produce X-ray radiation capable of reaching areas far beyond the galaxies,
altering its thermal and ionization state. Moreover, Ly$\alpha$ radiation is
expected to fill the neutral medium, of great importance for the 21 cm signal.
In the following, we comment on the heating, ionization and Ly$\alpha$ flux
present in the mostly neutral IGM.

#####  5.2.1 Heating of the IGM

From energy conservation, one can derive a differential equation which
determines the kinetic temperature of the gas, $T_{K}$. The thermal evolution
equation reads (e.g., Ref. [562])

$\frac{dT_{K}}{dt}+2HT_{K}-\frac{2}{3}\frac{T_{K}}{(1+\delta)}\frac{d\delta}{dt}+\frac{T_{K}}{1+x_{e}}\frac{dx_{e}}{dt}=\frac{2\mathcal{Q}}{3\,n_{b}(1+x_{e})}\leavevmode\nobreak\
,$ (5.18)

where $\mathcal{Q}$ is the total heating rate per unit volume. The second and
third terms of the left-hand side account for adiabatic expansion, the fourth
one induces changes in the temperature given by the ionization variation,
whereas the right-hand side term includes any possible interacting process
which can cool down or heat up the gas. Note that, in the absence of
heating/cooling $\mathcal{Q}$ terms, the above equation can be anallitically
solved, obtaining $T_{K}\propto a^{-2}(1+\delta)^{2/3}/(1+x_{e})$. In the
homogeneous, post-recombination limit, we have the well-known adiabatic
cooling evolution $T_{K}\propto a^{-2}$ as a consequence of the expansion of
the universe.

There are several possibles sources of heating and cooling in the IGM which
contribute to $\mathcal{Q}$. In the following, we shall comment on the most
relevant ones to take into account. While in this section, only the standard
processes are considered, further heating terms may arise from exotic DM
scenarios, such as accretion in PBHs. The heating effects in that case are
widely discussed in Part II of this thesis [6].

######  5.2.1.1 Compton cooling

Compton scattering $e+\gamma\leftrightarrow e+\gamma$ is the responsible of
the coupling between matter and radiation in the early universe. Even after
the decoupling of the CMB, the gas is still affected by scatterings with
photons, due to the large number of photons compared to baryons, enforcing the
gas temperature $T_{K}$ to be equal to the CMB temperature $T_{\gamma}$. The
corresponding cooling rate can be derived from the Boltzmann equation for the
Compton scattering, after performing a Fokker-Planck expansion (similar to Eq.
(4.35)) and obtaining the so-called Kompaneets equation [563]. Integrating it,
the Compton cooling rate is obtained [164],

$\mathcal{Q}|_{Compton}=\frac{4\sigma_{T}n_{e}\rho_{\gamma}}{m_{e}c}\left(T_{\gamma}-T_{K}\right),$
(5.19)

where $\rho_{\gamma}=(\pi^{2}/15)T_{\gamma}^{4}$ is the CMB energy density.
This result was first derived in Ref. [564]. At high redshifts, when Compton
cooling is very efficient and there are no other relevant sources of heating,
the gas temperature is coupled to the CMB temperature $T_{K}\rightarrow
T_{\gamma}\propto(1+z)$. However, due to the expansion of the universe, there
is a moment when the adiabatic cooling term becomes dominant, and the Compton
interactions no longer can couple both temperatures. This thermal decoupling
happens at a redshift $z_{D}$ when both rates (adiabatic and Compton) are
similar

$\displaystyle\left.2HT_{K}\right|_{z_{D}}\sim\left.\frac{2\mathcal{Q}|_{Compton}}{3k_{B}n_{b}(z_{D})(1+x_{e}^{\infty})}\right|_{z_{D}}\;$
(5.20)
$\displaystyle\Rightarrow\;1+z_{D}\simeq\left(\frac{45m_{e}H_{0}\sqrt{\Omega_{m}}}{4\pi^{2}\sigma_{T}T_{\gamma,0}x_{e}^{\infty}}\right)\simeq
150,$ (5.21)

where we have employed $x_{e}^{\infty}\simeq 3\times 10^{-4}$ as the leftover
free electron fraction from Recombination, which is roughly constant, as
estimated in Sec. 5.2.2 or computed via publicly available Recombination codes
such as RECFAST [565]. Therefore, we expect an adiabatic evolution of the
kinetic temperature $T_{K}\propto(1+z)^{2}$ for redshifts below $z_{D}\sim
150$, as long as no other heating sources are present.

######  5.2.1.2 X-ray heating

After the onset of star formation, some astrophysical objects such as high
mass X-ray binaries in starbust galaxies, supernova remnants or black holes in
mini-quasars may emit a strong background radiation at UV and X-ray
frequencies [566, 353]. Given that the ionization cross section can be
approximated as $\sigma(\nu)\propto\nu^{-3}$ (see, e.g., Ref. [30]), the mean
free path of photons goes as $\lambda_{mfp}\sim\nu^{3}$. While UV photons have
a lower mean free path and are absorbed in the neighborhood of astrophysical
sources, ionizing their environment, the more energetic X-rays can travel much
farther, being able to reach the mostly neutral IGM, distant from galaxies.
Photons with $\sim$ keV energies can therefore be absorbed, transferring their
energy either to ionization, heating or atomic transitions.

Solving the radiative transfer equation for X-rays, Eq. (4.1), it is possible
to obtain obtain the specific X-ray intensity (by number) $J_{\rm X}$ (in
units of s-1 keV-1 cm-2 sr-1),666Note that other convention in the literature
denotes with $J_{\rm X}$ the flux by energy, rather than by number, having
thus units of erg s-1 keV-1 cm-2, and implying an extra factor of $h\nu$. For
that quantity, we employ $I$. which reads

$J_{\rm
X}(z,\nu)=\frac{(1+z)^{2}}{4\pi}\int_{z}^{\infty}dz^{\prime}\frac{dt}{dz^{\prime}}\,(1+z^{\prime})\varepsilon_{\rm
X}(\nu^{\prime})\,e^{-\tau_{\rm X}(\nu,z,z^{\prime})}\leavevmode\nobreak\ ,$
(5.22)

where $\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$ is the comoving frequency at
emission, and $\tau_{\rm X}(\alpha,m_{e})$ is the IGM optical depth from
$z^{\prime}$ to $z$. This factor accounts for the attenuation of the X-ray
flux through the IGM. In analogy to Eq. (5.7), the comoving specific X-ray
emissivity $\epsilon_{\rm X}(z,\nu)$ is modeled as

$\varepsilon_{\rm X}(\nu)=\frac{L_{\rm X}(\nu)}{h\nu\,{\rm
SFR}}\dot{\rho}_{*}=\frac{L_{\rm X}(\nu)}{h\nu\,{\rm
SFR}}f_{*}\bar{\rho}_{b,0}\frac{df_{coll}(z,M_{\rm vir}^{\rm min})}{dt},$
(5.23)

where $L_{\rm X}/{\rm SFR}$ is the specific differential X-ray luminosity
normalized to the star formation rate escaping the host galaxies (in units of
erg s-1 Hz-1 $M_{\odot}^{-1}$ yr), whose spectral shape can be assumed to be a
power law, $L_{\rm X}\propto(\nu/\nu_{0})^{-\gamma_{\rm X}}$, with $\nu_{0}$ a
threshold frequency. The slope takes values around $\gamma_{\rm X}\simeq 1$,
ranging between $0.8$ and $1.5$, depending on the specific X-ray source [566,
494, 567]. High-redshift X-ray emitters are not directly observed yet, and
thus we shall consider the luminosity normalization as a free parameter, given
by the integrated luminosity (per SFR) over the energy band
$(2-10\leavevmode\nobreak\ {\rm keV})$, $L_{{\rm X}[2-10\leavevmode\nobreak\
{\rm keV}]}$/SFR.777There are other parameterizations to account for the
normalization, such as those followed in the early versions of 21cmFAST code
[378]. Instead of the luminosity, the heating efficiency $\xi_{X}$ was chosen,
defined as the number of X-ray photons per solar mass in stars, which can be
easily related to $L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm keV}]}/{\rm SFR}$.
We shall take as a fiducial value the extrapolation from nearby starbusts
galaxies, $L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm keV}]}/{\rm SFR}=3.4\times
10^{40}$ erg s-1 $M_{\odot}^{-1}$ yr [566, 494]. More recent estimates agree
with that default value, showing luminosities from X-ray compact binaries
between $10^{39}$ and $10^{40}$ erg s-1 $M_{\odot}^{-1}$ yr between redshifts
0 and 4 [568]. The threshold frequency $\nu_{0}$ is tipically assumed to lie
between $h\nu_{0}=300$ and $500$ eV, but its specific value is less relevant
given is degeneracy with $L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm keV}]}/{\rm
SFR}$ and the slope $\gamma_{X}$. Frequencies below $\nu_{0}$ are locally
absorbed in the nearby region of the emitting source, and do not reach the
IGM. Comparing to Eq. (5.23), one can see that the X-ray escape fraction is
effectively set to 1, since most of the X-rays are able to get away galaxies
due to the high mean free path.

The heating rate by X-ray sources is thus given by

$\mathcal{Q}|_{X}(z)=4\pi\,n_{b}(z)\int_{\nu_{0}}^{\infty}d\nu\,J_{\rm
X}(z,\nu)\,\sum_{j}(h\nu-E_{j}^{\rm th})f_{\rm
heat}\,f_{j}\,x_{j}\,\sigma_{j}(\nu)\leavevmode\nobreak\ ,$ (5.24)

where $\sigma_{j}$ are the photoionization cross-sections of species $j$, and
$f_{\rm heat}$ accounts for the fraction of energy from X-rays that is
converted into heat in the IGM, using the fit from Ref. [569]. Since
evaluating the integrals of Eqs. (5.22) and (5.24) together with the optical
depth can be quite time-consuming, some numerical implementations (such as the
one followed in the 21cmFAST code [570, 378, 485]) approximate the absorption
factor $e^{-\tau_{\rm X}}$ to speed-up the computation. Given its exponential
form, we can assume that all photons with optical depth large enough,
$\tau_{\rm X}>1$, are easily absorbed in the medium, whereas photons with
$\tau_{\rm X}\leq 1$ are able to travel without being absorbed. It is
equivalent to approximate the absorption factor as a Heaviside step function,
$e^{-\tau_{\rm X}}\simeq\theta\left[1-\tau_{\rm X}(\nu,z,z^{\prime})\right]$,
avoiding to integrate over the optical depth,

$\begin{split}\mathcal{Q}|_{X}(z)&=(1+z)^{2}n_{b}(z)\,\int_{z}^{\infty}dz^{\prime}\frac{dt}{dz^{\prime}}(1+z^{\prime})\,\int_{Max[\nu_{0},\nu_{\tau=1}]}^{\infty}d\nu\,\sum_{j}(h\nu-
E_{j}^{\rm th})f_{\rm heat}\\\
&\times\,f_{j}\,x_{j}\,\sigma_{j}(\nu)\varepsilon_{\rm
X}(\nu^{\prime})\leavevmode\nobreak\ ,\end{split}$ (5.25)

where $\nu_{\tau=1}$ is the frequency where the optical depth equals unity.

The function $f_{\rm heat}$ accounts for the fraction of X-ray energy
deposited as heating. Other amounts of this energy are absorbed in form of
ionization and excitation, having therefore the corresponding fractions
$f_{\rm ion}$ and $f_{\rm exc}$. Although these functions depend on the X-ray
energy and the ionization fraction $x_{e}$, they saturate at energies $\gtrsim
100$ eV, having a function of only $x_{e}$ [571, 569]. A rough approximation
gives [572]

$f_{\rm heat}\sim\frac{1+2x_{e}}{3},\;\;\;f_{\rm ion}\sim f_{\rm
exc}\sim\frac{1-x_{e}}{3}.$ (5.26)

Then, in mostly neutral media, the energy splits equally in the three
channels, while it is mainly deposited as heating when the ionization is high.

To gain insight about the impact of X-rays in the temperature evolution, we
can estimate the total heating by taking the on-the-spot approximation in the
flux (i.e., employing Eq. (5.4) in Eq. (5.24)), getting a simpler heating rate
(e.g., [494]),

$\displaystyle\mathcal{Q}|_{X}(z)=\frac{L_{{\rm X}[2-10\leavevmode\nobreak\
{\rm keV}]}}{{\rm SFR}}\,f_{*}f_{\rm
heat}\,\rho_{b}(z)\,\frac{df_{coll}(z,M_{\rm vir}^{\rm min})}{dt}.$ (5.27)

Thus, the increase in temperature per baryon in a Hubble time is

$\begin{split}\frac{2\,\mathcal{Q}|_{X}}{3\,H\,n_{b}}&\simeq 10\,{\rm
K}\;\left(\frac{-(1+z)df_{coll}(z,M_{*})/dz}{10^{-3}}\right)\\\
&\times\left(\frac{L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm keV}]}/{\rm
SFR}}{3.4\times 10^{40}erg\cdot
s^{-1}M_{\odot}^{-1}yr}\right)\left(\frac{f_{*}}{0.01}\right)\left(\frac{f_{\rm
heat}}{0.2}\right),\end{split}$ (5.28)

assuming $df_{coll}(z,M_{\rm vir}^{\rm min})/dt/H(z)=-(1+z)df_{coll}(z,M_{\rm
vir}^{\rm min})/dz\simeq 10^{-3}$ as a fiducial value at $z\simeq 20$. Note
that, given the growth of the collapse factor, the above formula increases by
an order of magnitude at $z\simeq 15$, as can be noticed in Fig. 2.1. As a
further estimation, ignoring adiabatic terms and changes in the ionization
fraction in the temperature equation, Eq. (5.18), we can integrate the
temperature taking Eq. (5.27) as the source in order to estimate the order of
magnitude increase in temperature generated by X-rays, $\Delta T|_{X}(z)$,

$\displaystyle\Delta T|_{X}(z)$ $\displaystyle=\frac{2}{3}\frac{L_{{\rm
X}[2-10\leavevmode\nobreak\ {\rm keV}]}}{{\rm SFR}}\,f_{*}f_{\rm heat}\,\mu
m_{p}\,\Delta f_{coll}(z,M_{*})$ (5.29) $\displaystyle\begin{split}&\simeq
10^{4}{\rm K}\;\Delta f_{coll}(z,M_{\rm vir}^{\rm min})\left(\frac{L_{{\rm
X}[2-10\leavevmode\nobreak\ {\rm keV}]}/{\rm SFR}}{3.4\times 10^{40}erg\cdot
s^{-1}M_{\odot}^{-1}yr}\right)\\\
&\times\left(\frac{f_{*}}{0.01}\right)\left(\frac{f_{\rm
heat}}{0.2}\right),\end{split}$ (5.30)

where we have assumed the fiducial values $L_{{\rm X}[2-10\leavevmode\nobreak\
{\rm keV}]}/{\rm SFR}\simeq 3.4\times 10^{40}ergs^{-1}M_{\odot}^{-1}yr$,
$f_{\rm heat}\simeq 0.2$ and $f_{*}\simeq 0.01$ [494]. Therefore, an increase
in the collapse fraction of $\Delta f_{coll}(z,M_{\rm vir}^{\rm min})\sim 1$
would lead to heat up the medium by $\sim 10^{4}$ K. Note however, that the
on-the-spot scenario would only be accurate for low redshifts, when the mean
free path is much shorter than the Hubble length, and for quantitative
computations, Eq. (5.24) (or Eq. (5.25)) must be used. Figure 5.2 shows the
temperature prediction for different astrophysical models, computed with
21cmFAST. As can be seen, the minimum virial temperature sets the timing of
the onset of heating, being earlier for lower threshold masses, as lower mass
halos are allowed to host star formation. On the other hand, the X-ray
luminosity determines the strenght of heating. The spin temperature, treated
in Sec. 4.2, is also shown, being coupled to the kinetic one by means of
scattering of Hydrogen atoms by Ly$\alpha$ radiation, whose flux is commented
in Sec. 5.2.3.

Figure 5.2: Evolution of the relevant temperatures as a function of redshift,
for different values of the minimum virial mass and X-ray luminosity per star
formation rate (in units of erg s${}^{-1}M_{\odot}^{-1}$yr). Thin lines stand
for the kinetic temperature of the gas $T_{K}$, thick lines for the spin
temperature $T_{S}$, while the blue line corresponds to the CMB temperature
$T_{\gamma}$.

######  5.2.1.3 21 cm heating

Recently, a new heating channel was proposed based on CMB photons interacting
via the 21 cm line [382]. The absorption and emission of CMB photons by the
hyperfine levels in Hydrogen atoms can also transfer some thermal energy in
the process. Its corresponding heating rate can be written as [382]

$\mathcal{Q}|_{21}=\frac{3}{4}\;n_{H}\;x_{\rm
HI}\;x_{CMB}\;A_{10}\;T_{*}\;\left(\frac{T_{\gamma}}{T_{S}}-1\right),$ (5.31)

where $x_{CMB}=(1-e^{-\tau_{21}})/\tau_{21}\simeq 1$, where the 21 cm optical
depth $\tau_{21}$ is small. Note that this heating implicitly depends on the
Ly$\alpha$ flux, and therefore on the astrophysical parameters, through
$T_{S}$. Given the factor $T_{\gamma}/T_{S}-1$, this heating rate turns on
when $T_{S}$ deviates from $T_{\gamma}$, during the Ly$\alpha$ coupling
period, vanishing however when the hyperfine transitions are driven by CMB
photons. Since this epoch (partly) coincides with the onset of X-ray heating,
it has a subdominant impact on the thermal history of the IGM, unless X-ray
emitters are quite faint [9]. Note that in the limit when the medium is very
hot, and the CMB is cooler than the IGM, $T_{\gamma}/T_{K}<1$, and cooling
instead of heating takes place. The maximum possible cooling from this source
appears in the limit when $T_{\gamma}/T_{K}$ can be neglected, and thus
$T_{\gamma}/T_{S}-1\simeq-1$.

One can easily estimate quantitatively this effect. Approximating the color
temperature to the kinetic one (which is a fair assumption, as commented in
Sec. 4.3.3), we can write

$\frac{T_{\gamma}}{T_{S}}-1\simeq\frac{x_{\alpha}}{1+x_{\alpha}}\left(\frac{T_{\gamma}}{T_{K}}-1\right).$
(5.32)

In Sec. 5.2.3.1, the Ly$\alpha$ coupling coefficient is estimated as
$x_{\alpha}\sim 0.5$ at $z\sim 20$ in the fiducial model. The largest possible
values for $T_{\gamma}/T_{K}$ are reached when heating sources are not
efficient enough and the kinetic temperature still cools adiabatically, having
$T_{\gamma}/T_{K}=(1+z_{D})/(1+z)$. Then, at $z\sim 20$, $T_{\gamma}/T_{K}\sim
7$ and $T_{\gamma}/T_{S}-1\sim 2$. In real scenarios, heating sources could be
already important at these times to produce temperature deviations from the
adiabatic evolution, lowering this value. Therefore, it must be taken as an
upper limit. Hence, the heat per baryon and per Hubble time reads

$\displaystyle\frac{2\,\mathcal{Q}|_{21}}{3\,H\,n_{b}}$
$\displaystyle\simeq\frac{1}{2}\;\mathfrak{f}_{H}\;\frac{A_{10}\;T_{*}}{H}\;\left(\frac{T_{\gamma}}{T_{S}}-1\right)$
(5.33) $\displaystyle\simeq 3\;{\rm
K}\;\left(\frac{21}{1+z}\right)^{3/2}\;\left(\frac{T_{\gamma}/T_{S}-1}{2}\right),$
(5.34)

which is roughly an order of magnitude lower than the corresponding X-ray heat
from Eq. (5.28). For lower redshifts, this difference even increases: while
the X-ray heating increases by an order of magnitude at $z\simeq 15$, Eq.
(5.34) would reach as much as $\sim 20$ K at the same epoch, if the assumption
of adiabatic evolution is still valid (which should not be true, since heating
sources could already be important) and the Ly$\alpha$ coupling is completely
efficient. These considerations make the 21 cm heating term negligible in most
of astrophysical scenarios, and consequently, is usually not considered, as
happens in 21cmFAST. It must be remarked, however, that if X-ray luminosities
are much lower than expected, this would be the dominant heating channel at
high redshifts. Furthermore, it could also be relevant for even higher
redshifts, deep in the Dark Ages, when collisional coupling becomes efficient,
allowing the factor $T_{\gamma}/T_{S}-1$ to deviate from 0.

######  5.2.1.4 Ly$\alpha$ heating

Ly$\alpha$ scattering with atoms can transfer a fraction of energy from
radiation to the gas, with a contribution from continuum and injected photons.
The heating rate, already computed in Sec. 4.3.4, reads

$\mathcal{Q}|_{\rm Ly\alpha,k}=\frac{4\pi
Hh\nu_{\alpha}\Delta\nu_{D}}{c}J_{\infty,k}I_{k}=n_{b}h\Delta\nu_{D}H\frac{J_{\infty,k}}{J_{0}}I_{k}\leavevmode\nobreak\
,$ (5.35)

where $k$ stands for $k=c$, continuum, or $k=i$, injected photons, and
$J_{0}(z)=c\bar{n}_{b}(z)/(4\pi\nu_{\alpha})$ the intensity corresponding to
one photon per baryon [370]. The integrals $I_{c}$ and $I_{i}$ encode the
details of the scattering effects and depend on the temperature and the Gunn-
Peterson optical depth $\tau_{GP}$ (the optical depth of Ly$\alpha$ photons
redshifting through the IGM). For these quantities, we use the fits provided
by Ref. [372].

It is straightforward to estimate the efficiency of Ly$\alpha$ photons to heat
the medium. For the ranges of interest, epochs between $z\sim 30$ and $z\sim
10$, and for temperatures between $\sim 10$ K and $\sim 100$ K, the heating
integrals $|I_{k}|$ are of the order of $\sim 10$ or lower [372]. The rate is
proportional to the energy $h\Delta\nu_{D}\simeq 0.5\,{\rm
K}\,\sqrt{T_{K}/(10^{2}{\rm K})}$. Employing the estimate of
$J_{\infty}/J_{0}\sim 3\times 10^{-2}$ at $z\sim 20$ from Sec. 5.2.3.1, we
obtain a heating per baryon in a Hubble time of

$\displaystyle\frac{2\,\mathcal{Q}|_{\rm Ly\alpha,k}}{3\,H\,n_{b}}$
$\displaystyle=\frac{2}{3}h\Delta\nu_{D}\frac{J_{\infty,k}}{J_{0}}I_{k}$
(5.36) $\displaystyle\simeq 0.1\;{\rm K}\;\left(\frac{T_{K}}{10^{2}{\rm
K}}\right)^{1/2}\left(\frac{J_{\infty}/J_{0}}{3\times
10^{-2}}\right)\left(\frac{I_{k}}{10}\right)\leavevmode\nobreak\ ,$ (5.37)

Note that this is two orders of magnitude lower than the estimation from
X-rays in Eq. (5.28), and also lower than the 21 cm contribution, Eq. (5.34).
This difference even increases for lower redshifts, when X-rays are more
efficient. Moreover, the term of injected photons cool the gas rather than
heating it at these temperatures and redshifts, which counteracts in some
amount the continuum term. Therefore, Ly$\alpha$ scatterings only mildly heat
up the gas compared to absorption of X-rays, see Refs. [382, 9] for further
comparisons. As happens with the 21 cm heating rate, with realistic X-ray
luminosities this heating channel becomes subdominant and in practice (as done
in 21cmFAST) is usually neglected.

#####  5.2.2 Ionization evolution

The evolution of the local ionized fraction of the neutral IGM, $x_{e}({\bf
x},z)\equiv n_{e}/n_{b}$, is given by the interplay between ionizations of H
and He atoms, and recombinations of positive ions. Part of the discussion is
similar to that of Sec. 5.1.1, with few differences. Since the neutral IGM is
a thick medium where most of ionizing photons are fastly absorbed, the
recombinations to the ground state do not effectively account for the net
ionized fraction, and then case B recombination must be employed, rather than
case A [354]. Moreover, the ionization of the neutral IGM is produced by X-ray
photons with a mean free path long enough to escape beyond the galaxy
environment, instead of ionizing UV photons locally absorbed. For that reason,
the on-the-spot approximation is not as good estimate as in Sec. 5.1.1, and
the proper ionization rate should be computed. Considering Eq. (5.1) and the
equivalent ones for HeI and HeIII, and taking into account that conservation
of charge implies $n_{e}=n_{\rm HII}+n_{\rm HeII}+2n_{\rm HeIII}$, we can
derive the evolution equation for $x_{e}$:

$\frac{dx_{e}}{dt}=\Lambda_{\rm ion}-\alpha_{\rm
B}\,C\,x_{e}^{2}\,n_{b}\,\mathfrak{f}_{\rm H}\leavevmode\nobreak\ ,$ (5.38)

where $\alpha_{\rm B}\simeq 2.6\times
10^{-13}(T_{K}/10^{4}K)^{-0.7}cm^{3}s^{-1}$ is the case-B recombination, and
$C$ is the clumping factor (evaluated at the scale on the scale of the
simulation cell). The total ionization rate $\Lambda_{\rm ion}$ is given by
the joint contribution of the ionizing rates $\Gamma_{j}$ of each species $j$,
with $n_{b}\Lambda_{\rm ion}=\sum_{j}\Gamma_{j}n_{j}$, and can be written as

$\Lambda_{\rm ion}=\sum_{j}4\pi
x_{j}\mathfrak{f}_{j}\int^{\infty}_{\nu_{0}}d\nu\,\sigma_{j}(\nu)\,F_{j}\,J_{X}(\nu),$
(5.39)

where $J_{X}$ is the X-ray flux, given by Eq. (5.22), $F_{j}$ accounts for the
contribution from primary and secondary ionizations, and the summation goes
over the different degrees of ionization $j={\rm HI,HeI,HeII}$. In the neutral
IGM, we can assume the same degree of ionization of the species HII and HeII,
thus $n_{\rm HII}/n_{\rm H}\simeq n_{\rm HeII}/n_{\rm He}\simeq x_{e}$, the
fraction of HeIII being negligible. Therefore, in the equation above we can
write $x_{\rm HI}=x_{\rm HeI}=1-x_{e}$ and $x_{\rm HeII}=x_{e}$. The factors
$F_{j}$ read [378]:

$F_{j}=1+(h\nu-E^{th}_{j})\left(\frac{f_{\rm ion,HI}}{E^{th}_{\rm
HI}}+\frac{f_{\rm ion,HeI}}{E^{th}_{\rm HeI}}+\frac{f_{\rm
ion,HeII}}{E^{th}_{\rm HeII}}\right),$ (5.40)

being $E^{th}_{j}$ the threshold energy for each species, and $f_{\rm ion,j}$
the fraction of the electron’s energy going into secondary ionizations of
species $j$. We employ the fits derived in Ref. [571] for the functions
$f_{\rm ion,j}$, which can be roughly approximated by Eq. (5.26) for HI, being
suppressed by an order of magnitude for Helium. The unity term in the equation
above stands for primary ionizations, while the second term accounts for
secondary ones. These secondary ionizations are produced by electrons which,
after receiving energy from a X-ray photon, become very fast and non-thermal,
colliding and distributing their energy in several channels until
thermalizing.

It is possible to sketch the $x_{e}$ evolution in some approximate limits.
Before the formation of X-ray sources and after Recombination, only a leftover
fraction remained ionized. One can easily estimate the evolution of this
remnant by neglecting ionizing sources in Eq. (5.38). In this simple limit,
the differential equation can be easily integrated, obtaining a slowly
decaying function of redshift towards a constant value. For the recombination
coefficient, it can be assumed that $T_{K}\simeq T_{\gamma,0}(1+z)$ (with
$T_{\gamma,0}=2.73$ K) due to Compton cooling, valid for high redshifts (as
shown in Sec. 5.2.1.1). Since Recombination takes place at $z_{\rm R}\sim
10^{3}$, when the universe was still completely ionized, $x_{e}(z_{\rm R})\sim
1$, hence, the remnant ionized fraction for redshifts much lower than $z_{\rm
R}$ is given by

$x_{e}\simeq\frac{H_{0}\sqrt{\Omega_{m}}}{\alpha_{B}(T_{\gamma,0})n_{\rm
H,0}(1+z_{\rm R})^{0.8}}\simeq 3\times 10^{-4},$ (5.41)

which remarkably coincides with the precise estimation from Recombination
codes such as RECFAST [565].888It has to be emphasized that at $z\sim 150$,
the Compton coupling is not efficient anymore, and the kinetic temperature
deviates from the CMB one, as commented above, but taking this into account
would not alter the order-of-magnitude estimate. On the other hand, this
fraction grows again when starbust galaxies or other sources start injecting
X-ray photons to the medium. Following the same procedure of Sec. 5.2.1.2, we
can treat the ionization rate using the on-the-spot approximation, to rouhgly
estimate the impact of X-rays on $x_{e}$. Similarly to Eq. (5.27), and
neglecting secondary ionizations, we get

$\Lambda_{\rm ion}\simeq\frac{L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm
keV}]}}{{\rm SFR}}\,f_{*}\frac{\mu
m_{p}}{h\nu_{0}}\,\frac{df_{coll}(z,M_{*})}{dt}.$ (5.42)

Neglecting recombinations, we can solve Eq. (5.38), obtaining

$\begin{split}\Delta x_{e}(z)&\sim 10^{-2}\Delta
f_{coll}(z,M_{*})\left(\frac{L_{{\rm X}[2-10\leavevmode\nobreak\ {\rm
keV}]}/{\rm SFR}}{3.4\times 10^{40}erg\cdot s^{-1}M_{\odot}^{-1}yr}\right)\\\
&\times\left(\frac{f_{*}}{0.01}\right)\left(\frac{500eV}{h\nu_{0}}\right).\end{split}$
(5.43)

Note that for the fiducial values considered, the increase in $x_{e}$ would be
somewhat mild. This fact, compared to the discussion in Sec. 5.1, indicates
that the progress of reionization processes must be mostly driven by the
expansion and overlap of HII bubbles due to UV radiation, rather than the
effect of X-rays in the mostly neutral regions, which has a moderate impact.
For more quantitative computations, however, a full and precise determination
of $x_{e}$ is required and must be computed properly.

#####  5.2.3 Ly$\alpha$ flux

One of the most relevant quantities for the determination of the spin
temperature is the Ly$\alpha$ flux, $J_{\alpha}$, responsible of the WF effect
commented in Sec. 4.2.2. There are two main sources which are able to produce
a relevant flux in the IGM at the redshifts of interest: the redshifted
radiation from stellar emission, providing a flux $J_{\alpha,*}$, and X-ray
excitation, $J_{\alpha,X}$. The total flux is then given by the sum
$J_{\alpha}=J_{\alpha,*}+J_{\alpha,X}$. These mechanims are reviewed in the
following.

######  5.2.3.1 Stellar emission

Once the first stars in the universe are formed, they release radiation
background between the Lyman limit and the Ly$\alpha$ frequency. Due to the
expansion of the universe, these photons can redshift to Ly-$n$ resonances,
and eventually lead to a $2p\rightarrow 1s$ transition through a decaying
cascade. Due to the high optical depth of the IGM, photons redshifting to
Lyman resonances are always absorbed. When a photon reaches the Ly-$n$
resonance at redshift $z$, it can be emitted at $z_{max,n}>z$, which is the
maximum emission redshift so that a photon reaches the Ly-n resonance at
redshift z:

$1+z_{max,n}=(1+z)\,\frac{\nu_{n+1}}{\nu_{n}}=(1+z)\,\frac{1-\frac{1}{(n+1)^{2}}}{1-\frac{1}{n^{2}}},$
(5.44)

where the Ly-$n$ frequency is $\nu_{n}=\nu_{LL}(1-n^{-2})$, with
$\nu_{LL}=13.6$ eV$/h$ the Lyman limit frequency. As stated in Sec. 4.3.1,
some radiative cascades do not contribute to the Ly$\alpha$ flux, since they
terminate in the 2s $\rightarrow$ 1s transition. We thus define the
probability to generate a Ly$\alpha$ photon after the absorption by the $n$
level as the recycled fraction of level $n$, $f_{rec}(n)$ [377]. This quantity
can be easily computed by means of an iterative algorithm from the selection
rules and the decay rates. The total Ly$\alpha$ flux can be written in terms
of the photon comoving emissivity $\varepsilon_{\alpha}(\nu,z)$ as [377, 573]

$J_{\alpha,*}=\frac{c(1+z)^{2}}{4\pi}\sum_{n=2}^{n_{max}}f_{rec}(n)\int_{z}^{z_{max,n}}dz^{\prime}\frac{\varepsilon_{\alpha}(\nu^{\prime}_{n},z^{\prime})}{H(z^{\prime})}\leavevmode\nobreak\
,$ (5.45)

where the emission frequency is
$\nu^{\prime}_{n}=\nu_{n}(1+z^{\prime})/(1+z)$.

Following the procedure of Sec. 5.1, we can employ a simple model for the
emissivity, proportional to the comoving star formation rate density
$\dot{\rho}_{*}(z)$ [574],

$\varepsilon_{\alpha}(\nu,z)=\mathcal{N}(\nu)\frac{\dot{\rho}_{*}(z)}{\mu
m_{p}}=\mathcal{N}(\nu)f_{*}\bar{n}_{b,0}\frac{df_{coll}(z,M_{*})}{dt}\leavevmode\nobreak\
,$ (5.46)

where $\mathcal{N}(\nu)$ is the the number of emitted photons per frequency.
The radiation spectral distribution in frequencies can be taken as a piecewise
power law between each $n$ and $n+1$ Lyman levels,

$\mathcal{N}(\nu)=N_{n}\frac{(\beta_{n}+1)\nu_{\alpha}^{\beta_{n}}}{(\nu_{n+1}^{\beta_{n}+1}-\nu_{n}^{\beta_{n}+1})}\left(\frac{\nu}{\nu_{\alpha}}\right)^{\beta_{n}}\leavevmode\nobreak\
,\;\;\;\;\nu_{n}\leq\nu\leq\nu_{n+1},$ (5.47)

with $\nu_{\alpha}=(3/4)\nu_{LL}$ the Ly$\alpha$ frequency, $N_{n}$ the number
of photons emitted between the $n$ and $n+1$ resonances and $\beta_{n}$ the
spectral index, whose numerical values can be derived from stellar synthesis
models [575, 574]. The function above is normalized as
$\int_{\nu_{n}}^{\nu_{n+1}}d\nu\mathcal{N}(\nu)=N_{n}$, while the
normalization of the total spectrum gives the total number of photons emitted
between the Ly$\alpha$ and the Lyman limit
$\int_{\nu_{\alpha}}^{\nu_{LL}}d\nu\mathcal{N}(\nu)=\sum_{n}N_{n}=N_{\rm
tot}$. For Population II (or low-metallicity) stars, $N_{\rm tot}\simeq 9690$,
out of which 6520 photons are emitted between Ly$\alpha$ and Ly$\beta$ lines
[575, 574]. For very massive and metal-free Population III stars, $N_{\rm
tot}\simeq 4800$, with 2670 photons between Ly$\alpha$ and Ly$\beta$ [495].
However, given the lack of observations of high-redshifts galaxies, the mass
distribution, metallicity and therefore the spectrum of these sources is still
highly uncertain. Hence, $N_{\rm tot}$ may be considered as a free parameter,
which is usually proportional to the number of ionizing photons
$N_{\gamma/b}$. In 21cmFAST, the fiducial scenario assumes Population II
stars, which have $N_{\gamma/b}=4\times 10^{3}$, that, with the spectrum of
Eq. (5.47), corresponds to a value of $N_{tot}\simeq 9690$.

As a side note, an emissivity proportional to $\nu^{-1}$ (normalized as
$\mathcal{N}(\nu)=N_{tot}/{\rm ln}(4/3)\,\nu^{-1}$) allows to solve
analytically the integral of Eq. (5.45), obtaining the following Ly$\alpha$
flux [9]:

$J_{\alpha,*}(z)=\frac{c\bar{n}_{b}(z)f_{*}N_{tot}}{4\pi\nu_{\alpha}ln(4/3)}\sum_{n=2}^{n_{max}}f_{rec}(n)\frac{\nu_{\alpha}}{\nu_{n}}\left(f_{coll}(z,M_{*})-f_{coll}(z_{max,n},M_{*})\right)\leavevmode\nobreak\
.$ (5.48)

This flux differs by less than 10% from the one computed with the piecewise
power-law spectrum presented above, and given its analytical form, can provide
insight about the flux evolution. Profiting from the formula above, it is
possible to estimate the production of Ly$\alpha$ photons in the high redshift
universe. Writing the intensity corresponding to one photon per baryon
$J_{0}(z)=c\bar{n}_{b}(z)/(4\pi\nu_{\alpha})$ as a reference value [370], it
is is possible to give an order-of-magnitude approximation for the Ly$\alpha$
flux as

$J_{\alpha,*}(z)\sim 3\times
10^{-2}J_{0}(z)\left(\frac{N_{tot}}{9690}\right)\left(\frac{f_{*}}{0.01}\right)\left(\frac{f_{coll}(z,M_{*})}{10^{-4}}\right),$
(5.49)

where the fiducial value taken for $N_{tot}$ corresponds to Pop II stars, as
commented above, and $f_{coll}(z,M_{*})\simeq 10^{-4}$ for a minimum virial
temperature of $10^{4}$ K at $z\simeq 20$ (see Fig. 2.1). With this result,
the Ly$\alpha$ coupling coefficient from Eq. (4.46) reads

$\displaystyle x_{\alpha}$ $\displaystyle\simeq 14.5\times
S_{\alpha}\frac{J_{\alpha,*}(z)}{J_{0}(z)}\left(\frac{1+z}{21}\right)^{2}$
(5.50) $\displaystyle\sim
0.5\;S_{\alpha}\left(\frac{1+z}{21}\right)^{2}\left(\frac{N_{tot}}{9690}\right)\left(\frac{f_{*}}{0.01}\right)\left(\frac{f_{coll}(z,M_{*})}{10^{-4}}\right),$
(5.51)

where, as stated in Sec. 4.3.2, $S_{\alpha}$ is an order 1 function accounting
for the details of radiative transfer. Therefore, we see that, for these
fiducial astrophysical parameters, the WF coefficient approaches
$x_{\alpha}\sim 1$ at $z\simeq 20$, implying an efficient WF effect and thus a
strong coupling between the spin temperature and the kinetic one from there
on.

######  5.2.3.2 X-ray excitation

Similarly to what happens with secondary ionizations, absorption of energetic
X-ray photons can lead to collisional excitation by fast electrons of Hydrogen
atoms, producing a Ly$\alpha$ photon at the de-excitation [384]. The
corresponding flux of this mechanism can be written as

$J_{\alpha,X}=\frac{cn_{b}(z)}{4\pi
h\nu_{\alpha}H\nu_{\alpha}}4\pi\int_{\nu_{0}}^{\infty}d\nu\,J_{\rm
X}(z,\nu)\,\sum_{j}(h\nu-E_{j}^{\rm th})f_{\rm
Ly\alpha}\,f_{j}\,x_{j}\,\sigma_{j}(\nu),$ (5.52)

where $f_{\rm Ly\alpha}$ is the fraction of electron energy which is deposited
as Ly$\alpha$ excitations, whose value is close to the total excitations
fraction, $f_{\rm exc}$, as $f_{\rm Ly\alpha}\simeq 0.8f_{\rm exc}$ [569]. In
the high energy limit, where $f_{\rm exc}$ saturates to constant values of
energy, the expression above can be approximated as

$J_{\alpha,X}\simeq\frac{c}{4\pi h\nu_{\alpha}H\nu_{\alpha}}\frac{f_{\rm
Ly\alpha}}{f_{\rm heat}}\mathcal{Q}|_{X}(z)=J_{0}(z)\frac{f_{\rm
Ly\alpha}}{f_{\rm heat}}\frac{\mathcal{Q}|_{X}(z)}{Hh\nu_{\alpha}n_{b}(z)}$
(5.53)

with $f_{\rm Ly\alpha}\simeq f_{\rm heat}$ at high energies [569]. As we have
done for the heating and ionization in the previous subsections, we can apply