Epoch of reionization power spectrum limits from Murchison Widefield Array data targeted at EoR1 field

ABSTRACT

Current attempts to measure the 21 cm power spectrum of neutral hydrogen during the epoch of reionization (EoR) are limited by systematics that produce measured upper limits above both the thermal noise and the expected cosmological signal. These systematics arise from a combination of observational, instrumental, and analysis effects. In order to further understand and mitigate these effects, it is instructive to explore different aspects of existing data sets. One such aspect is the choice of observing field. To date, MWA EoR observations have largely focused on the EoR0 field. In this work, we present a new detailed analysis of the EoR1 field. The EoR1 field is one of the coldest regions of the southern radio sky, but contains the very bright radio galaxy Fornax-A. The presence of this bright extended source in the primary beam of the interferometer makes the calibration and analysis of EoR1 particularly challenging. We demonstrate the effectiveness of a recently developed shapelet model of Fornax-A in improving the results from this field. We also describe and apply a series of data quality metrics that identify and remove systematically contaminated data. With substantially improved source models, upgraded analysis algorithms and enhanced data quality metrics, we determine EoR power spectrum upper limits based on analysis of the best ∼14-h data observed during 2015 and 2014 at redshifts 6.5, 6.8, and 7.1, with the lowest 2σ upper limit at z = 6.5 of Δ² ≤ (73.78 mK)² at k = 0.13 h Mpc⁻¹, improving on previous EoR1 measurement results.

methods: data analysis, methods: observational, techniques: interferometric, (cosmology:) dark ages, reionization, first stars

1 INTRODUCTION

Over cosmic time, hydrogen has gone through three major transitions: ionized hydrogen post-nucleosynthesis in the primordial plasma, neutral hydrogen as the universe cooled by expansion, and finally reionized hydrogen as a result of the formation of the first stars and galaxies (Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012). Our goal is to measure the H|$\rm {I}$| spin temperature and neutral fraction via its 21-cm signal to follow the evolution of neutral gas in the IGM during the epoch of reionization (EoR).

The GMRT (Giant Metrewave Radio Telescope; Paciga et al. 2013), LOFAR (LOw Frequency Array¹; Yatawatta et al. 2013), PAPER (Precision Array for Probing the Epoch of Reionization²; Parsons et al. 2010), LWA (Long Wavelength Array³; Ellingson et al. 2009), and the MWA (Murchison Widefield Array⁴; Bowman et al. 2013; Tingay et al. 2013) are first generation instruments with the primary goal of detecting the EoR signal.

In reality, detection of such a faint signal (⁠|${\sim}10\, \mathrm{mK}$|⁠) is complicated; i.e. foregrounds, which are radio emission from galactic and extragalactic sources, are up to five orders of magnitude brighter than the EoR signal (Jelić et al. 2008). They inherently have smooth spectra as opposed to the cosmological 21-cm signal which varies rapidly in frequency. Therefore, by developing suitable analysis approaches, they can be separated from the EoR signal in power space (Dillon et al. 2014, 2015a; Trott et al. 2016). The other major challenge is the systematic errors in calibration due to instrumental effects such as beam models (Joseph et al. 2020) or spectral behaviour of electric components (Trott & Wayth 2016), imperfect sky models (Barry et al. 2016; Byrne et al. 2019; Zhang et al. 2020), radio frequency interference (RFI; Offringa et al. 2015), and ionospheric distortions (Jordan et al. 2017; de Gasperin et al. 2018; Trott et al. 2018). In addition, due to low SNR, we need to integrate over long observation times to average down the noise. Over the last years, many efforts have been taken to develop novel analysis methodologies to overcome these challenges and detect the signal (Morales et al. 2019).

Due to the low sensitivity of first generation telescopes, the only effective measurement of the H i signal is via statistical measurements like the power spectrum. Therefore, we quantify the statistics of the brightness temperature contrast between the 21-cm signal and the CMB and take advantage of the isotropy and homogeneity of the universe to average the 3D k-space power spectrum over spherical shells of constant spatial scale |$\vert \boldsymbol{k} \vert \mathrm{Mpc^{-1}}$|⁠. The current experiments placed upper limits on the EoR power spectrum, such as MWA (Dillon et al. 2014, 2015b; Beardsley et al. 2016; Barry et al. 2019; Li et al. 2019; Trott et al. 2020), GMRT(Paciga et al. 2013), LOFAR (Patil et al. 2017; Mertens et al. 2020), and PAPER (Kolopanis et al. 2019). Due to high contamination by systematics, they have not detected the signal yet. However, in addition to learning invaluable lessons in design and analysis, the results are getting informative enough to constrain the astrophysical models (Greig et al. 2020, 2021).

The next generation of instruments, such as SKA-Low (Square Kilometer Array; Koopmans et al. 2014), HERA (Hydrogen Epoch of Reionization Array; DeBoer et al. 2016), and NenuFAR (New Extension in Nançy Upgrading LOFAR; Zarka et al. 2015), are under design/designed with substantially larger collecting area and hence higher sensitivity with respect to their first-generation counterparts. Therefore, beyond detection on limited redshift or spatial scales, they will be able to characterize the EoR power spectrum with more details and also potentially map the ionized bubbles of hydrogen during the EoR.

Efforts to calculate upper limits for the cosmological H i power spectra, using first generation of instruments, with on-going efforts to improve the precision and quality of the algorithms and their implementation (Paciga et al. 2013; Dillon et al. 2015b; Beardsley et al. 2016; Patil et al. 2017; Barry et al. 2019; Li et al. 2019), have produced the most recent 2σ upper limit of |$1.8 \times 10^3\, \mathrm{mK^2}$| at |$k=0.14\, \mathrm{h\, Mpc^{-1}}$| and z = 6.5 from 110 h of MWA field EoR0 high-band data (Trott et al. 2020).

The very bright (⁠|${\sim}100s\, \mathrm{Jy}$|⁠) and extended radio galaxy Fornax-A is located in the EoR1 field. Therefore, this field requires a more detailed sky model than the EoR0 field. Trott et al. (2020) measured the limits on the power spectrum of deep observations over three MWA EoR fields, and found that EoR1 produced higher limits than the other two fields at each studied redshift. In this work, ∼14 h of Phase-I MWA data, targeting the EoR1 field and observed during 2014 and 2015, is analysed using the RTS/CHIPS [Real Time System, Mitchell et al. (2008), Ord et al. (2010). The Cosmological H|$\rm{\small I}$| Power Spectrum, Trott et al. (2016)] pipeline to measure the upper limits at three redshifts: 6.5, 6.8, and 7.1. The calibration procedure including the source models were improved relative to the previous work (Trott et al. 2020) and customized to obtain more precise calibration solutions. In addition, the application of well-established data-quality metrics such as window power (Beardsley et al. 2016), removed contaminated data from the analysis data set. By applying new metrics, new systematic signatures in data were examined. This provided a multilayered understanding of data anomalies and hence allowed the development of well-targeted approaches for systematic mitigation, either by resolving shortcomings in the analysis or by more effective data selection methods. The measured upper limits over 14 h of data improve on the recent limits from EoR1 published at Trott et al. (2020).

The rest of this paper is organized as follows. In Section 2, we describe the instrument and the observational data we processed in this study and our initial understanding of the EoR1 field. In Section 3, we introduce the analysis pipeline and sky model that were applied to the data. Section 4 compares the data sets observed in 2015/2014, targeted at the EoR1 field, with the corresponding data observed in 2013 or targeted at the EoR0 field. We also compare data variations across different nights and pointings, and finally present the results. Lastly, in Section 5 the analysis results and future possible analysis improvements are discussed.

2 INSTRUMENT AND OBSERVATIONS

2.1 The Murchison Widefield Array

The Murchison Widefield Array is a radio interferometer located at the Murchison Radio Observatory, a radio-quiet site in Western Australia. While supporting several science programmes, it was primarily built with the goal of detecting the EoR signal. The technical design of MWA and its science capabilities are discussed in Tingay et al. (2013) and Bowman et al. (2013), respectively. MWA has been later upgraded in several phases (Wayth et al. 2018; Beardsley et al. 2019).

