# The 21-cm signal during the end stages of reionization

Sambit K. Giri,1, Michele Bianco2,3,5, Timothée Schaeffer4, Ilian T. Iliev5,
Garrelt Mellema6, and Aurel Schneider4
1 Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes
Alfv́ens väg 12, SE-106 91 Stockholm, Sweden
2 Laboratoire d’Astrophysique, Ecole Polytechnique Federale de Lausanne
(EPFL), Observatoire de Sauverny, Versoix 1290, Switzerland
3 Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-
Str 27, 8093 Zurich, Switzerland
4 Department of Astrophysics, University of Zurich, Winterthurerstrasse 190,
8057 Zurich, Switzerland.
5 Astronomy Centre, Department of Physics & Astronomy, Pevensey III Building,
University of Sussex, Falmer, Brighton, BN1 9QH, United Kingdom
6 Department of Astronomy and Oskar Klein Centre, AlbaNova, Stockholm
University, SE-10691 Stockholm, Sweden E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ, Report Number NORDITA
2024-003)

###### Abstract

During the epoch of reionization (EoR), the 21-cm signal allows direct
observation of the neutral hydrogen (H  I ) in the intergalactic medium (IGM).
In the post-reionization era, this signal instead probes H  I in galaxies,
which traces the dark matter density distribution. With new numerical
simulations, we investigated the end stages of reionization to elucidate the
transition of our Universe into the post-reionization era. Our models are
consistent with the latest high-redshift measurements, including ultraviolet
(UV) luminosity functions. Notably, these models consistently reproduced the
evolution of the UV photon background, which is constrained from
Lyman-$\alpha$ absorption spectra. We studied the dependence of this
background on the nature of photon sinks in the IGM, requiring mean free path
of UV photons to be $\sim$10 comoving-megaparsecs (cMpc) during the EoR that
increases gradually with time during late stages ($z\lesssim 6$). Our models
revealed that the reionization of the IGM transitioned from an inside-out to
an outside-in process when the Universe is less than 0.01 per cent neutral.
During this epoch, the 21-cm signal also shifted from probing predominantly
the H  I in the IGM to that in galaxies. Furthermore, we identified a
statistically significant number of large neutral islands (with sizes up to 40
cMpc) persisting until very late stages ($5\lesssim z\lesssim 6$) that can
imprint features in Lyman-$\alpha$ absorption spectra and also produce a knee-
like feature in the 21-cm power spectrum.

###### keywords:

techniques: interferometric, cosmology: theory, reionization, first stars,
early Universe, radio lines: galaxies

††pubyear: 2024††pagerange: The 21-cm signal during the end stages of
reionization–A

##  1 Introduction

Recent cutting-edge observational facilities, such as the James Webb Space
Telescope (JWST), have substantially improved our capabilities to explore the
high-redshift ($z\gtrsim 5$) Universe. In recent years, the dataset of early
galaxies and quasars has been quickly growing (e.g. Oesch et al., 2016; Atek
et al., 2018, 2023; Bakx et al., 2023; Robertson et al., 2023; Bunker et al.,
2023). These observations are enhancing our understanding of both cosmology
and astrophysics during the epoch of reionization (EoR), which is the period
when early photon sources formed and ionized the gas in the intergalactic
medium (IGM)(e.g. Sun et al., 2023; Schaeffer et al., 2023; Hassan et al.,
2023; Dayal & Giri, 2024; Lovell, 2024; Yung et al., 2024). We refer
interested readers to Morales & Wyithe (2010) and Dayal & Ferrara (2018) for
recent reviews. This work focuses on the end stages of reionization when the
IGM has been predominately heated and ionized.

Previous models of reionization ended at $z\approx 6$, which was a long-
standing observational constraint (e.g. Fan et al., 2006). Current
observations updated our understanding of the end stages of reionization,
including hints that it prolonged below $z\approx 6$ (e.g., Kulkarni et al.,
2019; Bosman et al., 2022). Observations have also suggested a substantial
evolution in global quantities, such as the mean ionizing (UV) photon
background (e.g., Wyithe & Bolton, 2011; Calverley et al., 2011) and the mean
free path (MFP) of ionizing photons (e.g., Worseck et al., 2014; Becker et
al., 2021), which has been challenging to reproduce in simulations. Recent
studies have developed models focused on understanding these quantities (e.g.,
D’Aloisio et al., 2018; Cain et al., 2021; Puchwein et al., 2019; Chardin et
al., 2016). The models generated in this work are consistent with these latest
measurements during the late stages.

Hydrogen, the most abundant element in our Universe, can be probed by using
the neutral hydrogen (H  I ) 21-cm line emission, which can be measured with
radio telescopes (e.g., Pritchard & Loeb, 2012). During the EoR, this signal
can probe the evolution of the IGM, which depends on cosmological structure
formation and the astrophysical processes of galaxy formation (e.g. Iliev et
al., 2012; Dixon et al., 2016; Schaeffer et al., 2023; Schneider et al.,
2023). Previous studies have demonstrated its capability to constrain not only
the properties of these early galaxies (e.g., Greig & Mesinger, 2015; Mirocha
& Furlanetto, 2019; Park et al., 2019; Qin et al., 2021a) but also aspects of
cosmology (e.g., Sitwell et al., 2014; Kern et al., 2017; Schneider, 2018; Liu
& Slatyer, 2018; Nebrin et al., 2019; Giri & Schneider, 2022).

In the post-reionization era, most H  I resides inside galaxies, shielded from
UV photons (e.g., Villaescusa-Navarro et al., 2014, 2018). Unlike during the
EoR, the 21-cm signal during this phase traces the large-scale distribution of
galaxies. Therefore, this signal serves as a tracer of the large-scale
structure of our Universe and has been proposed as a new cosmological probe
(e.g. Camera et al., 2013; Xu et al., 2016; Carucci et al., 2017; Wu & Zhang,
2022). However, during the end stages of reionization, the 21-cm signal traces
both the H  I in galaxies and the IGM (e.g. Wyithe & Loeb, 2008; Xu et al.,
2019). Previous studies have explored this using semi-analytical frameworks
and found that the residual H  I in the IGM can impact the signal at high-
density regions (Miralda-Escudé et al., 2000; Wyithe & Loeb, 2008). These
studies inferred an ‘outside-in’ nature of reionization during the end stages.
This suggests the ionization fraction is anti-correlated with matter density
due to higher recombination rates in dense regions. We will test this
hypothesis using fully numerical simulations.

Ongoing radio experiments, such as the Low-frequency Array (LOFAR; Mertens et
al., 2020), Murchison Wide-field Array (MWA; Trott et al., 2020) and Hydrogen
Reionization Array (HERA; The HERA Collaboration et al., 2022b, 2023), are
improving the upper limits on the power spectrum of the 21-cm signal. These
observations have provided valuable constraints on the early Universe physics
(e.g. Ghara et al., 2020, 2021; Mondal et al., 2020; Greig et al., 2021b;
Greig et al., 2021a; The HERA Collaboration et al., 2022a). Among these
measurements, Trott et al. (2020) went down to $z\approx 6.5$, closest to the
end stages of reionization. Additionally, the 21-cm signal has been detected
during the post-reionization era in cross-correlation (The CHIME Collaboration
et al., 2023). Cross-correlating the post-reionization 21-cm signal with other
tracers of cosmological matter distribution can further constrain cosmology
(e.g., Padmanabhan et al., 2020).

In the near future, we expect to detect the post reionization 21-cm power
spectrum with advanced radio experiments, including Baryon Acoustic
Oscillations In Neutral Gas Observations (BINGO; Abdalla et al., 2022), Five-
hundred-meter Aperture Spherical radio Telescope (FAST; Bigot-Sazy et al.,
2015), and Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX;
Newburgh et al., 2016). The Square Kilometre Array (SKA) is currently under
construction and will be powerful enough to go beyond the power spectrum,
providing images of the H  I distribution (Mellema et al., 2015). SKA will
consist of two frequency components, SKA-Mid and SKA-Low, focusing on
observing post-reionization and the EoR, respectively. The initial phase of
SKA-Low will cover a frequency range from about 50 to 350 MHz111The latest
details about the SKA can be found at
https://www.skao.int/en/explore/telescopes., corresponding to redshifts
between 30 and 3 (Koopmans et al., 2015). Thus, this component will be
instrumental in probing the end stages of reionization. The suite of
simulations developed in this work will play a vital role in understanding
these measurements.

Modelling the IGM during these end stages is challenging because almost all
the ionized bubbles are overlapping, making large distances in the IGM
transparent to UV photons (e.g. Xu et al., 2017; Giri et al., 2019).
Additionally, we need to model the astrophysical processes inside galaxies to
study the amount of H  I shielded from UV photons. In principle, a hydro-
dynamical N-body framework with radiative transfer calculation would be
capable of fully depicting the transition from EoR to Post-EoR. However, such
simulations are computationally challenging due to the large dynamical range
required to study the H  I in the IGM and inside galaxies. Computationally
cheap semi-numerical frameworks become inaccurate during these times as bubble
overlaps lead to photon non-conservation (e.g. Choudhury & Paranjape, 2018).
In Giri et al. (2019), we developed a suite of large-scale radiative transfer
simulations to study the IGM during these end stages. We will extend this
exploration with a new simulation suite tuned to the latest measurements and
including the post-reionization H  I distribution.

This paper is structured as follows: in the next section, we describe our
simulation framework. We present our models in Section 3 and discuss the
simulated 21-cm signal in Section 4. Finally, we summarise our findings in
Section 5.

##  2 Simulation framework

We first describe the $N$-body simulation used in this study to model the
photon sources during the epoch of reionization. To explore the evolution of
the 21-cm signal during the end stages of reionization, we require modelling
of both the H  I distribution in the IGM and the H  I inside galaxies that are
self-shielded from UV photons. Therefore, in subsection 2.2, we define the
sources of reionization and the method to model unresolved sinks in the IGM.
Finally, in subsection 2.3, we explain our method for modelling the H  I
content inside galaxies.

###  2.1 Cosmological structure formation

We model the cosmological structure formation by running an $N$-body code,
Pkdgrav3 (Potter et al., 2017), with $2048^{3}$ particles in a simulation box
of length $298$ cMpc in each direction. We chose this box size to reduce the
sample variance at large scales (i.e. $k\approx 0.1\,\mathrm{Mpc}^{-1}$; Iliev
et al., 2014; Giri et al., 2023), where the current radio experiments provided
the best upper limits. This simulation setup gives us a particle mass of
$m_{\rm part}=1.2\times 10^{8}M_{\odot}$. We identify dark matter haloes using
a friends-of-friend halo finder implemented in Pkdgrav3 with a minimum of 10
dark matter particles, which helped resolve haloes down to the high-mass
atomically cooling haloes (HMACHs) defined in Dixon et al. (2016). Though this
minimum number of particles is particularly low, we found a good match with
the halo mass function calculated from extended Press-Schechter formalism. We
show this validation test in Appendix A. We assume that star formation in
haloes smaller than HMACHs is suppressed, especially during the end stages of
reionization, due to radiative feedback (e.g. Nebrin et al., 2023).

