A new MWA limit on the 21 cm power spectrum at redshifts ∼13–17

ABSTRACT

Observations in the lowest Murchison Widefield Array (MWA) band between 75 and 100 MHz have the potential to constrain the distribution of neutral hydrogen in the intergalactic medium at redshift ∼13–17. Using 15 h of MWA data, we analyse systematics in this band such as radio-frequency interference (RFI), ionospheric and wide field effects. By updating the position of point sources, we mitigate the direction-independent calibration error due to ionospheric offsets. Our calibration strategy is optimized for the lowest frequency bands by reducing the number of direction-dependent calibrators and taking into account radio sources within a wider field of view. We remove data polluted by systematics based on the RFI occupancy and ionospheric conditions, finally selecting 5.5 h of the cleanest data. Using these data, we obtain 2σ upper limits on the 21 cm power spectrum in the range of |$0.1~ h~{\mathrm{ Mpc}}^{-1}\lessapprox k \lessapprox 1 ~ ~h~{\mathrm{ Mpc}}^{-1}$| and at z = 14.2, 15.2, and 16.5, with the lowest limit being |$6.3\times 10^6 ~\rm mK^2$| at |$k=0.14 ~h~{\mathrm{ Mpc}}^{-1}$| and at z = 15.2 with a possibility of a few  per cent of signal loss due to direction-independent calibration.

methods: data analysis, dark ages, reionization, first stars

1 INTRODUCTION

The 21 cm hyperfine transition traces the presence of neutral hydrogen throughout the Universe. In the period between recombination, z ∼ 1100, and the end of reionization, z ∼ 6, the redshifted 21 cm line, observed by the radio telescopes at 1420/(1 + z) MHz, encodes both the density and physical state of the intergalactic medium (IGM), as well as the properties and distribution of the earliest luminous sources (Furlanetto, Oh & Briggs 2006). Thus, detection of the 21 cm line is a powerful probe of the high-z Universe such as the Cosmic Dawn, when the first stars are formed, and the Epoch of Reionization (EoR), when the UV photons emitted from early galaxies ionized the entire neutral hydrogen in the IGM.

Attempts to detect this signal have included both interferometric attempts to measure fluctuations in the 21 cm brightness temperature during the EoR [Murchison Widefield Array (MWA; Tingay et al. 2013; Beardsley et al. 2016; Wayth et al. 2018a), LOFAR (van Haarlem et al. 2013), PAPER (Parsons et al. 2010), and GMRT (Paciga et al. 2013)] and global signal experiments that attempt to measure the evolution of the sky-averaged 21 cm brightness temperature [EDGES (Monsalve et al. 2017), SARAS (Singh et al. 2018), and BIGHORNS (Sokolowski et al. 2015)].

Notable progress in 21 cm line fluctuation analysis has been made at the higher frequency range of |$\rm 100~MHz\lt \nu \lt 200~MHz$| in this decade, for example using the data from MWA Phase I (Dillon et al. 2015; Beardsley et al. 2016; Barry et al. 2019), MWA Phase II (Li et al. 2019; Patwa et al.2021), LOFAR (Patil et al. 2017), and PAPER (Kolopanis et al. 2019). Latest upper limits on the 21 cm line signal reached |$1.8\times 10^3 ~\rm mK^2$| at z = 6.5 (Trott et al. 2020) and |$5.3\times 10^3 ~\rm mK^2$| at z = 9.1 (Mertens et al. 2020). These limits are starting to constrain the astrophysics of reionization (Greig et al. 2021a, b). However, the 21 cm power spectrum at the EoR has not been detected due to various systematics. For example, one of the major difficulties is to correct for the Earth’s ionosphere (Jordan et al. 2017; de Gasperin et al. 2018) and which can bias power spectrum estimation (Trott et al. 2018). Radio frequency interference (RFI), which can bias calibration, has to be flagged (Offringa et al. 2015; Morales et al. 2019). Most importantly, the bright foregrounds are a few orders of magnitude larger than the expected 21 cm signal (Jelić et al. 2008) and become further complicated by interacting with the instrumental response (Morales et al. 2019). The power spectrum measurement can also be biased by instrumental systematics such as the beam (Joseph et al. 2020) and spectral smoothness (Trott & Wayth 2016). Incomplete sky models can introduce calibration errors which can prevent detection of the cosmological power spectrum (Barry et al. 2016; Byrne et al. 2019; Zhang et al. 2020).

Further interest in 21 cm cosmology was generated by the detection of a broad absorption feature centred at 78 MHz reported by Bowman et al. (2018). The unexpected depth (∼500 mK) and profile of the EDGES detection suggest that fiducial expectations (e.g. Pritchard & Loeb 2012) regarding physical conditions at high redshift may need to be substantially revised (Barkana 2018; Fialkov & Barkana 2019), although the cosmological origin of this feature remains uncertain (Hills et al. 2018; Singh & Subrahmanyan 2019). Interferometric (non)detection of the larger than expected spatial fluctuations in the 21 cm signal which should accompany a large global feature could confirm or constrain the EDGES result. Previous interferometric constraints on the 21 cm fluctuations prior to reionization have been reported by MWA (hereafter, EW16 Ewall-Wice et al. 2016), LOFAR (Gehlot et al. 2019), OVRO-LWA (Eastwood et al. 2019), and LOFAR-AARTFAAC (Gehlot et al. 2020) although none of these results have yet reached the sensitivity required to substantially constrain theoretical models which might predict an EDGES-like global signal. Whether or not the EDGES signal is ultimately confirmed, these investigations emphasize the unique potential of the 21 cm line to probe the z > 15 Universe and the necessity of preparing for observations with upcoming facilities such as HERA (DeBoer et al. 2017), NenuFAR (Zarka et al. 2015), and SKA1-LOW (Koopmans et al. 2014).

Unfortunately, all of the same major difficulties which complicate a detection of the EoR signal apply at lower frequencies. Bright radio foregrounds, RFI, ionospheric distortions, instrumental systematics and difficulty of calibration are not only present but invariably more prominent at lower frequencies. In this work, we report on further efforts to measure the 21 cm power spectrum at z ∼13–17 using the MWA. We analyse ultralow frequency data, 75 MHz < ν < 100 MHz, using the real-time system (RTS; Mitchell et al. 2008) and the Cosmological HI Power Spectrum estimator (CHIPS; Trott et al. 2016) for the first time. We mainly investigate the above difficulties and attempt to optimize the calibration setting to mitigate the systematics.

This paper is organized as follows: in Section 2, we describe the MWA instrument, observations, and our processing pipelines; in Section 3, we describe the effect of several important sources of systematic error and mitigation strategies that we have employed to improve our data calibration; and in Section 4, we show power spectrum results using the best data sets; and we conclude in Section 5.

2 OBSERVATIONS AND DATA ANALYSIS

2.1 The MWA

The MWA (Tingay et al. 2013; Wayth et al. 2018a) is a low-frequency radio interferometer located at the Murchison Radio Observatory (MRO) in Western Australia . The detection of the EoR is one of the primary science goals of the MWA, along with Solar and Heliospheric radio emission, Galactic and Extragalactic radio astronomy, and time domain studies (Bowman et al. 2013; Beardsley et al. 2019).

In this work, we consider data from the Phase I MWA which consisted of 128 tiles distributed semirandomly over ∼3 km with a densely concentrated central core. Each tile is a 4 × 4 grid of crossed dipoles 1 m apart on top of a 25 m² mesh ground plane. The signals from 16 dipoles in each tile are combined in an analogue beamformer which can introduce discrete delays for each dipole hence allowing the primary beam of the tile to be electronically steered to different positions on the sky. The MWA receiver system digitizes, channelizes and selects a 30.72 MHz instantaneous bandwidth (Prabu et al. 2015) which is cross-correlated and averaged in the Graphics Processing Unit (GPU) based MWA correlator (Ord et al. 2015). For these data, the visibilities are recorded at 2 s and 40 kHz resolution.