The MWA Phase-I configuration, which was used to obtain the data in this work, consisted of 128 antenna-tiles each comprising 16 dual polarization dipoles placed on a regular grid. These tiles were distributed within a radius of 1.5 km, with a dense core of 50 m for the purpose of EoR experiments. Each dipole is sensitive to the entire sky and is optimized to operate in the |$80{-}300\, \mathrm{MHz}$| frequency band, out of which a processing bandwidth of |$30.72\, \mathrm{MHz}$| is selected. The radio signals from the dipoles of a tile are combined in an analogue beamformer. For EoR experiments, the beamformer steers the targeted direction of observation in steps of 6.8^○, each taking 27 min. The primary beam for each direction, or ‘pointing’, provides a wide field-of-view of about 25^○ at |$150\, \mathrm{MHz}$|⁠. In this work, we analyse the five central pointings, including two pointings before and after the zenith transit, inclusive. They are labelled in Local Sidereal Time (LST) order as Minus2, Minus1, Zenith, Plus1, and Plus2, respectively. The radio signal is then transmitted to receivers in which they are digitized and filtered into 24 coarse channels of |$1.28\, \mathrm{MHz}$| width, the edges of which need to be flagged prior to our analysis to avoid the aliasing effect introduced by channelization. Then they are passed to the correlator. The correlator then cross-multiplies the signals between all pairs of antennas at 10 kHz resolution. These cross-correlations are called visibilities, which are averaged over frequency and time and then transmitted to be archived. For the data included here, the correlation output resolution is 2s in time and 40 kHz in frequency. The averaged data are written as snapshots of cadence of 112 s, named as one single observation.

2.2 EoR observing programme

MWA EoR observations are collected over three frequency bands: high (⁠|$167{-}197\, \mathrm{MHz}$|⁠) centred at redshift z ≈ 6.8; low (⁠|$139{-}167\, \mathrm{MHz}$|⁠) centred at z ≈ 8.2, and ultra-low (⁠|$75{-}100\, \mathrm{MHz}$|⁠) centred at z ≈ 16. There are three targeted fields for EoR observations that are relatively devoid of galactic emission and extragalactic bright sources, named EoR0, EoR1, and EoR2, as shown in Fig. 1(a). Their coordinates are presented in Table 1.

Figure 1.

(a) MWA EoR fields, (b) EoR1 field, and (c) Fornax-A image.

Table 1.

MWA EoR fields. EoR0 and EoR1 transit over zenith at the MWA latitude.

	Centred at
Field	RA	Dec.
EoR0	0h	−27○
EoR1	4h	−27○
EoR2	10.33h	−10○

The physical characteristics of the field of interest for this work, EoR1, and its comparison to the widely studied one, EoR0, are discussed in Section 4.2. With MWA phase-I configuration, thousands of hours of EoR data were collected from mid-2013 to mid-2016.

2.3 Data

This work aims to analyse and measure the power spectrum limit of the EoR signal from the EoR1 field, observed in High band by the MWA Phase-I configuration. Trott et al. (2020) compared the power spectra from integrations of 300, 500, 680, and 820 EoR1 observations. They showed that adding data beyond 300 observations (until ∼800 observations) did not improve the power spectrum, due to the saturation by systematics (for more discussion, refer to Section 3). Therefore, partially motivated by finite computing resources, we aimed to demonstrate improvements on a final data set of 400–500 observations. Initially, we calibrated a pilot data set spanning all of the observation nights of 2015 (September to December). Based on the results of that processing, we excluded from further analysis nights with a relatively high rate of calibration failure, a high ionospheric metric, a high number of dead tiles or with data format issues. We also included the zenith pointing observations from 2014 (spanning August to September). Therefore, we began our fiducial analysis with 767 successfully calibrated snapshots from 2015 and 2014. This data set is then refined as described in Section 3.4. The final upper limits presented in Section 4.5 are obtained from 459 selected observations (∼14-h integration), covering five central pointings (including the zenith). There are 171 (⁠|${\sim}37{{\ \rm per\ cent}}$| of total) zenith observations in the final set.

2.4 Target field of EoR1

The targeted fields for EoR studies with MWA have been chosen to be in relatively cold parts of sky, but the presence of some bright radio sources such as Fornax-A in EoR1, shown in Fig. 1(b), is unavoidable. The bright and complex radio emission of Fornax-A [total flux density of |${\sim}500\, \mathrm{Jy}$| at |$189\, \mathrm{MHz}$| (Bernardi et al. 2013)], shown in Fig. 1(c), makes the process of calibration and foreground removal more complicated. Previous MWA studies of this field include McKinley et al. (2014) who modelled the spectral energy distribution of Fornax-A and Procopio et al. (2017) who improved the source catalogue in the EoR1 field by combining both GLEAM and GMRT |$150\, \mathrm{MHz}$| data. Since Fornax-A is an extended source with spatial and spectral structure, it is essential to develop a model which effectively represents its complicated morphology. This is further discussed in Section 3.3.

In contrast, the EoR0 field does not have anything as bright or extended as Fornax-A, and thus has been used more in MWA EoR analysis (Dillon et al. 2015b; Beardsley et al. 2016; Barry et al. 2019; Li et al. 2019; Trott et al. 2020). We compare and contrast analysis of the two fields in Section 4.2.

3 DATA CALIBRATION AND POWER SPECTRUM CALCULATION

Beardsley et al. (2013) showed that a full season of MWA observation on two fields, yielding an integration time of |$900{+}700\, \mathrm{hour}$|⁠, can potentially detect the EoR signal with an S/N of 14 on the amplitude. Were the measurements only noise-dominated, the larger integrations would average down the uncorrelated noise and lead to detection. However, in practice, the systematic errors, originated both from astrophysical and non-astrophysical sources, largely contribute to the power, exceeding the thermal noise level. There have been many efforts in the community to understand about different systematics and the effective strategies to mitigate them; e.g. the residual power due to imperfect calibration models (Offringa et al. 2015; Patil et al. 2017; Procopio et al. 2017; Li et al. 2019; Morales et al. 2019) and the spectral contamination from faint unmodelled sources (Barry et al. 2016), excess power due to imperfect telescope beam models (Beardsley et al. 2016; Barry et al. 2019; Li et al. 2019; Joseph et al. 2020), the contamination due to RFI (Offringa et al. 2015) and following possible contamination due to RFI excision (Offringa, Mertens & Koopmans 2019), the introduced spectral structure due to instrumental chromaticity of interferometers (Trott & Wayth 2016) and the power bias due to ionospheric distortion effects (Jordan et al. 2017; Trott et al. 2018). Trott et al. (2020) found that integrating beyond 300 observations (up to 820) did not significantly improve their results, implying the existence of a systematic floor of less than |${\sim}26\, \mathrm{h}$|⁠.

While systematic contamination can be mitigated through different techniques employed in analysis pipelines; it is also possible to employ data-quality metrics that detect and exclude highly contaminated data as described in Section 4.

Analysis of a data set to measure the power spectrum of the EoR signal requires three significant processes: (i) flagging of data affected by RFI, poor observing conditions or other anomalies, (ii) calibration and (iii) calculation of the power spectrum. We implement the RTS/CHIPS processing pipeline to analyse the data in this work (Jacobs et al. 2016). In this section, the current state of the pipeline is described, schematically shown in the diagram 2.

Figure 2.

The data flow from correlator to software pipeline: The visibilities are passed to AOFLAGGER to be RFI-flagged, then to RTS to be calibrated and then to CHIPS to measure the power spectrum.

3.1 RFI flagging

It is essential to mitigate the effect of RFI by detecting the contaminated data prior to processing. AOFLAGGER⁵ Offringa et al. (2010), Offringa, Van De Gronde & Roerdink (2012) is the primary algorithm applied to detect and flag the RFI-contaminated data. In addition, the algorithm also flags five 40 kHz channels in each coarse band: two 40 kHz channels on either side of the band to avoid the aliasing effect introduced by the channelizing procedure, and the middle 40 kHz band that has a DC offset (Offringa et al. 2015). Practically, the RFI occupancy for most channels is either very low (⁠|${\lt}5{{\ \rm per\ cent}}$|⁠) or very high (⁠|${\gt}90{{\ \rm per\ cent}}$|⁠) in case of contaminated data. Therefore, at the end, we fully flag any channel that contains more than 80 per cent RFI occupancy to conservatively exclude any highly contaminated data.

