Studying the morphology of reionisation with the triangle correlation function of phases

TL;DR

This paper introduces the triangle correlation function (TCF), a new phase-based statistical tool derived from the three-point correlation function, to analyze ionized region scales during reionization from 21cm data, offering advantages over traditional power spectrum methods.

Contribution

The paper presents the TCF as a novel phase-sensitive statistic for probing ionized bubble scales, demonstrating its effectiveness and robustness on simulated and semi-realistic reionization fields.

Findings

TCF peaks at characteristic bubble radius in simple models

Robustness of TCF under observational noise and resolution effects

Ability to identify multiple characteristic scales in complex fields

Abstract

We present a new statistical tool, called the triangle correlation function (TCF), inspired by the earlier work of Obreschkow et al. It is derived from the three-point correlation function and aims to probe the characteristic scale of ionized regions during the epoch of reionization from 21cm interferometric observations. Unlike most works, which focus on power spectrum, i.e. amplitude information, our statistic is based on the information we can extract from the phases of the Fourier transform of the ionization field. In this perspective, it may benefit from the well-known interferometric concept of closure phases. We find that this statistical estimator performs very well on simple ionization fields. For example, with well-defined fully ionized discs, there is a peaking scale, which we can relate to the radius of the ionized bubbles. We explore the robustness of the TCF when observational effects such as angular resolution and noise are considered. We also get interesting results on fields generated by more elaborate simulations such as 21CMFAST. Although the variety of sources and ionized morphologies in the early stages of the process make its interpretation more challenging, the nature of the signal can tell us about the stage of reionization. Finally, and in contrast to other bubble size distribution algorithms, we show that the TCF can resolve two different characteristic scales in a given map.

Tables

Table 1. Table 1: Peaking scale (in Mpc) for different methods to obtain the bubble size distribution of realisations of the 21CMFAST simulation at different redshifts.

$z$	$x_{\text{H }\text{ii}}$	TCF	RMFP	SPA
13.0	0.036	< 8.9	2.01	< 5.9
12.0	0.062	< 10.8	2.34	< 7.4
11.0	0.010	< 11.2	3.12	< 8.9
6.4	0.971	6.2 & 20.2	7.0	13.9
6.2	0.994	7.7	7.0	9.4

Table 2. Table 2: Peaking scales (Mpc) for a toy model with different types of observational effects: smoothing due to angular resolution and noise for an integration time of 1000 ​ h 1000 h 1000~{}\mathrm{h} . Actual radius size is 7.8 ​ Mpc 7.8 Mpc 7.8~{}\mathrm{Mpc} .

	Clean	Smoothed	Smoothing	Smoothing
	signal	signal	+ 1000h	scale
LOFAR	6.7	6.6	7.4	4.5
SKA1-Low core	6.7	8.2	–	17.0
SKA1-Low central	6.7	6.7	6.9	5.3

Equations

$$x(\mathbf{r})=\sum \hat{x}(\mathbf{k})e{}^{i\mathbf{k}\cdot \mathbf{r}}.$$

$$\widehat{ϵ}(\mathbf{k})=\frac{\hat{x}(\mathbf{k})}{|\hat{x}(\mathbf{k})|}={\mathrm{e}}^{i{\Phi }_{\mathbf{k}}}.$$

$$\delta T_{\mathrm{b}}$$

$$\simeq 27\mathrm{mK}{x}_{\text{Hi}}\left(1+{\delta }_b\right)\frac{{\Omega }_{\mathrm{b}}{h}^2}{0.023}\sqrt{\frac{0.15}{{\Omega }_{\mathrm{m}}{h}^2}}\sqrt{\frac{1+z}{10}}$$

$${\Xi }_n({\mathbf{r}}_1,...,{\mathbf{r}}_{n-1})=\frac{1}{V}{\int }_V{\mathrm{d}}^{D}r\prod _{j=1}^{n-1}x(\mathbf{r}+{\mathbf{r}}_j),$$

$$P_n\left({\mathbf{k}}_1,\dots {\mathbf{k}}_{n-1}\right)=\hat{x}({\mathbf{k}}_1)\dots \hat{x}({\mathbf{k}}_{n-1})\hat{x}\left(-\Sigma {\mathbf{k}}_j\right)$$

$${\Xi }_n\left({\mathbf{r}}_1,\dots {\mathbf{r}}_{n-1}\right)={\left[\frac{V}{(2\pi {)}^D}\right]}^{n-1}\int {\mathrm{d}}^D{k}_1{\mathrm{e}}^{i{\mathbf{k}}_1\cdot {\mathbf{r}}_1}\dots$$

$$\times \int {\mathrm{d}}^D{k}_{n-1}{\mathrm{e}}^{i{\mathbf{k}}_{n-1}\cdot {\mathbf{r}}_{n-1}}{P}_n\left({\mathbf{k}}_1,\dots {\mathbf{k}}_{n-1}\right).$$

$$\mathcal{P}(\mathbf{k})=\hat{x}(\mathbf{k})\hat{x}(-\mathbf{k})=|\hat{x}(\mathbf{k}){|}^2,$$

$$\mathcal{B}(\mathbf{k},\mathbf{q})=\hat{x}(\mathbf{k})\hat{x}(\mathbf{q})\hat{x}(-\mathbf{k}-\mathbf{q}),$$

$${\Xi }_3\left(\mathbf{r},\mathbf{s}\right)=\frac{V^2}{(2\pi {)}^{2D}}\iint {\mathrm{d}}^{D}k{\mathrm{d}}^{D}q{\mathrm{e}}^{i(\mathbf{k}\cdot \mathbf{r}+\mathbf{q}\cdot \mathbf{s})}\mathcal{B}(\mathbf{k},\mathbf{q}).$$

$$\{\begin{array}{rl}& s_x=\frac{1}{2}{r}_x-\frac{\sqrt{3}}{2}{r}_y\\ & s_y=\frac{\sqrt{3}}{2}{r}_x+\frac{1}{2}{r}_y\\ & s_z=r_z\end{array}$$

$$\mathbf{p}=\left(\begin{array}{c}{k}_x+\frac{1}{2}{q}_x+\frac{\sqrt{3}}{2}{q}_y\\ k_y-\frac{\sqrt{3}}{2}{q}_x+\frac{1}{2}{q}_y\\ k_z+q_z\end{array}\right),$$

$${\Xi }_3^{\ast }(\mathbf{r})=\frac{V^2}{(2\pi {)}^{2D}}\iint {\mathrm{d}}^{D}k{\mathrm{d}}^{D}q{\mathrm{e}}^{i\mathbf{p}\cdot \mathbf{r}}\frac{\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)}{|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|}.$$

$${\xi }_3^{\ast }(r)\equiv \frac{V^2}{(2\pi {)}^{2D}}\iint {\mathrm{d}}^{D}k{\mathrm{d}}^{D}q{\omega }_D\left(pr\right)\frac{\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)}{|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|},$$

$${\omega }_D(x)=\{\begin{array}{rlr}& \frac{\sin \left(x\right)}{x}& \mathrm{if}D=3,\\ & J_0(x)& \mathrm{if}D=2,\end{array}$$

$${\xi }_3^{\ast }(r)=\sum _{\mathbf{k}}\sum _{\mathbf{q}}{\omega }_D\left(pr\right)\frac{\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)}{|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|},$$

$$s(r)={\left(\frac{r}{L}\right)}^{3D/2}\sum _{k,q\le \pi /r}{\omega }_D\left(pr\right)\frac{\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)}{|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|}.$$

$$r_{\mathrm{peak}}=\left(0.838\pm 0.015\right)R+\left(-0.028\pm 0.099\right).$$

$$\ell (r)={\left(\frac{r}{L}\right)}^{3D/2}\sum _{k,q\le \pi /r}{\omega }_D\left(|\mathbf{k}-\mathbf{q}|r\right)\frac{\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)}{|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|},$$

$$A(l,m)I_{\nu }(l,m)=\sqrt{1-l^2-m^2}\iint V(u,v){\mathrm{e}}^{2i\pi (ul+vm)}\mathrm{d}u\mathrm{d}v,$$

$$V(u,v)={\hat{A}}^{\star }(u,v)\ast {\hat{I}}_{\nu }(u,v)$$

$$V_{ij}^{\mathrm{m}}$$

$$=|G_i{G}_j|{\mathrm{e}}^{i({\varphi }_i-{\varphi }_j)}{V}_{ij}^{\mathrm{true}},$$

$${\mathcal{B}}_{ijk}$$

$$=|G_i{G}_j{G}_k{|}^2{\mathrm{e}}^{i({\varphi }_i-{\varphi }_j)}{\mathrm{e}}^{i({\varphi }_j-{\varphi }_k)}{\mathrm{e}}^{i({\varphi }_k-{\varphi }_i)}{V}_{ij}^{\mathrm{true}}{V}_{jk}^{\mathrm{true}}{V}_{ki}^{\mathrm{true}}$$