2.2 Ultralow observations

Observations of the MWA EoR target fields in the lowest MWA band were initiated by EW16. We term these observations ultralow in contrast to the main EoR low (138–168 MHz) and high (168–197 MHz) bands. For the observations considered in this work, data were taken in a split band between 75–100 MHz as well as 152–156 MHz. Using information from the higher frequency sub-band offers the possibility of improved ionospheric measurements but here we only consider processing of the contiguous 75–100 MHz data. The EoR target fields are chosen to be cold regions of the radio sky and are designated EoR0 (centred at 0 h RA, −27 deg DEC) and EoR1 (4 h, −30 deg). The discrete nature of the beamformer delays means that the MWA does not track the target field but rather updates the pointing centre approximately each 30 min and allows the target field to drift through. Each such drift scan is called a pointing. We also restrict our analysis to the central three pointings and where noted focus on the zenith pointing to reduce variations due to different primary beam shapes. Observations used in this work are summarized in Table 1.

Table 1.

MWA ultralow observations including zenith and ±1 off zenith observations. The observation performed from 2014 September 22 to 2014 October 25.

Total data volume EoR0	7.4	h
Total data volume EoR1	6.4	h
Phase centre EoR0	0h–27°
Phase centre EoR1	4h–30°
Frequency range	74.88 <ν < 100.48	MHz
Coarse band channel	1.28	MHz
Frequency resolution	40	kHz
Time per snapshot	112	s
Time resolution	2	s
MWA configuration	Phase I

Fig. 1 shows our sky models for the EoR0 and EoR1 fields. The background diffuse emission is the improved Haslam map at 408 MHz¹ (Remazeilles et al. 2015), and also we show radio sources brighter than 10 Jy at 87 MHz as example, and the beam response of MWA at 87 MHz. Due to the very wide field of view of the MWA at the ultralow frequencies, the bright ‘A-team’ sources Fornax A and Pictor A are present in both the EoR0 ane EoR1 observations. In both sets of observations these sources are far from the phase centre, however in EoR1 they are in the main lobe of the primary beam. While outside of the primary beam in EoR0, Fornax A is still the brightest source above the horizon and Pictor A is the 680th brightest source at 81 MHz, at 1 apparent Jy and close to horizon. Note that removing Pictor A from the visibilities does not improve our power spectrum results for a 10 min integration for EoR0.

Figure 1.

Sky map of EoR0 and EoR1. The EoR0 and EoR1 are centred at (RA, Dec = 0 h, −27 deg) and (RA, Dec  = 4 h, −30 deg). The contours are MWA’s beam response at 87 MHz. As an example, coloured circles indicate the radio sources brighter than 10 Jy at 87 MHz and listed in our calibration source list. The redder (bluer) colour indicates bright (faint) source. The background diffuse emission is improved Haslam map at 408 MHz. The five star marks indicate examples of direction-dependent calibrators for zenith observation.

2.3 RTS and CHIPS

2.3.1 Calibration

We processed data through an analysis pipeline designed to measure the 21 cm power spectrum from MWA EoR data. Raw visibilities from the MWA archive are processed into calibrated residual visibilities though the MWA RTS (Mitchell et al. 2008). We here describe briefly the RTS calibration procedure. Pre-processing is performed by using COTTER which performs the RFI flagging using AOFLAGGER package (Offringa et al. 2015). The data are recorded as observations which are each 112 s of data with 2 s time resolution and 40 kHz frequency resolution. Each such observation is independently processed through the RTS.

The RTS, a GPU based calibration software, applies a sky-based calibration procedure. The sky model is generated from an all sky catalogue which is created by combining multiple observations using PUMA (Line et al. 2017). The base source catalogue includes GLEAM (Wayth et al. 2015), TGSS images (Procopio et al. 2017), and MWA extended baseline data to model bright complex sources using shapelets (Line et al. 2020). The sources are ranked by beam weighted flux density, where the beam weights are evaluated for all sources at central frequency of entire bandwidth. In this work, the 1000-2000 highest ranked sources within a certain radius from the pointing centre are selected. For the highband, radius of 20 deg is chosen. However, for ultralow frequency, the beam response is considerably wider, and we empirically set a radius of 60 deg.

The initial calibration is a direction-independent (DI) calibration which combines the 1000–2000 selected sources into a single compound calibrator and solves for each complex antenna gains (Jones matrix). These solutions are performed independently for each coarse channel in parallel and are stored for later processing. We exclude baselines shorter than 20λ and taper baselines shorter than 40 λ because they are dominated by Galactic diffuse emission which is not modelled in our catalogue.

Once the DI solution is obtained, the RTS performs direction-dependent (DD) calibration. The procedure is called peeling and operated within a calibration measurement loop. In order to obtain the Jones matrix and/or to correct for ionospheric effects, the RTS reads a peeling source list which also includes 1000–2000 sources. Using the peeling source list and DI calibration result, all modelled sources are subtracted from the observed visibilities. The model subtracted visibilities are phase shifted to the position of the desired calibrator. Then, the subtracted source model of the calibrator is added to the subtracted visibilities, and the RTS solves the Jones matrix for the DD calibrator for each coarse band. Typically, the five brightest sources are selected as DD calibrators, however, based on the results of Section 3.3.2, we used zero DD calibrators for the processing presented in our final limits. For the rest of the calibrator sources, instead of full Jones matrix, a correction of ionospheric displacement and an amplitude scaling are obtained using the entire frequency band. We list our fiducial RTS setting in Table. 2.

Table 2.

Our fiducial RTS parameters used for investigating data quality. A DD calibration is performed towards the number of sources listed as DD calibration. Ionospheric phase correction is performed towards the number of sources listed as ionosphere correction and this is also the number of sources to be subtracted. Baselines less than the minimum baseline length are excluded from the calibration. We down weight baselines shorter than the tapered baseline length in the antenna gain calibration.

Number of sources for DD calibration	5
Number of sources for ionosphere correction	1000
Amplitude correction	on
Number of sources to be subtracted	1000
Beam model	Analytic
Minimum baseline length	20 λ
Tapered baseline length	40 λ

The phase error for the visibility observed by a baseline u due to ionosphere is described as exp {−iαuλ²}, where αλ² is the offset in a direction and λ is observed wavelength. During the peeling step, the ionospheric offset is measured in the direction of each source in the linear regime that αuλ² ≪ 1, and is corrected for before source subtraction. This peeling-based ionospheric calibration has been tested at higher frequencies (Intema et al. 2009; Jordan et al. 2017, hereafter, CJ17). However, at ultralow frequencies, the linear assumption might be not valid for longer baselines. For a 3 km baseline, the assumption might bias the correction for ionosphere distortion of αλ² > 1^′. The typical ionospheric median offset is approximately 0.15’ at 200 MHz for quiet ionosphere data (CJ17). The median value corresponds to approximately 0.3^′ at 85 MHz, and thus the ionosphere correction should work well for a quiet ionosphere.

Using the full Jones matrix or ionospheric calibration, the calibrator is subtracted again, and the RTS continues the peeling for the next source. This DD calibration is operated with 8 s time resolution. As for the DI calibration, we exclude and taper short baselines to avoid signal loss due to unmodelled diffuse emission. Once the DD calibration is completed for all calibrators and all time steps, the RTS outputs source subtracted residual visibilities which are used for power spectrum estimation as described below.

2.3.2 Power spectrum estimation

The cylindrical power spectrum (2D power spectrum) is described in 2D Fourier space (⁠|$\rm k_{ \parallel }$|⁠, |$\rm k_{ \perp }$|⁠), where |$\rm k_{ \perp }$| corresponds to angular scales perpendicular to the line of sight and |$\rm k_{\parallel }$| to those the line of sight. Since the foreground has smooth spectrum, the foreground contamination is effectively reduced at higher |$\rm k_{\parallel }$| modes while mode mixing caused from interferometer’s chromaticity propagates foreground contamination to higher |$\rm k_{\parallel }$| called foreground (FG) wedge. In the rest of scale at high |$\rm k_{\parallel }$| and low |$\rm k_{\perp }$|⁠, the 21 cm signal has possibility to dominate the foreground contamination and the scales are called the EoR window Datta, Bowman & Carilli (2010), Morales et al. (2012), Trott, Wayth & Tingay (2012). The spherically averaged power spectrum (1D power spectrum) at k is evaluated by averaging the 2D power spectrum within the corresponding k-bins.