Throughout this study, we have considered a flat cold dark matter cosmological
model with parameters aligned with the Planck 2018 results (Planck
Collaboration et al., 2020). These parameters include a matter abundance of
$\Omega_{\rm m}=0.32$, a baryon abundance of $\Omega_{\rm b}=0.044$, a
dimensionless Hubble constant of $h=0.67$, and a standard deviation of matter
perturbations at the 8$h^{-1}$ cMpc scale denoted as $\sigma_{\rm 8}=0.83$. We
consider a primordial gas with Helium abundance mass factor of $Y_{p}=0.249$,
which gives a molecular weight of $\mu=(1-Y_{p})^{-1}=1.33$. We initialised
our $N$-body simulation at $z=150$ with second-order Lagrangian perturbation
theory (e.g. Bertschinger, 1998) using a transfer function generated with the
Boltzmann code, class (Lesgourgues, 2011; Blas et al., 2011). The $N$-body
simulation was run down to $z\approx 4$ while saving snapshots at every
$\sim$10.86 Myrs. We used the Piecewise Cubic Spline mass assignment scheme
(Sefusatti et al., 2016) to generate matter density fields on $256^{3}$ grids.

###  2.2 Reionization of the intergalactic medium

We simulate the evolution of the IGM using our state-of-the-art radiative
transfer simulation code pyC2Ray (Hirling, Bianco et al. 2023), which is an
upgraded version of C2Ray (Mellema et al., 2006). This code has a user-
friendly Python interface and it can leverage Graphics processing units (GPUs)
to reduce the computing time of solving the three-dimensional radiative
transfer equation in a cosmological simulation volume.

To leverage the latest measurements of the high-redshift sources, we have
implemented a flexible parameterization designed to identify viable ionizing
source models. We detail this parameterization in subsection 2.2.1. During
late times, the ionization state of the IGM may depend on small-scale
absorbers. In subsection 2.2.2, we describe how we deal with the absorption of
photons by unresolved density fluctuations.

####  2.2.1 Source model

We populate the dark matter haloes provided by our $N$-body simulation with
galaxies of stellar mass
$M_{\star}=f_{\star}(M_{h})(\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}})M_{h}$,
where $f_{\star}$ describes the fraction of baryonic mass falling into a halo
of mass $M_{h}$ that gets converted to stars. We parameterize this star
formation inside the haloes using the following stellar-to-halo relation,

$\displaystyle
f_{\star}(M_{\mathrm{h}})\equiv\frac{M_{\star}}{M_{\mathrm{h}}}=f_{\star,0}\left(\frac{M_{\mathrm{h}}}{10^{10}M_{\odot}}\right)^{\alpha_{\star}}\
,$ (1)

where $f_{\star,0}$ and $\alpha_{\star}$ are the normalization constant and
index of the power-law, respectively. We should note that our simulations
framework, pyC2Ray, implements a more generic parametrization provided in
Schneider et al. (2021, 2023), which is needed for studying non-cold dark
matter models (Giri & Schneider, 2022). In the current work, we consider the
reduced form given above that is enough to model the ultraviolet luminosity
functions (UVLFs; e.g. Gillet et al., 2020; Park et al., 2020).

We require the photons produced inside the sources to model the gas properties
in the IGM. We model the rate of the number of UV photons ($\dot{N}_{\gamma}$)
escaping into the IGM with the following relation,

$\displaystyle\dot{N}_{\gamma}(M_{h})=f_{\mathrm{esc}}(M_{h})\zeta_{\mathrm{ion}}\dot{M}_{\star}(M_{h})\,$
(2)

where $f_{\mathrm{esc}}$, $\zeta_{\mathrm{ion}}$ and $\dot{M}_{\star}$ are the
escape fraction, ionizing photon production efficiency and stellar mass
growth, respectively. $\zeta_{\mathrm{ion}}$ is defined as the number of
ionizing photons produced inside a source per unit proton mass. The mass
dependence of the escape fraction can be parameterized with a power-law given
as (e.g. Park et al., 2019),

$\displaystyle
f_{\mathrm{esc}}(M_{h})=f_{\mathrm{esc,0}}\left(\frac{M_{h}}{10^{10}M_{\odot}}\right)^{\alpha_{\mathrm{esc}}}\
,$ (3)

where $f_{\mathrm{esc,0}}$ and $\alpha_{\mathrm{esc}}$ are free parameters.

We model the star formation rate as
$\dot{M}_{\star}=f_{\star}(M_{h})(\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}})\dot{M}_{h}$,
where $\dot{M}_{h}$ is the halo accretion rate. We assume an exponential mass
growth model222A similar accretion model is implemented in the 21cmFAST
framework, but their timescales differ by a factor of $(1+z)$. See eq. 3 in
Park et al. (2019) for comparison. See Schneider et al. (2021) for a
comparison of several mass growth models. where the accretion rate is given
as,

$\displaystyle\dot{M}_{h}(z)=\frac{M_{h}(z)}{[\alpha_{h}(1+z)H(z)]^{-1}}\ .$
(4)

$H(z)$ is the Hubble parameter and $\alpha_{h}$ is a free parameter. We set
$\alpha_{h}=0.79$ to obtain a good match to detailed simulations from Behroozi
et al. (2020). This mass growth model incorporates both the gradual halo
growth and mergers (e.g. McBride et al., 2009; Dekel et al., 2013; Trac et
al., 2015).

For modelling the reionization of the IGM, pyC2Ray reads the halo catalogue
from the $N$-body simulation and calculates the corresponding
$\dot{N}_{\gamma}$ using Eq. 2. This quantity is then assigned to the relevant
grid cells in the gridded matter density fields. The transfer of radiation in
the IGM is modelled keeping $\dot{N}_{\gamma}$ fixed between successive time
steps of the $N$-body snapshots.

####  2.2.2 Sink model

Large-scale simulation volumes such as ours cannot resolve the small-scale
density structures relevant for absorbing UV photons. Several physical
interpretations exist of these unresolved sinks, including self-shielded
systems (i.e. Lyman-limit systems) and small-scale density fluctuations that
can increase recombination. Modelling these structures is crucial for this
work as the IGM in the late stages are sensitive to the properties of the
sinks (e.g. Cain et al., 2021; Gaikwad et al., 2023). Previous authors
exploring this problem have suggested different methods to model them (e.g.
Miralda-Escudé et al., 2000; Ciardi et al., 2006; Choudhury et al., 2009;
Sobacchi & Mesinger, 2014; Shukla et al., 2016; Mao et al., 2020; Bianco et
al., 2021a).

We have implemented multiple approaches to account for the effect of
unresolved sinks in pyC2Ray, such as through the clumping factor $\mathcal{C}$
and mean free path $R_{\mathrm{mfp}}$. The former approach boosts the
recombinations by increasing the ratio of variance of the matter density and
the square of the mean density. The latter approach defines a maximum allowed
distance for the photons to travel in the IGM, assuming these photons get
absorbed by unresolved sinks. We should note that the ionized regions in the
IGM can have sizes larger than $R_{\mathrm{mfp}}$ due to ionized bubble
overlaps.

In our framework, we have implemented a redshift evolving $R_{\mathrm{mfp}}$,
which is parameterised as,

$\displaystyle
R_{\mathrm{mfp}}(z)=A_{\mathrm{mfp}}\left(\frac{1+z}{5}\right)^{\eta_{\mathrm{mfp,0}}}\left[1+\left(\frac{1+z}{1+z_{\mathrm{mfp}}}\right)^{\eta_{\mathrm{mfp,1}}}\right]^{\eta_{\mathrm{mfp,2}}}\
,$ (5)

where $A_{\mathrm{mfp}}$, $z_{\mathrm{mfp}}$, $\eta_{\mathrm{mfp,0}}$,
$\eta_{\mathrm{mfp,1}}$ and $\eta_{\mathrm{mfp,2}}$ are free parameters. The
above expression is a modified form of the fit provided in Worseck et al.
(2014), which can be reconstructed by setting $\eta_{\mathrm{mfp,2}}=0$. A
non-evolving $R_{\mathrm{mfp}}$ can be modelled by setting both
$\eta_{\mathrm{mfp,1}}$ and $\eta_{\mathrm{mfp,2}}$ to zero. We introduced
this modification to incorporate recent measurements of $R_{\mathrm{mfp}}$,
which will be presented in Section 3.1.3.

###  2.3 Neutral hydrogen inside haloes

After the reionization of the IGM, the 21-cm signal is produced by H  I that
persists within massive galaxies, having self-shielded against the UV
background. Direct observation of this H  I in galaxies beyond $z\sim 0.4$
through 21-cm emission lines is quite challenging, and thus, observing the
integrated 21-cm line of unresolved galaxies over a wide sky area can provide
a more efficient and cost-effective probe for cosmological matter distribution
(Chang et al., 2010; Bagla et al., 2010). Conversely, simulating the H  I
content in galaxies is extremely complex, requiring high-resolution hydro-
dynamical N-body simulations coupled with radiative feedback and cooling
mechanisms (Bharadwaj & Srikant, 2004; Bird et al., 2014). However, a
reasonable approximation of the averaged H  I content can be obtained based on
the hosting halo mass, $M_{\mathrm{h}}$, and several models have been proposed
in the literature (see, e.g., Villaescusa-Navarro et al., 2014, 2018;
Padmanabhan et al., 2017; Modi et al., 2019; Spinelli et al., 2020).

Therefore, we adopt a straightforward approach and assign the mass of H  I
($M_{\mathrm{HI}}$) to the haloes in our simulation using the
$M_{\mathrm{HI}}-M_{\mathrm{h}}$ relation from Padmanabhan et al. (2017),
given as

$\displaystyle
M_{\text{HI}}(M_{\mathrm{h}})=\alpha\,f_{H,c}\,M_{\mathrm{h}}\,\left(\frac{M_{\mathrm{h}}}{10^{11}h^{-1}M}\right)^{\beta}\exp\left[-\left(\frac{v_{c,0}}{v_{c}(M_{\mathrm{h}})}\right)^{3}\right]\
.$ (6)

Here, $\alpha$ is a normalisation factor and represents the fraction of H  I ,
relative to cosmic hydrogen mass fraction,
$f_{\mathrm{H,c}}=(1-Y_{p})\Omega_{b}/\Omega_{m}$, associated with a dark
matter halo of mass $M_{\mathrm{h}}$. $\beta$ is the logarithmic slope of the
$M_{\mathrm{HI}}-M_{\mathrm{h}}$ correlation. $v_{c,0}$ is the virial velocity
below which the neutral hydrogen is supposed to be suppressed. We choose
$\alpha=0.9$, $\beta=-0.58$ and $v_{c,0}=10^{1.56}\,\rm km/s$, as used by
Padmanabhan et al. (2017), which was based on constraints from observational
data of the abundance and clustering of H  I in galaxies between redshift
$z\sim 0$ and $5$. We calculate the virial velocity of the halos, $v_{c}(M)$,
based on the relation:

$v_{c}(M_{\mathrm{h}})=96.6\,{\rm
km\,s^{-1}}\,\left(\frac{\Delta_{c}\,\Omega_{m}}{54.4}\right)^{\frac{1}{6}}\,\left(\frac{1+z}{3.3}\right)^{\frac{1}{2}}\left(\frac{M_{\mathrm{h}}}{10^{11}\,M_{\odot}}\right)^{\frac{1}{3}}$
(7)

where $\Delta_{c}=178$ is the mean over-density of the halo (Maller & Bullock,
2004; Barnes & Haehnelt, 2014). This paper will explore the evolution of the
21-cm signal using gridded simulation volumes conducted on a mesh with a low
spatial resolution ($\Delta r\sim 1.17$ cMpc). Therefore, after assigning a
value of $M_{\mathrm{HI}}$ to the haloes at different redshifts, we place them
on the same grids of the matter density distribution used to model the
reionization of the IGM.

##  3 Reionization models

Table 1: A list of reionization models along with the corresponding parameters. These parameters are grouped into source and sink parameters. We also mention the Thompson optical depth ($\tau_{\mathrm{e}}$) for our models that are consistent with the Planck 2018 constraints (Planck Collaboration et al., 2020). Model Name | Source | Sink | Thompson optical depth
---|---|---|---
| $\zeta_{\mathrm{ion}}\,f_{\star,0}\,f_{\mathrm{esc,0}}$ | $\alpha_{\star}$ | $\alpha_{\mathrm{esc}}$ | $A_{\mathrm{mfp}}$ | $z_{\mathrm{mfp}}$ | $\eta_{\mathrm{mfp,0}}$ | $\eta_{\mathrm{mfp,1}}$ | $\eta_{\mathrm{mfp,2}}$ | $\tau_{\mathrm{e}}$
Source1_SinkA | 4 | 0.3 | 0.0 | 210 | 6.0 | -9 | 9 | 1 | 0.049
Source2_SinkA | 4 | 0.3 | -0.3 | 210 | 6.0 | -9 | 9 | 1 | 0.051
Source3_SinkA | 4 | 0.3 | -0.5 | 210 | 6.0 | -9 | 9 | 1 | 0.053
Source1_SinkB | 4 | 0.3 | 0.0 | 210 | 5.5 | -9 | 9 | 1 | 0.049
Figure 1: The ultraviolet luminosity function (UVLF) gives the number density
of photon sources of different brightness ($M_{\mathrm{UV}}$). We plot several
measurements (Atek et al., 2018; Bouwens et al., 2021; Bowler et al., 2020;
Donnan et al., 2023; Ishigaki et al., 2018; Livermore et al., 2017; McLure et
al., 2013; Rojas-Ruiz et al., 2020; Stefanon et al., 2019) binned at three
different redshifts ($z=6,7$ and 8). The solid line represents the best-fit
model of the stellar-to-halo relation (Eq. 1) against these measurements.

The reionization of the IGM is less well constrained and, therefore, we
develop multiple models with properties consistent with available
measurements. First, we describe the source and sink parameter values of our
cosmic reionization models in subsection 3.1. In subsection 3.2, we describe
the reionization history of the IGM in our models. Later, in subsection 3.3
and 3.4, we study the evolution of the topology of neutral islands
distribution in the IGM and the UV background, respectively, in these models.

###  3.1 Model parameters

Here we will use the latest observations to find relevant UV photon source and
sink model parameters.

####  3.1.1 Stellar mass

We begin by finding suitable intrinsic source properties for our simulations
by using the ultraviolet luminosity function (UVLF) measurements, which give
the number density of sources at various UV brightness. We follow Park et al.
(2019) to express this quantity as

$\displaystyle\phi_{\rm
UV}=f_{\mathrm{duty}}(M_{h})\frac{\mathrm{d}n}{\mathrm{d}M_{h}}\frac{\mathrm{d}M_{h}}{\mathrm{d}M_{\rm
UV}}\ ,$ (8)

where $({\mathrm{d}n}/{\mathrm{d}M_{h}})$ and $M_{\mathrm{UV}}$ are the halo
mass function (HMF) and the absolute UV magnitude, respectively.
$f_{\mathrm{duty}}=\exp(-M_{\mathrm{turn}}/M_{h})$ is the duty cycle that
models the suppression of star formation in haloes below a mass of
$M_{\mathrm{turn}}$. $M_{\mathrm{UV}}$ quantifies the brightness of sources
that are given as

$M_{\mathrm{UV}}=M_{0}-2.5\biggr{[}\log_{\mathrm{10}}f_{\star}(M_{h},z)+\log_{\mathrm{10}}\frac{\Omega_{b}}{\Omega_{m}}\\\
+\log_{\mathrm{10}}\frac{\dot{M}_{h}(M_{h},z)}{\rm
M_{\odot}yr^{-1}}-\log_{\mathrm{10}}\frac{\kappa}{\rm
M_{\odot}yr^{-1}}\biggr{]}\ ,$ (9)

where $M_{0}=51.6$ corresponding to AB magnitude system (Oke, 1974). We use
$\kappa=1.15\times 10^{-28}~{}{\rm M_{\odot}yr^{-1}/(erg~{}s^{-1}Hz^{-1}})$,
calibrated for 1500 Å dust-corrected rest-frame UV luminosity. This
calibration assumes continuous star formation and a Salpeter stellar initial
mass function (Madau & Dickinson, 2014).

Figure 1 shows measurements of UVLFs that are binned at three redshifts
($z=6$, 7 and 8). Using values of $f_{\star,0}=0.1$ and $\alpha_{\star}=0.3$,
we find good fits to these measurements and adopt these values for our source
model. The truncation at the faint end is caused by the duty cycle with
$M_{\mathrm{turn}}=10^{9}M_{\odot}$, which is set by the smallest halo mass in
our simulations. We should note that $f_{\star,0}$ degenerates with
$\zeta_{\mathrm{ion}}$, which is currently loosely constrained by
measurements. The parameter $\alpha_{\star}$ models the mass-dependence of our
source model and allows the construction of star formation models in which
either low-mass or high-mass halos dominate the process. In all our source
models, we keep the parameters of the $f_{\star}(M_{h})$ relation fixed as the
UV photon budget will include the escape fraction model, which will be
discussed next.

####  3.1.2 ionizing photon budget

Figure 2: The Ly$\alpha$ escape fraction plotted against various stellar
masses. We plot three models (black lines) with different mass dependence of
the escape fraction that are consistent with estimations from CEERS (circles),
JADES (diamonds) and FRESCO (stars) survey data (Lin et al., 2024; Napolitano
et al., 2024). We should note that these data are provided at different
redshifts between $4.5\lesssim z\lesssim 8.5$. We estimate the UV escape
fraction through the relation $f_{\mathrm{esc}}\approx
0.15f_{\mathrm{esc,Ly\alpha}}$ (Begley et al., 2024).

In our framework, the total rate of escaping UV photons from a source is
defined by Eq. 2. While the star formation rate is well constrained by
observations of the UVLF (discussed in the previous subsection), the photon
escape rate is much harder to constrain. We have some indirect constraints on
the escape fraction of UV photons through the Ly$\alpha$ photon escape, given
the strong correlation between Ly$\alpha$ and UV escape (Chisholm et al.,
2018). For the comparison here, we assumed a relation for UV escape fraction
$f_{\mathrm{esc}}\approx 0.15f_{\mathrm{esc,Ly\alpha}}$ (Begley et al., 2024),
where $f_{\mathrm{esc,Ly\alpha}}$ is the Ly$\alpha$ escape fraction.

The Ly$\alpha$ photon escape can be estimated using the features in the high
redshift spectroscopic data, such as the UV slope and line strengths (e.g.
Zackrisson et al., 2013, 2017; Jensen et al., 2016; Jaskot et al., 2019; Giri
et al., 2020b; Begley et al., 2022). In Figure 2, we present some of the
latest constraints on $f_{\mathrm{esc,Ly\alpha}}$ at several stellar masses
$M_{\star}$. Though these measurements have a wide redshift range
($4.5\lesssim z\lesssim 8.5$), we consider them together. With the
$M_{\star}(M_{h})$ relation fixed against the UVLFs, we determine three models
for the UV escape fraction, $f_{\mathrm{esc}}(M_{h})$. We first set
$f_{\mathrm{esc,0}}$ to 0.02, and consider models with $\alpha_{\mathrm{esc}}$
varied as 0 (Source1), -0.3 (Source2) and -0.5 (Source3). These models yield
three different mass dependencies for the UV escape fraction. We list the
above three distinct source models as the first three entries in Table 1.

In all cases, the intrinsic ionizing efficiencies of the sources denoted as
$\zeta_{\mathrm{ion}}$, are kept fixed at 2000, which gives
$\zeta_{\mathrm{ion}}f_{\mathrm{\star,0}}f_{\mathrm{esc,0}}=4$. We made this
choice to obtain a consistent reionization history, which will be shown in
Section 3.2. Notably, the three UV source models, $\dot{N}_{\gamma}(M_{h})$,
exhibit significantly different mass dependencies. The final mass dependence
of $\dot{N}_{\gamma}$ is given by $\alpha_{\star}+\alpha_{\mathrm{esc}}$.
While the source model in ‘Source2’ is independent of halo mass, the
$\dot{N}_{\gamma}$ in ‘Source1’ and ‘Source3’ are correlated and anti-
correlated, respectively, with halo masses. We will study the impact of these
mass dependencies in the next section. The fourth model listed in Table 1 has
the same source property as ‘Source1’ but with a different sink model. We will
discuss the different sink models in the following subsection.

####  3.1.3 Unresolved absorbers

Figure 3: The mean free path models (depicted as lines) compared against
various measurements (Worseck et al., 2014; Becker et al., 2021; Zhu et al.,
2023; Davies et al., 2023). The fit provided in Worseck et al. (2014), shown
with a dashed line, does not agree with recent constraints. The ‘SinkA’ and
‘SinkB’ models considered here using the modified parameterization (Eq. 5)
have a better agreement with these constraints and converge to 10 and 20 cMpc
at high redshifts.

In this study, we set $\mathcal{C}=1$ and explore the impact of sink models
using the redshift dependent evolution of $R_{\mathrm{mfp}}$ defined in Eq. 5.
Previous studies have used the fit from Worseck et al. (2014) in their
reionization simulations (e.g. Shukla et al., 2016; Choudhury et al., 2021).
We can obtain this fit by setting $A_{\mathrm{mfp}}=175$ cMpc,
$\eta_{\mathrm{mfp,0}}=-4.4$ and $\eta_{\mathrm{mfp,2}}=0$, which match the
measurements at $z\approx 5$. However, Figure 3 indicates that this fit
overestimates the mean free path at $z\gtrsim 5$ compared to the latest
constraints on $R_{\mathrm{mfp}}$. These data points, suggesting a stronger
redshift dependence, lead us to adopt $A_{\mathrm{mfp}}=210$ cMpc and
$\eta_{\mathrm{mfp,0}}=-9$ for a better match.