$$=|G_i{G}_j{G}_k{|}^2{V}_{ij}^{\mathrm{true}}{V}_{jk}^{\mathrm{true}}{V}_{ki}^{\mathrm{true}},$$

$$\frac{\left(\genfrac{}{}{0ex}{}{N-1}{2}\right)}{\left(\genfrac{}{}{0ex}{}{N}{2}\right)}=\frac{(N-1)(N-2)}{2}\times \frac{2}{N(N-1)}=\frac{N-2}{N}=1-\frac{2}{N}.$$

$${\mathcal{B}}_{ijk}=\left(\sum _m{V}_{ij,{\tau }_m}\right)\times \left(\sum _m{V}_{jk,{\tau }_m}\right)\times \left(\sum _m{V}_{ki,{\tau }_m}\right),$$

$${\theta }_{\mathrm{AR}}=\frac{\lambda }{b_{\max }},$$

$${\Xi }_3^{(1)}(r)={\left(\frac{r}{L}\right)}^{3D/2}\sum _{k,q\le \pi /r}{\omega }_D\left(pr\right)\mathcal{B}\left(\mathbf{k},\mathbf{q}\right).$$

$${\Xi }_3^{(2)}(r)={\left(\frac{r}{L}\right)}^{3D/2}\sum _{k,q\le \pi /r}{\omega }_D\left(pr\right)|\mathcal{B}\left(\mathbf{k},\mathbf{q}\right)|.$$

$$x_{\text{Hii}}(\mathbf{r})=\sum _{i=1}^n\Theta \left(\frac{|\mathbf{r}-{\mathbf{a}}_i|}{R}\right),$$

Full text

††pubyear: 2019††pagerange: Studying the morphology of reionisation with the triangle correlation function of phases–13

1 Introduction

During the Epoch of Reionisation (EoR), from $z\sim 20$ to $z\sim 6$, early light sources ionised the hydrogen and helium atoms of the Intergalactic Medium (IGM). The information currently available on this period comes from indirect observations such as the redshift evolution of the density of star-forming galaxies, thought to be a major source of reionisation and estimated through UV luminosity densities (Bouwens et al., 2015; Ishigaki et al., 2015; Robertson et al., 2015); the integrated Thomson optical depth from Cosmic Microwave Background (CMB) observations (Planck Collaboration et al., 2016b, 2018); and from estimates of the averaged neutral fraction of the IGM obtained via the damping wings of quasars (Greig et al., 2017; Greig et al., 2018), surveys of Lyman-$\alpha$ emitters (Konno et al., 2014; Schenker et al., 2014; Mason et al., 2018) or gamma-ray bursts afterglows (Totani et al., 2014). Although these observations keep improving, in terms of both redshift and precision, they are still not sufficient to draw a precise history of reionisation. Many uncertainties remain on various aspects of the reionisation process such as the emissivity of early galaxies or the level of clumpiness in the IGM (Bouwens et al., 2017; Gorce et al., 2018). Interferometric measurements of the 21cm signal will potentially allow for maps of H i regions in the sky to give us a sense of both the topology and morphology of the reionisation process. To optimise the signal-to-noise ratio of current radio interferometers such as the Murchison Widefield Array (MWA)111http://www.mwatelescope.org, the Precision Array to Probe the Epoch of Reionization (PAPER)222http://eor.berkeley.edu and the Low Frequency Array (LOFAR)333http://www.lofar.org, many works focus on statistical estimators such as the $n$-point correlation functions ($n$-PCF) and in particular the power spectrum ($n=2$). To make the most out of upcoming observations, especially of 21cm tomographic imaging – a key science goal of the future Square Kilometer Array (SKA)444http://www.skatelescope.org/, it is useful to look at higher-order statistics. Previous works have focused on the use of the bispectrum, i.e. the Fourier transform of the 3-PCF to learn about the non-Gaussianity of the reionisation signal (Watkinson et al., 2017; Majumdar et al., 2018; Watkinson et al., 2019; Trott et al., 2019). Because they include both amplitude and phase information, the $n$-PCF ($n\geq 3$) will not complement the information carried by the 2-PCF as well as a correlation function solely based on the Fourier phases, such as the one we describe in this paper.

1.1 Phase information

Consider an ionisation field $x(\bm{r})$ with values ranging from 0 (neutral) to 1 (fully ionised). The visibilities observed by radio interferometers give information in Fourier space, so it is useful to expand $x(\bm{r})$ with a Fourier series,

[TABLE]

Each $\hat{x}(\bm{k})$ is a complex number which can be decomposed into an amplitude $|\hat{x}(\bm{k})|$ and a phase term $\Phi_{\bm{k}}$, so that $\hat{x}(\bm{k})=|\hat{x}(\bm{k})|\,\mathrm{e}^{i\,\Phi_{\bm{k}}}$. When the power spectrum $\mathcal{P}(k)=|\hat{x}(\bm{k})|^{2}$ of the field is considered, all the information contained in the phases is lost. If $x(\bm{r})$ is a pure Gaussian random field (GRF), then the phases are uniformly distributed on the interval $[0,2\pi]$ (Watts et al., 2003) and the power spectrum is sufficient to fully describe the field. Conversely, non-Gaussianity is traceable to the phases. In this work, we estimate the amount of information that can be extracted from phases to add to the usual power spectrum studies. Eventually, we want to use phases to learn about the structure of an ionisation field and extract a characteristic length, such as the average radius of ionised bubbles. In this perspective, we develop a new statistical estimator, derived from the 3-point correlation function and solely based on phases, called the triangle correlation function (TCF). It is similar to the line correlation function introduced by Obreschkow et al. (2013) to study elongated structures in dark matter fields. To consider phase information only, we compute statistics from the phase factor $\hat{\epsilon}(\bm{k})$ defined as

[TABLE]

By construction, the phase factor has an amplitude of one, so that its 2-PCF and power spectrum will vanish. The simplest correlation function related to $\hat{\epsilon}(\bm{k})$ is then its 3-PCF.

The top left panel of Fig. 1 shows an ionisation field made of randomly distributed ionised discs on a neutral background. The top right panel corresponds to the inverse Fourier transform of its phase factor $\hat{\epsilon}(\bm{k})$, a tracer of phase information in real space. We see that phases mostly preserve the edges of the ionised regions: they will provide information about the structure of the field. The same exercise is done in the lower panels of Fig. 1 but for an ionisation field extracted from the semi-numerical 21CMFAST simulation (Mesinger & Furlanetto, 2007) at $z=7.8$ and global ionised fraction $x_{\text{H\,{ii}}}=0.54$ (details on the simulation parameters can be found in Section 4.3). We see that the more complex structure of the ionisation field in this simulation is reflected in a more complex phase map. The boundaries between neutral and ionised regions are not as definite as they were for the toy model above. It will therefore be more difficult to extract information about the structure of this field from its phases.

To this day, little work has been done regarding phase information in the EoR, mostly because of the lack of a solid statistical framework to do so. Thyagarajan et al. (2018) have shown that the 21cm signal from EoR can be detected in bispectrum phase spectra, and mention a potential use of their results for H i intensity mapping experiments. Most of the existing work is related to the study of galaxy clustering: although the initial conditions to matter distribution in the Universe are Gaussian, the growth of cosmic structures will lead to some non-linearity and the distributions of both amplitudes and phases will be altered (Watts et al., 2003; Levrier et al., 2006). By characterising this alteration, one hopes to learn about the formation of cosmic structures and possibly about the properties of dark matter (Obreschkow et al. 2013, see also follow-up papers Wolstenhulme et al. 2015; Eggemeier et al. 2015; Eggemeier & Smith 2017). We will discuss the benefits of phases compared to amplitude information in more details in Section 6.

1.2 Current estimators of the size of ionised regions during reionisation

