Deep multi-redshift limits on Epoch of Reionisation 21 cm Power Spectra from Four Seasons of Murchison Widefield Array Observations

The spatial power spectrum quantifies the signal power as a function of spatial scale, $k$ ($h$Mpc^-1). It is the Fourier Transform of the two-point spatial correlation function, and can be estimated from the volume-normalised Fourier-Transformed brightness temperature field:

In radio interferometry, the angular scales are related to the Fourier modes of the measured interferometric visibilities, and the line-of-sight modes can be mapped with spectral channels (for a resonant line signal): $(u,v)\rightarrow k_{\bot}$, $\mathcal{F}(f)=\eta\rightarrow k_{\parallel}$. As such, we extract angular modes directly from the measured visibilities, such that:

where $A$ is the angular dimension and $\Omega$ is the observing volume. The brightness temperature and source flux density are related via the equation for the specific intensity, which is linear in the radio regime:

where $k_{B}$ is the Boltzmann constant. During the Fourier Transform, the area dimension is collected to yield power spectral flux density (Jansky).

Algorithmically, the power spectrum is formed from a data cube with dimensions $(u,v,f)$, where the angular Fourier Transform to ($u,v$) is performed natively by the interferometer, and the spectral channels are used to map the line-of-sight. Data are observed over time and integrated together coherently by gridding measured visibilities onto a common discretized $uv$-plane. The final steps are to Fourier Transform along frequency in each cell, and then to square to arrive at the unnormalised power. A power spectrum formed in this way may be used for cosmological measurements, because the three $k$-vectors are orthogonal and can easily be mapped to spherical $k$. In this work, the 2D (cylindrically-averaged) and 1D (spherically-averaged) power spectrum are used. The former principally acts to identify foreground leakage into the parameter space used for EoR analysis. The latter provides the cosmologically-relevant measurements. In 2D, the angular and line-of-sight modes are separated, and denoted $k_{\bot},k_{\parallel}$. In 1D, these are averaged to $k^{2}=k_{\bot}^{2}+k_{\parallel}^{2}$. We present the dimensionless power spectrum in 1D:

The alternative approach to approximating the power spectrum is the delay transform (Parsons et al., 2012; Jacobs et al., 2015). Here, each individual baseline’s data are Fourier Transformed along its frequency channels, and the power formed through their squared quantities. This approach is straight forward, but can be difficult to interpret cosmologically, because the $k$-vectors are non-orthogonal except on short baselines (the line-of-sight transform evolves with frequency) and there is no angular correlation encoded between baselines with similar lengths and orientations. Nonetheless, the delay transform is very useful for quality assurance to ensure the visibility data are not corrupted or contaminated in the regions of parameter space used for EoR analysis (i.e., the ‘EoR Window’, a region of $k_{\bot},k_{\parallel}$ parameter space outside of the region expected to be dominated by smooth-spectrum foregrounds).

CHIPS – the Cosmological HI Power Spectrum estimator (Trott et al., 2016) – is one of the signal processing pipelines used by the MWA EoR collaboration for taking calibrated data and processing them to output power spectra, with associated uncertainties. In its original form, it was intended to undertake a full thermal noise plus residual foreground signal inverse covariance weighting, to optimally extract cosmological information. There are difficulties with this approach, and these were explored in the literature at the time, and have more recently been demonstrated in the retracted results from the PAPER collaboration (Cheng et al., 2018). Inaccurate residual foreground models, and failure to fully and independently understand their internal covariance and covariance with the signal, can easily lead to signal loss. As such, CHIPS is used primarily (and entirely in this work), as an inverse variance estimator, where the baseline weighting is used for sampling.

