Gridded and direct Epoch of Reionisation bispectrum estimates using the Murchison Widefield Array

Abstract

We apply two methods to estimate the 21-cm bispectrum from data taken within the Epoch of Reionisation (EoR) project of the Murchison Widefield Array (MWA). Using data acquired with the Phase II compact array allows a direct bispectrum estimate to be undertaken on the multiple redundantly spaced triangles of antenna tiles, as well as an estimate based on data gridded to the uv-plane. The direct and gridded bispectrum estimators are applied to 21 h of high-band (167–197 MHz; z = 6.2–7.5) data from the 2016 and 2017 observing seasons. Analytic predictions for the bispectrum bias and variance for point-source foregrounds are derived. We compare the output of these approaches, the foreground contribution to the signal, and future prospects for measuring the bispectra with redundant and non-redundant arrays. We find that some triangle configurations yield bispectrum estimates that are consistent with the expected noise level after 10 h, while equilateral configurations are strongly foreground-dominated. Careful choice of triangle configurations may be made to reduce foreground bias that hinders power spectrum estimators, and the 21-cm bispectrum may be accessible in less time than the 21-cm power spectrum for some wave modes, with detections in hundreds of hours.

1. Introduction

Exploration of the growth of structure in the first billion years of the Universe is a key observational driver for many experiments. One tracer of the conditions within the early Universe is the 21-cm spectral line of neutral hydrogen, which encodes in its brightness temperature distribution details of the radiation field and gas properties in the intergalactic medium permeating the cosmos (Furlanetto, Oh, & Briggs 2006; Pritchard & Loeb 2008). Redshifted to low frequencies, the 21-cm line is accessible with radio telescopes (υ < 300 MHz), including current and future instruments. These include the Murchison Widefield Array, MWA^a (Bowman et al. 2013; Tingay et al. 2013; Jacobs et al. 2016); the Precision Array for Probing the Epoch of Reionisation, PAPER^b (Parsons et al. 2010); the LOw-Frequency ARray, LOFAR^c (van Haarlem et al. 2013; Patil et al. 2016); the Long Wavelength Array, LWA^d (Ellingson et al. 2009), and the future HERA (DeBoer et al. 2016) and SKA-Low (Koopmans et al. 2015).

The weakness of the signal, combined with the expectation that most of its information content is contained in the second moment (Wyithe & Morales 2007), which is uncorrelated across spatial Fourier wave mode, motivates the use of the power spectrum as a statistical tool for detecting and characterising the cosmological signal. Despite the ease with which the power spectrum can be computed from radio interferometric data, the presence of strong, spectrally structured residual foreground sources (Trott, Wayth, & Tingay 2012; Datta, Bowman, & Carilli 2010; Vedantham, Udaya Shankar, & Subrahmanyan 2012; Thyagarajan et al. 2015), complex instrumentation (Trott & Wayth 2016), and imperfect calibration (Patil et al. 2014; Barry et al. 2016), yield power spectra that are dominated by systematics. Thus far, a detection of signal from the Epoch of Reionisation (EoR) has not been achieved (Patil et al. 2016; Beardsley et al. 2016; Trott et al. 2016; Cheng et al. 2018). These systematics, combined with the expectation that non-Gaussian information can be extracted usefully from cosmological data, lead the discussion for other statistics. The bispectrum, as a measure of signal non-Gaussianity, is one such statistic that contains cosmologically relevant information (Bharadwaj & Pandey 2005; Majumdar et al. 2018; Watkinson et al. 2018), while being relatively straightforward to compute with interferometric data (Shimabukuro et al. 2017).

The bispectrum is the Fourier Transform of the three-point correlation function, and extracts higher-order correlations between different spatial scales. Its spatial and redshift evolution can be used to place different constraints on the underlying processes that set the 21-cm brightness temperature, and therefore it provides complementary information to the power spectrum. In an early paper exploring the use of the bispectrum for a model EoR signal, and radio interferometers, Bharadwaj & Pandey (2005) demonstrated that a strong non-Gaussian signal is produced by the presence of ionised regions, and discussed the behaviour of the power spectrum and bispectrum signals as a function of frequency channel separation, although they only consider non-Gaussianity due to the ionisation field modelled as non-overlapping randomly placed spherical ionised regions. Some recent work has explored the combination of bispectrum with other tracers (CII spectral features) to extract clean cosmological information (Beane & Lidz 2018). The bispectrum has also been used in the single-frequency (angular) case in the CMB community, where non-Gaussianities can be contaminated by structured foregrounds (Jung, Racine, & van Tent 2018).

Majumdar et al. (2018) explore the ability of the bispectrum to discriminate fluctuations in the matter density distribution from those of the hydrogen neutral fraction, reporting that for some triangle configurations the sign of the bispectrum is a marker for which of these processes is dominating the bispectrum. They show output bispectra for equilateral and isosceles configurations over a range of wavemodes and redshifts, including parameters of relevance to current low-frequency 21-cm experiments (z < 9, 0.1 < k < 1.0). For modes relevant to the MWA, the bispectrum amplitude fluctuates in sign with wavenumber and triangle geometry (stretched → equilateral → squeezed) with a range spanning 10³ − 10⁹ mK³h ⁻⁶ Mpc⁶. This range of potential signs and amplitudes in measurable modes and redshifts motivates us to study this signal in MWA data.

Watkinson et al. (2018) provide a useful tool for visualising the correspondence of real-space structures and bispectrum. They highlight that equilateral k-vector configurations probe above-average signal concentrated in filaments with a circular cross-section (their Figure 1). Stretched (flattened) k-vector triangle configurations (with one k-mode larger than the other two), by extension, probe above-average signal concentrated in filaments with ellipsoidal cross-sections (at the extreme these filaments tend towards planes). Finally, squeezed k-vector triangle configurations (with one k-mode smaller than the other two) correspond to a modulation of a large-scale mode over small-scale plane-wave concentrations of above-average signal, and therefore measure the correlation of the small-scale power spectrum with large-scale modes. Notably, they introduce and explore other bispectrum normalisations that are found to be more stable to parameter fluctuations. In this work, we discuss the relative merits of different bispectrum statistics for use with real data in the presence of real systematics.

Crucially, the switch to positive bispectrum at the end of reionisation occurs as we reach regimes/scales at which the concentration of above-average signal drive the non-Gaussianity. This will occur before the EoR (on scales where the density field is the dominant driver of the temperature fluctuations, or, if the spin temperature is not yet saturated during this phase, when heated regions are driving the non-Gaussianity) and towards the end of reionisation (when islands of 21-cm signal drive the non-Gaussianity).

Conversely, a negative-valued bispectrum will be unique to the phase when ionised regions drive the non-Gaussianity. In general, foreground astrophysical processes are not expected to produce a negative bispectrum, because they are associated with overdensities in the brightness temperature distribution (Lewis 2011; Watkinson & Pritchard 2014). These factors may play a future important role in discriminating real cosmological non-Gaussianity from contaminants.

Despite some work studying the sensitivity of current and future experiments for measuring the bispectrum (Shimabukuro et al. 2017; Yoshiura et al. 2015), these have used idealised scenarios that omit any residual foreground signal and systematics introduced by the instrument. Bharadwaj & Pandey (2005) discuss foreground fitting tools using frequency separation to study the bispectrum over visibility correlations across frequency, but this method breaks down for large field-of-view instruments where the interferometric response affects the foreground smoothness (Morales et al. 2012). Further, no 21-cm interferometric data have been used to estimate the bispectrum. In this work, we address both of these by presenting bispectrum estimators that can use real datasets, computing the expected impact of foregrounds measured by the instrument, and applying the estimators to 21 h of MWA EoR data.