Characterising the morphology and topology of H ii regions at different times through 21cm tomography is essential to the study of reionisation. For now, two types of metrics have been proposed to do so: Minkowski functionals (Gleser et al., 2006; Lee et al., 2008; Bag et al., 2018; Chen et al., 2018) and bubble size distributions (BSD). Using Minkowski functionals on tomographic images, one can learn about the topology of the reionisation process at a given time, for example the level of percolation or the shape of ionised regions. When the ionised regions largely overlap, there will likely form a unique extended percolated zone pierced by neutral tunnels (Bag et al., 2018; Elbers & van de Weygaert, 2018). Bubble size distributions allow for more quantitative results as they give a direct estimate of the size of ionised regions. To derive BSD, three main methods can be found in the literature: the friend-of-friends algorithm (FOF), based on the connections between ionised regions (Iliev et al., 2006); the spherical average method (SPA), which looks for the largest spherical volume which can cover an ionised region and exceed a given ionisation threshold (Zahn et al., 2007); and the mean-free-path method (RMFP), based on Monte-Carlo Markov Chain (MCMC) techniques and implemented in the semi-numerical simulation 21CMFAST (Mesinger & Furlanetto, 2007; Mesinger et al., 2011). The RMFP algorithm proceeds as follows: it picks a random ionised pixel, stores its distance to the closest neutral pixel in a random direction and stack these distances, measured for many pixels, in a histogram. A thorough comparison of these methods is made in Lin et al. (2016) for simulated clean 21cm maps and in Giri et al. (2018a) for maps including observational effects. A new method, which the authors call granulometry, has been recently presented in Kakiichi et al. (2017). The idea is to successively sieve a binary field with a spherical hole of increasing radius $R$. Depending on how many ionised pixels are sieved, one can build the probability distribution of H ii regions sizes. In their work, the authors prove that granulometry should perform well on future SKA observations, as long as the correct observing strategies are chosen. Lin et al. (2016) have also extended to 3D the well-known watershed algorithm: it connects same value pixels to find isodensity lines which they interpret as the edges of characteristic structures. Note that all these methods require real space data whereas the TCF can be used directly in Fourier space.

We first describe in Section 2 the 21cm signal, on which the simulated observations we use are based. In Section 3, we develop the mathematical formalism leading to the definition of the triangle correlation function, and in particular the definitions of $n$-point correlation functions and poly-spectra. In Section 4, we apply the TCF to simulated ionisation fields,from toy models to outputs of the 21CMFAST simulation. We relate our method to observations, and in particular to the closure phase relation, in Section 5. Finally, we discuss various aspects of our results: we first give evidence for the benefits of phase information in Section 6 and then consider computational performance in Section 7. Conclusions can be found in Section 8.

2 The 21cm signal

The neutral hydrogen 21cm line corresponds to the spin-flip transition of an electron between the two hyperfine levels of the ground state of the H i atom. As a tracer of neutral hydrogen, it is naturally a very interesting observable to learn about the EoR. For the last decades, many efforts have been made to design experiments capable of detecting this signal, despite its low amplitude, the presence of strong foregrounds and huge calibration challenges. Seen with respect to a radio background, the CMB, the evolution of the observed differential 21cm brightness temperature writes (Pritchard & Loeb, 2012):

[TABLE]

where $\delta_{b}$ is the baryon overdensity and $T_{\mathrm{S}}$ the spin temperature, which characterises the relative populations of the two spin states. During the late stages of reionisation, it is expected that the spin temperature will dominate the CMB temperature because of its coupling to the kinetic temperature of the gas, and we can make the approximation of the second line. With this approximation, we can map the 21cm signal on the sky and, thanks to Eq. 3, interpret cold spots, i.e. regions where the signal is weaker than average, either as ionised or underdense regions; and hot spots as neutral or overdense regions. Being able to differentiate underdense from ionised regions in brightness temperature measurements is a major challenge of 21cm tomography (Giri et al. (2018a), see also Giri et al. (2018b) for an efficient way to do so). In this work however, we consider segmented data, converted into binary ionisation fields. Because of redshift, a photon emitted with a wavelength of 21cm during the EoR ($6\lesssim z\lesssim 35$) will reach us today with a frequency $40\lesssim\nu\lesssim 200~{}\mathrm{MHz}$ i.e. the frequency range that the new generation of radio interferometers, such as MWA, PAPER, LOFAR and SKA, are expected to probe.

3 The triangle correlation function of phases

For a real ionisation field $x(\bm{r})$ of volume $V$ in dimension $D$, the $n$-point correlation function ($n$-PCF) measures the correlations between $n$ points, described as $n-1$ vectors $\{\bm{r}_{i}\},~{}i=1,\ldots n-1$. It is defined as:

[TABLE]

where $\bm{r}_{n}=0$. If we average this definition over all possible rotations, we obtain the isotropic $n$-PCF $\xi_{n}$. As can be seen in Eq. 4, the $n$-PCFs are convolutions and it will be easier to compute them in Fourier space. The Fourier transform of the $n$-PCF is called the $n^{\mathrm{th}}$ poly-spectrum $P_{n}$ such that

[TABLE]

and

[TABLE]

For $n=2$, $P_{2}$ is called the power spectrum $\mathcal{P}(\bm{k})$ and for $n=3$, we have the bispectrum $\mathcal{B}(\bm{k},\bm{q})$ such that

[TABLE]

where in the first equation we have used the fact that, because $x\left(\bm{r}\right)$ is a real field, its Fourier transform verifies $\hat{x}(-\bm{k})=\hat{x}^{*}(\bm{k})$. Note that we can also write the bispectrum in terms of three $k$-vectors $(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ forming a closed triangle i.e $\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}=0$.

Consider the 3-point correlation function (3-PCF), i.e. the inverse Fourier transform of the bispectrum:

[TABLE]

To study filamentary structures in matter fields, Obreschkow et al. (2013, hereafter O13) consider the 3-PCF in Eq. 8 for $\bm{r}=-\bm{s}$ i.e. three equidistant points forming a straight line. Because here we look for spherical structures, we will consider two vectors $\bm{r}$ and $\bm{s}$ forming an equilateral triangle, as it is the three-point shape closest to a sphere. In this case, $\bm{s}$ is just $\bm{r}$ rotated by an angle $\pi/3$:

[TABLE]

where the 2D case is limited to the first two equations555Note that in the actual computation of the triangle correlation function, we will consider rotations not only around the $z$-axis but also around the other two.. Let

[TABLE]

then we can write $\bm{k}\cdot\bm{r}+\bm{q}\cdot\bm{s}=\bm{p}\cdot\bm{r}$ in Eq. 8. If we choose to consider the phase factor of the bispectrum rather than its full form (see Eq. 2), we find the modified 3-PCF in dimension $D$

[TABLE]

Note that, according to Eggemeier et al. (2015), the phase factor of the bispectrum can be directly related to the normalised bispectrum defined in Watkinson et al. (2017) and later used in Watkinson et al. (2019) and Trott et al. (2019). Assuming ergodicity, a rotational average of the above expression gives the isotropic modified 3-PCF:

[TABLE]

where $\omega_{D}$ is the window function defined by

[TABLE]

for $J_{0}(x)$ the Bessel function of the first kind and order 0. Numerically, we will need to discretise the integral as

[TABLE]

for a periodic box with physical side length $L$, divided into $N^{D}$ cells. Each cell has a width $\Delta x=L/N$. In Fourier space, this box transforms into a box of same dimensions ($N^{D}$) but of side length $2\pi N/L$ and spacing $\Delta k=2\pi/L$. The largest mode $k=2\pi/\Delta x$ has the smallest wavelength. If we assume phases are uncorrelated below a given scale, then the modes whose wavelengths are smaller than this scale will have random phases. When the resolution is improved, i.e. $\Delta x$ reduced, the number of such modes increases and random phase terms are added to Eq. 14 so that the signal eventually diverges. Following O13, we introduce a cut-off $k\leq\pi/r$ on the sums of Eq. 14 to limit this number. Similarly, when increasing the size of the box $L$, we add modes with wavelength larger than the largest correlation scale within the box, and therefore add random phase terms to the sum. Because $\Delta k=2\pi/L$, the number of modes scales as $L^{D}$ and $s(\bm{r})$ diverges as $L^{3D/2}$. As discussed in O13, we introduce the pre-factor $(r/L)^{3D/2}$ to Eq. 14 to remove this divergence. These changes applied to the isotropic modified 3-PCF $\xi^{*}_{3}(r)$ define the triangle correlation function of phases:

[TABLE]

The imaginary part of $s(\bm{r})$ is zero therefore we only consider its real part in the following sections.

All these derivations have been done with $x(\bm{r})$ being the ionisation field $x_{\text{H\,{ii}}}$. If we consider the neutral field $x_{\text{H\,{i}}}$, we have $x_{\text{H\,{ii}}}=1-x_{\text{H\,{i}}}$ and $\hat{x}_{\text{H\,{ii}}}(\bm{k})=-\hat{x}_{\text{H\,{i}}}(\bm{k})$ which we can plug back into the equations above to find $s_{\text{H\,{ii}}}(r)=-s_{\text{H\,{i}}}(r)$. Therefore when applied to a mostly ionised field containing a few remote neutral islands, the correlations dominating the signal will be related to H i regions and the signal will be negative. In particular, because 21cm interferometric images trace the neutral gas in the sky, if we apply our method on 21cm brightness temperature maps, we obtain the exact same signal as for the corresponding ionisation field, but with a reversed sign. This also implies that during the middle stages of the reionisation process, positive and negative (respectively, H ii and H i) correlations overlap, and the TCF flattens down.

4 Application to simulated ionisation fields

