The 21-cm bispectrum from neutral hydrogen islands at z < 6

Spatial variations in the Lyman-α forest opacity at z < 6 seem to require a late end to cosmic reionization. In this picture, the universe contains neutral hydrogen ‘islands’ of up to 100 cMpc /h in extent down to redshifts as low as z ∼ 5 . 3 . This delayed end to reionization also seems to be corroborated by various other observables. An implication of this scenario is that the power spectrum of the cosmological 21-cm signal at z < 6 is enhanced relative to conventional reionization models by orders of magnitude. However, these neutral hydrogen islands are also predicted to be at the locations of the deepest voids in the cosmological large-scale structure. As a result, the distribution of the 21-cm signal from them is highly non-Gaussian. We derive the 21-cm bispectrum signal from these regions using high-dynamic-range radiative transfer simulations of reionization. We find that relative to conventional models in which reionization is complete at z > 6 , our model has a significantly larger value of the 21-cm bispectrum. The neutral islands also imprint a feature in the isosceles bispectrum at a characteristic scale of ∼ 1 cMpc − 1 . We also study the 21-cm bispectrum for general triangle configuration by defining a triangle index. It should be possible to detect the 21-cm bispectrum signal at ν ≳ 200 MHz using SKA1-LOW for 1080 hours of observation, assuming optimistic foreground removal


