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. Reference Tingay2013; Bowman et al. Reference Bowman, Cairns, Kaplan, Murphy and Oberoi2013). 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. Reference Wayth2018). 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. Reference Li2018, Joseph et al. Reference Joseph, Trott and Wayth2018). 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 uncorrelatedFootnote ^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 approximationFootnote ^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. Reference Datta, Bowman and Carilli2010):

(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. Reference Parsons, Pober, McQuinn, Jacobs and Aguirre2012; Ali et al. Reference Ali2015; Thyagarajan et al. Reference Thyagarajan2015), 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. Reference Morales, Hazelton, Sullivan and Beardsley2012, 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. (Reference Watkinson, Giri, Ross, Dixon, Iliev, Mellama and Pritchard2018) 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 Reference Eggemeier and Smith2017; Brillinger & Rosenblatt Reference Brillinger, Rosenblatt and Bernard1967). 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.

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.Footnote ^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)Footnote ^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. Reference Watkinson, Giri, Ross, Dixon, Iliev, Mellama and Pritchard2018) 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. Reference Jacobs2016). These 629 calibrated and peeled observations were used for bispectrum estimation.

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. Reference Intema, van Weeren, Röttgering and Lal2011; Franzen et al. Reference Franzen2016):

(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 (Reference Murray, Trott and Jordan2017). 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 (Reference Watkinson and Trott2019) 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. Reference Majumdar, Pritchard, Mondal, Watkinson, Bharadwaj and Mellema2018).

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.Footnote ⁱ

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 Reference Mesinger, Furlanetto and Cen2011). 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. Reference Beardsley2016; Wayth et al. Reference Wayth2018)]. Future work presented in Watkinson & Trott (Reference Watkinson and Trott2019) 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 (Reference Bharadwaj and Pandey2005), 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. Reference Watkinson, Giri, Ross, Dixon, Iliev, Mellama and Pritchard2018). 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 (Reference Watkinson and Trott2019) and demonstrated previously in Lewis (Reference Lewis2011). 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. Reference Majumdar, Pritchard, Mondal, Watkinson, Bharadwaj and Mellema2018). 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. (Reference Yoshiura, Shimabukuro, Takahashi, Momose, Nakanishi and Imai2015). 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 Reference Watkinson and Trott2019), 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.