We generate boxes of $N^{2}$ cells filled with randomly distributed ionised disks. Although the derivations performed in Section 3 are valid for a two or three-dimensional ionisation field, we now limit our work to 2D boxes, as they are closer to what one would observe with a radiotelescope. These boxes are binary: an ionised region will have pixels of value 1, whereas a neutral zone will be filled by zero pixels. Because the UV photons emitted by early galaxies have a very short mean free path in the surrounding neutral IGM, the boundary between ionised and neutral regions is expected to be sharp, and such a binary model is acceptable (Furlanetto et al., 2004a). Each simulated box has periodic boundary conditions and ionised regions are allowed to overlap666When they overlap, the maximum pixel value is set to 1, since it is a proxy for ionised level.. In a more realistic field, we would expect the ionised bubbles to be clustered around overdense regions rather than randomly distributed but we choose to first ignore this effect, as we will later consider it when applying our method to the 21CMFAST simulation. For each box, we compute the triangle correlation function as defined in Eq. 15 for a range of correlation scales $r$. We compare the resulting plot with the known size of the ionised bubbles filling the box. The programme we developed to compute these correlations is publicly available online777https://github.com/adeliegorce/Triangle_correlations.

4.1 Picking up characteristic scales

Fig. 2 shows the TCF computed for 30 realisations of a 2D box of $512^{2}$ pixels and side length $L=400~{}\mathrm{Mpc}$ filled with 70 binary bubbles of radius $R=10$ (in pixel units) i.e. about $7.8~{}\mathrm{Mpc}$. The solid line is the mean value for all 30 boxes, and the shaded area represents the $95\%$ confidence interval. Although the signal is very similar for all realisations of the box at small correlation scales, there is more variance at larger scales, as $r$ gets closer to $L$. To see if this scattering is a numerical or physical effect, we compute the triangle correlations of 30 boxes of same dimensions, but with random Fourier phases; the result is added in red to the figure. In this case, the mean signal is close to zero but variance can still be seen at large scales. Interestingly, Eggemeier & Smith (2017) note the same scatter in the line correlation function $\ell(\bm{r})$ of density fields. The line correlation function is the 3-PCF of Fourier phases for three equidistant points aligned in real space (Obreschkow et al., 2013). Through its analytical derivation, they trace this scatter back to the Gaussian part of the covariance matrix of $\ell(\bm{r})$. From another point of view, we see that the mean signal flattens out when many realisations of the same box are considered. This suggests that the power seen on large scales corresponds to correlations between three points located in separate bubbles and so traces a form of correlation in the bubble locations. Because our simulation randomly distributes bubbles in the box, it is then only natural that these correlations would average out. The remaining non-zero value on large scales decreases as we increase the number of bubbles in the box and therefore can be assimilated to shot noise. Note that the oscillations on the left side of the peak are due to the sharp edges of the bubbles: when symmetrical Gaussian distributions are used instead of binary bubbles, or once the field is smoothed by a given angular resolution, they vanish. We refer the interested reader to Watkinson et al. (2017) and Majumdar et al. (2018), where the authors discuss the various problems induced by solid spheres in the bispectrum. Finally, a comparison of numerical results with an analytical derivation of the TCF for a similar toy model can eb found in Appendix A.

The triangle correlation function probes equilateral triangles inscribed in the ionised discs of the box, and simple geometry shows that these have a side length of $\sqrt{3}R$, a scale at which we would expect $s(r)$ to peak. On Fig. 2, there is a clear peak in the signal, but at $r\simeq 6.7~{}\mathrm{Mpc}<\sqrt{3}R$. To check the physical meaning of this peak, we generate many boxes of same dimensions and fill each one of them with 50 binary bubbles of a given radius. Results are gathered in the upper panel of Fig. 3: the TCF peak shifts to larger scales as bubbles increase in size. We plot the peaking scale $r_{\mathrm{peak}}$ as a function of the bubble radius $R$ in the lower panel of Fig. 3, where the error bars correspond to the $95\%$ confidence interval for 10 realisations of the same box. Points are closely aligned, and a MCMC fit to a linear relation gives the following results, where both $r_{\mathrm{peak}}$ and $R$ are in $\mathrm{Mpc}$ and the uncertainty corresponds to the standard deviation on sampled parameters:

[TABLE]

We trace this shift between the peaking scale and the expected peak ($\sqrt{3}R$) back to the window function: when we take $\omega_{D}(x)=1$ in Eq. 15 to compute the TCF, the signal has a shape similar to what is found for the correct window function, but is stretched over the $x$-axis, so that it peaks at a scale $R\leq r_{\mathrm{peak}}\leq\sqrt{3}R$.

Here, we chose to work with a binary field, with no partially ionised regions where $0<x_{\mathrm{H_{II}}}<1$. To see if our results hold when such regions exist, as there would be if X-ray photons significantly contributed to reionisation (Visbal & Loeb, 2012) or if some observed regions are unresolved, we perform the same test but for symmetric 2D Gaussian distributions centred on the bubble radius $R$. The box dimensions are the same as before. As for binary bubbles, we observe a peak in the triangle correlation function at scales slightly smaller than $\sqrt{3}R$. A linear relation can still well describe the correlation between $r_{\mathrm{peak}}$ and $R$ but there is more scatter around the model. This is mostly due to the peaks being difficult to locate: on average, the signal is much flatter for Gaussian than for binary bubbles.

Note that we checked that these results are independent of the number of bubbles in the box. Similarly to what O13 find, the number of bubbles (or filling fraction) will only impact the amplitude of the signal, not its shape. Indeed, for a given radius, the amplitude of the signal will decrease as the number of bubbles – and so the filling fraction of the box, increases. This is likely due to the fact that increasing the number of bubbles increases the number of different translations between objects within the box and so randomises the phase terms (Eggemeier et al., 2015). In general, we find that as long as the filling fraction does not exceed $60\%$, there is still a clear peak in the signal and we can infer a characteristic scale. This corresponds to the fact that in the late stages of the reionisation process, the ionised regions have largely overlapped and they have no characteristic shape and size anymore. We can compare these results to those of Bharadwaj & Pandey (2005) who, through a more observational approach, relate the bispectrum of H i fluctuations to the correlations between the visibilities measured at three different baselines of an interferometer. They generate the same kind of toy models as us and find that their signal increases with the size of the ionised regions, but overlaps prevent their method from being used in the mid-stages of reionisation ($x_{\text{H\,{ii}}}\geq 0.5$).

4.2 Further tests

The TCF seems to correctly pick up the characteristic scale of spherical regions in an ionisation field. Because we want to use this method directly on Fourier data, for which no real space image will be available, we need to ensure that the peaking scale corresponds to spherical regions only. In their work, O13 define the line correlation function (LCF), which is analogous to our triangle correlation function but for two vectors $\bm{r}$ and $\bm{s}$ aligned ($\bm{s}=-\bm{r}$). Therefore $\ell(\bm{r})$ probes elongated rather than spherical structures and writes:

[TABLE]

where the window function has not changed compared to our expression but is now a function of $|\bm{k}-\bm{q}|r$. O13 find the LCF to be close to flat when there are only spherical structures in the field considered (see their Figs. 3 and 4). To check the significance of each correlation function, we compute both of them for fields filled with bubbles of various sizes and compare the results. We find that the LCF also picks up the bubble size, with a peak located at roughly the same correlation scales as for triangle correlations, but with a much weaker amplitude. To push the comparison further, we compute the TCF and the LCF for 2D boxes filled with more or less stretched binary ellipses. We define the prolateness $q$ of an ellipse as the ratio of its semi-minor axis to its semi-major axis so that an ellipse with prolateness $q=1$ is a flat disk. Fig. 4 shows the evolution of the integral of the two correlation functions as a function of the prolateness of the objects in the box 888Note that we integrate only on the correlation range $0.5~{}\mathrm{Mpc}\leq r\leq 20~{}\mathrm{Mpc}$ where the signal is more reliable (see Sec. 4.1).. As $q\rightarrow 1$, the triangle correlations signal increases, while the line correlations signal decreases, even reaching negative values for $q>0.4$. On the contrary, for $q=0.1$, the TCF is close to zero and the LCF shows very strong signal. This confirms that the TCF mostly picks up spherical structures and therefore can safely be used to infer characteristic bubble sizes directly from Fourier data.

