# The Impact of ionization Morphology and X-ray Heating on the Cosmological
21cm Skew Spectrum

J. H. Cook1,2, S. Balu2,3, B. Greig2,3,4, C. M. Trott1,2, J. L. B. Line1,2, Y.
Qin2,3, J. S. B. Wyithe2,4
1International Centre for Radio Astronomy Research, Curtin University, Perth,
Australia
2ARC Centre of Excellence for All Sky Astrophysics in 3D (ASTRO 3D)
3School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
4Research School of Astronomy & Astrophysics, Australian National University,
Canberra, ACT 2611, Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The cosmological 21cm signal offers a potential probe of the early Universe
and the first ionizing sources. Current experiments probe the spatially-
dependent variance (Gaussianity) of the signal through the power spectrum
(PS). The signal however is expected to be highly non-Gaussian due to the
complex topology of reionization and X-ray heating. We investigate the non-
Gaussianities of X-ray heating and reionization, by calculating the skew
spectrum (SS) of the 21cm signal using meraxes, which couples a semi-analytic
galaxy population with semi-numerical reionization simulations. The SS is the
cross-spectrum of the quadratic temperature brightness field with itself. We
generate a set of seven simulations from $z=30$ to $z=5$, varying the halo
mass threshold for hosting star-formation, the X-ray luminosity per star-
formation rate, and the minimum X-ray energy escaping host galaxies. We find
the SS is predominantly negative as a function of redshift, transitioning to
positive towards the start of reionization, and peaking during the midpoint of
reionization. We do not see a negative dip in the SS during reionization,
likely due to the specifics of modelling ionization sources. We normalise the
SS by the PS during reionization isolating the non-Gaussianities. We find a
trough ($k\sim\,0.1\,\textrm{Mpc}^{-1}$) and peak
($k\sim\,0.4-1\,\textrm{Mpc}^{-1}$) in the normalised SS during the mid to
late periods of reionization. These correlate to the ionization topology, and
neutral islands in the IGM. We calculate the cosmic variance of the normalised
SS, and find these features are detectable in the absence of foregrounds with
the SKA_LOW.

###### keywords:

cosmology: dark ages, reionization, first stars – methods: statistical

††pubyear: 2024††pagerange: The Impact of ionization Morphology and X-ray
Heating on the Cosmological 21cm Skew Spectrum–A

##  1 Introduction

The cosmological 21cm neutral hydrogen line promises to be an insightful probe
of the first luminous sources and the structure of the Universe during early
cosmic time. The first luminous sources (stars, galaxies, compact objects)
heat and ionize the surrounding intergalactic medium (IGM), through the
cumulative emission of ultraviolet (UV) and X-ray photons (see the following
review papers: Barkana & Loeb, 2001; Morales & Wyithe, 2010; Pritchard & Loeb,
2012; Furlanetto, 2016).111Reionization predominately occurs due to UV
photons, with some contribution from X-ray emission (up to $\sim 10$ percent
Mesinger et al. (2013)).. These ionized bubbles grow and eventually overlap,
culminating in the end of reionization by redshift $\sim 5.3$ (Bosman et al.,
2022). These bubbles encode information about these sources onto the
cosmological 21cm temperature brightness signal (Furlanetto & Oh, 2005). These
luminous sources also heat the neutral hydrogen medium through X-ray emission,
which encodes additional information about these sources (Pritchard &
Furlanetto, 2007; Furlanetto, 2016). The cosmological 21cm signal is measured
relative to the cosmic microwave background (CMB), and can be either in
relative emission or absorption. By measuring the 21cm signal we can construct
the spatial, and line of sight distributions of neutral hydrogen. This will
allow for the properties of the first luminous sources to be inferred through
their influence on the cosmological 21cm signal.

Most of the focus in the 21cm cosmological community has been on measuring
either the one or two-point statistics of the signal. The one point statistic
experiments determine the sky averaged quantities (for example the global mean
temperature). For example: The Shaped Antenna measurement of the background
Radio Spectrum 3 telescope (SARAS3, Nambissan et al., 2021); the Experiment to
Detect the Global EoR Signature (EDGES, Bowman et al., 2018). The two point
statistic experiments are primarily measured by radio interferometers. The
current generation of radio interferometers includes the Murchison Widefield
Array (MWA, Tingay et al., 2013; Wayth et al., 2018); Low-Frequency Array
(LOFAR, van Haarlem, M. P. et al., 2013); Hydrogen Epoch of Reionization Array
(HERA, DeBoer et al., 2017); the New extension in Nancay upgrading LOFAR
(NenuFAR, Zarka et al., 2012).

The two point statistic experiments calculate the PS of the 21cm signal which
is the Fourier transform of the two point correlation function. This measures
the Gaussianity or the variance of the signal as a function of comoving
spatial scale. If the signal is entirely Gaussian this would capture all the
information about the 21cm signal within cosmic variance222Since we cannot
truly measure the ensemble average power spectrum, we can only estimate it
over some volume. Each independent realisation therefore is a random sample of
the true PS with some cosmic variance.. The signal is however expected to be
highly non-Gaussian as it evolves during the Epoch of Heating (EoH) and the
Epoch of reionization (EoR) (Wyithe & Morales, 2007; Lidz et al., 2007). In
the former the non-Gaussianities are driven by the appearance of the first
luminous sources which heat the neutral IGM primarily through X-ray emission
(Furlanetto, 2006). During this period, the strong emission from the first
luminous sources drive above average temperature contrasts relative to the
IGM. Eventually as X-ray heating progresses the medium saturates driving the
non-Gaussianities to the matter density (Watkinson et al., 2018). During the
latter stages of reionization the spin temperature of the neutral hydrogen is
expected to be saturated $(T_{S}\gg T_{\rm{CMB}})$, and therefore the non-
Gaussianities are largely driven by the ionization topology (Hutter et al.,
2019). Analytical estimates of the characteristic size of the ionization
topology (bubbles) around individual luminous sources are $\sim
10\,\textrm{cMpc}$ during the late time period of reionsation (Wyithe & Loeb,
2004; Furlanetto & Oh, 2005; Zahn et al., 2007). However, Lin et al. (2016)
showed that the characteristic size is underestimated and is closer to $\sim
20-100\,\textrm{cMpc}$. Lin et al. (2016) and Giri et al. (2017) demonstrate
the difficulty of determining the characteristic size from the complex 3D
ionization topology during reionization.

Non-Gaussianity has been shown to be important in constraining the 21cm
signal, in particular during reionization, and could be important for
confirming a detection of the 21cm signal (Shimabukuro et al., 2017). Non-
Gaussianity has primarily been investigated by calculating the expected 21cm
bispectrum. The bispectrum is the Fourier transform of the three point
correlation function (Peebles, 1981), like the PS it probes the central third
order moment as a function of spatial scales. The bispectrum offers a
complementary picture of the 21cm signal especially during the EoR (Bharadwaj
& Pandey, 2005; Majumdar et al., 2018). Watkinson et al. (2018) investigated
the non-Gaussianity due to X-ray heating from stellar sources and high mass
X-ray binaries. Hutter et al. (2019) investigated the ionization morphology
and the effects on the non-Gaussianity of the 21cm signal during reionization.
Numerous studies have been conducted on the bispectrum and its sensitivity
during reionization (Bharadwaj & Pandey, 2005; Yoshiura et al., 2015;
Shimabukuro et al., 2016, 2017; Watkinson et al., 2017; Majumdar et al., 2018;
Mondal et al., 2021; Majumdar et al., 2020). Trott et al. (2019) measure the
bispectrum of MWA data, looking at a gridded and non-gridded estimator,
establishing upper limits. Watkinson et al. (2020) looked at the expected
foreground bispectrum, commenting on the detectability of the 21cm bispectrum
in the presence of foreground systematics. More recently Tiwari et al. (2022)
showed that the bispectrum can help constrain reionization parameters.

The bispectrum has low signal to noise relative to the PS, and is
computationally intensive to measure, even with the fast Fourier transform
method of Watkinson et al. (2017). To mitigate the difficulties related to
computation and sensitivity, much focus has been on calculating the
equilateral bispectrum (for example: Bharadwaj & Pandey, 2005; Yoshiura et
al., 2015; Watkinson et al., 2018), as well as the squeezed bispectrum, which
compresses one of the triangle mode sides (Chiang et al., 2014; Mondal et al.,
2021). Other work has investigated the higher order one point statistics of
the simulated 21cm signal, due to X-ray heating and Ly-$\alpha$ coupling (for
example: Watkinson & Pritchard, 2015; Ross et al., 2019). The CMB cosmology
community has investigated alternatives that probe the non-Gaussianity through
the cross spectrum of quadratic temperature fields with the temperature field
(Cooray, 2001). This is called the skew spectrum (SS), and is a collapsed form
of the bispectrum, compressing the information into a pseudo PS as a function
of one wavenumber (Fourier modes ($k$)) (Regan, 2017). Generalised in Szapudi
& Szalay (1997); Munshi et al. (1998) and first used by Cooray (2001), it is
now gaining interest in the 21cm community with the release of Ma & Peng
(2023, hereafter MP23) at the time of writing this paper. Again drawing on the
CMB cosmology community for inspiration, Dai et al. (2020) investigated what
information can be gained by combining the PS and the SS. They found that the
SS in conjunction with the PS offered increased constraints on cosmological
parameters. The SS promises to have better signal to noise than the bispectrum
because it integrates over all bispectrum triangle configurations for a given
Fourier mode $k$. Additionally the SS can be directly compared to the PS
because they can be measured at the same Fourier modes. The quadratic field
cross correlation approach also makes it easy to measure the SS from
simulations without having to first calculate the bispectrum.

In this work we use the updated version of meraxes, which couples a semi-
analytic galaxy formation model with a semi-numerical reionization simulation
to provide a realistic population of galaxies which can interact with the IGM
through a variety of feedback effects. These feedback effects include
supernovae, AGN, and photoheating, along with the infall/accretion of gas
(Mutch et al., 2016; Qin et al., 2017a; Qiu et al., 2019). Balu et al. (2023)
updated meraxes to include X-ray heating and spin temperature evolution for
the semi-analytic galaxy formation model. Additionally in Balu et al. (2023),
the halo merger trees for the $210\,h^{-1}\,\textrm{Mpc}$ simulations were
augmented to include all atomically cooled galaxies out to $z=20$ $(\sim
2\times 10^{7}\,h^{-1}\,\textrm{M}_{\odot})$. We build on the earlier work of
Balu et al. (2023) by performing additional simulations varying the ionization
morphology through changing the minimum mass threshold for galaxies hosting
star formation, the X-ray luminosity and minimum energy threshold for X-rays
escaping their host galaxies. Combined, these simulations enable the
exploration of the cosmic evolution of the 21cm signal using a realistic
population of galaxies from the cosmic dawn down to the completion of the EoR.
These are ideal for studying the non-Gaussianity of the EoH and the EoR using
the SS.

The paper is outlined as follows; in section 2 we define the PS and SS. In
section 3 we briefly describe the simulations performed in this work. Section
4 presents the thermal and ionization history of each simulation as well as
the statistics as a function of redshift. Section 5 presents the PS, the SS
and the normalised SS during the EoR for each simulation. Section 6 discusses
the detectability of the normalised SS for the future SKA_LOW radio
interferometer. We discuss and draw conclusions from the results in Section 7.
The cosmology used in this work is defined by Planck Collaboration et al.
(2021): $h=0.68$, $\Omega_{m}=0.31$, $\Omega_{b}=0.048$,
$\Omega_{\Lambda}=0.69$, $\sigma_{8}=0.81$, and $n_{s}=0.96$. All cosmological
scales are in comoving units.

##  2 Power Spectrum and Skew Spectrum

In this section we review the PS and the SS, and how they are calculated from
simulation volumes.

###  2.1 Power Spectrum

The PS is the Fourier transform of the two point correlation function, and
probes the Gaussianity of a random field as a function of spatial scale $k$
(Peebles, 1981):

$\langle\tilde{\delta}_{T}(\mathbf{k})\tilde{\delta}_{T}(\mathbf{k^{\prime}})\rangle\equiv(2\pi)^{3}\delta_{D}(\mathbf{k}+\mathbf{k^{\prime}})P(k).$
(1)

In Equation 1 the angular brackets $\langle\rangle$ denote the ensemble
average over different realisations of the Universe. The Dirac delta
$\delta_{D}(\mathbf{k}+\mathbf{k^{\prime}})$ restricts the average to
uncorrelated modes, and $P(k)$ is the spherically averaged power spectrum.
$\delta_{T}(\mathbf{x})$ is the brightness temperature (or density) contrast
$\delta_{T}(\mathbf{x})=(T(\mathbf{x})-\bar{T}(\mathbf{x}))/\bar{T}(\mathbf{x})$,
and $\bar{T}(\mathbf{x})$ is the mean temperature.
$\tilde{\delta}_{T}(\mathbf{k})$ is the three dimensional Fourier transform of
$\delta_{T}(\mathbf{x})$, defined by333In practice we perform the Fourier
transform over $\delta T_{b}(x)$ not $\delta_{T}(x)$.:

$\tilde{\delta}_{T}(\mathbf{k})=\frac{V}{N_{\rm{pix}}}\sum\delta_{T}(\mathbf{x})\,e^{-i\mathbf{k}\cdot\mathbf{x}}.$
(2)