Line-of-sight modes for computing the power spectrum are extracted from spectral sampling of the data. The MWA’s signal processing chain contains filterbanks that yield 24 coarse channels of 1.28 MHz over the full 30.72 MHz band. Within these coarse channels, the native spectral resolution is 10 kHz, but EoR data are observed at 40 kHz resolution and processed at 80 kHz resolution. The shape of the fine polyphase filterbank leads to poor bandpass characteristics at the coarse channel edges, and as such, a single 80 kHz channel is flagged at each end of each coarse channel. This yields regularly-spaced missing channels in the final output visibilities. A Fourier Transform over the data to retrieve the line-of-sight spatial scales will contain a comb shape due to these missing channels, where the $k_{\parallel}=0$ mode is copied in harmonics of the coarse channel separation. There are several ways to handle this; in this implementation of CHIPS, we use an ordinary kriging (a Gaussian Process Regression) to provide an interpolate of these data that uses the covariance structure of the data (Rasmussen & Nickisch, 2010; Wackernagel, 2003). Kriging has been used to fit for foregrounds in LOFAR EoR datasets, using an optimised set of hyperparameters (Mertens et al., 2018). The kriging kernel (variogram) is estimated conservatively to contain a noise-like variance and a frequency covariance which decays smoothly across several megahertz. Kriging applies an interpolation over unsampled data using linear weights of sampled data. The weights are computed to minimise the mean-squared error and yield an unbiased interpolate, subject to a specified data covariance matrix. In this work, we assume a covariance matrix with a nugget (variance term) provided by the thermal noise variance, and a Gaussian-shaped spectral covariance for the foregrounds with a characteristic length of 50 channels. The functional form of the kriging kernel is kept consistent across all datasets to avoid biasing results with fine-tuning, and is given by:

where the relative scaling of the nugget to the squared-exponential is appropriate for the relative amplitude of the thermal noise and foregrounds. The same kriging kernel is applied to all datasets, and is applied to the real and imaginary components of each $uv$-cell along the frequency distribution in the gridded $(u,v,\nu)$ data cube. The results are not strongly-dependent on the choice of spectral scale, with statistically-similar results occurring for values of 2–5 MHz. Testing demonstrates that the application of kriging does not bias results in the modes used for measurement (see Section 2.4), but can at higher-order modes. It was initially implemented to access more modes close to the coarse channel harmonics, but this has provided only limited improvement for some datasets. Nonetheless, it does yield improved results in those modes used for limits in this work ($k<0.4$ $h$Mpc^-1) compared with omitting the kriging. In future work, we will either (1) not apply kriging at all; (2) invest more effort to understand and refine it so that we are confident that is unbiased across the full range of $k$-modes. In this work, we retain it, because it is infeasible to re-process all of these data and we report results at unbiased wave modes only.

Data were observed with the Murchison Widefield Array, a general-purpose low-frequency radio interferometer operating at the Murchison Radioastronomy Observatory in Western Australia (Tingay et al., 2013; Bowman et al., 2013). Phase I of the array comprised a pseudo-random core for EoR science, surrounded by sparser remote tiles for angular resolution. In Phase II of the array (Wayth et al., 2018; Beardsley et al., 2019), the telescope consists of 256 tiles of 16 dual-polarisation dipole antennas in a regular 4$\times$4 grid, spread over $\sim$5 km. Only 128 tiles can be connected to the signal chain at one time, and the telescope operates in an "Extended" (survey science; long baselines) or "Compact" (high surface brightness sensitivity; short baselines) configuration. The EoR experiment uses the latter configuration. The sky model for instrument calibration and source subtraction are formed from the Phase I configuration, and auxiliary data from other telescopes.

The EoR experiment observes in three bands (Jacobs et al., 2016): ultralow-band (75–100 MHz), low-band (139–167 MHz), and high-band (167–197 MHz), with more than ninety percent of data observed in low- and high-bands. The data for this work are only taken from the upper two bands (139–197 MHz), consistent with the reionisation epoch when the signal is expected to be in emission (as compared with the Cosmic Dawn). The primary EoR experiment observes data from three observing fields, chosen to minimise sky temperature (away from the Galactic plane) and containing bright calibration sources. They are EoR0 (RA=0h, Dec=$-$27^o), EoR1 (RA=4h, Dec=$-$27^o), and EoR2 (RA=10.3h, Dec=$-$10^o). EoR1 contains Fornax A, an extended ($\sim$1 degree) double-lobed radio source part way down the main primary beam, and EoR2 contains Hydra A. These two sources both help and hinder data calibration and need to be subtracted with high precision for EoR science (Procopio et al., 2017; Trott & Wayth, 2016). EoR0 contains the setting Galactic plane in early (pre-zenith) pointings, yielding power from the horizon in power spectra. Of the data presented in this paper, 50% are EoR0, 28% are EoR1, and 22% are EoR2. This is a function of the data that have been calibrated, and not a reflection on the overall contributions of each. However, in general EoR2 is observable at the end of the season, and there are fewer hours available for it than for EoR0 and EoR1.