2. MWA Phase II Array

The MWA is a 256-tile low-frequency radio interferometer located in the Western Australian desert, on the future site of the Square Kilometre Array (SKA) (Tingay et al. 2013; Bowman et al. 2013). The telescope operates from 80 to 300 MHz with antennas spread over a 5-km diameter. Its primary science areas include exploration of the EoR, radio transients, solar and heliospheric studies, study of pulsars and fast transients, and the production of a full-sky low-frequency extragalactic catalogue. In 2016 it underwent an upgrade from 128 to 256 antenna tiles (Wayth et al. 2018). At any time, 128 of the tiles can be connected to the signal processing system. The array operates in a ‘compact’ configuration, utilising redundant spacings and short baselines for EoR science, or an ‘extended’ configuration, maximising angular resolution and instantaneous uv-coverage. The compact configuration is employed in this work.

Figure 1. Zoomed MWA compact configuration layout showing the two hexagonal sub-arrays of 36 tiles each, with redundant tile spacings. These short redundant baselines are used in this work to form equilateral and isosceles triangle bispectra with high sensitivity. Some of the longer baseline tiles of the MWA are not shown.

The compact configuration has a maximum baseline of 500 m and is optimised for EoR science. Figure 1 shows the tile layout, including the two 36-tile hexagonal subarrays of redundantly spaced tiles. The minimum redundant spacing is 14 m. The primary motivations for the hexagons are twofold: (1) to increase the sensitivity to angular scales of relevance for the EoR, allowing coherent addition of measurements from redundant baselines, and (2) enabling additional methods for calibrating the array (redundant calibration, Li et al. 2018, Joseph et al. 2018). For the bispectrum, there is an additional advantage of multiple, redundant equilateral triangle baselines being formed from the short spacings. These can be added coherently to study the bispectrum signal on particular scales, and allows for a direct bispectrum measurement (perfectly defined triangles formed from discrete baselines). These direct bispectrum results can be compared to a more general gridded bispectrum, whereby all baselines formed by an irregularly spaced array (such as MWA Phase I, or the non-hexagon tiles of Phase II compact) can be gridded onto the Fourier (uv-) plane, using a gridding kernel that represents the Fourier response function of the telescope (in this case, the Fourier Transform of the primary beam response to the sky). These estimators will both be explored in this work.

3. Power spectrum

We briefly review the power spectrum as the primary estimator for studying the EoR with 21-cm observations. The power spectrum is typically used to describe radio interferometer observations from the EoR, and contains all of the Gaussian-distributed fluctuation information. The power spectrum is the power spectral density of the spatial fluctuations in the 21-cm brightness temperature field. It is used because it encodes the fluctuation variance (where most of the EoR signal is expected to reside), and sums signal from across the observing volume to increase sensitivity. It is defined as

(1) $$P(\vec k) = {\delta _D}(\vec k - {\vec k^{'}}){1 \over {{\Omega _V}}}\langle {V^{ *} }(\vec k)V({\vec k^{'}})\rangle ,$$

where $V(\vec k) = V(u,v,\eta ) = {\cal F}{\cal T}\left( {V(u,v,\nu )} \right)$ is the measured interferometric visibility (Jansky), Fourier transformed along frequency (υ) to map frequency to line-of-sight spatial scales (Jy Hz) at a given point in the Fourier (uv-) angular plane (u, v); 〈〉 encode an ensemble average over different realisations of the Universe, and the δ_D-function ensures that we are expecting to measure a Gaussian random field where the different modes are uncorrelated^e. Further assuming spatial isotropy allows us to average incoherently in spherical shells, where $\vec k = k$. Ω_V provides the volume normalisation, where Ω_V = (BW) Ω is the product of the observing bandwidth and angular field of view. Converting from measured to physical units maps Jy² Hz² to mK² h ⁻⁶ Mpc⁶. After volume normalisation this becomes, mK² h ⁻³ Mpc³.

3.1. Power spectra with radio interferometric data

The power spectrum can be produced naturally with interferometric data. Unlike optical telescopes that produce images of the sky, or single-dish radio telescopes that acquire a single sky power, a radio interferometer visibility (Jy) directly measures Fourier representations of the sky brightness distribution at the projected baseline location (u = Δx/λ; v = Δy/λ). In the flat-sky approximation^f:

(4) $$V(u,v,\nu ) = \int_\Omega A(l,m,\nu )S(l,m,\nu )\exp ( - 2\pi i(ul + vm))dldm,$$

where A(l, m, υ) is the instrument primary beam response to the sky at position (l, m) from the phase centre and frequency υ, S(l, m, υ) is the corresponding sky brightness (Jy/sr, which is proportional to temperature), and the exponential encodes the Fourier kernel. The physical correspondence of sky projected on to the tile locations yields a fixed set of discrete but incomplete Fourier modes to be measured. This incompleteness leads to parts of the Fourier plane where there is no information. The line-of-sight spatial scales are obtained by Fourier Transform of visibilities measured at different frequencies, along frequency to map υ to η:

(5) $$V(\eta (k)) = {\cal F}{\cal T}(V(\nu )) = {{\Delta \nu } \over {{N_{{\rm{c}}h}}}}\sum\limits_{j = 1}^{{N_{{\rm{c}}h}}} V(\nu )\exp \left( { - {{2\pi ijk} \over {{N_{{\rm{c}}h}}}}} \right),$$

where N _ch is the number of spectral channels, Δυ is the spectral resolution, and j and k index frequency and spatial mode (Hz⁻¹, or seconds).

The attenuation of the sky due to the primary beam (and general sky finiteness) alters the complete continuous Fourier Transform to a windowed transform, whereby the primary beam response leaks signal into adjacent Fourier modes, as can be seen using the convolution theorem:

(6) \begin{equation} V(u,v,\eta) = {\kern2pt\tilde{\kern-2pt A}}(u,v,\eta) \circledast \tilde{S}(u,v,\eta), \end{equation}

where the true sky brightness distribution is convolved with the Fourier Transform of the primary beam response. This leakage implies that the visibility measured by a discrete baseline actually contains signal from a region of the Fourier plane, as described by the Fourier beam kernel, Ã (u, v, η).

In general, to compute the power spectrum from a large amount of data, we are motivated by the signal weakness to add the data coherently; i.e., we sum complex visibilities directly that contribute signal to the same point in the Fourier uv-plane. To do this, the measurement from each baseline is convolved with the Fourier beam kernel and ‘gridded’ (added with a weight) onto a common two-dimensional plane. Signal will add coherently, while noise adds as the square-root (because the thermal noise is uncorrelated between measurements). The weights for each measurement are also gridded with the kernel onto a similar plane. After addition of all the data, the signal uv-plane is divided by the weights to yield the optimal-weighted average signal at each point. The resulting cube resides in (u, v, υ) space, and can be Fourier transformed along frequency to obtain a cube in (u, v, η) space. The power spectrum can then be formed by squaring and normalising the cube, and averaging incoherently (in power) in spherical shells:

(7) $$P(|\vec k|) = {{\sum\limits_{i \in k} V_i^ * (\vec k){V_i}(\vec k){W_i}(\vec k)} \over {\sum\limits_{i \in k} {W_i}(\vec k)}},$$

where W are the weights and $|\vec k| = |({k_u},{k_v},{k_\eta })| = \sqrt {k_u^2 + k_v^2 + k_\eta ^2} $.

As an intermediate step, the cylindrically averaged power spectrum can be formed (e.g., Datta et al. 2010):