In reality, the reionisation process will not be as homogeneous as our toy models: it is unlikely that all ionised bubbles should have the same radius at a given redshift. It would therefore be useful if the triangle correlation function could differentiate between bubbles of different radii. To test for this, we generate a box of $512^{2}$ pixels and side length $L=400~{}\mathrm{Mpc}$, filled with 20 binary disks of radius $14~{}\mathrm{px}$ and 60 of radius $4~{}\mathrm{px}$. Results are displayed on Fig. 5a for 20 realisations of this box. We can clearly distinguish two peaks, which seem to match the bubble radii (dash-dotted vertical lines), once scaled by the coefficients of Eq. 16 (dotted vertical lines). Note that this result holds when we replace the binary bubbles by Gaussian disks. To ensure that the scales picked up by our function correspond to two different sizes of bubbles and not to the length and width of elongated structures, we compute the triangle correlations of a box filled with 60 binary ellipses of the same dimensions as the bubbles above i.e. a semi-major axis $\alpha=14~{}\mathrm{px}$ and a semi-minor axis $\beta=4~{}\mathrm{px}$. The resulting TCF, shown as an orange dashed line on the figure, differs largely from the TCF of discs. There is only one clear peak, corresponding to scales close to the semi-minor axis: because it probes equilateral triangles, the function seems to pick up elongated ellipses as rows of discs of radius $\beta$. The possibility to discern two peaks in the TCF in the presence of different bubble sizes will be limited by two factors. First, the number of bubbles, as we have seen that too many or too few objects leads to a weak signal. We find that when the ratio of large over small bubble numbers exceeds 10 or drops below 0.1, one of the two peaks flattens out. Second, the separation between the two radii: the TCF peaks of bubbles with $R=6~{}\mathrm{px}$ and $R=10~{}\mathrm{px}$ will overlap into one single wide peak and we will not be able to differentiate them anymore.

Analytic derivations predict a log-normal distribution for the bubble sizes during reionisation, which becomes increasingly peaked as bubble grow and merge (Furlanetto et al., 2004b). This first result has been confirmed by many authors looking at simulations (McQuinn et al., 2007; Zahn et al., 2007; Lin et al., 2016). Let two boxes where the bubble centres are randomly distributed, but the radius sizes are sampled from a log-normal distribution. Because we generate the box, we know the exact distribution of bubble radii within and we want the TCF to trace it. On Fig. 5b, we compare the TCF (solid line) with the actual radius distribution (histogram) for 2 different bubble size distributions: on centred on $R=5~{}\mathrm{px}$ with variance $2~{}\mathrm{px}$, corresponding to the early stages of EoR, and one with $R=11~{}\mathrm{px}$ with variance $1~{}\mathrm{px}$, corresponding to later stages. The number of bubbles is adjusted to reach a filling fraction $x_{\text{H\,{ii}}}=0.09$. Overall, there is a reasonably good match between the two. We see that as the mean radius increases, the discrepancy between peaking scale and real bubble size grows in a way that the linear relation of Eq. 16 cannot fully compensate for. Note that results are similar for a Gaussian distribution of the bubble sizes.

4.3 Results on 21CMFAST simulation

Now that our method has been tested and characterised on toy models, we apply it to more physical ionisation fields, extracted from the semi-numerical simulation of reionisation 21CMFAST999http://github.com/andreimesinger/21cmFAST (Mesinger & Furlanetto, 2007; Mesinger et al., 2011). We choose 21CMFAST because, thanks to relatively short executing times, it gives a large flexibility in parametrisation. It also includes an implementation of the random mean free path algorithm to estimate the bubble size distribution of data cubes, which we can compare to our estimator. Note that an important difference with the toy models is that the bubble locations are no longer random, and the TCF might pick up correlations between the bubble centres on large scales. The 21CMFAST code first generates a high-redshift linear density field which is then evolved to lower redshifts using linear theory and the Zel’dovich approximation. The ionisation field is extracted from this density field using excursion-set theory for haloes with virial temperature $T>10^{4}~{}\mathrm{K}$. We use this code to generate a box with sufficient resolution to obtain a H i field of $512^{3}$ pixels and side length $L=400~{}\mathrm{Mpc}$ for the following cosmology: $\Omega_{m}=0.308$, $\Omega_{b}=0.049$ and $H_{0}=67.74~{}\mathrm{km.s}^{-1}\mathrm{Mpc}^{-1}$ (Planck Collaboration et al., 2016a). The resulting reionisation history has its midpoint at $z=7.9$ for a duration of $\Delta z=z\left(x_{\text{H\,{ii}}}=0.10\right)-z\left(x_{\text{H\,{ii}}}=0.99\right)=4.9$ and gives an integrated Thomson optical depth of $\tau=0.067$ (see Fig. 7).

The output of the simulation is a 3D H i field but we choose to analyse 2D slices to be closer to actual observations. We also convert the given H i field into a binary H ii field in order to get positive triangle correlations when most of the field is ionised and negative correlations when the sky is mainly neutral (see Eq. 7b for $x_{\text{H\,{i}}}=1-x_{\text{H\,{ii}}}$ and discussion at the end of Sec. 3), in continuity with previous sections. Fig. 7 presents the TCF of 2D slices at redshifts $z=12.0,11.0,10.0,9.0,8.0,7.0,6.4$ and $6.2$ corresponding to global ionisation levels of $x_{\text{H\,{ii}}}=0.09,0.10,0.18,0.29,0.49,0.80,0.97$ and $0.99$ respectively, along with a picture of the corresponding H ii field. For each simulation, results are shown with variance estimated from the TCF of 20 Gaussian random fields (GRF) of same dimensions. This indicates that the signal is more significant at small scales, a point already mentioned in Section 4.1 and gives an idea of how well this field differs from a GRF. However, these errors will not take into account the non-Gaussian nature of the field if the 2D slice considered is not representative of the overall field, especially at high redshifts where the signal is highly non Gaussian. Indeed, consider the $z=6.2$ slice on Fig. 7: the small ionised regions seen in the field are isolated features. If we compute triangle correlations for more realisations of this simulation, it is unlikely that the exact same feature will appear in the field and so in the TCF. Therefore there will be more variance on small scales than the error bars on the Figure. Ideally, one would find a theoretical expression for the covariance of the TCF, but we keep this derivation for future work.

At high redshift, in what is often referred to as the pre-overlap phase, we see on the TCFs of Fig. 7 that there is more power at small scales but no clear peak. This is likely due to the variety of H ii regions: many small ionised regions around young sources parcel the neutral background out. However, by looking at the signal-to-noise ratio (SNR), we can infer an upper limit on the sizes of the ionised regions: at $z=13$, scales smaller than $8.9~{}\mathrm{Mpc}$ (after conversion with Eq. 16) contribute for $80\%$ of the cumulative SNR; at $z=12$ and $z=11$, this upper limit increases respectively to $10.8~{}\mathrm{Mpc}$ and $11.2~{}\mathrm{Mpc}$. The structure of the ionisation field at high redshift directly relates to the way the 21CMFAST algorithm is constructed and particularly to the use of excursion set theory. In excursion-set theory, when the average density of a region exceeds a given threshold, it collapses to form a halo. When applied to reionisation, the threshold additionally considers ionising photons production: if the region has produced a sufficient number of ionising photons with respect to its mass and volume, then it is considered ionised. Because of this, the ionisation field in the early stages of reionisation will have many very small sources rather than a few efficient sources. It would be interesting to compute the TCF of other reionisation simulations, in particular simulations with more efficient escape fractions (Seiler et al., 2019). This also includes simulations where reionisation is led by Active Galactic Nuclei (AGN): ionising sources are more scarce and have a better ionising efficiency, leading to a topology more similar to the toy models used in Sec. 4.1. In this perspective, our method may be able to differentiate between different reionisation scenarios. Such a study goes beyond the scope of this work but we refer the interested reader to Watkinson et al. (2019), where the authors use the bispectrum as a probe for non-Gaussianities due to X-ray heating.

Later on, when the global ionisation fraction reaches values between $25\%$ and $75\%$, negative signal coming from neutral regions overlaps with the positive signal from ionised regions. The TCF flattens, and therefore cannot give information about the morphology of the field. However, measuring a flat signal from actual data could be interpreted as the reionisation process being in its middle stages. For $z>7$, most of the sky is ionised and only a few remote neutral islands remain. This is the post-overlap phase. We see on Fig. 7 that the sizes of these neutral islands are efficiently picked up by the TCF. For $z=6.2$, there is a very clear negative peak at scales $r=7.7~{}\mathrm{Mpc}$ which correspond to a radius size of $r^{\prime}=9.2~{}\mathrm{Mpc}$ once the linear relation in Eq. 16 is applied as a correction. If we roughly estimate the size of the neutral zones, we find that they are about $\sim 11~{}\mathrm{px}\sim 8.6~{}\mathrm{Mpc}$ in radius, which is very close to the estimation given by triangle correlations. For $z=6.4$, we see two clear peaks in the signal, corresponding to the two sizes of ionised regions seen in the corresponding real space field. The first peak is spread over scales $4.6<r<7.8~{}\mathrm{Mpc}$ i.e. $5.5<r^{\prime}<9.3~{}\mathrm{Mpc}$, while the second one is more narrow, centred around $20.2~{}\mathrm{Mpc}$. For comparison, we find that we can sieve the two larger neutral islands in the field with disks of radius $\sim 23~{}\mathrm{px}\sim 18~{}\mathrm{Mpc}$; whereas the smaller ones can fit in disks whose radii range from $7~{}\mathrm{Mpc}$ to $12~{}\mathrm{Mpc}$. This motivates a further study of how well the TCF performs compared to common BSD algorithms.