The data used in this work span Phase I and II of the array. Despite redundancy being available in Phase II, and used in other MWA pipelines such as Fast Holographic Deconvolution (Barry et al., 2019a, FHD,) plus Omnical (Li et al., 2018, 2019; Zhang et al., 2020), the RTS calibration currently only performs sky-based calibration. Given the interest in understanding the utility of hybrid arrays with redundant and non-redundant baselines, it would be useful to compare the results for the different configurations. Unfortunately, the data in this work are 92% Phase I and 8% Phase II, with no individual Phase II set exceeding 5 hours. As such, comparison of datasets suffers from the small and concentrated number of observations, and is not useful.

Data used in this work were observed over five pointings, separated by $\sim$6.8 degrees on the sky (27 minutes per pointing), corresponding to the beamformer analogue delay settings that produce a consistent primary beam response. Of all of the pointings, the five central pointings (including zenith) are found to have the least contaminated power (Beardsley et al., 2016), and have the most well-behaved beam patterns, and are used exclusively in this work. Off-zenith pointings are only pursued for EoR0 and EoR2, because the zenith results from these fields are sufficiently-interesting to warrant coherently combining data from differing pointings. Table 1 lists the Alt/Az and name for each pointing used here.

The data for this work were observed in the 2013B, 2014B, 2015B, 2016B observing seasons, with a range of August 2013 – January 2017, spanning Phase I and Phase II configurations.

The data were selected by extracting all of the observations available in the MWA database (hosted by the Pawsey Supercomputing Centre) from 2013 to 2017, which had been successfully calibrated with the RTS (output files were produced with finite values, and bandpass and phase plots looked reasonable), had complete calibrated visibility files with the standard temporal (8 second) and spectral (80 kHz) resolution, and had an ionospheric activity metric value associated with them. In all cases, the most recently processed calibration was used, with processing dates ranging from 2017–2019. The same RTS version was used for all processing. Most data that satisfy the requirements are from 2013–2017. This search yielded 13,591 2-minute observations.

There are four main quality assurance metrics applied to refine the dataset to be selected for power spectrum analysis: (1) calibration success (no errors; all frequencies present with finite-valued data), (2) ionospheric activity, (3) delay-space EoR Window Power, and (4) delay-space EoR Wedge Power. Here we describe these metrics and the order in which they are applied.

1. 1.

   Calibration: Data are calibrated using the MWA Real-Time System (RTS, Mitchell et al., 2008; Jacobs et al., 2016; Trott et al., 2016); the RTS performs direction-independent and direction-dependent (DD) calibration. It uses the primary beam model from Sokolowski et al. (2017) to select the 1000 apparent brightest sources for each snapshot pointing for calibration. The brightest five are used for DD calibration with a full Jones matrix solution for each source. These DD solutions are applied to the full set of 1000, effectively peeling five and directly subtracting the remaining 995. These calibrated data are fed through a validation process that confirms existence of all spectral channels with finite-valued data;

2. 2.

   Ionospheric Activity: we use the ionospheric activity metric developed by Jordan et al. (2017), which uses the measured versus expected source positions of 1000 point sources in the field-of-view to estimate an ionospheric phase screen, and derive an activity metric that combines median source offset with source-offset anisotropy. Different thresholds are set for Phase I and Phase II datasets, where the reduced angular resolution of the latter array configuration yields a higher base activity level;

3. 3.

   Window Power: we use the delay transform power spectrum estimator to compute the power in the EoR Window, below the MWA’s first coarse channel harmonic. Power is computed incoherently across baselines. Window Power is computed in the region bounded by the main beam lobe and the first coarse channel harmonic, for baselines with length $<100\lambda$ ($0<|\vec{u}|<100\lambda$, $k_{\bot}<k_{\parallel}<0.4$);

4. 4.

   Wedge Power: we use the delay transform power spectrum estimator to compute the power below the EoR Window, in the primary beam main lobe wedge. Power is computed incoherently across baselines. The Wedge Power is computed for a region bounded by $k_{\parallel}=0$ and the main beam lobe wedge, also for baselines with length $<100\lambda$ ($0<|\vec{u}|<100\lambda$, $k_{\bot}>k_{\parallel}>0$). Both are normalised by the number of contributing cells to yield an average power per cell (Jy²).