(8) $$P({k_ \bot },{k_\parallel }) = {{\sum\limits_{i \in {k_ \bot }} {V^ * }(\vec k)V(\vec k)W(\vec k)} \over {\sum\limits_{i \in {k_ \bot }} W(\vec k)}},$$

and ${k_ \bot } = \sqrt {k_u^2 + k_v^2}, \ k_\parallel=k_\eta$. This is a useful estimator for discriminating contaminating foregrounds (continuum sources with power concentrated at small k _||) from 21-cm signal. Herein we will refer to this power spectrum, and its bispectrum analog, as the ‘gridded power spectrum’ and ‘gridded bispectrum’, respectively.

Alternatively, one can take the baselines themselves, and their visibilities measured along frequency, and take the Fourier Transform directly along the frequency axis. This ‘delay spectrum’ approach is utilised by some experiments with short baselines (Parsons et al. 2012; Ali et al. 2015; Thyagarajan et al. 2015), both to increase sensitivity when there are redundant spacings, and to work as a diagnostic. The frequency and η axes are not parallel, except at zero-length baseline. Because an interferometer is formed instantaneously from antennas with a fixed spatial offset, the baseline length in Fourier space (e.g., u) evolves with frequency as u = Δxυ/c, and this evolution is therefore increased for larger bandwidths and for longer baselines. For the short spacings of interest to the EoR, the correspondence is good, and the delay transform can be used to mimic the direct k _|| transform of gridded data (see Figure 1 of Morales et al. 2012, for a visual explanation). In general, ‘imaging’ arrays with many non-redundant spacings are suited to gridded power spectra, whereas redundant arrays, with a lesser number of multiply-sampled modes, are suited to delay power spectra. For the Phase II compact MWA, the two hexagonal subarrays have these short-spaced redundant baselines, and the ‘delay power spectrum’ and its bispectrum analog can also be used effectively. In general, we would not suggest use of the delay spectrum to undertake EoR science, because of the limitations discussed, but it is the appropriate analogue for the direct bispectrum estimator, and is therefore pertinent for the normalised bispectrum analysis.

4. Bispectrum

The bispectrum is the Fourier Transform of the three-point correlation function. Akin to the two-point correlation function (the Fourier dual of which is the power spectrum), the three-point correlation function measures the excess signal over that of a Gaussian random field distribution measured at three spatial locations, averaged over the volume. For a field with Fourier Transform denoted by $\Delta (\vec k)$, the bispectrum is formed over closed triangles of k vectors in Fourier space:

(9) $$\langle \Delta ({\vec k_1})\Delta ({\vec k_2})\Delta ({\vec k_3})\rangle = {\delta _D}({\vec k_1},{\vec k_2},{\vec k_3}){\rm{B}}\,({{\rm{\tilde k}}_1},{{\rm{\tilde k}}_2},{{\rm{\tilde k}}_3}).$$

Here the δ_D-function ensures closure in Fourier space. It has units of mK³ h ⁻⁶Mpc⁶ after volume normalisation. The bispectrum is often applied to matter density fields, where $\Delta (\vec k)$ is the Fourier Transform of matter overdensity, $\delta (\vec x) = {{\rho (\vec x)} \over {\overline \rho }} - 1$. In radio interferometric measurements, the coherence of the wavefront (the visibilities obtained by cross-correlating voltages from individual antennas) represents the Fourier Transform of the sky brightness temperature distribution, measured in Jansky.

As discussed earlier, this bispectrum estimator can be unstable, with cosmological simulations showing rapid fluctuations between positive and negative values as non-Gaussianity becomes negligible but the amplitude is still large. As such, Watkinson et al. (2018) suggest the normalised bispectrum as a more stable statistic:

(10) $${\cal B}({\vec k_1},{\vec k_2},{\vec k_3}) = {{{\rm{B}}\,({{{\rm{\tilde k}}}_1},{{{\rm{\tilde k}}}_2},{{{\rm{\tilde k}}}_3})\sqrt {{{\rm{k}}_1}{{\rm{k}}_2}{{\rm{k}}_3}} } \over {\sqrt {P({{\vec k}_1})P({{\vec k}_2})P({{\vec k}_3})} }},$$

where $P(|\vec k|)$ is the three-dimensional power spectrum, which describes the volume-normalised variance on a given spatial scale, and is the Fourier Transform of the two-point correlation function (Eggemeier & Smith 2017; Brillinger & Rosenblatt 1967). This normalisation isolates the contribution from the non-Gaussianity to the bispectrum, by normalising out the amplitude part of the statistic. It is akin to normalising the third central moment by σ ³ to calculate the skewness.

4.1. Bispectrum with radio interferometric data

Because the bispectrum is formed from the triple product of a triangle of wavespace measurements, it can be formed directly through the product of three interferometric visibilities. In the limit where the array has perfect (complete) uv-sampling, individual measurements of signal on triangles of baselines can be multiplied to form the bispectrum estimate. In the more general case, where an interferometer has instantaneously incomplete, but well-sampled baselines, there are two options for extracting the triangles of signal measurements: direct (via multiplication of measurements from three tiles forming a triangle of baselines), or gridded, where each uv-measurement is gridded onto the uv-plane (with its corresponding Fourier beam kernel and weights), and the final bispectra are computed from the fully-integrated and gridded data.

Direct bispectrum estimators can be applied to specific triangles according to the array layout, but these are usually unique, with irregular configurations (all three internal angles are distinct), leading to difficult cosmological interpretation and poor sensitivity. These issues arise for imaging-like arrays with pseudo-random layouts, but are alleviated for redundant arrays, where regular triangles (isosceles and equilateral) exist and are instantaneously available in many copies in the array. These features make interpretation more straightforward and increase sensitivity to these bispectrum modes.

Gridded bispectrum estimators can be applied to any array, yield improved sensitivity by coherent gridding of data, and may allow for a wider range of triangles to be probed. Nonetheless, they suffer from the increased difficulty of extracting robust estimates that correctly account for the correlation of data in uv-space.

With the benefit of having a redundant array, we will apply both sets of estimators to our data.

5. The gridded estimator

Each measured visibility encodes information about a small range of Fourier modes of the sky brightness distribution. Although each baseline is usually reported as a single number representing the antenna separation measured between antenna centres, the baselines actually measure a range of separations when accounting for the actual physical size.^g This translates to a range of Fourier modes being measured by a given baseline, and is equivalent to the statement that a finite primary beam response to the sky mixes Fourier modes through spectral leakage (effectively a taper on the continuous Fourier transform). Thus, when measurements from different baselines are combined coherently (with phase information) onto a uv-plane, they can be gridded with a kernel that is the Fourier Transform of the primary beam response to the sky. Such a gridding kernel captures the degree of spectral leakage introduced by the antenna response, and means that baselines of similar length and orientation have some shared information. The gridding kernel is represented by à in Equation (6). With a single defined visibility phase centre, all visibility measurements can be added with this beam kernel onto a single plane (for each frequency channel), along with their associated weights, to form a coherently averaged estimate for the Fourier representation of the sky brightness temperature:

(11) $${\hat V_{uv}} = {{\sum\limits_i V({u_i},{v_i})\tilde A({u_i},{v_i})W({u_i},{v_i})} \over {\sum\limits_i \tilde A({u_i},{v_i})W({u_i},{v_i})}},$$

where i indexes measurement and W is the weight associated with each. The bispectrum is then estimated as the sum over the beam-weighted gridded visibilities:

(12) $${\hat B_{123}} = {{\sum\limits_{j \in \tilde A} {{\hat V}_{j1}}{{\hat V}_{j2}}{{\hat V}_{j3}}{W_{j1}}{W_{j2}}{W_{j3}}} \over {\sum\limits_j {W_{j1}}{W_{j2}}{W_{j3}}}},$$

where

(13) $${W_{1j}} = {W_1}{\tilde A_{1j}},$$

are the beam-gridded measurement weights.

5.1. Gridded estimator noise

The gridded bispectrum estimator is formed from coherent addition of visibilities over all observations. As such, if a given visibility has thermal noise level σ _therm (Jy Hz)^h, the uncertainty on the bispectrum is

(14) \begin{equation} \Delta\hat{B}_{\rm TOT} = \frac{\sqrt{3}\sigma_{\rm therm}^3}{\sqrt{\displaystyle\sum_{{\kern2pt\tilde{\kern-2pt A}}} W_{1}W_{2}W_{3}}}, \end{equation}

where the denominator is the sum over the gridding kernel of the weights triplets.

The uncertainty on the normalised bispectrum is then:

(15) $$\Delta {\hat {\cal B}_{123,{\rm{TOT}}}} = {\cal B}\sqrt {{{\Delta {B^2}} \over {{B^2}}} + {{\Delta P_1^2} \over {4P_1^2}} + {{\Delta P_2^2} \over {4P_2^2}} + {{\Delta P_3^2} \over {4P_3^2}}} ,$$

where the uncertainties can contain both thermal noise and noise-like uncertainty from residual foregrounds.

For 300 2-min observations, and 24 triangles per 28 m baseline triad group, the expected thermal noise level for a complete dataset is

(16) $$\Delta B = 4.2 \times {10^{10}}{\rm{m}}{K^3}{h^{ - 6}}{\rm{M}}p{c^6}.$$

The presence of residual foregrounds will be studied in Section 10.2.

6. The direct estimator

As an alternate approach to the gridded estimator, visibilities are Fourier-transformed along the frequency direction to compute the delay transform, and closed bispectrum triangles formed from the closed redundant triads of antennas. This approach does not use the primary beam and ignores the local spatial correlations generated by the primary beam spatial taper. It also transforms along a dimension that changes angle with respect to k _|| as a function of baseline length, but approximates a k _|| Fourier Transform for small u (small angle).

The bispectrum for a given observation is the weighted average over all triads:

(17) $${\hat B_{123}} = {{(\sum\limits_i {V_{1i}}{W_{1i}})(\sum\limits_i {V_{2i}}{W_{2i}})(\sum\limits_i {V_{3i}}{W_{3i}})} \over {(\sum\limits_i {W_{1i}})(\sum\limits_i {W_{2i}})(\sum\limits_i {W_{3i}})}},$$

where i indexes over redundant triangles (triads). The final bispectrum estimate then performs a weighted average over observations, such that:

(18) $${\hat B_{123,{\rm{TOT}}}} = {{\sum\limits_j {{\hat B}_{123,j}}{W_j}} \over {\sum\limits_j {W_j}}},$$

where

(19) $${W_j} = (\sum\limits_i {W_{1i}})(\sum\limits_i {W_{2i}})(\sum\limits_i {W_{3i}}).$$

Figure 2. Schematic of how a stretched isosceles triangle configuration is extracted from redundant angularly-equilateral triangle baselines of the MWA Phase II hexagons.

6.1. Direct Estimator noise

The direct bispectrum estimator is formed from coherent addition of baseline triplets for a given observation, which are then averaged with relative weights to the final estimate. As such, if a given visibility has thermal noise level σ _therm (Jy Hz), the uncertainty on the bispectrum is

(20) $$\Delta {\hat B_{{\rm{TOT}}}} = {{\sqrt 3 \sigma _{{\rm{therm}}}^3} \over {\sqrt {\sum\limits_j {W_j}} }}.$$

The uncertainty on the normalised bispectrum is then given by the same expression as for the Gridded Estimator [Equation (15)].

For 300 observations, and 24 triangles per 28 m baseline triad group, the expected noise level is

(21) $$\Delta B = 7.1 \times {10^{11}}{\rm{m}}{K^3}{h^{ - 6}}{\rm{M}}p{c^6}.$$

7. Triangles considered for estimation

Unlike bispectrum estimates that can be obtained from Phase I data, where the array is in an imaging configuration with no redundant triangles, we aim to take advantage of the 72 redundant tiles in the hexagonal sub-arrays, afforded by the Phase II layout. This allows for both direct and gridded bispectrum estimators to be applied to matched observations with matched data calibration.

The most numerous (highest sensitivity) groups of redundant triangles are the angularly-equilateral configurations of the 14 m and 28 m baselines (48 and 24 sets, respectively). For these triangles, the equilateral configurations exist only for the η = 0 (k _|| = 0) line-of-sight mode. Other configurations of these closed angular triangles are isosceles or irregular triangles, depending on the η values chosen; however, the closed triangle requirement of the bispectrum demands that:

(22) $${\eta _1} + {\eta _2} + {\eta _3} = 0,$$

in addition to the angular components of the vectors summing to zero (as is enforced by choosing the closed triangle baselines).

For comparison with theoretical predictions, we will focus on equilateral and isosceles triangles. The 14 m and 28 m baselines are very short, corresponding to cosmological scales of k _⊥ ≃ 0.01hMpc⁻¹ at z = 9. Thus, although the equilateral configuration is cosmologically relevant and the easiest to interpret, these modes are expected to be heavily foreground dominated (i.e., they correspond to the line-of-sight DC mode, and the large angular scales of diffuse and point-source foreground emission). We consider them for completeness, but will show them to be cosmologically irrelevant from an observational perspective when computed this way. These same angularly-equilateral triangle configurations will, however, be used to form relevant isosceles configurations with η ₁ = η ₂ and η ₃ = −2 η ₁. Given that we aim to sample modes where foregrounds are not dominant in our power spectra, these isosceles configurations form ‘stretched’ (also referred to ‘flattened’ in Watkinson et al. 2018) configurations (k _|| > > k _⊥). Figure 2 shows how the stretched isosceles configurations are extracted from the data with a redundant baseline triad. Figure 3 then shows schematically the approximate vectors for two of the four isosceles configurations considered here.

8. Observations

The direct and gridded estimators are applied to 21.0 h of Phase II high-band zenith-pointed data, comprising 10.7 h (320 observations) on the EoR0 field (RA = 0 h, Dec. = −27°) and 10.3 h (309 observations) on the EoR1 field (RA = 4 h, Dec. = −30°). We observe 30.72 MHz in 384 contiguous 80 kHz channels, with a base frequency of 167.035 MHz. Approximately 15% of the observations were obtained from drift-scan data, where the telescope remains pointed at zenith for many hours and the sky drifts through. For consistency with the drift-n-shift data, we chose drift scan data observed with the field phase centres within 3° of zenith. The data were observed over 5 weeks from 2016 October 15 to November 28, and 1 week in 2017 July. Because the delay spectrum is used as part of the power spectrum estimator for the direct bispectrum, each observation was individually inspected for poor calibration or data quality, and bad observations excised from the dataset. The excised observations comprised ∼ 5% of the dataset, and primarily were due to poor calibration solutions over sets of data contiguous in time due to poor instrument conditions (e.g., many flagged tiles or spectral channels).

Figure 3. Schematic of how isosceles triangle vectors are extracted, overlaid on a power spectrum. We aim to choose triangles with vectors that reside in noise-like regions of the delay spectrum.