Let’s consider the spherical average (SPA) and the random mean free path (RMFP) methods. Note that these two algorithms need to be applied to real space images whereas the TCF is used on Fourier data. We choose not to compare our results to the friend-of-friends algorithm (Iliev et al., 2006) because it considers overlapping bubbles as a unique large ionised region, which is fundamentally different from our approach. We use versions of SPA and RMFP implemented by the authors of Giri et al. (2018a)101010tools21cm are found on https://github.com/sambit-giri/tools21cm.. Recall that the SPA algorithm looks for the largest sphere around each ionised region whose ionisation level exceeds a given threshold, here chosen to be $x_{\text{H\,{ii}}}=0.5$. RMFP, on the other side, relies on MCMC: it looks for the first neutral cell encountered in a random direction starting from an ionised pixel and records the length of the ray to build a histogram. The SPA and RMFP outputs are probability distributions of radii.Results can be seen on Fig. 8 for $z=6.2,6.4,11.0$ and 12.0 and are detailed in Table 1. For clarity, we plot $-s(r)$ in the first two panels to have a positive signal for all redshifts. At low redshift, our method mostly agree with the other BSDs, although the TCF exhibits a narrower peak; we are also able to clearly distinguish between two characteristic scales for $z=6.4$, whereas SPA et RMFP only give one most likely scale, apparently corresponding to the smaller neutral regions. At higher redshift, the SPA distributions tend to $r\sim 0~{}\mathrm{Mpc}$ and we choose to compare values of quantiles rather than maximum likelihood values: at $z=11$, according to SPA, there is a $80\%$ probability that the ionised regions have a radius smaller than $13.3~{}\mathrm{Mpc}$; similarly, scales $<11.2~{}\mathrm{Mpc}$ represent $80\%$ of the cumulative SNR for the TCF. These upper limits, computed for both $z=11$ and $z=12$, are shown as vertical lines on Fig. 8. RMFP gives a maximum likelihood radius of $3.12\pm 0.78~{}\mathrm{Mpc}$ (the error corresponding to the bin size), therefore the most quantitative result. For all three methods, the upper limit increases as redshift decreases, which corresponds to ionised regions growing with time. Note that the upper limit given by SPA is the lowest: our results corroborate those of Lin et al. (2016), who argue that SPA is heavily biased towards smaller bubble radii than the actual bubble sizes. Indeed, they find that for a single bubble of radius $R$, the SPA probability distribution peaks at $R/3$. On the contrary, they find the RMFP method to be unbiased and to peak at the correct bubble size. Note that on the two right panels of Fig. 8, we see some signal at large scales. This can be related to statistical noise on one side – we analyse a unique slice of the simulation; and to the potential correlations between separate bubbles, which are not randomly distributed anymore.

Overall, our estimator performs well with respect to comparable BSD algorithms. However, for this method to be an useful tool in the analysis of upcoming radio observations, a deeper analysis, including the study of different types of reionisation simulations, would be required. In practise, we would also use this method as a forward modelling process, and compare triangle correlations from observations to a set of predicted signals for different types of ionisation fields. We keep the construction of such a database for future work.

5 Relation to observations

5.1 Visibilities and closure phases

When gathering data with an interferometer, the signal will be measured for pairs of antennae separated by a baseline $D$. The measurement is made in terms of a complex visibility $V(u,v)$, where $u$ and $v$ are the projection of the baseline in wavelength units on the plane perpendicular to the vector $\bm{s}_{0}$ pointing to the phase reference position, i.e. the centre of the field to be imaged. In Eq. 3, the brightness temperature depends on the position in the sky $\bm{r}$ where the signal is observed: $\delta T_{\mathrm{b}}(\bm{r},z)\propto\,x_{\text{H\,{i}}}\,(\bm{r},z)$. $\delta T_{\mathrm{b}}$ is the intensity $I_{\nu}$ of the 21cm signal observed through the interferometer for a given frequency $\nu$ (i.e. for a given redshift, since $\nu=\nu_{21\mathrm{cm}}/(1+z)$)111111We assume the instrument to have narrow bandpass filters and so to probe exactly $\nu$ rather than a frequency band centred on $\nu$.. In what follows, we will use $I_{\nu}(\bm{r})=\bar{I}_{\nu}\,x_{\text{H\,{i}}}(\bm{r},z)$ and $\delta T_{\mathrm{b}}(\bm{r},z)$ interchangebly. If the interferometer probes a sufficiently small region of the sky compared to the beam width of the antennae, we can approximate this region by a flat plane. Then the source intensity distribution is a function of two real spatial variables $(l,m)$ and the van-Cittert theorem tells us that the complex visibility is the 2D inverse Fourier transform of the intensity, corrected by the normalised average effective collecting area $A(l,m)$ of the two antennae (Thompson et al., 2017):

[TABLE]

that we can re-write as

[TABLE]

where $*$ denotes a convolution, $\hat{A}$ and $\hat{I}_{\nu}$ are respectively the 2D Fourier transforms of $A$ and $I_{\nu}$ and we define $A^{\star}=A/\sqrt{1-l^{2}-m^{2}}$. Note that this approximation introduces a phase error that may need to be taken into account when applying our method, based on phase information. Because $A(l,m)$ is a measurable instrumental characteristic, if we know our instrument sufficiently well, we can deconvolve the measured complex visibilities by $\hat{A}^{\star}$ and obtain, after correcting for the different pre-factors, the $\hat{x}(\bm{k})$ terms to compute the triangle correlation function.

One of the interests methods based on phase information is that the phases of the measured visibilities, combined in a bispectrum, will not be sensitive to errors related to calibration (Jennison, 1958; Monnier, 2007). Consider the visibility $V_{ij}^{\mathrm{m}}$ measured between two antennae $i$ and $j$ at a given frequency $\nu$. It will have contributions from the true visibility $V_{ij}^{\mathrm{true}}$, coming from the cosmological signal, but also from the amplitude and phase errors of each antenna, modelled by a complex gain $G_{i}=|G_{i}|\,\mathrm{e}^{i\phi_{i}}$, such that:

[TABLE]

where $*$ denotes the complex conjugate. The amplitude of this gain will come from beam specific effects such as mirror reflectivity, detector sensitivity or local scintillations whereas the phase term can originate either from telescope errors or from outside effects such as atmospheric turbulence (Levrier et al., 2006; Monnier, 2007). If we combine the signal from three antennae forming a closed triangle, we can avoid this phase error and we will be left only with what is called the closure phase. Indeed, consider three baselines $ij$, $jk$ and $ki$ observing at the same given frequency $\nu$. We can write the bispectrum of their complex visibilities as

[TABLE]

where in the last step, the different phase terms cancel each other out: the phase of the measured bispectrum is the phase of the true bispectrum. By construction, the bispectrum $\mathcal{B}(\bm{k},\bm{q})$ in Eq. 7b considers a closed triangle configuration, for three vectors $\bm{k}$, $\bm{q}$ and $-\bm{k}-\bm{q}$. Because we choose to work in two dimensions, the three vectors lie in the same plane on the sky, perpendicular to the line-of-sight, and are measured for the same frequency so that the closure relation holds. We will then be able to use our method on observational data without worrying about calibration errors. For an example of the use of bispectrum closure phases in interferometry, we refer the reader to Thyagarajan et al. (2018). In this work, the authors compare the bispectrum phase spectra, i.e. the phases of the bispectrum $\mathcal{B}_{ijk}$ in Eq. 21, coming from different components of a simulated signal, i.e. a single point source, diffuse foregrounds and H i fluctuations from the EoR and demonstrate that a quantitative relationship exists between the EoR signal strength and the whole bispectrum phase power spectra.

These properties are a major benefit of our technique compared to other solutions found in the literature to estimate the characteristic size of ionised regions from interferometric data as these always require to reconstruct the real-space image corresponding to observations. Indeed, our method can directly use as input the complex visibilities observed by an interferometer, and if we choose to use closure phases, results will be independent of antenna-based calibration and calibration errors. However, in practise, there are some limitations to the use of the closure relation.

First, only a limited number of triangles can be constructed from the array of antennae of a telescope, therefore some information will be lost compared to simple baseline measurements. Monnier (2007) count that for an array made of $N$ antennae, there are $n=\binom{N}{3}=N(N-1)(N-2)/6$ possible closed triangles, $\binom{N}{2}$ independent Fourier phases and $\binom{N-1}{2}$ independent closure phases. Therefore the amount of phase information recovered from closure phases is

[TABLE]

