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

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.


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). 1 .These ionized bubbles grow and eventually overlap, culminating in the end of reionization by redshift ∼ 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 ★ E-mail: Jaiden.cook@curtin.edu.au 1 Reionization predominately occurs due to UV photons, with some contribution from X-ray emission (up to ∼ 10 percent Mesinger et al. (2013)).
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 variance 2 .The signal is however expected to 2 Since we cannot truly measure the ensemble average power spectrum, we 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 (  ≫  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 ∼ 10 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 ∼ 20 − 100 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. 2016Shimabukuro et al. , 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- 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 can only estimate it over some volume.Each independent realisation therefore is a random sample of the true PS with some cosmic variance.
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 ()) (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 .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 ℎ −1 Mpc simulations were augmented to include all atomically cooled galaxies out to  = 20 (∼ 2×10 7 ℎ −1 M ⊙ ).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): ℎ = 0.68, Ω  = 0.31, Ω  = 0.048, Ω Λ = 0.69,  8 = 0.81, and   = 0.96.All cosmological scales are in comoving units.

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

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  (Peebles 1981): (1) In Equation 1the angular brackets ⟨⟩ denote the ensemble average over different realisations of the Universe.The Dirac delta   (k+k ′ ) restricts the average to uncorrelated modes, and () is the spherically averaged power spectrum.  (x) is the brightness temperature (or density) contrast   (x) = ( (x) − T (x))/ T (x), and T (x) is the mean temperature.δ (k) is the three dimensional Fourier transform of   (x), defined by3 : We define the dimensionless PS: where the angular brackets now denote the incoherent average in a spherical shell of width Δ log  = 0.173, and δ *  (k, ) is the conjugate transpose of δ (k, ).

Bispectrum and Skew Spectrum
The bispectrum of the 21cm brightness temperature fluctuations is defined as: 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 k 1 +k 2 +k 3 = 0. Without loss of generality we let k 1 = k, k 2 = q, and k 3 = − (k + q).The SS is the integral of the bispectrum (, , |k + q|) over all possible triangle configurations, for a fixed triangle side : (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: 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 δ 2 (k) below: Similar to Equation 3 we can define the dimensionless SS: 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 Δ log  = 0.173.For consistency and for comparison we use the same bins for calculating the SS and the PS throughout this work.

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.

A Realistic Galaxy Population
We make use of the L210_N4320 dark matter-only simulation of the genesis suite of -body simulations (Power et al. in prep).L210_N4320 has 4320 3 dark matter particles in a cubical volume of side length  = 210ℎ −1 Mpc achieving a halo mass resolution of ∼ 5 × 10 8  ⊙ .
The halo merger trees from L210_N4320 were further 'augmented' to a halo mass resolution of ∼ 3×10 7  ⊙ , the atomic cooling limit at  = 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  = 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 ∼ 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 selfconsistently 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  esc (≤ 1) of the galaxies is a power-law in redshift  (also see Section 3.4.1): where  esc,0 is the escape fraction normalisation and  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  ∼ 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).

IGM Ionization State
meraxes computes the IGM ionization, following the seminumerical code 21cmfast (Mesinger et al. 2011), via an excursionset 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 ∼ 0.2 ℎ −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.
where  b * (, |) is the number of stellar baryons in a sphere of radius  centred at  and redshift ,   = 4000 is the number of UV ionizing photons per baryon (Barkana & Loeb 2007),  atom is the number of neutral H i in the same volume, nrec is the average number of recombinations in the IGM (Sobacchi & Mesinger 2014), and x 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  = 50 Mpc and decrease it successively down to the size of a cell.

21cm Signal
The differential brightness temperature of the 21cm emission from a cloud of H i gas illuminated by CMB radiation of temperature  CMB is given by: where  S is the IGM spin temperature which determines the energy level populations of the H i hyperfine states,   0 is the optical depth,  nl ≡ / ρ − 1 is defined as the evolved Eulerian density contrast ( is the density),  () is the Hubble parameter,   / is the line-ofsight co-moving velocity gradient, and  H i is the neutral fraction.meraxes sources the density and the velocity fields from the N-body simulations and creates self-consistent  S and  H i fields.

