# Probing the intergalactic medium during the Epoch of Reionization using 21 cm signal power spectra

R. Ghara<sup>1,2,3,4</sup>①, A. K. Shaw<sup>2,5</sup>①, S. Zaroubi<sup>1,6,7</sup>①, B. Ciardi<sup>7</sup>, G. Mellema<sup>8</sup>①, L. V. E. Koopmans<sup>1</sup>①, A. Acharya<sup>7</sup>①, M. Choudhury<sup>2,9</sup>①, S. K. Giri<sup>10</sup>①, I. T. Iliev<sup>11</sup>①, Q. Ma<sup>12</sup>, and F. G. Mertens<sup>13</sup>

Received 31 January 2024 / Accepted 16 April 2024

## ABSTRACT

**Context.** The redshifted 21 cm signal from the Epoch of Reionization (EoR) directly probes the ionization and thermal states of the intergalactic medium during that period. In particular, the distribution of the ionized regions around the radiating sources during EoR introduces scale-dependent features in the spherically averaged EoR 21 cm signal power spectrum.

**Aims.** The goal is to study these scale-dependent features at different stages of reionization using numerical simulations and to build a source model-independent framework to probe the properties of the intergalactic medium using EoR 21 cm signal power spectrum measurements.

**Methods.** Under the assumption of high spin temperature, we modeled the redshift evolution of the ratio of the EoR 21 cm brightness temperature power spectrum to the corresponding density power spectrum using an ansatz consisting of a set of redshift and scale-independent parameters. This set of eight parameters probes the redshift evolution of the average ionization fraction and the quantities related to the morphology of the ionized regions.

**Results.** We tested this ansatz on different reionization scenarios generated using different simulation algorithms and found that it is able to recover the redshift evolution of the average neutral fraction within an absolute deviation $\le 0.1$.

**Conclusions.** Our framework allows us to interpret 21 cm signal power spectra in terms of parameters related to the state of the IGM. This source model-independent framework is able to efficiently constrain reionization scenarios using multi-redshift power spectrum measurements with ongoing and future radio telescopes such as LOFAR, MWA, HERA, and SKA. This will add independent information regarding the EoR IGM properties.

**Key words.** radiative transfer – galaxies: formation – galaxies: high-redshift – intergalactic medium – cosmology: theory – dark ages, reionization, first stars

## 1. Introduction

The formation of the first sources of radiation at the end of the Universe's Dark Age is one of the landmark events in Cosmic history. During the first billion years, radiation from the first stars, galaxies, quasars (QSOs), and high-mass X-ray binaries (HMXBs) permanently changed the ionization and thermal state of the Universe. It is expected that radiation from early X-ray sources such as HMXBs and mini-QSOs changed the thermal state of the cold intergalactic medium (IGM) much before the IGM became highly ionized (see, e.g., Pritchard & Furlanetto 2007; Thomas & Zaroubi 2011;

Mesinger et al. 2011; Islam et al. 2019; Ross et al. 2019; Eide et al. 2020). The onset of the first sources that changed the IGM's thermal state is known as the Cosmic Dawn (CD). The subsequent period when the IGM's atomic neutral hydrogen (Hα) became ionized is known as the Epoch of Reionization (EoR). A few indirect probes such as the observations of Gunn-Peterson optical depth in $z \ge 6$ QSO spectra and Thomson scattering optical depth of the cosmic microwave background (CMB) photons provide us with useful information about the rough timing and duration of the EoR (see, e.g., Fan et al. 2006; McGreer et al. 2015; Bañados et al. 2018; Planck Collaboration VI 2020; Mitra et al. 2015).

<footnote n="1">¹ Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
e-mail: ghara.raghunath@gmail.com</footnote>

<footnote n="2">² ARCO (Astrophysics Research Center), Department of Natural Sciences, The Open University of Israel, 1 University Road, PO Box 808, Ra'anana 4353701, Israel
e-mail: abinashkumarshaw@gmail.com</footnote>

<footnote n="3">³ Haverford College, 370 Lancaster Ave, Haverford, PA 19041, USA</footnote>

<footnote n="4">⁴ Center for Particle Cosmology, Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="5">⁵ Department of Computer Science, University of Nevada, Las Vegas, Nevada 89154, USA</footnote>

<footnote n="6">⁶ Department of Natural Sciences, The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel</footnote>

<footnote n="7">⁷ Max-Planck-Institute for Astrophysics, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany</footnote>

<footnote n="8">⁸ The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova 10691, Stockholm, Sweden</footnote>

<footnote n="9">⁹ Center for Fundamental Physics of the Universe, Department of Physics, Brown University, Providence 02914, RI, USA</footnote>

<footnote n="10">¹⁰ Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvén's väg 12, 106 91 Stockholm, Sweden</footnote>

<footnote n="11">¹¹ Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Pevensey II Building, Brighton BN1 9QH, UK</footnote>

<footnote n="12">¹² School of Physics and Electronic Science, Guizhou Normal University, Guiyang 550001, PR China</footnote>

<footnote n="13">¹³ LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, 75014 Paris, France</footnote>
However, many details about these epochs, such as the exact timing, properties of the sources and their evolution, feedback mechanisms, and morphology of the ionized and heated regions, are still unknown.

Observations of the redshifted 21 cm radiation produced by HI in the IGM can provide us with information related to the timing, the morphology of the ionized and heated regions, and properties of the ionizing and heating sources (see, e.g., Pritchard & Loeb 2012; Zaroubi 2013; Shaw et al. 2023a; Ghara et al. 2024, for reviews). Many of the world's large radio observation facilities have aimed for measuring the brightness temperature of this redshifted HI 21 cm radiation (hereafter 21 cm signal) from the CD and EoR. Radio observations using a single antenna, such as EDGES2 (Bowman et al. 2018), SARAS2 (Singh et al. 2017), REACH (de Lera Acedo et al. 2022), and LEDA (Price et al. 2018), aim to measure the redshift evolution of the sky-averaged 21 cm signal. However, observing the morphological distribution of the 21 cm signal in the sky is expected to tell us more about these epochs. Radio interferometers, such as the Low-Frequency Array (LOFAR)<fnref n="1" /> (van Haarlem et al. 2013; Patil et al. 2017), the New Extension in Nançay Upgrading LOFAR (NenuFAR)<fnref n="2" /> (Munshi et al. 2024), the Amsterdam ASTRON Radio Transients Facility and Analysis Center (AARTFAAC) (Gehlot et al. 2022), the Precision Array for Probing the Epoch of Reionization (PAPER)<fnref n="3" /> (Parsons et al. 2014; Kolopanis et al. 2019), the Murchison Widefield Array (MWA)<fnref n="4" /> (e.g. Tingay et al. 2013; Wayth et al. 2018), and the Hydrogen Epoch of Reionization Array (HERA)<fnref n="5" /> (DeBoer et al. 2017), have been commissioned to measure the spatial fluctuations in the HI 21 cm signal at different stages of the CD and EoR.

Due to limited sensitivity, the radio interferometer-based observations aim to detect this signal in terms of the statistical quantities such as the spherically averaged power spectrum ($\Delta^2_{\delta T_b}(k, z)$) of the differential brightness temperature ($\delta T_b$) of the HI signal at different redshifts ($z$) and wave-numbers ($k$). The upcoming Square Kilometre Array (SKA)<fnref n="6" /> will be more sensitive and will also produce tomographic images of the CD and EoR 21 cm signal (Mellema et al. 2015; Ghara et al. 2017).

Observing the 21 cm signal from CD and EoR is very challenging and it has remained undetected by the radio observations to date. The measured HI signal is severely contaminated by the galactic and extra-galactic foregrounds. While the foregrounds are more substantial than the expected CD and EoR HI signal by several orders of magnitude (see, e.g., Ghosh et al. 2012), their smooth frequency dependence allows them to be either subtracted (Harker et al. 2009; Bonaldi & Brown 2015; Chapman et al. 2016; Mertens et al. 2018; Hothi et al. 2021), avoided (Datta et al. 2010; Liu et al. 2014) or suppressed (Datta et al. 2007; Ghara et al. 2016). These observations also face severe challenges at the calibration step of the data analysis process. Nevertheless, recent improvements in the calibration methods (see, e.g., Kern et al. 2019, 2020; Mevius et al. 2022; Gan et al. 2022, 2023), and the mitigation of the foregrounds (e.g., Mertens et al. 2018; Liu et al. 2014) made it possible to obtain noise dominated upper limits of $\Delta^2_{\delta T_b}(k, z)$. For example, $\Delta^2_{\delta T_b}(k = 0.14 \, h \, \text{Mpc}^{-1}, z = 6.5) \approx (43)^2 \, \text{mK}^2$,

$\Delta^2_{\delta T_b}(k = 0.075 \, h \, \text{Mpc}^{-1}, z = 9.1) \approx (73)^2 \, \text{mK}^2$, and $\Delta^2_{\delta T_b}(k = 0.34 \, h \, \text{Mpc}^{-1}, z = 7.9) \approx (21.4)^2 \, \text{mK}^2$ are the best upper limits obtained from MWA (Trott et al. 2020), LOFAR (Mertens et al. 2020), and HERA (Abdurashidova et al. 2023) EoR observations, respectively.

These recent upper limits have started to rule out CD and EoR scenarios including those which do not require either an unconventional cooling mechanism or the presence of a strong radio background in addition to the CMB (e.g., Ghara et al. 2020; Greig et al. 2021; Mondal et al. 2020; Abdurashidova et al. 2022a). For example, the recent HERA EoR observation results as reported in Abdurashidova et al. (2022a,b) show that the IGM temperature must be higher than the adiabatic cooling threshold by redshift 8, while the soft band X-ray luminosity per star formation rate of the first galaxies are constrained ($1\sigma$ level) to $[10^{40.2}-10^{41.9}] \, \text{erg s}^{-1}/(\text{M}_\odot \, \text{yr}^{-1})$. In addition, the recent results from the global HI 21 cm signal observations, such as SARAS and EDGES, have also started ruling out EoR and CD scenarios and putting constraints on the properties of the early sources, models of dark matter, and level of radio backgrounds (e.g., Barkana 2018; Fialkov et al. 2018; Muñoz & Loeb 2018; Nebrin et al. 2019; Chatterjee et al. 2019, 2020; Ghara & Mellema 2020; Ghara et al. 2022; Bera et al. 2023).

These previous studies have put constraints mainly on the astrophysical source properties using either Bayesian inference techniques (e.g., Park et al. 2019; Cohen et al. 2020) or Fisher matrices (e.g., Ewall-Wice et al. 2016; Shaw et al. 2020). The main reason behind this is the fact that 21 cm signal simulation codes take the source parameters as input. However, it should be realized that the 21 cm signal measurements do not probe the astrophysical sources directly. In addition, the inference on the properties of the astrophysical sources is limited by the ambiguity of the source model used in the inference framework. The observed 21 cm signal, on the other hand, directly probes the ionization and the thermal states of the IGM. Therefore, we emphatically aim to constrain the IGM properties rather than the astrophysical source parameters.