In addition, as will be discussed in Section 4.3, we apply the technique of Sky-Subtracted Incoherent Noise Spectra (SSINS; Wilensky et al. 2019), to mitigate the ultra-faint RFI contamination. SSINS estimates the noise by taking advantage of the negligible sky variance over short-time intervals of the interferometer measurements; and hence subtracts two time-adjacent visibility measurements and then averages them for all the baselines in the observation. This tool provides more sensitivity to interference, such as digital TV (DTV), hidden below the thermal noise level of individual baselines.

3.2 RTS calibration pipeline

The data is calibrated using the Real Time System (RTS; Mitchell et al. 2008; Ord et al. 2010). The RTS was initially developed for real-time calibration and imaging of MWA data, but is now employed as an offline pipeline.

The RTS has a number of features aimed specifically at dealing with large low-frequency data including: measurement and correction of ionospheric effects, direction-dependant beam models, baseline-dependant averaging, use of point and extended source models, as well as parallelization over frequency.

The RTS may be used in one of two modes, producing either calibrated residual visibilities or images. In this work, the first mode is used to create calibrated visibilities, which will be added together on a uvw-grid in CHIPS to form the power spectrum.

A detailed description of the algorithm is discussed in Mitchell et al. (2008), which is reviewed here. The RTS accomplishes the calibration process through two major steps;

• Step 1: For each observation, a compound calibrator source is created, made up of the 1000 apparently brightest sources available in the source catalogue. The source catalogue is created by cross-matching with PUMA (Line et al. 2017). Full details of PUMA and this cross-matched catalogue can be found in Line et al. (2020), but we include a brief description in Section 3.3. The single compound calibrator provides a model for the visibilities. By fitting the observed visibilities to this model, one can determine the direction-independent complex voltage gain or Jones matrices for each antenna.

• Step 2: Bright sources, in comparison to faint sources, are less affected by flux contamination via the point-spread function of the remaining sources in the field, and thus should be peeled first to reduce subsequent error. Thus, all sources are ranked in order of apparent brightness. Then the model visibilities of the 1000 brightest sources are subtracted from the visibility set and then each source is individually peeled/calibrated in a Calibration Measurement Loop. The Calibration Measurement Loop performs these steps for each source individually in subsequent order: (a) the model visibilities for the source are added back to the visibility set and the phase centre of the observation is rotated toward the actual position of the source; (b) the gain and positional offset of the source due to ionospheric refractive effects are measured and corrected; (c) the full direction-dependent calibration is computed and applied, but only for the five strongest sources due to computational limitations. For the rest of the sources, only ionospheric corrections are computed and applied; (d) once the measurements converge, they are used to subtract the source from the visibility set, thus improving the model visibility iteratively. This step is performed on an 8s cadence. Averaging the visibilities to 8s resolution increases the SNR while provides the required time resolution for sampling ionospheric fluctuations which is less than a minute. Once the calibrated residual visibilities are created, they are averaged to 80 kHz in frequency for power estimation. This again provides higher SNR while resolves in the expected spectral fluctuations of EoR.

3.3 Source catalogue and extended source modelling

Three sky models, which we label puma2017a, puma2017b, and puma2020, are used in this work to measure the effect of improved sky models on the power spectrum limits. All three used the GLEAM catalogue (Hurley-Walker et al. 2017; Wayth et al. 2015) as a base, cross-matched using PUMA to other catalogues to provide additional positional and spectral information (for details see Line et al. 2020). The first sky model, puma2017a, is based on the work of Procopio et al. (2017), which used TGSS GMRT (Intema et al. 2017) images to create higher resolution source models within ∼10° of the EoR1 field centre. The GLEAM catalogue was still being finalised at the time of this publication. Once the GLEAM catalogue had been finalised, the extended models from Procopio et al. (2017) were retained, and all other sources were replaced with updated models based on the full release version of GLEAM. We label this model puma2017b. RTS has the capability of digesting models of extended sources such as shapelet-based ones. A shapelet-model for Fornax-A is included in the sky model which is used in both calibration and subtraction. For both puma2017 sky models, the same model of Fornax-A was used. This was a shapelet model, based on phase I MWA data, as detailed in Riding, Mitchell & Webster (2017).

The puma2020 sky model includes an updated Fornax-A shaplet model, which utilizes a combination of MWA Phases I and II data to increase angular resolution. For details on this model, please see Line et al. (2020). In addition, an in situ adjustment was made to the flux density and spectral index of Fornax-A in order to better match the observed flux density and spectral index in a sample of zenith data.

The EoR1 limits published in Trott et al. (2020) used previously calibrated and peeled data products that were gradually processed over 2 yr, using whichever sky model was available at time of processing. A re-processing of that entire data set would require a large investment of computational resources. We estimate 50 per cent of the data in the Trott et al. (2020) EoR1 limit used the puma2017a model, and 50 per cent the puma2017b model.

All of the observations used in the calculation of final limits in this paper are processed with the puma2020 sky model.

3.4 Data metrics

The RTS calibration process also provides a number of metrics which we use to evaluate data quality and determine the suitability of each processed observation for power spectrum estimation. The three most important metrics are as follows:

• Calibration-based metrics: RTS provides the amplitude and phase of complex gain solutions for each tile. These solutions are a useful tool to evaluate the calibration quality, since they reflect how smoothly the calibration procedure is performed over the frequency for each tile, both in shape and amplitude. In addition, we can diagnose any tile with non-consistent profile relative to others.

• IonoQA: The state of the ionosphere during the observation can be measured as a single value for ionospheric activity by combining the median source offset and amount of offset directionality (Jordan et al. 2017). This ‘IonoQA’ metric value varies from 2 (low ionospheric activity) to multiple 10s (high ionospheric activity). Trott et al. (2018) demonstrated the need to exclude ionospherically active data to avoid bias in cosmological power-spectrum measurements.

• Visibility noise rms: Since calibrated visibilities are averaged to 8-s resolution, each of 112-s observations consists of 14 measurement sets or equivalently 7 pairs of opposing even/odd cadences.⁶ Visibility noise rms is calculated as the standard deviation of differences between two sets of visibilities in each pair and is denoted as σ_vis. The reported rms values throughout this work are calculated, while the first and last time-samples are excluded due to some extra known errors, introduced by partial flagging or missing data samples. In addition, three 80 kHz channels at the beginning, centre, and end of each coarse band are excluded once calculating the rms to avoid the discontinuity due to the effect of either aliasing or DC offsets. The related frequency gaps in the x-axis of rms plots, presented in Section 4, are also removed.

3.5 Power measurement – CHIPS

In this work, CHIPS (the Cosmological H i Power Spectrum Estimator; Trott et al. 2016) is used to measure the power spectrum. It is an inverse variance estimator that computes the power from gridded, calibrated visibilities.

As mentioned before, we are attempting to study the physical properties of the hydrogen gas during EoR by probing the statistics of its temperature field, which are affected by the reionization process. By denoting temperature fluctuation (relative to the mean) as |$T(\boldsymbol{r})$| at vector position |$\boldsymbol{r}$|⁠, we represent the autocorrelation function of the temperature fluctuations signal as |$\xi (\boldsymbol{r})$|

$$\begin{eqnarray} \xi (\boldsymbol{r}) = \langle T(\boldsymbol{r_0}) T(\boldsymbol{r_0} + \boldsymbol{r}) \rangle _{\boldsymbol{r_0}} \end{eqnarray}$$

(1)

in which <.> denotes the averaging over the sample volume, Vol. The Gaussian statistics of a signal is well described by its power spectrum, |$P(\boldsymbol{k})$|⁠, in k-space. Under the assumption of the isotropic 21-cm signal, it provides the signal power spectrum over different spatial scales of k (⁠|${=}|\boldsymbol{k}|$|⁠).

The power spectrum is defined as the Fourier Transform of the autocorrelation function:

$$\begin{eqnarray} P(k) = \int _{\rm Vol} \xi (\boldsymbol{r}){\rm exp}(-2\pi i \boldsymbol{k}\,\cdot\,\boldsymbol{r}) {\rm d}\boldsymbol{r} . \end{eqnarray}$$

(2)

By replacing the autocorrelation function, the power spectrum equals to the spatial covariance of the signal, integrated over the spatial volume Vol:

$$\begin{eqnarray} P(k) = \frac{1}{\rm Vol}\langle \tilde {T}(\boldsymbol{k})\tilde {T}^{*}(\boldsymbol{k}) \rangle . \end{eqnarray}$$

