# Reionization Parameter Inference from 3D Minkowski Functionals of the 21 cm
Signals

Kangning Diao Department of Astronomy, Tsinghua University, Beijing 100084,
China  Zhaoting Chen Jodrell Bank Centre for Astrophysics, School of Physics
and Astronomy, The University of Manchester, Manchester M13 9PL, UK  Institute
for Astronomy, The University of Edinburgh, Royal Observatory, Edinburgh EH9
3HJ, UK  Xuelei Chen Department of Physics, College of Sciences, Northeastern
University, Shenyang 110819, China  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100101, China  Key Laboratory of Radio
Astronomy and Technology, Chinese Academy of Sciences, Beijing 100101, China
Yi Mao Department of Astronomy, Tsinghua University, Beijing 100084, China

(Received XXX; Revised YYY; Accepted ZZZ)

###### Abstract

The Minkowski Functionals (MFs), a set of topological summary statistics, have
emerged as a powerful tool for extracting non-Gaussian information. We
investigate the prospect of constraining the reionization parameters using the
MFs of the 21 cm brightness temperature field from the epoch of reionization
(EoR). Realistic effects, including thermal noise, synthesized beam, and
foreground avoidance, are applied to the mock observations from the radio
interferometric array experiments such as the Hydrogen Epoch of Reionization
Array (HERA) and the Square Kilometre Array (SKA). We demonstrate that the MFs
of the 21 cm signal measured with SKA-Low can be used to distinguish different
reionization models, whereas the MF measurement with a HERA-like array cannot
be made accurately enough. We further forecast the accuracies with which the
MF measurements can place constraints on reionization parameters, using the
standard MCMC analysis for parameter inference based on forward modeling. We
find that for SKA-Low observation, MFs provide unbiased estimations of the
reionization parameters with accuracies comparable to the power spectrum (PS)
analysis. Furthermore, joint constraints using both MFs and PS can improve the
constraint accuracies by up to $30\%$ compared to those with the PS alone.
Nevertheless, the constraint accuracies can be degraded if the EoR window is
shrunk with strong foreground avoidance. Our analysis demonstrates the promise
of MFs as a set of summary statistics that extract complementary information
from the 21 cm EoR field to the two-point statistics, which suggests a strong
motivation for incorporating the MFs into the data analysis of future 21 cm
observations.

Cosmology(343) — Reionization(1383) — 21 cm lines(690) — Markov chain Monte
Carlo(1889)

††journal: ApJ††software: 21cmFAST(Mesinger et al., 2011; Murray et al.,
2020), 21CMMC(Greig & Mesinger, 2015, 2017, 2018), 21cmSense(Pober et al.,
2013, 2014), Tools21cm(Giri et al., 2020), OSCAR, HERA-sim.

##  1 Introduction

The formation and evolution of the very early galaxies is a crucial part of
cosmic evolution history. With the unprecedented sensitivity of instruments
such as the James Webb Space Telescope (JWST) (Gardner et al., 2023), for the
first time a large population of very early galaxies at $z\gtrsim 10$ are
observed with great details (see, e.g. Castellano et al. 2022; Finkelstein et
al. 2023; Naidu et al. 2022; Adams et al. 2023; Austin et al. 2023). These
observations draw attention to the lack of understanding in early structure
formation, as the mass and the luminosity of the early galaxies show potential
inconsistencies with predictions from theory (Harikane et al., 2023b, 2022, a;
Bouwens et al., 2023). Observations of the very early Universe can therefore
serve as a powerful probe for cosmology and are indispensable to advance the
theory of cosmic structure formation.

The formation of early galaxies leads to the cosmic epoch of reionization
(EoR; see Furlanetto et al. 2006 for a review). The baryonic matter in the
Universe after recombination mainly consists of neutral hydrogen atoms (H i),
taking up $74\%$ of the mass density of baryons (Dodelson & Schmidt, 2020).
After the formation of the luminous objects, ultraviolet (UV) photons are
copiously produced, ionizing the H i in the intergalactic medium (IGM). As a
result, H i in the IGM serves as a tracer of the EoR (e.g. Gnedin & Ostriker,
1997; Gnedin & Shaver, 2004; Furlanetto et al., 2006; Morales & Wyithe, 2010;
Pritchard & Loeb, 2012), and is directly linked to the luminosity and spatial
distribution of the early galaxies (Muñoz et al., 2022).

The redshifted 21 cm line, arising from the hyperfine transition of atomic
hydrogen (Hellwig et al., 1970), can be observed by radio telescopes. While
individual galaxies can not be resolved, high-redshift 21 cm experiments are
capable of covering large cosmological volumes to the Gigaparsec scale, ideal
for cosmological constraints. Experiments targeting the 21 cm line during EoR
include the Low Frequency Array (LOFAR; van Haarlem et al., 2013; Yatawatta et
al., 2013), Murchison Widefield Array (MWA; Tingay et al., 2013), Precision
Array for Probing the Epoch of Reionization (PAPER; Parsons et al., 2010), and
Hydrogen Epoch of Reionization Array (HERA; DeBoer et al., 2017). In the
future, the low-frequency array of the Square Kilometre Array (SKA; Braun et
al., 2019) will be capable of producing high-fidelity images for EoR
measurements (Mellema et al., 2013).

These 21 cm experiments will likely measure the H i distribution in terms of
two-point statistics and several more sophisticated statistical properties
extracted from the tomographic 21 cm images. The 21 cm power spectrum (PS),
which is well-studied both in observation (see, e.g. Cheng et al., 2018;
Lanman & Pober, 2019; Trott et al., 2020; Abdurashidova et al., 2022) and in
theory (see, e.g. Greig & Mesinger, 2015; Li et al., 2019; Pagano & Liu,
2020), traces the underlying dark matter distribution at large scales (see,
e.g. Xu et al. 2019; Georgiev et al. 2022). The signals directly measured by
radio interferometers are visibilities in Fourier space (Smirnov, 2011),
thereby naturally linked to the PS (e.g. Morales & Wyithe 2010). However, the
growth of the ionized bubbles induces non-Gaussianity into the 21 cm
distribution (see, e.g. Ciardi et al., 2003; Mellema et al., 2006; Shimabukuro
et al., 2015; Watkinson & Pritchard, 2015; Mondal et al., 2016; Hutter et al.,
2020; Shaw et al., 2020) due to reionization patchiness that can be as large
as the scale of tens of Megaparsec, as almost no 21 cm signal is emitted from
the ionized regions. Thus, non-Gaussian statistics are required to give the
full description of the clustering in the IGM during reionization.

Limited by the foreground and the principle of interferometry, the 21 cm
tomographic data measures the fluctuations, not the total signal — the
positive signals in the data indicate the neutral structures and the negative
signals trace the ionized bubbles. Thus the non-Gaussianity affects the
clustering of the 21 cm signal through two aspects — morphology and topology
of ionized bubbles. Non-Gaussian morphology can be characterized in terms of
higher-order correlation functions such as three-point correlation function
(Hoffmann et al., 2019) in configuration space and bispectrum (e.g. Bharadwaj
& Ali, 2004; Yoshiura et al., 2015; Shimabukuro et al., 2017; Hutter et al.,
2020; Watkinson et al., 2022; Tiwari et al., 2022) and trispectrum (e.g.
Lewis, 2011; Shaw et al., 2019) in Fourier space. While these statistics have
been explored in cosmological theories (e.g. Adshead et al., 2013), these
correlation functions do not fully characterize the field distribution in a
highly non-Gaussian case (Carron, 2011).

The size statistics also characterize the morphological information. These
include the bubble size distribution (Kakiichi et al., 2017; Giri et al.,
2018a; Busch et al., 2020; Shimabukuro et al., 2022) and neutral island
distributions (Xu et al., 2014, 2017; Giri et al., 2019a), or feature in these
distributions such as the largest cluster statistics(Pathak et al., 2022;
Dasgupta et al., 2023). Previous works have shown from simulations that the
EoR undergoes a phase transition (Furlanetto & Oh, 2016) during which the
ionized bubbles are changed from being small and isolated initially to being
connected and percolated as the bubbles expand and merge.

On the other hand, unlike the conventional Fourier-space polyspectra methods,
topological statistics, which are intrinsically defined in configuration space
and measured from tomographic data, provide complementary information to the
morphology (e.g. Pratten & Munshi, 2012). Topological properties of the field
have been characterized in terms of the Minkowski Functionals (MFs; e.g.
Gleser et al., 2006; Wang et al., 2015; Yoshiura et al., 2017; Kapahtia et
al., 2019), shape finders (Bag et al., 2019) and Betti numbers (Giri &
Mellema, 2021). The EoR can be divided into five stages by quantifying the
percolation of ionized bubbles and neutral islands with these statistics (Chen
et al., 2019).

It is proved and well known that, for an excursion set field, the MFs contain
the complete global topological information (Hadwiger, 1959). As such, the MFs
are powerful tools for characterizing the topology of cosmological signals.
Thus the MFs have been applied to the data analysis in cosmic microwave
background (CMB) (e.g. Novikov et al., 1999; Buchert et al., 2017) and the
large-scale structure surveys (e.g. Hikage et al., 2006; Appleby et al.,
2022).

In this work, we explore the detectability of MFs, taking into account the
resolution and instrumental effects of interferometric array experiments. We
consider the configurations of a redundant, HERA-like array as well as the
SKA. We study the evolution of the observed MFs during different stages of
EoR. As a demonstration of concept, we extend the
21CMMC111https://github.com/21cmfast/21CMMC code (Greig & Mesinger, 2015,
2017, 2018) as an MCMC sampler to infer the EoR parameters from the MFs of the
21 cm brightness temperature field. Two reionization models, the faint
galaxies and bright galaxies of Greig & Mesinger (2017, 2018), are used to
study parameter inferences with MFs in different scenarios.

In previous work, Kapahtia et al. (2021) exploited topological statistics,
specifically Contour Minkowski Tensors and Betti numbers in 2D images, to
constrain astrophysical parameters from the EoR. In comparison, our work
includes observational effects such as the foreground wedge cut to perform a
more realistic forecast. In addition, our analysis is not only based on the
information in 3D images that contain the information along the LOS but also
exploits the MFs as a function of different thresholds as opposed to averaging
over threshold values as in Kapahtia et al. (2021). More modes are therefore
included in our analysis.