INTRODUCTION
Observations of the Lyman-α forest point to a late end of reionization (Kulkarni et al. 2019;Keating et al. 2020;Bosman et al. 2022).In our previous work, we explored the implications of a reionization model that agrees with these observational constraints at redshifts 5-8 for the 21-cm power spectrum (Raste et al. 2021).We found that given the late end of reionization, the power spectrum of the 21-cm signal at redshifts z = 5-6 is orders of magnitude higher than previous estimates.This signal should be detectable by hera and ska1-low in ∼ 1000 hours of observations, assuming optimistic foreground subtraction (Raste et al. 2021).
However, models predict that the large islands of neutral hydrogen that persist till redshift z ∼ 5.5 in our reionization models are in the deepest density voids in the universe.As a result, the 21-cm signal from them should be significantly non-Gaussian, which should lead to a large bispectrum signal.Furthermore, these neutral islands have highly irregular shapes that might hold clues about the galaxies that contributed to reionization.While the power spectrum is sensi-⋆ E-mail: janakee@ncra.tifr.res.intive to the size of the ionized regions, the bispectrum is sensitive to their shapes, which makes it a promising tool to study the topology of reionization (Hutter et al. 2020).Unlike the power spectrum, the value of the bispectrum signal for different triangle configurations can be negative as well as positive and this can encode information on various features of the reionization.Simulations have shown that the modelling parameters of density, ionization, X-ray heating and Lyman-α coupling drive the non-Gaussianity at various scales.These processes determine the shape, magnitude, peak and sign of the 21-cm bispectrum as a function of redshift (Shimabukuro et al. 2016;Majumdar et al. 2018;Watkinson et al. 2019;Hutter et al. 2020;Kamran et al. 2021b;Ma et al. 2021;Shaw et al. 2021;Kamran et al. 2021a).The 21-cm bispectrum is also a function of the redshift-space distortions and light-cone anisotropy (Bharadwaj et al. 2020;Majumdar et al. 2020;Kamran et al. 2021b;Shaw et al. 2021;Mondal et al. 2021).It has been consistently shown by multiple authors that observing the 21-cm bispectrum together with the power spectrum can reduce the inferred uncertainty on reionization parameters (Shimabukuro et al. 2016(Shimabukuro et al. , 2017;;Watkinson et al. 2022;Tiwari et al. 2022).
The 21-cm bispectrum signal can be observed by correlat- Figure 1.The bottom panel shows the 21-cm brightness temperature, ∆T b , in a reionization model that is consistent with the Lyα forest at z ∼ 6.Large neutral hydrogen 'islands' are seen at z < 6.In order to understand the 21-cm bispectrum signal from these neutral hydrogen islands, we contrast this reionization model with a more conventional one in which reionization finishes early, at z > 6.This model is shown in the top panel.
In this paper we compute the bispectrum of the ionized hydrogen fraction, the gas density and the 21-cm brightness temperature for our late reionization model.In Section 2, we briefly describe our simulation.We also discuss in this section the calculation of bispectra using the fast FFT code BiFFT presented by Watkinson et al. (2017) and various normalisations of the bispectrum.We present our results in Section 3 for equilateral, isosceles and scalene triangle configurations.Finally, we discuss the prospects of detecting the 21-cm bispectrum with ska1-low in Section 4, and conclude in Section 5.

METHODS
We have discussed our simulation in detail in Kulkarni et al. (2019) and Raste et al. (2021).Here we repeat only the essential details.To obtain the gas density and velocity fields, we have used the p-gadget-3 code, a modified version of the gadget-2 code (Springel et al. 2001;Springel 2005).This simulation is similar to the simulations from the Sherwood Simulation Suite (Bolton et al. 2017) with their 160-2048 initial conditions, containing 2048 3 gas and dark matter particles and 160 h −1 cMpc box length with periodic boundary conditions.For further processing, the gas density is gridded by projecting the smooth particle hydrodynamic (SPH) kernel in our simulation onto a Cartesian grid of size 2048 3 .This gives a grid resolution of 78.125 h −1 ckpc.The ionization and temperature fields are computed using the aton code (Aubert & Teyssier 2008, 2010), which solves the radiative transfer equation by using the M1 approximation for the first moment (Aubert & Teyssier 2008;Levermore 1984;González et al. 2008).
We calculate the differential brightness temperature (∆T b ) box using the density, ionization, and peculiar velocity boxes and assuming TS ≫ TCMB, We take the z-axis of the simulation box as the line-of-sight direction to calculate the peculiar velocity gradient dvp/ds.The above expression assumes dvp/ds ≪ H, which is generally a valid assumption.However, for the very few simulation cells in our computation that have that have dvp/ds ∼ H, we follow the standard practice and enforce a ceiling of |dvp/ds| < 0.5H(z) (Santos et al. 2010;Mesinger et al. 2011;Mao et al. 2012).
The reionization in our 'late reionization' simulation model ends at redshift at z ∼ 5.3, and the midpoint of reionization occurs at redshift z ∼ 7.1.We also study an 'early reionization' model, in which the evolution of the volume-averaged ionized hydrogen fraction is calibrated to match the evolution in the Haardt & Madau (2012) model of reionization.In this model, reionization is complete at z ∼ 6.7.The two simulations are identical in all aspects apart from the source emissivity.Figure 1 compares the 21-cm brightness temperature in the two models.

Computing the Bispectrum
The bispectrum B of a field F (r) is defined as, where F (k) is the Fourier transform of F (r).The Dirac delta δD term requires that the wavevectors k1, k2 and k3 form a closed triangle in the Fourier space.
We follow Scoccimarro (2015), Sefusatti et al. (2016), and Watkinson et al. (2017) to efficiently compute the bispectrum without using multiple nested loops through the Fourier box 1 .This algorithm is described in detail by Watkinson et al. (2017) and is implemented by these authors in their publicly available code, BiFFT 2 .We calculate bispectrum using a modified Python version of BiFFT.
1 In this work, we do not subtract the mean from the simulation box before calculating bispectra.All the information about the mean value is located only in the real part of the k = 0 mode of the Fourier box.This mode is not used while calculating power spectrum or bispectrum.We confirm that our results do not change by subtracting the mean from the simulation box. 2 https://bitbucket.org/caw11/bifft/

Triangle Configuration
For any triangle formed by the wavevectors k1, k2 and k3, where θ is the angle between k1 and k2 arms of the triangle.

Density Bispectrum
Figure 2 shows the evolution of isosceles bispectra of the gas density from redshift z = 9.02 to 5.26.The isosceles triangle configuration is used for four representative values of k1 = 0.2, 0.5, 0.75, and 1 cMpc −1 .For each value of k1 (= k2), the figure shows a range of values of k3 available in the simulation box, between 0.08 and 2 cMpc −1 , depending on k1.We have normalised the bispectra using Equation 5 and the box has been reduced to resolution of 256 3 for computa-tional ease (see Appendix A for a comparison with results from the higher resolution box).The normalised density bispectrum is of the order of a few units, and grows with time due to an increase in the non-Gaussianity induced by structure formation.The normalised bispectrum also increases in amplitude with increasing k1, as the small scales have more non-Gaussianity compared with larger scales.

Neutral Hydrogen Fraction Bispectrum
Figures 3 and 4 show, respectively, the unnormalised and normalised bispectrum of the neutral hydrogen fraction at various redshifts between z = 9.02 and 5.26, at k1 = 0.2, 0.5, 0.75 and 1 cMpc −1 .Each of the two figures show the bispectrum for both of our reionization models in separate panels.Apart from the difference in the redshift of reionization, the bispectrum in the two reionization models are qualitatively similar, only shifted in redshift (time).
Recall that cos θ = 0.5 denotes the equilateral triangle configuration.Values of cos θ from −1 to 0.5 correspond to the stretched limit, and those from 0.5 to 1 correspond to the squeezed limit.A large positive value of the equilateral bispectrum indicates an overabundance of roughly spherical structures of higher-than-average values of xHI (neutral regions) embedded in the background of lower-than-average values of xHI (ionized regions) (Lewis 2011;Hutter et al. 2020).A large negative value indicates an over-abundance of below-average structures embedded in the above-average background.This allows us to interpret the evolution that we see in Figure 3.
At high redshifts, the xHI distribution has 'holes' of belowaverage values (the ionized regions), yielding a negative values for almost all k modes and triangle configurations.However, the bispectra at very small scale (large k) stretched limit triangles are positive.These triangles correspond to overabundance of small-scale above-average filamentary structure.This is perhaps the small scale neutral 'valleys' between spherical ionized bubbles.
With the progress of reionization, the ionized regions become larger and start merging.As reionization crosses the half-way point, the distribution of xHI now has 'islands' of above-average values, yielding a positive value of the bispectrum.At the mid-point of reionization, even if roughly half of the volume is occupied by ionized IGM and half by neutral, the shapes of ionized regions are roughly spherical compared to neutral regions, which have more irregular shapes.Therefore, we see that the stretched and squeezed limit bispectrum start becoming positive at lower redshifts, however small scale (large k) equilateral bispectra still remain negative.They are signature of small scale spherical ionized regions embedded in neutral IGM.This signature persists even during the later half of reionization.At large scales (small k), on the other hand, the bispectrum is positive for all configurations during the later part of reionization.This suggests that, on large scales, now there are above-average neutral regions embedded in below-average ionized IGM.
In the stretched limit, cos θ → −1, the bispectrum measures the probability of the neutral hydrogen pancakes that demarcate ionized bubbles.Such pancake-like boundary surfaces are the over-abundant structure at all redshifts.Consequently, the bispectrum for cos θ → −1 is positive at all redshifts after the start of reionization.
In the squeezed limit, cos θ → 1, a positive value of the bispectrum indicates an overabundance of small-scale positive perturbations in the xHI distribution and a large-scale modulation of these perturbations.As soon as the ionized regions grow to a reasonable size, the distribution of xHI inside the large ionized regions is trivially uniform, but that in the leftover neutral regions has small-scale perturbations due to smaller ionized bubbles.The bispectrum turns positive when this happens.
At the end of reionization, at redshifts of z ∼ 5.4, the bispectrum shows complex features that are sensitive to the morphology of the residual xHI islands.In Figure 4, the normalisation accentuates the features in the bispectrum that we see in Figure 3.This is most pronounced for the oscillating features at redshift z ∼ 5.4.These oscillations are a qualitatively distinct signature that seems to mark the end phases of reionization in both early and late reionization models, and so could be a useful smoking gun for identifying the redshift of reionization from observations.

21-cm Bispectrum
Figures 5 and 6 respectively show the unnormalised and normalised isosceles bispectrum of the 21-cm brightness temperature for the same redshift range and the k values as in Figures 3 and 4. For approximately equilateral triangles, the bispectrum is negative in the early stages of reionization, changes sign when reionization is half complete, and then is positive at lower redshift, before settling to zero in the post-reionization era.The bispectrum is positive at almost all redshifts in the stretched limit.Similarly, the bispectrum in the squeezed limit also stays positive at all but the very early stages of reionization.At the end of reionization, for z ∼ 5.4, the bispectrum shows a set of complex features that map to the 21-cm brightness structure of the residual neutral hydrogen islands in the voids.
It is noteworthy that, similar to the 21-cm power spectrum explored in our previous work (Raste et al. 2021), the predicted bispectrum is very large at z < 6 in the late reionization model, while the bispectrum at these redshifts is zero in the fiducial early reionization model.This is due to the persistence of neutral hydrogen islands at these redshift in the late reionization model.This signal in the 21-cm bispectrum is directly induced by the opaque regions seen in the Lyα forest at these redshifts.
In the absence of spin temperature fluctuations, most of the brightness temperature bispectrum is induced by the fluctuations in xHI in the range of redshifts and wavenumbers con-  B1).
sidered here.Therefore, we see that brightness temperature bispectrum follows a very similar trend as the bispectrum of the neutral hydrogen fraction.To study this effect more quantitatively, we break down our ∆T b bispectra in various autoand cross-bispectra components of density and neutral fraction.Ignoring the redshift space distortion, the 21-cm signal at any redshift can be written as, where, T0(z) is the base 21-cm signal at redshift z and there are spatial fluctuations due to density (δD) and neutral H i fraction (xHI).The average 21-cm signal is, And its fluctuation is, Defining δHI = (xHI − ⟨xHI⟩)/⟨xHI⟩ and δD,HI = (δDxHI − ⟨δDxHI⟩)/⟨δDxHI⟩, we have, Notice that the density bispectrum does not contribute directly to the ∆T b bispectra.We should also emphasise that while calculating cross bispectra, the order of the fields is important in a non-equilateral configuration.For example, ⟨δDδHIδHI⟩ might not be the same as ⟨δHIδDδHI⟩ if k1 ̸ = k2.
In Figure 7 we show each of the cross bispectra (unnormalised) component, multiplied by its respective weight, from Eq 9 for the equilateral bispectrum configuration (left) and the isosceles configuration at k1 = 0.5 cMpc −1 , for redshift z = 5.95.The gray dashed curve shows the weighted sum of these individual components.Comparing it with the ∆T b bispectra, we see that this predicted bispectrum is slightly different from the computed bispectra, since we have ignored  here the effect of velocity fluctuations in Eq 6.But this effect is small.We also note that the bispectra of various cross components fluctuate around zero a lot more than the auto bispectra.These fluctuations occur where the non-Gaussianity is close to 0. However, the sum of these cross components do not show these fluctuations, as their effects average out.Finally, note that the ∆T b bispectra have shapes very similar to the neutral fraction bispectra.Hence, in absence of spin temperature fluctuations, the neutral fraction fluctuations dominate the 21-cm fluctuations.
Other than a few k-modes, which show the positive fluctuations, the equilateral bispectrum of ∆T b is negative at intermediate k-modes and postive at very large and very small k-modes.Also notice that for the equilateral configurations ⟨δHIδHIδD,HI⟩ and ⟨δHIδD,HIδHI⟩, as well as ⟨δD,HIδD,HIδHI⟩ and ⟨δHIδD,HIδD,HI⟩ have very similar shapes, other than few fluctuations.However, their shapes are very different for the isosceles configurations.