It is important to note that with the power-law dependence of
$R_{\mathrm{mfp}}$ provided in Worseck et al. (2014) will approach zero at
high redshifts ($z\gtrsim 7$), potentially hindering the formation of large
ionized bubbles and violating reionization history constraints. We address
this issue by introducing the modification shown in Eq. 5, which has three
additional parameters, $z_{\mathrm{mfp}}$, $\eta_{\mathrm{mfp,1}}$ and
$\eta_{\mathrm{mfp,2}}$. We set $\eta_{\mathrm{mfp,1}}=9$ and
$\eta_{\mathrm{mfp,2}}=1$, exploring two models with $z_{\mathrm{mfp}}=6$ and
$z_{\mathrm{mfp}}=5.5$. In these models, $R_{\mathrm{mfp}}$ converges to 10
and 20 cMpc, respectively, as evident in Figure 3. These models are marked as
‘SinkA’ and ‘SinkB’. We have listed the sink parameter values in Table 1.

###  3.2 Reionization history

Figure 4: Top-left panel: The intergalactic medium (IGM) reionization history,
quantified by the volume-averaged neutral fraction $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$, of our models compared with constraints
from several measurements. These measurements include Ly$\alpha$ dark pixel
fraction (circles; McGreer et al., 2015; Jin et al., 2023), Ly$\alpha$
equivalent width (stars; Mason et al., 2018; Mason et al., 2019; Hoag et al.,
2019; Bolan et al., 2022; Bruton et al., 2023; Morishita et al., 2023; Jones
et al., 2023), Ly$\alpha$ emitter clustering (diamonds; Ouchi et al., 2010),
quasar damping wing (squares; Davies et al., 2018b; Wang et al., 2020; Greig
et al., 2022) and galaxy damping wing (pentagons; Hsiao et al., 2023; Umeda et
al., 2023). Top-right panel: The volume-averaged ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$; solid) and mass-averaged neutral
fraction ($\langle x_{\mathrm{HI}}\rangle_{\mathrm{m}}$; dashed) at the end
stages of reionization ($z\lesssim 7$). Our models show a transition from
inside-out to outside-in nature of reionization when $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 5\times 10^{-5}$. Bottom-left
panel: The plot shows Euler characteristics ($\chi$), which measures the
topology of the distribution neutral islands in the IGM, evolving with
$\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$. All the reionization models
have a distinct evolution, suggesting that the source and sink models affect
the topology. Bottom-right panel: The redshift evolution of the ionizing (UV)
background compared against several constraints from quasar observations
(Wyithe & Bolton, 2011; Calverley et al., 2011; D’Aloisio et al., 2018; Davies
et al., 2018a, 2023). This background grows over time until $z\approx 7$. Post
this epoch, the value is governed by the choice of $R_{\mathrm{mfp}}$, which
can be inferred from the slower growth at low redshifts. The model with
‘SinkB’ has a background more significant than the models with ‘SinkA’ due to
larger $R_{\mathrm{mfp}}$ values at each redshift (see Figure 3).

First, we investigate the reionization of the IGM in the different models
simulated in this work. We show the volume-averaged neutral fraction ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$) of the IGM in the top-left panel of
Figure 4. Our models are carefully calibrated with the choice of
$\zeta_{\mathrm{ion}}$ to to be consistent with a range of observational
constraints on the reionization history. These constraints are represented
with scattered points in the panel. The end of reionization is delayed in our
models and falls well below $z\approx 6$, which is in line with the recent
findings (e.g. Kulkarni et al., 2019; Bosman et al., 2022).

The progress of reionization depends on our choice of source model. In the
case of the ‘Source1’ model ($\alpha_{\star}+\alpha_{\mathrm{esc}}>0$), the
process is primarily driven by large-mass haloes. Conversely, in the ‘Source3’
model ($\alpha_{\star}+\alpha_{\mathrm{esc}}<0$), small-mass haloes play a
dominant role in driving reionization. The reionization history of models
using ‘Source1’ shows a delay compared to those using ‘Source3’ due to the
higher abundance of small-mass haloes in the latter. The model employing
‘Source2’ with $\alpha_{\star}+\alpha_{\mathrm{esc}}=0$ exhibits a
reionization history that falls in between the previous two scenarios.

To understand the impact of the mean free path, we conducted simulations
called ‘Source1_SinkA’ and ‘Source1_SinkB,’ which share the same source model
but differ in the mean free path values. Initially, the reionization histories
of these two models are identical. The $R_{\mathrm{max}}$ value in ‘SinkA’
model is smaller than that in ‘SinkB’ model at all times (see Figure 3).
However, in the middle stages, the reionization process accelerates in
‘Source1_SinkB’ compared to the other model, which can be understood from the
larger value for $R_{\mathrm{max}}$. Therefore, the mean free path for
ionising photons has a substantial impact on the later stages of reionization.

In the top-right panel of Figure 4, we focus on the end stages of reionization
($z\lesssim 6.8$), plotting both volume-averaged ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$) along with the mass-averaged ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{m}}$) neutral fraction. Throughout most of the
reionization process, we observe $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{m}}$$\lesssim$$\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$, indicative of the inside-out nature of
reionization (Iliev et al., 2006). However, as reionization nears completion
($\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 5\times 10^{-5}$), we
observe a transition to an outside-in reionization pattern. During these very
late stages, high recombination rates in dense regions dominate the
reionization process, a phenomenon assumed in several previous studies (e.g.,
Miralda-Escudé et al., 2000; Wyithe & Loeb, 2008).

Additionally, we estimated the Thompson scattering optical depth
($\tau_{\mathrm{e}}$) derived from our models, which serves as a crucial
constraint of the reionization history. We have listed the values in Table 1.
The Planck mission has precisely constrained this parameter to $0.054\pm
0.007$ at the 68 per cent confidence level (Planck Collaboration et al.,
2020), representing a reduction compared to previous WMAP results (Hinshaw et
al., 2013). This decline suggests a potential delay in the end of reionization
(e.g. Mitra et al., 2015). Our models agree with these latest
$\tau_{\mathrm{e}}$ constraints.

###  3.3 Topology of neutral islands

To study the topology of the neutral island distribution in the IGM, we
estimated the Euler characteristics ($\chi$) of our simulations, which
describes the evolution of the connectivity of these islands due to the
reionization process. We refer interested readers to Friedrich et al. (2011)
and Giri & Mellema (2021) for more discussion about topological evolution
during reionization. The Euler characteristics is defined as

$\displaystyle\chi=N_{\rm bubbles}-N_{\rm tunnels}+N_{\rm islands}\ ,$ (10)

where $N_{\rm bubbles}$, $N_{\rm tunnels}$ and $N_{\rm islands}$ are the
number of ionized bubbles, tunnels and neutral islands, respectively.

In the bottom-left panel of Figure 4 , we present $\chi$ for all our models at
different epochs of reionization. At early times, $\chi$ increases due to the
increasing number of ionized bubbles in the IGM. When small mass haloes drive
reionization, a larger number of these ionized bubbles form. Therefore,
‘Source3’ has the most prominent peak value of $\chi$. The sink models do not
affect the topology at these early stages, as revealed by the overlapping
evolution of the two sink models with the same source property. Thus, we infer
that the $\chi$ at early stages of reionization is useful in distinguishing
between source models.

Over time, the ionized bubbles merge, creating tunnels that decrease $\chi$
and lead to a negative minimum during the intermediate stages. We observe
differences in the behaviour of the two sink models during these stages. In
the case of the ‘SinkB’ model, which allows larger ionized bubbles to form for
each source, this merger happens more quickly compared to the ‘SinkA’ models.
When connections occur with fewer number of bubbles, we get less number of
tunnels, which is the case for ‘SinkB’. Consequently, the ‘Source1_SinkB’
model exhibits a minimum with the highest value. This distinctive evolution of
$\chi$ can help distinguish between different sink models.

In the late stages of reionization, $\chi$ depended on the number of neutral
islands. Due to the small number of these islands, the magnitude of $\chi$
remains low during this period. We observed that the topology of the neutral
islands tends to converge across all models at late times. Although the
numerical count of these neutral islands is relatively small, their sizes
remain notably significant, which will be further explored in Section 4.

###  3.4 Ionizing photon background

We now study the mean ionization rate ($\Gamma$) within the IGM provided by
our simulations. The bottom-right panel in Figure 4 shows our results compared
to observational constraints. Notably, these constraints predominantly apply
to lower redshifts ($z\lesssim 6.5$). Replicating these constraints has posed
a significant challenge in the field, often requiring abrupt changes in source
properties towards the end stages of reionization (Chardin et al., 2016;
Puchwein et al., 2019). Our simulations show a smooth evolution of $\Gamma$,
produced with our model for unresolved sinks through a smooth evolution of
$R_{\mathrm{mfp}}$.

For a particular source model, the growth of $\Gamma$ is sensitive to the
expansion of ionized bubbles merging together, rendering the IGM transparent
to UV photons. During the final stages of reionization, most parts of the IGM
become transparent to UV photons emitted by sources throughout the simulation
volume. Therefore, we can consider the evolution of $\Gamma\propto
R_{\mathrm{mfp}}$ (e.g., Haardt & Madau, 2012) during this phase. In previous
models (e.g., Dixon et al., 2016; Giri et al., 2019), we assumed a non-
evolving $R_{\mathrm{mfp}}$ that resulted in the flattening of $\Gamma$ during
late times. In this study, however, we assumed an increasing
$R_{\mathrm{mfp}}$ during late times (see Section 2.2.2), causing $\Gamma$ to
increase as reionization proceeds during the final stages. Furthermore, the
evolution of $\Gamma$ observed in ‘SinkA’ models closely resembles that seen
in several hydro-dynamical simulations, such as Thesan (Garaldi et al., 2022)
and CODA-III (Lewis et al., 2022).

##  4 The 21-cm signal

This section presents the evolution of the 21-cm signal in our simulations. In
subection 4.1, we begin by describing the quantity measured by radio
telescopes and study the large-scale distribution of H  I during the end
stages of reionization. In subection 4.2, we discuss the existence of large
neutral islands still remaining during this stage. Lastly, we present the
power spectrum of the 21-cm signal in the concluding subsection.

Figure 5: Slices from our model at redshift $z\approx 5.436$. Top panels: We
show the snapshot from our $N$-body run in volumes with a box length of 298
cMpc. The left and right panels depict the dark matter over-density ($\delta$)
and cumulative stellar mass on the grid for our source model that was fixed
using the UVLFs (see Section 3.1.1) on the same grid, respectively. Bottom
panels: We illustrate the brightness temperature ($\delta T_{\mathrm{b}}$)
corresponding to the 21-cm signal measured by radio telescopes. The left and
right panels show this signal produced by the H