(3)

On the other hand, the interferometer measures the sky brightness as visibilities in u, v, w domain. Given the flat-sky approximation, the measured visibilities, V(u,v) are the double Fourier Transform of the sky brightness, I(l, m). This can be represented as

$$\begin{eqnarray} V(u,v) = \int \int _{\rm source} I(l,m) {\rm e}^{-2\pi i (ul + vm)}{\rm d}l {\rm d}m, \end{eqnarray}$$

(4)

where u, v are Fourier modes of the measured visibility and l, m are the corresponding angular scales, which are represented by modes perpendicular to the line of sight, or the angular modes,k_⊥, in k-space. Additionally, the modes parallel to the line of sight, k_∥, may be mapped with spectral channels of measurement, while |$k^2 = k_{\perp }^2 + k_{\parallel }^2$|⁠. The sky brightness is linearly related to brightness temperature at radio regime and, hence, the Power Spectrum can be obtained as

$$\begin{eqnarray} P(k) = \frac{1}{\rm Vol} \langle \tilde {V}(\boldsymbol{k})\tilde {V}^{*}(\boldsymbol{k}) \rangle . \end{eqnarray}$$

(5)

CHIPS starts off by performing the coherent integration of calibrated visibilities via gridding them all on to uv-plane per w-snapshot. This procedure creates the data cube in (u, v, f). Next, each cell is Fourier Transformed along frequency and finally squared.

Measured power spectrum and upper limits are reported as the 1D power spectrum which is defined as

$$\begin{eqnarray} \Delta ^2(k) \equiv \frac{k^3}{2\pi ^2}P(k) , \end{eqnarray}$$

(6)

and is the integrated total power over spatial scale k in mK².

3.6 Power spectrum averaging scheme

While the spherically averaged 1D power spectrum (over 3D spatial scales of |$k=|\boldsymbol{k}|$|⁠) provides the required cosmological measurements, the 2D cylindrically averaged power spectrum is an informative tool to learn about the power at different regions of 2D space. Fig. 3 shows an example 2D power spectrum. It provides the power distribution over angular (k_⊥) and line-of-sight (k_∥) modes and will be used throughout this work for diagnostics. Since foregrounds vary smoothly as a function of frequency, they contaminate the low k_∥ modes. However, the instrument chromaticity couples the foregrounds into higher k_∥ modes which shapes a characteristic wedge in 2D space, called foreground wedge. We can perform our measurement in the region above foreground wedge which is called the EoR window.

Figure 3.

A typical 2D power spectrum for an integration of 20 observations. The region that is marked with the dashed line shows the selected range of (k_⊥, k_∥) of which the measurements are included in averaging from 3D power spectrum to 1D.

As discussed in Barry et al. (2019), Li et al. (2019), and Trott et al. (2020), cuts in k-space need to be implemented when averaging to the 1D power spectrum to exclude regions contaminated by foreground and instrumental systematics. Specifically these are as follows:

• 12λ > k_⊥ > 50λ; to avoid contaminated low k_⊥ and poor uv-coverage at high k_⊥.

• |$k_{\parallel } \gt 0.1\, \mathrm{h\, Mpc^{-1}}$|⁠; to avoid the still present leakage from the foreground wedge at low k_⊥.

• k_∥ > 3.5k_⊥; to shift the horizon limit slightly higher preventing sub-horizon leakage.

The selected k-space region is marked by a dashed line in Fig. 3. The same cuts and binning scheme are applied throughout the paper to average the 3D k-space to 1D.

3.7 Improvements to calibration pipeline

In this section, we demonstrate the major changes to the calibration pipeline, including those of the source catalogues and the peeling algorithm, and their effect in 1D power space.

We start off by showing the 1D PS for a single typical observation processed with the puma2017a source catalogue (as described in Section 3.3). This processing is similar to that performed by Trott et al. (2020) and is shown as the blue line in Fig. 4.

Figure 4.

The improvement of data reduction as shown by comparing 1D power spectrum of one typical zenith observation processed by different versions of the source catalogue and calibration pipeline. The dashed line shows the thermal noise curve.

The first improvement is due to the peeling algorithm. In order to effectively remove foreground sources, we set the RTS to peel the 1000 brightest sources in the field, in the order of their apparent brightness. However, due to a numerical issue that only affected peeling of faint sources, peeling fewer sources (orange power spectrum) would result in less contamination in comparison to complete 1000-source peeling (blue power spectrum). This difference is due to the fact that the peeling process with fewer sources does not go through the procedure of faint source peeling. Previous EoR1 processing in Trott et al. (2020) had some observations with fewer number of sources peeled.

Having resolved the numerical issue with the peeling of faint sources, the effective peeling of 1000 sources is reflected by a less contaminated power measurement shown by the green spectrum. Also, RTS processing typically treats the five brightest sources as full direction-dependant (DD) calibrators. Although reducing the number of DD calibrators has been shown to improve the results for ultra-low data (Yoshiura et al., submitted), the results were inconclusive for high-band data.

Another improvement is based on the availability of dead dipole information from 2014 onward. Dead dipoles occur when the low-noise amplifier on one of the cross dipoles fails. RTS forms the beam model for each tile individually. So the statistics of dead dipoles of each single tile provides us with the information to adjust the beam profile of the tile and the visibility weights. By updating the beam model for each tile, we take into account the effect of any dead dipoles, if it contains any. Without this information, additional sources of error are included (Joseph et al. 2020). Therefore, by flagging the dead dipoles within each tile, we obtain a less contaminated power spectrum as shown by red line.

Finally, upgrading the sky models further improved the results by mitigating the foreground contamination. The purple line is the result of processing with puma2017b that includes updated sky model, and the brown line is the result of processing with the puma2020 catalogue that includes the upgraded shapelet model of Fornax-A based on its Phase-I/II image (refer to Section 3.3). Puma2020 catalogue is the basis of this work. There is ∼6 times improvement in power measurement at |$P(k=0.1\, \mathrm{h\, Mpc^{-1}})$| due to improvements in the source catalogue and peeling process.

4 RESULTS

In this section, we demonstrate some of the prominent differences between the MWA EoR0 and EoR1 fields, as well between the 2015 data analysed here and the 2013 data analysed in some previous MWA results. We then present a series of data quality metrics used to identify contaminated data and select a best sample for power estimation. Our final power spectrum results are presented in Section 4.5.

4.1 Comparison of 2013 and 2015 Data

Previous MWA EoR data analysis works such as Dillon et al. (2014, 2015a), Beardsley et al. (2016), and Barry et al. (2019) have been focused on data observed during 2013, while Li et al. (2019) and Trott et al. (2020) studied 2016 and 2013–2016 data, respectively. In this work, we aim to process the data observed during 2015 (and 2014). Therefore, it is important to investigate the variability of data between the 2013 and 2015 seasons. We compare two data sets including nine zenith-pointing observations of a single night, observed in each of these two years. Both data sets are selected according to the metrics elaborated in Section 4.3.

Fig. 5(a) illustrates the visibility noise rms as a function of frequency of the two data sets, indicating that 2013 visibilities are more contaminated than 2015 ones. While the correlator temporal resolution of 2013 and 2015 data is |$0.5\, \mathrm{s}$| and |$2\, \mathrm{s}$|⁠, respectively, the main difference between the data of these two seasons is due to the recording of dead dipoles (as discussed in Section 3.7). There are ∼50 tiles with one dead dipole in 2015 data and the tiles with more than two dead dipoles are flagged (Joseph et al. 2020). For the 2015 observations, the beam shape and visibility weights are adjusted if there is a dead dipole recorded, but that is not the case with 2013 observations. Moreover, there is one local jump in visibility noise rms values in 2013 observations which is due to a discontinuity in applied digital gains. This jump is not noticeable for the 2015 data where smoother digital gains are applied across the band. Also, increased rms values are noticeable at the coarse band edges where digital non-linearities produce varying coarse band shapes. To investigate the effect of the application of dead dipole information, we disabled the feature of dead-dipole flagging in the 2015 data. In that case, as it is shown in Fig. 5b, the rms values of 2013 and 2015 observations generally agree with each other, though there are still some differences due to dissimilar dead dipoles in the two data sets.

Figure 5.

Comparison of visibility noise rms of 2015 and 2013 data.