The remainder of this paper is organized as follows. In Section 2, we briefly
summarize the theoretical framework to generate mock observations, including
the 21 cm brightness field and the instrumental effects. In Section 3, we
discuss the evolution and underlying information of the observed MFs. We then
present the parameter fitting results with the MFs coupled with different
instrumental effects in Section 4 and make concluding remarks in Section 5. We
adopt the values of cosmological parameters in this paper as follows,
($\Omega_{\Lambda}$, $\Omega_{\rm M}$, $\Omega_{b}$, $n$, $\sigma_{8}$, $h$) =
(0.69, 0.31, 0.048, 0.97, 0.81, 0.68), in consistency with the Planck 2015
results (Planck Collaboration et al., 2016).

##  2 Simulations and Mock Observations

###  2.1 Simulations of the 21 cm signals

The fluctuation of 21 cm emission is represented by its flux density, or
equivalently the 21 cm brightness temperature $T_{b}$. In the optically thin
approximation, it can be written as (e.g. Furlanetto et al., 2006; Mellema et
al., 2013)

$\displaystyle\delta T_{b}(\nu)=27x_{{\rm
HI}}(1+\delta_{m})\bigg{(}1-\frac{T_{\gamma}}{T_{{\rm S}}}\bigg{)}$ (1)
$\displaystyle\times\bigg{(}\frac{1+z}{10}\frac{\Omega_{m}}{0.27}\bigg{)}^{1/2}\bigg{(}\frac{\Omega_{b}h}{0.044\times
0.7}\bigg{)}$
$\displaystyle\times\bigg{(}1+\frac{1}{a\,H(z)}\frac{\mathrm{d}v_{\|}}{\mathrm{d}r_{\|}}\bigg{)}^{-1}\quad\mathrm{mK},$

where $x_{\rm HI}$ is the neutral fraction of the IGM,
$\delta_{m}({\mathbf{x}},z)\equiv\rho/\bar{\rho}-1$ is matter density
fluctuations, $T_{\gamma}$ is the CMB temperature and $T_{\rm S}$ is the spin
temperature. Here we focus on the stage where the gas has been thoroughly
heated so that $T_{\rm S}\gg T_{\gamma}$. The signal is also affected by the
proper gradient along the line of sight (LOS) of the peculiar velocity,
$\mathrm{d}v_{\|}/\mathrm{d}r_{\|}$.

Figure 1: Illustrations of the instrumental effects in mock SKA and HERA-like
observations. The top left panel shows a slice of the 21 cm signal at $z=7$
(corresponding to $\bar{x}_{\rm HI}=0.32$) in the faint galaxies model. Panels
in the middle (bottom) row show the 21 cm signals that are applied with the
instrumental effects with the HERA-like (SKA) observations. In each of the
middle and bottom rows, panels from the left to the right column show the 21
cm signals that are convolved with the PSF (first column), a realization of
thermal noise simulations (second column), the PSF-convolved 21 cm signals
with thermal noise applied (third column), the PSF-convolved 21 cm signals
that are applied with thermal noise and foreground wedge cut (fourth column),
respectively.

We use the publicly available code
21cmFAST222https://github.com/21cmfast/21cmFAST (Mesinger et al., 2011; Murray
et al., 2020)v3.0.3 to simulate the 21 cm $T_{b}$ field, and the main
procedure can be briefly summarized as follows. First, a high-resolution
initial condition (IC) is generated through linear perturbation theory. Then
the second-order Lagrangian perturbation theory (LPT, e.g. Scoccimarro, 1998)
is applied to evolve IC to the density and velocity fields at lower redshifts,
then smoothing them to lower resolution. Next, following the excursion set
theory (Furlanetto et al., 2004), 21cmFAST estimates the ionization field by
comparing the reciprocal of the ionizing efficiency $\zeta$ and the collapse
fraction $f_{\rm coll}(\mathbf{x},R,z)$ within a spherical volume. Here,
$f_{\rm coll}$ is the fraction of mass collapsed into halos, and $\zeta$ is an
aggregation of several physical parameters to describe the photon production
rate per baryon in a dark matter halo as described in Mesinger et al. (2011).

A spherical region around the position $\mathbf{x}$ with the radius $R$ at the
redshift $z$ is fully ionized if the condition,

$f_{\rm coll}(\mathbf{x},R,z)\geq\zeta^{-1}\,,$ (2)

is satisfied. Following this recipe, 21cmFAST iterates through every cell in
the simulation. For each cell, $R$ is initially set to $R_{\rm mfp}$, the mean
free path of the ionizing photons or technically the largest possible radius
of ionized bubbles. The criterion in Equation (2) is then repeatedly checked
from large to small radii until a radius $R$ is found to meet the criterion,
in which case the central cell of the spherical region is set to be fully
ionized. If no radius greater than the cell size $R_{\rm cell}$ meets the
criterion, then the cell is partially ionized with the ionized fraction set to
be $f_{\rm coll}(\mathbf{x},R_{\rm cell},z)\,\zeta$.

The reionization process can then be simulated by modeling $\zeta$ and $f_{\rm
coll}$. In this paper, we focus on a two-parameter model, with the minimum
virial temperature ($T_{\rm vir}$) and the ionizing efficiency ($\zeta$). The
parameter $T_{\rm vir}$ represents the threshold temperature required for a
halo to host ionizing sources, effectively tuning the collapse function
$f_{\rm coll}$. Higher values of $T_{\rm vir}$ indicate a reduced number of
ionizing sources, while larger $\zeta$ values lead to a faster expansion of
ionized bubbles. Thus, there exists a degeneracy between $T_{\rm vir}$ and
$\zeta$ in terms of their impact on the ionization history. However, it has
been noted by Yoshiura et al. (2015) that despite this degeneracy, the
morphology of the ionized regions exhibits distinct characteristics associated
with different parameter combinations, as manifested in the MFs.

Table 1: Telescope parameters at redshift $z=7$ to generate mock observations. Parameters |  HERA-like |  SKA-Low
---|---|---
Antenna number |  331 |  512
System temperature ($T_{\rm sys}$) [K] |  $100+60\left(\frac{\nu}{300\rm MHz}\right)^{-2.55}$ |  $100+60\left(\frac{\nu}{300\rm MHz}\right)^{-2.55}$
Dish Area ($A_{\rm Dish}$) [m2] |  153.93 |  962.11
Integration time ($t_{\rm int}$) [s] |  10 |  10
Declination |  -38∘ |  -30∘
Critical frequency ($\nu_{\rm c}$) [MHz] |  110 |  110
Frequency resolution ($\Delta\nu$) [MHz] |  0.097 |  0.097
Default wedge slope ($m$) |  0.5 |  0.35
PSF FWHM [Mpc] |  25 |  3
Averaged thermal noise erroraaThis is the rms of smoothed thermal noise with 1000 hours observation at $z=7$. The SKA-Low result is consistent with other estimations (e.g. Giri et al., 2018b). ($\sqrt{\left\langle N^{2}\right\rangle}$) [mK] |  1.37 |  1.01

In this work, all our semi-numerical simulations have a comoving size of
$512\times 512\times 128$ Mpc3 with $\rm(2\,Mpc)^{3}$ cells. The length of
$512$ Mpc in the transverse plane corresponds to an image size of $\sim
3^{\circ}$, matching the field-of-view (FoV) of the instruments. The length of
128 Mpc along the line-of-sight corresponds to a redshift interval of $\Delta
z\sim 0.5$, representing the observations in a frequency sub-band. In Sections
2 and 3, we adopt a choice of parameters $\zeta=30$ and $\rm{log}[T_{\rm
vir}/K]=4.70$, corresponding to the faint galaxies model in Greig & Mesinger
(2018), as the benchmark model that shows the characteristics of MFs. The top
left panel of Figure 1 shows a slice of simulated brightness temperature at
$z=7$. In Section 4 and afterward, different choices of model parameters using
both faint galaxies and bright galaxies models are discussed.

###  2.2 Mock Observations

In this section, we use the simulated 21 cm signal to generate mock
observations with instrumental effects. The minimum scales probed are limited
by the point spread function (PSF), which for interferometric observation is
the synthesized beam of the weighted distribution of interferometric baselines
(Condon & Ransom, 2016). The observed signal also has noise, which is assumed
to be thermal here. The resulting image will be dominated by the radio
foregrounds, which can be avoided by applying a foreground wedge cut (Liu et
al., 2014).

####  2.2.1 Thermal Noise

In an interferometer array, each pair of elements forms a baseline, denoted by
its vector components in units of wavelength $(u,v)$, and the interferometric
visibility data is recorded independently, contributing to a point in the $uv$
plane. The thermal noise error can be drawn from a Gaussian distribution with
the following root-mean-square (rms) for every single baseline (e.g.
Furlanetto et al., 2006; Giri et al., 2018b)

$\sigma_{N}^{\rm single}=\frac{\lambda^{2}T_{\rm sys}}{\epsilon\Omega A_{\rm
Dish}\sqrt{2\Delta\nu t_{\rm int}}},$ (3)

where $\Omega\approx\theta^{2}$ is the beam area, $\lambda$ is the signal
wavelength, $T_{\rm sys}$ is the system temperature of the telescope,
$\Delta\nu$ is the frequency resolution, $t_{\rm int}$ is the integration
time, $A_{\rm Dish}$ is the physical area of station and $\epsilon$ is an
efficiency factor which follows