The 2-min observations were each calibrated through the MWA Real-Time System (RTS; Mitchell et al. 2008), as is routinely performed for MWA EoR data, and 1 000 of the brightest (apparent) sources peeled from the dataset (Jacobs et al. 2016). These 629 calibrated and peeled observations were used for bispectrum estimation.

9. Results

We begin by reporting the bispectrum estimates for the two methods and fields, and then report the normalised bispectra, which incorporate the power spectrum estimates. Table 1 shows the bispectrum estimates and their one-sigma uncertainties (thermal noise) for the direct and gridded estimators, both observing fields and for different triangle configurations. Bold-faced results indicate bispectrum estimates that are consistent with thermal noise. These tend to be those that are extremely stretched isosceles configurations (cos θ ∼ 1), with estimates that sit well outside the primary foreground contamination parts of parameter space. Conversely, the equilateral triangle configurations that use the k _|| = 0 mode exclusively show extremely large detections. There is no suggestion that these are 21-cm cosmological bispectrum detections, but rather are foreground contaminants. This will be explored more fully in Section 10. Note also that the thermal noise levels reported here are a factor of a few larger than the theoretical expectation derived in Sections 5 and 6. This is due to the fraction of data with weights that are less than unity, indicating flagged baselines and spectral channels.

Also listed in Table 1 is the expected bispectrum values from simulations that assume either bright or faint galaxies drive reionisation (Greig & Mesinger 2017). The largest amplitudes are for the smallest k-modes, which also tend to be more foreground dominated.

The normalised bispectrum, $\cal {B}$, is normalised by the power spectra at each of the k modes forming the triangles. Figure 4 shows the power spectra for the EoR1 and EoR0 fields for the full datasets as used in the gridded estimator. These have been processed through the CHIPS power spectrum estimator (Trott et al. 2016). Figures 5 and 6 show the corresponding delay spectra, as used in the direct bispectrum. There are small differences between the two power spectrum estimators, as is expected given that delay spectra do not grid with primary beams, and Fourier Transform along frequency, yielding different results for longer baselines. The signature of Galactic emission from close to the horizon is evident in the EoR0 power spectra, while it is less structured in EoR1, where the Galactic Centre has set. Most notably, the delay spectra show large foreground leakage into the EoR window (k _|| < 0.4), yielding large power spectrum denominator values for the normalised bispectrum.

Table 1. Bispectrum estimates and one-sigma uncertainties for the direct and gridded bispectra for each observing field and triangle type. Bold-faced values indicate bispectrum estimates that are consistent with thermal noise. The right-hand column lists expected bispectrum values from simulation for faint and bright galaxies driving reionisation. k modes are comoving and measured in h Mpc⁻¹.

Figure 4. Gridded power spectra for the 21 h of observations on two fields used in this work, as processed through the CHIPS estimator (Trott et al. 2016).

Figure 5. Delay transform power spectra for the EoR1 field for the data used in this analysis. Note the large leakage into the EoR window, which yields large denominators for the normalised bispectrum.

Figure 6. Delay transform power spectra for the EoR0 field for the data used in this analysis.

Using these data, Table 2 describes the normalised bispectrum. Bold-faced results are broadly consistent with thermal noise (< 5σ), again reflecting the modes that are least affected by foregrounds. The difference between the dimensional and reduced bispectrum results is due to the different power spectral estimators. Also notable is the difference in amplitude of the gridded and direct normalised bispectrum estimates. Due to the division by the power spectrum, the normalised bispectrum is heavily dependent on the details of the power spectrum estimates, which fluctuate substantially in foreground-affected regions. The delay-space power spectra show increased foreground power in the EoR window, and this is reflected in a larger power spectrum estimate, and therefore a lower normalised bispectrum. This reliance highlights the complexity for interpreting the normalised bispectrum with foreground-affected data.

10. Bispectrum signature of foregrounds

Estimates of bispectrum sensitivity for operational and future 21-cm experiments are incomplete without a treatment of foregrounds. Despite the expectation that point source, continuum foregrounds only impact a region of the three dimensional EoR parameter space (k_x, k_y, k _||), in reality the details of the instruments, complexity of extragalactic and Galactic emission, limited bandwidth and calibration errors leave residual contaminating signal throughout the full parameter space. Although these methods perform very well to remove such signal, the extreme dynamic range demanded by this experiment translate to bias that exceeds the expected cosmological signal strength. The results presented here are clearly foreground-dominated, particularly for the equilateral triangle configuration.

As such, the bispectrum signature of foregrounds can be computed for a simple point-source foreground model. We first consider the expected foreground bispectrum, which quantifies the bias in the measurement, and then turn to the variance of the foreground bispectrum, which quantifies the additional noise term.

We employ a model where the sky is populated with a random distribution of unresolved extragalactic point sources that follow a low-frequency number counts distribution (Intema et al. 2011; Franzen et al. 2016):

(23) $${{dN} \over {dS}} = \alpha {S^\beta }{\kern 1pt} {\kern 1pt} {\rm{J}}{y^{ - 1}}{\rm{s}}{{\rm{r}}^{ - 1}},$$

where α ≃ 3900 and β = −1.59 for sources with flux density at 150 MHz of less than 1 Jansky. We assume there is no source clustering and spectral dependence, yielding a Poisson-distributed number of sources in each differential sky area.

The clustering of point sources in the power spectrum has been studied by Murray, Trott, & Jordan (2017). They find that source clustering will be unimportant for the MWA (unless the clustering is extreme, which is not measured), but may be important for the SKA, which can clean to deeper source levels. Nonetheless, the structure due to clustered point-source foregrounds only changes the amplitude of the foreground structure in the EoR wedge as a function of angular scale (k _⊥). Because the line-of-sight spectral component is unaffected, the signature of clustered foregrounds in the EoR Window is mostly unchanged. These more realistic point-source foregrounds will be considered in the simulations of Watkinson & Trott (2019) and here we retain the analytic signature of the Poisson foregrounds.

We further assume that the primary beam can be approximated by a frequency-dependent Gaussian:

(24) $$A(l,m,{\nu _0}) = \exp \left( { - {{({l^2} + {m^2})\nu _0^2{{\rm{A}}_{\rm{eff}}}} \over {2{c^2}{\epsilon^2}}}} \right),$$

Table 2 Normalised bispectrum estimates, ${\cal B}$, and one-sigma uncertainties for the direct and gridded bispectra for each observing field and triangle type. Bold-faced values indicate bispectrum estimates that are consistent with thermal noise.

where A_eff is the tile effective area and ϵ encodes the conversion from an Airy disc to a Gaussian.

The visibility is given by Equation (4) for frequency υ. To compute the line-of-sight component to the visibility, we Fourier Transform over frequency channels, after employing a frequency taper (window function) to reduce spectral leakage from the finite bandwidth:

(25) \begin{array}{*{20}{l}} {V(u,v,\eta )}& = &{\int d ldmS(l,m,{\nu _0})A(l,m,{\nu _0})}\\ {}&{}&{ \times \int d \nu \Upsilon (\nu )\exp \left( { - 2\pi i\frac{{\nu (xl + ym)}}{c}} \right)\exp \left( { - 2\pi i\nu \eta } \right)}\\ {}& = &{\int d ldmS(l,m,{\nu _0})A(l,m,{\nu _0})}\\ {}&{}&{ \times \int d \nu \Upsilon (\nu )\exp \left( { - 2\pi i\nu (xl/c + ym/c + \eta )} \right)} \end{array}

(26) \begin{array}{*{20}{l}} = &{\int d ldmS(l,m)A(l,m)}\\ {}&{ \times \;\tilde \Upsilon (xl/c + ym/c + \eta )\,{\rm{J}}yHz,} \end{array}