I

in the intergalactic medium (IGM) of our Source1_SinkA model and shielded H

I

inside galaxies residing in the dark matter haloes. In all the panels, we zoom
on to a region with a visibly large neutral island. Comparing these zoomed
regions in the bottom panels, we infer that this neutral island is situated in
a region with fewer haloes or the cosmic void.

###  4.1 Differential brightness temperature

The 21-cm signal produced by H  I during reionization can be found at radio
frequencies below 235 MHz. The interferometry-based radio telescopes record
the differential brightness temperature corresponding to this signal, which is
given as (e.g. Mellema et al., 2013)

$\displaystyle\delta T_{b}(\mathbf{x},z)=T_{0}(z)x_{\rm
HI}(\mathbf{x},z)\left[1+\delta_{b}(\mathbf{x},z)\right]\left(1-\frac{T_{\rm
cmb}(z)}{T_{\rm s}(\mathbf{x},z)}\right),$ (11)

where $T_{\rm cmb}$ is the cosmic microwave background (CMB) temperature, and
$x_{\rm HI}$ is the neutral hydrogen fraction. $\delta_{b}$ describes the H  I
gas perturbation field that is assumed to follow the dark matter. The factor
$T_{0}$ depends only on cosmology that is given by,

$\displaystyle T_{0}(z)\approx 22\left(\frac{\Omega_{\rm
b}h^{2}}{0.022}\right)\left(\frac{0.144}{\Omega_{\rm
m}h^{2}}\frac{1+z}{7}\right)^{1/2}{\rm mK}\ .$ (12)

$T_{\mathrm{s}}$ is the spin temperature of the gas, which we assume to be
much higher than $T_{\mathrm{CMB}}$. Studies have shown that this spin
temperature saturation occurs by the early stages of reionization (e.g. Ghara
et al., 2015; Ross et al., 2017, 2021). Therefore it remains a reasonable
assumption during the later stages of reionization and the post-reionization
era.

Figure 6: The size distribution of the neutral islands in ‘Source1_SinkA’
(blue) and ‘Source1_SinkB’ (red) models at three different volume-averaged
neutral fractions $\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$. These
distributions show that significantly large neutral islands remain until the
very late stages of reionization. The distributions from ‘Source1_SinkB’ have
the peaks positioned at larger sizes than ‘Source1_SinkA’ during each epoch.

We first inspect the large-scale distribution of H  I in our simulation suite.
Figure 5 illustrates slices from different fields at $z\approx 5.436$. The
top-left panel shows the dark matter distribution. We assume the baryons to
follow this distribution in IGM. In the top-right panel, we show the
distribution of stellar masses. These gridded distributions are estimated by
assigning stellar mass to dark matter haloes and summing them up for each grid
point. The stellar masses were calculated using Eq. 1 with the free parameters
fixed against the UVLF observations (see Section 3.1.1). We observe that the
distribution of the source masses follows the large-scale distribution of
matter, indicating that our sources reside in the dense cosmic filaments.

In the bottom panels of Figure 5, we present the $\delta T_{\mathrm{b}}$
slices from our ‘Source1_SinkA’ model. At $z\approx 5.436$, the volume-
averaged neutral fraction is $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-4}$, and so it corresponds
to the final stages of reionization. We chose this epoch to understand the
phase when reionization transitions from inside-out to outside-in, which was
identified in Section 3.2. The bottom-left panel displays the $\delta
T^{\mathrm{IGM}}_{\mathrm{b}}$, which corresponds to the signal from H  I in
the IGM. We observe that the patterns in this signal follow the dense
structures in the matter distribution (top-left panel). This correlation is
caused due to the higher recombination rates in the dense regions.

While most regions exhibit a very low signal strength ($\delta
T^{\mathrm{IGM}}_{\mathrm{b}}\lesssim 10^{-2}\,\rm mK$), there are a few areas
with high values, $\gtrsim 10\,\rm mK$. These regions correspond to neutral
islands that remain shielded from UV photons. These islands are located at
large enough distance from photon sources such that the recombination rate,
although low, is enough to counterbalance the incoming UV photons. In the
inset, we show a zoom in on one such a neutral island. This inset correspond
to a size of $30\,\rm cMpc$ per side. Our simulations contain numerous such
islands in other slices, and we will conduct a statistical investigation of
these islands in the next subsection.

In the bottom-right panel, we show the $\delta T^{\mathrm{halo}}_{\mathrm{b}}$
produced by H  I inside galaxies residing in haloes, which we assigned by
using the method described in Section 2.3. This signal follows the
distribution of sources (top-right panel) and the dense cosmic filaments (top-
left panel) by design. At this very late stage, the strength of this signal is
comparable to the signal produced in the IGM.Within the neutral island region
shown in the zoom in, the $\delta T^{\mathrm{halo}}_{\mathrm{b}}$ is
predominantly low ($\lesssim 10^{-2}\,\rm mK$), where the neutral island is
located. Therefore, the resulting 21-cm signal ($\delta
T^{\mathrm{halo}}_{\mathrm{b}}+\delta T^{\mathrm{IGM}}_{\mathrm{b}}$) from
this region will follow the under-dense matter distribution. We will discuss
the impact of these islands in the statistical measurements of the 21-cm
signal in subsection 4.3.

###  4.2 Neutral islands during the very late stages of reionization

In Giri et al. (2019), we found few but very large (with lengths $\sim 100$
cMpc) neutral islands at $\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx
10^{-2}$. With the models in this study, which extend to even later epochs
($z\lesssim 6$), we continue to observe notably large neutral islands, as
discussed earlier. While Figure 5 visually presents one such large neutral
island in the slice, the entire simulation volume contains more instances. We
employed the mean-free-path size distribution algorithm implemented in
Tools21cm (Giri et al., 2020a) to investigate these islands further. See Giri
et al. (2018a) for more information about this algorithm.

In Figure 6, we present the island size distribution (ISD), illustrating the
probability of sizes for neutral islands in two sink models with the same
sources (‘Source1’) at different epochs. The ‘SinkA’ model is represented in
blue, while the ‘SinkB’ model is shown in red. We see that significantly large
neutral islands exist in our model that can impact the statistical
measurements. For the case at $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-4}$ (dot-dashed line),
‘SinkA’ and ‘SinkB’ have an average value of island sizes $\bar{R}\simeq
3.2\,\rm cMpc$ and $3.6\,\rm cMpc$, respectively. Additionally, we observe
neutral islands as large as $\sim$20 cMpc in both models. These large islands
can be detected in the upcoming image data from SKA-Low using the structure
identification approaches developed in Giri et al. (2018b) and Bianco et al.
(2021b); Bianco et al. (2024). At $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-3}$, the neutral island
sizes can reach $\sim$40 cMpc. We will discuss their impact on the power
spectra in the next subsection.

The two sink models show distinct ISDs at all the epochs in Figure 6. The
island sizes in the ‘SinkB’ model is larger than that of ‘SinkA’ at all
epochs. In Giri et al. (2019), we developed a model based on packing spheres
in simulation volumes to comprehend the distribution of neutral islands. This
model suggested that the empty spaces between spheres of large sizes are
large. Since the ‘SinkB’ model permits the formation of larger bubbles around
every source compared to ‘SinkA’, the resulting neutral islands are more
prominent, as revealed in the figure.

###  4.3 21-cm power spectrum

Figure 7: The 21-cm power spectra
$\Delta_{\mathrm{21}}=k^{3}P_{\mathrm{21}}/(2\pi^{2})$ for our models
(Source1_SinkA: blue, Source2_SinkA: orange, Source3_SinkA: green and
Source1_SinkB: red). The solid lines represent $\Delta_{\mathrm{21}}$,
assuming the contribution from H

I

in both the IGM and inside the haloes. The dashed and dotted lines represent
$\Delta_{\mathrm{21}}$ considering only H

I

in the IGM and inside haloes, respectively. Top panel: We show the redshift
evolution of $\Delta_{\mathrm{21}}$ at $k=0.1$ (left), 0.5 (middle) and 1.0
(right) Mpc-1. To focus on the end stages of reionization, the insets of the
panels zoom into the region where $z\lesssim 6.5$. Until very late stages
($\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\lesssim 10^{-4}$) of
reionization, the contribution from the H

I

in the IGM dominates. Bottom panel: We compare the contribution of the H

I

in the IGM and inside haloes to $\Delta_{\mathrm{21}}$ at three very late
stages where $\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$ $\approx 10^{-3}$
(left), $10^{-4}$ (middle) and $10^{-5}$ (right). The contribution from the
latter initially dominates on small scales and evolves to all scales over
time.

The spatial characteristics of the 21-cm signal can be probed with the power
spectrum $(P_{\mathrm{21}}(k)=\langle\delta T^{*}_{\mathrm{b}}(k)\delta
T_{\mathrm{b}}(k)\rangle)$ that the radio telescopes are attempting to detect.
Post-reionization H  I follows the galaxy distribution and, therefore, the
matter distribution. Hence, most studies model the post-reionization 21-cm
power spectrum as a biased version of the matter power spectrum (e.g. Santos
et al., 2015; Padmanabhan et al., 2017; Obuljen et al., 2018). However, we do
not make this assumption and estimate the power spectrum directly from our
simulation volumes where the H  I both in the IGM and inside galaxies are
included.

Figure 7 illustrates the dimensionless power spectrum
$\Delta_{\mathrm{21}}(k)=k^{3}P_{\mathrm{21}}(k)/(2\pi^{2})$ of our models.
The top panels show the redshift evolution at three different scales
($k\approx 0.1,0.5,1.0~{}\mathrm{Mpc}^{-1}$ ). While the solid lines represent
the $\Delta_{\mathrm{21}}$ of the signal produced by H  I from both haloes and
IGM (total), the dashed lines represent the signal from the IGM alone. In
these panels, we do not include the $\Delta_{\mathrm{21}}$ for the H  I inside
galaxies because their magnitude is small for most redshifts. However, the
noticeable difference between the solid and dashed line is caused by the
inclusion of this H  I .

In all our models, this power spectrum is dominated by the H  I in the IGM
until very late stages of reionization ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-3}$). After this epoch, the
contribution from the H  I inside galaxies begin to become more critical. The
transition of the 21-cm signal to probing the galaxies occurs when the
$\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$ is less than $10^{-4}$. In the
top panels, this epoch is marked by a flattening of the power spectra in all
models and at all scales. By detecting this unique signature, future radio
experiments would be able to provide evidence of the end of reionization. The
epoch of the transition from IGM to galaxies coincides with the time when
$\langle x_{\mathrm{HI}}\rangle_{\mathrm{m}}$ becomes larger than $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$. As discussed in Section 3.2, the
distribution of H  I in the IGM transitions from an inside-out to outside-in
nature in our simulations.