Using these metrics, cuts are made as follows, with the intent that unusually high-valued, or unusually low-valued snapshots are omitted. We take the full dataset and compute the average and standard deviation wedge and window power metrics for each data-type (Phase I or II, and high- and low-band). We omit snapshots that have wedge and window power values that fall outside of the primary mode of the distribution of values, and also snapshots that show consistent wedge power but high window power. The ionospheric cuts are made based on the metric values that are expected to produce biased results according to the analysis of Trott et al. (2017). After the cuts are applied, datasets are then ordered by ionospheric metric value only. The distributions of Window Power are not Gaussian, and are generally multi-modal. Figure 1 displays the distribution for EoR0 high-band, with blue denoting included observations, and red denoting omitted. There is a clear separation of the two clusters, suggesting that calibration errors, and not statistical fluctuations, are the primary cause. This distinct behaviour is also observed for the other datasets. The data cuts are conservative insofar as all observations in the primary mode are included for all datasets.

Table 2 describes the components of each dataset, and the changes after each stage of assessment.

Figures 2 – 5 show examples of power in the EoR window versus ionospheric activity for different fields and frequencies. The omitted observations plotted in these figures are all from Phase I; Phase II observations extracted from the database were consistent with a quiet ionosphere.

Blue and green filled circles denote observations that meet the assessment criteria for Phase I and II, respectively. Open red circles are observations that are omitted due to high Window power or high ionospheric activity. These cuts remove 10–20 percent of the data. The majority of the removed data contain particular tiles or baselines with very poor calibration, leading to excess power in those modes.

Correctly calculating the normalisation from measurement units (Jansky and Hertz) to cosmological units (megaparsecs) seems trivial, but is complicated by the choices made during the analysis pipeline (e.g., cross-multiplying even and odd samples for the cross power spectrum, and the definition of Stokes I with respect to polarisation axes). It is important to ensure that the normalisation is correct, and that signal is not being lost due to coherence of coherently-gridded data. The former can be ensured via matching of different approaches to calculating the noise, and to internal consistency between independent MWA analysis pipelines (Barry et al., 2019a). The latter can be achieved using a signal simulation; ensuring that the input power spectrum is recovered after passing through the CHIPS pipeline.

We performed a large simulation using a modified version of 21cmFAST (Mesinger et al., 2011; Greig & others, 2020), tuned to the MWA’s large primary beam size and frequency resolution. A box with 6,400 pixels on each side was produced, with 7,500 cMpc on angular sides, and 1.172 cMpc line-of-sight resolution. This large simulation is important for creating simulations that can be directly applied to MWA, without need for padding or interpolation. The data are produced as a light-cone, with signal evolution as a function of redshift. We apply a beam model to the slices, perform an angular Fourier Transform, and convert the brightness temperature units to Jansky per steradian.

Figure 6 shows the input EoR 21cm power spectrum over 384 channels centred at $z=6.8$ (green) and the signal power recovered through CHIPS (red). This demonstrates good consistency of the input and output signal levels in both shape and amplitude.

Given that we do not perform any post-calibration subtraction of signal from the data, these results provide confidence that signal loss is not occurring in the CHIPS pipeline.

In addition to the regular simulation, we also perform the same operation but with the missing channels corresponding to those in the actual data. We apply the kriging to the mitigate the missing channels and check that the output power levels are still consistent. This procedure demonstrates that (1) the kriging is not biasing the results, and (2) that kriging is offering no benefit at small $k$, but does close to the coarse channel harmonics. Note that we do not use those harmonic modes in our measurements. Figure 7 shows the 1D power spectra from the same simulation with the kriging applied (blue) and the full dataset (red). The contaminated modes, where there is a discrepancy, correspond to the location of the coarse channel harmonic. The primary results shown in this work are for $k<0.4~{}h$Mpc^-1.

We also compare simulation outputs through RTS-CHIPS with those through FHD-$\epsilon$ppsilon (Barry et al., 2019b), to ensure power and noise-level consistency. A similar EoR simulation was produced and passed through both pipelines, yielding consistent results (N. Barry, in prep.). Although this deviates somewhat from previous MWA EoR papers where the actual data sets were processed through both pipelines (Barry et al., 2019b; Li et al., 2019; Beardsley et al., 2016, Figure 7,), this retains the same philosophy of ensuring robustness of results using independent calibration and analysis methods.

We present a census of the data from four years of the MWA EoR experiment; the cleanest subsets of data assessed in this work, taken from three observing fields, two broad observing spectral bands, and multiple telescope pointings (for EoR0 and EoR2).

Ultimately, the aim is to combine data from individual pointings coherently. The data have been processed using the same phase centre, with a common framework, and with an optimal weighting to allow for coherent addition of data (i.e., the weights are carried through the analysis such that separate datasets can be added using an optimal weighting). However, the primary beam response of the telescope changes appreciably between pointings, and there is scope for decoherence if the telescope response is not correctly modelled, which would lead to undesirable signal loss.