We define the dimensionless PS:

$\Delta^{2}_{T}(k,z)\equiv\frac{k^{3}}{(2\pi^{2})V}\,\mskip
1.5mu\overline{\mskip-1.5mu{\delta T}\mskip-1.5mu}\mskip
1.5mu_{b}^{2}\langle\tilde{\delta}_{T}(\mathbf{k},z)\tilde{\delta}^{*}_{T}(\mathbf{k},z)\rangle\>\rm{mK^{2}},$
(3)

where the angular brackets now denote the incoherent average in a spherical
shell of width $\Delta\log{k}=0.173$, and
$\tilde{\delta}^{*}_{T}(\mathbf{k},z)$ is the conjugate transpose of
$\tilde{\delta}_{T}(\mathbf{k},z)$.

###  2.2 Bispectrum and Skew Spectrum

The bispectrum of the 21cm brightness temperature fluctuations is defined as:

$\langle\tilde{\delta}_{T}(\mathbf{k}_{1})\tilde{\delta}_{T}(\mathbf{k}_{2})\tilde{\delta}_{T}(\mathbf{k}_{3})\rangle\equiv(2\pi)^{3}\delta_{D}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})B(k_{1},k_{2},k_{3}),$
(4)

again the angular brackets and the Dirac delta function denote the ensemble
average. For the bispectrum the ensemble average is over all triplet values
that satisfy the closed triangle condition
$\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=0$. Without loss of generality
we let $\mathbf{k}_{1}=\mathbf{k}$, $\mathbf{k}_{2}=\mathbf{q}$, and
$\mathbf{k}_{3}=-\left(\mathbf{k}+\mathbf{q}\right)$. The SS is the integral
of the bispectrum $B(k,q,|\mathbf{k}+\mathbf{q}|)$ over all possible triangle
configurations, for a fixed triangle side $k$:

$S_{\gamma}(k)=\frac{1}{(2\pi)^{3}}\int_{\mathbb{R}}d^{3}q\,B(k,q,|\mathbf{k}+\mathbf{q}|).$
(5)

It can be shown that Equation 5 is equivalent to the cross spectrum of the
mean subtracted squared temperature field, to the temperature field:

$\langle\tilde{\delta}_{T^{2}}(\mathbf{k})\tilde{\delta}_{T}(\mathbf{k^{\prime}})\rangle\equiv(2\pi)^{3}\delta_{D}(\mathbf{k}+\mathbf{k^{\prime}})S_{\gamma}(k),$
(6)

Similarly to the PS definition in Equation 1 the delta function and the
angular brackets denote the ensemble average. We define the Fourier transform
of the squared temperature field $\tilde{\delta}_{T^{2}}(\mathbf{k})$ below:

$\tilde{\delta}_{T^{2}}(\mathbf{k})=\frac{V}{N_{\rm{pix}}}\sum\left(\delta_{T}(\mathbf{x})\right)^{2}\,e^{-i\mathbf{k}\cdot\mathbf{x}}.$
(7)

Similar to Equation 3 we can define the dimensionless SS:

$\Delta^{2}_{T^{2},T}(k,z)\equiv\frac{k^{3}}{(2\pi^{2})V}\,\mskip
1.5mu\overline{\mskip-1.5mu{\delta T}\mskip-1.5mu}\mskip
1.5mu_{b}^{3}\langle\tilde{\delta}_{T^{2}}(\mathbf{k},z)\tilde{\delta}^{*}_{T}(\mathbf{k},z)\rangle.$
(8)

We calculate the SS by performing a three dimensional Fourier transform of the
mean subtracted and squared temperature field. We then take the product of
this with the conjugate of the Fourier transform of the temperature field, and
then average in spherical shells of width $\Delta\log{k}=0.173$. For
consistency and for comparison we use the same bins for calculating the SS and
the PS throughout this work.

##  3 MERAXES

In order to simulate the cosmic 21cm signal, we use the meraxes (Mutch et al.,
2016) semi-analytical galaxy formation and evolution model. In this section,
we give a brief summary of meraxes and refer the reader to other relevant
works for further details.

###  3.1 A Realistic Galaxy Population

We make use of the L210_N4320 dark matter-only simulation of the genesis suite
of $N$-body simulations (Power et al. in prep). L210_N4320 has $4320^{3}$ dark
matter particles in a cubical volume of side length $L=210h^{-1}\,\rm{Mpc}$
achieving a halo mass resolution of $\sim 5\times 10^{8}~{}M_{\odot}$.

The halo merger trees from L210_N4320 were further ‘augmented’ to a halo mass
resolution of $\sim 3\times 10^{7}~{}M_{\odot}$, the atomic cooling limit at
$z=20$, using the Monte-Carlo algorithm code DarkForest (Qiu et al., 2021).
This is achieved by sampling low-mass haloes from a conditional halo mass
function that is based on the extended Press-Schechter theory (Bond et al.,
1991; Bower, 1991; Lacey & Cole, 1993), after modifications to match the
N-body simulations’ halo mass functions (HMFs). These haloes are then
‘grafted’ onto the L210_N4320 merger tree in a manner such that the final
augmented HMFs agree with those from high resolution N-body simulations (see
Figure 2 of Balu et al., 2023). DarkForest also assigns and evolves the
positions of the newly added haloes using the local halo density field and the
linear continuity equation (see Qiu et al., 2021; Balu et al., 2023, for
further details). We therefore effectively create an N-body simulation that
has a statistically complete galaxy population down to redshift $z=20$; we
deploy meraxes on this augmented simulation.

The goal of meraxes is to simulate the growth and evolution of galaxies during
the EoR in a self-consistent manner. This is achieved through detailed and
physically motivated prescriptions for varied astrophysical phenomena such as
radiative cooling of gas, star formation, supernovae (SNe) and active galactic
nuclei feedback, and mergers (Mutch et al., 2016; Qin et al., 2017b; Qin et
al., 2017a; Qiu et al., 2019).

For each simulation snapshot, a dark matter halo increases its baryonic mass
in proportion to the universal cosmic baryonic fraction. This mass is added to
a ‘hot gas’ reservoir from where it can cool down to form a ‘cold gas disk’.
Following a star formation prescription based on the Kennicutt-Schmidt law
(Kennicutt, 1998), stars are created out of this cold disk when a cold mass
threshold is reached. The cadence of our simulation is constructed so that the
longest time-step is $\sim 16$ Myr. Hence newly formed massive stars can go
SNe in the same time-step and less massive stars can survive for a few
snapshots. meraxes therefore, has implementations for both instantaneous and
delayed SNe feedback. The primary impact of SNe is to heat up the cold gas in
a galaxy. SNe therefore move a portion of the cold gas to the hot halo and in
very extreme energetic cases can even remove the gas from the galaxy
altogether.

The amount of stellar mass in a galaxy fixes the amount of ionizing UV and
X-ray photons that it produces. Once the local environment of a galaxy is
ionized, the cooling properties of the IGM are affected. meraxes couples
reionization feedback and galaxy growth by self-consistently modifying the
amount of gas that is accreted onto a galaxy depending on the local UV
background and the local IGM ionization state. The UV escape fraction $f_{\rm
esc}(\leq 1)$ of the galaxies is a power-law in redshift $z$ (also see Section
3.4.1):

$f_{\rm esc}=f_{\rm esc,0}\bigg{(}\frac{1+z}{6}\bigg{)}^{\alpha_{\rm esc}},$
(9)

where $f_{\rm esc,0}$ is the escape fraction normalisation and $\alpha_{\rm
esc}=0.20$ is the power-law index.

In this work, we adopt the same fiducial simulation as Balu et al. (2023),
L210_AUG (hereon labelled Fiducial). This simulation has been calibrated with
respect to the UV luminosity functions and the colour-magnitude relation
within a rigorous Bayesian framework (Qiu et al., 2019), as well as the
stellar mass functions (Balu et al., 2023), at $z\sim 4-10$. Reionization
parameters were tuned such that the reionization history is consistent with
existing measurements of the IGM neutral fraction and the CMB optical depth
(see Balu et al., 2023, in particular Figure 3 and Table 2).

###  3.2 IGM Ionization State

meraxes computes the IGM ionization, following the semi-numerical code
21cmfast (Mesinger et al., 2011), via an excursion-set formalism (Furlanetto
et al., 2004). First, we grid the simulation volume and assign galaxies to the
voxels based on their positions. We subdivide our simulation volume into
$1024^{3}$ cells, corresponding to a cell size of $\sim 0.2~{}h^{-1}$ Mpc. In
spheres of decreasing radii, we compare the number of ionizing photons from
the stellar baryons and the total baryons in the IGM. After accounting for
recombinations that can happen in the densest parts of the IGM, we flag a cell
as ionized when the number of ionizing photons is higher than that of the
neutral baryons.

$N_{\rm b*}(\boldsymbol{x},z|R)N_{\gamma}f_{\rm esc}\geq N_{\rm
atom}(\boldsymbol{x},z|R)(1+\bar{n}_{\rm rec})(1-\bar{x}_{e}),$ (10)

where $N_{\rm b*}(\boldsymbol{x},z|R)$ is the number of stellar baryons in a
sphere of radius $R$ centred at $\boldsymbol{x}$ and redshift $z$,
$N_{\gamma}=4000$ is the number of UV ionizing photons per baryon (Barkana &
Loeb, 2007), $N_{\rm atom}$ is the number of neutral H i in the same volume,
$\bar{n}_{\rm rec}$ is the average number of recombinations in the IGM
(Sobacchi & Mesinger, 2014), and $\bar{x}_{e}$ is the mean electron fraction
accounting for the secondary ionizations caused by X-ray photons. Motivated by
the mean-free path of a typical UV photon in the IGM (Songaila & Cowie, 2010),
we set the maximum of $R=50$ Mpc and decrease it successively down to the size
of a cell.

###  3.3 21cm Signal

The differential brightness temperature of the 21cm emission from a cloud of H
i gas illuminated by CMB radiation of temperature $T_{\rm{CMB}}$ is given by:

$\begin{split}\delta T_{b}(\nu)&=\dfrac{T_{\rm
S}-T_{\rm{CMB}}}{1+z}(1-e^{-\tau_{\nu_{0}}})\\\ &\approx
27x_{\textsc{\text{H\,{i}}}}(1+\delta_{\rm nl})\left(\dfrac{H}{dv_{\rm
r}/dr+H}\right)\left(1-\frac{T_{\rm{CMB}}}{T_{\rm S}}\right)\\\
&\mathrm{\hskip 8.5359pt}\times\left(\frac{1+z}{10}\frac{0.15}{\Omega_{\rm
M}h^{2}}\right)^{1/2}\left(\dfrac{\Omega_{\rm
b}h^{2}}{0.023}\right)\,\mathrm{mK},\end{split}$ (11)

where $T_{\rm S}$ is the IGM spin temperature which determines the energy
level populations of the H i hyperfine states, $\tau_{\nu_{0}}$ is the optical
depth, $\delta_{\rm nl}\equiv\rho/\bar{\rho}-1$ is defined as the evolved
Eulerian density contrast ($\rho$ is the density), $H(z)$ is the Hubble
parameter, $dv_{r}/dr$ is the line-of-sight co-moving velocity gradient, and
$x_{\text{H\,{i}}}$ is the neutral fraction. meraxes sources the density and
the velocity fields from the N-body simulations and creates self-consistent
$T_{\rm S}$ and $x_{\text{H\,{i}}}$ fields.

####  3.3.1 Spin Temperature

As can be seen from the Equation (11), the spin temperature $T_{\rm S}$ of the
IGM plays a major role in the 21cm signal. The level populations of the H i
hyperfine states depend on a number of physical processes in the IGM,
including the amount and the energy of the UV and X-ray photons. $T_{\rm S}$
is related to the UV and X-ray emission via:

$T_{\rm S}^{-1}=\dfrac{T_{\rm CMB}^{-1}+x_{\alpha}T_{\alpha}^{-1}+x_{c}T_{\rm
K}^{-1}}{1+x_{\alpha}+x_{\rm c}},$ (12)

where $x_{\alpha}$ and $x_{\rm c}$ are the Wouthuysen-Field coupling
(Wouthuysen, 1952; Field, 1958) and the collisional coupling coefficients
respectively. $x_{\rm c}$ is computed by taking into account the collisions of
H i atoms amongst themselves as well as with free electrons and protons in the
IGM (Zygelman, 2005; Furlanetto & Furlanetto, 2006). $x_{\alpha}$ depends on
the local Ly $\alpha$ background flux and closely follows the implementation
in Mesinger et al. (2011). $T_{\alpha}$ is the ‘colour’ temperature, $T_{\rm
K}$ is the kinetic temperature of the IGM, and we assume $T_{\alpha}=T_{\rm
K}$ (Field, 1959).

The spin temperature field is therefore very sensitive to the $T_{\rm K}$,
which is impacted by X-ray heating. The evolution of the $T_{\rm K}$ depends
on the angle-averaged X-ray intensity $J(\boldsymbol{x},E,z)$ which is
computed as a function of the position $\boldsymbol{x}$, X-ray photon energy
$E$, and redshift $z$:

$J(\boldsymbol{x},E,z)=\dfrac{(1+z)^{3}}{4\pi}\int_{z}^{\infty}dz^{\prime}\frac{cdt}{dz^{{}^{\prime}}}\epsilon_{\rm
X}e^{-\tau},$ (13)

where we have integrated the comoving X-ray emissivity $\epsilon_{\rm X}$ back
along the light cone, and $e^{-\tau}$ accounts for the probability that an
X-ray photon emitted at redshift $z^{\prime}$ survives till $z$. We compute
$\epsilon_{\rm X}$ as a function of the position $\boldsymbol{x}$, X-ray
photon energy $E_{e}=E(1+z^{\prime})/(1+z)$ at the emitted redshift
$z^{\prime}$:

$\epsilon_{X}(\boldsymbol{x},E_{e},z^{\prime})=\dfrac{L_{X}}{\rm
SFR}\times\rm{SFRD}(\boldsymbol{x},E_{e},z^{\prime}),$ (14)

where $L_{X}/{\rm SFR}$ is the galaxies’ specific X-ray luminosity per unit
star formation rate (SFR), and SFRD is the star formation rate density. We
assume a power-law in X-ray photon energy $E$, $L_{X}/{\rm SFR}\propto
E^{-\alpha_{X}}$, with $\alpha_{X}=1$ which is consistent with observations of
high mass X-ray binaries in the local Universe (Mineo et al., 2012; Fragos et
al., 2013; Pacucci et al., 2014), and is normalised:

$L_{X<2~{}\rm{keV}}/{\rm
SFR}=\int_{E_{0}}^{2~{}\rm{keV}}dE_{e}~{}L_{X}/{SFR},$ (15)

where $L_{X<2~{}\rm{keV}}/{\rm SFR}$ is the soft-band X-ray luminosity per SFR
in units of (erg s-1 M${}_{\odot}^{-1}$ yr), and $E_{0}$ fixes the minimum
energy for an X-ray photon so that it is not absorbed within the galaxy.

###  3.4 Simulations

Table 1: Astrophysical parameter summary for the seven meraxes simulations. See text for details. Name | Minimum Halo Mass ($M_{\rm thresh}$) | $f_{\rm esc,0}$ | X-ray Luminosity ($L_{X<2\,\rm{keV}}/\rm{SFR}$) | X-ray threshold ($E_{0}$) | Comments
---|---|---|---|---|---
| $[M_{\sun}]$ |  | $[\rm{erg\,s^{-1}\,M^{-1}_{\sun}\,yr}]$ | $[\rm{keV}]$ |
Fiducial | $3\times 10^{7}$ | 0.14 | $3.16\times 10^{40}$ | $0.5$ | Fiducial simulation
mid_$M_{h}$ | $10^{9}$ | 0.25 | $3.16\times 10^{40}$ | $0.5$ | Intermediate halo mass threshold simulation
high_$M_{h}$ | $10^{10}$ | 0.45 | $3.16\times 10^{40}$ | $0.5$ | High halo mass threshold simulation
low_$L_{X}$ | $10^{8}$ | 0.14 | $3.16\times 10^{38}$ | $0.5$ | Low X-ray luminosity simulation
high_$L_{X}$ | $10^{8}$ | 0.14 | $3.16\times 10^{42}$ | $0.5$ | High X-ray luminosity simulation
low_$E_{0}$ | $10^{8}$ | 0.14 | $3.16\times 10^{40}$ | $0.2$ | Low X-ray energy threshold simulation
high_$E_{0}$ | $10^{8}$ | 0.14 | $3.16\times 10^{40}$ | $1$ | High X-ray energy threshold simulation

To aid our physical interpretation of the features present in the SS we run a
further set of six simulations with meraxes, in addition to our fiducial
simulation.

In particular, we are interested in the physical processes which impact the
morphology of the 21cm signal. To explore the impact of the ionization
morphology we vary the minimum mass for halos hosting star formation. Setting
the UV escape fraction to zero in galaxies below a given mass threshold alters
the size and distribution of the ionized regions (i.e. produces larger, more
isolated bubbles for an increasing mass threshold). With regard to the heating
morphology, we vary the X-ray luminosity and the minimum energy threshold for
X-rays escaping their host environment. Increasing the X-ray energy threshold
decreases the prevalence of bubbles of heated IGM gas transitioning toward an
effective uniform background of IGM heating. Table 1 summarises the
simulations used in this work along with the values of the parameters that are
varied.

####  3.4.1 Halo Mass Threshold

To explore the impact of the minimum halo mass on the EoR morphology, we
modify the $f_{\rm esc}$ (also see Equation 9 and Section 3.2) prescription as
follows:

$f_{\rm esc}=\begin{cases}f_{\rm
esc,0}\bigg{(}\frac{1+z}{6}\bigg{)}^{\alpha_{\rm esc}},&M_{\rm halo}\geq
M_{\rm thresh}\\\ \\\ 0,&M_{\rm halo}<M_{\rm thresh}.\\\ \end{cases}$ (16)

We run simulations with mass thresholds $M_{\rm thresh}=10^{9}$ and
$10^{10}\,M_{\odot}$, and label them mid_$M_{h}$ and high_$M_{h}$ respectively
(we point out that our Fiducial simulation contains all haloes down to
$3\times 10^{7}\,M_{\odot}$). By considering an increasing halo mass
threshold, we effectively decrease the total number of galaxies capable of
contributing to reionization. To compensate for the loss of ionizing sources
we increase the UV escape fraction of those remaining star-forming galaxies to
ensure a reionization history consistent with our observational constraints.
We therefore, increase the UV escape fraction normalisation $f_{\rm esc,0}$
(see Equation 9) to $[0.25,0.45]$ respectively (see Fig 2). These simulations
can thus probe the impact of ionization morphology in the PS and SS. An
increasing mass threshold should result in larger ionized regions (changing
the physical location of features in the PS/SS). We emphasise that we still
populate and evolve the galaxies in the haloes below the mass threshold, and
these galaxies can start contributing to the UV ionization budget when their
host halo mass passes the threshold. We point out that we do not suppress
emission of X-ray photons by these galaxies. In this manner, we fix the X-ray
background across these simulations (Fiducial, mid_$M_{h}$ and high_$M_{h}$ )
to be the same. This was a deliberate design choice to isolate the impact of
the EoR morphology on the 21-cm statistics.

####  3.4.2 X-Ray Luminosity

We also consider two simulations with a lower and higher X-ray luminosity
$L_{X<2\,\rm{keV}}/\rm{SFR}=[3.16\times 10^{38},3.16\times
10^{42}]\>\rm{ergs\,s^{-1}\,M_{\odot}^{-1}\,yr}$ as compared to the fiducial
value of $3.16\times 10^{40}\>\rm{ergs\,s^{-1}\,M_{\odot}^{-1}\,yr}$. These
two simulations, labelled low_$L_{X}$ and high_$L_{X}$ correspond to the
L210_AUG_LOWX and L210_AUG_HIGHX simulations in Balu et al. (2023). These
simulations cover a range of X-ray luminosities per SFR, one order of
magnitude broader than what is observed in the local soft band X-ray
luminosity (Mineo et al., 2012; Fialkov et al., 2016) based on the range
adopted in Greig & Mesinger (2017).

For low_$L_{X}$ , X-ray heating is inefficient and the IGM ionizes before it
is heated (21cm signal always remains in absorption). This produces large
temperature contrasts between the ionized and neutral regions resulting in
much higher amplitudes for the 21-cm statistics.

####  3.4.3 X-ray Energy Threshold

We explore the impact of the X-ray photon energy threshold $E_{0}$ by
producing two simulations with $E_{0}=0.2\,\rm{keV}$ and $E_{0}=1\,\rm{keV}$
compared to $E_{0}=0.5\,\rm{keV}$ for the Fiducial simulation. Decreasing the
energy threshold, coupled with our power-law X-ray spectral energy
distribution, results in a higher fraction of softer X-ray photons. As softer
photons have shorter mean free paths, more heat energy is deposited closer to
the host galaxies resulting in more prevalent bubbles of heating around the
first galaxies. Increasing the energy threshold removes this heating
morphology as the X-ray photons now penetrate much deeper into the IGM before
depositing their heat energy resulting in an effective uniform background of
heating. In effect, varying this energy threshold will alter the amplitude of
the 21-cm statistics during the heating epoch (see e.g. Pacucci et al., 2014;
Greig & Mesinger, 2017).

####  3.4.4 Fiducial no Spin Temperature

Figure 1: 21cm differential brightness temperature light cone slices for each
simulation volume as a function of redshift from $z=20-5$. The top slice is
the Fiducial simulation (labelled), and the bottom slice is the high_$E_{0}$
simulation. The colour bar is a symmetric log-scale, where blue indicates
absorption relative to the CMB, red indicates emission relative to the CMB,
and gray indicates either ionization or zero signal. Each lightcone slice is
fixed to the same temperature scale.

As a comparison we create an additional simulation of the Fiducial signal in
the spin temperature saturation limit where $T_{S}\gg T_{\mathrm{CMB}}$ (from
hereon Fid_no$T_{S}$ ). This approximation effectively sets the temperature
contrast independent of the spin temperature during reionization. In this
limit the temperature field is proportional to the matter density
$\delta_{\rho}$ and ionization field $x_{\mathrm{HI}}$ (Cooray, 2005;
Furlanetto et al., 2006; Lidz et al., 2007). Shimabukuro et al. (2017) and
Majumdar et al. (2018) both explore the matter density and ionization
bispectrum components to the total 21cm bispectrum in this limit. Majumdar et
al. (2018) in particular finds that the negative sign of the bispectrum might
be an important indicator of the ionization topology during reionization.

##  4 Results

In this section, we analyse the thermal cosmological history of the 21cm
brightness temperature signal for the simulation sets. We plot 2D slices of
the lightcone boxes as a function of redshift for each simulation. We also
calculate the neutral fraction $(\bar{x}_{\textrm{HI}})$ for each simulation
coeval box as a function of redshift, and compare the results of each
simulation. We then calculate the PS and SS for each simulation as a function
of redshift for large $(k\sim 0.1\,\mathrm{Mpc}^{-1})$ and small $(k\sim
1\,\mathrm{Mpc}^{-1})$ spatial scales.

###  4.1 21cm Lightcones

Fig 1 shows a lightcone slice of each simulation as a function redshift. In
descending order the panels show the Fiducial, mid_$M_{h}$ , high_$M_{h}$ ,
low_$L_{X}$ , high_$L_{X}$ , low_$E_{0}$ , and high_$E_{0}$ simulations. The
colour bar is a log symmetric colour map, where blue indicates absorption, and
red indicates emission relative to the CMB. Gray indicates zero temperature
difference. During reionization $\delta T_{b}=0$ is typically associated with
regions that are ionized.

The same $N$-body dark matter particle genesis simulations are used to
generate each of the different meraxes simulations. Therefore each simulation
has the same dark matter halo distribution. In Fig 1 this is evident at high
redshifts ($z\lesssim 20$) and during reionization ($z\lesssim 8$) in the
location and approximate size of the first ionization regions. There are some
obvious differences in the temperature contrast due to the different X-ray
heating parameters. Of note, we see that the low_$L_{X}$ is always in
absorption, even during reionization, and the high_$L_{X}$ simulation is
heavily preheated at high redshift. The low_$E_{0}$ simulation has small
regions of localised heated gas that appear in emission at $z>15$. In
contrast, the high_$E_{0}$ simulation results in a more uniform heating of the
IGM and thus the brightness temperature is relatively featureless. There is
also a clear difference in the size of ionization regions between the
Fiducial, mid_$M_{h}$ , and high_$M_{h}$ simulations, with the size increasing
from Fiducial to high_$M_{h}$ at fixed redshift. We discuss these features in
the context of the statistics in the following subsections.

###  4.2 Neutral Fraction and Ionization

Figure 2: The ionization history of the average neutral hydrogen (neutral
fraction $\bar{x}_{\textrm{HI}}$) IGM calculated for each simulation
(labelled).

In Fig 2 we show the average neutral fraction ($\bar{x}_{\textrm{HI}}$)
calculated for each simulation coeval box as a function of redshift. The
Fiducial model was calibrated to match existing observational constraints (see
Balu et al. (2023) for details). The halo mass threshold simulations have
increased UV escape fractions as a function of halo mass to ensure similar
reionization histories for easier comparison of the ionization morphology. The
solid black line is the average neutral fraction for the Fiducial simulation.
By $z=10$ the Fiducial simulation is already partially ionized at the $\sim 5$
percent level. We note that the ionization history for the low_$E_{0}$ (purple
dashed line), the high_$E_{0}$ (light purple dash dotted line), and the
low_$L_{X}$ (light green dash dotted line) simulations are effectively
identical to the Fiducial case. This is expected, since ionization is
predominantly driven by UV photons not X-ray emission, additionally these
simulations also have the same halo mass thresholds and escape fraction
$(f_{\rm esc,0})$ as the Fiducial simulation (see Table 1). In contrast the
high_$L_{X}$ (dark green dashed line) simulation undergoes reionization early
relative to the Fiducial, due to the increase in the number of ionizations
following secondary collisions of the X-ray photons. The high_$L_{X}$
simulation has a ionization fraction of $\sim 10$ percent reionization by
$z=10$. This is not unexpected since X-ray emission can be responsible for at
most $10$ percent of ionization (Mesinger et al., 2013).