We selected three epochs from our simulations to focus on the transition
process and plot them in the bottom panels of Figure 7. Along with the
$\Delta_{\mathrm{21}}$ of the total (solid lines) and IGM (dashed lines)
fields, we include the field that corresponds to the H  I in the haloes only
(dotted lines). The left panel displays the $\Delta_{\mathrm{21}}$ when
$\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$ is approximately $10^{-3}$,
where the power from the IGM is seen to be $\sim 2$ orders of magnitude higher
than that from the galaxies. Therefore we can neglect still neglect the latter
contribution at this time.

As time progresses, the small-scale signal transitions to probe the H  I
inside haloes, as evident in the middle panel, when $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-4}$. This scale dependent
transition is consistent with the findings in Xu et al. (2019). As discussed
in Section 4.2, a significant number of large neutral islands remains during
this epoch. We observe that all models have an IGM $\Delta_{\mathrm{21}}$ with
a knee-like feature at wave-number $k_{\mathrm{knee}}$ between 0.4 and 0.6
$\mathrm{Mpc}^{-1}$. This scale approximately corresponds to the average scale
of the neutral islands during this time as $k_{\mathrm{knee}}\sim
2/\bar{R}$333In Georgiev et al. (2022), a similar relation was identified, but
with the sizes of ionized bubbles. This $k$ value determined the scale beyond
which the 21-cm power spectrum became a biased tracer of matter distribution
during the early stages of reionization.. Due to this feature, we cannot
assume $\Delta_{\mathrm{21}}$ to be a biased version of the matter power
spectrum.

In the right-hand panel, we show $\Delta_{\mathrm{21}}$ when $\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-5}$, which is the time by
when the contribution from the IGM has become negligible. Almost all the
signal comes from the H  I inside dark matter haloes. As the
$\Delta_{\mathrm{21}}$ is probing the haloes during this time, it can be
modelled as a biased version of the matter power spectrum. Therefore,
$\Delta_{\mathrm{21}}$ can be used to study cosmological models after this
transition era.

##  5 Conclusion

In this study, we enhanced our fully numerical cosmic reionization simulation
framework, pyC2Ray (Hirling, Bianco et al. 2023), enabling it to incorporate
various high-redshift observations, including the ultraviolet luminosity
functions (UVLFs) and Lyman-$\alpha$ (Ly$\alpha$) absorption spectra.
Utilizing this framework, we investigated the evolution of the cosmological
21-cm signal, particularly emphasizing the very late stages of reionization.
Our simulation suite was constructed based on early galaxy formation models
aligned with the latest high-redshift measurements. We calibrated our model
parameters for star formation within dark matter haloes using available UVLF
measurements at high redshifts. Additionally, we considered models for UV
photon escape that are consistent with constraints for the Ly$\alpha$ escape
fraction.

In simulation volumes such as those utilized here (with a box length of 298
cMpc), the small scale absorbers are not easily resolved. Therefore, we
implemented a sink model that limits the propagation of UV photons beyond a
certain distance $R_{\mathrm{mfp}}$. Previous studies have demonstrated the
effectiveness of this sink model in explaining observations of Ly$\alpha$
absorption spectra (e.g., Cain et al., 2021; Bosman et al., 2022; Davies et
al., 2023). A comparison of several methods for modelling these sinks will be
presented in an upcoming paper (Georgiev et al., in prep). We adopted two
redshift dependencies for $R_{\mathrm{mfp}}$, calibrated to constraints from
relevant studies. These models yield a distinct evolution of reionization
during the later stages ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\lesssim 0.5$). Most frameworks for
reionization modelling typically assume a fixed $R_{\mathrm{mfp}}$ when
interpreting observations (e.g. Mondal et al., 2020; Greig et al., 2021b; Qin
et al., 2021b), which may introduce biases on constraints when analyzing the
final stages of reionization.

We also studied the evolution of the UV background ($\Gamma$) in our
simulations, which is an output of the radiative transfer solver in pyC2Ray.
This background is modulated by the sink model and, consequently, influenced
by the choice of $R_{\mathrm{mfp}}$, consistent with findings in previous
studies (e.g., Haardt & Madau, 2012; Sobacchi & Mesinger, 2014; Becker et al.,
2021; Gaikwad et al., 2023). It is important to note that although we enforce
an $R_{\mathrm{mfp}}$ for the radiative transfer around each source, the sizes
of ionized regions in our simulations can be much larger due to overlaps. Our
models consistently align with constraints on $R_{\mathrm{mfp}}$ and $\Gamma$.
This improvement is a significant enhancement compared to our previous 21-cm
signal simulations (e.g. Dixon et al., 2016; Giri et al., 2019).

A caveat of our study of the UV background is that we have considered only one
method of modelling the sinks. For example, a clumping factor evolving with
redshift or density-dependent clumping factor (e.g. Mao et al., 2020; Bianco
et al., 2021a) can affect the IGM reionization during the end stages and,
consequently, affect the evolution of the $R_{\mathrm{mfp}}$ and $\Gamma$.
Previous studies have also attributed the measured evolution of the UV
background to the uniqueness of the UV emissivity during late stages (e.g.
Chardin et al., 2016; Puchwein et al., 2019). We defer the exploration of such
sink and source models to future work.

The 21-cm signal serves as a valuable probe for understanding the topology of
the large-scale distribution of H  I in our Universe (Giri & Mellema, 2021;
Kapahtia et al., 2021; Elbers & van de Weygaert, 2023). In the post-
reionization era, this signal primarily probes the H  I content within dark
matter haloes, becoming a biased tracer of the matter power spectrum. During
this period, the topology of H  I distribution is typically assumed to be
outside-in in nature (e.g., Miralda-Escudé et al., 2000; Wyithe & Loeb, 2008).
In this topology, dense regions exhibit a higher neutral fraction than under-
dense regions due to increased recombination rates (Finlator et al., 2009;
Choudhury et al., 2009). We investigated this transition phase and found the
topology to remain inside-out until the Universe reaches a neutral fraction of
$\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-4}$, a
significantly later stage compared to previous models (e.g., Finlator et al.,
2009).

Our reionization models exhibit a distinct evolution of the 21-cm power
spectra, highlighting the potential of this signal to constrain both source
and sink models. Focusing on the very late stages of reionization ($\langle
x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\lesssim 10^{-3}$), we studied the
dependence of the signal to the H  I content in both the IGM and dark matter
haloes. The relevance of H  I inside dark matter haloes becomes more
pronounced after $\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx
10^{-4}$, coinciding with the phase when the topology of reionization
undergoes a transition in our models. The redshift evolution of the 21-cm
power spectra reveals a distinct flattening at this epoch, potentially serving
as a marker for the 21-cm signal entering the post-reionization era and the H
I topology transitioning to outside-in. Furthermore, we found that the
transition process of the signal from probing the IGM to the haloes is a scale
dependent, which begins at small-scale and transfer to large-scales. This
finding is consistent with the study in Xu et al. (2019).

In Giri et al. (2019), we investigated the existence of large neutral islands
during the end stages of reionization, detectable in the image data from the
SKA-Low. In this study, our simulations reveal a statistically significant
number of such neutral islands during the very late stages of reionization
($\langle x_{\mathrm{HI}}\rangle_{\mathrm{v}}$$\approx 10^{-4}$). The sizes of
these islands were as large as 20 cMpc, which can be detected in the image
data with our feature identification framework (Giri et al., 2018b; Bianco et
al., 2021b; Bianco et al., 2024). We studied the distribution of the sizes of
these islands, demonstrating that they are significant enough to impact the
large-scale fluctuations of the 21-cm signal. Notably, we observe a knee-like
feature in the 21-cm power spectra corresponding to these islands, persisting
until the period when the signal of H  I in haloes starts to dominate.

In modelling the H  I inside dark matter haloes, we do not consider the in-
and out-flow of gas between the IGM and galaxies, as this requires
sophisticated hydro-dynamical simulation. Additionally, we follow a
straightforward approach using the $M_{\mathrm{HI}}-M_{\mathrm{h}}$ relation
from Padmanabhan et al. (2017) that assumes that the H  I in massive galaxies,
residing in high mass halo $M_{\mathrm{h}}\gtrsim 10^{10}\,M_{\odot}$, can
self-shield against the UV photons. However, recent H  I surveys provided
observational evidence that suggests a more complex scenario with Ultra Faint
Dwarf Galaxy able to keep most of their reservoir of H  I throughout EoR (see
Irwin et al., 2007; Saul et al., 2012; Giovanelli et al., 2013; Janesh et al.,
2019). These low-mass galaxies can impact our study of the end stages of
reionization, and we will explore their impact in future work.

Upcoming observations of the 21-cm signal will help study these late stages of
reionization ($z\lesssim 6$) in more detail. As the low redshifts are expected
to have a greater signal-to-noise ratio, we expect better constraining
capability during these stages. Therefore, our models can be helpful in
interpreting the signal during these late epochs.

## Acknowledgements

We acknowledge Benoit Semelin and Yves Revaz for their helpful discussions.
Nordita is supported in part by NordForsk. GM is supported by the Swedish
Research Council project grant 2020-04691_VR and AS is supported by the Swiss
National Science Foundation (SNF) via the grant PCEFP2_181157. We acknowledge
the allocation of computing resources provided by the National Academic
Infrastructure for Supercomputing in Sweden (NAISS) and the Swiss National
Supercomputing Centre (CSCS) under the SKA share with the project ID sk015. We
have utilised the following Python packages for manipulating the simulation
outputs and plotting results: numpy (Van Der Walt et al., 2011), scipy
(Virtanen et al., 2020), and matplotlib (Hunter, 2007).

## Data Availability