Previously, Mirocha et al. (2013) considered the features of the redshift evolution of the sky-averaged brightness temperature curves within a simplified global HI signal framework which does not invoke any astrophysical sources and attempted to constrain physical properties of the IGM in terms of Ly$\alpha$ background, overall heat deposition, mean ionization fraction, and their time derivatives. In the context of 21 cm signal power spectrum, studies such as Ghara et al. (2020, 2021) used the recently obtained upper limits from LOFAR (Mertens et al. 2020) and MWA (Trott et al. 2020) to constrain the properties of the IGM at different stages of the EoR. These studies use the outputs from GRIZZLY (Ghara et al. 2015a) simulations and characterize the IGM in terms of quantities such as the sky-averaged ionization fraction, average gas temperature, sky-averaged brightness temperature, the volume fraction of the “heated regions” in the IGM with its brightness temperature $T_b$ larger than the background CMB temperature $T_\gamma$, and the characteristic size of these heated regions. For example, using the recent upper limits from LOFAR (Mertens et al. 2020), Ghara et al. (2020) ruled out reionization scenarios at redshift 9.1 where heating of the gas is negligible and the IGM is characterized by ionized fraction $\ge 0.13$, a distribution of the ionized regions with a characteristic size $\ge 8 \, h^{-1} \, \text{Mpc}$, and a full width at half-maximum $\ge 16 \, h^{-1} \, \text{Mpc}$. In an alternative approach, Shimabukuro et al. (2022) used artificial neural networks to build a framework that estimates the size

<footnote n="1">¹ http://www.lofar.org/</footnote>

<footnote n="2">² https://nenufar OBS-Nancay.fr/en/homepage-en/</footnote>

<footnote n="3">³ http://eor.berkeley.edu/</footnote>

<footnote n="4">⁴ http://www.mwatelescope.org/</footnote>

<footnote n="5">⁵ https://reionization.org/</footnote>

<footnote n="6">⁶ http://www.skatelescope.org/</footnote>
distribution of the ionized regions using the EoR 21 cm power spectrum.

Our previous studies such as Ghara et al. (2020, 2021), which aim at constraining the properties of the CD and EoR IGM parameters, use source-parameter dependent GRIZZLY simulations. The inputs of their framework are a set of source parameters, such as the ionization efficiency, the minimum mass of dark matter halos that host UV emitting sources, the X-ray emission efficiency, and the minimum mass of dark matter halos that host X-ray emitting sources. The framework provides a set of derived IGM parameters in addition to the 21 cm signal observable. It is not straightforward to build a mathematical framework that directly connects the complex morphology of the IGM to the 21 cm signal observable by skipping the source-parameter dependence. Recently, Mirocha et al. (2022) have attempted to build such a galaxy-free phenomenological model for the EoR 21 cm signal power spectra. The model assumes uniform $T_S$, spherical ionized bubbles and a binary ionization field. While the model efficiently predicts the 21 cm signal power spectrum for volume average neutral fraction $\bar{x}_{HI} \ge 0.8$, the prediction accuracy rapidly drops for reionization stages with $\bar{x}_{HI} \le 0.8$ which shows the complexity level of the problem.

Unlike our aforementioned IGM inference framework, the main goal of this work is to develop a source parameter-free phenomenological model of EoR 21 cm signal power spectra in terms of quantities related to the IGM. We keep our model simple by ignoring the effect of spin-temperature fluctuations and targeting the IGM only during the EoR. The amplitude and the shape of $\Delta_{\delta T_b}^2(k, z)$ as a function of $k$ during different stages of the EoR depend on the ionization fraction and the complex morphology of the ionized regions at that period. The aim here is to use the multi-redshift measurements of the EoR 21 cm signal power spectra to constrain the IGM properties during the EoR.