The mid_$M_{h}$ (dark blue dashed line) and the high_$M_{h}$ (light blue dash
dotted line) simulations begin reionization later than the Fiducial
simulation, as it takes longer for haloes to gravitationally grow in excess of
their respective mass thresholds to emit ionizing UV photons. Nevertheless,
the mid_$M_{h}$ and the high_$M_{h}$ simulations neutral fraction profiles
result in a similar ionization history to the Fiducial. This is due to the
effective parameterisation of the escape fraction relative to the halo mass
threshold. To compensate for the loss of ionizing sources as a function of
increasing the halo mass threshold, the UV escape fraction was increased
proportional to the halo mass threshold. The mid_$M_{h}$ and high_$M_{h}$
simulations have UV $f_{\rm esc,0}$ values of $0.25$, and $0.45$, compared to
the Fiducial with $0.14$. Therefore, more ionizing UV photons escape per unit
mass from the same higher mass halos in high_$M_{h}$ and mid_$M_{h}$ , for the
same amount of star formation, compared to the Fiducial simulation. As
reionsation progresses, star formation increases, this results in a more rapid
(sharper) reionization relative to the Fiducial simulation.

###  4.3 21cm Statistics

####  4.3.1 Mean Brightness Temperature

Figure 3: Mean brightness temperature for all simulations, calculated from the
coeval boxes as a function of redshift.

Fig 3 shows the mean temperature for each simulation as a function of
redshift. As previously mentioned, the Fiducial model, low_$L_{X}$ and
high_$L_{X}$ simulations are taken from Balu et al. (2023). We see that the
mean temperature for these three simulations agree with those shown in Figure
8 of Balu et al. (2023), for more details on the low_$L_{X}$ and high_$L_{X}$
simulations we refer the reader to this work. For all simulations we see a
characteristic absorption feature which occurs when Ly-$\alpha$ emission
couples the spin temperature to the gas temperature. As the gas expands
adiabatically it cools relative to the CMB, increasing the relative
absorption. For most simulations this absorption trough occurs at
approximately $z\sim 15$ (with the exception of the high_$L_{X}$ simulation).
The timing of the trough depends on X-ray heating which eventually drives the
IGM into emission (with the exception of the low_$L_{X}$ simulation),
occurring during reionization at $z\gtrsim 10$. This culminates in a peak
roughly at the midpoint of reionization in the redshift range of $z\sim 6-8$.

Figure 4: The PS as a function of redshift for all simulations at a fixed
spatial scale. Panel (a) shows the redshift evolution at $k\sim
0.1\,\rm{Mpc}^{-1}$, and panel (b) shows the redshift evolution at $k\sim
1\,\rm{Mpc}^{-1}$.

As mentioned in the previous section the halo mass simulations both start
reionization later, but conclude earlier than the Fiducial. This delay results
in a higher heating and reionization rates for the halo mass simulations. The
higher ionization rate and heating being directly associated with the larger
escape fraction. This directly affects the mean temperature and the amplitude
of the PS (as seen in Fig 4).

The X-ray energy threshold simulations follow a similar evolution to the
Fiducial. The biggest difference is between the amplitude and onset of X-ray
heating between $10<z<15$. Here the low_$E_{0}$ (purple dashed line) undergoes
heating earlier than the Fiducial and the high_$E_{0}$ (light purple dash
dotted line) simulations. This should be expected since there are relatively
more lower energy soft X-rays available to heat the IGM. Additionally the
timing of heating should be earlier for the low_$E_{0}$ simulation because the
mean free path of X-ray photons is proportional to their energy, meaning
softer X-rays deposit their energy into the IGM before harder X-rays.

####  4.3.2 21cm Power Spectrum

Fig 4 shows the PS calculated at spatial scales $k\sim 0.1\,\mathrm{Mpc}^{-1}$
(panel (a)), and $k\sim 1\,\mathrm{Mpc}^{-1}$ (panel (b)), as a function of
redshift for each coeval box, and each simulation. The lines in Fig 4
correspond to the same simulations in Fig 3. The features discussed in Fig 3
are broadly mirrored in Fig 4 (a) for most simulations. We find the Fiducial
simulation in panel (a) of Fig 4 at $k\sim 0.1\,\rm{Mpc}^{-1}$ has a peak at
$z\sim 16$ during Ly-$\alpha$ pumping (Wouthuysen, 1952; Field, 1958), and a
peak at $z\sim 6.5$ during the mid point of reionization. Most simulations do
not have the characteristic three peak structure seen in (Furlanetto, 2006;
Mesinger et al., 2013), with the exception of the low_$E_{0}$ simulation and
the high_$L_{X}$ simulation, the latter which has four peaks. The peak during
the EoH for the Fiducial simulation has mostly merged with the peak during the
EoR due to the delayed heating.

We also include the PS as a function of redshift for the Fid_no$T_{S}$
simulation (dashed black line). We only show the Fid_no$T_{S}$ simulation up
to $z=12$, since the approximation is only valid during reionization. We find
good agreement between the Fiducial simulation and the Fid_no$T_{S}$
simulation during the middle and late periods of reionization at $k\sim
0.1\,\rm{Mpc}^{-1}$. At the larger scales the ionization topology is the
dominant component in the PS. At smaller scales X-ray heating and other
coupling effects are important for setting the spin temperature. These
differences result in disagreement at $k\sim 1\,\rm{Mpc}^{-1}$ since the spin
temperature calculation is omitted in the Fid_no$T_{S}$ simulation.

The mid_$M_{h}$ and the high_$M_{h}$ simulations in panel (a) have broadly the
same PS, differing at most in amplitude during the EoR, with the high_$M_{h}$
simulation peaking earlier than the mid_$M_{h}$ simulation. The mid_$M_{h}$
and high_$M_{h}$ simulations compared to the Fiducial have reionization
topologies driven by larger (and more biased) sources. This results in an
increase in the amplitude of the 21cm PS. The PS of both simulations are
broadly the same as the Fiducial simulation, only significantly deviating at
$z\sim 10$ when reionization begins to become more significant for the
Fiducial simulation.

For the low_$L_{X}$ and high_$L_{X}$ simulations, we find the same features at
$k\sim 0.1\,\rm{Mpc}^{-1}$ and $k\sim 1\,\rm{Mpc}^{-1}$ in Fig 4 as seen in
Balu et al. (2023) for the L210_AUG_lowX and L210_AUG_highX simulations444The
low_$L_{X}$ and high_$L_{X}$ simulations are the L210_AUG_lowX and
L210_AUG_highX from Balu et al. (2023).. The high amplitude of the low_$L_{X}$
simulation is due to the higher temperature contrast that results from the
colder IGM (lower X-ray luminosity). The lower amplitude of the high_$L_{X}$
simulation is due to the higher X-ray luminosity heating, this reduces the
temperature fluctuations on all scales. The inefficient heating due to lower
X-ray luminosity in the low_$L_{X}$ simulation results in the merging of the
EoH and EoR peaks. In the high_$L_{X}$ simulation we have four peaks, an early
peak at $z\sim 22$ coincident with the relatively weak absorption trough
(Ly-$\alpha$ pumping), a peak at $z\sim 18$ which corresponds to the EoH
heating, with the IGM being in emission at this stage. There is a peak during
the midpoint of reionization at $z\sim 7$. In addition to the expected three
peaks there is an additional peak at $z\sim 10$. This is due to the first
ionization sources, with the high_$L_{X}$ simulation having $\sim 10$ percent
ionization by $z=10$.

For the high_$E_{0}$ case at $k\sim 0.1\,\rm{Mpc}^{-1}$ the first peak occurs
at $z\sim 14$, which is roughly coincident with the absorption trough in the
mean temperature brightness in Fig 3 which reaches a minimum at $z\lesssim
14$. The second less prominent peak occurs during the EoR ($z\sim 6$). The
lower amplitude and later occurrence of the peak during the EoR are a direct
result of the higher X-ray energy threshold, which produce less structure in
the IGM due to the more uniform heating of the longer mean-free path X-ray
photons. The first patches of emission that are correlated on the largest
scales appear later relative to the other simulations, we see similar
behaviour in Figure 1 of Greig & Mesinger (2017) (bottom row of Figure 1).
This delayed and more uniform heating means the Ly-$\alpha$ pumping dominates
the amplitude of the PS for a longer period.

Figure 5: The SS as a function of redshift for all simulations at a fixed
spatial scale. Panel (a) shows the redshift evolution at $k\sim
0.1\,\textrm{Mpc}^{-1}$, panel (b) shows the redshift evolution at $k\sim
1\,\textrm{Mpc}^{-1}$.

For the low_$E_{0}$ simulation at fixed scale $k\sim 0.1\,\rm{Mpc}^{-1}$ we
see the three peaked structure. The peak due to Ly-$\alpha$ occurs early at
$z\sim 17$, with the peak during the EoH occurring at $z\sim 12$ which
coincides with the pockets of emission seen in Fig 1. The amplitude of the EoH
is greater for the low_$E_{0}$ simulation due to the lower X-ray energy
threshold. These lower energy X-rays efficiently heat the local medium around
the first luminous sources, producing higher temperature contrasts in the IGM
due to the more inhomogeneous heating. This produce more inhomogeneous
structures, increasing the overall power. These correlate on the largest
scales resulting in a strong peak. The high temperature contrast regions are
the first to ionize, thus during the EoR, the relative amplitude of the EoR
peak returns to that of the Fiducial model.

####  4.3.3 21cm Skew Spectrum

Fig 5 shows the SS calculated at $k\sim 0.1\,\rm{Mpc}^{-1}$ (panel (a)), and
$k\sim 1.0\,\rm{Mpc}^{-1}$ (panel (b)), as a function of redshift for each
coeval box, and each simulation. The lines correspond to the same simulations
in Fig 3 and 4.

The SS for each simulation, broadly mirrors the mean brightness temperature
for both small (panel (b)) and large (panel (a)) spatial scales in Fig 5. The
transition to a positive SS occurs rapidly within one coeval redshift box (a
cosmic blink). This happens in all simulations except low_$L_{X}$ which is
always in absorption555We do not include the Fid_no$T_{S}$ in this paragraph,
since it is only relevant during the EoR.. The transition to positive SS for
almost all simulations occurs earlier ($\Delta z<1$) and more rapidly than the
mean temperature brightness. Looking at Fig 1 there is a clear explanation.
During the transition period between the EoH and the EoR ($z\sim 10$), the
star formation rate increases. The X-ray luminosity is proportional to the
star formation rate. This results in the appearance of heated islands that are
in emission relative to the rest of the IGM which is undergoing more uniform
and less efficient heating. These heated islands skew the temperature
distribution towards positive values, resulting in a positive SS. This
eventually tapers as the IGM saturates and the heated regions overlap,
completing the transition from absorption to emission.

This transition happens significantly earlier for the low_$E_{0}$ simulation
($\Delta z\sim 2.5$). Due to the lower X-ray energy threshold, the low_$E_{0}$
simulation undergoes earlier and more intense local heating. The amplitude of
the effect appears to be scale dependent, with the transition in panel (b) at
smaller scales occurring much closer in redshift to the transition in the mean
temperature brightness. Conversely, for the high_$E_{0}$ simulation we see a
$\Delta z\sim 1$ lag in this transition at large and small scales. This is due
to the delayed and more uniform heating from the high X-ray energy threshold.

During reionization, the ionization morphology appears to dominate the signal
in the SS, with a peak associated with the midpoint of reionization, with the
exception of the low_$L_{X}$ simulation which has a trough (effectively
mirrored about the x-axis). As the neutral fraction drops below $50$ percent
most of the medium is ionized, and this tapers the SS, much like the PS,
driving the signal to zero as reionization progresses. Like the PS, the SS
amplitude is largest for the high_$M_{h}$ simulation at the largest scale due
to the fact that reionization is driven by larger, more biased galaxies.
However, the high_$L_{X}$ simulation has the highest overall amplitude at all
scales in Fig 5. The efficient X-ray heating in the high_$L_{X}$ simulation
drives the non-Gaussianities, and results in an earlier reionization, and
therefore a peak at ($z\sim 7$) compared to the other simulations.

Figure 6: The dimensionless SS plotted against the dimensionless PS as a
function of redshift for each simulation at a fixed spatial scale $k\sim
0.1\,\rm{Mpc}^{-1}$. Each curve indicates the co-evolution of the PS and SS
amplitude as a function of redshift, with the right pointed triangle,
indicating the start point at $z=20$, and the left pointing triangle
indicating the endpoint at $z=5$.

For contrast we calculate the SS for the Fid_no$T_{S}$ simulation for all
redshifts up to $z=12$ (black dashed line). At large scales ($k\sim
0.1\,\mathrm{Mpc}^{-1}$) we find a dip in the SS from positive to negative.
This feature agrees with similar work performed by Majumdar et al. (2018) who
look at different bispectrum of the expected 21cm signal in the saturation
limit ($T_{S}\gg T_{\mathrm{CMB}}$). This transition occurs because the first
ionization regions appear, which have zero signal and skew the distribution in
the negative direction (growing peak at zero in an otherwise positive
distribution). This result clearly contrasts with the other simulations, of
particular note is the difference in amplitude with Fid_no$T_{S}$ having $2-3$
orders of magnitude lower SS amplitude in comparison. This illustrates the
importance of the spin temperature in the higher order statistics.