The source code used for the simulations of this study, the Pkdgrav3
(https://bitbucket.org/dpotter/pkdgrav3/) and pyC2Ray
(https://github.com/cosmic-reionization/pyC2Ray), are publicly available. We
have made the 21-cm signal simulation data freely available at
https://doi.org/10.5281/zenodo.10785609. Our data can be read and manipulated
using our public tool, Tools21cm (https://github.com/sambit-giri/tools21cm).

## References

  * Abdalla et al. (2022) Abdalla E., et al., 2022, Astronomy & Astrophysics, 664, A14
  * Atek et al. (2018) Atek H., Richard J., Kneib J.-P., Schaerer D., 2018, Monthly Notices of the Royal Astronomical Society, 479, 5184
  * Atek et al. (2023) Atek H., et al., 2023, Monthly Notices of the Royal Astronomical Society, 519, 1201
  * Bagla et al. (2010) Bagla J., Khandai N., Datta K. K., 2010, Monthly Notices of the Royal Astronomical Society, 407, 567
  * Bakx et al. (2023) Bakx T. J., et al., 2023, Monthly Notices of the Royal Astronomical Society, 519, 5076
  * Barnes & Haehnelt (2014) Barnes L. A., Haehnelt M. G., 2014, Monthly Notices of the Royal Astronomical Society, 440, 2313–2321
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, Monthly Notices of the Royal Astronomical Society, 508, 1853
  * Begley et al. (2022) Begley R., et al., 2022, Monthly Notices of the Royal Astronomical Society, 513, 3510
  * Begley et al. (2024) Begley R., et al., 2024, Monthly Notices of the Royal Astronomical Society, 527, 4040
  * Behroozi et al. (2020) Behroozi P., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5702
  * Bertschinger (1998) Bertschinger E., 1998, Annual Review of Astronomy and Astrophysics, 36, 599
  * Bharadwaj & Srikant (2004) Bharadwaj S., Srikant P. S., 2004, Journal of Astrophysics and Astronomy, 25, 67–80
  * Bianco et al. (2021a) Bianco M., Iliev I. T., Ahn K., Giri S. K., Mao Y., Park H., Shapiro P. R., 2021a, Monthly Notices of the Royal Astronomical Society, 504, 2443
  * Bianco et al. (2021b) Bianco M., Giri S. K., Iliev I. T., Mellema G., 2021b, Monthly Notices of the Royal Astronomical Society, 505, 3982
  * Bianco et al. (2024) Bianco M., Giri S. K., Prelogović D., Chen T., Mertens F. G., Tolley E., Mesinger A., Kneib J.-P., 2024, Monthly Notices of the Royal Astronomical Society, 528, 5212
  * Bigot-Sazy et al. (2015) Bigot-Sazy M.-A., et al., 2015, arXiv preprint arXiv:1511.03006
  * Bird et al. (2014) Bird S., Vogelsberger M., Haehnelt M., Sijacki D., Genel S., Torrey P., Springel V., Hernquist L., 2014, Monthly Notices of the Royal Astronomical Society, 445, 2313–2324
  * Blas et al. (2011) Blas D., Lesgourgues J., Tram T., 2011, Journal of Cosmology and Astroparticle Physics, 2011, 034
  * Bolan et al. (2022) Bolan P., et al., 2022, Monthly Notices of the Royal Astronomical Society, 517, 3263
  * Bosman et al. (2022) Bosman S. E., et al., 2022, Monthly Notices of the Royal Astronomical Society, 514, 55
  * Bouwens et al. (2021) Bouwens R., et al., 2021, The Astronomical Journal, 162, 47
  * Bowler et al. (2020) Bowler R., Jarvis M., Dunlop J., McLure R., McLeod D., Adams N., Milvang-Jensen B., McCracken H., 2020, Monthly Notices of the Royal Astronomical Society, 493, 2059
  * Bruton et al. (2023) Bruton S., Lin Y.-H., Scarlata C., Hayes M. J., 2023, The Astrophysical Journal Letters, 949, L40
  * Bunker et al. (2023) Bunker A. J., et al., 2023, arXiv preprint arXiv:2302.07256
  * Cain et al. (2021) Cain C., D’Aloisio A., Gangolli N., Becker G. D., 2021, The Astrophysical Journal Letters, 917, L37
  * Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, Monthly Notices of the Royal Astronomical Society, 412, 2543
  * Camera et al. (2013) Camera S., Santos M. G., Ferreira P. G., Ferramacho L., 2013, Physical Review Letters, 111, 171302
  * Carucci et al. (2017) Carucci I. P., Corasaniti P.-S., Viel M., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 018
  * Chang et al. (2010) Chang T.-C., Pen U.-L., Bandura K., Peterson J. B., 2010, Nature, 466, 463
  * Chardin et al. (2016) Chardin J., Puchwein E., Haehnelt M. G., 2016, Monthly Notices of the Royal Astronomical Society, 465, 3429
  * Chisholm et al. (2018) Chisholm J., et al., 2018, Astronomy & Astrophysics, 616, A30
  * Choudhury & Paranjape (2018) Choudhury T. R., Paranjape A., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3821
  * Choudhury et al. (2009) Choudhury T. R., Haehnelt M. G., Regan J., 2009, Monthly Notices of the Royal Astronomical Society, 394, 960
  * Choudhury et al. (2021) Choudhury T. R., Paranjape A., Bosman S. E., 2021, Monthly Notices of the Royal Astronomical Society, 501, 5782
  * Ciardi et al. (2006) Ciardi B., Scannapieco E., Stoehr F., Ferrara A., Iliev I., Shapiro P., 2006, Monthly Notices of the Royal Astronomical Society, 366, 689
  * D’Aloisio et al. (2018) D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018, Monthly Notices of the Royal Astronomical Society, 473, 560
  * Davies et al. (2018a) Davies F. B., Hennawi J. F., Eilers A.-C., Lukić Z., 2018a, The Astrophysical Journal, 855, 106
  * Davies et al. (2018b) Davies F. B., et al., 2018b, The Astrophysical Journal, 864, 142
  * Davies et al. (2023) Davies F. B., et al., 2023, arXiv preprint arXiv:2312.08464
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Physics Reports, 780, 1
  * Dayal & Giri (2024) Dayal P., Giri S. K., 2024, Monthly Notices of the Royal Astronomical Society, 528, 2784
  * Dekel et al. (2013) Dekel A., Zolotov A., Tweed D., Cacciato M., Ceverino D., Primack J., 2013, Monthly Notices of the Royal Astronomical Society, 435, 999
  * Dixon et al. (2016) Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, Monthly Notices of the Royal Astronomical Society, 456, 3011
  * Donnan et al. (2023) Donnan C., et al., 2023, Monthly Notices of the Royal Astronomical Society, 518, 6011
  * Elbers & van de Weygaert (2023) Elbers W., van de Weygaert R., 2023, Monthly Notices of the Royal Astronomical Society, 520, 2709
  * Fan et al. (2006) Fan X., Carilli C., Keating B., 2006, Annu. Rev. Astron. Astrophys., 44, 415
  * Finlator et al. (2009) Finlator K., Özel F., Davé R., Oppenheimer B. D., 2009, Monthly Notices of the Royal Astronomical Society, 400, 1049
  * Friedrich et al. (2011) Friedrich M. M., Mellema G., Alvarez M. A., Shapiro P. R., Iliev I. T., 2011, Monthly Notices of the Royal Astronomical Society, 413, 1353
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, MNRAS, 525, 4093
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, Monthly Notices of the Royal Astronomical Society, 512, 4909
  * Georgiev et al. (2022) Georgiev I., et al., 2022, MNRAS, 513, 5109
  * Ghara et al. (2015) Ghara R., Choudhury T. R., Datta K. K., 2015, Monthly Notices of the Royal Astronomical Society, 447, 1806
  * Ghara et al. (2020) Ghara R., et al., 2020, MNRAS, 493, 4728
  * Ghara et al. (2021) Ghara R., et al., 2021, MNRAS, 503, 4551
  * Gillet et al. (2020) Gillet N. J., Mesinger A., Park J., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1980
  * Giovanelli et al. (2013) Giovanelli R., et al., 2013, The Astronomical Journal, 146, 15
  * Giri & Mellema (2021) Giri S. K., Mellema G., 2021, MNRAS, 505, 1863
  * Giri & Schneider (2022) Giri S. K., Schneider A., 2022, Phys. Rev. D, 105, 083011
  * Giri et al. (2018a) Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018a, Monthly Notices of the Royal Astronomical Society, 473, 2949
  * Giri et al. (2018b) Giri S. K., Mellema G., Ghara R., 2018b, Monthly Notices of the Royal Astronomical Society, 479, 5596
  * Giri et al. (2019) Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T., 2019, Monthly Notices of the Royal Astronomical Society, 489, 1590
  * Giri et al. (2020a) Giri S., et al., 2020a, JOSS, 5, 2363
  * Giri et al. (2020b) Giri S. K., Zackrisson E., Binggeli C., Pelckmans K., Cubo R., 2020b, Monthly Notices of the Royal Astronomical Society, 491, 5277
  * Giri et al. (2023) Giri S. K., Schneider A., Maion F., Angulo R. E., 2023, Astronomy & Astrophysics, 669, A6
  * Greig & Mesinger (2015) Greig B., Mesinger A., 2015, MNRAS, 449, 4246
  * Greig et al. (2021a) Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021a, Monthly Notices of the Royal Astronomical Society, 500, 5322
  * Greig et al. (2021b) Greig B., et al., 2021b, Monthly Notices of the Royal Astronomical Society, 501, 1
  * Greig et al. (2022) Greig B., Mesinger A., Davies F. B., Wang F., Yang J., Hennawi J. F., 2022, Monthly Notices of the Royal Astronomical Society, 512, 5390
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, The Astrophysical Journal, 746, 125
  * Hassan et al. (2023) Hassan S., et al., 2023, ApJ, 958, L3
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, The Astrophysical Journal Supplement Series, 208, 19
  * Hirling et al. (2023) Hirling P., Bianco M., Giri S. K., Iliev I. T., Mellema G., Kneib J.-P., 2023, arXiv preprint arXiv:2311.01492
  * Hoag et al. (2019) Hoag A., et al., 2019, The Astrophysical Journal, 878, 12
  * Hsiao et al. (2023) Hsiao T. Y.-Y., et al., 2023, arXiv preprint arXiv:2305.03042
  * Hunter (2007) Hunter J. D., 2007, Comput. Sci. Eng., 9, 90
  * Iliev et al. (2006) Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, Monthly Notices of the Royal Astronomical Society, 369, 1625
  * Iliev et al. (2012) Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., Ahn K., 2012, Monthly Notices of the Royal Astronomical Society, 423, 2222
  * Iliev et al. (2014) Iliev I. T., et al., 2014, MNRAS, 439, 725
  * Irwin et al. (2007) Irwin M. J., et al., 2007, The Astrophysical Journal, 656, L13
  * Ishigaki et al. (2018) Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, The Astrophysical Journal, 854, 73
  * Janesh et al. (2019) Janesh W., Rhode K. L., Salzer J. J., Janowiecki S., Adams E. A. K., Haynes M. P., Giovanelli R., Cannon J. M., 2019, The Astronomical Journal, 157, 183
  * Jaskot et al. (2019) Jaskot A. E., Dowd T., Oey M., Scarlata C., McKinney J., 2019, The Astrophysical Journal, 885, 96
  * Jensen et al. (2016) Jensen H., Zackrisson E., Pelckmans K., Binggeli C., Ausmees K., Lundholm U., 2016, The Astrophysical Journal, 827, 5
  * Jin et al. (2023) Jin X., et al., 2023, The Astrophysical Journal, 942, 59
  * Jones et al. (2023) Jones G. C., et al., 2023, arXiv preprint arXiv:2306.02471
  * Kapahtia et al. (2021) Kapahtia A., Chingangbam P., Ghara R., Appleby S., Choudhury T. R., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 026
  * Kern et al. (2017) Kern N. S., Liu A., Parsons A. R., Mesinger A., Greig B., 2017, The Astrophysical Journal, 848, 23
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in AASKA14. http://pos.sissa.it/
  * Kulkarni et al. (2019) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E., Puchwein E., Chardin J., Aubert D., 2019, Monthly Notices of the Royal Astronomical Society: Letters, 485, L24
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv preprint arXiv:1104.2932
  * Lewis et al. (2022) Lewis J. S., et al., 2022, Monthly Notices of the Royal Astronomical Society, 516, 3389
  * Lin et al. (2024) Lin X., et al., 2024, arXiv preprint arXiv:2401.09532
  * Liu & Slatyer (2018) Liu H., Slatyer T. R., 2018, Physical Review D, 98, 023501
  * Livermore et al. (2017) Livermore R., Finkelstein S., Lotz J., 2017, The Astrophysical Journal, 835, 113
  * Lovell (2024) Lovell M. R., 2024, Monthly Notices of the Royal Astronomical Society, 527, 3029
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Maller & Bullock (2004) Maller A. H., Bullock J. S., 2004, MNRAS, 355, 694
  * Mao et al. (2020) Mao Y., Koda J., Shapiro P. R., Iliev I. T., Mellema G., Park H., Ahn K., Bianco M., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1600
  * Mason et al. (2018) Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., De Barros S., Vanzella E., 2018, The Astrophysical Journal, 856, 2
  * Mason et al. (2019) Mason C. A., et al., 2019, Monthly Notices of the Royal Astronomical Society, 485, 3947
  * McBride et al. (2009) McBride J., Fakhouri O., Ma C.-P., 2009, Monthly Notices of the Royal Astronomical Society, 398, 1858
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, Monthly Notices of the Royal Astronomical Society, 447, 499
  * McLure et al. (2013) McLure R., et al., 2013, Monthly Notices of the Royal Astronomical Society, 432, 2696
  * Mellema et al. (2006) Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, New Astronomy, 11, 374
  * Mellema et al. (2013) Mellema G., et al., 2013, Exp. Astron., 36, 235
  * Mellema et al. (2015) Mellema G., Koopmans L., Shukla H., Datta K. K., Mesinger A., Majumdar S., 2015, arXiv preprint arXiv:1501.04203
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Miralda-Escudé et al. (2000) Miralda-Escudé J., Haehnelt M., Rees M. J., 2000, The Astrophysical Journal, 530, 1
  * Mirocha & Furlanetto (2019) Mirocha J., Furlanetto S. R., 2019, Monthly Notices of the Royal Astronomical Society, 483, 1980
  * Mitra et al. (2015) Mitra S., Choudhury T. R., Ferrara A., 2015, Monthly Notices of the Royal Astronomical Society: Letters, 454, L76
  * Modi et al. (2019) Modi C., Castorina E., Feng Y., White M., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 024–024
  * Mondal et al. (2020) Mondal R., et al., 2020, MNRAS, 498, 4178
  * Morales & Wyithe (2010) Morales M. F., Wyithe J. S. B., 2010, Annual review of astronomy and astrophysics, 48, 127
  * Morishita et al. (2023) Morishita T., et al., 2023, The Astrophysical Journal Letters, 947, L24
  * Napolitano et al. (2024) Napolitano L., et al., 2024, arXiv e-prints, p. arXiv:2402.11220
  * Nebrin et al. (2019) Nebrin O., Ghara R., Mellema G., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 051
  * Nebrin et al. (2023) Nebrin O., Giri S. K., Mellema G., 2023, Monthly Notices of the Royal Astronomical Society, p. stad1852
  * Newburgh et al. (2016) Newburgh L., et al., 2016, in Ground-based and Airborne Telescopes VI. pp 2039–2049
  * Obuljen et al. (2018) Obuljen A., Castorina E., Villaescusa-Navarro F., Viel M., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 004
  * Oesch et al. (2016) Oesch P., et al., 2016, The Astrophysical Journal, 819, 129
  * Oke (1974) Oke J. B., 1974, ApJS, 27, 21
  * Ouchi et al. (2010) Ouchi M., et al., 2010, The Astrophysical Journal, 723, 869
  * Padmanabhan et al. (2017) Padmanabhan H., Refregier A., Amara A., 2017, Monthly Notices of the Royal Astronomical Society, 469, 2323
  * Padmanabhan et al. (2020) Padmanabhan H., Refregier A., Amara A., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3935
  * Park et al. (2019) Park J., et al., 2019, MNRAS, 484, 933
  * Park et al. (2020) Park J., Gillet N., Mesinger A., Greig B., 2020, Monthly Notices of the Royal Astronomical Society, 491, 3891
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, Astronomy & Astrophysics, 641, A6
  * Potter et al. (2017) Potter D., Stadel J., Teyssier R., 2017, Computational Astrophysics and Cosmology, 4, 2
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys., 75, 086901
  * Puchwein et al. (2019) Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2019, Monthly Notices of the Royal Astronomical Society, 485, 47
  * Qin et al. (2021a) Qin Y., Mesinger A., Greig B., Park J., 2021a, Monthly Notices of the Royal Astronomical Society, 501, 4748
  * Qin et al. (2021b) Qin Y., Mesinger A., Bosman S. E., Viel M., 2021b, Monthly Notices of the Royal Astronomical Society, 506, 2390
  * Robertson et al. (2023) Robertson B., et al., 2023, Nature Astronomy
  * Rojas-Ruiz et al. (2020) Rojas-Ruiz S., Finkelstein S. L., Bagley M. B., Stevans M., Finkelstein K. D., Larson R., Mechtley M., Diekmann J., 2020, The Astrophysical Journal, 891, 146
  * Ross et al. (2017) Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, Monthly Notices of the Royal Astronomical Society, 468, 3785
  * Ross et al. (2021) Ross H. E., et al., 2021, MNRAS, 506, 3717
  * Santos et al. (2015) Santos M. G., et al., 2015, arXiv preprint arXiv:1501.03989
  * Saul et al. (2012) Saul D. R., et al., 2012, ApJ, 758, 44
  * Schaeffer et al. (2023) Schaeffer T., Giri S. K., Schneider A., 2023, MNRAS, 526, 2942
  * Schneider (2018) Schneider A., 2018, Phys. Rev. D, 98, 063021
  * Schneider et al. (2021) Schneider A., et al., 2021, Phys. Rev. D, 103, 083025
  * Schneider et al. (2023) Schneider A., Schaeffer T., Giri S. K., 2023, Physical Review D, 108, 043030
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M., 2016, Monthly Notices of the Royal Astronomical Society, 460, 3624
  * Shukla et al. (2016) Shukla H., Mellema G., Iliev I. T., Shapiro P. R., 2016, Monthly Notices of the Royal Astronomical Society, 458, 135
  * Sitwell et al. (2014) Sitwell M., Mesinger A., Ma Y.-Z., Sigurdson K., 2014, Monthly Notices of the Royal Astronomical Society, 438, 2664
  * Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, Monthly Notices of the Royal Astronomical Society, 440, 1662
  * Spinelli et al. (2020) Spinelli M., Zoldan A., De Lucia G., Xie L., Viel M., 2020, Monthly Notices of the Royal Astronomical Society, 493, 5434–5455
  * Stefanon et al. (2019) Stefanon M., et al., 2019, The Astrophysical Journal, 883, 99
  * Sun et al. (2023) Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., Wetzel A., Cochrane R. K., 2023, The Astrophysical Journal Letters, 955, L35
  * The CHIME Collaboration et al. (2023) The CHIME Collaboration et al., 2023, The Astrophysical Journal, 947, 16
  * The HERA Collaboration et al. (2022a) The HERA Collaboration et al., 2022a, The Astrophysical Journal, 924, 51
  * The HERA Collaboration et al. (2022b) The HERA Collaboration et al., 2022b, The Astrophysical Journal, 925, 221
  * The HERA Collaboration et al. (2023) The HERA Collaboration et al., 2023, The Astrophysical Journal, 945, 124
  * Trac et al. (2015) Trac H., Cen R., Mansfield P., 2015, The Astrophysical Journal, 813, 54
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Umeda et al. (2023) Umeda H., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Isobe Y., Zhang Y., 2023, arXiv preprint arXiv:2306.00487
  * Van Der Walt et al. (2011) Van Der Walt S., et al., 2011, Comput. Sci. Eng., 13, 22
  * Villaescusa-Navarro et al. (2014) Villaescusa-Navarro F., Viel M., Datta K. K., Choudhury T. R., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 050
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro F., et al., 2018, The Astrophysical Journal, 866, 135
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature methods, 17, 261
  * Wang et al. (2020) Wang F., et al., 2020, The Astrophysical Journal, 896, 23
  * Worseck et al. (2014) Worseck G., et al., 2014, Monthly Notices of the Royal Astronomical Society, 445, 1745
  * Wu & Zhang (2022) Wu P.-J., Zhang X., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 060
  * Wyithe & Bolton (2011) Wyithe J. S. B., Bolton J. S., 2011, Monthly Notices of the Royal Astronomical Society, 412, 1926
  * Wyithe & Loeb (2008) Wyithe J. S. B., Loeb A., 2008, Monthly Notices of the Royal Astronomical Society, 383, 606
  * Xu et al. (2016) Xu Y., Hamann J., Chen X., 2016, Physical Review D, 94, 123518
  * Xu et al. (2017) Xu Y., Yue B., Chen X., 2017, The Astrophysical Journal, 844, 117
  * Xu et al. (2019) Xu W., Xu Y., Yue B., Iliev I. T., Trac H., Gao L., Chen X., 2019, Monthly Notices of the Royal Astronomical Society, 490, 5739
  * Yung et al. (2024) Yung L. A., Somerville R. S., Finkelstein S. L., Wilkins S. M., Gardner J. P., 2024, Monthly Notices of the Royal Astronomical Society, 527, 5929
  * Zackrisson et al. (2013) Zackrisson E., Inoue A. K., Jensen H., 2013, The Astrophysical Journal, 777, 39
  * Zackrisson et al. (2017) Zackrisson E., et al., 2017, The Astrophysical Journal, 836, 78
  * Zhu et al. (2023) Zhu Y., et al., 2023, The Astrophysical Journal, 955, 115

##  Appendix A Validating the halo catalogue

Our simulation framework, pyC2Ray, utilizes the dark matter haloes to model
photon sources in our cosmological simulations. In this study, we have used
the friend-of-friends algorithm implemented in Pkdgrav3 to find dark matter
haloes in snapshots from $N$-body simulation. As we wanted to model
reionization at large scales ($\gtrsim$100 Mpc) caused by sources residing in
haloes as small as $\sim$10${}^{9}M_{\odot}$, we had to find a balance between
the box size and the number of particles. To ensure this halo catalogue, we
compare the corresponding halo mass function (HMF) against the ones modelled
with extended Press-Schechter (EPS) calculation.

Figure 8: The halo mass function calculated from the halo catalogue of our
$N$-body simulations (points) and modelled using extended Press-Schechter
formalism (solid lines). The redshift values are marked with the colour of the
lines.

We follow the EPS formalism described in Schneider (2018) and refer interested
readers to this work and the references therein. In Figure 8, we observe a
generally robust agreement among the HMFs across all redshifts. However, some
discrepancies are noticeable, particularly at the high-mass end. This
discrepancy is attributed to the substantial Poisson noise inherent in
modelling the massive haloes within the selected simulation volume.