Bispectrum for generic triangle configurations
Moving from the isosceles triangle to more general triangle configurations, in Figure 8, we show the evolution of the normalized bispectra for available triangle configurations for the 256 3 resolution cube of neutral fraction (top panels) and brightness temperature (bottom panels), in the late (right) and early (left) reionization models.
The mapping between the triangle index and k values is as follows: Of the (k1, k2, k3) triplet, two values are taken from array (0.1, 0.2, 0.5, 0.75, 1.0 cMpc −1 ), with k1 ≥ k2.We choose the third k value with θ/π from the array [0.01, 0.05, 0.1, 0.2, 0.33, 0.4, 0.5, 0.6, 0.7, 0.85, 0.95].Then, (k1, k2, k3) triplets are sorted in increasing order (ka, k b , kc) with ka ≥ k b ≥ kc.The triangle index of a given triplet is its rank in the sorted sequence, with ka moving fastest.We tabulate the mapping between the triangle index and the (k1, k2, k3) triplets in Appendix B. Indices > 100 roughly correspond to k > 0.5 cMpc −1 .The fluctuations at higher triangle indices are due to the triangle configuration fluctuating between stretched limit triangles, which correspond to the large positive bispectra values and other configurations, including equilateral, which correspond to negative bispectra values.
We can see that for early stages of reionization, the bispectra of neutral fraction for various triangle configurations is negative.With the evolution of reionization, they become positive.Towards the end of reionization, the amplitude of the fluctuations increase.The post-reionization bispectra have small amplitude for small triangles (small k-modes, large scale), but they become larger for larger triangle configurations (large k-modes, small scale).The brightness temperature bispectra have very similar behaviour, however, their amplitude is slightly more positive, reflecting the effect of density bispectra.The post-reionization normalised ∆T b bispectra are usually larger than the reionization bispectra at 2.0 1.8 1.6 1.4 1.2 0.8 0.2 Figure 10.Unnormalised isosceles bispectra of x HI (top) and ∆T b (bottom) for the late reionization model at k = 1 cMpc −1 for boxes 0 to 5 and redshifts 7.14, 5.95, 5.41 and 5.26 from left to right.Boxes 0 to 3 have very slight differences, which suggests that perhaps, the bispectra are more dependent on the neutral vs ionized state of the IGM, and the partially ionized regions do not affect the bispectra significantly.Box 4 has large, flat, positive ∆T b bispectra at all redshifts, which matches with the post-reionization bispectra in the last panel.In this box, we have removed the effect of neutral regions.Therefore, the post-EoR bispectra is an effect of residual neutral fraction, as expected.∆T b bispectrum of Box 5 is in essence the gas density bispectrum.
SKA1-LOW, 1080 hours Late Reionization 1.0 0.9 0.8 0.7 0.6 0.4 0.1 Figure 11.ska1-low sensitivity (red curves) for 1080 hours of tracking mode observation with optimistic foreground removal at redshift z = 5.96 (ν = 204 MHz) and k 1 = 0.2, 0.5, 0.75 and 1.0 cMpc −1 from left to right.We compare it with our late reionization model (black curves).The positive and negative sensitivity curves are the 1-σ upper and lower bounds for the real part of noise bispectra.
all k-mode triangles and do not show strong fluctuations with triangle index.