whee γ(υ) is the spectral taper, Δυ. is the channel resolution, and we have performed the Fourier Transform over frequency. For analytic tractability, in this work we use a Gaussian taper, with a characteristic width, Σ ≃BW/7, such that the edges of the band are consistent with zero and it is well-matched to a Blackman–Harris taper:

(27) $$\Upsilon (\nu ) = {1 \over {\sqrt {2\pi {\Sigma ^2}} }}\exp - {{{\nu ^2}} \over {2{\Sigma ^2}}},$$

with corresponding Fourier Transform,

(28) $$\tilde \Upsilon (\eta ) = \exp - 2{\pi ^2}{\Sigma ^2}{\eta ^2}.$$

The bispectrum is formed from the triple product of visibilities. Accounting for the fact that the point sources are only correlated locally (δ_D(l ₁ + l ₂ + l ₃ = 0)), its expected value with respect to foregrounds is

(29) $$\langle {V_1}{V_2}{V_3}\rangle = \int dldm\langle {S^3}(l,m)\rangle {A^3}(l,m)\exp \left( { - 2{\pi ^2}{\Sigma ^2}{T^2}} \right)$$

where

(30) $$\matrix{ {{T^2} = {{\left( {{{{x_1}l} \over c} + {{{y_1}m} \over c} - {\eta _1}} \right)}^2}} \hfill \cr { + {{\left( {{{{x_2}l} \over c} + {{{y_2}m} \over c} - {\eta _2}} \right)}^2} + {{\left( {{{{x_3}l} \over c} + {{{y_3}m} \over c} - {\eta _3}} \right)}^2}.} \hfill \cr } $$

Here, the source counts have been separated from the spatial integral. This is a general expression for a triplet of baselines. We can now simplify this for triangles, particularly those with isosceles configurations (where the equilateral is a single case of an isosceles).

Closed triangles follow the relations:

(31) $${x_1} + {x_2} = - {x_3}$$

(32) $${y_1} + {y_2} = - {y_3}$$

(33) $${\eta_1} + {\eta_2} = - {\eta_3}$$

and we define, without loss of generality, the following relations for the isosceles configurations considered in this work:

(34) \begin{array}{*{20}{l}} {{x_1}}& = &{ - 2{x_2}}\\ {{x_2}}& = &{{x_3}}\\ {{y_1}}& = &0\\ {{y_2}}& = &{ - {y_3}}\\ {{y_2}}& = &{{x_1}\cos \pi /6 = 2{x_2}\cos \pi /6 = \sqrt 3 {x_2}}\\ {2{\eta _2}}& = &{2{\eta _3} = - {\eta _1}.} \end{array}

Making these substitutions in Equation (29), completing the squares and collecting terms, we find:

(35) \begin{array}{*{20}{c}} {\langle {V_1}{V_2}{V_3}\rangle = \int d ldm\langle {S^3}(l,m)\rangle {A^3}(l)}\\ { \times \exp \left( { - 12{\pi ^2}{\Sigma ^2}\left( {x_2^2/{c^2}({l^2} + {m^2}) + \eta _2^2} \right)} \right).} \end{array}

The source count expectation value uses the source number counts distribution and the fact that the number of sources at any sky location is Poisson-distributed to find:

(36) \begin{array}{*{20}{l}} {\langle {S^3}(l,m)\rangle }& = &{\int_S {{S^3}} ({\nu _0})\frac{{dN}}{{dS}}}\\ {}& = &{\frac{\alpha }{{4 + \beta }}S_{{\rm{m}}ax}^{4 + \beta }\,{\rm{J}}{y^3}s{r^{ - 1}}} \end{array}

Incorporating the primary beam from Equation 24, moving to polar coordinates, and performing the integral over (l, m), we find for the expected foreground bispectrum bias:

(37) \begin{array}{*{20}{c}} {\langle {V_1}{V_2}{V_3}\rangle = \frac{\alpha }{{4 + \beta }}S_{{\rm{m}}ax}^{4 + \beta }}\\ { \times \frac{\pi }{\theta }\exp \left( { - \frac{{{\pi ^2}{\rm{B}}{{\rm{W}}^2}\eta _2^2}}{{25}}} \right)\,{\rm{J}}{y^3}H{z^3}} \end{array}

where

(38) $$\theta = {{3{{\rm{A}}_{\rm{eff}}}\nu _0^2} \over {{c^2}}} + {{{\pi ^2}{\rm{B}}{W^2}u_2^2} \over {25\nu _0^2}},$$

and BW is the experiment bandwidth. This factor combines the primary beam (spatial taper) and spectral taper components into a single factor.

The equilateral configuration can be derived from this expression with η ₂ = 0. For the 28-m baselines, a maximum source flux density of 1 Jy and A_eff = 21 m², and performing the cosmological conversions, we expect a bispectrum estimate of:

(39) $$B(x = 28) \simeq 8.6 \times {10^{19}}{\rm{m}}{K^3}{h^{ - 6}}{\rm{M}}p{c^6},$$

which is comparable to the estimates found in Section 9. For 14-m baselines, B(x = 14) ≃ 1.0 ×10²⁰mK³h ⁻⁶Mpc⁶.

The isosceles configurations incorporate the η term. For k _|| > 0.1, this term decays to below the noise, which is consistent with that observed in the data. The signature of this isosceles foreground dimensional bispectrum in k _⊥ - k _|| space is shown in Figure 7. For the k ₁ = 0.1 h Mpc⁻¹ stretched configuration, we expect for 28 m (14 m):

(40) $$B \simeq 1.7 \times {10^{12}}{\kern 1pt} {\kern 1pt} (1.0 \times {10^{12}}){\rm{m}}{K^3}{h^{ - 6}}{\rm{M}}p{c^6}.$$

The squeezed configurations of large k _⊥ combined with small k _|| might be interesting for future studies, depending on the expected cosmological signal on these scales. Given that the power spectrum is expected to be small on these combination of line-of-sight and angular scales, most EoR experiments are not designed for high sensitivity here (k _⊥ = 0.1 corresponds to 200-m baselines).

Figure 7. Point-source-foreground dimensional bispectrum signature of isosceles triangle vectors in k _⊥ − k _||-space (note the stretched logarithmic colour bar). In this model, the expected foreground signal has fallen to below the expected cosmological signal value by k _|| ≥ 0.12.

Interestingly, the expected foreground bispectrum signal is positive, due to its constituent astrophysical sources being associated with overdensities. Conversely, the stretched isosceles 21-cm bispectrum from the cosmological signal will be negative on many scales during reionisation (Majumdar et al. 2018).

10.1. Normalised foreground bispectrum

The normalised bispectrum also contains the expected power spectrum values for a foreground model. In line with the methodology developed in the previous section, we can write the expected power spectrum at (u, v, η) as

(41) \begin{array}{*{20}{l}} {P(u,v,\eta )}& = &{\langle {V^*}(u,v,\eta )V(u,v,\eta )\rangle }\\ {}& = &{\frac{\alpha }{{3 + \beta }}S_{{\rm{m}}ax}^{3 + \beta }} \end{array}

(42) $${ \times \left( {{\rm{erf}}\left( {\frac{{b + 2a}}{{\sqrt 2 a}}} \right) - {\rm{erf}}\left( {\frac{{b - 2a}}{{\sqrt 2 a}}} \right)} \right)}$$