Spin Temperature
As can be seen from the Equation (11), the spin temperature  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. S is related to the UV and X-ray emission via: where   and  c are the Wouthuysen-Field coupling (Wouthuysen 1952;Field 1958) and the collisional coupling coefficients respectively. 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).  depends on the local Ly  background flux and closely follows the implementation in Mesinger et al. (2011).  is the 'colour' temperature,  K is the kinetic temperature of the IGM, and we assume   =  K (Field 1959).
The spin temperature field is therefore very sensitive to the  K , which is impacted by X-ray heating.The evolution of the  K depends on the angle-averaged X-ray intensity  (, , ) which is computed as a function of the position , X-ray photon energy , and redshift : where we have integrated the comoving X-ray emissivity  X back along the light cone, and  −  accounts for the probability that an Xray photon emitted at redshift  ′ survives till .We compute  X as a function of the position , X-ray photon energy   =  (1+ ′ )/(1+) at the emitted redshift  ′ : where   /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 ,   /SFR ∝  −   , with   = 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: where  <2 keV /SFR is the soft-band X-ray luminosity per SFR in units of (erg s −1 M −1 ⊙ yr), and  0 fixes the minimum energy for an X-ray photon so that it is not absorbed within the galaxy.

Simulations
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.

Halo Mass Threshold
To explore the impact of the minimum halo mass on the EoR morphology, we modify the  esc (also see Equation 9and Section 3.2) prescription as follows: We run simulations with mass thresholds  thresh = 10 9 and 10 10  ⊙ , and label them mid_ ℎ and high_ ℎ respectively (we point out that our Fiducial simulation contains all haloes down to 3 × 10 7  ⊙ ).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  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_ ℎ and high_ ℎ ) to be the same.This was a deliberate design choice to isolate the impact of the EoR morphology on the 21-cm statistics.

X-Ray Luminosity
We also consider two simulations with a lower and higher X-ray luminosity  <2 keV /SFR = [3.16×10 38, 3.16×10 42 ] ergs s −1 M −1 ⊙ yr as compared to the fiducial value of 3.16 × 10 40 ergs s −1 M −1 ⊙ yr.These two simulations, labelled low_  and high_  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_  , 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.

X-ray Energy Threshold
We explore the impact of the X-ray photon energy threshold  0 by producing two simulations with  0 = 0.2 keV and  0 = 1 keV compared to  0 = 0.5 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).

Fiducial no Spin Temperature
As a comparison we create an additional simulation of the Fiducial signal in the spin temperature saturation limit where   ≫  CMB (from hereon Fid_no  ).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   and ionization field  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.

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 ( xHI ) 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 ( ∼ 0.1 Mpc −1 ) and small ( ∼ 1 Mpc −1 ) spatial scales.