With as little as 40 antennae, we are able to recover $95\%$ of the phase information but 1000 antennae are not enough to reach $99.9\%$. Note that for the 296 antennae of the SKA1-Low central array and the 48 antennae of LOFAR, we recover respectively $99.3\%$ and $95.8\%$ of the phase information from closure phases. However, even if almost all the phase information is recovered, Readhead et al. (1988) find that the noise level of a phase-only observation will still be at least twice higher than a map made from full visibility data, because when ignoring amplitudes, half of the signal is lost. Additionally, Readhead et al. (1988) show that because the bispectrum is a triple product, there will be new sources of noise compared to single baseline observations. First, if the same signal is measured on different time intervals, the observed bispectrum will not only be the product of three complex numbers $V_{ij}$ anymore, but of the sum of each observation on each time interval:

[TABLE]

where $V_{ij,\tau_{m}}$ is the visibility measured for baseline $ij$ on time interval $\tau_{m}$. Then the cross terms combining signals integrated on different time intervals (e.g. $V_{ij,\tau_{m}}$ and $V_{jk,\tau_{l}}$) will give incoherent phase terms that can be assimilated to noise. The second potential source of noise mentioned by Readhead et al. (1988) corresponds to the same kind of reasoning, but in the spatial domain: if we consider the triangle formed by three baselines $(ij,ik,kl)$, then on a redundant array there will be contributions not only from the $ij+jk+kl=0$ triangles but also from identical baselines who are not part of a triangle. See the example on Fig. 9 of a redundant array: all numbered vectors correspond to the same baseline but only 1, 2 and 4 form triangles. When probing $ijk$ triangles, the measured bispectrum will take the sum of the four signals as the visibility corresponding to this baseline (e.g. $V_{ij}=V_{1}+V_{2}+V_{3}+V_{4}$) and additional cross terms, irrelevant to closure phases because they include $V_{3}$, will arise.

Readhead et al. (1988) show however than these two types of noise can be reduced to a sensible signal-to-noise ratio if enough frames are used in the integration. Similarly, Carilli et al. (2018) mention that some instrumental effects such as polarisation leakage or cross coupling of antennae, called "closure errors" can lead to a departure from the closure relation.

Finally, the limited number of triangles one can construct from $N$ given antennae will limit the sampling of $(u,v)$ space and worsen the sparsity of observations. Applying our method to sparse data goes beyond the scope of this work, but it would be interesting to see how triangle correlations perform with this additional difficulty. For now, we will limit ourselves to noisy data sets. We refer the interested reader to Trott et al. (2019), where the authors compute the normalised bispectrum defined in Watkinson et al. (2019) from the closed triangle visibilities of MWA Phase II.

5.2 Instrumental effects

To see how our method performs when applied to actual observations, we now add instrumental effects to our ionisation maps: beam smoothing and noise corresponding to observations by the core of SKA1-Low, its central area, and by LOFAR. We pick three simulated ionisation fields: the first comes from the toy models described in Section 4. It is made of 70 bubbles of radius $R=10~{}\mathrm{px}$, and assumed to correspond to a redshift $z=9$. The two others are extracted from 21CMFAST at redshifts $z=6.2$ and 6.4. They are shown, in this order, on the left panels of Fig. 10. We first convert each one of them into a brightness temperature map according to Eq. 3 with the following cosmology: $\Omega_{m}=0.309$, $\Omega_{b}=0.049$ and $H_{0}=67.74~{}\mathrm{km.s}^{-1}\mathrm{Mpc}^{-1}$.

The resolution of 21-cm tomographic data will be first limited by the angular resolution of the interferometer considered. The full width at half maximum (FWHM) of an interferometer is given by (in radians):

[TABLE]

where $\lambda$ is the redshifted 21cm wavelength i.e. $\lambda=21\mathrm{cm}\times(1+z)$ and $b_{\mathrm{max}}$ is the maximum baseline of the interferometer. We have $b_{\mathrm{max}}=3500$, $1000$ and $3400~{}\mathrm{m}$ for LOFAR, SKA1-Low core and SKA1-Low central respectively. To account for this effect, we convolve the $\delta T_{b}$ map with a Gaussian kernel of FWHM $\theta_{\mathrm{AR}}d_{\mathrm{c}}(z)$, with $d_{\mathrm{c}}(z)$ the comoving distance at the redshift (i.e. frequency) considered. For SKA central, we get the second-to-left panels on the figure – note that because they have similar maximum baselines, the angular resolution of SKA central is close the one of LOFAR. For SKA core, because $b_{\mathrm{max}}$ is much smaller, the smoothing blurs the shape of the ionised bubbles to an extreme point, and our method will perform poorly. Finally, we simulate realistic instrumental noise by using the measurement equation software OSKAR121212https://github.com/OxfordSKA/OSKAR.

Fig. 10 presents the results: the first column corresponds to triangle correlations for a clean field. We identify the peaking scale for this clean field in order to compare it to the scale picked on corrupted data later on and indicate it as a dotted blue line on each plot. Note that for the third field (21CMFAST at $z=6.4$), there are two scales picked up. The larger one is well defined, whereas we use an interval $6.2\pm 1.6~{}\mathrm{Mpc}$ for the smaller one. On the figure, the second column is for the field smoothed with the angular resolution of SKA1-Low central; and the last three for the smoothed signal with instrumental noise from each of the three experiments considered (respectively SKA1-Low core, central and LOFAR) and 1000 hours of integration time. We choose 1000h as it is the most commonly used value in the literature, but recent works have shown that as few as 324h of observations with SKA can be sufficient to differentiate between different reionisation models (Binnie & Pritchard, 2019). The dash-dotted lines on each plot mark the smoothing scale of the corresponding experiment, whose values are given in Table 2. One can see that results with SKA core are not satisfying because of a too low angular resolution that blurs the edges of the ionised regions: the signal is mostly flat or seems to pick up the smoothing scale. However, results for SKA central and LOFAR prove very satisfactory on the toy model (first row), with a clear peak at a scale close to the one found with clean signal (see Table 2 for details): applying the linear relation of Eq. 16 to the peaking scale could give a good approximation of the actual size of ionised regions in the image. When the integration time is increased, we get for both LOFAR and SKA1-Low central a peaking scale even closer to the one found for clean signal. Unfortunately, the LOFAR sensitivity does not allow to extract much information from the 21CMFAST maps: at $z=6.2$ the signal is mostly flat, and although there is a clear peak for the $z=6.4$ map, it cannot resolve the two characteristic scales, however picked up when accounting for SKA central sensitivity. Note that for each plot, the error bars correspond to the variance of the TCF computed for a GRF of same dimensions. We have also tried our statistic on 21CMFAST boxes at higher redshifts, and we find that, because there is no clear characteristic scale in the field (see Section 4.3), the TCF peaks at the smoothing scale corresponding to the telescope considered (listed in Table 2). There is however no risk of misinterpretation of this peak since real data would be deconvolved from telescope properties such as angular resolution before being analysed.

In all the work above, we assumed that foreground pollution was completely removed from the data cubes: Chapman et al. (2015) presented efficient foreground removal techniques that will produce good quality 21cm maps. However, foreground residuals could still impact our results. A precise instrumental calibration is usually required to separate the 21cm signal from foregrounds and is one of the main challenges of upcoming experiments although as mentioned before, it seems that bispectrum phases and therefore our results would be unaffected by calibration errors (Thyagarajan et al., 2018).

6 Why phase information?

In order to see if our phase-only estimator truly supplements power spectrum information, we generate a 2D ionisation field made of 70 binary disks of radius $R=10~{}\mathrm{px}=7.8~{}\mathrm{Mpc}$ and compute its 2-point correlation function and its TCF. We then shuffle its Fourier phases – replace them by random phases ranging from 0 to $2\pi$, and compute the corresponding 2-PCF and TCF. Results can be seen on the upper and lower panels of Fig. 11 respectively. On the left panels, we see that reshuffling the phases has made the field lose all its structure: there are no bubbles any more. Because we kept the absolute value of the field unchanged, the 2-PCFs in the second column are exactly identical. However, the extra information carried by phases is clearly visible when comparing the TCFs on the right-hand panels: there is almost no signal for the field with random phases whereas the field with structured phases exhibits a clear peak.

For completeness, we compute three different types of triangle correlations encompassing more or less phase information and compare results on Fig. 12. We do so for three boxes filled with 70 binary bubbles: the first (upper panel) with bubbles of radius $R=7$, the second (middle panel) $R=9$ and the third (lower panel) $R=12$.

The first type of triangle correlations considered is the one we have used so far, defined in Eq. 15 as a modified inverse Fourier transform of the phase factor of the bispectrum (blue solid line on Fig. 12). 2. 2.

The second uses the full bispectrum, amplitude included (purple dashed line on Fig. 12):

[TABLE] 3. 3.

The third only uses the amplitude of the bispectrum (green dotted line on Fig. 12):

[TABLE]

We see that the main advantage of the phase-only TCF is that, compared to statistics using bispectrum amplitude, the peaking scale is much better defined. It is also less sensitive to the filling fraction than the two other functions: in the bottom panel, where bubbles are larger and the filling fraction higher, phase correlations still exhibit a well-defined peak whereas the other two flatten.