(43) $${ \times \exp - 4{\pi ^2}{\Sigma ^2}{\eta ^2}\sqrt {\frac{\pi }{{4a}}} \exp \frac{{{b^2}}}{{4a}},}$$

where

(44) $$a = \frac{{2\pi {c^2}}}{{\nu _0^2{{\rm{A}}_{\rm{eff}}}/{\varepsilon ^2}}} + \frac{{4{\Sigma ^2}|x{|^2}}}{{{c^2}}}$$

(45) $$b = \frac{{8{\Sigma ^2}|x|\eta }}{c},$$

encode the spatial and spectral tapers, and |x|² = x ² + y ² (without loss of generality). This expression is derived from the Fourier Transform over Gaussians, and then the integral over dldm.ⁱ

When η = 0 and for the 28-m baseline triangles, the expected bispectrum normalisation is

(46) $$\sqrt {P{{(u,v,\eta )}^3}} /V = 2.8 \times {10^{21}}{\kern 1pt} {\kern 1pt} {\rm{m}}{\rm{K}^3}{h^{ - 6}}{\rm{M}}p{c^6}.$$

For the 14-m triangles, we find $\sqrt {P{{(u,v,\eta )}^3}} /V = 3.3 \times {10^{21}}{\kern 1pt} {\kern 1pt} {\rm{m}}{\rm{K}^3}{h^{ - 6}}{\rm{M}}p{c^6}$. When compared with the expected bispectrum value, we find that (28 m):

(47) $$\langle {\cal B}\rangle = 1.7,$$

and $\left\langle \cal B \right\rangle = 0.6$ for the 14-m baselines, which exceed the equilateral triangle configuration estimates from the MWA data. As with the bispectrum estimate, the isosceles configurations have expected power values that fall rapidly with η, and are less comparable to the data in these idealised scenarios. However, for the k ₁ = 0.1h Mpc ⁻¹ stretched configuration, we expect for 28 m (14 m):

(48) $$\langle {\cal B}\rangle = 4.0{\kern 1pt} {\kern 1pt} (240,000).$$

These values are ratios of very small numbers, and therefore are highly dependent on numerical details and are not representative. However, they may lend support to the idea that the normalised bispectrum is difficult to interpret because it relies on foreground details in both the bispectrum and power spectrum.

Alternatively, the combination of power spectrum, dimensional bispectrum and normalised bispectrum may help to shed additional light on whether data are really foreground free. Given the different behaviour of foregrounds in these statistics, this information may be used to discriminate cosmological information from foregrounds, or to help to design some iterative foreground cleaning algorithm, taking into consideration their behaviour in cosmological simulations. In this scenario, the normalised bispectrum may provide useful information.

Figure 8 displays the normalised foreground bispectrum for isosceles configurations in k _⊥ − k _|| space. For the point-source foregrounds, the power spectrum denominator dominates over the expected bispectrum signal, yielding values < 10⁻³ across all of parameter space. This presents an interesting divergence from the usual expectation of foreground bias in the power spectrum, where foregrounds add overall signal. In this case, a large measurement that exceeds the thermal noise is consistent with a cosmological origin, and not with residual foregrounds.

10.2. Foreground bispectrum error

We now turn our attention to consideration of the signal variance due to residual foregrounds, $\langle B_{{\rm{F}}G}^2\rangle $, such that [cf., Equation (14)]:

(49) $$\Delta {B^2} = {{3\sigma _{{\rm{therm}}}^6} \over {\sum\limits_j {W_j}}} + \langle B_{{\rm{F}}G}^2\rangle ,$$

and

(50) $$\langle B_{{\rm{F}}G}^2\rangle = \langle V_1^ * V_2^ * V_3^ * {V_1}{V_2}{V_3}\rangle .$$

This reduces to a relatively simple expression for the simple point-source case, due to the cancelling of complex components (this is not generally true for the covariance). Using the same formalism as earlier, and again considering the Poisson-distributed nature of the flux density of the sources, we find:

(51) \begin{array}{*{20}{l}} {\langle B_{{\rm{F}}G}^2\rangle }& = &{\int {{S^6}} \frac{{dN}}{{dS}}dS\int {{A^6}} (l,m)dldm}\\ {}&{}&{ \times \int {\vec \Upsilon } {{\rm{e}}^{( - 2\pi i({\eta _1}\Delta {\nu _{12}} + {\eta _2}\Delta {\nu _{34}} + {\eta _3}\Delta {\nu _{56}})}}d\vec \nu }\\ {}& = &{\alpha \frac{{S_{{\rm{m}}ax}^{7 + \beta }}}{{7 + \beta }}\left( {12\pi {\rm{ }}\frac{{{c^2}{\varepsilon ^2}}}{{\nu _0^2{A_{{\rm{eff}}}}}}} \right){{\rm{e}}^{ - 4{\pi ^2}{\Sigma ^2}(\eta _1^2 + \eta _2^2 + \eta _3^2)}},}\\ {}& = &{\alpha \frac{{S_{{\rm{m}}ax}^{7 + \beta }}}{{7 + \beta }}\left( {12\pi \frac{{{c^2}{\varepsilon ^2}}}{{\nu _0^2{A_{{\rm{eff}}}}}}} \right){{\rm{e}}^{ - 6{\pi ^2}{\Sigma ^2}\eta _1^2}},} \end{array}

where $\vec \Upsilon \equiv \Upsilon ({\nu _1})\Upsilon ({\nu _2})\Upsilon ({\nu _3})\Upsilon ({\nu _4})\Upsilon ({\nu _5})\Upsilon ({\nu _6})$. This expression is flat in angular modes, and decays rapidly in line-of-sight modes.

Comparing this with the expected value of the foreground bispectrum, Equation (37), we can form the foreground bispectrum signal-to-error ratio (Figure 9). The signal bias exceeds the uncertainty for small scales, but on the larger scales of interest for EoR, the uncertainty dominates. Nonetheless, there is no line-of-sight dependence, demonstrating that the foreground bias and uncertainty both drop rapidly and are negligible for k _|| > 0.12hMpc⁻¹, implying that for larger k _|| scales, point-source foregrounds are not significant in the signal or noise budget.

Figure 8. Point-source-normalised foreground bispectrum signature of isosceles triangle vectors in k _⊥ − k _|| -space (note the stretched logarithmic colour bar). In this model, the expected foreground signal has fallen to below the expected cosmological signal value by k _|| ≥ 0.12.

One can also now compare the foreground uncertainty to the expected thermal noise level. For the EoR1 field data, the measured uncertainty for the direct bispectrum estimator was 6.7 × 10¹² mK³ Mpc⁶. Figure 10 shows this level (green line) compared with the foreground bispectrum error (red line) as a function of line-of-sight scale. (The gridded estimator has slightly lower noise level, but the distinction is not significant when compared to the large gradient of the foreground contribution.) As with the foreground bias, the error induced by residual foregrounds drops steeply beyond k _|| = 0.12hMpc⁻¹, and falls below the thermal noise (even in this case with a small dataset for the EoR1 field).

Figure 9. Ratio of point-source foreground bispectrum bias to uncertainty, for isosceles triangle vectors in k _⊥ − k _|| -space (note linear plot). The bias exceeds the uncertainty at large angular modes, but rapidly falls below for larger k _⊥, with no dependence on line-of-sight scale.

Figure 10. Errors for point-source-normalised foreground bispectrum (red) and thermal noise (green), for isosceles triangle vectors in k _⊥ − k _|| -space and 300 observations used in this work for the EoR1 field.

As a final step to assessing the advantages of the bispectrum compared with the power spectrum to detect the cosmological signal, we divide the expected 21-cm bispectrum values for the faint and bright galaxies, presented in Table 1 by the foreground error for these modes, and compare with that for the power spectrum. The power spectrum values for faint and bright galaxies are taken from the same underlying dataset generated by 21-cm FAST (Mesinger, Furlanetto, & Cen 2011). For all but the k ₁ = 0.1hMpc^-1 mode, the foreground uncertainty is negligible, and the ratio is uninteresting. For k ₁ = 0.1hMpc^-1, k ₂ = k ₃ = 0.2hMpc^-1:

(52) $$\matrix{ {{P_{21}}} \hfill & = \hfill & {1.3 \times {{10}^4}m{K^2}Mp{c^3}} \hfill \cr {\Delta {P_{FG}}} \hfill & = \hfill & {0.2 \times {{10}^{ - 1}}m{K^2}Mp{c^3}} \hfill \cr {{B_{21}}} \hfill & = \hfill & {4.4 \times {{10}^9}m{K^3}Mp{c^6}} \hfill \cr {\Delta {B_{FG}}} \hfill & = \hfill & {4.6 \times {{10}^{10}}m{K^3}Mp{c^6},} \hfill \cr } $$

yielding better performance for the power spectrum, within the foreground dominated region.

However, outside of the foreground ‘wedge’, which exists in both power spectrum and bispectrum space, the data uncertainty is limited by the thermal noise, and here the bispectrum achieves higher signal-to-noise ratio for a set observation time:

(53) $$\matrix{ {{P_{21}}} \hfill & = \hfill & {1.3 \times {{10}^4}m{K^2}Mp{c^3}} \hfill \cr {\Delta {P_{{\rm{therm}}}}} \hfill & = \hfill & {5.7 \times {{10}^6}m{K^2}Mp{c^3}} \hfill \cr {{B_{21}}} \hfill & = \hfill & {4.4 \times {{10}^9}m{K^3}Mp{c^6}} \hfill \cr {\Delta {B_{{\rm{therm}}}}} \hfill & = \hfill & {5.3 \times {{10}^{11}}m{K^3}Mp{c^6}.} \hfill \cr } $$

Taking the ratios we find,

$$\matrix{ {{P_{21}}/\Delta {P_{{\rm{t}}herm}}} \hfill & = \hfill & {0.002} \hfill \cr {{B_{21}}/\Delta {B_{{\rm{t}}herm}}} \hfill & = \hfill & {0.008.} \hfill \cr } $$

Accounting for the fact that the gridded bispectrum averages down with t ^1.5 while the gridded power spectrum averages with t, the observing time multiple (above 10 h) for a detection (SNR = 1) is

$$\matrix{ {{t_P}} \hfill & = \hfill & {500 \times } \hfill \cr {{t_B}} \hfill & = \hfill & {25 \times .} \hfill \cr } $$

Therefore, the bispectrum detection can theoretically be achieved in a fraction of the time of the power spectrum detection, for thermal noise-limited modes close to the EoR wedge. An SNR = 1 detection level could potentially be reached in 250 h, for this wave mode. This conclusion is relevant for the MWA, where the excellent instantaneous uv-coverage allows for rapid observation of triangle configurations. [We note that the power spectrum SNR shown here is not inconsistent with previous expectations for the performance of the MWA, because it applies only to this single mode (Beardsley et al. 2016; Wayth et al. 2018)]. Future work presented in Watkinson & Trott (2019) will explore a more full range of triangle configurations and foreground bias and error.

11. Discussion and conclusions

As discussed, the model used for the foreground bispectrum signal predicts that isosceles configurations have amplitudes that fall rapidly with non-zero k _||. Nevertheless, we find that the numerator and denominator of the normalised bispectrum scale such that its amplitude in this model increases as a function of η. For the 28-m baselines, the ratio doubles by k _|| > 0.014h Mpc ⁻¹. Over the same k range, the dimensional bispectrum is rapidly decaying, losing seven orders of magnitude from the k _|| = 0 mode. In this model, the expected foreground signal has fallen to below the expected cosmological signal value by k _|| ≥ 0.12. The results from the MWA datasets have some modes thermal noise-limited at 10 h, and only the k _|| = 0 mode is clearly foreground dominated for all experiments. In line with the discussion of Bharadwaj & Pandey (2005), it is possible that the dimensional bispectrum is less affected by foregrounds than the power spectrum. However, the normalised bispectrum is more difficult to interpret, given the different observational foreground effects on the bispectrum numerator and power spectra denominator.

Reduction of foreground contamination is an active field of research in 21-cm EoR experiments, and primary motivator for testing statistics other than the power spectrum. Despite the normalised bispectrum providing a cosmologically stable and robust estimate of non-Gaussianity compared with the dimensional bispectrum, the expected foreground value is difficult to discriminate from the expected signal value (Watkinson et al. 2018). Conversely, the dimensional bispectrum yields values that are highly significant detections, showing clear foreground contamination. Thus, the normalised bispectrum may not be the best discriminant in real EoR experiments. For future experiments, with higher sensitivity, exploration of modes with negative bispectra may help discrimination from foreground contamination, where the bispectrum is expected to be positive This is explored further in Watkinson & Trott (2019) and demonstrated previously in Lewis (2011). It would also be interesting to study the signature of calibration errors in radio data on bispectrum estimates, to explore whether they have an imprint that can be discriminated from the cosmological signal.

The thermal noise levels, as are achieved in these 10 h datasets for large k isosceles configurations, are 3–4 orders of magnitude larger than the expected bispectrum value for these configurations at low redshifts (Majumdar et al. 2018). The gridded bispectrum noise scales with observation time to the power of 1.5, requiring a 1 000 h observation with the MWA to achieve a cosmological detection. This estimate is in line with predictions from Yoshiura et al. (2015). Further advantage may be gained from incoherent addition of isosceles triangle configurations with similar vector lengths, where the bispectrum is expected to vary slowly with changing parameters. An initial test of this for the k ₁ = 0.1h Mpc⁻¹ mode shows an improvement in sensitivity by a factor of ten for the gridded estimator, yielding a theoretical detection of the signal with 150 h of data. The direct estimator scales incoherently with time (t ^0.75), due to the incoherent addition of triangles from different observations, but does utilise coherent addition of instantaneously-redundant triads.

We have presented the first effort to estimate the cosmological bispectrum from the EoR with 21 h of MWA data, and have shown the parts of parameter space that are consistent with thermal noise at this level, using two types of bispectrum estimator. These two approaches are presented in order to demonstrate the practicalities of estimation of the bispectrum with real radio interferometer data. We have also derived a form for the expected bispectrum signature of point-source foregrounds for equilateral and isosceles configurations, and demonstrated broad consistency between the analytic model and the estimates obtained from the data.

By considering the foreground bispectrum variance in the noise estimation, we have demonstrated that both the foreground bias and variance are insignificant for k _|| < 0.12h Mpc⁻¹, allowing these regions of parameter space to be probed with dominant thermal noise. Due to the ability of the gridded bispectrum estimator to reduce thermal noise proportional to t ^1.5, unlike the power spectrum which reduces with t, the 21-cm cosmological bispectrum may be detectable with fewer observing hours than the power spectrum for arrays with excellent instantaneous uv- coverage (i.e., with well-sampled baselines). This insight makes observational pursuit of the bispectrum worthwhile for some current instruments.

In a companion paper (Watkinson & Trott 2019), we explore optimal triangles to study from a signal and foreground contamination ratio perspective. Future work can also study the signature of calibration errors on the bispectrum. This work helps to define the optimal observational strategy and approach to bispectrum studies.