The power spectrum is calculated by using the Cosmological HI Power Spectrum estimator (CHIPS; Trott et al. 2016, hereafter, CT16). CHIPS has been used in various previous works to analyse MWA data at higher frequencies (e.g. Trott et al. 2020, hereafter, CT20). We mention some recent modifications here. The inverse co-variance weighting method is replaced by an inverse variance estimator as explained in CT20. This implementation is introduced to avoid signal loss due to mis-modelling of foregrounds. Frequency modes were transformed by regular fast Fourier transformation (FFT). Although a Least Squares Spectral Analysis (LSSA) was suggested for the frequency transform in CT16, more recently CHIPS omits LSSA to avoid signal loss.

Due to aliasing from the polyphase filterbank used in the signal processing chain of MWA, 40 kHz and 80 kHz fine channels at centre and edges, respectively, of each 1.28 MHz coarse bands are always flagged. Thus, the spectral behaviour of the output visibilities has a comb structure, and this results in foreground contamination seen as spectral harmonics in the EoR window. To mitigate this effect, a Gaussian processing based kriging method was introduced in CHIPS (CT20). However, we do not use kriging in this work because the kriging does not show clear improvements on the power spectrum at large scales in CT20.

3 SYSTEMATIC ERRORS AND MITIGATION

In order to assess their quality, we examined several variable aspects of our data. Three main sources of systematic error, namely RFI, the ionosphere and the wide field of view, were evaluated. Here we describe the data quality thresholds implemented based on these evaluations, and the modifications to the standard MWA EoR calibration strategy that were used to improve our power spectrum results.

3.1 RFI

Although the MRO is an extremely radio-quiet site, RFI is unavoidable in the ultralow data, particularly in the FM band (87–108 MHz) where bright terrestrial signals can come over the horizon, be ducted through the troposphere, or be reflected off aircraft or even the Moon (McKinley et al. 2013). We investigated the RFI environment in the ultralow band by considering the flagging statistics of AOFlagger as well as SSINS (Wilensky et al. 2019), a different RFI detection algorithm developed for EoR data.

3.1.1 AOFlagger occupancy

AOFlagger returns a boolean flagging value for every time, frequency and baseline sample to indicate the suspected presence of RFI. For each MWA observation we average these flags in time to determine the flagging occupancy as a function of observing frequency. Figs 2 and 3 show the mean and median occupancy of the 2014 EoR0 and EoR1 data, respectively. Unsurprisingly, there is frequent contamination in the FM band with many channels showing persistent contamination above 10 per cent. Below the FM band, the occupancy is typically at the ∼0.5 per cent level and clearly shows evidence of the coarse band shape modulating the rate of flagging false positives. Spectral structure introduced by such systematic flagging effects may eventually limit 21 cm observations (Offringa, Mertens & Koopmans 2019). Also notable is the clearly higher ensemble occupancy of EoR1 with respect to EoR0; possibly owing to a change in the underlying visibility statistics caused by the brightness of Fornax A. In future, fine-tuning of the flagging algorithm is likely warranted but here we note that trials of ignoring the AOFlagger flags below the FM band did not produce improved results indicating that these effects are not limiting our current analysis.

Figure 2.

Percentage of samples flagged by AOFlagger as a function of observational frequency for EoR0. Bottom panel shows median occupancy, middle panel shows mean occupancy, and top panel focuses on median occupancy below 89 MHz.

Figure 3.

Percentage of samples flagged by AOFlagger as a function of observational frequency for EoR1. Bottom panel shows median occupancy, middle panel shows mean occupancy, and top panel focuses on median occupancy below 89 MHz.

We further considered the flagging occupancy of individual observations so as to select data with low RFI for power spectrum measurement. Figs 4 and 5 show the mean AOFlagger occupancy averaged over the top and bottom half of our band for each observation of EoR0 and EoR1, respectively. There is no obvious correlation between the RFI below (75−87 MHz) and in the FM band (87–99 MHz). Based on these distributions we selected low RFI data as having occupancy <3 per cent in the FM band and <1 per cent below.

Figure 4.

Percentage of samples flagged by AOFlagger for each EoR0 observation. Dotted lines indicate the selection criteria described in the text. All EoR0 observations pass the cut on RFI below the FM band. Colours represent observations from the same night.

Figure 5.

Percentage of samples flagged by AOFlagger for each EoR1 observation. The colours in this do not correspond to the same night of observation as for EoR0.

3.1.2 SSINS occupancy

We also searched for RFI using SSINS. Unlike AOFlagger which operates on individual baselines, SSINS considers the statistics of an incoherent sum over baselines as a function of time and frequency. In the EoR high band (168–198 MHz), SSINS has proven to be particularly sensitive to faint broad-band intererence such as digital television (DTV). Barry et al. (2019) found that using a SSINS-based matched filter to excise observations contaminated by DTV led to a significant improvement in the measured power spectrum limit. Although there is no known equivalent of DTV in the ultralow band, we applied SSINS to our data to search for unknown unknowns. Fig. 6 compares the RFI occupancy as measured by SSINS and AOFlagger. The occupancies are highly correlated and there is no evidence of SSINS-detected RFI which has been missed by AOFlagger. Consequently, we did not add any additional SSINS-based selection to our data.

Figure 6.

Percentage of samples flagged by AOFlagger versus flagged by SSINS.

3.2 Ionosphere

3.2.1 Ionospheric conditions

Earth’s ionosphere is a partially ionized gas at altitude of 50 km to 1000 km. The ionosphere refracts incoming radio waves and shift their phase (Thompson, Moran & Swenson 2017). As a consequence of ionospheric refraction, point source subtraction can leave large residuals since the apparent positions of extra-galactic radio sources are shifted from their true position stored in the source catalogue. Since the positional offset due to the ionosphere is proportional to ν⁻², the ionospheric conditions becomes critical at ultralow frequency. de Gasperin et al. (2018) investigated the effect of the ionosphere at 60 MHz. At the low frequencies, 2nd order terms caused by Faraday rotation become non-negligible because the term is proportional to ν⁻³. The error due to Faraday rotation is less than a few per cent compared to the first order. Thus, we ignore the Faraday rotation effect since our frequencies of interest are higher than 75 MHz, but introducing a method to correct these effects might be required to achieve the detection of 21 cm line at lower frequencies.

Using the RTS-measured ionospheric offset at each source direction, we can explore the ionospheric conditions for each 2 min snapshot. CJ17 developed the cthulhu² package for investigating the ionosphere.Cthulhu defines a metric of ionospheric activity based on the median offset and the anisotropy of ionospheric structure, where the anisotropy is quantified as eigenvalue evaluated using principal component analysis (PCA) to the source offsets. This metric has already been used in previous high-band data analysis (e.g. CT20).

Fig. 7 shows ionospheric conditions at ultralow frequency. In CJ17, data within the range of median |$\lt 0.15 ~\rm arcmin$| and PCA <70 are classified as ionospherically quiet. Note that the median offset is calculated from the measured offsets in the direction of 1000 sources at 75 MHz and the median value is normalized to 200 MHz with (75 MHz/200 MHz)². For the EoR1 field, 61 per cent of data are ionospherically quiet and the figure is 60 per cent for the EoR0 field. We note, however, that there is no quantitative threshold guaranteed not to contaminate the power spectrum for the ultralow data.

Figure 7.

Ionosphere condition at ultralow frequency. The PCA value and the median offset are calculated by using cthulhu from log files of RTS with fiducial setting. Colours indicate different dates of the observation. We set the median offset of 0.15 arcmin as the threshold of ionosphere activity, indicating with vertical dashed line.