The effect of inclusion of dead-dipole information in 2D and 1D power space is shown in Fig. 6. The 2D power spectra comparison of 2013 and (dead-dipole flagged) 2015, shown in Fig. 6(a), generally indicates less contaminated power measurements for 2015. Moreover according to Fig. 6(b), while the 1D power spectrum of 2015 without dead-dipole flagging mostly agrees with that of 2013, it is noticeably (1.3 times at lowest k-modes) more contaminated than that of dead-dipole flagged data from 2015.

Figure 6.

Power spectrum comparison between 2015 and 2013 data.

Both the visibility noise rms and the power spectrum confirm there is less contamination in 2015 data relative to 2013 data, thereby reducing the systematic errors in this analysis of the 2015 data in comparison.

4.2 Comparison of EoR0 and EoR1 Fields

The MWA EoR0 field has been broadly analysed and studied in previous works such as Dillon et al. (2015a), Beardsley et al. (2016), Barry et al. (2019), Li et al. (2019), and Trott et al. (2020). However, the EoR1 data has only been studied by Paul et al. (2016) and Trott et al. (2020), with the former analysing a 3-h integration with the CASA/Delay-Spectrum pipeline and the latter analysing ∼19 h of zenith-pointing observations, using the RTS/CHIPS pipeline. Zhang et al. (2020) also studied the performance of different calibration methods on MWA Phase-II data from EoR0 and EoR1 fields. In this work, we aim to analyse EoR1 data with the RTS/CHIPS pipeline over a wider range of beam pointings with an improved version of the RTS, higher fidelity source catalogues and improved data-quality metrics. Before focusing on EoR1, we consider some of the main differences between these fields. The sky temperature of the two fields are also compared in Appendix  A.

4.2.1 Per pointing comparison

One of the underlying assumptions of all EoR experiments is that the calibrated visibilities are stable in time, allowing for long integrations which can attain the sensitivity required to detect the faint cosmological signal. However, as shown in Beardsley et al. (2016) and Barry et al. (2019), there are prominent differences visible between the different MWA pointings for EoR0. In this section, we present differences visible between different pointings both for EoR0 and EoR1, and across different nights for EoR1. We illustrate these differences through the median of the visibility noise rms (rms hereafter) and the IonoQA (defined in Section 3.4).

Fig. 7 displays the per-observation rms and IonoQA, as well as the per-pointing 1D power spectra of three representative nights. Two nights (2015 September 24 and 12) from EoR1 are presented in the first two rows, and the third row presents a night (2015 September 24) from EoR0.

Figure 7.

The columns show the median rms, IonoQA, and 1D power spectrum, respectively. The first and second rows show the EoR1 observation nights 2015 September 24 and 12, respectively, while the third row shows the EoR0 observation night 2015 September 24. Each data point indicates one observation and is colour-coded by pointing: Minus2 (blue), Minus1 (orange), Zenith (green), Plus1 (red), and Plus2 (purple). The dashed lines in power spectra plots represent the thermal noise.

Comparing the EoR1 nights, September 24 (first row) has a linear and predictable variation of rms within each pointing transition. We interpret this as the primary beam repeatedly passing over regions of varying galactic emission. For the September 12 data (second row) the behaviour is more complicated; the linear, periodic structure still appears but some unknown effect is also introducing additional structure within each pointing. There appears to be no significant difference in the IonoQA values for these nights.

For the EoR0 observations from September 24 (third row), we observe an overall decline in the rms values over the course of the night. We interpret this as lower sky temperature associated with the setting of the Galactic Centre, as predicted in Fig. A1. The overall rms for EoR0 observations is also higher than for EoR1, as discussed in Section  A. For this night, the IonoQA also appears to show some larger variations than seen in the EoR1 data later in the night.

A small number of observations with IonoQA value of greater than 5 do not show up in the IonoQA plot as it is limited to quality assurance criteria 4.3.

We note that the 1D power spectra of different pointings over these typical nights do not show a systematic intra-pointing variability either on different nights or on different spatial scales.

Fig. 8 shows the rms across the whole frequency band for three observations from each pointing of EoR1 (top panel) and EoR0 (bottom panel). The rms generally declines as we move towards higher frequencies as expected from the lower sky temperature. However, in the case of off-zenith pointings, the rms does not follow a uniform trend and tends to slightly increase at high frequency. The slope increases for observations further from the zenith, even more so for the Minus3 and Minus4 pointings (not presented in this figure). As discussed in Beardsley et al. (2016), although the MWA EoR fields are primarily devoid of Galactic plane contamination, the sidelobes of the beam are highly contaminated by Galactic emission when it rises above the horizon, one of the complications of wide-field instruments. This condition is significantly worse in before-zenith pointing transitions such as Minus4 to Minus1. EoR1 shows widefield contamination on either side of zenith, albeit at a lower level.

Figure 8.

RMS over frequency per pointing of (a) EoR1 and (b) EoR0.

Per-pointing analysis of EoR1 data shows that the five central pointings are behaving consistently and thus approved for further analysis to measure the limits.

4.2.2 Integrated 2D power spectrum comparison

The power spectra of a ∼1-h integration of EoR1 and EoR0 (6 observations from each of five central pointings) are compared in Fig. 9. The 2D power spectrum comparison (top panel) indicates higher contamination in EoR window and so generally higher 1D power (bottom panel) for EoR1. However, the power in the wedge is higher for EoR0 and is particularly noticeable along the horizon in the ratio plot (middle panel). Thyagarajan et al. (2015) demonstrated that the wide-field measurements are highly sensitive to diffuse emission at the horizon. This feature is due to the shortening of baselines near the horizon that makes them more sensitive to larger angular scales. They also subtend larger solid angles near the horizon and so capture a larger integrated emission. Therefore, we expect Galactic emission, which is mostly present at off-zenith transitions of EoR0 (Trott et al. 2020), to introduce higher contamination if the foreground modelling and removal is not perfect. This excess power is leaking into large spatial scales or equivalently small k_⊥ areas of the horizon in the ratio plot, and is captured by lowest bin of 1D power spectrum. This difference shows up as lower EoR1 power at lowest spatial modes, or equivalently large-scale measurements in 1D power spectrum (bottom panel).

Figure 9.

EoR1 and EoR0 power spectrum (top panel) and their comparison in 2D (middle panel) and 1D (bottom panel) power space.

Having improved the analysis tools that would be promising to perform further deep measurements of EoR1 field and validate its capability to less contaminated limits in comparison to EoR0, at least at large scales.

4.3 Data quality metrics

All 767 candidate observations were successfully calibrated. In order to select the best data set for power spectrum analysis, we need to filter out any highly contaminated observations. We combine a number of criteria applied by previous MWA EoR analyses together with novel methods of identifying anomalous data. Specifically, our selection is performed over six stages, as follows: (i) autocorrelation amplitudes, (ii) calibration-based metrics, (iii) window/wedge power metric, (iv) individual 1D power spectrum, (v) ultra-fine RFI, (vi) phase component of gain profile. Each one is described below in more detail and summarized in Table 2.

• Autocorrelation profiles

  Examining the (time-averaged) autocorrelation amplitudes of individual tiles readily identifies features such as faulty unflagged tiles or digital TV (DTV) interference. Since the MWA is fairly stable with well-behaved autocorrelations, inferences about extreme outliers and large changes in frequency can be made even prior to calibration. The autocorrelation amplitudes generally vary smoothly over frequency. For all tiles with the same cable-group,⁷ the amplitudes are distributed with a maximum deviation of roughly |$20{{\ \rm per\ cent}}$| around the mean value. However, some tiles are extreme outliers or have large deviations over short frequency-scales as high as |$25-30{{\ \rm per\ cent}}$| of the mean value. Thus, the associated observations were removed from the data set.

  To learn more about the large variations appearing in some tiles, we measured the standard deviation of the autocorrelation amplitudes over time as a function of frequency. We noticed that some tiles, specifically 55, 116, 125, and 163, show large deviations in some frequency channels. The detected standard deviation values might be as large as those produced by DTV, as illustrated in Fig. 10, and are caused by a sharp peak at a specific time/frequency of the amplitude profile. Since such features are not flagged by AOFlagger, we probed their effect on power measurements by flagging these tiles. Their effect on the measured power spectra is not noticeable at this stage; implying that the level of their introduced systematics is below the sensitivity of current analysis. However, with future development of higher precision techniques and larger data sets, removing these observations from the data might reduce the systematic floor.