Untangling the bispectrum
To understand which features of our simulation box correspond to which details in the 21-cm bispectrum, we take our original ionization fraction box (Box 0) and construct several different modified boxes using following prescription: In Figure 9, we show slices of brightness temperature boxes created using these modified ionization boxes at redshift z = 5.96.We show unnormalised isosceles bispectra for these boxes at k1 = 1 cMpc −1 , for various redshifts in Figure 10.
Boxes 0 to 3 have very slight differences.Specifically, we can see that box 0 and box 3 have very similar shape and magnitude for almost all redshifts and almost all k-modes.This suggests that perhaps, the bispectra are more dependent on the neutral vs ionized state of the IGM, and the partially ionized regions do not affect the bispectra significantly.However, when we remove these ionized shapes and regions from the simulation box in Box 4, the shape of the bispectrum changes completely.Boxes 4 and 5 are mostly ionized and neutral boxes, respectively, with only some partially ionized regions left.Boxes 4 and 5 provide useful extreme cases with which to compare and contrast boxes 0-3.Box 5 represents a nearly completely neutral volume; only small partially ionized regions exist in an otherwise neutral IGM.Similarly, box 4 represents a nearly completely ionized volume; only small partially neutral regions are present in an otherwise ionized box.A comparison of such configurations with boxes 0-3 isolates contribution to the bispectrum of the most visible component of boxes 0-3, namely the large residual neutral regions.Indeed, the bispectra of boxes 4 and 5 are qualitatively different from the bispectra in boxes 0-3.The ionization fraction bispectra for Box 4 and Box 5 have similar shape with the sign inverted, as we see in top panels of Figure 10.The bispectra also have a smaller amplitude.The 21-cm bispectra are positive in boxes 4 and 5.In box 5, the bispectrum is set by density fluctuations.These do not evolve significantly from redshift 7 to 5, leading to the almost negligible change with redshift in the bispectrum in the bottom panels of Figure 10.Box 4 on the other hand shows a bispectrum that is set by the small number of partially ionized cells.The resultant bispectra matches the post-reionization values that we see in Figure 5.The redshift evolution seem in the 21cm bispectrum in box 4 is due to evolution in the number of partially neutral cells, which affect the 21-cm in spite of their small volume fraction as the rest of the volume has zero brightness.The flatness of the bispectra in boxes 4 and 5 clearly highlight the features introduced in the bispectrum by the neutral hydrogen regions required by the Lyα Forest.In future 21-cm measurements, it would be valuable to correlate such features with the Lyα Forest data.