Strong ionospheric contamination can degrade the RTS calibration results. For example, as shown in Fig. 8, the amplitude of calibrated visibilities correlates with the median source offset, where the visibility amplitude is calculated as an average over all baselines. Ideally, the visibility amplitudes should not depend on ionospheric conditions. The opposite trend is found in the averaged amplitude of the DI calibration complex gains which are calculated using all tiles and all channels. These results indicate the ionospheric error biases the gains due to effective errors of the patch source model. This figure motivates removing data with median offsets of >0.15 arcmin.

Figure 8.

The visibility amplitude versus the ionosphere median offset. Different colours indicate different dates of observation. Top and bottom panels show the result of EoR0 and EoR1. The offset correlates with visibility amplitude when the offset is larger than 0.15 arcmin indicated as vertical dashed line.

3.2.2 Updating the calibration model

Fig. 8 indicates that even for the DI calibration (patch calibration), correction for the ionosphere is essential. In order to improve the calibration, we update the patch catalogue based on measured ionospheric offsets. We first measure the ionosphere offsets at each source using full band data via the peeling calibration with 64 s integration time. Next, we shift the position of each source based on the measured offset at 88 MHz. Using the updated source model, we repeat the usual patch (DI) calibration with the RTS for each obsid.

Fig. 9 shows the ratio between the median value of the DI gain result with ionosphere patch calibration and without the update. While the ratio is consistent with one for the data with quiet ionosphere, the gain tends to be underestimated for data with active ionosphere. This result indicates that due to the offset of the apparent position of the brightest sources, the patch catalogue cannot model the sky correctly, and the error of the calibration model results in biased gains.

Figure 9.

Ratio of gain amplitude with ionosphere patch calibration update and the gain without the update. We calculate median gain amplitude among all tiles at 75 MHz. The ratio shows clear correlation with median ionosphere offset. These are all zenith observations. Threshold of the ionosphere median offset is indicated as vertical dashed line.

The 2D power spectrum is useful to evaluate the effect of ionosphere correction in 2D Fourier space. In the power spectrum regime, we compare three strategies of ionospheric calibration : (i) fiducial RTS ionosphere correction with original source list, (ii) without ionosphere correction, (iii) RTS with updating the source list. In Fig. 10, we plot the logarithm of ratio of each power spectrum. Compared to fiducial correction, peeling without ionosphere correction leaves residuals at all scales as shown in left panel of Fig. 10. The right panel of Fig. 10 shows that the updating source list method slightly improves the power spectrum in the EoR window even in data with good ionospheric conditions. This result indicates that the error of DI calibration results in unwanted frequency structure in gain which propagates the foreground power into the EoR window. The effect is more prominent for data with active ionosphere, but we omit such data from deep integration to reduce systematic errors as similar to previous works (e.g. CT20).

Figure 10.

Comparison of power spectrum of different ionosphere correction strategies for the EoR0 field. Left-hand panel shows the logarithm of the ratio of no ionosphere correction to fiducial run. Blue indicates power spectrum without ionosphere correction has larger power than fiducial. Right-hand panel shows the logarithm of the ratio of ionospheric updated source catalogue to fiducial. Red indicates that power spectrum with updated source is less than fiducual.

We mention that introducing ionospheric offsets in DI calibration is a new procedure in this work. While we show the method can improve the calibration at ultralow frequency, this strategy might be not effective at high band. This is because the improvements are not prominent for quiet data even at the ultralow frequency and ionosphere offset at the high bands will be 4 times lower than the offset at ultralow frequency. This work focus on ultralow frequency, and we leave more quantitative discussion at higher band for future works.

3.3 Wide Field of view

The primary beam of the radio telescope becomes wider with ν⁻², and the MWA’s beam is tremendously large at ultralow frequencies. For example, the beam FWHM is 42 deg at 87 MHz. Thus, the wide field effects become even more prominent than the EoR high band. In this section, we present analyses related to the beam model, DD calibration, contamination from Galactic plane and point source subtraction.

3.3.1 Beam model

In the fiducial run, we use an analytic MWA beam model. While the analytic model is fairly accurate, the true beam pattern does differ from the analytic model mainly due to mutual coupling effects (Ung et al. 2020). As the MWA beam model is not well tested in the ultralow band, we performed tests comparing an analystical implementation of the beam model in the RTS (which is much faster), to the Fully Embedded Element (FEE) model (Sutinjo et al. 2015), which has been shown to be more accurate at higher frequencies (Sokolowski et al. (2017), Line et al. (2018), Chokshi et al. (2021)).

The beam model is used to create the calibration model. As the result of the DI calibration, the averaged gain amplitude of the analytic beam model is roughly 1 per cent larger than that of the FEE model. We also compare the power spectrum with the analytic beam model and the FEE beam model. While the effect of the beam model is negligible within the EoR window, the power spectrum with the analytic model is roughly 5 per cent smaller than that with the FEE model in the foreground wedge. Although it is difficult to guarantee the accuracy of the FEE model at the ultralow frequencies, the FEE model is the most reliable beam model according to recent beam modelling experiments. Thus we use the FEE model in the calibration of best data sets.

3.3.2 DD calibration

The RTS typically performs DD calibration for a few of the brightest sources. In the process of DD calibration, the Jones matrix for the DD calibrator is solved per each tile independently for each coarse channel. This process is intended to reduce the systematic errors from bright sources which can dominate the visibilities. However, for accurate DD calibration, a high signal-to-noise ratio is essential, and the number of sources to be calibrated direction dependently is a free parameter. We note that the result of the DD calibration is only used to subtract the selected DD calibrators. For the rest of sources, the result of DI calibration and ionosphere correction is applied. Further, reducing the number of DD calibrators reduces the number of degrees of freedom used in calibration and hence should only reduce any possible signal loss.

Using 25 min of integrated zenith observation for each the EoR0 and the EoR1 fields, we here compare the 2D power spectrum with five, one, and zero DD sources. The apparent flux density of the five DD sources for the zenith observation of the EoR0 (EoR1) field are in range from 23 (34) apparent Jy to 35 (737) apparent Jy at 88 MHz. These five sources are indicated with star markers in Fig. 1. We note that depending on the observation, RTS can select different DD calibrator since the beam response varies with time and pointings. In left panels of Fig. 11, we compare the 2D power spectrum with five DD sources and one DD source. An improvement is found in the EoR window between FG wedge and first coarse band harmonic. This red region means that the power spectrum with one DD source is roughly 1.5 times lower than that with 5 DD sources. This indicates that the DD calibration for the fainter 4 sources leaves additional spectral structure in the residuals. In the right panels of Fig. 11, we compare the 2D power spectrum with one DD source and without any DD sources. The power spectrum without DD calibrator shows further improvements in the EoR window especially for the EoR1 field. This indicates the spectrally non-smooth residual of Fornax A, which is the primary DD calibrator of the EoR1, dominates these scales. We will use zero DD calibrators in our final run since the power spectrum shows the lowest power in the EoR window.

Figure 11.

Comparison of power spectrum of DD calibration for five sources (DD5), one source (DD1), and zero sources (DD0). Left- and right-hand panels are the logarithm of the ratios of DD1 to DD5 and of DD0 to DD1. Top row is for EoR0 and bottom row is for EoR1. Red indicates that the power spectrum with less DD calibrators has lower power in the EoR window.

3.3.3 Galactic plane

The primary beam of the MWA is very large and the sidelobes can often produce significant response to the Galactic plane, especially at the ultralow frequencies. As the ultralow observation strategy is a drift-and-shift tracking towards the desired field, the beam response varies with pointing every 30 min. The different pointings are denoted by different numbers. The zenith observation is ‘0’, and the pointings before and after zenith are ‘-1’ and ‘ + 1’ (Beardsley et al. 2016). Here we investigate the contamination of Galactic plane and pointing dependence of the power spectrum.