• Calibration-based metrics

  At this stage, all of observations are probed in terms of the following:

  • Calibration solutions

    The RTS-derived bandpass and gain solutions for each observation are examined. In general, the amplitude and phase components of these complex quantities are consistent for all the tiles within one observation. Observations with at least one tile with an abnormal profile, typically large or sudden variations in their gain amplitudes are excluded. In addition, observations with more than six flagged tiles are removed.

  • Ionospheric activity

    Based on the studies by Jordan et al. (2017), Trott et al. (2018), the ionospherically active observations with large offset and anisotropy measure should be avoided. So we exclude observations with IonoQA greater than 5.

  • Visibility noise rms

    Visibility rms (discussed in Section 3.4) generally varies smoothly over frequency. However, highly contaminated observations whose profile fluctuation over at least one fine channel was larger than |$40{{\ \rm per\ cent}}$| of the overall rms change over the whole band were excluded.

• Window/wedge power metric

  Next, we apply the window/wedge metric as developed by Beardsley et al. (2016) and employed in previous studies such as Trott et al. (2020). This is measured by a delay transform power spectrum estimator, computed incoherently across all short baselines and treats each baseline separately. Window power estimates the power in a selected region of the EoR window; between the main beam lobe and the first coarse channel harmonic. Wedge power measures the power in the primary beam main lobe wedge. The cuts on window power pick values within 2σ range of the distribution, as shown in the Fig. 11(a). The cut is also shown in the scatter plot of window power versus rms in Fig. 11(b). We also apply 2σ cut on wedge power as shown in Fig. 11(c).

• Individual 1D power spectra

  Having examined the incoherent delay spectrum of observations in the last step, we need to verify their power spectra which is a coherent measurement of the gridded visibilities. Hence, the 1D power spectrum of each observation is individually measured, for consistency with previous work (Trott et al. 2020). Fig. 12 shows the histogram of P(lowest-bin) for candidate observations. The long tail of this histogram shows the clearly contaminated observations with high power measurements raised by some systematics. We apply Chauvenet’s Criterion to detect the outliers. This method cuts the outliers with P(k = 0.1) > 1.35 × 10⁵mK² in the long tail, as shown with dashed red line in the figure. We examined the effect of their inclusion in the whole data set power spectrum and see a drop in the total power implying they are clearly contaminated and need to be excluded from the total set.

  Having cleaned the data set by applying the last four metrics, 459 observations are selected. We can be even more conservative by considering two more subtle metrics, as follows:

• Ultra-faint RFI

  As described in Section 3.1, standard MWA RFI flagging is handled by AOFlagger in this work. However, we are able to further detect significant ultra-faint RFI contamination occurrences in the data using SSINS. While we expect a Gaussian distribution for the mean-subtracted noise spectrum of a non-contaminated observation, as illustrated in Fig. 13(a), the most repeated occurrences of contaminated cases are either bright DTV interference as in Fig. 13(b) or a time-variant source of interference present over the whole observation time lapse and frequency band as shown in Fig. 13(c). In total, 120 observations with abnormal SSINS profiles are identified by visual inspection. The abnormal cases include those illustrated in Figs 13(b) and (c) and presence of narrow-band RFI over time or frequency.

• Phase component of gain profile

  There are some observations that show excessive phase wrapping on the tiles associated with one or more receiver. This issue is detected in 118 observations.

Figure 10.

The standard deviation of the autocorrelation amplitude over time as a function of frequency for tiles 152, 163, 125, and 55 (respectively, from top to bottom), colour-coded for 11 observations. Tile 152 is an example of a healthy tile, while the other three tiles show spikes at some frequency channels. The observation plotted in blue is DTV interference.

Figure 11.

Window/wedge power cut-off.

Figure 12.

Histogram of P(lowest bin) for candidate observations.

Figure 13.

Sky-subtracted Noise Spectrum of three observations with normal (top), DTV-interfered (middle), and biased (bottom) noise spectrum. Western Australia channel 6 (⁠|$177.5\, \mathrm{MHz}$|⁠), channel 7 (⁠|$184.5\, \mathrm{MHz}$|⁠), and channel 8 (⁠|$191.6\, \mathrm{MHz}$|⁠) fall within MWA high-band frequency range.

Table 2.

Data quality metrics that were applied to the data set in order and the associated number of passed observations. The final power spectrum is obtained from 459 selected observations that are selected by first four metrics. The last two metrics of (v) and (vi) are applied for further systematic analysis purpose.

Metric	Description	Number of passed observation
	Successful calibration	767
i	Auto-correlation	758
ii	Calibration-based metrics:	522
	solutions, IonoQA, rms
iii	Window/wedge power	494
iv	Individual power spectrum	459
v	SSINS	339
vii	Phase component of gain	221

Having applied all six metrics, 221 observations are selected.

As the inclusion of the data detected by metric (v) and (vi) in the data set does not make any noticeable change in the power spectrum, we could not verify if they are definitely a source of systematic contamination in the measurement at this stage. However, they provide a cautiously refined data set for further precise analysis of systematic contamination as discussed in next Section 4.4.

4.4 Data selection based on visibility rms

The main limitation on detecting the EoR signal at this time are the (un)known systematics. Visibility noise rms could represent a combination of different systematics, as seen in the different receiver temperatures inferred in Section  A. So we explore the residual systematic errors based on this metric. To make sure that our data set is homogeneous, we perform this test only for the 221 observations that have been filtered by all the aforementioned six metrics. As shown in Fig. 11(b), for the 459 observations used in our final results, the median rms values are centred at |${\sim}14.8\, \mathrm{Jy}$| within a range of ∼3 Jy. By splitting the data into two groups of low-rms and high-rms with median rms value of lower and higher than |$14.8\, \mathrm{Jy}$|⁠, respectively, we see significantly different results for their power measurements (Fig. 14). Two data sets with the same number of observations from each rms group were selected. The 1D power spectrum of low-rms data set is significantly lower than the high-rms one, shown in Fig. 14(a). To investigate the effect of presence of high-rms data in the whole data set, we also compared the power spectrum of a data set including equal amounts of low-rms and high-rms data with another data set containing only low-rms data, shown in Figs 14(b) and (c). This demonstrates the increased level of contamination in the case of high-rms visibilities.

Figure 14.

Power spectrum comparison between low- and high-rmsbservations.

There are 167 low-rms observations in the best 221-observation data set. The lowest upper limits obtained from the low-rms data is |${7.34\times 10^3}\, \mathrm{mK^2}$| at |$k=0.13\, \mathrm{h\, Mpc^{-1}}$| and z = 6.5. Among the 221-observation clean data set, the high-rms observations mostly belong to a particular part of the observing season. For example, |$100{{\ \rm per\ cent}}$| of 2015 September observations exhibit of low rms while all 2014 September and August observations are all of high rms. As the high-rms observations are concentrated in well-defined time-frames, they might be affected by some hardware artefact that occurred and lasted during that time. Having learned about these effects, a future approach might be to diagnose and avoid the observations of those faulty nights prior to mass processing.

4.5 Final power-spectrum measurements

To avoid the cosmological evolution of the signal over the bandwidth of interest, we need to measure it within an interval that assures the signal ergodicity. We pick sub-bands of width |$15.36\, \mathrm{MHz}$| that will be tapered to an effective width of |${\sim}10\, \mathrm{MHz}$| once a Blackman–Harris window function is applied. Hence, we will have three sub-bands with central redshifts of 6.5, 6.8, and 7.1, each including 192 fine frequency channels.

In the search for data which is only noise-dominated, we explored six metrics to exclude highly contaminated data. While the first four metrics are required to assure the data is not contaminated, we also conservatively applied the last two metrics to eliminate the possible finer signatures of systematics in the data. 459 observations passed the first four metrics but only 221 pass also the metrics (v) and (vi). Due to the relatively low integration time of selected or finely selected data sets, the possible effects of detected systematics would not be evident in the final power measurement. However, the lowest 2σ upper limits for 221-observation data set is |${7.29\times 10^3}\, \mathrm{mK^2}$| at |$k=0.13\, \mathrm{h\, Mpc^{-1}}$| and z = 6.5.