21cm Lightcones
Fig 1 shows a lightcone slice of each simulation as a function redshift.In descending order the panels show the Fiducial, mid_ ℎ , high_ ℎ , low_  , high_  , low_ 0 , and high_ 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   = 0 is typically associated with regions that are ionized.
The same -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 ( ≲ 20) and during reionization ( ≲ 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_  is always in absorption, even during reionization, and the high_  simulation is heavily preheated at high redshift.The low_ 0 simulation has small regions of localised heated gas that appear in emission at  > 15.In contrast, the high_ 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_ ℎ , and high_ ℎ simulations, with the size increasing from Fiducial to high_ ℎ at fixed redshift.We discuss these features in the context of the statistics in the following subsections.

Neutral Fraction and Ionization
In Fig 2 we show the average neutral fraction ( xHI ) 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  = 10 the Fiducial simulation is already partially ionized at the ∼ 5 percent level.We note that the ionization history for the low_ 0 (purple dashed line), the high_ 0 (light purple dash dotted line), and the low_  (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 (  esc,0 ) as the Fiducial simulation (see Table 1).In contrast the high_  (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_  simulation has a ionization fraction of ∼ 10 percent reionization by  = 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_ ℎ (dark blue dashed line) and the high_ ℎ (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_ ℎ and the high_ ℎ 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_ ℎ and high_ ℎ simulations have UV  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_ ℎ and mid_ ℎ , 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.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_  and high_  simulations we refer the reader to this work.For all simulations we see a characteristic absorption feature which occurs when Ly- 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  ∼ 15 (with the exception of the high_  simulation).The timing of the trough depends on Xray heating which eventually drives the IGM into emission (with the exception of the low_  simulation), occurring during reionization at  ≳ 10.This culminates in a peak roughly at the midpoint of reionization in the redshift range of  ∼ 6 − 8.

Mean Brightness Temperature
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 <  < 15.Here the low_ 0 (purple dashed line) undergoes heating earlier than the Fiducial and the high_ 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_ 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.(Wouthuysen 1952;Field 1958), and a peak at  ∼ 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_ 0 simulation and the high_  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.

21cm Power Spectrum
We also include the PS as a function of redshift for the Fid_no  simulation (dashed black line).We only show the Fid_no  simulation up to  = 12, since the approximation is only valid during reionization.We find good agreement between the Fiducial simulation and the Fid_no  simulation during the middle and late periods of reionization at  ∼ 0.1 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  ∼ 1 Mpc −1 since the spin temperature calculation is omitted in the Fid_no  simulation.
The mid_ ℎ and the high_ ℎ simulations in panel (a) have broadly the same PS, differing at most in amplitude during the EoR, with the high_ ℎ simulation peaking earlier than the mid_ ℎ simulation.The mid_ ℎ and high_ ℎ 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  ∼ 10 when reionization begins to become more significant for the Fiducial simulation.
For the low_  and high_  simulations, we find the same features at  ∼ 0.1 Mpc   4 .The high amplitude of the low_  simulation is due to the higher temperature contrast that results from the colder IGM (lower X-ray luminosity).The lower amplitude of the high_  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_  simulation results in the merging of the EoH and EoR peaks.In the high_  simulation we have four peaks, an early peak at  ∼ 22 coincident with the relatively weak absorption trough (Ly- pumping), a peak at  ∼ 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  ∼ 7.In addition to the expected three peaks there is an additional peak at  ∼ 10.This is due to the first ionization sources, with the high_  simulation having ∼ 10 percent ionization by  = 10.
For the high_ 0 case at  ∼ 0.1 Mpc −1 the first peak occurs at  ∼ 14, which is roughly coincident with the absorption trough in the mean temperature brightness in Fig 3 which reaches a minimum at  ≲ 14.The second less prominent peak occurs during the EoR ( ∼ 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- pumping dominates the amplitude of the PS for a longer period.
For the low_ 0 simulation at fixed scale  ∼ 0.1 Mpc −1 we see the three peaked structure.The peak due to Ly- occurs early at  ∼ 17, with the peak during the EoH occurring at  ∼ 12 which coincides with the pockets of emission seen in Fig 1 .The amplitude of the EoH is greater for the low_ 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.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_  which is always in absorption5 .The transition to positive SS for almost all simulations occurs earlier (Δ < 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 ( ∼ 10), the star formation rate increases.The Xray 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.

21cm Skew Spectrum
This transition happens significantly earlier for the low_ 0 simulation (Δ ∼ 2.5).Due to the lower X-ray energy threshold, the low_ 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_ 0 simulation we see a Δ ∼ 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_  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_ ℎ simulation at the largest scale due to the fact that reionization is driven by larger, more biased galaxies.However, the high_  simulation has the highest overall amplitude at all scales in Fig 5 .The efficient X-ray heating in the high_  simulation drives the non-Gaussianities, and results in an earlier reionization, and therefore a peak at ( ∼ 7) compared to the other simulations.
For contrast we calculate the SS for the Fid_no  simulation for all redshifts up to  = 12 (black dashed line).At large scales ( ∼ 0.1 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 (  ≫  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  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  = 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  ∼ 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_ ℎ , high_ ℎ and high_ 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  = 0.1 Mpc −1 for the low_ 0 , the high_  and the low_  simulations we do find some weaker peak and trough features towards the earlier stages of reionization ( ∼ 7.5), when the Fid_no  simulation SS is negative.The absolute change in amplitude as a result of these features are on the order of 10 3 − 10 5 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  , 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.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_  and high_  simulations.

Skew Spectrum and Power Spectrum Co-evolution
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_ ℎ and high_ ℎ 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_ ℎ and high_ ℎ simulations compared to the Fiducial.
For the low_  simulation the SS is always negative, there is also a fairly flat SS from  = 10 to  = 6, with a small peak.These corresponds to the turning points in the curve, similar to the Fiducial simulation.The high_  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_  simulation which are correlated with the two peaks seen between 10 <  < 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.
For the low_ 0 and high_ 0 simulations, the evolution is notably different.For the low_ 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_ 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_ 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_ 0 on other simulations.

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.

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 -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 (), referred to as the normalised SS: In Equation 17Δ 2  and Δ 2  2  are the dimensionless spherically averaged PS and SS, which have units of mK 2 and mK 3 respectively.The subscripts  and  2  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.

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  HI ∼ 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 (()) from Equation 17 This is typically related to the ionization morphology (Zaldarriaga et al. 2004).Notably, the PS and SS are very similar as a function of , 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 ().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 () 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 xHI ∼ 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_ ℎ and high_ ℎ 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 xHI ∼ 0.25 the maxima is significantly larger for the Fiducial simulation relative to the mid_ ℎ and high_ ℎ 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_ ℎ with its much larger escape fraction ionizes a larger volume per ionizing source, preventing the appearance of neutral islands.Thus the high_ ℎ simulation does not exhibit a strong local maxima.
The location of the local minimum for the Fiducial simulation changes with decreasing neutral fraction from  ∼ 0.2 Mpc −1 to ∼ 0.08 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 xHI ∼ 0.5 appears at  ∼ 0.1 Mpc −1 for all simulations.The expected size of the ionizing regions ranges from 20−100 Mpc (Wyithe & Loeb 2004; Figure 9. PS (row one), SS (row two) and normalised SS (row three), of the Fiducial (solid black line), low_ 0 (double dotted dash line), and the high_ 0 (dash dotted line).Each figure from left to right is the xHI ∼ 0.75, 0.5, and 0.25 for each simulation.Zaldarriaga et al. 2004;Lin et al. 2016), with the minima corresponding to size scales from ∼ 31 − 79 Mpc.The local maxima on the other hand appears at  ∼ 1 Mpc −1 for the Fiducial,  ∼ 0.6 Mpc −1 for the mid_ ℎ and  ∼ 0.4 Mpc −1 for the high_ ℎ simulation.These correspond to sizes 6.3 Mpc, 3.8 Mpc, and 1.9 Mpc for the Fiducial, mid_ ℎ and the high_ ℎ simulations respectively.Here the low_  simulation power spectra differs from the Fiducial and high_  at all three ionization states.This difference is driven by the larger absolute temperature difference of the low_  simulation compared to the other two simulations.The flattening of the low_  PS at  < 1 Mpc −1 is still prevalent, since reionization still proceeds the same for the low_  , high_  and Fiducial simulation.The differences at high  are due to the signal temperature on small scales, in the Fiducial and the high_  simulations the signal is in emission, with the high_  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_  simulation which remains in absorption as reionization progresses. 6hen we look at the SS we see a similar trend, with the ex-ception that low_  is negative compared to the Fiducial and high_  scenarios.The shape of low_  relative to the other simulations is broadly mirrored about the -axis, and shows a similar morphology to the PS, as is similar with the Fiducial and the high_  simulations.The normalised SS shows similar features to the Fiducial for both the low_  (albeit negative) and for the high_  .There is a peak/trough at  ∼ 0.1 Mpc −1 , and a subsequent peak/trough at  ∼ 1 Mpc −1 .However, this is less pronounced in the low_  case, and this feature has more of a flattening from  = 1 − 10 Mpc −1 , which is also seen in the PS.Due to the similarity in the ionization morphology between the Fiducial, low_  and the high_  simulations especially at  HI ∼ 0.5, and 0.25, the amplitudes of the higher modes ( > 1 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_ 0 simulation, and high_ 0 is the dashed dotted line.In all cases the SS and the PS of the Fiducial and the low_ 0 simulations are almost identical during reionization, differing at most at  > 1 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_ 0 simulation, increasing  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 Xray 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 ( < 1 Mpc −1 ) and by the gas density on small scales ( ≥ 1 Mpc −1 ) (Watkinson et al. 2018;Majumdar et al. 2018;Ma & Peng 2023).

DETECTABILITY OF THE NORMALISED SKEW SPECTRUM
There is a clear potential advantage to the normalised SS (()) compared to the SS.The trough at  ∼ 0.1 Mpc −1 and the peak at  ∼ 1 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 (), 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: 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 Δ 7 .The number of modes is proportional to the sampled comoving volume of space  used to calculate the spectrum.Equation 18assumes 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(Mondal et al. , 2016)): The non-Gaussianities in the PS are proportional to the Trispectrum  (, ) and  −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.
7 Δ is not a constant since logarithmic bins are typically used to calculate the PS.The logarithmic bin width log Δ = 0.173.

Cosmic Variance of the Power Spectrum and the Skew Spectrum
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 ℎ −1 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 -mode.Subfigures 10a and 10b 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  HI ∼ 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. (2015bMondal et al. ( , 2016)); Greig et al. (2022), where we find an earlier transition to non-Gaussianities in Subfigure 10a in agreement with that seen by Balu et al. (2023).The uncertainties at large scales are dominated by the Gaussian component at ( ≲ 0.3Mpc −1 ) for the PS and ( ≲ 0.2Mpc −1 ) for the SS.This transition is seen at ( ≲ 0.5Mpc −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.

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: (() =   ()/ ()) is the relative error for either the PS or the SS, labelled with the subscripts  and  2  respectively.() 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.
Fig 11 shows the correlation for the xHI ∼ 0.75, 0.5, and the 0.25 Fiducial simulation coeval boxes as a function of .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   and   2  is zero ( (, ) = 0).In this case the relative errors for the PS and SS are equal (  =   2  ) and only depend on the volume  and the shell volume  2 Δ.Equation 20 therefore simplifies to: We use Equation 21 as a model for the Gaussian component of the uncertainties in ().

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 xHI = 0.25, xHI = 0.5, and xHI = 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 (S/N() = 1/  ()).In this case we assume the full relative uncertainty for (), which includes the non-Gaussianities in Equation 20.Notably,   ∝  −1/2 (S/N ∝  1/2 ) therefore, to estimate the S/N for an SKA_LOW observation, we can replace the To determine the SKA_LOW comoving volume for an observation, we first need to know the field of view Ω  and the observing bandwidth Δ.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  ∼ /, where  is the observing wavelength, and  is the station diameter which we assume is 35 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 xHI respectively.Finally assuming an observing bandwidth of ∼ 30 MHz for each neutral fraction, we then determine the comoving volume using Equation 28from Hogg (1999).We then calculate the S/N for each neutral fraction for the Fiducial simula-  (2015a).We see that for all neutral fractions for all spatial scales the S/N > 10, with a max signal to noise of ∼ 300 for xHI = 0.75.Thus, in the absence of thermal noise, foregrounds and systematics, the features present in () should be detectable.MP23 perform their own S/N analysis for the SKA_LOW with the addition of 1000h of thermal noise.They find a S/N of ∼ 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  = 0.1 Mpc −1 and peak at  ≤ 1 Mpc −1 .

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 ( ∼ 0.1 Mpc −1 ) and small ( ∼ 1 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 xHI ∼ 0.75, 0.5 and 0.25.In all simulations we find a local minimum at  ∼ 0.1 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  = 0.1 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.

Sign of the SS During Reionization
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  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 nonlinear.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.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  ∼ 12 ( HI ∼ 0.98).Specifically, we apply varying levels of ionization morphology from xHI ∼ 0.98 to xHI ∼ 0.75 to the  ∼ 12 Fiducial simulation volume, to create pseudo-ionization states at  ∼ 12.
A 2D slice of the Fiducial simulation at  ∼ 12 is shown in Subfig 14a, and in Subfig 14b we show the same slice with the xHI ∼ 0.75 ( ∼ 7.24) ionization field applied.We calculate the differential temperature brightness distribution for each pseudo-ionization state, these are displayed in Subfig 14c, 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 14d.At large scales ( ∼ 0.1 Mpc −1 ) for ionization states of xHI < 0.9 (greater than 10% ionization) the SS is positive, compared to the initial negative SS of the Fiducial simulation at  ∼ 12.This is similar to the results seen in Figure 3 of MP23.

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 observations8 .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.
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.
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.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.21cm differential brightness temperature light cone slices for each simulation volume as a function of redshift from  = 20 − 5.The top slice is the Fiducial simulation (labelled), and the bottom slice is the high_ 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.

Figure 2 .
Figure 2. The ionization history of the average neutral hydrogen (neutral fraction xHI ) IGM calculated for each simulation (labelled).

Figure 3 .
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_  and high_  simulations are taken fromBalu et al. (2023).We see that the mean temperature for these three simulations agree with those shown in Figure8ofBalu et al. (2023), for more details on the low_  and high_  simulations we refer the reader to this work.For all simulations we see a characteristic absorption feature which occurs when Ly- 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  ∼ 15 (with the exception of the high_  simulation).The timing of the trough depends on Xray heating which eventually drives the IGM into emission (with the exception of the low_  simulation), occurring during reionization at  ≳ 10.This culminates in a peak roughly at the midpoint of reionization in the redshift range of  ∼ 6 − 8.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 Fig4).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 <  < 15.Here the low_ 0 (purple dashed line) undergoes heating earlier than the Fiducial and the high_ 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_ 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.

Fig 4
Fig 4 shows the PS calculated at spatial scales  ∼ 0.1 Mpc −1 (panel (a)), and  ∼ 1 Mpc −1 (panel (b)), as a function of redshift for each coeval box, and each simulation.The lines in Fig 4 correspond

Figure 4 .
Figure 4.The PS as a function of redshift for all simulations at a fixed spatial scale.Panel (a) shows the redshift evolution at  ∼ 0.1 Mpc −1 , and panel (b) shows the redshift evolution at  ∼ 1 Mpc −1 .

Figure 5 .
Figure 5.The SS as a function of redshift for all simulations at a fixed spatial scale.Panel (a) shows the redshift evolution at  ∼ 0.1 Mpc −1 , panel (b) shows the redshift evolution at  ∼ 1 Mpc −1 .

Fig 5
Fig 5 shows the SS calculated at  ∼ 0.1 Mpc −1 (panel (a)), and  ∼ 1.0 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_  which is always in absorption 5 .The transition to positive SS for almost all simulations occurs earlier (Δ < 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 ( ∼ 10), the star formation rate increases.The Xray 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

Figure 6 .
Figure 6.The dimensionless SS plotted against the dimensionless PS as a function of redshift for each simulation at a fixed spatial scale  ∼ 0.1 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  = 20, and the left pointing triangle indicating the endpoint at  = 5.

Fig 6
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_ ℎ , and high_ ℎ in panel (a), the low_  and high_  in panel (b), and the low_ 0 , and high_ 0 simulations in

Figure 7 .
Figure 7. PS (row one), SS (row two) and normalised SS (row three), of the Fiducial (solid black line), mid_ ℎ (double dotted dash line), and the high_ ℎ (dash dotted line).Each figure from left to right is the xHI ∼ 0.75, 0.5, and 0.25 for each simulation.

Figure 8 .
Figure 8.PS (row one), SS (row two) and normalised SS (row three), of the Fiducial (solid black line), low_  (double dotted dash line), and the high_  (dash dotted line).Each figure from left to right is the xHI ∼ 0.75, 0.5, and 0.25 for each simulation.
for each simulation.() 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 () (third row) as a function of spatial scale, with the different neutral fractions of ∼ 0.75, ∼ 0.5 and ∼ 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_ ℎ the double dot dashed line, and the high_ ℎ is the dashed dot line.For both the PS and the SS we see a flattening of the spectra ( < 1 Mpc −1 ) as reionization progresses.
Fig 8 has the same layout as Fig 7, here low_  is the double dot dashed line, and the dash dot line is the high_  simulation.

Figure 10 .
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).

Figure 11 .
Figure 11.Correlation coefficient between the Fiducial PS and SS for the  HI ∼ 0.75 (solid), the 0.5 (dash dotted line), and the 0.25 (double dot dashed line) neutral fractions.

Figure 12 .
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 xHI = 0.5.

Figure 13 .
Figure 13.Signal to noise ratio of  ( ) for the Fiducial simulation scaled by the expected SKA_LOW observing comoving volume.The solid curve corresponds to the Fiducial neutral fractions  HI ∼ 0.75, 0.5 for the dash dotted line, and 0.25 for the double dot dashed line.
Figure 14.Subfigure (a) shows a 2D slice through the meraxes  ∼ 12 coeval box, with  HI ∼ 0.98.Subfigure (b) shows the same slice, with a the  HI ∼ 0.75 ionization field taken from  ∼ 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.

Table 1 .
Astrophysical parameter summary for the seven meraxes simulations.See text for details.