MP23 investigate the SS during the EoH and the EoR with 21cmfast . In Figure 3
and 4 of MP23, they display the SS as a function of redshift at different
fixed scales (small and large) for several different simulations. In general
they find two positive peaks in the SS across all spatial scales for all
simulations. MP23 associate the first peak with X-ray heating which couples
the spin temperature to the matter density. However, the transition from
negative SS to positive occurs at $z=14-12$ for the simulations presented in
MP23. We find this transition is likely a result of the appearance of the
first ionization regions, since at this redshift range the simulations are at
the minimum temperature (maximal absorption relative to the CMB). The second
transition from positive to negative occurs at $z\sim 10$ when the neutral IGM
goes from absorption to emission, effectively flipping the temperature
distribution and its asymmetry about the y-axis. The final transition happens
when the IGM is more than $50$ percent ionized, in this case the remaining
neutral medium (which is in emission) skews the temperature distribution
towards positive values.

The simulations in MP23 have effectively the same X-ray parameters as the
Fiducial simulation in our work. However, we do not find a negative SS at the
start of reionization, instead we see a single peak and trough for the
Fiducial, mid_$M_{h}$ , high_$M_{h}$ and high_$E_{0}$ simulations.
Furthermore, we note that the absolute amplitude of our SS for all simulations
is $10-10^{2}$ times large than those in MP23. At the largest scale
$k=0.1\,\mathrm{Mpc}^{-1}$ for the low_$E_{0}$ , the high_$L_{X}$ and the
low_$L_{X}$ simulations we do find some weaker peak and trough features
towards the earlier stages of reionization $(z\sim 7.5)$, when the
Fid_no$T_{S}$ simulation SS is negative. The absolute change in amplitude as a
result of these features are on the order of $10^{3}-10^{5}\,\mathrm{mK}^{3}$.
This is comparable to the absolute change in MP23 that results in the a sign
transition. The lack of a sign transition for the meraxes simulations, in
contrast to MP23 and the Fid_no$T_{S}$ , is a clear indicator, that the spin
temperature is inherently far more skewed than in the latter cases. We further
discuss the differences, and the sources of this larger skewness in meraxes in
Section 7.1.

####  4.3.4 Skew Spectrum and Power Spectrum Co-evolution

Fig 6 shows the dimensionless SS and PS of each simulation as a function of
redshift. Each simulation set is separated into individual panels, with the
mid_$M_{h}$ , and high_$M_{h}$ in panel (a), the low_$L_{X}$ and high_$L_{X}$
in panel (b), and the low_$E_{0}$ , and high_$E_{0}$ simulations in panel (c).
Each curve is taken at $k\sim 0.1\,\rm{Mpc}^{-1}$, the colour coded triangles
show the direction the curve travels as a function of redshift. The start
point is indicated by the right facing triangle (relative to the peak), and
the endpoint being the left facing triangle.

The curves in Fig 6 demonstrate how the SS and PS co-evolve as a function of
redshift. For the Fiducial simulation we see the PS increases as the SS
becomes more negative, this occurs during absorption at early redshifts. Then
as the IGM heats due to X-ray emission, the PS amplitude decreases, and the SS
amplitude increases towards zero, until eventually becoming positive at a
fixed PS amplitude. The PS and SS then both increase as reionization and
heating occur in tandem, eventually culminating in a downward trend towards
zero after the midpoint of reionization. This evolution is broadly mirrored by
the other simulations with the exception of the low_$L_{X}$ and high_$L_{X}$
simulations.

For the halo mass simulations, the co-evolution of the SS and PS as a function
of redshift is practically identical compared to the Fiducial simulation. The
deviations occur due to the differences between the PS and SS amplitudes.
During the EoH since all simulations have the same X-ray background, the main
difference is due to heating from UV emission. For the Fiducial simulation
this can occur at lower halo masses. The heating leads to a reduction in the
PS amplitude for the Fiducial simulation when the IGM is in absorption; this
effect is however minimal compared to X-ray heating (Furlanetto et al., 2006).
Once reionization starts, the mid_$M_{h}$ and high_$M_{h}$ simulations can
emit more UV photons per unit mass per unit star formation, this results in
relatively more UV heating. When the IGM is in emission, this leads to higher
SS and PS amplitudes for the mid_$M_{h}$ and high_$M_{h}$ simulations compared
to the Fiducial.

For the low_$L_{X}$ simulation the SS is always negative, there is also a
fairly flat SS from $z=10$ to $z=6$, with a small peak. These corresponds to
the turning points in the curve, similar to the Fiducial simulation. The
high_$L_{X}$ simulation demonstrates deviation from the Fiducial evolution,
with the turning point from negative to positive SS amplitude happening at low
PS amplitude. We see more turning points with increasing PS amplitude for the
high_$L_{X}$ simulation which are correlated with the two peaks seen between
$10<z<20$ in Fig 4 panel (b). These finally culminate in a turning point at
high SS and PS amplitude during the peak of reionization, transitioning to the
zero amplitude for both.

Figure 7: PS (row one), SS (row two) and normalised SS (row three), of the
Fiducial (solid black line), mid_$M_{h}$ (double dotted dash line), and the
high_$M_{h}$ (dash dotted line). Each figure from left to right is the
$\bar{x}_{\textrm{HI}}\sim 0.75,\>0.5,$ and $0.25$ for each simulation.

For the low_$E_{0}$ and high_$E_{0}$ simulations, the evolution is notably
different. For the low_$E_{0}$ simulation, the structure mirrors the Fiducial
simulation, with the sharp transition in amplitude for the SS happening
earlier and at higher PS amplitudes due to the increased heating morphology.
The more uniform heating of the IGM in the case of the high_$E_{0}$ simulation
leads to less variation, and thus to less rapid growth of the SS and PS
amplitude with respect to the other simulations. Notably, the transition from
negative to positive SS amplitude occurs over a longer period for high_$E_{0}$
. The PS amplitude in this case varies as the SS amplitude transitions from
negative to positive; this is a key difference between the high_$E_{0}$ on
other simulations.

##  5 Power Spectrum and Skew Spectrum During reionization

In this section we calculate the spherically averaged PS and the spherically
averaged SS for each simulation during the epoch of reionization. We
investigate what information we are sensitive to during reionization, and what
measuring the PS and the SS together can reveal about the physics of
reionization. We investigate this by measuring the normalised SS, discussed in
the following section 5.1.

###  5.1 Normalised Skew Spectrum

The SS like the bispectrum is a measure of the non-Gaussianity through the
central third order moment statistics of the temperature field. Watkinson et
al. (2018) found significant fluctuations in the expected bispectrum around
the zero point. These large fluctuations are due to the PS amplitude present
in the statistic. Inspired by Eggemeier & Smith (2016), Watkinson et al.
(2018) normalises the bispectrum to a unitless ‘normalised’ bispectrum, which
is normalised by the PS and the $k$-modes. This normalisation is akin to
measuring skewness, which is the central third order moment normalised by the
cube of the standard deviation (the variance to the power of $3/2$). In this
work we perform a similar normalisation to remove the Gaussian component of
the amplitude in the SS. We normalise the SS by the PS taken to the power of
$3/2$, providing a unitless quantity $\gamma(k)$, referred to as the
normalised SS:

$\gamma(k)=\frac{\Delta^{2}_{T^{2}T}}{\left(\Delta^{2}_{T}\right)^{3/2}}.$
(17)

Figure 8: PS (row one), SS (row two) and normalised SS (row three), of the
Fiducial (solid black line), low_$L_{X}$ (double dotted dash line), and the
high_$L_{X}$ (dash dotted line). Each figure from left to right is the
$\bar{x}_{\textrm{HI}}\sim 0.75,\>0.5,$ and $0.25$ for each simulation.

In Equation 17 $\Delta^{2}_{T}$ and $\Delta^{2}_{T^{2}T}$ are the
dimensionless spherically averaged PS and SS, which have units of
$\textrm{mK}^{2}$ and $\textrm{mK}^{3}$ respectively. The subscripts $T$ and
$T^{2}T$ indicate the dimensionless PS and SS respectively. Deviations from a
flat distribution as a function of spatial scale will be indicative of when
the PS of the temperature fluctuations is more or less significant relative to
the SS.

Figure 9: PS (row one), SS (row two) and normalised SS (row three), of the
Fiducial (solid black line), low_$E_{0}$ (double dotted dash line), and the
high_$E_{0}$ (dash dotted line). Each figure from left to right is the
$\bar{x}_{\textrm{HI}}\sim 0.75,\>0.5,$ and $0.25$ for each simulation.

###  5.2 Power Spectrum, Skew Spectrum, and the Normalised Skew Spectrum

In this section we focus our investigation on what the PS and SS look like for
each simulation at different stages of reionization, specifically at neutral
fractions of $x_{\rm{HI}}\sim 0.75,\,0.5$, and $0.25$. This allows for a
direct comparison of the state of the IGM for each simulation set as compared
to the Fiducial simulation. In addition to the PS and SS at these states of
the IGM, we also calculate the normalised SS $(\gamma(k))$ from Equation 17
for each simulation. $\gamma(k)$ allows for the isolation of the non-
Gaussianity of the signal by normalising out the PS amplitude. In Figs 7, 8
and 9 we show the spherically averaged PS (first row), the spherically
averaged SS (second row), and $\gamma(k)$ (third row) as a function of spatial
scale, with the different neutral fractions of $\sim 0.75$, $\sim 0.5$ and
$\sim 0.25$ representing the first, second and third columns of each figure.

In Fig 7 we look at the PS, SS and the normalised SS of the halo mass
simulation set. In Fig 7 the Fiducial simulation is the solid black line, the
mid_$M_{h}$ the double dot dashed line, and the high_$M_{h}$ is the dashed dot
line. For both the PS and the SS we see a flattening of the spectra
($k<1\,\rm{Mpc}^{-1}$) as reionization progresses. This is typically related
to the ionization morphology (Zaldarriaga et al., 2004). Notably, the PS and
SS are very similar as a function of $k$, this is due to the Gaussian
component present in the SS amplitude.

In the normalised SS we see an interesting feature develop in all three
simulations, there is a local minima and local maxima in $\gamma(k)$. The
local minima is caused by the flattening of the PS due to the ionization
morphology. The flattening is more prevalent in the PS than the SS, and
therefore a minima appears in $\gamma(k)$ at the characteristic scales of
ionizing regions. The minima only becomes prevalent during the mid to late
stages of reionization once the ionized regions have percolated. The local
maxima on the other hand is due to small scale structures. At
$\bar{x}_{\rm{HI}}\sim 0.75$ the amplitude is tied to the non-linear
clustering of the ionized sources and their individual ionized bubbles. This
is evident due to the increased amplitude for the Fiducial simulation relative
to the mid_$M_{h}$ and high_$M_{h}$ simulations. The Fiducial simulation
contains many more smaller mass galaxies increasing the non-linear amplitude
on these scales. As reionization proceeds, the maxima grows in amplitude. At
$\bar{x}_{\rm{HI}}\sim 0.25$ the maxima is significantly larger for the
Fiducial simulation relative to the mid_$M_{h}$ and high_$M_{h}$ simulations.
The peak is now driven by the prevalence of neutral islands in the IGM (see
e.g. Fig 1). Due to the smaller escape fraction, the Fiducial simulation
contains the largest number of neutral islands. On the other hand, the
high_$M_{h}$ with its much larger escape fraction ionizes a larger volume per
ionizing source, preventing the appearance of neutral islands. Thus the
high_$M_{h}$ simulation does not exhibit a strong local maxima.

The location of the local minimum for the Fiducial simulation changes with
decreasing neutral fraction from $k\sim 0.2\,\rm{Mpc}^{-1}$ to $\sim
0.08\,\rm{Mpc}^{-1}$. The characteristic scales of ionizing regions is
expected to increase with redshift and decreasing (increasing) neutral
(ionization) fraction (Giri et al., 2017). The local minima at
$\bar{x}_{\rm{HI}}\sim 0.5$ appears at $k\sim 0.1\,\rm{Mpc}^{-1}$ for all
simulations. The expected size of the ionizing regions ranges from
$20-100\,\rm{Mpc}$ (Wyithe & Loeb, 2004; Zaldarriaga et al., 2004; Lin et al.,
2016), with the minima corresponding to size scales from $\sim
31-79\,\rm{Mpc}$. The local maxima on the other hand appears at $k\sim
1\,\rm{Mpc}^{-1}$ for the Fiducial, $k\sim 0.6\,\rm{Mpc}^{-1}$ for the
mid_$M_{h}$ and $k\sim 0.4\,\rm{Mpc}^{-1}$ for the high_$M_{h}$ simulation.
These correspond to sizes $6.3\,\rm{Mpc}$, $3.8\,\rm{Mpc}$, and
$1.9\,\rm{Mpc}$ for the Fiducial, mid_$M_{h}$ and the high_$M_{h}$ simulations
respectively.