We confirmed the accuracy of the power normalization through a set of simulations. As is shown in Appendix  B, the simulated EoR signal is consistently recovered through RTS/CHIPS pipeline. The measured power and 2σ upper limits of 459 selected observations are presented in Table 3 and Fig. 15. This provides us with the lowest limit of |$\Delta ^2 \le {5.44\times 10^3}\, \mathrm{mK^2}$| at |$k=0.13\, \mathrm{h\, Mpc^{-1}}$| and z = 6.5.

Figure 15.

Measured 1D power spectrum (solid lines) and 2σ upper limits (dashed) for the 459 observation at z = 6.5 (red), z = 6.8 (green), and z = 7.1 (blue).

Table 3.

Measured power (Δ²), 2σ upper limit (⁠|$\Delta _U^2$|⁠), and the square root of the upper limit (Δ_U) for three redshifts. The reported k values point to the bin centre.

	z = 6.5			z= 6.8			z = 7.1
k	|$\boldsymbol{\Delta ^2}$|	|$\boldsymbol{\Delta _U^2}$|	|$\boldsymbol{\Delta _U}$|	|$\boldsymbol{\Delta ^2}$|	|$\boldsymbol{\Delta _U^2}$|	|$\boldsymbol{\Delta _U}$|	|$\boldsymbol{\Delta ^2}$|	|$\boldsymbol{\Delta _U^2}$|	|$\boldsymbol{\Delta _U}$|
|$(\mathrm{h\, Mpc^{-1}})$|	(mK2)	(mK2)	(mK)	(mK2)	(mK2)	(mK)	(mK2)	(mK2)	(mK)
0.135	5.12E+03	5.44E+03	73.78	8.73E+03	9.28E+03	96.34	1.43E+04	1.52E+04	123.25
0.177	2.05E+04	2.13E+04	145.97	2.95E+04	3.06E+04	175.05	4.54E+04	4.71E+04	216.91
0.248	4.05E+04	4.23E+04	205.55	6.36E+04	6.61E+04	257.19	1.07E+05	1.11E+05	333.57
0.319	1.37E+05	1.44E+05	379.12	1.51E+06	1.55E+06	1246.86	2.33E+06	2.40E+06	1550.42
0.389	2.10E+07	2.16E+07	4649.95	8.76E+07	9.02E+07	9498.26	1.34E+08	1.38E+08	11758.27
0.460	1.04E+08	1.07E+08	10327.4	5.13E+07	5.29E+07	7273.08	8.05E+07	8.29E+07	9106.82
0.531	4.93E+06	5.08E+06	2254.25	3.72E+05	3.95E+05	628.23	6.53E+05	6.84E+05	827.07
0.602	3.58E+05	3.76E+05	612.94	6.54E+05	6.90E+05	830.54	1.33E+06	1.39E+06	1177.4
0.672	1.65E+06	1.71E+06	1305.86	2.20E+06	2.28E+06	1511.57	3.43E+06	3.56E+06	1886.78
0.743	1.72E+06	1.78E+06	1336.03	2.35E+06	2.45E+06	1564.64	3.26E+07	3.36E+07	5800.56
0.814	2.06E+08	2.12E+08	14555.44	2.94E+08	3.03E+08	17404.05	1.28E+09	1.32E+09	36371.56
0.885	7.52E+08	7.75E+08	27841.21	1.08E+09	1.12E+09	33406.26	5.76E+08	5.94E+08	24369.03
0.956	2.77E+07	2.86E+07	5348.69	3.95E+07	4.07E+07	6382.98	2.80E+06	2.95E+06	1718.88
1.026	9.88E+05	1.08E+06	1041.16	2.45E+06	2.59E+06	1607.81	5.22E+06	5.46E+06	2337.5

4.6 Comparison with previous results

The MWA high-band EoR1 data processing by RTS/CHIPS pipeline has only been previously studied by Trott et al. (2020). In this section, we discuss the results and address the differences in the two analysis approaches, as summarized in Tables 4 and 5.

Table 4.

Comparison of analysis approach of this work and Trott et al. (2020).

Feature	This work	Trott et al. (2020)
No. of selected observations	459	600
Observation time	2015/2014	2013-2016
Pointings	Five central	zenith
MWA phase of operation	Phase-I	phase-I/II
Source catalogue	puma2020	puma2017a/b
Shapelet modelling	New model based on phaseI/II image of Fornax	old model
IonoQA metric cut-off	5	30
Other applied metrics	Auto-correlation	N/A
	rms trend
	Calibration solution profiles

Table 5.

Comparison of the upper limits between this work and Trott et al. (2020) at |$k=0.14 \, \mathrm{h\, Mpc^{-1}}$| which is the lowest common k.

	|$\boldsymbol{\Delta _U}(mK)$| at k = 0.14 (⁠|$\mathrm{h\, Mpc^{-1}}$|⁠)
	z = 6.5	z = 6.8	z = 7.1
This work (459-observation)	145.970	175.05	216.91
Trott et al. (2020) (600-observation)	183.8	199.9	305.0

As discussed in Section 3.7, the upgraded version of RTS and source catalogues have significantly improved data reduction and power measurement results of EoR1 data. The Trott et al. (2020) study also includes 2013 data that suffers from lack of dead dipole information. With more conservative data refinement strategies, a less systematic-contaminated power measurement has been obtained. In Fig. 16, the 1D power spectra of a sample of 20 observations reduced with both versions of the pipelines are compared.

Figure 16.

Comparison of 1D power spectra of this work (blue) and Trott et al. (2020) (red) for 20 observations.

4.7 Discussion

Currently, the MWA EoR power spectrum measurements do not reach the expected thermal noise and hence are limited by systematics. In this work, a 14-h integration using an improved (RTS+CHIPS) analysis pipeline has been demonstrated the possibility of lowering this systematic floor. By deeper integration of observations and further mitigation of systematic errors, we will potentially be able to reach the noise-dominated level of signal power and ultimately detect the EoR signal. In this regime, the assurance of signal preservation will be critical and requires a full end-to-end simulation which is beyond the scope of this work and is a part of the collaboration research program. However, in the case of this work in which the integration time is short and measurements are higher than the noise level, the possible signal loss does not largely interfere with our results. Nevertheless the lack of signal loss issue is (in)directly validated by different probes: As studied in the recent deep measurements of EoR upper limits with RTS+CHIPS, Trott et al. (2020) showed that there is no signal loss occurring in the CHIPS pipeline. It is also shown that there is no signal loss occurring in RTS+CHIPS due to imperfect telescope modelling. In the measurements by the the parallel MWA EoR analysis pipeline (FHD+ εppsilon), lack of signal loss is also validated via simulation (Barry et al. 2016; Barry et al. 2019; Li et al. 2019). The consistent results from cross-checking the two independent parallel pipelines, assures the validity of signal loss test results for both pipelines (Beardsley et al. 2016; Jacobs et al. 2016; Barry et al. 2019; Li et al. 2019).

5 CONCLUSIONS

This work has presented the analysis of Phase-I MWA EoR data targeted at the EoR1 field using the RTS/CHIPS pipeline. This data includes five central pointing observations of the field in 2015/2014. By applying more cautious data cleaning strategies, while improving the analysis pipeline and source modelling, we measure the upper limit of the power spectrum at three redshift bands centred at 6.5, 6.8, and 7.1, with the lowest measurement obtained at z = 6.5 of |$\Delta ^2 \le (73.78\, \mathrm{mK})^2$| at |$k=0.13\, \mathrm{h\, Mpc^{-1}}$|⁠. The upper limits are the result of integrating |${\sim}14\, \mathrm{hr}$| of data, which despite shorter integration time is 1.26 times lower (at |$k=0.14 \mathrm{h\, Mpc^{-1}}$|⁠) than the previous EoR1 analysis results by Trott et al. (2020) which is a |${\sim}19\, \mathrm{hr}$| integration. This improvement confirms the necessity of improving the modelling of bright and extended sources with significant spectral structure such as Fornax-A to achieve more robust foreground mitigation, in addition to removing potentially contaminated data via visibility rms metrics. Although the cuts we used are motivated by removing systematics, they can also lead to signal loss if the data is not systematically limited. In future, we need to expand our signal loss simulation to a more representative sample to show their effect on the cosmological signal.