Fig. 12 shows examples of the beam response superimposed on the Haslam map. For the EoR0 field, the first side lobe overlaps with the Galactic plane especially for the −1 pointing. For the EoR1 field, the Galactic plane at 75 MHz overlaps with the edge of the main lobe.

Figure 12.

Beam response at 81 MHz over the 408 MHz improved version of Haslam map for EoR0 and EoR1. The contour indicates beam responses of 0.001, 0.1, and 0.5. Solid line (0 pointing) is the zenith observation, the dashed (−1 pointing) and dot–dashed (+1 pointing) lines are the off-zenith observations.

The contamination from Galactic plane can be seen in the 2D power spectrum as shown in Fig. 13. We select 10 min of integrated observation with the quiet ionosphere for each pointing.

Figure 13.

2D power spectrum of EoR0 (top row) and EoR1 (bottom row). Compared different grid points (pointings = 0, −1, +1). Each plot includes 10 min of integrated observations. These are N–S power spectrum. For the EoR0 of −1 pointing, strong Galactic plane contamination is found at horizon limit. The right two panels are the logarithmic ratio of 0 pointing to −1 pointing and 0 pointing to  +1 pointing. Red indicates that power spectrum of 0 pointing has less power than others.

Top row panels of Fig. 13 show the power spectrum for the EoR0 field. Since the Galactic plane is just above the horizon, the difference between pointings are visible within the horizon limit. The side lobe propagates the Galactic plane into the FG wedge for −1 pointing.

The power spectrum of EoR1 field is shown in bottom row panels. The difference between pointings are not as prominent as the EoR0. The Galactic plane becomes strong for −1 pointing. This is because the primary beam overlaps the Galactic anti-centre as shown in Fig. 12.

Note that the N–S polarization has powerful Galactic plane contamination compared to the E–W polarization. This is due to the MWA’s beam responses being different for each polarization. Thus, in the Fig. 13, we only plot N–S polarization But even for the E–W polarization, the 2D power spectrum shows similar behaviour, but the contamination is not as powerful as for the N–S.

However, the ratio shows noise like behaviour in the EoR window. Therefore, we do not remove any pointing from our final integration for both the EoR fields. On the other hand, the power spectrum of −1 pointing shows ∼10 times stronger contamination from Galactic plane compared to the zenith observation especially for the EoR0 field. In later section, we will investigate the effect of Galactic plane in our final upper limit by comparing power spectrum with −1 pointing and without it.

3.3.4 Point source subtraction

The source list used for patch calibration and peeling are made from a large catalogue built using PUMA (Line et al. 2017). Most sources are described as a point source or multiple Gaussian sources. Some of the brightest extended sources, including Fornax A and Pictor A, are modelled using shapelets (Line et al. 2020). For creating the source list, the base frequency of beam model and maximum radius from the pointing centre are free parameters. Sources are ranked by the beam weighted flux density and only sources within the radius are selected. For the highband, radius of 20 deg and central frequency of entire frequency range are chosen. However, for ultralow frequency, the optimized parameter is not known, and therefore we investigate the effect of the parameters.

First we demonstrate the effect of the radius. Two source lists for radius of 30 and 60 deg are used RTS calibration separately, and the power spectrum comparison is shown in Fig. 14. The difference between power spectrum with radius of 30 and 60 shows clear effect at the FG wedge. This is not surprising because a number of sources outside of the 30 deg radius are brighter than many sources within. We also compare the source list with radius of 75 deg, but the difference with radius of 60 deg are negligible because these two sources lists are almost identical to each other.

Figure 14.

Difference plot of 2D power spectrum for EoR0 E–W with foreground subtraction within the radius of 30 deg |$P_{\rm k}(\rm 30deg)$| from the pointing centre and 60 deg |$P_{\rm k}(\rm 60\,deg)$|⁠. Blue indicates |$P_{\rm k}(\rm 30\,deg)$| has large power than |$P_{\rm k}(\rm 60\,deg)$|⁠.

Even for the EoR0 field, sources listed in the peeling model with radius 90 deg can pick up Pictor A. Although Pictor A is far from the phase centre, it is bright at 80 MHz. However, even if we subtract the Pictor A, we cannot find any improvement for 20 min integration. The subtraction of Pictor A should be important for further calibration processing, but the source is in the edge of primary beam and the calibration has some risks due to beam correction. Thus, we do not include the Pictor A in the peeling source list for the EoR0 data.

The base frequency has no-prominent effect compared to the radius. We make source lists with radius of 60 deg and base frequency of 75, 80, 85 and 90 MHz. The difference of the power spectrum for each of these source list only shows noise-like difference. Thus, we use the central frequency of the full band as the base frequency of beam model.

Due to the incompleteness of the radio catalogue, no sources are listed in (RA, Dec)≈(-1 h, 0). At these frequencies, existing catalogues have suffered from noise and ionosphere. Thus, developing a new deep catalogue will be important for future calibration especially before the upcoming SKA1-Low observation. New surveys such as LoBES (Lynch et al. in preparation) and GLEAM-X (Wayth et al. 2018b; Hurley-Walker et al. in preparation) will fill in the remaining gaps in our sky model.

4 POWER SPECTRUM RESULTS

Consolidating the results of the above discussions, we attempt to diagnose a clean data set for our final power spectrum calculation. There are 446 snapshots observed in 2014. The RTS fails to generate DI calibration result for 21 observations. This might be due to some tiles not working properly. By visually inspecting the RTS DI gain phases, we omit a further 128 snapshots because they show rapid variation in gain phase spectrum.

We then further calibrate remaining data with source lists including 2000 objects. The source lists are built with a radius of 60 deg using the beam model at 90 MHz. Regarding the source lists, we update the position of each source using the measured ionosphere offsets as in Section 3.2.2. The RTS DI calibration is performed with the FEE beam model. We use zero sources for the full DD calibration (i.e. we only use the antenna gain solution obtained in the DI calibration). The point source subtraction and the ionosphere correction within the RTS are performed for 1000 sources. We list the RTS setting for the final run in Table 3.

Table 3.

Same as the Table 2, but for the final upper limits.

Number of sources for DD calibration	0
Number of sources for ionosphere correction	1000
Amplitude correction	on
Number of sources to be peeled	1000
Beam model	FEE model
Minimum baseline length	20 λ
Tapered baseline length	40 λ

In order to avoid the contamination from both the ionosphere and RFI, we separate the data into clean and polluted data sets. As the threshold, we regard the data with either the ionosphere median offset of more than 0.15 arcmin and averaged RFI occupancy of more than 3 per cent at higher coarse bands as the polluted data.

Fig. 15 shows the 1D power spectrum of the clean data, polluted data and integrating both data. For the EoR0 field, the clean data are clearly better than the polluted data at |$k\sim 0.1~h~{\rm Mpc} ^{-1}$|⁠. This indicates the data threshold works well to improve the integration. However, the clean data has twice as many observations as the polluted data. Thus, this comparison might not be fair. But, even if we integrate the clean and polluted data, the power spectrum does not change at the scale.

Figure 15.

Comparing the measured 1D power spectrum for the clean data, polluted data, and integrating both data. We only use zenith observations and E–W polarization. The lower straight lines are the 2σ uncertainty.

On the other hand, for the EoR1 field, the threshold does not improve the power spectrum. Rather the clean data shows slightly higher power than the polluted data, although the ratio of clean data volume to polluted data is only 0.75 in the EoR1. This seems to indicate that the power spectrum of the EoR1 field is limited by other systematics rather than the ionosphere and the RFI.

We omit the systematically polluted data from our further integration because integrating polluted data does not show improvement. Furthermore, the data with high RFI occupancy can generate unwanted bias in the power and the ionospherically active data has biased gain value even after updating the source list. Therefore, omitting the polluted data should be worth to avoid unknown biases. Although the RFI and ionosphere might not be dominant systematics for the EoR1 field, we remove the polluted data from further integration for a conservative analysis.