Fig 8 has the same layout as Fig 7, here low_$L_{X}$ is the double dot dashed
line, and the dash dot line is the high_$L_{X}$ simulation. Here the
low_$L_{X}$ simulation power spectra differs from the Fiducial and
high_$L_{X}$ at all three ionization states. This difference is driven by the
larger absolute temperature difference of the low_$L_{X}$ simulation compared
to the other two simulations. The flattening of the low_$L_{X}$ PS at
$k<1\,\rm{Mpc}^{-1}$ is still prevalent, since reionization still proceeds the
same for the low_$L_{X}$ , high_$L_{X}$ and Fiducial simulation. The
differences at high $k$ are due to the signal temperature on small scales, in
the Fiducial and the high_$L_{X}$ simulations the signal is in emission, with
the high_$L_{X}$ simulation having a larger emission temperature which causes
larger temperature offsets with the ionized regions and thus produces a
higher-amplitude 21cm PS. This is not the case for the low_$L_{X}$ simulation
which remains in absorption as reionization progresses.666It should be noted
that at smaller scales $k>1\,\rm{Mpc}^{-1}$ approximations in the construction
of the simulation yield numerical artefacts that could potentially impact the
results.

When we look at the SS we see a similar trend, with the exception that
low_$L_{X}$ is negative compared to the Fiducial and high_$L_{X}$ scenarios.
The shape of low_$L_{X}$ relative to the other simulations is broadly mirrored
about the $k$-axis, and shows a similar morphology to the PS, as is similar
with the Fiducial and the high_$L_{X}$ simulations. The normalised SS shows
similar features to the Fiducial for both the low_$L_{X}$ (albeit negative)
and for the high_$L_{X}$ . There is a peak/trough at $k\sim
0.1\,\rm{Mpc}^{-1}$, and a subsequent peak/trough at $k\sim 1\,\rm{Mpc}^{-1}$.
However, this is less pronounced in the low_$L_{X}$ case, and this feature has
more of a flattening from $k=1-10\,\rm{Mpc}^{-1}$, which is also seen in the
PS. Due to the similarity in the ionization morphology between the Fiducial,
low_$L_{X}$ and the high_$L_{X}$ simulations especially at $x_{\rm{HI}}\sim
0.5$, and $0.25$, the amplitudes of the higher modes
($k>1\,\textrm{Mpc}^{-1})$ appear to be affected more by the X-ray heating,
than the ionization morphology itself.

Fig 9 shows the X-ray threshold simulation PS, SS and normalised SS. Here the
double dotted dashed line is the low_$E_{0}$ simulation, and high_$E_{0}$ is
the dashed dotted line. In all cases the SS and the PS of the Fiducial and the
low_$E_{0}$ simulations are almost identical during reionization, differing at
most at $k>1\,\rm{Mpc}^{-1}$. Again, the ionization morphology for all three
simulations is practically identical because they contain the same halo mass
distribution and threshold.

In the case of the high_$E_{0}$ simulation, increasing $E_{0}$ removes the
lowest energy photons, this moves the peak of the X-ray distribution to higher
energies. This results in a higher relative fraction of harder X-ray photons.
Since harder X-rays heat on larger scales (more uniform heating), this appears
to reduce the amplitude on all scales for both the PS and SS respectively.
However the reduction seems to impact the SS more than the PS. Fluctuations
and non-Gaussianities during reionization are largely driven by ionization
morphology at large scales ($k<1\,\textrm{Mpc}^{-1}$) and by the gas density
on small scales ($k\geq 1\,\textrm{Mpc}^{-1}$) (Watkinson et al., 2018;
Majumdar et al., 2018; Ma & Peng, 2023).

##  6 Detectability of the Normalised Skew Spectrum

There is a clear potential advantage to the normalised SS ($\gamma(k)$)
compared to the SS. The trough at $k\sim 0.1\,\rm{Mpc}^{-1}$ and the peak at
$k\sim 1\,\rm{Mpc}^{-1}$ demonstrate the sensitivity of the normalised SS to
non-Gaussianity in the 21cm signal. Additionally, calculating the SS
essentially requires the calculation of the PS, it is therefore
straightforward to construct the normalised SS by dividing out the Gaussian
amplitude. Naturally the detectability of the features present in $\gamma(k)$,
particularly the trough, is important to estimate if there is any potential
for its use as a probe of reionization. In this section we investigate the
statistical uncertainty of the PS and the SS, where we use these errors to
propagate the expected statistical uncertainty on the normalised SS. We save
the discussion of system noise from interferometers and foregrounds for future
work; these effects deserve independent consideration.

First, we consider the statistical uncertainty of the 21cm PS, also known as
the cosmic variance. Therefore of a random Gaussian field, the PS and by
extension the variance, describe all the information contained in the field.
In this case the uncertainty on the PS in the absence of thermal or
instrumental noise is determined by the Poisson sampling:

$\sigma(\Delta^{2}_{T})(k)=\Delta^{2}_{T}(k)\sqrt{\frac{(2\pi)^{2}}{Vk^{2}\Delta
k}}.$ (18)

The Poisson sampling error is proportional to one over the square root of the
number of Fourier modes for a given spherical shell with width $\Delta
k$777$\Delta k$ is not a constant since logarithmic bins are typically used to
calculate the PS. The logarithmic bin width $\log\Delta k=0.173.$. The number
of modes is proportional to the sampled co-moving volume of space $V$ used to
calculate the spectrum. Equation 18 assumes Gaussianity, however the signal
becomes non-Gaussian as reionization progress (Cooray, 2005; Furlanetto, 2006;
Wyithe & Morales, 2007). These non-Gaussianities correlate the signal on
different Fourier modes (Mondal et al., 2015a), which in turn introduces a
noise floor to the expected PS cosmic variance (Mondal et al., 2015b, 2016):

$\sigma(\Delta^{2}_{T})(k)=\Delta^{2}_{T}(k)\sqrt{\frac{(2\pi)^{2}}{Vk^{2}\Delta
k}}+\sqrt{\frac{T(k,k)}{V}}$ (19)

The non-Gaussianities in the PS are proportional to the Trispectrum $T(k,k)$
and $V^{-1/2}$. The effect of these non-Gaussianities is to flatten the signal
to noise providing a fundamental limit to the signal detectability (Mondal et
al., 2015a). A calculation of the analytic cosmic variance of the SS is
outside of the scope of this work. In this section we assume that the relative
uncertainties in the skew-spectrum follow a similar relationship to the PS.
Since the SS contains a Gaussian amplitude component, this seems reasonable,
and we will demonstrate this in the following sections.

###  6.1 Cosmic Variance of the Power Spectrum and the Skew Spectrum

(a) $\Delta^{2}_{T}(k)$ Uncertainty

(b) $\Delta^{2}_{T^{2}T}(k)$ Uncertainty

Figure 10: The numerically estimated statistical uncertainties on PS (a) and
the SS (b) (solid line), compared to the expected Gaussian uncertainties
(dashed line), and the estimated non-Gaussian component (dotted line).

To investigate the cosmic variance during the EoR, we follow the same method
outlined in Balu et al. (2023). We split each of the simulation volumes in
Figs 7, 8 and 9 into $27$ sub volumes each with side length
$70\,h^{-1}\,\rm{Mpc}$. We then calculate the PS and SS for each sub volume.
The cosmic variance for the PS and the SS is numerically estimated by
calculating the variance with respect to the mean power and SS for each
$k$-mode. Subfigures 10(a) and 10(b) show the sample error (solid black line)
for the PS and the SS, compared to the expected Gaussian uncertainty (dashed
line) for the Fiducial model at a neutral fraction of $x_{\textrm{HI}}\sim
0.5$. We also estimate the Non-Gaussian component (dotted line) by subtracting
the expected Gaussian errors from the sampled errors, these results are
similar to those shown in Figure 3 of Greig et al. (2022). We find qualitative
agreement with Mondal et al. (2015b, 2016); Greig et al. (2022), where we find
an earlier transition to non-Gaussianities in Subfigure 10(a) in agreement
with that seen by Balu et al. (2023).

The uncertainties at large scales are dominated by the Gaussian component at
$(k\lesssim 0.3\textrm{Mpc}^{-1})$ for the PS and $(k\lesssim
0.2\textrm{Mpc}^{-1})$ for the SS. This transition is seen at $(k\lesssim
0.5\textrm{Mpc}^{-1})$ in Mondal et al. (2015a) and Greig et al. (2022). Balu
et al. (2023) attributes the earlier transition to the spin temperature
evolution and detailed physical prescriptions for the higher non-Gaussianity
in the L210 box. Overall, we find the assumption that the uncertainties in the
SS have a similar form as the PS to be a good one.

###  6.2 Uncertainty in the Normalised Skew Spectrum

To estimate the uncertainty in the normalised SS we calculate the linear first
order error propagation of equation 17:

$R^{2}_{\gamma}(k)=\frac{9}{4}R^{2}_{T}(k)+R^{2}_{T^{2}T}(k)-3R_{T}(k)R_{T^{2}T}(k)\rho(k),$
(20)

$(R(k)=\sigma_{x}(k)/X(k))$ is the relative error for either the PS or the SS,
labelled with the subscripts $T$ and $T^{2}T$ respectively. $\rho(k)$ is the
Pearson correlation coefficient of the dimensionless PS and SS as a function
of spatial scale. We find significant correlation between the PS and the SS as
calculated from the sub volumes in the previous section.

Figure 11: Correlation coefficient between the Fiducial PS and SS for the
$x_{\textrm{HI}}\sim 0.75$ (solid), the $0.5$ (dash dotted line), and the
$0.25$ (double dot dashed line) neutral fractions.

Fig 11 shows the correlation for the $\bar{x}_{\textrm{HI}}\sim 0.75$, $0.5$,
and the $0.25$ Fiducial simulation coeval boxes as a function of $k$. The
average correlation is $0.82$, $0.88$ and $0.92$ respectively for each of the
coeval boxes. This result is not surprising and demonstrates that the SS
during reionization has significant contribution from the PS amplitude, which
is evident from the morphological similarity of the PS and the SS.

We can derive an expression for the Gaussian component errors in the
normalised skew spectrum if we assume that Trispectrum component of $R_{T}$
and $R_{T^{2}T}$ is zero $(T(k,k)=0)$. In this case the relative errors for
the PS and SS are equal $(R_{T}=R_{T^{2}T})$ and only depend on the volume $V$
and the shell volume $k^{2}\Delta k$. Equation 20 therefore simplifies to:

$R_{\gamma}(k)=\frac{R_{T}(k)}{2}\sqrt{13-12\rho(k)}.$ (21)

Figure 12: The first order propagated uncertainties on the normalised SS
(solid black line), compared to the Gaussian model uncertainties (dashed
line), and the estimated non-Gaussian component (dotted line), for the
Fiducial simulation at $\bar{x}_{\textrm{HI}}=0.5$.

We use Equation 21 as a model for the Gaussian component of the uncertainties
in $\gamma(k)$. Fig 12 shows the sample estimated errors (solid black line)
calculated from Equation 20. The red double dot dashed line shows the
estimated Gaussian uncertainties estimated from Equation 21. Finally the non-
Gaussian component (dotted line) was likewise estimated by subtracting the
Gaussian uncertainty model at $\bar{x}_{\textrm{HI}}=0.5$ from the full
uncertainty estimation. The transition from Gaussian dominated to non-Gaussian
dominated uncertainties occurs at $k\lesssim 0.3\textrm{Mpc}^{-1}$.
Interestingly, we find that the Non-Gaussian uncertainties are roughly flat as
a function of $k$.

####  6.2.1 Detection Predictions

In this section we perform a rudimentary signal to noise (S/N) estimate for
the normalised SS, for future SKA_LOW observations. For the estimate we assume
the Fiducial simulation as the 21cm signal, and we consider the neutral
fractions $\bar{x}_{\textrm{HI}}=0.25$, $\bar{x}_{\textrm{HI}}=0.5$, and
$\bar{x}_{\textrm{HI}}=0.75$, and their respective redshifts $7.4$, $6.6$ and
$6.1$.

The signal to noise ratio can be defined as the inverse of the relative
uncertainty ($\textrm{S/N}(k)=1/R_{\gamma}(k)$). In this case we assume the
full relative uncertainty for $\gamma(k)$, which includes the non-
Gaussianities in Equation 20. Notably, $R_{\gamma}\propto V^{-1/2}$
$(\textrm{S/N}\propto V^{1/2})$ therefore, to estimate the S/N for an SKA_LOW
observation, we can replace the simulation comoving volume
($V^{1/2}_{\rm{sub}}$) by the SKA_LOW comoving volume ($V^{1/2}_{\rm{SKA}}$).
This is done by first calculating the S/N for the Fiducial simulation, then
normalising out $V^{1/2}_{\rm{sub}}$, and finally scaling by
$V^{1/2}_{\rm{SKA}}$.

To determine the SKA_LOW comoving volume for an observation, we first need to
know the field of view $\Omega_{f}$ and the observing bandwidth $\Delta\nu$.
From these values we can determine the comoving volume for each redshift
(neutral fraction) (Hogg, 1999). The angular width of the main lobe of the
primary beam for an interferometer is approximately given by
$\theta\sim\lambda/D$, where $\lambda$ is the observing wavelength, and $D$ is
the station diameter which we assume is $35\,\rm{m}$ (Turner, 2015). Using the
observing wavelength for each neutral fraction we calculate the field of view
to be $2.87^{2}\,\deg^{2}$, $2.6^{2}\,\deg^{2}$, and $2.44^{2}\,\deg^{2}$, for
each $\bar{x}_{\textrm{HI}}$ respectively. Finally assuming an observing
bandwidth of $\sim 30\,\rm{MHz}$ for each neutral fraction, we then determine
the comoving volume using Equation 28 from Hogg (1999). We then calculate the
S/N for each neutral fraction for the Fiducial simulation by taking the ratio
of $\gamma(k)/\sigma(\gamma(k))$. The S/N is then scaled by
$\sqrt{V_{\rm{SKA}}/V_{\rm{sub}}}$.