This paper is structured as follows. In Sect. 2, we describe the basic methodology of our framework. We present our results in Sect. 3, before concluding in Sect. 4. The cosmological parameters used throughout this study are the same as the N-body simulations employed here (i.e., $\Omega_m = 0.27$, $\Omega_\Lambda = 0.73$, $\Omega_B = 0.044$, $h = 0.7$ (Wilkinson Microwave Anisotropy Probe (WMAP); Hinshaw et al. 2013).

## 2. Framework

### 2.1. The EoR 21 cm signal

The differential brightness temperature ($\delta T_b$) of the 21 cm signal from a region at angular position $\mathbf{x}$ and redshift $z$ can be expressed as (see, e.g., Madau et al. 1997; Furlanetto et al. 2006),

$$ \delta T_b(\mathbf{x}, z) = 27 x_{HI}(\mathbf{x}, z) \left[ 1 + \delta_B(\mathbf{x}, z) \left( \frac{\Omega_B h^2}{0.023} \right) \times \left( \frac{0.15}{\Omega_m h^2} \frac{1+z}{10} \right)^{1/2} \left[ 1 - \frac{T_\gamma(z)}{T_S(\mathbf{x}, z)} \right] \right] \text{mK.} \quad (1) $$

Here $T_S$, $x_{HI}$, and $\delta_B$ are respectively the spin temperature of HI, the neutral hydrogen fraction, and the baryonic density contrast of the region located at $(\mathbf{x}, z)$. The quantity $T_\gamma$ is the radio background temperature at 21 cm wavelength for redshift $z$. In this study, we assume a high spin temperature limit, i.e., $T_S \gg T_\gamma$. This is expected to be the case in the presence of efficient X-ray heating.

Here, we use the GRIZZLY code (Ghara et al. 2015a,b) to generate brightness temperature maps during the EoR. The inputs

for this code are the uniformly gridded dark-matter density and velocity field cubes and the corresponding dark-matter halo list. The GRIZZLY simulations considered in this study use the dark-matter fields and the corresponding halo lists within comoving cubes of side $500 h^{-1}$ Mpc, produced from the PRACE$^7$ project PRACE4LOFAR N-body simulations (see e.g., Giri et al. 2019; Kamran et al. 2021, for the details of the simulation). The simulation assumes that all dark-matter halos with masses larger than $M_{min}$ contribute to reionization. The stellar mass $M_\star$ inside a halo of mass $M_{halo}$ is assumed to be $M_\star \propto M_{halo}^{\alpha_s}$. We choose the ionization efficiency ($\zeta$) so that the reionization process ends roughly at $z \sim 6.5^8$. The reionization models considered in this study are inside-out in nature where the very dense regions around the sources get ionized first. We refer to Ghara et al. (2015a, 2020) for the details of the method and the source parameters. Our fiducial GRIZZLY model, as shown in Fig. 1, corresponds to a choice of $M_{min} = 10^9 M_\odot$ and $\alpha_s = 1$ and spans from redshift 6.9 to 15.5. We also produce 23 more reionization scenarios by choosing different combinations of [$M_{min}, \alpha_s$] where we vary $M_{min}$ between $10^9-10^{11} M_\odot$ and $\alpha_s$ between 0.3–2. Smaller values of $M_{min}$ and larger values of $\alpha_s$ will create a more patchy reionization scenario. We use all these reionization scenarios for building and testing our model of the $\delta T_b$ power spectrum. We note that all our simulations include the redshift-space distortion effects based on the cell moving method (Ghara et al. 2015a; Ross et al. 2021).

Figure 1 shows a slice through a simulated light-cone of the EoR 21 cm signal $\delta T_b$ (Eq. (1)). The figure represents how the fluctuations $\delta T_b$ in the sky (shown by the vertical axis) evolve with redshift or distance from the observer (represented by the horizontal axis). Our assumption of $T_S \gg T_\gamma$ makes the HI 21 cm signal $\delta T_b$ positive in the neutral regions while the ionized regions are represented by $\delta T_b = 0$. We note that ionized regions in the IGM are absent around $z = 15$ where the fluctuations in $\delta T_b$ are governed by the density fluctuations only (Eq. (1)). Small isolated ionized regions gradually appear around the high-density peaks in $\delta_B$. Over time, the isolated ionized regions grow in size and overlap with each other. This overlap can occur as early as when the IGM volume is ionized by a few tens of percent depending upon the reionization history. These overlaps eventually create complex percolated structures of the ionized regions which grow in volume over time as the reionization progresses. For $\bar{x}_{HI} \ge 0.5$ (i.e., $z \ge 8$ in Fig. 1), the sizes of the ionized regions are smaller than the neutral regions. Visually, the ionized regions are embedded into the neutral regions. It becomes the opposite at reionization stages with $\bar{x}_{HI} \le 0.5$. At these stages, the distribution of the neutral regions is more meaningful compared to the distribution of the ionized regions.

### 2.2. The EoR 21 cm signal power spectrum

This study is based on the k-dependent features of the dimensionless power spectrum of coeval $\delta T_b$ cube at redshift $z$, i.e., $\Delta_{\delta T_b}^2(z, k) = k^3 P_{\delta T_b}(z, k)/2\pi^2$ during different stages of the EoR. Assuming statistical homogeneity of the signal, one can define the 3D power spectrum for a coeval signal volume $V$ as $\delta_K(k -$

<footnote>$^7$ Partnership for Advanced Computing in Europe: http://www.prace-ri.eu/</footnote>

<footnote>$^8$ Some probes, such as the Ly$\alpha$ forest observations at $z \approx 5.5$, suggest a late reionization compared to our fiducial reionization model (see, e.g., Becker et al. 2018; Eilers et al. 2018). However, the exact end of the EoR is still debated and hence, we choose our fiducial simulation from our earlier works (e.g., Shaw et al. 2023b; Ghara et al. 2024).</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Light-cone of EoR 21 cm signal. This shows the redshift ($z$) evolution of the corresponding differential brightness temperature ($\delta T_b$) from $z \approx 15.5$ to 6.9. We assume a high spin temperature limit, i.e., $T_S \gg T_\gamma$. This light cone is generated using the GRIZZLY code and represents our fiducial EoR scenario.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Power spectra of EoR 21 cm signal brightness temperature and their different components. The a–d panels show $\Delta_{\delta T_b}^2$, $\Delta_{\delta\delta}^2$, $\Delta_{x_{HI}x_{HI}'}^2$, and $|\Delta_{x_{HI}\delta}^2|$ respectively as a function of $k$ at different stages of reionization. The e panel shows the redshift evolution of the 21 cm power spectrum at different scales. The f panel compares the redshift evolution of $\Delta_{x_{HI}x_{HI}'}^2$, $\Delta_{\delta\delta}^2$, and $|\Delta_{x_{HI}\delta}^2|$ at $k = 0.1 \text{ h Mpc}^{-1}$. The power spectra correspond to our fiducial GRIZZLY EoR scenario as shown in Fig. 1.</figcaption></figure>

$k')P_{\delta T_b}(z, k) = V^{-1}\langle\delta\hat{T}_b(z, k) \delta\hat{T}_b^*(z, k')\rangle$, where $\delta_K(k - k')$ denotes the 3D Kroneker's delta function and $\delta\hat{T}_b(z, k)$ is the Fourier transform of the EoR 21 cm signal $\delta T_b(x, z)$. Here we use the spherically averaged power spectrum $P_{\delta T_b}(z, k)$ which is computed by averaging $P_{\delta T_b}(z, k)$ within spherical shells of certain widths in 3D Fourier space. According to Eq. (1), the EoR 21 cm power spectrum $\Delta_{\delta T_b}^2(z, k)$ depends on the power spectra of the density and neutral fraction fields ($\Delta_{\delta\delta}^2$ and $\Delta_{x_{HI}x_{HI}'}^2$ respectively), and their cross power spectrum $\Delta_{x_{HI}\delta}^2$. Here, $\Delta_{\delta\delta}^2$ is associated with a field given by Eq. (1) for $x_{HI}(x, z) = 1$, therefore powered by the density fluctuations only. On the other hand, the field associated with $\Delta_{x_{HI}x_{HI}'}^2$ assumes $\delta_B(x, z) = 0$ in Eq. (1), and thus is independent of the density fluctuations and only depends on the neutral fraction fluctuations (Lidz et al. 2007; Georgiev et al. 2022)<fnref n="9" />. $\Delta_{x_{HI}\delta}^2$ is the cross-power spectrum of the fields associated with $\Delta_{\delta\delta}^2$ and $\Delta_{x_{HI}x_{HI}'}^2$.

Figure 2 shows the evolution of the EoR 21 cm signal power spectrum $\Delta_{\delta T_b}^2$ as well as the power spectra of the density field $\Delta_{\delta\delta}^2$, neutral fraction field $\Delta_{x_{HI}x_{HI}'}^2$, and their cross-power spectrum $\Delta_{x_{HI}\delta}^2$. The power spectra correspond to our fiducial GRIZZLY EoR scenario as presented in Fig. 1. The a–d panels of Fig. 2 show $\Delta_{\delta T_b}^2$, $\Delta_{\delta\delta}^2$, $\Delta_{x_{HI}x_{HI}'}^2$, and $|\Delta_{x_{HI}\delta}^2|$ respectively as a function of $k$ at different redshifts. The panel e shows the redshift evolution of $\Delta_{\delta T_b}^2$ for different scales while panel f compares the redshift evolution of the $\Delta_{x_{HI}x_{HI}'}^2$, $\Delta_{\delta\delta}^2$, and $|\Delta_{x_{HI}\delta}^2|$ for $k = 0.1 \text{ h Mpc}^{-1}$. The high-density regions get ionized first in this inside-out reionization model. This causes anti-correlation between $x_{HI}$ and $\delta$ and thus, negative values for the cross-power spectrum $\Delta_{x_{HI}\delta}^2$ which suppresses the large-scale $\delta T_b$ power spectrum at the initial stage of the EoR<fnref n="10" />. However, the suppression is less significant at the small-scales, causing a tilt in the $\Delta_{\delta T_b}^2$ compared to the $\Delta_{\delta\delta}^2$ (see

<footnote n="9">⁹ Lidz et al. (2007), Georgiev et al. (2022) use the field of fluctuations in $x_{HI}$, $\delta x_{HI}$. This leads to additional, higher-order, cross terms when composing the 21 cm power spectra in terms of the constituent fields $\delta$ and $\delta x_{HI}$.</footnote>

<footnote n="10">¹⁰ This phase of strong suppression of the large scale 21 cm power spectrum was first pointed out by Lidz et al. (2007) who called it “equilibration”. Georgiev et al. (2022) studied it in some more detail and showed that it is caused by the near cancellation of the positive and negative terms in the decomposition of the 21 cm power spectra.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 3. Ratio of $\Delta_{\delta T_b}^2$ to $\Delta_{\delta\delta}^2$ obtained from simulations during EoR. The top panel shows $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ as a function of $k$ at different redshifts. These correspond to our fiducial GRIZZLY model which is also shown in Fig. 2. The bottom panel shows evolution of $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ at $k = 0.05 \, h \, \text{Mpc}^{-1}$ as a function of neutral fraction. Different curves stand for different reionization scenarios generated using GRIZZLY by varying $\alpha_s$ and $M_{\text{min}}$. The black curve represents our fiducial GRIZZLY simulation which corresponds to $\alpha_s = 1$ and $M_{\text{min}} = 10^9 M_{\odot}$.</figcaption></figure>

panels a and b). For $z \ge 9$, $|\Delta_{x_{HI}\delta}^2|$ remains larger than $\Delta_{x_{HI}x_{HI}}^2$ (see panels d and f). For $z \le 9$, $\Delta_{x_{HI}x_{HI}}^2$ becomes the dominant term. The interplay between the $\Delta_{x_{HI}x_{HI}}^2$ and $\Delta_{x_{HI}\delta}^2$ contributions causes a minimum in the $\Delta_{\delta T_b}^2$ vs $z$ curves around $z \sim 10$ (see bottom panels). The large-scale $\delta T_b$ power spectrum increases as reionization progresses. For example, $\Delta_{\delta T_b}^2(k = 0.1 \, h \, \text{Mpc}^{-1})$ increases from $z = 9$ to $z = 7.3$ as $\Delta_{x_{HI}x_{HI}}^2$ becomes dominant compared to $\Delta_{x_{HI}\delta}^2$. For $z \le 7.3$, $\Delta_{\delta T_b}^2(k = 0.1 \, h \, \text{Mpc}^{-1})$ quickly drops as the majority of the IGM gets ionized. This causes $\Delta_{\delta T_b}^2(k = 0.1 \, h \, \text{Mpc}^{-1})$ to peak at redshift 7.3 with amplitude of $\approx 10 \, \text{mK}^2$. The peak amplitudes and the associated redshifts change with $k$ (see panel e of Fig. 2).

The top panel of Fig. 3 shows the ratio of $\Delta_{\delta T_b}^2$ to $\Delta_{\delta\delta}^2$ (also known as the 21 cm signal bias) as a function of $k$ at different redshifts for the fiducial GRIZZLY model. The bottom panel of Fig. 3 shows the evolution of $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ for $k = 0.05 \, h \, \text{Mpc}^{-1}$ as a function of $\bar{x}_{HI}$. The different curves in the bottom panels correspond to different GRIZZLY reionization models, with the thick black curve representing the fiducial one. This ratio is expected to be 1 for $\bar{x}_{HI} = 1$. The curves show that the ratio first decreases from 1 to a minimum at an early stage of reioniza-

tion. We denote $\bar{x}_{HI}$ at this stage as $\bar{x}_{HI,\text{min}}$. This minimum corresponds to the equilibration phase caused by the anti-correlation between the density and neutral fraction fields in our inside-out reionization model where $\Delta_{x_{HI}\delta}^2 < 0$. The ratio then increases with the increase of the size of the ionized regions and reaches a maximum and further decreases to zero as $\bar{x}_{HI}$ approaches zero towards the end of the EoR. We denote the ionization fraction at the stage when the maximum occurs as $\bar{x}_{HI,\star}$. The qualitative features of the different curves in the bottom panel of Fig. 3 are similar, although the stages when the minimum and maximum occur (i.e., the values of $\bar{x}_{HI,\text{min}}$ and $\bar{x}_{HI,\star}$) and the peak amplitude of the ratio $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ changes with the patchiness of the reionization scenarios. The curve which reaches the largest bias or ratio value (approximately 14) corresponds to the model which uses the largest values for both $M_{\text{min}}$ ($10^{11} M_{\odot}$) and $\alpha_s$ (2), and thus represents the most patchy reionization scenario among the considered models. Nevertheless, the evolutionary features of this ratio as a function of $\bar{x}_{HI}$ remain the same despite being an early or late reionization.

### 2.3. Modeling scale dependence of the EoR 21 cm signal power spectrum at a given redshift

In this section, we aim to model the complex scale dependence of the $\delta T_b$ power spectrum (e.g., see Figs. 2 and 3) as we described in the previous section. The $k$-dependence of $\Delta_{\delta T_b}^2$ evolves with redshift. It should be noted that this study considers features of the power spectrum for the range $0.05 \, h \, \text{Mpc}^{-1} \le k \le 0.6 \, h \, \text{Mpc}^{-1}$, which covers the scales probed by EoR observations such as LOFAR and MWA. The overall feature of the power spectrum, as we have seen in the top panel of Fig. 3 (see also, Xu et al. 2019; Georgiev et al. 2022), suggests that one possible ansatz to represent the $k$-dependence of the coeval EoR 21 cm signal power spectrum can be

$$ \Delta_{\delta T_b}^2(k, z) = A \frac{\left(\frac{k}{k_0}\right)^\eta \Delta_{\delta\delta}^2(k, z)}{1 + \left(\frac{k}{k_0}\right)^\eta} \quad (2) $$

Here $A, k_c, k_0, \gamma$, and $\eta$ are the parameters to fit $\Delta_{\delta T_b}^2$ at a particular redshift, considering $\Delta_{\delta\delta}^2$ to be known for the background cosmology. The form of the numerator in Eq. (2) is chosen to compensate for the difference in slope between $\Delta_{\delta T_b}^2$ and $\Delta_{\delta\delta}^2$ during the early stages of the EoR (see Fig. 2). On the other hand, the denominator accounts for the fall of $\Delta_{\delta T_b}^2$ relative to the corresponding $\Delta_{\delta\delta}^2$ at the small scales (e.g., for $k \ge 0.1 \, h \, \text{Mpc}^{-1}$) during the advanced stages of the EoR (e.g., for $\bar{x}_{HI} \le 0.7$). For $\bar{x}_{HI} = 1$, one expects $A = 1, \gamma = 0, k_0 \to \infty$ and $\eta$ to be positive. It is expected that the $k_0$ values are much larger compared to the smallest $k$ value achieved by the EoR HI observations. For example, LOFAR reaches $k \sim 0.05 \, h \, \text{Mpc}^{-1}$ as reported in papers such as Patil et al. (2017). We set $k_c = 0.05 \, h \, \text{Mpc}^{-1}$ which makes the $A$ parameter equal to the ratio of $\Delta_{\delta T_b}^2$ to $\Delta_{\delta\delta}^2$ at $k = 0.05 \, h \, \text{Mpc}^{-1}$ if $k_0$ is quite large compared to $0.05 \, h \, \text{Mpc}^{-1}$.

We first attempt to fit the power spectra at different redshifts using Eq. (2), separately, by varying $A, k_0, \gamma$, and $\eta$ on a grid for our fiducial GRIZZLY scenario. We explore the log($A$), $\gamma$, log($k_0$), and $\eta$ parameter space on equal-spaced grids. The chosen parameter ranges are $[-4, 2]$, $[-2, 2]$, $[-1, 0]$, and $[0, 3]$, respectively. We estimate the best-fit values of these four parameters at each redshift by minimizing the error

$$ \Psi_4^2(S: [z, A, \gamma, k_0, \eta]) = \frac{1}{n} \sum_{i=1,n} \left( \frac{\Delta_{\delta T_b, m}^2(S)}{\Delta_{\delta T_b, o}^2(z, k_i)} - 1 \right)^2. $$
<figure><img src="image_4.png" /><figcaption>Fig. 4. Outcome of fitting EoR power spectra using Eq. (2). Left to right panels show the evolution of the best-fit values of the parameters A, γ, k₀ and η and fitting error as a function of x̄ₖ for different reionization scenarios generated using GRIZZLY. Bold lines correspond to the fiducial GRIZZLY reionization scenario.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Outcome of fitting EoR power spectra using Eq. (3). Left to right panels show the evolution of the best-fit values of the parameters A and γ and fitting error as a function of x̄ₖ for different reionization scenarios. It should be noted that here we have used Eq. (3) while we used Eq. (2) for Fig. 4. The reionization scenarios are the same as in Fig. 4. Bold lines correspond to the fiducial GRIZZLY reionization scenario.</figcaption></figure>

We consider a *k* range in between 0.05 and 0.6 *h* Mpc⁻¹ and divide it into *n* = 9 log spaced bins. Here, Δ²ₖₜₘ represents the *modeled* power spectrum using Eq. (2) for a set of input parameters. Δ²ₖₜₐ is the *observed* or *input* power spectrum which we consider for fitting. The evolution of the best-fit parameter values as a function of x̄ₖ is shown in the different panels of Fig. 4. The rightmost panel, which shows the goodness of the fit, indicates that the maximum error in our fitting is well within 10%.

The left panel of Fig. 4 presents the best-fit values of A as a function of x̄ₖ. These roughly agree with the values presented in Fig. 3. We find that k₀ reaches the maximum value of the range or becomes unconstrained (simultaneously η becomes ill-defined) during the early stages of reionization (x̄ₖᵢ ≥ x̄ₖᵢ, min) before significant overlap between the isolated ionized regions occurs. Eventually, the formation of large overlapped ionized regions changes the *k* dependencies on the small-scale power spectra. With the growth of overlapping ionized regions, the characteristic bubble size increases, which implies a decrease in k₀ values while remaining much larger than 0.05 *h* Mpc⁻¹. We repeated the fitting of Eq. (2) for our different GRIZZLY scenarios and found a qualitatively similar dependence of k₀ and η on reionization history (see grey lines in Fig. 4).

It is expected from Eq. (2) that the parameters γ, k₀, and η might be degenerate. To reduce degenerate parameters and to minimize the number of parameters in our ansatz, we fix k₀ and η at typical values of the best-fit parameters throughout the reionization history. We fixed η = 1.5, while we fixed k₀ = 0.3 *h* Mpc⁻¹ for x̄ₖᵢ ≤ x̄ₖᵢ, min and infinity otherwise. This reduces the number of parameters and modifies Eq. (2) to

$$ \Delta_{\delta T_b}^2 = \begin{cases} \Delta_{\delta 0}^2 A \left( \frac{k}{0.05} \right)^{\gamma} \left( 1 + \left( \frac{k}{0.5} \right)^{1.5} \right), & \text{if } \bar{x}_{\text{HI}} \le \bar{x}_{\text{HI, min}} \\ \Delta_{\delta 0}^2 A \left( \frac{k}{0.05} \right)^{\gamma}, & \text{otherwise.} \end{cases} \quad (3) $$

As before, we fixed k_c to 0.05, which ensures that the parameter A represents the ratio Δ²ₖₜₘ/Δ²ₖₐ at *k* = 0.05 *h* Mpc⁻¹.

To check the performance of the form in Eq. (3), we vary the parameters A and γ and compare the fitted power spectrum

Δ²ₖₜₘ(z, k, A, γ) with the simulated input power spectrum Δ²ₖₜₐ(z, k) at each redshift independently. We explore the log(A) and γ parameter space on grids where we have chosen the same parameter ranges as above, i.e., [-4, 2] and [-2, 2], respectively. We estimate the best-fit values of A and γ at each redshifts by minimizing the error

$$ \Psi_2^2(z, A, \gamma) = \frac{1}{n} \sum_{i=1,n} \left( \frac{\Delta_{\delta T_b, m}^2(z, k_i, A, \gamma)}{\Delta_{\delta T_b, o}^2(z, k_i)} - 1 \right)^2 . $$

The two left panels of Fig. 5 present the evolution of the best-fit values of A and γ as a function of x̄ₖ for different GRIZZLY reionization scenarios. The right panel of the figure presents Ψ₂×100 as a function of x̄ₖ. The figure shows that, even with the simplified form of Δ²ₖₜₘ as used in Eq. (3), the fitting error Ψ₂ is ≤10% and not drastically different compared to Ψ₄. On the other hand, the evolution of the γ parameter is now smoother compared to the four-parameter case.

The top panel of Fig. 6 shows a comparison between the power spectrum of our fiducial reionization scenario and its best-fit power spectrum Δ²ₖₜₐ(z, k) obtained using Eq. (3). The bottom panel of Fig. 6 shows the percentage fitting error Err(Δ²ₖₜₘ) = (1 - Δ²ₖₜₐ/Δ²ₖₜₘ) × 100% as a function of *k* modes for the different redshifts. The plot shows that Eq. (3) can predict the power spectrum of the fiducial EoR scenario with error ≤10%. We roughly find similar results for other GRIZZLY reionization scenarios.

## 2.4. Modeling the redshift evolution of the EoR 21 cm signal power spectrum

We already find that the form of the 21 cm power spectrum as used in Eq. (3) works well for different stages of reionization. However, for any given reionization scenario, the best-fit values of the A and γ evolve with redshift or alternatively x̄ₖ.
<figure><img src="image_6.png" /><figcaption>Fig. 6. Comparing simulated and fitted EoR power spectra using Eq. (3). The top panel shows a comparison between the 21 cm signal power spectrum of our fiducial reionization scenario and its best-fit power spectrum obtained using Eq. (3). Different colors represent different redshifts. The lines correspond to the power spectrum $\Delta_{\delta T_b}^2$ from GRIZZLY simulations while the points stand for the fitted power spectrum $\Delta_{\delta T_b, f}^2$. The bottom panel shows the percentage fitting error $\text{Err}(\Delta_{\delta T_b}^2) = (1 - \frac{\Delta_{\delta T_b, f}^2}{\Delta_{\delta T_b}^2}) \times 100\%$.</figcaption></figure>

Thus, it is important to understand the behavior of A and $\gamma$ parameters as a function of $\bar{x}_{HI}$ to come up with a final set of redshift-independent parameters which can be constrained using multi-redshift EoR 21 cm signal power spectra measurements.

The left panel of Fig. 5 shows the dependencies of the best-fit values of A as a function of $\bar{x}_{HI}$. We note that both the bottom panel of Fig. 3 and the left panel of Fig. 5 represent the evolution of $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ at $k = 0.05 \text{ h Mpc}^{-1}$ as a function of $\bar{x}_{HI}$. While the curves in the bottom panel of Fig. 3 show the estimates directly from the GRIZZLY simulation, the curves in the left panel of Fig. 5 represent the best-fit value of parameter A when using Eq. (3) for modeling $\Delta_{\delta T_b}^2$. A = 1 when $\bar{x}_{HI} = 1$ and the 21 cm signal power spectrum is completely determined by density fluctuations. As reionization progresses (i.e., $\bar{x}_{HI}$ decreases), A first decreases and reaches its minimum value at $\bar{x}_{HI} = \bar{x}_{HI,\min}$. Thereafter, A increases until $\bar{x}_{HI} = \bar{x}_{HI,\star}$ where it reaches a maximum value, which we call $A_{\star}$. After this, A decreases and $A \to 0$ for $\bar{x}_{HI} \to 0$ when reionization ends.

We model the dependence of A on $\bar{x}_{HI}$ in the following way.

$$A = A_{\star} \left( \frac{\bar{x}_{HI}}{\bar{x}_{HI,\star}} \right)^{\alpha_A} \left( \frac{1 - \bar{x}_{HI}}{1 - \bar{x}_{HI,\star}} \right)^{\beta_A} + 10^{-5 \times (1 - \bar{x}_{HI})} \quad (4)$$

with $\beta_A = \alpha_A (1 - \bar{x}_{HI,\star})/\bar{x}_{HI,\star}$. In Eq. (4), $\bar{x}_{HI,\star}$, $A_{\star}$, and the slope $\alpha_A$ are free parameters that set the evolution of A as a function of $\bar{x}_{HI}$. The reionization scenarios considered in this study suggest $\bar{x}_{HI,\star} \le 0.5$. The first term of the right hand side in the above equation shows the $\bar{x}_{HI}$ dependence of $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ at $k = 0.05 \text{ h Mpc}^{-1}$ for stages when the ionization power spectrum ($\Delta_{\delta T_b, \gamma_{HI}}^2$) dominates the 21 cm signal power spectrum in comparison with the matter density. The second part of the right hand side of this equation shows the $\bar{x}_{HI}$ dependence during the initial stages of the EoR corresponding to $\bar{x}_{HI} > \bar{x}_{HI,\min}$ when the density power spectrum ($\Delta_{\delta\delta}^2$) and the anti-correlation between

the density and neutral fraction ($\Delta_{\delta T_b}^2$) are important (see, e.g., Fig. 2). We note that the second term rapidly decreases as $\bar{x}_{HI}$ decreases and is negligible at the stages when the peak of A vs $\bar{x}_{HI}$ curves occurs. Thus, we neglect the second term when we determine the values for $\beta_A$ in terms of $\alpha_A$ at the maximum. Deriving an analytical form for $\bar{x}_{HI,\min}$ is not straightforward, as we estimate $\bar{x}_{HI,\min}$ numerically for a set of our input parameters and thus $\bar{x}_{HI,\min}$ is a derived quantity for a given model.

Next, we check the accuracy of the fitting form of A as used in Eq. (4). We consider the evolution of $A = \Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ ($k = 0.05 \text{ h Mpc}^{-1}$) as a function of $\bar{x}_{HI}$ from different GRIZZLY reionization scenarios as inputs. We vary $A_{\star}$, $\bar{x}_{HI,\star}$, and $\alpha_A$ on regularly spaced grids respectively in the ranges [0, 20], [0, 1], and [0, 2], and determine their best-fit values. For a particular reionization scenario, we fit the evolution of A values corresponding to different stages of reionization together using Eq. (4). The top panel of Fig. 7 shows the best-fit values of the parameters $A_{\star}$, $\bar{x}_{HI,\star}$, and $\alpha_A$. Although we observe a clear correlation between $A_{\star}$ and $\bar{x}_{HI,\star}$, as A is sensitive to small changes in these parameters, we still keep all of them as independent parameters in our final ansatz. The bottom panel of Fig. 7 shows the evolution of $A - A_f$ as a function of $\bar{x}_{HI}$ for different GRIZZLY reionization scenarios. Here, $A_f$ represents the best-fit value of A obtained using Eq. (4).

Next, we consider the $\gamma$ parameter. The middle panel of Fig. 5 shows the dependencies of the best-fit values of $\gamma$ on $\bar{x}_{HI}$. As expected, $\gamma \to 0$ for $\bar{x}_{HI} \to 1$. For $\bar{x}_{HI} > \bar{x}_{HI,\min}$, $\gamma$ increases as $\bar{x}_{HI}$ decreases and roughly shows a power-law dependency on $\bar{x}_{HI}$. We modeled this part as $\gamma = 10^{-2} \bar{x}_{HI}^{\gamma_p}$ where $\gamma_p$ is the power-law index. For $\bar{x}_{HI} \le \bar{x}_{HI,\min}$, $\gamma$ decreases with $\bar{x}_{HI}$ and roughly shows an exponential drop while $\gamma$ reaches negative values for $\bar{x}_{HI} \to 0$. We modeled this drop of $\gamma$ with the decreases of $\bar{x}_{HI}$ as $\gamma = \gamma_1 \exp[\gamma_c \bar{x}_{HI}] + \gamma_0$. The parameter $\gamma_c$ controls the rate of decrease of $\gamma$ for $\bar{x}_{HI} \le \bar{x}_{HI,\min}$, while $\gamma_0 + \gamma_1$ represents $\gamma$ values when $\bar{x}_{HI} \to 0$. Therefore we choose the fitting form for $\gamma$ to be written as

$$\gamma = (\gamma_1 \exp[\gamma_c \bar{x}_{HI}] + \gamma_0) (1 - \mathcal{H}(\bar{x}_{HI,\min})) + 10^{-2} \bar{x}_{HI}^{\gamma_p} \mathcal{H}(\bar{x}_{HI,\min}), \quad (5)$$

where, $\mathcal{H}$ is the Heaviside step function. In order to reduce the number of parameters, we use a boundary condition for $\gamma$. As the two discontinuous functional forms of $\gamma$ for $\bar{x}_{HI} \le \bar{x}_{HI,\min}$ and $\bar{x}_{HI} > \bar{x}_{HI,\min}$ have the same value for $\bar{x}_{HI} = \bar{x}_{HI,\min}$, we can determine $\gamma_p$ as $\gamma_p = \log(\gamma_1 \exp[\gamma_c \bar{x}_{HI,\min}] + \gamma_0)/(0.01 \log \bar{x}_{HI,\min})$.

Next, we check the accuracy of the fitting form of $\gamma$ as defined in Eq. (5). We consider different evolutions of $\gamma$ as a function of $\bar{x}_{HI}$ from the middle panel of Fig. 5. We vary $\gamma_0$, $\log \gamma_1$, and $\gamma_c$ on regularly spaced grids respectively in the ranges [-1.5, 0], [-3, 0], and [0, 10] and determine their best-fit values. Note that the fitting considers different $\gamma$ values corresponding to different stages of a particular reionization history to obtain the best-fit values. The top panel of Fig. 8 shows the best-fit values of the parameters $\gamma_1$, $\gamma_c$, and $\gamma_0$. The bottom panel of Fig. 8 shows the deviation $\gamma - \gamma_f$ as a function of $\bar{x}_{HI}$ for different GRIZZLY reionization scenarios. Here, $\gamma_f$ represents the $\gamma$ value for the best-fit value of the parameters following Eq. (5). Similarly to the fit for A, here we find a prominent correlation between $\gamma_c$ and $\gamma_1$. As $\gamma_c$ and $\log \gamma_1$ behave roughly linearly, we consider $\gamma_1 = 10^{(1.3-\gamma_c)/(2.25)}$ which reduces the number of parameters to $\gamma_c$ and $\gamma_0$. Here, $\gamma_c$ accounts for the change in k-dependence of the bias $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ with $\bar{x}_{HI}$ for $\bar{x}_{HI} \le \bar{x}_{HI,\min}$. $\gamma_0$ accounts for the power-law dependence on k feature of $\Delta_{\delta T_b}^2/\Delta_{\delta\delta}^2$ in addition to small-scale feature $1/[1+(k/0.3)^{1.5}]$ at stages when $\bar{x}_{HI} \to 0$.
<figure><img src="image_7.png" /><figcaption>Fig. 7. Results of fitting the dependence of A on $\bar{x}_{HI}$ using Eq. (4). Top panel shows the best-fit values of the parameters $A_*$, $\bar{x}_{HI,*}$ and $\alpha_A$ when we fit different curves of the left panel of Fig. 5 using Eq. (4). Here, A represents the ratio $\Delta^2_{\delta T_b}/\Delta^2_{\delta\delta}$ at $k \sim 0.05 h \text{ Mpc}^{-1}$. The bottom panel shows the difference between the input/true A values and predicted (or fitted) $A_f$ values with the best-fit parameters as a function of average ionization fraction $\bar{x}_{HI}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. Results of fitting the dependence of $\gamma$ on $\bar{x}_{HI}$ using Eq. (5). Top panel shows the best-fit values of the parameters $\log \gamma_1$, $\gamma_c$, and $\gamma_0$ when we fit different curves in the middle panel of Fig. 5 using Eq. (5). The bottom panel shows the difference between the input or true $\gamma$ values and predicted/fitted $\gamma_f$ values with the best-fit parameters as a function of average ionization fraction $\bar{x}_{HI}$.</figcaption></figure>

Now, we are left with an ansatz which depends on five redshift-independent parameters $A_*$, $\bar{x}_{HI,*}$, $\alpha_A$, $\gamma_c$, and $\gamma_0$ to model the EoR 21 cm power spectrum as a function of k modes and the reionization history, which is parametrized by the globally averaged neutral fraction $\bar{x}_{HI}(z)$. In the following section, we introduce our modeling of reionization history.

## 2.5. Modeling the reionization history

The EoR 21 cm signal observations with radio interferometers are initially going to produce power spectrum $\Delta^2_{\delta T_b}(z, k)$ at different redshifts $z$. However, our ansatz (Eqs. (3)–(5)) predicts $\Delta^2_{\delta T_b}$ as a function of k and $\bar{x}_{HI}$. Therefore, it is necessary to model the reionization history $\bar{x}_{HI}$ as a function of redshift $z$.

The top panel of Fig. 9 shows the redshift evolution of the volume-averaged neutral fraction $\bar{x}_{HI}$ for all the GRIZZLY reionization scenarios considered in this work. We find $\bar{x}_{HI} \approx 1$ for $z \ge 11$ while the bulk of reionization occurs within a narrow redshift window ($\Delta z \le 4$) at $z \le 11$. We note that the reionization histories are asymmetric around $\bar{x}_{HI} = 0.5$. There is no well-established analytical form which accurately represents the redshift evolution of $\bar{x}_{HI}$. One possible analytical form to represent the asymmetric redshift evolution of $\bar{x}_{HI}$ is (e.g., Heinrich et al. 2017)

$$ \bar{x}_{HI}(z_0, \alpha_0, \Delta z, z) = \frac{1}{2} \left[ 1 - \tanh \left\{ \frac{y(z_0) - y(z)}{\Delta y} \right\} \right], \quad (6) $$ 
where $y(z) = (1 + z)^{\alpha_0}$, and $\Delta y = \alpha_0 (1 + z)^{\alpha_0 - 1} \times \Delta z$.

<figure><img src="image_9.png" /><figcaption>Fig. 9. Modeling of reionization history $\bar{x}_{HI}(z)$ using Eq. (6). The top panel shows the evolution of the neutral fraction as a function of redshift for different reionization scenarios generated using GRIZZLY code. The bold black curve corresponds to our GRIZZLY fiducial reionization scenario. The bottom panel shows fitting error $\Delta \bar{x}_{HI} = \bar{x}_{HI} - \bar{x}_{HI,f}$ on the redshift evolution of $\bar{x}_{HI}$. We have used an asymmetric tanh function (see Eq. (6)) to model the redshift evolution of $\bar{x}_{HI}$.</figcaption></figure>

Here $z_0$, $\Delta z$, and $\alpha_0$ are free parameters which govern the evolution of the mean neutral fraction $\bar{x}_{HI}$ with $z_0$ representing the redshift at which $\bar{x}_{HI} = 0.5$, $\Delta z$ the duration of the reionization, and $\alpha_0$ the asymmetry of reionization history around $z_0$. We note that
<table><caption>**Table 1.** List and description of the eight parameters used to model the redshift evolution of the EoR 21 cm brightness temperature power spectrum $\Delta_{\delta T_b}^2$.</caption><thead><tr><th>Parameters</th><th>Description</th></tr></thead><tbody><tr><td>$z_0$</td><td>Redshift corresponding to $\bar{x}_{HI} = 0.5$.</td></tr><tr><td>$\Delta z$</td><td>Redshift range of reionization in a tanh model.</td></tr><tr><td>$\alpha_0$</td><td>Asymmetry parameter around $\bar{x}_{HI} = 0.5$ in the redshift evolution of $\bar{x}_{HI}$.</td></tr><tr><td>$A_*$</td><td>Maximum value of the bias at $k = 0.05 \text{ h Mpc}^{-1}$.</td></tr><tr><td>$\bar{x}_{HI,*}$</td><td>Mean neutral fraction at the redshift when the bias at $k = 0.05 \text{ h Mpc}^{-1}$ gets the maxima.</td></tr><tr><td>$\alpha_A$</td><td>Power-law index on $\bar{x}_{HI}$ which accounts for the change of bias as a function of $\bar{x}_{HI}$ at $k = 0.05 \text{ h Mpc}^{-1}$.</td></tr><tr><td>$\gamma_c$</td><td>Account for the change in scale-dependence of bias with $\bar{x}_{HI}$.</td></tr><tr><td>$\gamma_0$</td><td>Account for the all-scale feature of bias in addition to small-scale feature $1/[1 + (k/0.3)^{1.5}]$ at stages with $\bar{x}_{HI} \to 0$.</td></tr></tbody></table>

**Notes.** Here bias represents the ratio of $\delta T_b$ to density power spectra $\Delta_{\delta T_b}^2 / \Delta_{\delta \delta}^2$.

the parametrization of $\bar{x}_{HI}$ as a function of redshift is not unique. Here, we have used an analytically simpler form for $\bar{x}_{HI}(z)$. An alternative parametrization of the reionization history such as the one used in Trac (2018) could also produce a good fit to the reionization scenarios used in this study.

The bottom panel of Fig. 9 shows the fitting error $\Delta \bar{x}_{HI} = (\bar{x}_{HI} - \bar{x}_{HI,f})$ between the simulated ($\bar{x}_{HI}$) and best-fit ($\bar{x}_{HI,f}$) values of the reionization histories. The best-fit ionization history for a given simulated reionization scenario is obtained by varying $z_0$, $\Delta z$, and $\alpha_0$ on uniformly spaced grids and minimizing the mean square error

$$E_{\text{tot}}(z_0, \Delta z, \alpha_0) = \frac{1}{n_z} \sum_{z < 15.5} [\bar{x}_{\text{HI}}(z_0, \alpha_0, \Delta z, z) - \bar{x}_{\text{HI}}(z)]^2.$$

Here $n_z$ is the number of redshifts considered for a reionization history, which can differ for different EoR models of GRIZZLY. We vary $z_0$, $\Delta z$, and $\alpha_0$ in the range [5, 15], [0, 3], and [0, 10] ranges, respectively. All the GRIZZLY reionization scenarios show a good fit with $\Delta \bar{x}_{HI} \lesssim 0.03$ for the majority of the reionization redshifts ranges. At the same time, the error increases up to ~0.08 near the end of EoR when $\bar{x}_{HI} \to 0$. This shows that a tanh form has some trouble capturing the fast drop in $\bar{x}_{HI}$ during the tail end of the reionization.

In Sect. 2.4, we used five redshift-independent parameters to model the scale dependence of $\Delta_{\delta T_b}^2$ as a function of $\bar{x}_{HI}$. In this section, we used three $z$ independent parameters to model the redshift evolution of $\bar{x}_{HI}$. Thus in the end, we are left with a set of eight redshift-independent free parameters $\theta = [z_0, \Delta z, \alpha_0, A_*, \bar{x}_{HI,*}, \alpha_A, \gamma_c, \gamma_0]$ to model both the scale dependence and redshift evolution of the EoR 21 cm power spectrum together with the reionization history. See Table 1 for a description of the parameters.

## 3. Results

We next apply our eight-parameter ansatz of the EoR 21 cm power spectra to various reionization scenarios obtained from different simulations. As inputs, we consider a C2RAY, a 21CMFAST, and all the GRIZZLY simulations as mentioned in the previous section. These three 21 cm simulation frameworks are based on very different source models, but to be consistent with our ansatz for $\Delta_{\delta T_b}^2$, we assume $T_S \gg T_\gamma$ in all of them. We refer the readers to Ghara et al. (2015a), Mellema et al. (2006), and Mesinger et al. (2011) for the details of the algorithms used in these three simulations.

From each reionization scenario, we extract power spectra for eight different redshifts covering the majority of the

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Three different input reionization histories and power spectra which are used to check the performance of our ansatz. The top, middle and bottom panels correspond to C2RAY, 21CMFAST, and GRIZZLY fiducial input reionization scenarios. The left panels show the redshift evolution of $\bar{x}_{HI}$ while the right panels show the set of input power spectra used in the MCMC analysis.</figcaption></figure>

EoR. The redshift ranges for C2RAY, 21CMFAST, and GRIZZLY scenarios are [6.1, 9], [6.4, 9.6], and [7.1, 11.1] respectively. Figure 10 shows the redshift evolution of $\bar{x}_{HI}$ (left column) and
the corresponding simulated EoR 21 cm signal power spectra (right column) for the C2RAY (top row), 21CMFAST (middle row), and the fiducial GRIZZLY scenario (bottom row). The redshifts of these power spectra span ~60 MHz of observational bandwidth while the frequency difference between two adjacent redshifts is Δν ~ 7 MHz. The latter is smaller than the typical bandwidth used in EoR 21 cm data analysis. For example, Mertens et al. (2020) used 12 MHz bandwidth to estimate the upper limits on the 21 cm power spectrum at redshift 9.1. The main motivation to use a smaller bandwidth is to resolve the fast evolving nature of the large-scale power spectrum around the peak. However, the choice of a smaller bandwidth will need more observation hours to reach the same signal-to-noise ratio.

Next we use these power spectra to constrain the eight parameters of our ansatz using Markov chain Monte Carlo (MCMC) based parameter estimation framework. We use the publicly available code COSMOMC<fnref n="11" /> (Lewis & Bridle 2002) for exploring the log-likelihood of these eight parameters θ. The log-likelihood in our MCMC algorithm is estimated as,

$$ \chi^2(\theta) = - \sum_{i,j} \left( \frac{\Delta^2_{\delta T_b, m}(z_i, k_j, \theta) - \Delta^2_{\delta T_b, o}(z_i, k_j)}{\Delta^2_{o, \text{err}}} \right)^2, $$

where Δ²_{δT_b,m} and Δ²_{δT_b,o} are the modeled and the simulated input power spectra, respectively. The index i runs over the eight input redshifts while the index j runs over k-bins with 0.05 h Mpc⁻¹ ≤ k ≤ 0.6 h Mpc⁻¹. Each MCMC analysis here is done with 8 independent walkers (sequences of parameter values in MCMC), each of which takes 10⁶ steps.

The quantity Δ²_{o, err} in the denominator is the error used in our MCMC analysis. In principle, this error should include measurement error, sample variance, and the imperfection of the ansatz. However, as the aim is to show a proof of concept of our ansatz, we consider the simple case in which the error Δ²_{o, err} = 1 mK² is k and z independent. In general, Δ²_{o, err} is expected to increase towards higher redshifts. Considering any observation, the scale dependence of Δ²_{o, err} changes with redshift as the uv coverage of an interferometric observation changes with observation frequency and also the sky noise varies with the frequency. Here, we do not consider such realistic situations and will address these issues in a follow-up work.

We run the MCMC analysis on the C2RAY, 21CMFAST, and GRIZZLY scenarios. The input power spectra for each scenario are presented in Fig. 10. The outcomes of the MCMC analysis are presented in the following subsections.

### 3.1. Scenario I: C2RAY

First, we consider a reionization scenario which is generated using the EoR 21 cm signal modeling code C2RAY (Mellema et al. 2006). This code uses the gridded density field and halo lists from N-body simulations and applies “Conservative Causal Ray-tracing method” based 3D radiative transfer to produce ionization fraction fields at different stages of reionization. The set of δT_b power spectra and the redshift evolution of $\bar{x}_{HI}$ of the input reionization history as obtained from a C2RAY simulation are shown in the top row of Fig. 10. This scenario uses the same dark-matter fields and halo list as used in the GRIZZLY simulations. C2RAY also considers contributions from all dark-matter halos with their masses larger than 10⁹ M_☉ and assumes the rate of production of the ionizing photons to be

1.3 × 10⁴² × M_halo/M_☉ s⁻¹. In this simulation, the mean-free-path length of the ionizing photons is chosen as 70 Mpc. In this C2RAY reionization scenario, $\bar{x}_{HI}$ decreases from 0.9 to 0 as reionization progressed between z ≈ 8 and 6. We find that the evolution of the ratio Δ²_{δT_b}/Δ²_{δδ} at k = 0.05 h Mpc⁻¹ reaches a maximum value of ≈5 at z ≈ 6.4 when $\bar{x}_{HI}$ ≈ 0.3.

Figure 11 shows the posteriors of the eight parameters θ of our power spectrum ansatz when we use the set of C2RAY power spectra as inputs to our MCMC framework. The off-diagonal panels show the joint probability distribution for a pair of parameters where 2D contours represent 1σ and 2σ confidence levels respectively. The curves in the diagonal panels represent the marginalized probability distributions of the individual ansatz parameters. The plot shows that while most of the parameters are well-constrained, some of them are not. One reason behind this might be the degeneracy of these parameters with the other parameters. The best-fit parameter values obtained from this analysis are z₀ = 6.9, Δz = 0.83, α₀ = 6.4, A_★ = 5.2, $\bar{x}_{HI,★}$ = 0.27, α_A = 1.8, γ_c = 2, and γ₀ = -1.2 (see also Table 2). The best-fit values of A_★ and $\bar{x}_{HI,★}$ agree well with the input reionization scenario which corresponds to A_★ = 5.01 and $\bar{x}_{HI,★}$ = 0.22. The comparison between input and best-predicted models for this reionization scenario is shown in the left column of Fig. 12. Here The top-left panel shows the comparison as a function of k at different redshifts while the middle-left panel shows the redshift evolution of the predicted and input ratio Δ²_{δT_b}/Δ²_{δδ} for different k-bins. The curves indicate that our ansatz performs well at different stages of the reionization and for the k range we considered here. A comparison between the input reionization history and the ansatz predictions (bottom-left panel of Fig. 12) suggests an excellent recovery of the redshift evolution of $\bar{x}_{HI}$ using our model. Figure 13 shows the corresponding deviation $\bar{x}_{HI} - \bar{x}_{HI,best-fit}$ where $\bar{x}_{HI,best-fit}$ represents the prediction using the MCMC best-fit parameters. Here, we find that the deviation remains between ±0.05 as indicated by the blue thick curve.

### 3.2. Scenario II: 21CMFAST

Our second input reionization scenario is generated using the publicly available semi-numerical 21 cm code 21CMFAST (Mesinger et al. 2011; Park et al. 2019). In this semi-numerical approach, the density fields are generated following the first-order perturbation theory (Zel'dovich 1970) while the ionization fields are produced using the excursion-set approach (Furlanetto et al. 2004). We assume the following parameters: fraction of galactic gas in stars for 10¹⁰ M_☉ halo $f_{★,10}$ = 0.05, the power-law index for star formation and halo mass relation α_★ = 0.5, the UV ionizing escape fraction for 10¹⁰ M_☉ halo $f_{esc,10}$ = 0.1, the power-law index for UV escape fraction and halo mass relation α_{esc} = -0.5, the characteristic mass scale for star formation suppression $M_{turn}$ = 5 × 10⁸ M_☉, star-formation timescale in units of the Hubble time t_★ = 0.5, number of ionizing photons per stellar baryon N_γ = 5000, and mean-free path of ionizing photons $R_{mfp}$ = 50 Mpc (for the details of these parameters, see Park et al. 2019). In this case, reionization ends at z ~ 6 while the majority of the ionization happens below z ≲ 10. The corresponding input reionization history and the input power spectra are shown in the middle row of Fig. 10. The input set of the 21 cm power spectra of the reionization scenario shows features similar to those of the fiducial GRIZZLY and C2RAY scenarios. We find that the ratio Δ²_{δT_b}/Δ²_{δδ} at k = 0.05 h Mpc⁻¹ reaches a maximum value of ≈1.3 at z ≈ 7.2 corresponding to $\bar{x}_{HI}$ = 0.35.

<footnote n="11">¹¹ https://cosmologist.info/cosmomc/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Posterior of the ansatz parameters (see Table 1). This shows the constraints of the fiducial EoR scenario obtained using the MCMC analysis in terms of the EoR power spectrum model parameters. The contours in the two-dimensional contour plots represent $1\sigma$ and $2\sigma$ confidence levels respectively. The curves in the diagonal panels are the marginalized probability distributions of the eight parameters. The MCMC analysis is performed on inputs from C2RAY.</figcaption></figure>

We repeat the MCMC analysis (as done for C2RAY) to obtain the posterior of our ansatz parameters. The middle row of Table 2 shows the posterior constraints on the eight parameters of our ansatz. The best-fit parameter values are $z_0 = 7.7$, $\Delta z = 1.7$, $\alpha_0 = 5.5$, $A_\star = 1.6$, $\bar{x}_{HI,\star} = 0.38$, $\alpha_A = 1.9$, $\gamma_c = 3.5$, and $\gamma_0 = -0.8$. The corresponding power spectra predictions at all the redshifts agree well with the input simulated power spectra (see the top and middle panels of the central column of Fig. 12). Similar to the previous case, the best-fit values of $A_\star$ and $\bar{x}_{HI,\star}$ agree well with the input 21CMFAST reionization scenario which has $A_\star = 1.3$ and $\bar{x}_{HI,\star} = 0.34$. We compare the simulated reionization history $\bar{x}_{HI}(z)$ with that predicted from our MCMC analysis in the bottom-central panel of Fig. 12. Similar to the C2RAY scenario, our framework works efficiently in this case as well and recovers reionization history $\bar{x}_{HI}$ within a

maximum deviation between $\pm 0.03$ as shown by the orange line in Fig. 13.

### 3.3. Scenario III: GRIZZLY

The simulated input power spectra and the reionization history for the fiducial GRIZZLY reionization scenario are shown in the bottom panels of Fig. 10. As described in Sect. 2.1, this corresponds to GRIZZLY parameters $\zeta = 1$, $M_{\min} = 10^9 M_{\odot}$, and $\alpha_s = 1$ (for details, see Ghara et al. 2020). This results in a reionization history where the IGM gets 10% ionized at $z \approx 11$ while the reionization ends around $z \approx 7$ (see top right panel of Fig. 12). Here, the ratio $\Delta^2_{\delta T_0}/\Delta^2_{\delta\delta}$ at $k = 0.05 \, h \, \text{Mpc}^{-1}$ reaches a maximum value of $\approx 6$ at $z \approx 7.3$ which corresponds to $\bar{x}_{HI} = 0.25$.
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Comparison between the input reionization scenario and recovered model from the MCMC analysis using EoR power spectra as input to the phenomenological model considered in this work. Left to right columns represent C2RAY, 21CMFAST, and GRIZZLY simulation respectively. From top to bottom we show, respectively, the ratio of 21 cm signal to density power spectrum ($\Delta_{\delta T_b}^2/\Delta_{\delta \delta}^2$) as a function of k-scales at different redshifts, the redshift evolution of $\Delta_{\delta T_b}^2/\Delta_{\delta \delta}^2$ for different scales and the corresponding reionization histories. The solid curves represent the inputs from simulations while the dotted curves stand for the MCMC best-fit prediction on the power spectrum model used in this work.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Redshift evolution of the difference between the input model neutral fraction and best-fit neutral fraction. The MCMC analysis here considers our eight-parameter model for the EoR power spectra. Different thin grey curves stand for different reionization scenarios from GRIZZLY simulation. The black curve represents our fiducial GRIZZLY simulation while the blue and orange curves present C2RAY and 21CMFAST input scenarios. The best-fit parameter values are obtained using input power spectra at eight redshifts for each reionization scenario.</figcaption></figure>

The posterior constraints on the ansatz parameters, obtained using MCMC analysis, are summarized in the bottom row of Table 2. The best-fit parameter values are $z_0 = 8.1$, $\Delta z = 1.2$, $\alpha_0 = 6.6$, $A_\star = 5.2$, $\bar{x}_{HI,\star} = 0.24$, $\alpha_A = 1.3$, $\gamma_c = 2.4$, and

$\gamma_0 = -0.88$. Similar to the previous cases, the predicted power spectra and the best-fit values of $A_\star$ and $\bar{x}_{HI,\star}$ are in good agreement with the simulated input values. The $A_\star$ and $\bar{x}_{HI,\star}$ values from the input reionization scenario are 5.8 and 0.22, respectively. The predicted evolution of $\bar{x}_{HI}$ as shown in the top right panel of Fig. 12 shows good agreement within the absolute deviations of 0.05.

We repeat our MCMC analysis for the other 23 GRIZZLY models as represented by the shaded grey lines in Figs. 9 and 3 considering the input power spectra at the same redshifts used for the fiducial GRIZZLY model. We plot the difference between the input and predicted neutral fraction for all these models in Fig. 13. These suggest that our ansatz represents the EoR 21 cm power spectrum close enough and can recover the reionization history $\bar{x}_{HI}(z)$ within an absolute error of ~0.1 (see the grey lines in Fig. 13). We note that the evolution of $\bar{x}_{HI}$ as a function of redshift is very steep around $z_0$ and thus even a small error in the estimate of $z_0$ will result in a large difference between the input (true) and predicted $\bar{x}_{HI}$. This is evident from the plots for various GRIZZLY models in Fig. 13 between redshift 7 and 8. The deviations could become larger for the models having steeper slopes around $z_0$.

# 4. Summary and conclusions

The redshifted 21 cm signal from the IGM during the EoR encodes unique information about that period. The 21 cm
<table><caption>**Table 2.** Summary of the outputs from the MCMC analysis and constraints on the eight parameter models of the EoR power spectra.</caption><thead><tr><th>Scenario</th><th>Parameters</th><th>Explored range</th><th>Mean</th><th>Standard deviation</th><th>Best fit</th><th>68% limits</th><th>95% limits</th></tr></thead><tbody><tr><td rowspan="8">C2RAY</td><td>z<sub>0</sub></td><td>[5, 15]</td><td>7.30</td><td>0.32</td><td>6.88</td><td>[6.9, 7.5]</td><td>[6.7, 7.98]</td></tr><tr><td>Δz</td><td>[0, 3]</td><td>1.24</td><td>0.56</td><td>0.83</td><td>[0.56, 1.40]</td><td>[0.42, 2.50]</td></tr><tr><td>α<sub>0</sub></td><td>[0, 10]</td><td>6.40</td><td>2.49</td><td>6.43</td><td>[5.03, 9.56]</td><td>[1.72, 10]</td></tr><tr><td>A<sub>★</sub></td><td>[0, 20]</td><td>4.18</td><td>1.61</td><td>5.17</td><td>[2.12, 5.18]</td><td>[1.22, 7.63]</td></tr><tr><td>x&#773;<sub>HI,★</sub></td><td>[0, 1]</td><td>0.23</td><td>0.07</td><td>0.27</td><td>[0.16, 0.30]</td><td>[0.09, 0.35]</td></tr><tr><td>α<sub>A</sub></td><td>[0.1, 2]</td><td>1.33</td><td>0.46</td><td>1.81</td><td>[1.07, 1.97]</td><td>[0.51, 2.0]</td></tr><tr><td>γ<sub>c</sub></td><td>[0, 10]</td><td>4.65</td><td>2.39</td><td>2.05</td><td>[1.4, 6.00]</td><td>[1.07, 9.07]</td></tr><tr><td>γ<sub>0</sub></td><td>[-3, 2]</td><td>-0.52</td><td>0.38</td><td>-1.20</td><td>[-0.66, -0.01]</td><td>[-1.31, 0.00]</td></tr><tr><td rowspan="8">21CMFAST</td><td>z<sub>0</sub></td><td>[5, 15]</td><td>8.16</td><td>0.42</td><td>7.70</td><td>[7.53, 8.70]</td><td>[6.92, 9.46]</td></tr><tr><td>Δz</td><td>[0, 3]</td><td>1.95</td><td>0.57</td><td>1.70</td><td>[1.31, 2.61]</td><td>[0.95, 2.94]</td></tr><tr><td>α<sub>0</sub></td><td>[0, 10]</td><td>5.91</td><td>2.47</td><td>5.50</td><td>[4.01, 9.21]</td><td>[1.56, 9.95]</td></tr><tr><td>A<sub>★</sub></td><td>[0, 20]</td><td>1.39</td><td>0.41</td><td>1.60</td><td>[0.86, 1.71]</td><td>[0.59, 2.30]</td></tr><tr><td>x&#773;<sub>HI,★</sub></td><td>[0, 1]</td><td>0.31</td><td>0.08</td><td>0.38</td><td>[0.23, 0.41]</td><td>[0.13, 0.47]</td></tr><tr><td>α<sub>A</sub></td><td>[0.1, 2]</td><td>1.27</td><td>0.46</td><td>1.90</td><td>[0.95, 1.95]</td><td>[0.43, 2.0]</td></tr><tr><td>γ<sub>c</sub></td><td>[0, 10]</td><td>3.78</td><td>1.80</td><td>3.50</td><td>[1.55, 4.62]</td><td>[0.94, 7.51]</td></tr><tr><td>γ<sub>0</sub></td><td>[-3, 2]</td><td>-0.68</td><td>0.43</td><td>-0.80</td><td>[-0.96, 0.05]</td><td>[-1.40, 0.0]</td></tr><tr><td rowspan="8">GRIZZLY</td><td>z<sub>0</sub></td><td>[5, 15]</td><td>8.03</td><td>0.21</td><td>8.10</td><td>[7.8, 8.2]</td><td>[7.60, 8.44]</td></tr><tr><td>Δz</td><td>[0, 3]</td><td>1.25</td><td>0.20</td><td>1.20</td><td>[1.11, 1.40]</td><td>[0.96, 1.56]</td></tr><tr><td>α<sub>0</sub></td><td>[0, 10]</td><td>5.61</td><td>2.30</td><td>6.60</td><td>[5.1, 9.8]</td><td>[2.0, 10]</td></tr><tr><td>A<sub>★</sub></td><td>[0, 20]</td><td>5.21</td><td>0.71</td><td>5.20</td><td>[4.41, 5.80]</td><td>[3.83, 6.61]</td></tr><tr><td>x&#773;<sub>HI,★</sub></td><td>[0, 1]</td><td>0.26</td><td>0.06</td><td>0.24</td><td>[0.22, 0.32]</td><td>[0.15, 0.40]</td></tr><tr><td>α<sub>A</sub></td><td>[0.1, 2]</td><td>1.46</td><td>0.38</td><td>1.30</td><td>[1.20, 1.94]</td><td>[0.80, 2.0]</td></tr><tr><td>γ<sub>c</sub></td><td>[0, 10]</td><td>2.41</td><td>0.90</td><td>2.40</td><td>[1.60, 3.20]</td><td>[0.76, 3.79]</td></tr><tr><td>γ<sub>0</sub></td><td>[-3, 2]</td><td>-1.05</td><td>0.51</td><td>-0.88</td><td>[-1.20, -0.50]</td><td>[-2.50, -0.25]</td></tr></tbody></table>

**Notes.** Top-to-bottom panels stand for the C2RAY, 21CMFAST, and GRIZZLY fiducial input scenarios. We use eight input power spectra for each of the reionization scenarios. The left to right columns represent the scenarios considered in this work, set of parameters, their ranges used in the MCMC analysis, mean parameter value, standard deviation, best-fit value, 68% and 95% credible limits obtained from the MCMC chains. The MCMC analysis is done using 8 chains each with $10^6$ steps.

observations indirectly probe the properties of the ionizing and heating sources and directly probes the ionization and thermal states of the IGM during the first billion years of our Universe. In this study, we focus on inferring the properties of the EoR IGM, rather than those of astrophysical sources, through 21 cm signal observations. The large-scale amplitude and scale-dependent features of the EoR 21 cm brightness temperature power spectrum depend on the ionization fraction of hydrogen and the morphology and distribution of the ionized or neutral regions. Our main aim is to develop a source parameter-free phenomenological model that constrains the properties of the EoR IGM using multi-redshift 21 cm power spectrum measurements. The framework constrains the redshift evolution of the average neutral fraction and a set of quantities related to the morphology and distribution of the ionized regions.

Using different GRIZZLY simulations, we study the scale-dependent features of the 21 cm power spectra at different stages of EoR. The quantity we aim to model is the ratio of the 21 cm brightness temperature ($\delta T_b$) to the density power spectrum also known as the 21 cm bias in the literature. We modeled this ratio as $\Delta^2_{\delta T_b}/\Delta^2_{\delta \delta} = A \frac{(\frac{\delta}{\delta \delta})^\gamma}{1+\left(\frac{\kappa}{\delta^3}\right)^{1.5}}$. Here, A represents bias at $k = 0.05 h \text{ Mpc}^{-1}$. We tested the goodness of fit of our ansatz at various stages of reionization using 24 different GRIZZLY scenarios. These tests suggest that the aforementioned functional form of the ratio of $\delta T_b$ to density power spectra efficiently reproduce

the EoR 21 cm power spectra for different reionization histories, accurately within $\le 10\%$ error (see Fig. 6).

As the A and $\gamma$ parameters in the above-mentioned ansatz evolve during reionization, we additionally model how these parameters evolve as a function of $\bar{x}_{HI}$. The model for A (Eq. (4)) uses three parameters, the maximum value of the ratio ($A_★$), the corresponding neutral fraction $\bar{x}_{HI,★}$, and a power-law index $\alpha_A$. The evolution of $\gamma$ can be described with two parameters (Eq. (5)): $\gamma_c$ which accounts for the change in scale-dependence, and $\gamma_0$ which accounts for the deviation of the scale dependence of $\Delta^2_{\delta T_b}/\Delta^2_{\delta \delta}$ from $1/[1+(k/0.3)^{1.5}]$ at small-scales (see Sect. 2.4 for details).

Using the GRIZZLY simulations, we fit the evolution of $\bar{x}_{HI}$ as a function of redshift using three parameters (see Eq. (6)). These are: redshift $z_0$ which corresponds to $\bar{x}_{HI} = 0.5$, redshift range of reionization $\Delta z$ in a tanh reionization model, and asymmetry parameter $\alpha_0$ to invoke asymmetry in history around $\bar{x}_{HI} = 0.5$. We tested the goodness of this form of $\bar{x}_{HI}(z)$ using 24 different GRIZZLY reionization models and found them to be consistent within $\Delta \bar{x}_{HI} = \pm 0.1$ (see Fig. 9).

In the end, we are left with a set of eight redshift and scale-independent parameters $\theta = [z_0, \Delta z, \alpha_0, A_★, \bar{x}_{HI,★}, \alpha_A, \gamma_c, \gamma_0]$ to jointly model the redshift evolution and scale-dependence of the ratio $\Delta^2_{\delta T_b}/\Delta^2_{\delta \delta}$. We demonstrate the performance of this ansatz with 24 GRIZZLY models, one reionization scenario from C2RAY and one from 21CMFAST. We use as an input the power spectra
simulated at eight redshifts within the interval $z = [7.1-11.1]$ (corresponds to ~60 MHz bandwidth) and perform a Bayesian MCMC analysis to constrain the ansatz parameters for each of the three reionization scenarios. All these tests collectively indicate that our ansatz reproduces the scale-dependence and the redshift evolution of the ratio $\Delta^2_{\delta T_b}/\Delta^2_{\delta \delta}$ reasonably well for a variety of reionization models considered here. The predicted redshift evolution of $\bar{x}_{HI}$, using the best-fit MCMC parameter set $\theta$, matches nicely with the input reionization history within $\Delta\bar{x}_{HI} = \pm 0.1$ (see Fig. 13). At the same time the constrained values for $A_*$, $\bar{x}_{HI,*}$ match closely with the reionization scenarios (see Fig. 12 and Table 2).

Our aforementioned approach is similar in spirit to a few previous studies. For example, Battaglia et al. (2013) simulate a reionization redshift field for a given density field (filtered at a particular scale). However, their approach fundamentally assumes a strong correlation between the density and reionization redshift fields. On the other hand, McQuinn & D’Aloisio (2018) take a perturbative approach to construct an EoR 21 cm signal field using the underlying density field. Their formalism assumes some source field and patchiness-dependent bias factors as well as a characteristic size of ionized regions to connect the 21 cm signal field with the density field. Unlike them, our approach is conceptually simpler and directly connects the EoR 21 cm signal power spectrum with the matter power spectrum. This phenomenological model directly exploits the features of multi-redshift EoR 21 cm signal power spectra for predicting the quantities related to the EoR IGM states. Thus, it makes our ansatz more flexible, computationally inexpensive and faster while exploring IGM parameters in the context of current and future observations. Additionally, our phenomenological ansatz is agnostic to the various methods of EoR simulations (3D and 1D radiative transfer and excursion-set based), whereas the analysis in Mirocha et al. (2022) is restricted to the excursion-set based simulations only.

The analysis presented in this paper is only based on the “inside-out” reionization scenario where the highly dense regions around the radiating sources become ionized first. We speculate that the same is also applicable for an “outside-in” reionization model as the scale-dependence of the EoR power spectra as well as the redshift evolution of the large-scale power spectrum are qualitatively similar to the “inside-out” case for the wavenumber range of our interest (see e.g., Figs. 2 and 3 of Pagano & Liu 2020). However, in the “outside-in” case, we do not expect a trough of the large-scale power spectra at $\bar{x}_{HI,min}$ and thus we expect an even simpler form of the ansatz for the EoR power spectra. We leave a detailed investigation for the “outside-in” reionization scenario for a future study.

The accuracy of the ansatz predictions of the EoR 21 cm power spectra suggests that it will be useful to constrain reionization scenarios using existing and upcoming measurements from LOFAR, MWA, HERA, and SKA. Here, we have used a simple-minded constant error on the input power spectra. It will be interesting to see the performance of this model when realistic errors are taken into account. We plan to address this in our future work.

Our model is based on $\delta T_b$ power spectra with high spin temperatures. This model also works if the gas temperature of the neutral IGM remains uniform. However, things get complicated when the presence of spin temperature fluctuations modifies the power spectra significantly. One expects a more complex evolution of the power spectrum in that case. Modeling such behaviors is out of the scope of this paper and will be addressed in a follow-up work.

Here we consider only the EoR 21 cm power spectra measurements to constrain the reionization history. In addition, information from several other probes such as the Thomson scattering optical depth measurement from CMB observation (e.g., Planck Collaboration VI 2020), the Gunn-Peterson trough in high-$z$ quasar spectra (e.g., Fan et al. 2006; Becker et al. 2015), observations of high-$z$ Ly-$\alpha$ emitters (e.g., Hu et al. 2010; Morales et al. 2021), and the Ly-$\alpha$ damping wings in high-$z$ quasar spectra (e.g., Bañados et al. 2018) can also be combined with the HI measurements from the EoR to tighten the constraint on the reionization history. We plan to address this in a future study using POLAR (Ma et al. 2023) algorithm which is based on GRIZZLY and the semi-analytical galaxy formation code L-Galaxies 2020 and self-consistently model the evolution of galaxy properties during the EoR.

Although our framework is efficient in recovering the redshift evolution of the average ionization fraction using the measurements of the EoR 21 cm signal power spectrum, there are several aspects that need improvements. While most of the parameters of this framework are easy to interpret, understanding the physical meaning of a few parameters such as $\gamma_c$ and $\gamma_0$ is non-trivial. These parameters are connected to the morphology and distribution of the ionized regions. It is important to understand how these quantities are linked to the distribution of morphological quantities such as volume, surface, and mean curvature (see, e.g., Giri et al. 2018; Giri & Mellema 2021; Ghara et al. 2024). This study is also beyond the scope of this work. We plan to address it in a future study.

**Acknowledgements.** R.G., A.K.S., and S.Z. acknowledge support from the Israel Science Foundation (grant no. 255/18). Furthermore, R.G. acknowledges support from the Kaufman Foundation (Gift no. GF01364). S.Z. also acknowledges Alexander von Humboldt Foundation for the Humboldt Research award. A.K.S. also acknowledges support from National Science Foundation (grant no. 2206602). G.M. acknowledges support by the Swedish Research Council grant 2020-04691. L.V.E.K. acknowledges the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 884760, “CoDEX”).

## References

* Abdurashidova, T. H. C. Z., Adams, T., Aguirre, J. E., et al. 2023, ApJ, 945, 124

* Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022a, ApJ, 924, 51

* Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022b, ApJ, 925, 221

* Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473

* Barkana, R. 2018, Nature, 555, 71

* Battaglia, N., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 81

* Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402

* Becker, G. D., Davies, F. B., Furlanetto, S. R., et al. 2018, ApJ, 863, 92

* Bera, A., Ghara, R., Chatterjee, A., Datta, K. K., & Samui, S. 2023, JApA, 44, 10

* Bonaldi, A., & Brown, M. L. 2015, MNRAS, 447, 1973

* Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67

* Chapman, E., Zaroubi, S., Abdalla, F. B., et al. 2016, MNRAS, 458, 2928

* Chatterjee, A., Dayal, P., Choudhury, T. R., & Hutter, A. 2019, MNRAS, 487, 3560

* Chatterjee, A., Dayal, P., Choudhury, T. R., & Schneider, R. 2020, MNRAS, 496, 1445

* Cohen, A., Fialkov, A., Barkana, R., & Monsalve, R. A. 2020, MNRAS, 495, 4845

* Datta, K. K., Bharadwaj, S., & Choudhury, T. R. 2007, MNRAS, 382, 809

* Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526

* de Lera Acedo, E., de Villiers, D. I. L., Razavi-Ghods, N., et al. 2022, Nat. Astron., 6, 984

* DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001

* Eide, M. B., Ciardi, B., Graziani, L., et al. 2020, MNRAS, 498, 6083

* Eilers, A.-C., Davies, F. B., & Hennawi, J. F. 2018, ApJ, 864, 53

* Ewall-Wice, A., Hewitt, J., Mesinger, A., et al. 2016, MNRAS, 458, 2710

* Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117
Fialkov, A., Barkana, R., & Cohen, A. 2018, Phys. Rev. Lett., 121, 011101

Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1

Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181

Gan, H., Koopmans, L. V. E., Mertens, F. G., et al. 2022, A&A, 663, A9

Gan, H., Mertens, F. G., Koopmans, L. V. E., et al. 2023, A&A, 669, A20

Gehlot, B. K., Koopmans, L. V. E., Offringa, A. R., et al. 2022, A&A, 662, A97

Georgiev, I., Mellema, G., Giri, S. K., & Mondal, R. 2022, MNRAS, 513, 5109

Ghara, R., & Mellema, G. 2020, MNRAS, 492, 634

Ghara, R., Choudhury, T. R., & Datta, K. K. 2015a, MNRAS, 447, 1806

Ghara, R., Datta, K. K., & Choudhury, T. R. 2015b, MNRAS, 453, 3143

Ghara, R., Choudhury, T. R., & Datta, K. K. 2016, MNRAS, 460, 827

Ghara, R., Choudhury, T. R., Datta, K. K., & Choudhuri, S. 2017, MNRAS, 464, 2234

Ghara, R., Giri, S. K., Mellema, G., et al. 2020, MNRAS, 493, 4728

Ghara, R., Giri, S. K., Ciardi, B., Mellema, G., & Zaroubi, S. 2021, MNRAS, 503, 4551

Ghara, R., Mellema, G., & Zaroubi, S. 2022, JCAP, 2022, 055

Ghara, R., Bag, S., Zaroubi, S., & Majumdar, S. 2024, MNRAS, 530, 191

Ghosh, A., Prasad, J., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2012, MNRAS, 426, 3295

Giri, S. K., & Mellema, G. 2021, MNRAS, 505, 1863

Giri, S. K., Mellema, G., & Ghara, R. 2018, MNRAS, 479, 5596

Giri, S. K., D'Aloisio, A., Mellema, G., et al. 2019, JCAP, 2019, 058

Greig, B., Mesinger, A., Koopmans, L. V. E., et al. 2021, MNRAS, 501, 1

Harker, G., Zaroubi, S., Bernardi, G., et al. 2009, MNRAS, 397, 1138

Heinrich, C. H., Miranda, V., & Hu, W. 2017, Phys. Rev. D, 95, 023513

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hothi, I., Chapman, E., Pritchard, J. R., et al. 2021, MNRAS, 500, 2264

Hu, E. M., Cowie, L. L., Barger, A. J., et al. 2010, ApJ, 725, 394

Islam, N., Ghara, R., Paul, B., Choudhury, T. R., & Nath, B. B. 2019, MNRAS, 487, 2785

Kamran, M., Ghara, R., Majumdar, S., et al. 2021, MNRAS, 502, 3800

Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2019, ApJ, 884, 105

Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2020, ApJ, 888, 70

Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133

Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

Lidz, A., Zahn, O., McQuinn, M., et al. 2007, ApJ, 659, 865

Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023019

Ma, Q.-B., Ghara, R., Ciardi, B., et al. 2023, MNRAS, 522, 3284

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429

McGreer, I. D., Mesinger, A., & D'Odorico, V. 2015, MNRAS, 447, 499

McQuinn, M., & D'Aloisio, A. 2018, JCAP, 2018, 016

Mellema, G., Iliev, I. T., Pen, U.-L., & Shapiro, P. R. 2006, MNRAS, 372, 679

Mellema, G., Koopmans, L., Shukla, H., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 10

Mertens, F. G., Ghosh, A., & Koopmans, L. V. E. 2018, MNRAS, 478, 3640

Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662

Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955

Mevius, M., Mertens, F., Koopmans, L. V. E., et al. 2022, MNRAS, 509, 3693

Mirocha, J., Harker, G. J. A., & Burns, J. O. 2013, ApJ, 777, 118

Mirocha, J., Muñoz, J. B., Furlanetto, S. R., Liu, A., & Mesinger, A. 2022, MNRAS, 514, 2010

Mitra, S., Choudhury, T. R., & Ferrara, A. 2015, MNRAS, 454, L76

Mondal, R., Fialkov, A., Fling, C., et al. 2020, MNRAS, 498, 4178

Morales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120

Muñoz, J. B., & Loeb, A. 2018, Nature, 557, 684

Munshi, S., Mertens, F. G., Koopmans, L. V. E., et al. 2024, A&A, 681, A62

Nebrin, O., Ghara, R., & Mellema, G. 2019, JCAP, 2019, 051

Pagano, M., & Liu, A. 2020, MNRAS, 498, 373

Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933

Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, ApJ, 788, 106

Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, ApJ, 838, 65

Planck Collaboration VI. 2020, A&A, 641, A6

Price, D. C., Greenhill, L. J., Fialkov, A., et al. 2018, MNRAS, 478, 4193

Pritchard, J. R., & Furlanetto, S. R. 2007, MNRAS, 376, 1680

Pritchard, J. R., & Loeb, A. 2012, Rep. Progr. Phys., 75, 086901

Ross, H. E., Dixon, K. L., Ghara, R., Iliev, I. T., & Mellema, G. 2019, MNRAS, 487, 1101

Ross, H. E., Giri, S. K., Mellema, G., et al. 2021, MNRAS, 506, 3717

Shaw, A. K., Bharadwaj, S., & Mondal, R. 2020, MNRAS, 498, 1480

Shaw, A. K., Chakraborty, A., Kamran, M., et al. 2023a, JApA, 44, 4

Shaw, A. K., Ghara, R., Zaroubi, S., et al. 2023b, MNRAS, 522, 2188

Shimabukuro, H., Mao, Y., & Tan, J. 2022, Res. Astron. Astrophys., 22, 035027

Singh, S., Subrahmanyan, R., Udaya Shankar, N., et al. 2017, ApJ, 845, L12

Thomas, R. M., & Zaroubi, S. 2011, MNRAS, 410, 1377

Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, 7

Trac, H. 2018, ApJ, 858, L11

Trott, C. M., Jordan, C. H., Midgley, S., et al. 2020, MNRAS, 493, 4711

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Wayth, R. B., Tingay, S. J., Trott, C. M., et al. 2018, PASA, 35, 33

Xu, W., Xu, Y., Yue, B., et al. 2019, MNRAS, 490, 5739

Zaroubi, S. 2013, in The First Galaxies, eds. T. Wiklind, B. Mobasher, & V. Bromm, Astrophys. Space Sci. Lib., 396, 45

Zel'dovich, Y. B. 1970, A&A, 5, 84