Consequently, 2.8 and 2.4 h of the best data are available for the EoR0 and the EoR1 fields, respectively. The integrated gridded visibilities are separated into three groups with respect to frequency: |$\rm 75\!-\!88~MHz$|⁠, |$\rm 81\!-\!94~MHz$|⁠, and |$\rm 89\!-\!101~MHz$|⁠. These separated visibility boxes are used for the power spectrum estimation at each redshift.

Fig. 16 shows the 2D power spectrum of the best data sets for the E–W polarization. The EoR0 field is better than the EoR1 field in the EoR window. Since the difference is prominent at higher |$\rm k_{\parallel }$|⁠, this indicates there is additional spectrally non-smooth contamination for the EoR1 field. This may be due to contamination from residuals of subtracted bright sources, as has been found for the EoR1 field at higher frequencies (Rahimi et al., in preparation).

Figure 16.

Measured two-dimensional power spectrum of the best clean data set for E–W polarization for the EoR0 (left) and the EoR1 (right). The frequency range is 75–87 MHz. Right-hand panel is the logarithm of the ratio between the EoR0 and the EoR1. Red indicates that the power spectrum for the EoR0 is lower than the power spectrum for the EoR1.

Fig. 17 shows the measured 1D power spectrum for each frequency group. The EoR0 for E–W polarization has the lowest power spectrum at |$k=0.1~h~{\mathrm{ Mpc}}^{-1}$|⁠. Due to the Galactic plane, the power spectrum of EoR0 for N–S polarization is larger than the E–W polarization. Note that the power spectrum is reduced by a factor of 1.15 at |$k\lt 0.07 ~h~{\mathrm{ Mpc}}^{-1}$| if we omit the data of −1 pointing for the EoR0 field. The power spectrum of EoR1 is roughly 3 times larger than the EoR0 at |$k\sim 0.1 ~h~{\mathrm{ Mpc}}^{-1}$| and at z = 16.5. In Table 4, we summarize the best 2σ upper limits for each field and polarization.

Figure 17.

Measured power spectrum of the best clean data sets for each redshift. The solid, dashed, dot–dashed and dotted lines are the EoR0 for E–W polarization, EoR0 N–S, EoR1 E–W and EoR1 N–S respectively. Thick lines show 2σ thermal uncertainties.

Table 4.

2σ upper limits on the amplitude of the 21 cm signal in units of |$\rm mK^2$| for each EoR field, polarization and redshift. The result of DI calibration is used, and the RTS setting for the limit is listed in Table 3. There is a possibility of a few per cent of signal loss due to DI calibration for the EoR0 field.

k (⁠|$~h~\rm Mpc^{-1}$|⁠)	Redshift		EoR0 E–W	EoR0 N–S	EoR1 E–W	EoR1 N–S
9.7e-02	z = 16.5	(15.2 < z < 17.9)	1.5e + 07	6.0e + 07	3.4e + 07	4.0e + 07
1.3e-01	z = 16.5	(15.2 < z < 17.9)	1.3e + 07	2.7e + 07	4.4e + 07	5.0e + 07
1.0e-01	z = 15.2	(14.1 < z < 16.4)	1.0e + 07	9.9e + 07	2.0e + 07	2.7e + 07
1.4e-01	z = 15.2	(14.1 < z < 16.4)	6.3e + 06	3.0e + 07	2.1e + 07	2.4e + 07
1.0e-01	z = 14.2	(13.1 < z < 15.2)	1.0e + 07	2.5e + 07	1.7e + 07	1.3e + 07
1.4e-01	z = 14.2	(13.1< z < 15.2)	1.0e + 07	1.4e + 07	3.4e + 07	1.2e + 07

We finally address the possibility of signal loss. DD calibration can lead significant signal loss at short baselines if Galactic diffuse emission is not included in the calibration model (e.g. Patil et al. 2016). One way to mitigate the signal loss is the exclusion of short baselines from the calibration although this cut can cause an enhancement in noise. We emphasize that the signal loss due to DD calibration should not affect our final upper limits because the DD Jones matrix calibration is not used. Furthermore, as shown in Section 3.3.2, the upper limits are enhanced as increasing the number of DD calibration sources. Thus, the signal loss associated with DD calibration is trivial in our final limits. On the other hand, however, the result of DI calibration applied to our final result can depend on the threshold wavelength as short baselines less than 20 λ and 40 λ are excluded and tapered in the DI calibration. To investigate this issue, we compare the angular power spectrum by varying the threshold value using 10 snapshots. Following (e.g. Patil et al. 2016), we exclude baselines shorter than 200 λ, and then the angular power spectrum at |u| < 200 shows an enhancement of less than 10 per cent for the EoR0 field. Further signal loss may not exist because the value does not change for higher threshold values. For the EoR1 field, on the other hand, we find roughly 5 per cent reduction in the angular power spectrum with high threshold wavelength at u < 50. The reduction might indicate a bias due to the error of the calibration catalogue because significantly bright sources dominate the diffuse emission even at u = 10 in the EoR1 field. It should be mentioned that because the sensitivity decreases when increasing the minimum baseline length, the antenna gain results have larger errors and become more susceptible to the systematics on long baselines, and hence it is hard to tell if the bias is due to signal loss or other systematic errors. For more concrete analysis, an end-to-end simulation including diffuse emission is required to probe it. Building such a large simulation is beyond the scope of this paper. Thus, we only caution there are a few  per cent of signal loss on our upper limits.

5 DISCUSSION AND CONCLUSION

By analysing 15 h of MWA Phase I ultralow data, we have investigated and mitigated against several systematic errors that dominate at these lower frequencies. Our findings can be summarized as follows: (i) While only 0.5 per cent of the lower bands (75–85 MHz) are contaminated by the RFI, the FM radio occupies more than 10 per cent of data at many fine channels. (ii) Since the ionospheric effect becomes stronger with ν⁻², the ionosphere makes modelling errors in the DI calibration. Updating the calibration model reduces the systematics, but further implementation will be required to remove the error completely. (iii) Due to the wide field of view, a calibration model has to include radio sources within a radius of 60 deg from pointing centre. (iv) Because the primary beam side lobe overlaps the Galactic plane for the EoR0, the power spectrum of N–S polarization shows prominent Galactic plane contamination at the horizon limit of the foreground wedge. (v) The DD calibration leaves spectrally non-smooth residuals after the point source subtraction as shown in Section 3.3.2 and therefore we used zero DD calibrators and the residual contamination is mitigated in the EoR window.

Based on the above jackknife analysis, we optimize the calibration of the RTS and remove systematic polluted data using thresholds on the RFI occupancy and ionosphere condition. We could use roughly 35 per cent of all data as clean data sets. By calculating the two dimensional power spectrum using the clean data, we find the EoR window of EoR1 is highly contaminated compared to the EoR0. The source of difference is not revealed but the two brightest sources (Fornax A and Pictor A) might contaminate the spectral fluctuation. For the spherical averaged 1 dimensional power spectrum, the systematics dominated best upper limits of roughly |$10^7 ~ \rm mK^2$| is obtained for the EoR0 and for the E–W polarization at 13 < z < 18 with a possibility of signal loss of a few  per cent.

We mention here that the evolution of 21 cm signal along line of sight within the band width which is called the light-cone effect. Ignoring this effect can bias the cosmological power spectrum by a factor of 2 at the cosmic dawn (z ≈18) with the redshift range Δz of 2.5 (Ghara, Datta & Choudhury 2015; Greig & Mesinger 2018). The effect is important at timing of a rapid evolution which may have happened at around z = 18 as indicated by the strong absorption in Bowman et al. (2018). Thus, power spectrum analysies at these low frequencies need to address this effect.

In Section 4, we calculated the power spectrum using gridded visibility data with bandwidth of 12.8 MHz. This frequency range is wider than that in Eastwood et al. (2019), they used only 2.6 MHz. Since we apply Blackman–Nuttall window (e.g. equation (17) in Thyagarajan et al. 2013), the effective bandwidth becomes roughly 4.65 MHz. Due to down weighting at the edges of the 12.8 MHz band, the effective redshift range is narrowed as well. While the exact redshift range is not clear due to the taper, it should be Δz ≈ 1.0 at z = 16.6.