PROSPECTS OF DETECTION
For a preliminary study of the prospects of detecting the 21-cm bispectrum modelled above, we use the 21cmSense3 (Pober et al. 2013(Pober et al. , 2014) ) and PyObs21 (Watkinson et al. 2022) codes to model the bispectrum covariance induced by instrument noise for ska1-low.
Our assumed telescope parameters are the same as in Raste et al. (2021).We have used the core configuration with 224 elements, with each element having diameter of 38 m and a field of view of about 12.5 deg 2 at 150 MHz (de Lera Acedo et al. 2020).The shortest and longest baselines for this configuration are 35.1 m and 887 m, respectively, resulting in an angular resolution of 10 ′ at 150 MHz.We assume a tracking mode of 6 hr, at z ∼ 6 (νo = 204 MHz), for 180 days.We include all k-modes in our calculation (using the 'opt' foreground mode, i.e., assuming that the foregrounds are some-how removed from the data).Using these noise estimates, we generate 100 noise boxes of 128 3 resolution elements with different seeds and calculate their bispectra.In Figure 11, we show the variance of these noise bispectra and compare it with the signal at redshift z = 5.96 for various k-modes4 .We see that the bispectrum sensitivity is better at lower kmodes (larger scale).These results suggest that ska1-low should be able to detect the 21-cm isosceles bispectra at small k-modes with ∼ 1000 hours of tracking mode observation, assuming optimistic foreground removal.These observations would help us understand the history of reionization, as well as the geometry of neutral islands at the end of reionization.