We also studied the difference between two MWA EoR fields, EoR0 and EoR1, each characterized by different foregrounds which changes the sky temperature. EoR0 is more contaminated by galactic emission on large scales, while EoR1 includes the bright radio source of Fornax-A. As shown in the comparison power spectra of the two fields in Fig. 9, while EoR1 is generally more contaminated by sources than EoR0, the strong presence of galactic emission in EoR0 along the horizon border impacts the wide-field measurements at large scale. The EoR1 analysis in this work and Trott et al. (2020) include ∼14 and |$19\, \mathrm{hour}$| integration, respectively. The improvements in the analysis tools warrant performing a deeper analysis similar to that carried out by Trott et al. (2020) for EoR0, using 100s of hours of data.

In this work, we explored various aspects of the data to reduce the final upper limit:

• Since the 2013 data does not contain information on the dead dipoles which affects the modelled beam response, additional systematics are introduced into the power measurement of any data set which includes 2013 data. Using 2015/2014 data in this work, the selected data set results in power measurements which are improved relative to previous published limits, despite the shorter integration time.

• As autocorrelations convey the same signature of error contamination as cross-correlations, they are a powerful tool for cleaning the data. In the case of MWA, they are also computationally efficient, as they can be implemented without the need for calibration.

• One of the effective approaches to mitigate systematics is to avoid all affected data once a hardware-related issue is recognized. As an example, the high-rms data studied in this work appeared to occur for a finite time duration on the scale of months for all LSTs, indicating a hardware-related issue. By flagging such data sets, we can prevent the introduction of identified systematics into the power spectrum.

Each new EoR data analysis can further illuminate the strengths and weaknesses of the existing EoR data sets and indicate new directions for the next generation of experiments. With the current development of SKA-Low and other new telescopes, these experiments can advise more effective and robust pipelines. This work again confirms how developing more precise sky models and upgrading analysis algorithms equips us for more effective mitigation of the systematic errors. However, even though many improvements have been made, we selected just |${\sim}60{{\ \rm per\ cent}}$| of data as a relatively reliable data set to obtain less contaminated measurements. The data loss was due to poor observational conditions, poor analysis outcomes or poor quality data. Particularly, large quantities of the poor quality data we excised seems to be related to hardware issues, specifically regarding a few tiles within each observation. Thus, along with the efforts toward more robust analysis algorithms and sky models, it is also essential to sustain routines in the hardware platform to lessen data failures, hence extracting the most out of observing and processing hours.

ACKNOWLEDGEMENTS

The authors wish to thank the anonymous referee for providing constructive suggestions that greatly improved the paper. This research was supported by Australian Government Research Training Program Scholarship and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

Computational resources used were awarded under Astronomy Australia Ltd’s ASTAC merit allocation scheme on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR programme receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government.

This work makes use of the Murchison Radio-astronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site.

DATA AVAILABILITY

The data used in this work is publicly available in MWA portal: https://asvo.mwatelescope.org/.

Footnotes

1

http://www.lofar.org

2

http://eor.berkeley.edu

3

http://lwa.unm.edu

4

http://www.mwatelescope.org

5

http://aoflagger.sourceforge.net

6

Each 112-s snapshot (observation) accommodates a series of 14 measurement sets, each covering 8-s integration time. We alternately label the 8-s cadences as even and odd. Therefore, there are 7 pairs of opposing even/odd cadences in each observation.

7

It refers to beamformer-to-receiver cable group.

APPENDIX A: COMPARISON OF SKY TEMPERATURE OF EOR0 AND EOR1 FIELDS

We compare the sky temperature over EoR0 and EoR1 patches by means of both the Global Sky Model (GSM; de Oliveira-Costa et al. 2008) and directly from our observational data.

In Fig. A1, the GSM prediction of the sky temperature in the EoR0 and EoR1 fields is illustrated. The GSM data has been weighted by the zenith-pointing MWA beam. Fig. A1a shows the sky temperature as a function of LST at a frequency of |$180\, \mathrm{MHz}$| which is the MWA high-band central frequency. EoR0 is observed at zenith at |$LST=0\, \mathrm{hr}$| and EoR1 at |$LST=4\, \mathrm{hr}$|⁠. The plot indicates the lower temperature of EoR1(⁠|${\sim}117\, \mathrm{K}$|⁠) compared to EoR0(⁠|${\sim}159\, \mathrm{K}$|⁠) at |$180\, \mathrm{MHz}$|⁠, while also showing the peak at |$LST=17\, \mathrm{hr}$| which corresponds to observing the Galactic Centre at MWA zenith. Fig. A1 b also illustrates the trend of sky temperature over the whole MWA frequency band for EoR0 and EoR1, confirming a theoretically colder sky at EoR1.

Figure A1.

GSM predicted sky temperature with MWA zenith-pointing beam, indicating EoR1 is a colder patch of sky comparing to EoR0.

To demonstrate these differences with our observational data, we take advantage of the temperature information embedded in the visibility noise rms. Following Thompson (1999), T_sys can be obtained as

$$\begin{eqnarray} T_{\rm sys} = \frac{A\sigma _{\rm vis,p}\sqrt{{\rm d}f {\rm d}t}}{k_{\rm B}} \end{eqnarray}$$

(A1)

in which A is the effective area per tile, σ_{vis, p} is the visibilities noise rms of a single polarization p, as defined in Section 4.1, df and dt are the frequency and time resolution of the observations and k_B is the Boltzmann constant. The exact same analysis was performed on the EoR1 field and the EoR0 field, and thus the main difference in the visibility rms should be related to the system temperature. In Fig. A2, the rms and corresponding system temperatures of both EoR1 and EoR0 are compared, giving a system temperature of |${\sim}240\, \mathrm{K}$| for EoR0 and ∼180K for EoR1 at |$180\, \mathrm{MHz}$|⁠. The observationally measured system temperature of a radio interferometer, T_sys, is the sum of the sky temperature and the receiver temperature. Thus, receiver temperatures of |${\sim}63\, \mathrm{K}$| and |${\sim}81\, \mathrm{K}$| for EoR1 and EoR0 respectively are deduced. The estimated receiver temperatures are broadly consistent with the reported temperature of |$50\, \mathrm{K}$| at |$150\, \mathrm{MHz}$| by Tingay et al. (2013) but do indicate some unidentified systematic effects which also differ between the fields. This discrepancy might be explained by the possible differences of the related components between the two fields; e.g. the basis of GSM model is the Galactic diffuse emission while not including the bright extragalactic sources such as Fornax-A (which is the dominant emission in case of EoR1). It should also be considered that the calculations of system temperature of the two fields are based on two typical individual observations which might not perfectly represent the actual system temperature of their field. The results are summarized in Table A1.

Figure A2.

System temperature (scaled on the left vertical axis) and visibility noise rms (scaled on the right vertical axis) for EoR0 (in dotted orange and solid red, respectively) and EoR1 (in dashed blue and solid aqua, respectively). EoR1 visibilities show less visibility noise which agrees with lower sky temperature over this field.

Table A1.

Sky/system temperature at MWA EoR fields.

	GSM predicted	Observationally driven
Field	Sky temperature (K)	System temperature (K) at 180 MHz.
EoR0	159	240
EoR1	117	180

APPENDIX B: RECOVERING SIMULATED EoR SIGNAL WITH RTS/CHIPS PIPELINE

In order to validate our normalization and provide a measure of confidence in our limits, we performed an end-to-end simulation to generate power spectra from a known input model. We chose a fully neutral IGM power spectrum model of the EoR from Furlanetto et al. (2006) as our input. In order to keep the simulation simple, we generated a measurement uv-plane with Gaussian noise which emulates the desired power spectrum. Then, through a process called degridding, we integrated the product of the uv-plane and the uv representation of the instrumental beam at every baseline location. This creates instrumental visibilities which are then analysed by the RTS/CHIPS pipeline to generate an output power spectrum. Calibration was not performed in this normalization validation.

Fig. B1 shows both the input power spectrum from Furlanetto et al. (2006) (purple) and the output power spectrum from simulated visibilities analysed by RTS and CHIPS (black). Given the agreement, we are confident in our overall normalization of our power spectrum limits. While there are some differences, this can be investigated in the future using image realizations of the EoR, foregrounds, and various calibration schemes.

Figure B1.

Input power spectrum of a theoretical EoR signal (purple) and the recovered output signal through the RTS/CHIPS pipeline (black). The normalization of the output power spectrum is consistent with the input, thus giving confidence to the normalization of our EoR upper limits in the absence of calibration.