Figure 13: Signal to noise ratio of $\gamma(k)$ for the Fiducial simulation
scaled by the expected SKA_LOW observing comoving volume. The solid curve
corresponds to the Fiducial neutral fractions $x_{\textrm{HI}}\sim 0.75$,
$0.5$ for the dash dotted line, and $0.25$ for the double dot dashed line.

Fig 13 shows the resulting estimated S/N for $\gamma(k)$ during reionization
for $\bar{x}_{\textrm{HI}}=0.25$ (dashed dot line),
$\bar{x}_{\textrm{HI}}=0.5$ (dash double dot line), and
$\bar{x}_{\textrm{HI}}=0.75$ (solid line). The tapering of the S/N as a
function of $k$ is characteristic of the non-Gaussian component in the cosmic
variance and is also seen in Figure 4 from Mondal et al. (2015a). We see that
for all neutral fractions for all spatial scales the $\textrm{S/N}>10$, with a
max signal to noise of $\sim 300$ for $\bar{x}_{\textrm{HI}}=0.75$. Thus, in
the absence of thermal noise, foregrounds and systematics, the features
present in $\gamma(k)$ should be detectable. MP23 perform their own S/N
analysis for the SKA_LOW with the addition of $1000$h of thermal noise. They
find a S/N of $\sim 20$ for their SS estimates. In the ideal case assuming
foregrounds and systematic noise can be removed, we should expect to be
sensitive to the normalised SS trough at $k=0.1\,\textrm{Mpc}^{-1}$ and peak
at $k\leq 1\,\textrm{Mpc}^{-1}$.

##  7 Discussion and Conclusion

We investigate the PS and SS as a function of redshift and spatial scale, for
a set of seven meraxes simulations. We vary the halo mass threshold, the X-ray
luminosity per star formation rate and the X-ray energy threshold. We find the
SS as a function of redshift broadly follows the mean differential temperature
brightness as a function of redshift, at large ($k\sim
0.1\,\mathrm{Mpc}^{-1}$) and small ($k\sim 1\,\mathrm{Mpc}^{-1}$) scales. We
do not see a negative sign for the SS amplitude during reionization as seen by
MP23; we discuss this in the following subsection 7.1 We further investigate
the ionization state and statistics of the IGM during the EoR. We look at the
spherically averaged PS, SS and the normalised SS for each simulation set
compared to the Fiducial at $\bar{x}_{\textrm{HI}}\sim 0.75,\>0.5$ and $0.25$.
In all simulations we find a local minimum at $k\sim 0.1\,\textrm{Mpc}^{-1}$
during the midpoint of reionization. This minimum corresponds to the
characteristic ionization topology (bubble) size during the EoR (Wyithe &
Loeb, 2004; Furlanetto & Oh, 2005; Lin et al., 2016), and we see the evolution
of the minima to larger scales with decreasing neutral fraction (increasing
ionization fraction). We expect that this feature should be detectable at
$k=0.1\,\textrm{Mpc}^{-1}$ in the absence of instrumental, thermal noise and
foreground contamination for the SKA_LOW and by extension current
interferometric experiments. We also see evidence of a local maxima in the
halo mass simulations, that shifts to larger scales as a function of the halo
mass threshold. This corresponds to small ionized or hot regions around these
halos. Our study highlights the importance of higher order statistics, and
what additional astrophysical information might be gained from calculating
both the SS and the PS.

The halo mass threshold simulations display the importance of ionization
topology as demonstrated in Fig 7. This clearly has the biggest impact on the
structure present in the normalised SS during reionization. There are however
some important limitations and caveats related to the halo mass simulations.
These simulations varied the ionization morphology by restricting the halo
mass threshold above which galaxies could emit ionizing UV photons, and
scaling the amount of UV emission as a function of the halo mass threshold.
Galaxies below the threshold however still produced X-ray emission, and thus
contribute to heating the IGM. Although this model is nonphysical, it allowed
for the halo mass threshold simulations to have a comparable X-ray background
to the Fiducial. This effectively isolated the impact of the ionization
topology on the 21cm statistics independent of X-ray heating. Changing the
X-ray emission in the same manner would delay heating, and have an undesirable
impact on the 21cm statistics.

###  7.1 Sign of the SS During Reionization

(a) meraxes slice $(x_{\mathrm{HI}}\sim 0.98)$

(b) meraxes slice $(x_{\mathrm{HI}}\sim 0.75)$

(c) meraxes PDFs

(d) meraxes SS

Figure 14: Subfigure (a) shows a 2D slice through the meraxes $z\sim 12$
coeval box, with $x_{\mathrm{HI}}\sim 0.98$. Subfigure (b) shows the same
slice, with a the $x_{\mathrm{HI}}\sim 0.75$ ionization field taken from
$z\sim 7.24$. Subfigure (c) shows the PDFs for both coeval boxes with the
different ionization fields, and subfigure (d) shows the resulting skew
spectra for both examples.

meraxes like 21cmfast calculates the temperature brightness field using the
excursion set formalism. However, even though our Fiducial simulation has the
same X-ray parameters as the simulations in MP23, we see some significant
differences between our results. Most notably the lack of a sign change
(positive to negative) during reionization in our SS for all simulations
(which include the spin temperature calculation). Additionally, we find the
absolute amplitude of our SS is approximately two orders of magnitude greater
than MP23, indicating that our differential temperature brightness
distribution is significantly more skewed on all scales. We do however see the
negative SS amplitude during reionization for the Fid_no$T_{S}$ simulation,
similar to Majumdar et al. (2018). We note that Majumdar et al. (2018) perform
a similar simulation with 21cmfast in the spin saturated limit. However, there
are limitations to this assumption and the spin temperature evolution as seen
in this work is important for understanding the higher order statistics of the
21cm signal.

The origin of the positive to negative sign transition (or vice versa), is a
direct result of the ionization morphology. For a mean zero differential
temperature brightness distribution, the ionized cells take the negative value
of the mean distribution (prior to mean subtraction). This acts to skew the
distribution in the opposite direction at the relevant scales for
reionization. The lack of a sign change in meraxes is due to the significantly
larger asymmetry (skewness) in the temperature distribution at all stages of
reionization and prior. This larger asymmetry is a direct result of the of the
spin temperature evolution in meraxes compared to 21cmfast . A complete
comparison of these two simulation packages is far beyond the scope of this
paper, here we summarise the key differences which result in the lack of a
sign change in our SS during reionization.

Two major difference between meraxes and 21cmfast are the halo mass function
(HMF) and the determination of the density fields. The HMF and the density
field of meraxes is determined from the input merger trees of N-body
simulations which are inherently non-linear. Additionally, the individual
haloes in meraxes are discrete and treated as independent galaxies which have
their own properties. 21cmfast on the other hand uses the Press-Schechter
conditional HMF renormalised to Sheth et al. (2001) , with a continuous
density field calculated by smoothing over the perturbed density field.
Furthermore, the Sheth et al. (2001) HMF is known to overproduce haloes by a
factor of two compared to N-body simulations (Watson et al., 2013; Murray et
al., 2013), thus there are more sources per unit mass when compared to
meraxes. Another significant difference is how star formation is handled in
meraxes compared to 21cmfast . In meraxes star formation for a galaxy can
cease due to feedback effects such as supernovae and AGN, this can lead to
intermittent star formation, and thus intermittent heating and UV production.
In contrast, in 21cmfast , once there is sufficient mass for star formation in
a given cell, star formation will continue unabated. Integrating this effect
over the time it takes to build the X-ray background, implies that 21cmfast
produces more sources, and more UV ionization per source compared to meraxes.
This difference is particularly notable when comparing the depth of the
absorption trough for the Fiducial simulation in Fig 3 to the absorption
trough in Figure 1 in MP23. The absorption trough for the Fiducial simulation
is almost a factor of two deeper than the MP23 simulations, moreover the
duration of the heating period is also more than a factor of two greater.

These effects, along with the inherently more skewed non-linear density in
meraxes from the N-body simulations, results in a more negatively skewed
temperature brightness distribution at all scales when compared to 21cmfast .
This asymmetry grows as the IGM adiabatically cools. When the IGM transitions
from negative to positive the highly skewed tail also transitions to positive.
The inherent asymmetry in the differential temperature distribution outweighs
that introduced by the ionization morphology on relevant scales.

The brightness distribution of the simulations in MP23 relative to meraxes is
not as strongly negatively skewed, and thus a small amount of ionization is
required to cause the initial sign transition. Therefore, meraxes with its
much larger negative asymmetry requires a significantly larger ionization
fraction to cause the equivalent sign flip. To better understand what level of
ionization is required to induce a sign transition with meraxes, we apply
varying levels of ionization fields from lower redshifts to the Fiducial
meraxes simulation at $z\sim 12$ ($x_{\mathrm{HI}}\sim 0.98$). Specifically,
we apply varying levels of ionization morphology from
$\bar{x}_{\mathrm{HI}}\sim 0.98$ to $\bar{x}_{\mathrm{HI}}\sim 0.75$ to the
$z\sim 12$ Fiducial simulation volume, to create pseudo-ionization states at
$z\sim 12$.

A 2D slice of the Fiducial simulation at $z\sim 12$ is shown in Subfig 14(a),
and in Subfig 14(b) we show the same slice with the $\bar{x}_{\mathrm{HI}}\sim
0.75$ ($z\sim 7.24$) ionization field applied. We calculate the differential
temperature brightness distribution for each pseudo-ionization state, these
are displayed in Subfig 14(c), with the original box in solid blue. The
decreasing ionization fraction results in shifting the peak of the
distribution closer to zero. We then calculate the SS for each pseudo-
ionization state, these are shown in Subfig 14(d). At large scales ($k\sim
0.1\,\mathrm{Mpc}^{-1}$) for ionization states of $\bar{x}_{\mathrm{HI}}<0.9$
(greater than $10\%$ ionization) the SS is positive, compared to the initial
negative SS of the Fiducial simulation at $z\sim 12$. This is similar to the
results seen in Figure 3 of MP23.

###  7.2 Future Work

Future work should consider how the ionization topology is quantitatively
linked to the trough feature seen in the halo mass simulations. This could be
investigated by measuring the characteristic scales of ionization for the
coeval boxes in this work using the methods outlined in Lin et al. (2016), in
particular the mean free path method. Understanding the characteristic
ionization scale and variance, and comparing them to the trough position and
width are important for quantitatively understanding how the ionization
topology affects the higher order statistics. This could additionally be
performed on a simple toy model similar to the one used in Majumdar et al.
(2018), where the characteristic scales can be varied to understand their
impact on the observed features in the normalised SS, independent of other
affects such as X-ray heating.

This work highlights the importance of the spin temperature in the higher
order statistics of the expected 21cm signal. Furthermore, this work
demonstrates the impact of the different implementations of star formation in
the 21cm signal between meraxes and 21cmfast . More work is required to
understand the significance of these differences and how they impact the
expected higher order statistics of the 21cm signal.

Future work will consider the practicality of calculating the SS with the
current and future generations of radio interferometric instruments. The
squared sky temperature can not be measured directly, so we must estimate it
from interferometric visibilities derived from observations888In Section A of
the appendix, we outline one method and its challenges for estimating the SS
from radio interferometric visibilities.. This ultimately involves a
convolution of the signal in Fourier space or a multiplication in image space
and subsequent inversion back to Fourier space. Each additional step in the
process propagates systematic and instrumental effects and spreads them across
different Fourier modes. These effects already impact the PS for numerous
experiments which are systematics limited rather than thermal noise limited.
Further investigation is needed to understand how these effects propagate
through the SS, and whether they render a realistic measurement impractical.

## Acknowledgements

We acknowledge the helpful discussion of Simon Mutch who assisted with this
work. This research was supported by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project
number CE170100013. JHC is supported by a Research Training Program
scholarship. CMT is supported by an ARC Future Fellowship under grant
FT180100321. The International Centre for Radio Astronomy Research (ICRAR) is
a Joint Venture of Curtin University and The University of Western Australia,
funded by the Western Australian State government. The MWA Phase II upgrade
project was supported by the Australian Research Council LIEF grant
LE160100031 and the Dunlap Institute for Astronomy and Astrophysics at the
University of Toronto. This scientific work makes use of the Murchison Radio-
astronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji
people as the traditional owners of the Observatory site. Support for the
operation of the MWA is provided by the Australian Government (NCRIS), under a
contract to Curtin University administered by Astronomy Australia Limited. We
acknowledge the Pawsey Supercomputing Centre which is supported by the Western
Australian and Australian Governments.

Part of this work was performed on the OzSTAR national facility at the
Swinburne University of Technology. The OzSTAR program partially receives
funding from the Astronomy National Collaborative Research Infrastructure
Strategy (NCRIS) allocation provided by the Australian Government. This
research was also undertaken with the assistance of resources from the
National Computational Infrastructure (NCI Australia), an NCRIS-enabled
capability supported by the Australian Government.