CONCLUSIONS
We have computed the 21-cm bispectrum in simulations of cosmic reionization that are consistent with the observed Lyman-α forest at z > 5. Our findings can be summarised as follows: • A late end of reionization causes the 21-cm bispectrum and the neutral hydrogen fraction bispectrum to have large values at z < 6 (frequencies greater than 200 MHz).This is in contrast to traditional reionization models, which typically predict zero bispectrum and power spectrum at these frequencies.
• The neutral fraction bispectrum is negative at all scales during the early stages of reionization.At later stages, the bispectrum starts to become positive at small scales (large k).Towards the end of reionization, the bispectra are positive at all scales for squeezed and stretched limit triangles.But they remain negative around the equilateral configuration at small scales (large k).This suggests that most of the non-spherical, "elongated" features have the geometry of above-average regions (neutral islands) embedded in the below-average background (ionized regions), but there are still some small spherical below-average regions (ionization bubbles) embedded within above-average regions (neutral islands).
• The 21-cm bispectrum follows the trends seen in the neutral fraction bispectrum.
• For generic triangle configurations, the normalised bispectra of neutral fraction are negative at high redshifts.They then turn positive with the progress of reionization.The ∆T b bispectra show very similar behaviour, however, their amplitude is slightly more positive, reflecting the effect of the density bispectra.This effect of density is most readily observed in post-reionization brightness temperature bispectra, which are usually larger than the reionization bispectra at all k-mode triangles and do not show strong fluctuations with triangle index.
• Partially ionised regions do not affect the shape of xHI or 21-cm bispectra significantly.However, removing large ionised regions or neutral islands will completely change the shape, amplitude and sign of bispectra at any redshift.
• With about a 1,000 hr of tracking mode observation and optimistic foreground removal, ska1-low should be able to detect the bispectrum at z ∼ 6 for small k-modes (large scale).
Overall, this work adds to the realization that statistics beyond the power spectrum can be very useful to understand reionization history from the future 21-cm observations.Similar to our previous work on the 21-cm power spectrum from the late reionization model, this work also highlights the importance of relatively higher frequency corresponding to redshifts z < 6 for the late reionization.Combining these 21-cm bispectrum and power spectrum results with the James Webb Space Telescope (JWST) and Nancy Grace Roman Space Telescope should lead to further insights into reionization physics.
APPENDIX A: RESOLUTION TESTS Our simulation boxes of the ionization, density and velocity have resolution of 2048 3 pixels, therefore, our brightness temperature box also has the same resolution.However, it is difficult to calculate bispectrum for the full resolution boxes.Therefore, in our work, we have presented bispectrum for 256 3 pixel box.We reduce the box resolution by averaging neighboring n 3 pixels into one large pixels, where n > 1.In Figures A1, we show unnormalised and normalised isosceles bispectra at k1 = 0.5 cMpc −1 and power spectra for 512 3 , 256 3 and 128 3 pixels resolutions at z = 5.96 .
As expected, the small scale (large k) power spectra diverge with resolution, as the fluctuations on these scales are averaged out while lowering the resolution.A similar trend can be seen for the unnormalised bispectra where the large k modes show larger divergence for 128 3 box.However, the normalised bispectra show very little fluctuations with box resolution as the effect of resolution on bispectrum and power spectrum cancel out due to normalisation.Note that the peaks in 128 3 (red) box in bottom middle and right panels occur only at kmodes where the bispectra have very small amplitude.Overall, the results presented in our work are robust with box resolution at relevant scales.