As an experiment, we calculated power spectrum using gridded visibility with frequency range from 78.1 MHz to 84.5 MHz. The bandwidth of 6.4 MHz should correspond to Δz ≈ 0.5 at z = 16.6. Fig. 18 shows that the power spectrum with bandwidth of 6.4 MHz is 10 times larger than 12.8 MHz. Since band width becomes half, the resolution in |$\rm k_{\parallel }$| becomes worse, and the coarse band channel harmonics cannot be well separated each other. Since the enhancement seems more drastic compared to the light-cone effect, we provided the upper limits with wider bandwidth.

Figure 18.

Measured power spectrum of the best clean data sets at z = 16.5 for the EoR0 field. The solid and dashed lines are the results with bandwidths of 12.8 and 6.4 MHz, respectively. Thick lines show 2σ thermal uncertainties.

The DD effect on the direction to sources omitted from DD calibrator should be solved. When we peel such sources, ionospheric effect is corrected and the amplitude of model is just scaled by a factor over the all coarse bands before the source subtraction. Therefore the DD calibration using a smooth polynomial function over the all coarse bands might be worth. For example, such analysis is done in Mertens et al. (2020). They solved DD gain calibration with a third-order Bernstein polynomial so that the solution is spectrally smooth. They expected that non-smooth instrumental response is well corrected by the result of DI calibration and also assume that DD effects (e.g. ionosphere and beam response) are spectrally smooth. These assumption should also be valid in our case. Developing new implementation of the RTS is beyond the purpose of this paper. Thus, we left such modification for future work.

Instrumental systematics and exotic physical models regarding the Cosmic Dawn have been suggested to explain the EDGES absorption feature. Our upper limit is still ∼1–2 orders of magnitude too high to strongly constrain such models. However, the upper limit obtained in this work is an order of magnitude lower than the previous analysis EW16. This improvement encourages further processing of ultralow data. Upcoming MWA Phase III observation performed with 256 tiles simultaneously, and the MWA Phase III has higher sensitivity and higher resolution compared to the MWA Phase I. Thus the future observation at ultralow frequencies would reduce the upper limits further by providing improved radio catalogues and reducing instrumental systematics, and the result will give us important insights on the EDGES results and the Cosmic Dawn. Furthermore, as we could find systematics such as the RFI condition, gain error due to ionosphere and wide field of view, continuing the data analysis at the lowest band would reveal unknown systematics and be worthwhile for future analysis of SKA1-Low data.

DATA AVAILABILTY

The data used in this work are publicly available from the MWA archive via https://asvo.mwatelescope.org/

ACKNOWLEDGEMENTS

We would like to thank to anonymous referee for his/her useful comments. We also thank to Dr. Padmanabhan for her comments on early version of this paper. This research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. SY is supported by JSPS Overseas Research Fellowships and this work was supported in part by JSPS KAKENHI, grants: 16H05999 (SY) and 21J00416 (SY). CMT is supported by an ARC Future Fellowship through project number FT180100321. 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 program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. We acknowledge the Pawsey Supercomputing Centre which is supported by the Western Australian and Australian Governments. This scientific 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. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center.

Footnotes

1

https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_info.cfm

2

https://gitlab.com/chjordan/cthulhu

APPENDIX A: CALCULATED POWER SPECTRUM AND THERMAL NOISE UNCERTAINTY

We here summarize all power spectrum at ultralow data for comparing in future analysis. In Tables A1, A2 and A3, we show the results at z = 16.5, 15.2 and 14.2. The list contains the measured power spectrum and 2σ thermal noise for each poralization and each EoR field. As presented in Section 4, the results are estimated from 2.8 to 2.4 h of the best clean data.

Table A1.

The measured power spectrum and 2σ thermal noise uncertainty in units of |$\rm mK^2$| for each EoR field, polarization at z = 16.5. Same as Table 4, the result of DI calibration is used, and the RTS setting for the result is listed in Table 3. There is a possibility of signal loss of a few per cent due to DI calibration.

z = 16.5 (15.2<z<17.9)	EoR0				EoR1
k |$~h~\rm Mpc^{-1}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|
1.8e-02	4.5e + 08	4.5e + 08	6.9e + 02	6.9e + 02	5.7e + 08	5.7e + 08	8.7e + 02	8.7e + 02
4.0e-02	9.1e + 08	1.0e + 09	4.7e + 03	4.7e + 03	1.2e + 09	1.1e + 09	5.9e + 03	5.9e + 03
6.7e-02	1.2e + 08	2.6e + 08	2.2e + 04	2.2e + 04	9.1e + 07	9.1e + 07	2.7e + 04	2.7e + 04
9.7e-02	1.5e + 07	6.0e + 07	6.7e + 04	6.7e + 04	3.4e + 07	4.0e + 07	8.5e + 04	8.5e + 04
1.3e-01	1.2e + 07	2.7e + 07	1.2e + 05	1.2e + 05	4.4e + 07	4.9e + 07	1.5e + 05	1.5e + 05
1.7e-01	6.4e + 07	7.2e + 07	3.4e + 05	3.4e + 05	1.1e + 08	9.2e + 07	4.4e + 05	4.4e + 05
2.1e-01	2.9e + 09	2.7e + 09	6.5e + 05	6.5e + 05	4.3e + 09	4.1e + 09	8.2e + 05	8.2e + 05
2.5e-01	6.9e + 10	7.0e + 10	7.9e + 05	7.9e + 05	1.1e + 11	1.0e + 11	1.0e + 06	1.0e + 06
2.9e-01	3.8e + 10	4.0e + 10	1.3e + 06	1.3e + 06	5.7e + 10	5.6e + 10	1.6e + 06	1.6e + 06
3.4e-01	2.4e + 08	7.2e + 08	2.8e + 06	2.8e + 06	3.6e + 08	4.0e + 08	3.6e + 06	3.6e + 06
3.8e-01	6.5e + 07	1.2e + 08	2.9e + 06	2.9e + 06	2.7e + 08	2.1e + 08	3.7e + 06	3.7e + 06
4.3e-01	3.2e + 08	3.1e + 08	4.2e + 06	4.2e + 06	1.2e + 09	8.5e + 08	5.4e + 06	5.4e + 06
4.8e-01	2.6e + 10	2.6e + 10	8.3e + 06	8.3e + 06	4.7e + 10	4.4e + 10	1.1e + 07	1.1e + 07
5.4e-01	5.6e + 11	5.6e + 11	8.0e + 06	8.0e + 06	9.5e + 11	9.1e + 11	1.0e + 07	1.0e + 07
5.9e-01	2.6e + 11	2.7e + 11	1.1e + 07	1.1e + 07	4.1e + 11	3.9e + 11	1.4e + 07	1.4e + 07
6.5e-01	1.3e + 09	2.0e + 09	1.4e + 07	1.4e + 07	1.4e + 09	1.4e + 09	1.8e + 07	1.8e + 07
7.0e-01	4.7e + 08	3.3e + 08	1.8e + 07	1.8e + 07	3.8e + 09	3.7e + 09	2.3e + 07	2.3e + 07
7.6e-01	4.6e + 10	4.4e + 10	2.3e + 07	2.3e + 07	7.7e + 10	7.2e + 10	2.9e + 07	2.9e + 07
8.2e-01	1.7e + 12	1.7e + 12	2.9e + 07	2.9e + 07	2.7e + 12	2.6e + 12	3.6e + 07	3.6e + 07
8.8e-01	7.4e + 11	7.7e + 11	3.6e + 07	3.6e + 07	1.1e + 12	1.1e + 12	4.5e + 07	4.5e + 07
9.4e-01	2.5e + 09	5.6e + 09	4.4e + 07	4.4e + 07	3.3e + 09	3.3e + 09	5.5e + 07	5.5e + 07

Table A2.