7 Computational performance

In Sec. 4.3 we compared the results of different methods to derive a bubble size distribution or equivalent, but we did not mention computational performance. Note that our code has been parallelised via OpenMP (OpenMP Architecture Review Board, 2013) in order to reduce computing times: because most of the code consists of sums, i.e. nested loops, the parallelisation is very efficient. Three parameters will play an essential role in determining the computing time necessary to evaluate triangle correlations:

The number $N$ of pixels in the box. 2. 2.

The range of correlations scales $r$ for which we compute the triangle correlations: because of the limits of our sum ($\sum_{k\leq\pi/r}$, see Eq. 15), there will be more terms to sum over for smaller scales. Therefore, in this work, we have chosen to compute triangle correlations for $0.5~{}\mathrm{Mpc}\leq r\leq 30~{}\mathrm{Mpc}$. 3. 3.

The box side length $L$, because it determines the sampling of $k$-cells: for smaller $\Delta k$, there will be more modes with norm smaller than $\pi/r$ and hence more terms to sum over.

Most of the boxes analysed in this work had $512^{2}$ pixels for a side length $L=400~{}\mathrm{Mpc}$ and triangles correlations were computed for $r/\mathrm{Mpc}\in[0.5,30]$, for an average computation time of 180 minutes on 20 cores. For comparison, the SPA algorithm, as implemented by Giri et al. (2018a), takes about half an hour to run on a 3D box with same side length and sampling; whereas RMFP, as implemented in 21CMFAST, only takes a couple minutes. Therefore more work is required to improve the computational efficiency of our algorithm further and so its usability. An approach similar to the one presented in the Appendix of Eggemeier et al. (2015), where the nested sums in Eq. 15 are replaced by rotational averages in real space, could greatly improve computation times.

8 Conclusions

Following the work of O13, we have constructed a new statistical tool, based on phase information only and called the triangle correlation function (TCF), which can be used to determine the characteristic scale of ionised regions on 21cm interferometric data from the Epoch of Reionisation. Indeed, if we plot this function over a range of correlation scales for simple fields, made of perfectly spherical fully ionised regions on a neutral background, we find a peaked signal, and the peaking scale can be directly related to the actual size of the bubbles. From such toy models, we have derived important properties for the TCF: it can differentiate between ionised regions of different sizes, and distinguish spherical from elongated structures. We have also found its results to be more reliable on scales smaller than about a twentieth of the physical side length of the field studied, as the finite size of the box implies more sample variance at larger scales. To see if our method can be applied to observational data, we have confronted the TCF with ionisation fields corrupted by instrumental noise or angular resolution from LOFAR or SKA, and we have found that it still performs well at giving the characteristic scale of ionised regions, as long as the integration time is sufficient – here, we have worked with 1000 hours. By comparing results for 3-PCF including amplitude information to $s(r)$, we also proved that a statistical tool using phase information only will be more efficient at picking up the characteristic scale of spherical structures in a field than a correlation function also using information from the bispectrum amplitude. Indeed, the phase-only function is more peaked, allowing a better identification of the characteristic scale. Moving on to the more elaborate reionisation simulation 21CMFAST (Mesinger & Furlanetto, 2007), we have found that our method gives a good estimation of the size of remote neutral islands at the very end of the reionisation process. In particular, and contrarily to other BSD algorithms such as RMFP and SPA, it is able to resolve two characteristic scales in a field. In the early stages of the simulation, because there are many very small ionised regions covering the neutral background, the signal is dominated by Poisson noise, and there is no clear characteristic scale to pick up. Therefore, for fields with a very low global ionised fraction, we rather use our method to infer an upper limit on the size of ionised regions. Note that during the overlap phase, positive correlations from ionised regions overlap with negative signal from neutral zones and the signal flattens out: an absence of signal can be interpreted as reionisation being in its middle stages and so we can learn about the duration of the process. In general, it will be more difficult to extract phase information from non binary ionisation field because of the more complex structure of the field, which is reverberated in the phase information: Fig. 1 compares the real-space phase information of one of our toy models with $\epsilon(\bm{r})$ for one of our 21CMFAST boxes. These results, on both toy models and more elaboration simulations, hold when instrumental effects such as telescope angular resolution and instrumental noise are added to clean maps.

Apart from noise and instrumental resolution, there will be some other important questions to solve before applying our method to true 21cm tomographic data. Kakiichi et al. (2017) proved that the cold spots in 21cm tomographic images trace H ii regions more efficiently at low redshift ($z\leq 7$) and high filling fraction ($x_{\text{H\,{ii}}}\geq 0.4$): earlier in the reionisation process, cold spots due to a local underdensity rather than ionisation are more frequent (see Sec. 2). Corroborated with the fact that the TCF performs well on remote neutral islands at the end of the reionisation process, we can expect to get good results on low redshift data. However, correctly identifying ionised regions in 21cm observations i.e. transforming a map of the differential brightness temperature $\delta T_{b}$ into a binary field made of fully ionised and fully neutral regions remains a challenge – for an extended discussion of this issue, see Giri et al. (2018a). Note that this is not an issue for the friends-of-friends (Iliev et al., 2006; Friedrich et al., 2011) and watershed methods (Lin et al., 2016) since they segment the data themselves. Additionally, it remains to be seen exactly how to exploit closure phases, and what is the effect of the quality of the data, especially its sparsity, on our results. We would also need to see how foreground removal techniques impact the non-Gaussianity of the signal, and therefore our results. Finally, the practical application of triangle correlations to observational data would require forward-modelling i.e. comparing measurements to the signal obtained for a number of simulations corresponding to various parametrisations of reionisation to infer the reionisation scenario corresponding to the observation. In this perspective, one would need to compute the triangle correlations of many more simulations.

Acknowledgements

The authors thank the referee for a useful report which helped improve the results of this paper. They also thank Catherine A. Watkinson, Marian Douspis and Benoît Semelin for useful discussions on various aspects of the analysis; Nithyanandan Thyagarajan for his helpful comments on a draft version of this paper; and Emma Chapman for kindly providing results from OSKAR runs to perform our noise analysis. AG and JRP acknowledge financial support from the European Research Council under ERC grant number 638743-FIRSTDAWN. AG’s work is supported by a PhD studentship from the UK Science and Technology Facilities Council (STFC).

This research made use of astropy, a community-developed core Python package for astronomy (Astropy Collaboration et al., 2013, 2018); matplotlib, a Python library for publication quality graphics (Hunter, 2007); scipy, a Python-based ecosystem of open-source software for mathematics, science, and engineering (Jones et al., 2001) – including numpy (Oliphant, 2006), and emcee, an implementation of the affine invariant MCMC ensemble sampler (Foreman-Mackey et al., 2013).

Appendix A Analytical derivation for a toy model

Consider a 2D box of volume $V=L^{2}$ filled with $n$ fully ionised bubbles of radius $R$, randomly distributed throughout the box so that their centres are located at $\bm{a}_{i}$ for $i\in\{1,n\}$. The ionisation field of the box is then

[TABLE]

where $\Theta\left(x\right)$ is the Heaviside step function worth 1 if $x\leq 1$ and 0 otherwise. We need to take the Fourier transform of this field, derived as follows:

[TABLE]

where to go from the first to the second line, we have used the change of variables $\bm{r}=\bm{r}-\bm{a}_{i}$. From the fourth to the fifth, we have let $z=kr$ to be able to use:

[TABLE]

If we define $W(y)=J_{1}(y)/y$, the Fourier transform of the top-hat window function in 2D, we have the final expression of the Fourier transform of our ionisation field:

[TABLE]

Indeed, if we let $\rho$ the mean number density of ionised bubbles such that $\rho=n/V$ and ignore overlapping131313This will be a reasonable assumption for the filling fractions considered., $\bar{x}_{\text{H\,{ii}}}=\pi R^{2}\rho$. From this we deduce analytic expressions for the power spectrum

[TABLE]

which is a real number; and the bispectrum of the box

[TABLE]

To find the analytic expression of our triangle correlations function for this toy model, we then plug Eq. 32 into Eq. 15.

In Fig. 13, we compare the result of a numerical integration of the triangle correlation function computed from the bispectrum in Eq. 7b with a version where we compute analytically the triangle correlation function, knowing the locations of the ionised bubbles and the size of the box only, according to Eq. 32 and detailed above. We consider 20 realisations of a box of dimensions $L=80~{}\mathrm{Mpc}$ for $100^{2}$ pixels, filled with 20 bubbles of radius $R=3~{}\mathrm{px}=2.4~{}\mathrm{Mpc}$141414The choices of box size and bubble number are limited by the computational cost of the analytic method.. We see that there is a good match between both computational methods: we have the confirmation that our estimator indeed traces the bubble distribution in the ionisation field considered.