We can test for this decoherence by comparing the regions of parameter space that should be consistent, e.g., the foreground-dominated EoR Wedge region at low $k_{\parallel}$. Because this is sky power, it should be retained between pointings and upon coherent addition of pointings. It will not be identical; the different sky response will mean that the power is in different locations, but it should retain an overall power level. Figure 20 plots the $k_{\parallel}=0$ mode power for the zenith (red), off-zenith (green) and combined (black) datasets for EoR0High at $z=6.5$. The power is retained during coherent addition. Figure 21 then displays the ratio (left) and fractional difference (right) of 2D power for the zenith pointing and combined pointings. A ratio consistent with unity and small fractional difference ($<10\%$) in the foreground wedge (low $k_{\parallel}$) demonstrates that combining pointings coherently is reasonable, because power is not being lost.

Having established that coherent addition of data from the same observing field and frequency range, is possible, we combine the best zenith-pointed observations with the four off-zenith pointings for EoR0 high- and low-bands, and the EoR2 high- and low-bands. Note that we do not present off-zenith pointings for the EoR1 field, due to the poor results from the zenith data. Figures 22 and 23 display EoR0 results, and Figures 24 and 25 display EoR2 results.

The best results at each redshift are reproduced in Table 4.

The results of combining data from different pointings for EoR0 and EoR2 demonstrate better performance in EoR0 at low redshift and EoR2 at high redshift. Given that the distributions of ionospheric activity and EoR Window Power are comparable between the fields, this is likely due to the different Galactic and extended structures drifting through the primary beam sidelobes as a function of frequency. The MWA primary beam introduces strong spectral gradients in the beam nulls, amplifying any effect of mis-modelling of the sources or primary beam in this regions, and potentially imprinting strong spectral structure on residual signals.

In EoR0, the Galaxy presents more prominently at low frequency due to the larger beam size, whereas the Puppis A, Centaurus A and Centaurus B sources in the EoR2 sidelobes may be better placed with respect to large spectral gradients in the primary beam. Figure 26 demonstrates this in the high-band 2D power spectrum, showing the ratio and difference of the power in the EoR0 to EoR2 fields. Aside from an overall decrement in the power in EoR0 (red), there is a power increase along the horizon modes in EoR0, indicative of the effect of the spectrally-smooth Galactic Plane. Conversely, the lack of any extended models for the spectrally complex Centaurus A and Puppis A sources leads to additional leaked power in the EoR Window on large scales.

We can study the instrumental response to the sky for the pointings where the results are best: EoR0 in high-band, and EoR2 in low-band. In Figures 27 and 28 we overlay the sky with the spectral gradient of the primary beam response, for the zenith pointing of EoR0 and EoR2, respectively (Cook & Seymour, 2020). In the right-hand plots, the beam gradient is shown separately for clarity.

In EoR0, the Galactic Plane is in the spectrally-flat horizon region, and the bright source 3C444 (80 Jy) is away from regions of large index. No bright extended sources reside in regions of large spectral gradient.

In EoR2, the Galactic Plane, with a large number of supernova remnants (e.g., Puppis A, $>$200 Jy), and extended radio galaxies (e.g., Centaurus A, $>$1000 Jy, and Centaurus B, $>$100 Jy) contains most of its structure in spectrally-flat regions, with high spectral indices mostly avoided by these complex sources.

For both of these fields, one can see that small changes in the location of the spectrally-steep beam nulls (as occurs when changing frequency bands) could lead to complex sources incurring large instrumental indices. Given that these bright sources near the field edges are rarely well-subtracted in the current sky models (if a good model even exists), there is potential for leakage into the EoR Window. Dead dipoles (i.e., those where the dipole is present but not delivering signal) within a tile will tend to smooth out these nulls, leading to additional complexity in the sampled signal. The impact of this for these fields is left for future work.

To study this, we plot the overlaid images for the EoR0 Low-band and EoR2 High-band (zenith) in Figure 29. In EoR0, 3C444 (50 Jy) resides in an area of large spectral index. As one of the brightest sources in the field, this may be having an impact if not adequately modelled. In EoR2, large areas of the Galactic Plane overlay regions of large spectral index. Without more in-depth modelling, which we leave to the next stage of this work, it is difficult to ascertain the direct impact of this complicated field.