$\displaystyle\begin{split}\epsilon=\left\\{\begin{array}[]{ll}1,&\nu<\nu_{\rm
c}\\\ \left(\frac{\nu_{\rm c}}{\nu}\right)^{2},&\nu\geq\nu_{\rm
c}.\end{array}\right.\end{split}$ ()

When calculating the thermal noise, we discretize the $uv$ space into 2D
grids. Assuming $N_{\rm uv}$ baselines fall in the same pixel in one night,
and the final visibility data is the coherent average of $N_{\rm d}$ nights,
then the noise rms of this pixel is calculated as

$\sigma_{N}^{\rm pixel}=\frac{\sigma_{N}^{\rm single}}{\sqrt{N_{\rm
d}N_{uv}}},$ (4)

where we assume that natural weighting is applied when producing the image
cube. Throughout this paper, we use natural weighting, i.e. the grid weights
are the number of baselines in each $u$-$v$ grid.

Here we consider two telescope configurations, HERA-like and SKA-Low, whose
parameters are summarized in Table 1. We calculate the SKA-Low thermal noise
through the publicly available code Tools21cm333https://github.com/sambit-
giri/tools21cm (Giri et al., 2020), following the procedure in Giri et al.
(2018b) and Hassan et al. (2019). We generate the HERA-like thermal noise in
the same manner, based on Tools21cm and 21cmSense444https://github.com/steven-
murray/21cmSense (Pober et al., 2013, 2014). We refer to the latter settings
as “HERA-like” for reasons discussed later in Section 2.2.2. Slices of the
simulated 1000 hours HERA-like and SKA-Low noise are shown in Figure 1.

####  2.2.2 Point Spread Function

The interferometer measures the Fourier modes of the sky signal. The
visibility data in the $u$-$v$ space are averaged into the $u$-$v$ grids and
then transformed back to the angular space for imaging. Since the gridded
visibility is the multiplication of the sky signal in Fourier space with the
grid weights, the images are then the convolution of the sky signal and the
synthesized beam, which is the Fourier transform of the grid weights (Condon &
Ransom, 2016). The synthesized beam is also referred to as the point spread
function (PSF). In the case of natural weighting, the PSF is the renormalized
Fourier transform of the baseline distribution.

Figure 2: The PSFs used for the mock observations at $z=8$ using SKA-Low
(left) and for HERA-like observations (right), respectively.

For SKA-Low, we use the v3 station layout (de Lera Acedo et al., 2020) to
simulate the baseline distribution of a typical tracking observation. We
choose the phase center to be the EoR0 field (Lynch et al., 2021) at ${\rm
RA}=0\,{\rm h}$, ${\rm Dec}=-30\,{\rm deg}$. The observation time of one track
is assumed to be a total of 4 hours with a time resolution of 10 seconds. The
coherent average of multiple nights is assumed to reach the total integration
time of 1000 hours. The baseline distribution is simulated using the
OSKAR555https://github.com/OxfordSKA/OSKAR code (Mort et al., 2010). The
baseline distribution is then gridded into $256\times 256$ u-v grids, with the
cell size corresponding to the 512 Mpc box size in the transverse plane. The
baseline distribution is then Fourier transformed to obtain the PSF. The PSF
at $z=8$ is shown in the left panel of Figure 2 for the purpose of
illustration, with the corresponding spatial resolution of the PSF listed in
Table 1.

HERA, on the other hand, is a redundant baseline array designed for redundant
calibration and power spectrum estimation (Dillon et al., 2020). Here, we do
not intend to simulate realistic drift scan observations of HERA, but to use
its array layout to generate PSFs that are qualitatively different from SKA-
Low observations for comparison. Throughout this paper, we refer to this
scenario as “HERA-like” observations. We use the
hera_sim666https://github.com/HERA-Team/hera_sim code to retrieve the station
layout to calculate the baseline coordinates. We then symmetrize the baseline
distribution, as a crude approximation of a coherent average of different
nights when the same patch of the sky passes through. The baseline
distribution is then gridded to obtain the PSF. The resulting PSF at $z=8$ is
shown in the right panel of Figure 2.

The simulated PSFs are then convolved with the signal cubes. Note that we do
not include primary beam attenuation in generating the mock observation, since
the $512$ Mpc box length is within the primary beam of the instruments. In our
case, the primary beam attenuation is simply multiplied in the signal, which
is then divided out. This step only increases the thermal noise level of the
pixels. Therefore, we note to the readers that the thermal noise simulated can
only be achieved with an integration time larger than 1000 hours in real
observations. The chromaticity of the primary beam may introduce additional
structure in the observed 21 cm signal, whose effects are beyond the scope of
this work.

####  2.2.3 Foreground Wedge

The major challenge for the EoR 21 cm line observation is the contamination of
the foregrounds, which is four to five orders of magnitude higher than the
signal. We follow the foreground avoidance method by filtering out the
foreground-dominated $k$-modes outside the EoR window. A thorough explanation
of the EoR window can be found in e.g. Liu & Shaw (2020), and we briefly
summarize the main idea below.

In the cylindrical 2D Fourier (i.e. $k_{\parallel}$-$k_{\perp}$) plane, where
$k_{\parallel}$ ($k_{\perp}$) is the comoving wave number parallel
(perpendicular) to the LOS, respectively, the region satisfying
$k_{\parallel}<mk_{\perp}$ is dominated by the foregrounds, but the regions
where $k_{\parallel}>mk_{\perp}$ (i.e. “EoR window”) is believed to have
negligible foreground contamination so it can be used to study the 21 cm
signal. Here m is the so-called “foreground wedge” slope

$m=\frac{DH_{0}E(z)\sin\theta}{c(1+z)},$ (5)

where $H_{0}$ is the Hubble constant, $D$ is the comoving distance,
$E(z)=H(z)/H_{0}$, the ratio of Hubble parameter at $z$ and Hubble constant at
present, and $\theta$ is the characteristic angular distance outside of which
the foreground emission is negligible due to power beam attenuation. In the
optimistic case, the foreground emission is contained within the primary beam
FoV of the instrument.

It is worth noting that, the feature of foreground wedge arises from the
chromaticity of the baselines, i.e. the fact that the angular scale measured
by a baseline changes with frequency (see e.g. Figure 5 of Murray & Trott
2018). For tomographic imaging where the visibilities are averaged in a map-
making process, the chromaticity no longer exists and there is no well-defined
foreground wedge (see e.g. Hothi et al. 2021). However, the residual
foreground cleaning process is expected to remove the foreground power at high
$k_{\parallel}$ and effectively recover an observation window above the wedge
(see e.g. Chen et al. 2023). Therefore, we expect a Fourier filter based on
the foreground wedge criterion can be used to extract the H i signal in the
images.

For EoR observations, telescopes with wide FoV are typically used. For
instruments like HERA and SKA-Low, the beam sidelobes are non-negligible up to
the physical horizon. The chromaticity of the sidelobes introduces foreground
contamination at high $k_{\parallel}$ modes, lifting the slope of the
foreground wedge (Pober et al., 2016). The sky signal is cut off at the
physical horizon, and therefore there is a large wide-field effect around the
horizon line $\sin\theta=1$ (Thyagarajan et al., 2015). These effects need to
be carefully studied and mitigated in real observations and are beyond the
scope of our work. We assume the optimistic scenario where the foreground
wedge is suppressed to the angular scale of the primary beam with the values
listed in Table 1. Larger values of $m$ are also explored later to examine the
impact of the position of the foreground wedge on MFs in Section 3.3. Finally,
we note that identification of ionized regions in the presence of foregrounds
may be achieved through machine learning techniques for future SKA-Low
observations (Bianco et al., 2024), which we leave for future work.

##  3 Minkowski Functionals of the 21 cm Brightness Temperature Field

###  3.1 Definition of Minkowski Functionals

Figure 3: The MFs around transition points with SKA-Low noise between
different stages. We show four MFs (from top to bottom) – $V_{0}$,
$V_{1}$,$V_{2}$, and $V_{3}$, respectively, around four different transition
points (from left to right). The grey-dotted, red solid, and blue dashed lines
represent the average MFs before, during, and after the transition point. The
1$\sigma$ errors of each line are shown in the shaded region with the same
color.

The MFs (Minkowski, 1903) are a set of functionals describing the morphology
of contour surfaces of a field. In three-dimensional space, MFs consist of
four independent functionals (Schmalzing & Buchert, 1997), which are defined
as follows:

$\displaystyle V_{0}(\nu)$ $\displaystyle=\frac{1}{V}\int_{V}\mathrm{d}^{3}x\
\Theta[u(\boldsymbol{x})-\nu\sigma]$ (6) $\displaystyle V_{1}(\nu)$
$\displaystyle=\frac{1}{6V}\int_{\partial F_{\nu}}\mathrm{d}s$ (7)
$\displaystyle V_{2}(\nu)$ $\displaystyle=\frac{1}{6\pi V}\int_{\partial
F_{\nu}}\mathrm{d}s\ [\kappa_{1}(\boldsymbol{x})+\kappa_{2}(\boldsymbol{x})]$
(8) $\displaystyle V_{3}(\nu)$ $\displaystyle=\frac{1}{4\pi V}\int_{\partial
F_{\nu}}\mathrm{d}s\ \kappa_{1}(\boldsymbol{x})\kappa_{2}(\boldsymbol{x})$ (9)

where $u(\boldsymbol{x})$ is the field being studied, $\sigma$ is the rms of
the field, $\nu\equiv u_{\rm thr}/\sigma$ represents the relative threshold of
the field, $V$ is the volume that contains the field, $\Theta$ represents the
Heaviside step function, $\partial F_{\nu}$ denotes the contour surface whose
threshold value is $u_{\rm thr}$, and $\kappa_{1}$ and $\kappa_{2}$ are the
principle curvatures of the $\partial F_{\nu}$. The morphological meanings of
$\\{V_{0},V_{1},V_{2},V_{3}\\}$ are the volume inside, surface area, mean
curvature and the Euler characteristics of $\partial F_{\nu}$. Here, “inside”
denotes the region where $u(\boldsymbol{x})>u_{\rm thr}$. The $V_{3}$ is
related to the genus number $g$ by $g=1-V_{3}/2$, which describes the homotopy
of the field by $gV=(\rm number\ of\ holes)-(\rm number\ of\ isolated\
islands)+1$.

We take the Kodenderink method (Koenderink, 1984) to compute MFs. In this
method, $\kappa_{1}+\kappa_{2}$ and $\kappa_{1}\kappa_{2}$ are represented by
the first- and second-order partial derivatives of the field, whereby these
integrals can be transformed into volume integrals. We developed a Python
package pyMIN777https://github.com/dkn16/pyMIN to perform the calculation.
Unless otherwise mentioned, we calculate the MFs from $-3\sigma$ to $3\sigma$
with a sampling interval of $0.1\sigma$, thus including 61 data points in
total for each functional.

A special case of the MFs is that of the Gaussian random field, which can be
analytically expressed as follows.

$\displaystyle V_{0}(\nu)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\Phi\left(\frac{\nu}{\sqrt{2}}\right)$
(10) $\displaystyle V_{1}(\nu)$
$\displaystyle=\frac{2}{3}\frac{\lambda}{\sqrt{2\pi}}e^{-\nu^{2}/2}$ (11)
$\displaystyle V_{2}(\nu)$
$\displaystyle=\frac{2}{3}\frac{\lambda^{2}}{\sqrt{2\pi}}\nu e^{-\nu^{2}/2}$
(12) $\displaystyle V_{3}(\nu)$
$\displaystyle=\frac{\lambda^{3}}{\sqrt{2\pi}}(\nu^{2}-1)e^{-\nu^{2}/2}$ (13)

where $\Phi$ is the Gaussian error function,
$\lambda\equiv\sqrt{\sigma_{1}/6\pi\sigma}$ and $\sigma_{1}$ is the rms of the
$|\nabla u|$ field. These functions provide a template for comparison. In
particular, MFs closer to the Gaussian analytical form represent a more
Gaussian field.

Throughout this paper, the 21 cm brightness temperature field is normalized by
its standard deviation before calculating the MFs, a different treatment from
that of Chen et al. (2019) where the MFs are calculated directly on the
temperatures. As we focus on mock observations with different levels of
thermal noise and reionization models, the normalization procedure sets the
comparison between different scenarios on a fairer basis.

###  3.2 Evolution during Reionization

We begin by investigating the evolution of MFs during the reionization process
through simulated mock observations. In our study, we employ a fiducial
observation scenario involving a 1000h integration time. In this section, we
assume no foregrounds as the most optimistic scenario and assume the
configuration of SKA-Low to generate the mock observations, but leave the
discussion of the impact of filtering out foreground wedge to Section 3.3.

To investigate the evolution of MFs, we examine the critical redshifts before,
during, and after each of the four transition points described in Chen et al.
(2019). The average neutral fractions of the IGM at the transition points are
$\bar{x}_{{\rm HI}}\approx 0.9$, 0.7, 0.3, and 0.15. For detailed discussions
of the different stages of EoR and the transitions, see Chen et al. (2019).
For every redshift, we construct ten independent mock observations and
calculate the corresponding mean and the sampling variance of the MFs.

To untangle the effects of noise and PSF, we start by removing the thermal
noise component in the mock observation and investigate the noise-free signal
box convolved with the PSF. We find that the noise-free scenario exhibits
consistent characteristics as in Chen et al. (2019), as the configuration of
SKA-Low is capable of resolving the small-scale structure of the 21 cm image
shown in the bottom left panel of Figure 1. To avoid redundancy, we choose not
to present our result of MFs in the noise-free scenario and refer the
interested readers to Figure 2 of Chen et al. (2019). We then proceed to
calculate the MFs convolved with PSF and assume the SKA-Low thermal noise with
1000h observation time, and show the results in Figure 3. Each line in the
figure represents the mean MFs along with the shaded regions being their
$1\sigma$ errors, derived from ten independent realizations. The presence of
noise impedes a direct correlation between the evolution of MFs and that of
ionized bubbles. Despite the presence of biases in all MFs, the relatively
small $1\sigma$ errors indicate that the noised MFs is dominated by the 21 cm
signal. Consequently, these biased MFs still exhibit distinct features that
trace the evolution of the EoR which we summarize as follows.

During the transition from the ionized bubble stage (Stage 1) to the ionized
fiber stage (Stage 2) ($\bar{x}_{\rm HI}\sim 0.90$), the reionization process
remains largely Gaussian, corresponding to relatively small bubbles of ionized
regions. As the temperature field is normalized, MFs at $\bar{x}_{\rm
HI}\gtrsim 0.90$ are almost identical. The application of PSF smoothing and
the presence of Gaussian thermal noise obscures the subtle non-Gaussian
signatures that start to emerge at $\bar{x}_{\rm HI}\sim 0.90$, originating
from the growth of the ionized bubbles. We conclude that the first transition
point of morphology described by the MFs can not be observed in realistic
observations.

As the reionization process progresses, distinct features begin to emerge
within the noised MFs. The transition from the ionized fiber stage to the
sponge stage (Stage 2 to Stage 3, $\bar{x}_{\rm HI}\sim 0.70$) coincides with
the emergence of large connected ionized regions. The MFs deviate from the
analytical Gaussian form and start to evolve with the evolution of neutral
fraction, demonstrating the non-Gaussianity of the EoR. Distinctive features
of the peaks and troughs in $V_{1,2,3}$ start to evolve. The merging of
ionized bubbles decreases the surface area of the ionizing front, as shown in
the decreasing of the peak in $V_{1}$ around $\delta T_{b}\sim 0$. Similarly,
a decrease of peak structures also appears in $V_{2}$ and $V_{3}$. From
$\bar{x}_{\rm HI}\sim 0.7$ onwards, the non-Gaussiantiy of the EoR starts to
dominate, which makes the MFs a powerful probe of its evolution.

During the transition from the sponge stage to the neutral fiber stage (Stage
3 to Stage 4, $\bar{x}_{\rm HI}\sim 0.32$), the originally decreasing peaks of
MFs start to increase as seen in the third column of Figure 3. It indicates
that the dominating morphology is no longer large neutral regions surrounding
the ionized bubbles and filament structures. Instead, the neutral IGM is
enclosed by the ionized region, forming small patches of neutral islands (Xu
et al., 2017) as reionization progresses. This shift in morphology can also be
seen by the fact that the position of the $\delta T_{b}\sim 0$ peak in $V_{1}$
and trough in $V_{3}$ shift from positive to negative values. The MFs trace
closely the evolution of the neutral islands, and the morphology of the
islands can be used to indicate the evolution of IGM towards the end of EoR
(Giri et al., 2019b).

During the transition from the neutral fiber stage to the neutral island stage
(Stage 4 to Stage 5, $\bar{x}_{\rm HI}\sim 0.16$), the majority of the IGM
become ionized, and the features of MFs are predominantly influenced by
thermal noise, mirroring the ionized bubbles at $\bar{x}_{\rm HI}\gtrsim 0.9$.
However, as seen in the fourth column of Figure 3, the MFs still demonstrate
sizeable non-Gaussian features. The neutral islands do not trace the luminous
objects and their host clusters, and therefore they are much more non-
Gaussian. The peak features in $V_{1,2,3}$ continue to grow, indicating that
the islands are increasingly disconnected as the Universe gets more ionized.

In conclusion, the evolution of the noised MFs shows mixed effects of Gaussian
smoothing from the thermal noise, as well as evolution in its features tracing
the morphology of EoR. For deep observations using SKA-Low, the signal-to-
noise level required for resolving the morphology can be reached and MFs can
be used to constrain the reionization, which we will discuss in Section 4.

###  3.3 The Impact of Foreground Wedge

Figure 4: The MF $V_{3}$ at $z=8$, assuming thermal noise with SKA-Low deep
(top) and HERA-like (bottom) observations, respectively. Here we assume
different scenarios of foreground cut — the no-foreground-cut scenario (blue
solid line), an optimistic (default) foreground wedge slope $m=0.35$ ($0.5$)
for SKA (HERA-like) observations (red solid line), and a pessimistic scenario
where the value of the wedge slope $m$ is doubled from the default value (cyan
solid line). All fields are normalized to have the same Gaussian template
(grey dotted line). We also show the 1$\sigma$ errors in shaded regions with
the same color.

The continuum foreground emission dominates the signals on large scales along
the LOS, thereby causing the loss of large-scale information. In this
subsection, we explore the impact of different observational configurations on
the observed MFs and test the robustness of MFs in the presence of
foregrounds. As mentioned in Section 2, we focus primarily on SKA-Low, while
only presenting a HERA-like case for extreme comparisons.

To avoid redundancy, in this subsection, we only show the comparisons of
$V_{3}$ in different settings of observations. The Euler characteristic
$V_{3}$ is related to the genus number of the EoR field, a well-studied
morphological indicator for EoR (e.g. Lee et al. 2008). The differences in
$V_{3}$ and their implications can be trivially extended to other MFs.

In Figure 4, we show the results of the MF $V_{3}$ at $z=8$ as an illustrative
example to highlight the characteristics of different scenarios. For a fair
comparison, we normalize each field with the standard deviation of the $m=0$
case instead of their respective standard deviations. The Gaussian template
with the unit standard deviation is also shown to demonstrate the
(non)-Gaussianity of $V_{3}$ in different scenarios, which we discuss in
detail below.

(1) No Foreground Cut Case. — We begin by considering the most optimistic
scenario where foreground contamination can be completely removed, as
indicated by the “$m=0$” case in Figure 4. In this case, the $V_{3}$ obtained
from SKA-Low exhibits notable deviations from the Gaussian template, while
$V_{3}$ with HERA-like observations is Gaussian-like.

Despite the resemblance with a Gaussian field, the amplitude of $V_{3}$ is
attenuated and the $1\sigma$ region is extended compared to the Gaussian case,
both of which arise as the consequences of PSF smoothing. For HERA-like
redundant layouts, the sky is only sampled at some fixed angular scales, with
the majority of the morphological information lost. Consequently, with less
remaining information, the field becomes more susceptible to thermal noise
fluctuations, leading to larger error bars.

(2) Foreground Effects. — We consider two scenarios of foreground
contamination. The first scenario — $m=0.35$ ($0.5$) for SKA (HERA-like)
observations, respectively — represents an optimistic case where foreground
removal algorithms successfully mitigate the foreground, bringing the wedge
down to the primary beam FoV as listed in Table 1. The second scenario
represents a pessimistic case where the foreground slope $m$ is magnified by a
factor of 2 due to incomplete foreground removal in beam sidelobes.

Across all cases with foregrounds, we observe that the loss of modes leads to
a decline in non-Gaussianity compared to the no-foreground-cut case. For the
SKA-Low configuration, the presence of the foregrounds results in fewer
features in the MFs. As expected, raising the foreground wedge results in more
Gaussian MFs with less information. Regarding the HERA-like configuration, we
observe an increase in the noise level and the peak structure, due to the loss
of information as well as the numerical artifacts from filtering out the
$k_{\parallel}$ modes.

In conclusion, we find that for SKA-Low observations, the non-Gaussian
features of MFs are robust against the presence of foregrounds, and can be
observed assuming a realistic level of foreground mitigation.

##  4 Reionization Parameter Inference with MFs

###  4.1 The Impact of Reionization Parameters

Figure 5: The average MFs and corresponding $1\sigma$ errors in shaded
regions, assuming SKA-Low noise. We consider the no foreground cut scenario
(left) and the foreground scenario with the wedge slope $m=0.7$ (right),
respectively. We show the results by varying four different sets of
reionization parameters that all yield $\bar{x}_{\rm HI}\approx 0.60$ at
$z=8$.

In this section, we investigate the impact of reionization parameters on the
observed MFs of the 21 cm signal and in turn, use the MFs to constrain the
EoR. For each scenario, we choose four sets of reionization parameters that
yield identical reionization history, with $\left(\zeta,{\rm
log_{10}[{T_{vir}/K}]}\right)=\left(15,4.29\right)$,
$\left(30,4.70\right),\left(50,4.94\right)$ and $\left(200,5.42\right)$. As
discussed in Section 2, increasing values of $\left(\zeta,{\rm
log_{10}[{T_{vir}/K}]}\right)$ lead to more inhomogeneous reionization, as the
reionization sources will be less numerous. We generate ten independent mock
realizations for each parameter set and compute the corresponding MFs. For
simplicity, we show the results at $z=8$ with $\bar{x}_{\rm HI}=0.60$, and the
conclusions can be trivially extended to other redshifts. For simplicity, in
the rest of this paper, we use $\rm{log}_{10}T_{vir}$ to denote ${\rm
log_{10}[{T_{vir}/K}]}$.

The left panel of Figure 5 shows the MFs corresponding to different parameter
values, assuming SKA-Low noise and no residual foreground. Despite the
presence of noise, the MFs successfully break the degeneracy among the
parameter sets, as evidenced by the varying shapes of the MFs in four cases.
Specifically, the MFs $\\{V_{1},V_{2},V_{3}\\}$ show substantial variations in
dependence on the parameter values, while $V_{0}$ exhibits negligible changes.
One prominent feature is that higher values of $\zeta$ and $T_{\rm vir}$ tend
to produce lower peaks in $\\{V_{1},V_{2},V_{3}\\}$. This phenomenon can be
attributed to the “faint” model, characterized by low $\zeta$ and $T_{\rm
vir}$, which generates numerous small ionized bubbles. Although some of these
small bubbles may be smoothed out due to limited resolution and thermal noise,
a substantial number remains in the signal. In contrast, the “bright” model,
characterized by higher $\zeta$ and $T_{\rm vir}$, leads to larger but fewer
ionized bubbles with less structure. Therefore, it can be seen that the MFs
are more sensitive to the reionization parameters when the reionization model
yields a larger population of small ionized bubbles.

To further test the robustness of MFs in distinguishing reionization models
with the foreground cut, we choose $m=0.7$, which represents a case of severe
foreground contamination, and calculate the MFs in the four sets of parameters
as shown in the right panel of Figure 5. As expected, the
$\\{V_{0},V_{1},V_{2},V_{3}\\}$ curves in this case appear slightly more
Gaussian compared to the no-foreground-cut scenario. Nevertheless, the
significant differences of the MFs $\\{V_{1},V_{2},V_{3}\\}$ between different
parameter sets still provide strong evidence for the robust constraining power
of MFs.

In conclusion, MFs exhibit considerable power in distinguishing different
reionization models, even in the presence of instrumental effects and
foreground cuts. This indicates that the MFs are likely to be exploited in
constraining reionization parameters — a topic we will discuss in the
remainder of this section.

###  4.2 21CMMC Setup

Now we quantify the constraining power of MFs. For the forecast study, we
extend the 21CMMC package (Greig & Mesinger, 2015) and use MFs as the summary
statistics for constraining EoR. We describe the procedure below.

Figure 6: The posterior distribution of EoR parameters using the SKA-Low
configuration with 1000 h observation for the faint galaxies (left) and bright
galaxies (right) model, respectively. The lower left plot in each panel shows
the $1\sigma$ (dark-shaded region) and $2\sigma$ region (light-shaded region),
with the results from MFs shown in red and PS shown in blue. The dashed lines
indicate the fiducial parameter values. The lower-right plot and upper-left
plot in each panel show the marginalized probability distribution, where
$1\sigma$ intervals are shown as the shaded regions.

Figure 7: Same as Figure 6 but for the HERA-like configuration with 1000 h
observation (and with different axis scales).

Our mock observations are $512\times 512\times 128$ Mpc3 coeval cubes on a
$256\times 256\times 64$ grid, the same as the coeval cube sampled within
MCMC. Then the coeval cube and noise cube are all smoothed to the same
telescope configuration before calculating MFs.

In this proof-of-concept work, we adopt a simple Gaussian likelihood. Ignoring
the correlation between any two points in MFs, the likelihood simply takes the
form
$\log\mathcal{L}(d|\theta)=-\sum_{ij}(f(\theta)_{ij}-d_{ij})^{2}/2\sigma_{ij}^{2}$,
where the indices $i$ and $j$ represent the $j$-th data point in the $i$-th
redshift bin, $d$ is the data to be fit, $f(\theta)$ denotes the data
simulated within MCMC of parameter $\theta$ and $\sigma$ is the standard
deviation of the data. In this work, we pre-calculate the $\sigma$ by forward
simulation of the fiducial model 1000 times, each time varying the initial
conditions and noise realizations. According to our test, the $\sigma$ starts
to converge after 100 realizations, and we only exploit the MFs
$\\{V_{1},V_{2},V_{3}\\}$ in calculating the likelihood, provided that their
changes caused by varying reionization parameters are more notable than
$V_{0}$.

Two fiducial models are considered, the faint galaxies model ($\zeta=30$,
$\log_{10}T_{\rm vir}=4.70$) and the bright galaxies model ($\zeta=200$,
$\log_{10}T_{\rm vir}=5.48$) as in Greig & Mesinger (2017). The parameters to
be constrained are $\zeta$ and $T_{\rm vir}$. We impose a flat prior for
$20\leq\zeta\leq 350$ and for $4\leq\rm\log_{10}T_{\rm vir}\leq 6$. For mock
observations, we again focus on the SKA-Low configuration, and also explore
HERA-like scenarios for an extreme comparison, with the exact specifications
listed in Table 1.

We simplify the fitting by using comoving boxes at $z=7$ and $z=8$ instead of
using lightcones. The comoving boxes with $512\times 512\times 128\,{\rm
Mpc^{3}}$ boxes at $z=7$ and $z=8$ represent an approximation of splitting the
lightcones between $z=7$ and $z=8$ into two sub-bands. As discussed in Section
3.2, the MFs between $z=7$ and $z=8$ are most sensitive to the EoR parameters,
so we expect that most constraining power comes from the $z=7-8$ sub-bands.
Our 21CMMC fitting uses 128 walkers with 100 iterations each, and the initial
conditions for all MC samples are fixed. To achieve a balance between
computational efficiency and inclusion of the cosmic variance, in each
scenario considered in this section, we conducted ten independent runs with
different initial conditions, resulting in $\sim 10^{5}$ samples in total.

We further investigate the parameter fitting with the 21 cm power spectrum
(PS) in the same settings, providing a benchmark for the EoR information
contained in MFs. The sensitivities for HERA and SKA-Low are calculated via
the 21cmSense code with the lower and upper limits of $k$ corresponding to the
box length and the resolution of the assumed telescope, respectively. We
account for three effects that affect the sensitivity: the thermal noise of
the telescope, the sampling variance of the mock observation, and the variance
introduced by the different initial conditions in 21cmFAST (i.e. theoretical
systematics). Joint constraints are also explored by simply combining the log-
likelihoods, to verify if the MFs provide any information beyond the 2-point
statistics.

###  4.3 Results with Telescope Noise

We start by exploring the scenario with thermal noise but no foreground
contamination, taking the fiducial parameter values of the faint galaxies and
bright galaxies models. According to the discussion in Section 4.1, the former
is expected for tighter constraints on reionization parameters than the
latter. Our parameter constraints for SKA-Low are shown in Figure 6 and listed
in Tables 2 and 3. We find that, with SKA-Low observations using MFs, both
$\zeta$ and $\rm T_{\rm vir}$ can be constrained with precise measurement
errors and are consistent with the fiducial values. For faint galaxies model
with SKA-Low, the constraining accuracies of $\zeta$ and $\log_{10}T_{\rm
vir}$ are $(1.23\%^{+2.30\%}_{-2.33\%},0.13\%^{+0.30\%}_{-0.30\%})$ from MFs,
and $(-2.13\%^{+1.70\%}_{-1.40\%},-0.30{\%}^{+0.19\%}_{-0.17\%})$ from PS. For
bright galaxies model with SKA-Low, the constraining accuracies of $\zeta$ and
$\log_{10}T_{\rm vir}$ are
$(7.13\%^{+11.39\%}_{-12.01\%},0.27{\%}^{+0.60\%}_{-0.57\%})$ from MFs, and
$(4.36\%^{+7.17\%}_{-7.51\%},0.20\%^{+0.38\%}_{-0.36\%})$ from PS. Therefore,
in the case of SKA-Low measurements for both faint galaxies and bright
galaxies models, MFs demonstrate comparable (yet slightly worse) constraining
powers to PS.

In comparison, the results with HERA-like configurations are shown in Figure 7
and listed in Tables 2 and 3. For faint galaxies model with HERA-like
observation, the constraining accuracies of $\zeta$ and $\log_{10}T_{\rm vir}$
are $(6.70\%^{+10.90\%}_{-11.07\%},1.11{\%}^{+1.19\%}_{-1.21\%})$ from MFs,
and $(-1.93\%^{+1.80\%}_{-1.87\%},-0.26{\%}^{+0.19\%}_{-0.19\%})$ from PS. For
bright galaxies model with HERA-like observation, the constraining accuracies
of $\zeta$ and $\log_{10}T_{\rm vir}$ are
$(7.18\%^{+23.49\%}_{-20.99\%},0.18{\%}^{+1.28\%}_{-1.08\%})$ from MFs, and
$(5.02\%^{+7.66\%}_{-7.22\%},0.26{\%}^{+0.36\%}_{-0.44\%})$ from PS. We find
that, with HERA observation, the constraints from MFs are much worse than
those from PS, due to the lack of long baselines to resolve small-scale
structures and insufficient sampling of different angular scales.

Figure 8: The bubble size distribution of the reionization field in the
simulation for the faint galaxies model at $z=7$ (solid yellow line) and $z=8$
(solid blue line). The FWHM of the PSF is shown for the HERA-like
configuration (dot-dashed lines) and SKA-Low (dashed lines), respectively.
Figure 9: The posterior distribution of EoR parameters using the SKA-Low
configuration with 1000 h observation and no foreground cut, for the faint
galaxies model. We show the constraint using only $V_{1}$ (red), $V_{2}$
(blue), $V_{3}$ (green), and using all three MFs jointly (black),
respectively.

To further compare the HERA-like configuration with SKA-Low, we calculate the
bubble size distribution (BSD) of the reionization field in our simulation for
the faint galaxies case and compare it with the angular resolution of the
instruments shown in Figure 8. We find that SKA-Low is capable of resolving
major structures in the 21 cm field, while HERA-like configuration can only
resolve the large-scale structures. Furthermore, the lack of $u$-$v$ sampling
from the redundant configuration means that the thermal noise in image space
is highly correlated as shown in Figure 1, thereby further smearing the
structure. The combination of the lack of long baselines and redundant layouts
means that a HERA-like instrument is not ideal for topological studies of EoR.

In addition, to make a diagnostic analysis regarding which MF is more
sensitive to reionization parameters than others, we perform the posterior
distribution using the information of the MF $V_{1}$, $V_{2}$, and $V_{3}$,
respectively, shown in Figure 9. We find that the inference result from each
separate MF is comparable, and the joint result with all three MFs is much
tighter than with any of each MF.

In conclusion, for SKA-Low observations in the case of no foreground
contamination, the MFs can be exploited as summary statistics to constrain
reionization parameters with accuracies of the same order of magnitude as the
power spectrum.

Table 2: Bayesian inference with MFs and/or PS for the faint galaxies model. Three scenarios are considered here — SKA-Low observation with thermal noise, HERA-like observation with thermal noise, and SKA-Low observation with both noise and foreground cut (“FG cut”). For each scenario, we show the inference results from the MFs and from the PS, respectively, and for the SKA (noise+FG cut) scenario, the result of joint inference from both MFs and PS. |  | SKA(Noise) | HERA-like(Noise) | SKA(Noise+FG Cut)
---|---|---|---|---
Parameter | True value | MFs | PS | MFs | PS | MFs | PS | Joint
$\log_{10}\left(T_{\rm vir}/{\rm K}\right)$ | 4.700 | $4.706^{+0.014}_{-0.014}$ | $4.686^{+0.009}_{-0.008}$ | $4.752^{+0.056}_{-0.057}$ | $4.688^{+0.009}_{-0.009}$ | $4.718^{+0.022}_{-0.024}$ | $4.689^{+0.010}_{-0.009}$ | $4.699^{+0.008}_{-0.009}$
$\zeta$ | 30 | $30.37^{+0.69}_{-0.70}$ | $29.36^{+0.51}_{-0.42}$ | $32.01^{+3.27}_{-3.32}$ | $29.42^{+0.56}_{-0.54}$ | $31.14^{+1.06}_{-1.16}$ | $29.46^{+0.55}_{-0.54}$ | $30.05^{+0.46}_{-0.49}$
Table 3: Same as Table 2 but for the bright galaxies model. |  | SKA(Noise) | HERA-like(Noise) | SKA(Noise+FG Cut)
---|---|---|---|---
Parameter | True value | MFs | PS | MFs | PS | MFs | PS | Joint
$\log_{10}\left(T_{\rm vir}/{\rm K}\right)$ | 5.480 | $5.495^{+0.033}_{-0.031}$ | $5.491^{+0.021}_{-0.020}$ | $5.481^{+0.070}_{-0.059}$ | $5.494^{+0.020}_{-0.024}$ | $5.495^{+0.044}_{-0.043}$ | $5.475^{+0.028}_{-0.042}$ | $5.477^{+0.024}_{-0.027}$
$\zeta$ | 200 | $214.26^{+22.78}_{-24.02}$ | $208.72^{+14.34}_{-15.02}$ | $214.36^{+46.97}_{-41.97}$ | $210.40^{+15.33}_{-14.43}$ | $214.65^{+30.64}_{-30.33}$ | $198.32^{+18.88}_{-26.41}$ | $199.83^{+16.14}_{-17.84}$

###  4.4 Results with Noise and Foreground Avoidance

Figure 10: Same as Figure 6 using the SKA-Low configuration but the foreground
cut is applied herein. We also show the joint constraint (green region).

Now we test the robustness of the constraining power of MFs in the presence of
the foregrounds, following the foreground avoidance recipe described in
Section 2.2.3. In Figure 10, we show the results with the foreground avoidance
technique. Here we only consider SKA-Low, because, for the HERA mock
observation, the information along LoS is almost lost, as we discussed in
Section 3.3. We adopt the foreground wedge slope $m=0.7$ herein, the
pessimistic scenario discussed in Sections 2.2.3 and 3.3.

Assuming the foreground avoidance, for faint galaxies model with SKA-Low
observation, the constraining accuracies of $\zeta$ and $\log_{10}T_{\rm vir}$
are $(3.80\%^{+3.53\%}_{-3.87\%},-0.38{\%}^{+0.47\%}_{-0.51\%})$ from MFs, and
$(-1.80\%^{+1.83\%}_{-1.80\%},-0.23\%^{+0.21\%}_{-0.18\%})$ from PS. For
bright galaxies model with SKA-Low observation, the constraining accuracies of
$\zeta$ and $\log_{10}T_{\rm vir}$ are
$(7.33\%^{+15.32\%}_{-15.17\%},0.27\%^{+0.80\%}_{-0.78\%})$ from MFs, and
$(-0.84\%^{+9.44\%}_{-13.21\%},-0.09\%^{+0.51\%}_{-0.77\%})$ from PS. This
indicates that after removing the information underneath the foreground wedge,
the constraining powers of both MFs and PS are degraded due to loss of
information. Comparing the change in the error bars of the parameter fittings,
MFs lose more information in the faint galaxies case than PS. On the other
hand, in the bright galaxies case, the fractional changes in error bars are
similar for MFs and PS. The larger fractional decrease in constraining power
for MFs in the faint galaxies model is because the MFs are more sensitive to
the topology of the reionization process, which has more impact in the faint
galaxies model that generates a larger population of small ionized bubbles.
Nevertheless, while the fractional signal loss is larger in the faint galaxies
case, the absolute constraining power of MFs is stronger in the faint galaxies
case than in the bright galaxies case.

We find that results for both PS and MFs show slight bi-modality in the 1D
posterior distribution, which is also reported in Watkinson et al. (2022). It
hints towards the fact that the fitting is susceptible to parameter
degeneracy. As MFs and PS probe different information about the reionization
process, we expect that the joint analysis can break the degeneracy and give
tighter constraints, which we will explore in Section 4.5.

In conclusion, we find that the MFs are robust against reasonable levels of
thermal noise and foreground cut and can be used to constrain EoR parameters.
For SKA-Low observation, MFs have comparable constraining powers to PS, making
them a set of promising summary statistics beyond the two-point correlation.

###  4.5 Joint Constraints using MFs and PS

Now, we further explore the constraining power of MFs, by studying the joint
analysis using MFs and PS in the parameter fitting. The 2D posterior is shown
in Figure 10 with the fitting results listed in the rightmost columns of
Tables 2 and 3.

The joint analysis gives tighter constraints for both bright galaxies and
faint galaxies cases than the constraints using MFs or PS alone. In
particular, for the faint galaxies case, compared to the results from only PS,
joint fittings improve about $\sim 10\%$ accuracies to both $\zeta$ and
$\log_{10}T_{\rm vir}$ measurements. For the bright galaxies case, the
improvements on the parameter accuracies increase to $\sim 20\%-30\%$.

Furthermore, as shown in Figure 10, the joint constraints also break the
parameter degeneracy. In the faint galaxies case, the 1$\sigma$ credible
region from the MFs and that from PS are slightly misaligned, with the former
significantly larger. As a result, the best fit of the joint analysis is more
precise about the true values, with the joint credible region further shrunk
from the PS credible region. In the bright galaxies case, the credible region
from PS finds a local (but not global) likelihood maximum at the lower values
of the reionization parameters than the true values, while the credible region
from the MFs is very elongated. In comparison, the joint analysis has only one
maximum of likelihood, and the joint credible region is much smaller than that
of the MFs or PS. In general, regarding the improvement of joint analysis upon
PS, we find that the improvement is more significant in the bright galaxies
case than in the faint galaxies case.

We conclude from these results that the MFs contain complementary information
about the EoR parameters to the PS.

###  4.6 Caveats

In this subsection, we discuss several points that may affect the results.

First of all, the morphological structure in the 21 cm brightness temperature
field may be oversimplified in this work because the field is simulated using
21cmFAST which is based on the excursion set model of reionization. Realistic
astrophysical processes, e.g. inclusion of inhomogeneous recombination (e.g.
Mao et al. 2020), might result in a more sophisticated morphological
structure, which might enhance the constraining power of MFs.

In addition, several observational effects such as redshift-space distortion
(e.g. Mao et al. 2012) and lightcone effect (e.g. Datta et al. 2012), which
are neglected in this work, can introduce large anisotropic effects in the 21
cm signal. Since the foreground wedge affects the sampling of the anisotropic
modes along the LoS, how MFs are sensitive to reionization parameters is
subject to further tests.

We also implicitly assume in this work that the foreground wedge corresponds
to the angular scale that is the instrument FoV of the primary beam and
neglects the effect of sidelobes. This assumption requires sufficient
mitigation of the wide-field effect in 21 cm observations. The mitigation of
the wide-field effect is performed in data calibration and imaging, which is
beyond the scope of this paper.

Finally, our posterior inference adopts the standard MCMC algorithm with a
simple Gaussian likelihood and, in particular, neglects the correlation
between any threshold $\delta T_{b}/\sigma_{T_{b}}$ bins and between any
redshift bins in the MFs. This treatment is sufficient for our forecast study,
but it will be necessary for future realistic analysis to include the complete
modeling of the full covariance matrix in the MCMC analysis (Prelogović &
Mesinger, 2023) or perform the simulation-based inference (SBI; Zhao et al.,
2022a, b) using the MFs.

##  5 Summary

In this paper, we explore the viability of using MFs as summary statistics for
21 cm observations during the EoR.

First, we test the robustness of MFs in the presence of observational effects
by generating mock observations that account for the effects of thermal noise,
synthesized beam, and foreground avoidance. The mock observations are
simulated using the configurations of SKA-Low, as well as a redundant array
layout similar to HERA.

We find that the synthesized beam of radio interferometers smooths the 21 cm
brightness temperature field, altering the structure of the MFs. The Gaussian
thermal noise obscures the underlying structure of reionization. Both effects
result in the loss of information in the non-Gaussian reionization process.

For redundant array configurations, the lack of long baselines, which results
in insufficient sampling in different angular scales, takes out the structures
of the MFs in the observed 21 cm signal. As a result, HERA-like configurations
are not ideal for constraining EoR models.

For the SKA-Low array, on the other hand, sufficient sensitivity and baseline
coverage help retain the information in the MFs. Assuming a deep total
integration time of 1000 hours, the topology of reionization can be seen in
the observed signal with its clear evolution, and as such, MFs can be used to
distinguish different stages of reionization as well as different models.

It is expected that foreground avoidance causes further loss of information in
the MFs, making them more Gaussian. However, assuming that the foreground
wedge corresponds to the primary beam of the instrument, we find that there is
still enough information in the observation window for MFs to capture the
topology of EoR.

Next, we forecast the prospects of constraining EoR parameters with MFs. For
this purpose, we build an MCMC parameter-fitting framework that uses forward
modeling. We find that for deep observations using SKA-Low, MFs can be used to
provide unbiased estimations of the reionization parameters, with the forecast
sensitivity comparable to the PS. For the HERA-like configurations,
nevertheless, the constraint accuracies of the reionization parameters using
MFs are about ten times worse than using the PS.

When the foreground cut is applied, the constraint accuracies of the
reionization parameters with the MFs are degraded, as we demonstrate with SKA-
Low observations. The presence of foregrounds leads to more information loss
in the MFs compared to the PS because the MFs are more sensitive to the
topology of the reionization fields.

Jointly fitting the parameters using both MFs and PS, we find that the joint
constraints are tighter than each separate fitting, which implies that the MFs
contain complementary information to the two-point statistics. The
improvements in the constraint accuracies depend on the reionization model.
Compared to the accuracy using only PS, the accuracies in the reionization
parameters are improved with $\sim 30\%$ ($10\%$) in the bright galaxies model
(faint galaxies model).

Our investigation highlights the prospect of constraining the reionization
parameters by exploiting the MFs as complementary summary statistics to the
PS. With the upcoming SKA-Low that will generate high-fidelity 21 cm images
with a large dynamical range, the full topological information in the 21 cm
image cube will be crucial in understanding the reionization and the driving
force behind it.

## acknowledgments

This work is supported by the National SKA Program of China (grant No.
2020SKA0110401), NSFC (grant No. 11821303, 12361141814), and the Chinese
Academy of Sciences grant ZDKYYQ20200008. ZC’s research is funded by a UKRI
Future Leaders Fellowship [grant MR/X005399/1; PI: Alkistis Pourtsidou]. KD
thanks Xiaosheng Zhao, Ce Sui, Bohua Li, and Meng Zhou for inspiring
discussions and help. We acknowledge the Tsinghua Astrophysics High-
Performance Computing platform at Tsinghua University for providing
computational and data storage resources that have contributed to the research
results reported within this paper.

## References

  * Abdurashidova et al. (2022) Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221, doi: 10.3847/1538-4357/ac1c78
  * Adams et al. (2023) Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755, doi: 10.1093/mnras/stac3347
  * Adshead et al. (2013) Adshead, P., Hu, W., & Miranda, V. 2013, Phys. Rev. D, 88, 023507, doi: 10.1103/PhysRevD.88.023507
  * Appleby et al. (2022) Appleby, S., Park, C., Pranav, P., et al. 2022, ApJ, 928, 108, doi: 10.3847/1538-4357/ac562a
  * Austin et al. (2023) Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJ, 952, L7, doi: 10.3847/2041-8213/ace18d
  * Bag et al. (2019) Bag, S., Mondal, R., Sarkar, P., et al. 2019, Monthly Notices of the Royal Astronomical Society, 485, 2235, doi: 10.1093/mnras/stz532
  * Bharadwaj & Ali (2004) Bharadwaj, S., & Ali, S. S. 2004, MNRAS, 352, 142, doi: 10.1111/j.1365-2966.2004.07907.x
  * Bianco et al. (2024) Bianco, M., Giri, S. K., Prelogović, D., et al. 2024, MNRAS, doi: 10.1093/mnras/stae257
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145
  * Braun et al. (2019) Braun, R., Bonaldi, A., Bourke, T., Keane, E., & Wagg, J. 2019, arXiv e-prints, arXiv:1912.12699, doi: 10.48550/arXiv.1912.12699
  * Buchert et al. (2017) Buchert, T., France, M. J., & Steiner, F. 2017, Classical and Quantum Gravity, 34, 094002, doi: 10.1088/1361-6382/aa5ce2
  * Busch et al. (2020) Busch, P., Eide, M. B., Ciardi, B., & Kakiichi, K. 2020, MNRAS, 498, 4533, doi: 10.1093/mnras/staa2599
  * Carron (2011) Carron, J. 2011, ApJ, 738, 86, doi: 10.1088/0004-637X/738/1/86
  * Castellano et al. (2022) Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJ, 938, L15, doi: 10.3847/2041-8213/ac94d0
  * Chen et al. (2023) Chen, Z., Chapman, E., Wolz, L., & Mazumder, A. 2023, MNRAS, 524, 3724, doi: 10.1093/mnras/stad2102
  * Chen et al. (2019) Chen, Z., Xu, Y., Wang, Y., & Chen, X. 2019, ApJ, 885, 23, doi: 10.3847/1538-4357/ab43e6
  * Cheng et al. (2018) Cheng, C., Parsons, A. R., Kolopanis, M., et al. 2018, ApJ, 868, 26, doi: 10.3847/1538-4357/aae833
  * Ciardi et al. (2003) Ciardi, B., Ferrara, A., & White, S. D. M. 2003, MNRAS, 344, L7, doi: 10.1046/j.1365-8711.2003.06976.x
  * Condon & Ransom (2016) Condon, J. J., & Ransom, S. M. 2016, Essential Radio Astronomy
  * Dasgupta et al. (2023) Dasgupta, S., Pal, S. K., Bag, S., et al. 2023, J. Cosmology Astropart. Phys, 2023, 014, doi: 10.1088/1475-7516/2023/05/014
  * Datta et al. (2012) Datta, K. K., Mellema, G., Mao, Y., et al. 2012, MNRAS, 424, 1877, doi: 10.1111/j.1365-2966.2012.21293.x
  * de Lera Acedo et al. (2020) de Lera Acedo, E., Pienaar, H., Razavi Ghods, N., et al. 2020, arXiv e-prints, arXiv:2003.12744, doi: 10.48550/arXiv.2003.12744
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001, doi: 10.1088/1538-3873/129/974/045001
  * Dillon et al. (2020) Dillon, J. S., Lee, M., Ali, Z. S., et al. 2020, MNRAS, 499, 5840, doi: 10.1093/mnras/staa3001
  * Dodelson & Schmidt (2020) Dodelson, S., & Schmidt, F. 2020, Modern Cosmology, doi: 10.1016/C2017-0-01943-2
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13, doi: 10.3847/2041-8213/acade4
  * Furlanetto & Oh (2016) Furlanetto, S. R., & Oh, S. P. 2016, MNRAS, 457, 1813, doi: 10.1093/mnras/stw104
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181, doi: 10.1016/j.physrep.2006.08.002
  * Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1, doi: 10.1086/423025
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Georgiev et al. (2022) Georgiev, I., Mellema, G., Giri, S. K., & Mondal, R. 2022, MNRAS, 513, 5109, doi: 10.1093/mnras/stac1230
  * Giri et al. (2020) Giri, S., Mellema, G., & Jensen, H. 2020, JOSS, 5, 2363, doi: 10.21105/joss.02363
  * Giri & Mellema (2021) Giri, S. K., & Mellema, G. 2021, MNRAS, 505, 1863, doi: 10.1093/mnras/stab1320
  * Giri et al. (2019a) Giri, S. K., Mellema, G., Aldheimer, T., Dixon, K. L., & Iliev, I. T. 2019a, MNRAS, 489, 1590, doi: 10.1093/mnras/stz2224
  * Giri et al. (2019b) —. 2019b, MNRAS, 489, 1590, doi: 10.1093/mnras/stz2224
  * Giri et al. (2018a) Giri, S. K., Mellema, G., Dixon, K. L., & Iliev, I. T. 2018a, MNRAS, 473, 2949, doi: 10.1093/mnras/stx2539
  * Giri et al. (2018b) Giri, S. K., Mellema, G., & Ghara, R. 2018b, MNRAS, 479, 5596, doi: 10.1093/mnras/sty1786
  * Gleser et al. (2006) Gleser, L., Nusser, A., Ciardi, B., & Desjacques, V. 2006, MNRAS, 370, 1329, doi: 10.1111/j.1365-2966.2006.10556.x
  * Gnedin & Ostriker (1997) Gnedin, N. Y., & Ostriker, J. P. 1997, ApJ, 486, 581, doi: 10.1086/304548
  * Gnedin & Shaver (2004) Gnedin, N. Y., & Shaver, P. A. 2004, ApJ, 608, 611, doi: 10.1086/420735
  * Greig & Mesinger (2015) Greig, B., & Mesinger, A. 2015, MNRAS, 449, 4246, doi: 10.1093/mnras/stv571
  * Greig & Mesinger (2017) —. 2017, MNRAS, 472, 2651, doi: 10.1093/mnras/stx2118
  * Greig & Mesinger (2018) —. 2018, MNRAS, 477, 3217, doi: 10.1093/mnras/sty796
  * Hadwiger (1959) Hadwiger, H. 1959, MatZe, 71, 124, doi: 10.1007/BF01181393
  * Harikane et al. (2023a) Harikane, Y., Nakajima, K., Ouchi, M., et al. 2023a, arXiv e-prints, arXiv:2304.06658, doi: 10.48550/arXiv.2304.06658
  * Harikane et al. (2022) Harikane, Y., Ono, Y., Ouchi, M., et al. 2022, ApJS, 259, 20, doi: 10.3847/1538-4365/ac3dfc
  * Harikane et al. (2023b) Harikane, Y., Ouchi, M., Oguri, M., et al. 2023b, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9
  * Hassan et al. (2019) Hassan, S., Liu, A., Kohn, S., & La Plante, P. 2019, MNRAS, 483, 2524, doi: 10.1093/mnras/sty3282
  * Hellwig et al. (1970) Hellwig, H., Vessot, R. F. C., Levine, M. W., et al. 1970, IEEE Transactions on Instrumentation Measurement, 19, 200, doi: 10.1109/TIM.1970.4313902
  * Hikage et al. (2006) Hikage, C., Komatsu, E., & Matsubara, T. 2006, ApJ, 653, 11, doi: 10.1086/508653
  * Hoffmann et al. (2019) Hoffmann, K., Mao, Y., Xu, J., Mo, H., & Wandelt, B. D. 2019, MNRAS, 487, 3050, doi: 10.1093/mnras/stz1472
  * Hothi et al. (2021) Hothi, I., Chapman, E., Pritchard, J. R., et al. 2021, MNRAS, 500, 2264, doi: 10.1093/mnras/staa3446
  * Hutter et al. (2020) Hutter, A., Watkinson, C. A., Seiler, J., et al. 2020, MNRAS, 492, 653, doi: 10.1093/mnras/stz3139
  * Kakiichi et al. (2017) Kakiichi, K., Majumdar, S., Mellema, G., et al. 2017, MNRAS, 471, 1936, doi: 10.1093/mnras/stx1568
  * Kapahtia et al. (2019) Kapahtia, A., Chingangbam, P., & Appleby, S. 2019, J. Cosmology Astropart. Phys, 2019, 053, doi: 10.1088/1475-7516/2019/09/053
  * Kapahtia et al. (2021) Kapahtia, A., Chingangbam, P., Ghara, R., Appleby, S., & Choudhury, T. R. 2021, J. Cosmology Astropart. Phys, 2021, 026, doi: 10.1088/1475-7516/2021/05/026
  * Koenderink (1984) Koenderink, J. J. 1984, Biological Cybernetics, 50, 363, doi: 10.1007/BF00336961
  * Lanman & Pober (2019) Lanman, A. E., & Pober, J. C. 2019, MNRAS, 487, 5840, doi: 10.1093/mnras/stz1639
  * Lee et al. (2008) Lee, K.-G., Cen, R., Gott, J. Richard, I., & Trac, H. 2008, ApJ, 675, 8, doi: 10.1086/525520
  * Lewis (2011) Lewis, A. 2011, J. Cosmology Astropart. Phys, 2011, 026, doi: 10.1088/1475-7516/2011/10/026
  * Li et al. (2019) Li, B., Chen, Z., Cai, Y.-F., & Mao, Y. 2019, MNRAS, 487, 5564, doi: 10.1093/mnras/stz1619
  * Liu et al. (2014) Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023018, doi: 10.1103/PhysRevD.90.023018
  * Liu & Shaw (2020) Liu, A., & Shaw, J. R. 2020, PASP, 132, 062001, doi: 10.1088/1538-3873/ab5bfd
  * Lynch et al. (2021) Lynch, C. R., Galvin, T. J., Line, J. L. B., et al. 2021, PASA, 38, e057, doi: 10.1017/pasa.2021.50
  * Mao et al. (2020) Mao, Y., Koda, J., Shapiro, P. R., et al. 2020, MNRAS, 491, 1600, doi: 10.1093/mnras/stz2986
  * Mao et al. (2012) Mao, Y., Shapiro, P. R., Mellema, G., et al. 2012, MNRAS, 422, 926, doi: 10.1111/j.1365-2966.2012.20471.x
  * Mellema et al. (2006) Mellema, G., Iliev, I. T., Pen, U.-L., & Shapiro, P. R. 2006, MNRAS, 372, 679, doi: 10.1111/j.1365-2966.2006.10919.x
  * Mellema et al. (2013) Mellema, G., Koopmans, L. V. E., Abdalla, F. A., et al. 2013, Experimental Astronomy, 36, 235, doi: 10.1007/s10686-013-9334-5
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955, doi: 10.1111/j.1365-2966.2010.17731.x
  * Minkowski (1903) Minkowski, H. 1903, MatAn, 57, 447, doi: 10.1007/BF01445180
  * Mondal et al. (2016) Mondal, R., Bharadwaj, S., & Majumdar, S. 2016, MNRAS, 456, 1936, doi: 10.1093/mnras/stv2772
  * Morales & Wyithe (2010) Morales, M. F., & Wyithe, J. S. B. 2010, ARA&A, 48, 127, doi: 10.1146/annurev-astro-081309-130936
  * Mort et al. (2010) Mort, B. J., Dulwich, F., Salvini, S., Adami, K. Z., & Jones, M. E. 2010, in 2010 IEEE International Symposium on Phased Array Systems and Technology, 690–694, doi: 10.1109/ARRAY.2010.5613289
  * Muñoz et al. (2022) Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657, doi: 10.1093/mnras/stac185
  * Murray et al. (2020) Murray, S., Greig, B., Mesinger, A., et al. 2020, JOSS, 5, 2582, doi: 10.21105/joss.02582
  * Murray & Trott (2018) Murray, S. G., & Trott, C. M. 2018, ApJ, 869, 25, doi: 10.3847/1538-4357/aaebfa
  * Naidu et al. (2022) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14, doi: 10.3847/2041-8213/ac9b22
  * Novikov et al. (1999) Novikov, D., Feldman, H. A., & Shandarin, S. F. 1999, International Journal of Modern Physics D, 8, 291, doi: 10.1142/S0218271899000225
  * Pagano & Liu (2020) Pagano, M., & Liu, A. 2020, MNRAS, 498, 373, doi: 10.1093/mnras/staa2118
  * Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, AJ, 139, 1468, doi: 10.1088/0004-6256/139/4/1468
  * Pathak et al. (2022) Pathak, A., Bag, S., Dasgupta, S., et al. 2022, J. Cosmology Astropart. Phys, 2022, 027, doi: 10.1088/1475-7516/2022/11/027
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Pober et al. (2013) Pober, J. C., Parsons, A. R., DeBoer, D. R., et al. 2013, AJ, 145, 65, doi: 10.1088/0004-6256/145/3/65
  * Pober et al. (2014) Pober, J. C., Liu, A., Dillon, J. S., et al. 2014, ApJ, 782, 66, doi: 10.1088/0004-637X/782/2/66
  * Pober et al. (2016) Pober, J. C., Hazelton, B. J., Beardsley, A. P., et al. 2016, ApJ, 819, 8, doi: 10.3847/0004-637X/819/1/8
  * Pratten & Munshi (2012) Pratten, G., & Munshi, D. 2012, MNRAS, 423, 3209, doi: 10.1111/j.1365-2966.2012.21103.x
  * Prelogović & Mesinger (2023) Prelogović, D., & Mesinger, A. 2023, MNRAS, 524, 4239, doi: 10.1093/mnras/stad2027
  * Pritchard & Loeb (2012) Pritchard, J. R., & Loeb, A. 2012, Reports on Progress in Physics, 75, 086901, doi: 10.1088/0034-4885/75/8/086901
  * Schmalzing & Buchert (1997) Schmalzing, J., & Buchert, T. 1997, ApJ, 482, L1, doi: 10.1086/310680
  * Scoccimarro (1998) Scoccimarro, R. 1998, MNRAS, 299, 1097, doi: 10.1046/j.1365-8711.1998.01845.x
  * Shaw et al. (2019) Shaw, A. K., Bharadwaj, S., & Mondal, R. 2019, MNRAS, 487, 4951, doi: 10.1093/mnras/stz1561
  * Shaw et al. (2020) —. 2020, MNRAS, 498, 1480, doi: 10.1093/mnras/staa2090
  * Shimabukuro et al. (2022) Shimabukuro, H., Mao, Y., & Tan, J. 2022, Research in Astronomy and Astrophysics, 22, 035027, doi: 10.1088/1674-4527/ac4ca3
  * Shimabukuro et al. (2015) Shimabukuro, H., Yoshiura, S., Takahashi, K., Yokoyama, S., & Ichiki, K. 2015, MNRAS, 451, 467, doi: 10.1093/mnras/stv965
  * Shimabukuro et al. (2017) —. 2017, MNRAS, 468, 1542, doi: 10.1093/mnras/stx530
  * Smirnov (2011) Smirnov, O. M. 2011, A&A, 527, A106, doi: 10.1051/0004-6361/201016082
  * Thyagarajan et al. (2015) Thyagarajan, N., Jacobs, D. C., Bowman, J. D., et al. 2015, ApJ, 804, 14, doi: 10.1088/0004-637X/804/1/14
  * Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007, doi: 10.1017/pasa.2012.007
  * Tiwari et al. (2022) Tiwari, H., Shaw, A. K., Majumdar, S., Kamran, M., & Choudhury, M. 2022, J. Cosmology Astropart. Phys, 2022, 045, doi: 10.1088/1475-7516/2022/04/045
  * Trott et al. (2020) Trott, C. M., Jordan, C. H., Midgley, S., et al. 2020, MNRAS, 493, 4711, doi: 10.1093/mnras/staa414
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2, doi: 10.1051/0004-6361/201220873
  * Wang et al. (2015) Wang, Y., Park, C., Xu, Y., Chen, X., & Kim, J. 2015, ApJ, 814, 6, doi: 10.1088/0004-637X/814/1/6
  * Watkinson et al. (2022) Watkinson, C. A., Greig, B., & Mesinger, A. 2022, MNRAS, 510, 3838, doi: 10.1093/mnras/stab3706
  * Watkinson & Pritchard (2015) Watkinson, C. A., & Pritchard, J. R. 2015, MNRAS, 454, 1416, doi: 10.1093/mnras/stv2010
  * Xu et al. (2019) Xu, W., Xu, Y., Yue, B., et al. 2019, Mon. Not. Roy. Astron. Soc., 490, 5739, doi: 10.1093/mnras/stz2926
  * Xu et al. (2017) Xu, Y., Yue, B., & Chen, X. 2017, ApJ, 844, 117, doi: 10.3847/1538-4357/aa7b71
  * Xu et al. (2014) Xu, Y., Yue, B., Su, M., Fan, Z., & Chen, X. 2014, Astrophys. J., 781, 97, doi: 10.1088/0004-637X/781/2/97
  * Yatawatta et al. (2013) Yatawatta, S., de Bruyn, A. G., Brentjens, M. A., et al. 2013, A&A, 550, A136, doi: 10.1051/0004-6361/201220874
  * Yoshiura et al. (2017) Yoshiura, S., Shimabukuro, H., Takahashi, K., & Matsubara, T. 2017, MNRAS, 465, 394, doi: 10.1093/mnras/stw2701
  * Yoshiura et al. (2015) Yoshiura, S., Shimabukuro, H., Takahashi, K., et al. 2015, MNRAS, 451, 266, doi: 10.1093/mnras/stv855
  * Zhao et al. (2022a) Zhao, X., Mao, Y., Cheng, C., & Wandelt, B. D. 2022a, ApJ, 926, 151, doi: 10.3847/1538-4357/ac457d
  * Zhao et al. (2022b) Zhao, X., Mao, Y., & Wandelt, B. D. 2022b, ApJ, 933, 236, doi: 10.3847/1538-4357/ac778e