## Data Availability

The simulated data referenced in this work, is available upon reasonable
request to the primary author.

## References

  * Balu et al. (2023) Balu S., Greig B., Qiu Y., Power C., Qin Y., Mutch S., Wyithe J. S. B., 2023, Monthly Notices of the Royal Astronomical Society, 520, 3368
  * Barkana & Loeb (2001) Barkana R., Loeb A., 2001, Physics Reports, 349, 125
  * Barkana & Loeb (2007) Barkana R., Loeb A., 2007, Reports on Progress in Physics, 70, 627
  * Bharadwaj & Pandey (2005) Bharadwaj S., Pandey S. K., 2005, MNRAS, 358, 968
  * Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  * Bosman et al. (2022) Bosman S. E. I., Davies F. B., Becker G. D., Keating L. C., Davies R. L., Zhu Y., Eilers A.-C., 2022, Monthly Notices of the Royal Astronomical Society, 514, 55
  * Bower (1991) Bower R. G., 1991, Monthly Notices of the Royal Astronomical Society, 248, 332
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Bridle & Schwab (1999) Bridle A. H., Schwab F. R., 1999, in Taylor G. B., Carilli C. L., Perley R. A., eds, Astronomical Society of the Pacific Conference Series Vol. 180, Synthesis Imaging in Radio Astronomy II. p. 371
  * Chiang et al. (2014) Chiang C.-T., Wagner C., Schmidt F., Komatsu E., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 048
  * Cooray (2001) Cooray A., 2001, Phys. Rev. D, 64, 043516
  * Cooray (2005) Cooray A., 2005, Monthly Notices of the Royal Astronomical Society, 363, 1049
  * Dai et al. (2020) Dai J.-P., Verde L., Xia J.-Q., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 007
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Eggemeier & Smith (2016) Eggemeier A., Smith R. E., 2016, Monthly Notices of the Royal Astronomical Society, 466, 2496
  * Fialkov et al. (2016) Fialkov A., Cohen A., Barkana R., Silk J., 2016, Monthly Notices of the Royal Astronomical Society, 464, 3498
  * Field (1958) Field G. B., 1958, Proceedings of the IRE, 46, 240
  * Field (1959) Field G. B., 1959, ApJ, 129, 536
  * Fragos et al. (2013) Fragos T., et al., 2013, The Astrophysical Journal, 764, 41
  * Furlanetto (2006) Furlanetto S. R., 2006, Monthly Notices of the Royal Astronomical Society, 371, 867
  * Furlanetto (2016) Furlanetto S. R., 2016, The 21-cm Line as a Probe of Reionization. Springer International Publishing, Cham, pp 247–280, doi:10.1007/978-3-319-21957-8_9, %****␣main.bbl␣Line␣125␣****https://doi.org/10.1007/978-3-319-21957-8_9
  * Furlanetto & Furlanetto (2006) Furlanetto S. R., Furlanetto M. R., 2006, Monthly Notices of the Royal Astronomical Society, 374, 547
  * Furlanetto & Oh (2005) Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031
  * Furlanetto et al. (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, The Astrophysical Journal, 613, 1
  * Furlanetto et al. (2006) Furlanetto S. R., Peng Oh S., Briggs F. H., 2006, Physics Reports, 433, 181
  * Giri et al. (2017) Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2017, Monthly Notices of the Royal Astronomical Society, 473, 2949
  * Greig & Mesinger (2017) Greig B., Mesinger A., 2017, Monthly Notices of the Royal Astronomical Society, 472, 2651
  * Greig et al. (2022) Greig B., Wyithe J. S. B., Murray S. G., Mutch S. J., Trott C. M., 2022, Monthly Notices of the Royal Astronomical Society, 516, 5588
  * Hogg (1999) Hogg D. W., 1999, arXiv e-prints, pp astro–ph/9905116
  * Hutter et al. (2019) Hutter A., Watkinson C. A., Seiler J., Dayal P., Sinha M., Croton D. J., 2019, Monthly Notices of the Royal Astronomical Society, 492, 653
  * Kennicutt (1998) Kennicutt R. C., 1998, The Astrophysical Journal, 498, 541
  * Lacey & Cole (1993) Lacey C., Cole S., 1993, Monthly Notices of the Royal Astronomical Society, 262, 627
  * Lidz et al. (2007) Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007, The Astrophysical Journal, 659, 865
  * Lin et al. (2016) Lin Y., Oh S. P., Furlanetto S. R., Sutter P. M., 2016, Monthly Notices of the Royal Astronomical Society, 461, 3361
  * Ma & Peng (2023) Ma Q.-B., Peng L., 2023, Monthly Notices of the Royal Astronomical Society, 523, 640
  * Majumdar et al. (2018) Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A., Bharadwaj S., Mellema G., 2018, Monthly Notices of the Royal Astronomical Society, 476, 4007
  * Majumdar et al. (2020) Majumdar S., Kamran M., Pritchard J. R., Mondal R., Mazumdar A., Bharadwaj S., Mellema G., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5090
  * Mesinger et al. (2011) Mesinger A., Furlanetto S., Cen R., 2011, Monthly Notices of the Royal Astronomical Society, 411, 955
  * Mesinger et al. (2013) Mesinger A., Ferrara A., Spiegel D. S., 2013, Monthly Notices of the Royal Astronomical Society, 431, 621
  * Mineo et al. (2012) Mineo S., Gilfanov M., Sunyaev R., 2012, MNRAS, 419, 2095
  * Mondal et al. (2015a) Mondal R., Bharadwaj S., Majumdar S., Bera A., Acharyya A., 2015a, Monthly Notices of the Royal Astronomical Society: Letters, 449, L41
  * Mondal et al. (2015b) Mondal R., Bharadwaj S., Majumdar S., 2015b, Monthly Notices of the Royal Astronomical Society, 456, 1936
  * Mondal et al. (2016) Mondal R., Bharadwaj S., Majumdar S., 2016, Monthly Notices of the Royal Astronomical Society, 464, 2992
  * Mondal et al. (2021) Mondal R., Mellema G., Shaw A. K., Kamran M., Majumdar S., 2021, Monthly Notices of the Royal Astronomical Society, 508, 3848
  * Morales & Wyithe (2010) Morales M. F., Wyithe J. S. B., 2010, Annual Review of Astronomy and Astrophysics, 48, 127
  * Munshi et al. (1998) Munshi D., Melott A. L., Coles P., 1998, Generalised Cumulant Correlators and Hierarchical Clustering, doi:10.48550/ARXIV.ASTRO-PH/9812271, https://arxiv.org/abs/astro-ph/9812271
  * Murray et al. (2013) Murray S., Power C., Robotham A., 2013, Astronomy and Computing, 3-4, 23
  * Mutch et al. (2016) Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, Monthly Notices of the Royal Astronomical Society, 462, 250
  * Nambissan et al. (2021) Nambissan J. T., et al., 2021, Experimental Astronomy, 51
  * Pacucci et al. (2014) Pacucci F., Mesinger A., Mineo S., Ferrara A., 2014, Monthly Notices of the Royal Astronomical Society, 443, 678
  * Peebles (1981) Peebles P. J. E., 1981, The Large-Scale Structure of the Universe. Princeton University Press, Princeton, doi:doi:10.1515/9780691206714, https://doi.org/10.1515/9780691206714
  * Planck Collaboration et al. (2021) Planck Collaboration et al., 2021, A&A, 652, C4
  * Pritchard & Furlanetto (2007) Pritchard J. R., Furlanetto S. R., 2007, Monthly Notices of the Royal Astronomical Society, 376, 1680
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Qin et al. (2017a) Qin Y., et al., 2017a, MNRAS, 472, 2009
  * Qin et al. (2017b) Qin Y., et al., 2017b, Monthly Notices of the Royal Astronomical Society, 472, 2009
  * Qiu et al. (2019) Qiu Y., Mutch S. J., da Cunha E., Poole G. B., Wyithe J. S. B., 2019, MNRAS, 489, 1357
  * Qiu et al. (2021) Qiu Y., Mutch S. J., Elahi P. J., Poulton R. J. J., Power C., Wyithe J. S. B., 2021, Monthly Notices of the Royal Astronomical Society, 500, 493
  * Regan (2017) Regan D., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 020
  * Ross et al. (2019) Ross H. E., Dixon K. L., Ghara R., Iliev I. T., Mellema G., 2019, Monthly Notices of the Royal Astronomical Society, 487, 1101
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, Monthly Notices of the Royal Astronomical Society, 323, 1
  * Shimabukuro et al. (2016) Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2016, Monthly Notices of the Royal Astronomical Society, 458, 3003
  * Shimabukuro et al. (2017) Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2017, Monthly Notices of the Royal Astronomical Society, 468, 1542
  * Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, Monthly Notices of the Royal Astronomical Society, 440, 1662
  * Songaila & Cowie (2010) Songaila A., Cowie L. L., 2010, The Astrophysical Journal, 721, 1448
  * Szapudi & Szalay (1997) Szapudi I., Szalay A. S., 1997, A New Class of Estimators for the N-point Correlations, doi:10.48550/ARXIV.ASTRO-PH/9704241, https://arxiv.org/abs/astro-ph/9704241
  * Tingay et al. (2013) Tingay S. J., et al., 2013, Publ. Astron. Soc. Australia, 30, e007
  * Tiwari et al. (2022) Tiwari H., Shaw A. K., Majumdar S., Kamran M., Choudhury M., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 045
  * Trott et al. (2019) Trott C. M., et al., 2019, Publications of the Astronomical Society of Australia, 36, e023
  * Turner (2015) Turner W., 2015, SKA-TEL-SKO-0000008, SKA Phase 1 System Level 1 Requirements, p. 40
  * Watkinson & Pritchard (2015) Watkinson C. A., Pritchard J. R., 2015, Monthly Notices of the Royal Astronomical Society, 454, 1416
  * Watkinson et al. (2017) Watkinson C. A., Majumdar S., Pritchard J. R., Mondal R., 2017, Monthly Notices of the Royal Astronomical Society, 472, 2436
  * Watkinson et al. (2018) Watkinson C. A., Giri S. K., Ross H. E., Dixon K. L., Iliev I. T., Mellema G., Pritchard J. R., 2018, Monthly Notices of the Royal Astronomical Society, 482, 2653
  * Watkinson et al. (2020) Watkinson C. A., Trott C. M., Hothi I., 2020, Monthly Notices of the Royal Astronomical Society, 501, 367
  * Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, Monthly Notices of the Royal Astronomical Society, 433, 1230
  * Wayth et al. (2018) Wayth R. B., et al., 2018, Publ. Astron. Soc. Australia, 35, e033
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Wyithe & Loeb (2004) Wyithe J. S. B., Loeb A., 2004, Nature, 432, 194
  * Wyithe & Morales (2007) Wyithe J. S. B., Morales M. F., 2007, Monthly Notices of the Royal Astronomical Society, 379, 1647
  * Yoshiura et al. (2015) Yoshiura S., Shimabukuro H., Takahashi K., Momose R., Nakanishi H., Imai H., 2015, Monthly Notices of the Royal Astronomical Society, 451, 266
  * Zahn et al. (2007) Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., Zaldarriaga M., Furlanetto S. R., 2007, The Astrophysical Journal, 654, 12
  * Zaldarriaga et al. (2004) Zaldarriaga M., Furlanetto S. R., Hernquist L., 2004, The Astrophysical Journal, 608, 622
  * Zarka et al. (2012) Zarka P., Girard J. N., Tagger M., Denis L., 2012, in Boissier S., de Laverny P., Nardetto N., Samadi R., Valls-Gabaud D., Wozniak H., eds, SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics.
  * Zygelman (2005) Zygelman B., 2005, ApJ, 622, 1356
  * van Haarlem, M. P. et al. (2013) van Haarlem, M. P. et al., 2013, A&A, 556, A2

##  Appendix A Estimating the SS From Visibilities

The most straightforward method for estimating the SS, is to calculate it from
image lightcones constructed form radio interferometric visibilities. In this
case, visibilities measured at different frequencies are used to constructed
image slices at each frequency. These images can then be mean subtracted,
squared and Fourier transformed. Multiplying the quadratic Fourier temperature
field by the conjugate of the Fourier temperature field, and averaging
spherically will estimate the SS. There are however numerous challenges. First
assume an ideal scenario where the foregrounds are perfectly subtracted, and
are not thermal noise limited. The image noise may not be entirely Gaussian
due to the spatial correlations introduced by the image PSF. Additionally,
squaring the image changes the noise statistics. Furthermore, primary beam
correcting the image would upscale the image noise. A solution might be to
consider smaller image subset where the primary beam is relatively constant.
This will limit the Fourier space resolution, but may be a desirable trade off
when considering instrumental effects introduced by the primary beam and the
image noise. This could be achieved with a spatial taper such as a circular
top hat to limit aliasing in Fourier space. The other limitation to this
method is the computational cost of imaging each channel, for low frequency
radio interferometers there are typically of order $100$ channels for a full
bandwidth observation. One solution could be to image only a subset of
channels, or to average channels together at the expense of bandwidth smearing
effects (Bridle & Schwab, 1999). Estimating the SS from radio interferometric
measurements, and the associated challenges will be focus of future work.