Same as A1, but for z = 15.2.

z = 15.2 (14.1<z<16.4)	EoR0				EoR1
k |$~h~\rm Mpc^{-1}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|
1.9e-02	1.9e + 08	2.1e + 08	3.9e + 02	3.9e + 02	3.1e + 08	3.1e + 08	4.9e + 02	4.9e + 02
4.1e-02	4.1e + 08	5.7e + 08	2.7e + 03	2.7e + 03	8.7e + 08	7.8e + 08	3.5e + 03	3.5e + 03
6.9e-02	5.5e + 07	2.2e + 08	1.3e + 04	1.3e + 04	8.2e + 07	6.4e + 07	1.6e + 04	1.6e + 04
1.0e-01	1.0e + 07	9.9e + 07	4.0e + 04	4.0e + 04	2.0e + 07	2.7e + 07	5.1e + 04	5.1e + 04
1.4e-01	6.2e + 06	3.0e + 07	6.9e + 04	6.9e + 04	2.1e + 07	2.3e + 07	8.8e + 04	8.8e + 04
1.7e-01	2.2e + 07	3.6e + 07	2.0e + 05	2.0e + 05	6.3e + 07	4.8e + 07	2.6e + 05	2.6e + 05
2.1e-01	1.2e + 09	1.2e + 09	3.8e + 05	3.8e + 05	2.7e + 09	2.4e + 09	4.9e + 05	4.9e + 05
2.6e-01	3.3e + 10	3.4e + 10	4.7e + 05	4.7e + 05	7.4e + 10	7.0e + 10	6.0e + 05	6.0e + 05
3.0e-01	1.8e + 10	2.1e + 10	7.6e + 05	7.6e + 05	4.1e + 10	3.9e + 10	9.7e + 05	9.7e + 05
3.5e-01	1.6e + 08	4.7e + 08	1.7e + 06	1.7e + 06	2.4e + 08	2.3e + 08	2.1e + 06	2.1e + 06
4.0e-01	4.4e + 07	8.9e + 07	1.8e + 06	1.8e + 06	5.0e + 08	3.7e + 08	2.2e + 06	2.2e + 06
4.5e-01	1.2e + 08	2.4e + 08	2.5e + 06	2.5e + 06	1.3e + 09	1.0e + 09	3.2e + 06	3.2e + 06
5.0e-01	1.1e + 10	1.2e + 10	5.0e + 06	5.0e + 06	3.0e + 10	2.5e + 10	6.3e + 06	6.3e + 06
5.6e-01	2.6e + 11	2.7e + 11	4.8e + 06	4.8e + 06	6.6e + 11	6.1e + 11	6.1e + 06	6.1e + 06
6.1e-01	1.2e + 11	1.4e + 11	6.4e + 06	6.4e + 06	3.1e + 11	2.8e + 11	8.1e + 06	8.1e + 06
6.7e-01	4.6e + 08	1.9e + 09	8.3e + 06	8.3e + 06	1.3e + 09	1.1e + 09	1.1e + 07	1.1e + 07
7.3e-01	2.6e + 08	4.3e + 08	1.1e + 07	1.1e + 07	2.6e + 09	2.8e + 09	1.4e + 07	1.4e + 07
7.9e-01	2.0e + 10	2.0e + 10	1.4e + 07	1.4e + 07	5.2e + 10	4.3e + 10	1.7e + 07	1.7e + 07
8.5e-01	8.1e + 11	8.4e + 11	1.7e + 07	1.7e + 07	1.9e + 12	1.8e + 12	2.2e + 07	2.2e + 07
9.1e-01	3.6e + 11	4.0e + 11	2.1e + 07	2.1e + 07	8.0e + 11	7.4e + 11	2.7e + 07	2.7e + 07
9.8e-01	1.7e + 09	5.3e + 09	2.6e + 07	2.6e + 07	2.9e + 09	1.8e + 09	3.3e + 07	3.3e + 07

Table A3.

Same as Table A3, but for z = 14.2.

z = 14.2 (13.1<z<15.2)	EoR0				EoR1
k |$~h~\rm Mpc^{-1}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|	|$\Delta ^2_{\rm N-S}$|	|$\Delta ^2_{\rm E-W}$|	|$\rm Thermal_{\rm N-S}$|	|$\rm Thermal_{\rm E-W}$|
1.9e-02	1.4e + 08	1.4e + 08	2.3e + 02	2.3e + 02	2.3e + 08	2.1e + 08	3.0e + 02	3.0e + 02
4.2e-02	2.6e + 08	3.4e + 08	1.7e + 03	1.7e + 03	6.0e + 08	5.3e + 08	2.2e + 03	2.2e + 03
7.1e-02	4.4e + 07	9.0e + 07	8.0e + 03	8.0e + 03	6.3e + 07	3.2e + 07	1.0e + 04	1.0e + 04
1.0e-01	1.0e + 07	2.5e + 07	2.5e + 04	2.5e + 04	1.7e + 07	1.3e + 07	3.2e + 04	3.2e + 04
1.4e-01	1.0e + 07	1.4e + 07	4.3e + 04	4.3e + 04	3.4e + 07	1.2e + 07	5.6e + 04	5.6e + 04
1.8e-01	2.9e + 07	2.6e + 07	1.3e + 05	1.3e + 05	8.5e + 07	3.2e + 07	1.6e + 05	1.6e + 05
2.2e-01	8.7e + 08	7.3e + 08	2.4e + 05	2.4e + 05	1.9e + 09	1.7e + 09	3.1e + 05	3.1e + 05
2.7e-01	2.0e + 10	2.1e + 10	3.0e + 05	3.0e + 05	5.1e + 10	5.0e + 10	3.8e + 05	3.8e + 05
3.1e-01	2.1e + 10	2.3e + 10	6.8e + 05	6.8e + 05	5.2e + 10	5.0e + 10	8.7e + 05	8.7e + 05
3.6e-01	1.7e + 09	2.2e + 09	7.4e + 05	7.4e + 05	4.1e + 09	3.3e + 09	9.5e + 05	9.5e + 05
4.1e-01	2.0e + 08	1.1e + 08	1.1e + 06	1.1e + 06	1.2e + 09	4.8e + 08	1.4e + 06	1.4e + 06
4.7e-01	3.7e + 08	1.8e + 08	1.6e + 06	1.6e + 06	1.4e + 09	9.0e + 08	2.0e + 06	2.0e + 06
5.2e-01	7.7e + 09	6.7e + 09	3.1e + 06	3.1e + 06	2.0e + 10	1.9e + 10	4.0e + 06	4.0e + 06
5.8e-01	1.6e + 11	1.7e + 11	3.0e + 06	3.0e + 06	4.5e + 11	4.4e + 11	3.9e + 06	3.9e + 06
6.3e-01	7.7e + 10	8.7e + 10	4.0e + 06	4.0e + 06	2.1e + 11	2.0e + 11	5.1e + 06	5.1e + 06
6.9e-01	7.9e + 08	1.3e + 09	5.2e + 06	5.2e + 06	2.0e + 09	1.3e + 09	6.7e + 06	6.7e + 06
7.6e-01	6.8e + 08	7.9e + 08	6.7e + 06	6.7e + 06	5.2e + 09	4.1e + 09	8.7e + 06	8.7e + 06
8.2e-01	1.4e + 10	1.2e + 10	8.6e + 06	8.6e + 06	3.5e + 10	3.0e + 10	1.1e + 07	1.1e + 07
8.8e-01	5.2e + 11	5.3e + 11	1.1e + 07	1.1e + 07	1.3e + 12	1.3e + 12	1.4e + 07	1.4e + 07
9.5e-01	2.3e + 11	2.5e + 11	1.3e + 07	1.3e + 07	5.5e + 11	5.4e + 11	1.7e + 07	1.7e + 07
1.0e + 00	4.1e + 09	3.1e + 09	1.6e + 07	1.6e + 07	3.6e + 09	2.8e + 09	2.1e + 07	2.1e + 07