Figure 2 .
Figure2.The normalised isosceles gas density bispectrum from redshift z = 9.02 to 5.26 for k 1 = 0.2 cMpc −1 , 0.5 cMpc −1 , 0.75 cMpc −1 and 1.0 cMpc −1 (left to right).Amplitude of density bispectrum grows with time due to the increasing non-Gaussianity with the formation of structures.This amplitude also increases with k 1 , as the small scales have more non-Gaussianity compared with larger scales.

Figure 7 .
Figure7.Unnormalised cross bispectra components, multiplied by their respective weights, from Eq 9 for equilateral (left) and isosceles (right, at k 1 = 0.5 cMpc −1 ) bispectrum configuration at z = 5.95 for the late reionization model.We compare the weighted sum (gray dashed curves) of these components with the ∆T b bispectra (thick black curves).

Figure 8 .
Figure 8. Normalized bispectrum for neutral fraction (top panel) and brightness temperature (bottom panel) for the late (right panel) and early (left panel) reionization models from redshift z = 9.02 to 5.26 for all triangle configurations (see TableB1).

Figure 9 .
Figure 9. ∆T b boxes 0 to 5 for redshift z = 5.96 for the late reionization model.These boxes were computed after modifying the ionization field by hand.Box 0 is the original simulations box, whereas in Box 1, x HII < 0.5 regions are set to x HII = 0, and in Box 2, x HII ≥ 0.5 regions are set to x HII = 1.Box 3 has both these approximations, essentially converting the box to a binary field of 0 and 1.Finally, in Box 4 x HII ≤ 0.5 regions are set to x HII = 1 and in Box 5 x HII ≥ 0.5 regions are set to x HII = 0, which respectively removes neutral and ionized regions from the simulation box.
5 is set to xHII = • Box 3: Both of the above (essentially converting the box to a binary field of 0 and 1) • Box 4: xHII ≤ 0.5 is set to xHII = • Box 5: xHII ≥ 0.5 is set to xHII =

Table B1 .
Values of (k 1 , k 2 , k 3 ) triplets for all triangle indices shown in Figure8.All k values are in units of cMpc −1 .