2.1 Murchison Widefield Array details

The MWA is a radio interferometer located at the Murchison Radio-astronomy Observatory, a radio quiet zone in Western Australia (Tingay et al. Reference Tingay2013). In this work, we are concerned with the phase II extended configuration (Wayth et al. Reference Wayth2018) of the telescope. This configuration contains 128 tiles with a maximum baseline of approximately 5.2 km. Its operating frequencies are from 70 to 300 MHz, with a field of view of approximately 26 degrees at 150 MHz. The MWA produces data at a spectral resolution of 10 kHz and a temporal resolution of 0.5 s. Each observation is typically 2 min in duration covering only a 30.72 MHz instantaneous bandwidth.

2.2 Data collection and pre-processing

Data were obtained from the MWA ASVO archiveFootnote ^a through the MWA TAP service. The TAP service was queried with the following settings: observations are within a circle of radius $5^{\circ}$ of the phase centre and centred on channel 144 (184.32 MHz, which corresponds to $z \sim 6.8$ ). With these settings, we obtain observations centred on one of two frequency bands designated for EoR. Furthermore, the settings allow us to obtain as many observations as we can without the sky changing significantly between observations. As mentioned earlier, these observations were made in the Phase II configuration, which allows us to extend upon EoR observations performed by Zheng et al. (Reference Zheng, Wu, Guo, Johnston-Hollitt, Shan, Duchesne and Li2020).

In this work, we inspect the fields used by MWA (EoR0 and EoR1, Barry et al. Reference Barry2019; Trott et al. Reference Trott2020; Rahimi et al. Reference Rahimi2021), fields which have been used by HERA (Abdurashidova et al. Reference Abdurashidova2022), and fields which were chosen by other metrics (prefixed with ‘SKAEOR’, Zheng et al. Reference Zheng, Wu, Guo, Johnston-Hollitt, Shan, Duchesne and Li2020). The phase centres and number of observations for all the fields investigated are given in Table 1. A figure highlighting the fields on the radio sky is given in Figure 1. Due to the use of analogue beamformers in the MWA, the interferometer is only capable of coarse pointing (Tingay et al. Reference Tingay2013). Thus, from the perspective of the MWA, the sky is broken into certain ‘grid numbers’ (also known as pointing numbers) which describe the pointing. For each field, numerous observations are used in order to separate field-based structure from poor observations (e.g. due to the health of the telescope itself).

Table 1. Table of fields which have been downloaded and processed, alongside the right ascension (deg) and declination (deg) of their phase centres. Listed are EoR0 and EoR1 (Lynch et al. Reference Lynch2021), fields used by HERA (Abdurashidova et al. Reference Abdurashidova2022), and fields chosen by other metrics (prefixed with ‘SKAEOR’, Zheng et al. Reference Zheng, Wu, Guo, Johnston-Hollitt, Shan, Duchesne and Li2020).

Figure 1. Fields of interest in this study highlighted on the 408 MHz all-sky map (Haslam et al. Reference Haslam, Salter, Stoffel and Wilson1982). Highlighted in green are the MWA EoR1 and EoR2 fields (Lynch et al. Reference Lynch2021). In pink are HERA fields (Abdurashidova et al. Reference Abdurashidova2022). In orange are fields selected by other metrics (Zheng et al. Reference Zheng, Wu, Guo, Johnston-Hollitt, Shan, Duchesne and Li2020).

Once the data had been downloaded, they were calibrated for direction-independent gains with HyperdriveFootnote ^b (Jordan et al. Reference Jordan2025). The sky model used is a combination of different surveys and models, including LoBES (Lynch et al. Reference Lynch2021) for the EoR 0 field, Procopio et al. (Reference Procopio2017) for the EoR1 field, Line et al. (Reference Line2020) for a shapelets model of Fornax A, (Cook, Trott, & Line Reference Cook, Trott and Line2022) for the Centaurus A model and Galactic Plane Supernova Remnants, and the GLEAM survey (Wayth et al. Reference Wayth2015; Hurley-Walker et al. Reference Hurley-Walker2016). Only Stokes I information were used for the sky model. We used the brightest 8000 sources. In addition, Hyperdrive uses a per-channel calibration at the resolution of 2 s/40 kHz. Furthermore, Hyperdrive considers the leakage terms in the Jones matrix of the beam response, utilises a simulated Full Embedded Element (FEE) beam model (Sokolowski et al. Reference Sokolowski2017) and solves for calibration solutions per frequency channel using a similar process described by Mitchell et al. (Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008). No antenna flagging algorithms were applied to the data.

Calibration solutions were assessed using the MWA’s Quality Analysis pipeline (Nunhokee et al. Reference Nunhokee, Null, Trott, Jordan and Line2024). Amplitude calibration solutions were normalised to the median, and the phase solutions were also unwrapped before calculating the metrics.

The calibrated visibilities were imaged with WSCLEAN (Offringa et al. Reference Offringa2014). Deconvolution (CLEAN, Högbom Reference Högbom1974) was applied until the data reached the $3\sigma$ noise threshold or reached the maximum number of iterations (set to 10 000). Since the field of view of the main lobe of the MWA beam is approximately 21 degrees at 184.32 MHz, we set scale of each pixel to 15 arcsec and the image size to $5\,064 \times 5\,064$ pixels.

3.1 Image metrics

Motivation: Residual signal in peeled images indicates the presence of unmodelled source sidelobe.

Expected behaviour: Root-mean-square noise should be minimised and dynamic range maximised.

In this work, deconvolved multi-frequency synthesis images are produced. Hence, residual signal in these images indicates the presence of unmodelled source sidelobes. The root-mean-square (RMS) and dynamic range (DR) are used to quantify the residual signal. The RMS should be minimised and DR maximised; a high RMS may completely obscure the 21 cm signal and may indicate significant contamination from source sidelobes, while a low DR indicates a lack of bright calibrators or, again, contamination from source sidelobes. The RMS and DR are calculated with the following equations

(1) \begin{equation} \sigma_{\text{rms}} = \sqrt{\frac{1}{N}\sum_i^N{x_i^2}} \end{equation}

and

(2) \begin{equation} \chi = \frac{X}{\sigma_{\text{rms}}} \end{equation}

respectively, where $x_i$ is the value of each pixel, N is the number of pixels, and X is the value of the brightest pixel.

Table 2. Table of metrics to be extracted from archival MWA data and their desired behaviour.

Deconvolving in the image plane is conceptually related to the process of peeling in the visibility plane: they both aim to remove unwanted sidelobes from the data. The peeling of visibilities is a common step in creating power spectra for EoR science; hence, the RMS of deconvolved images can still be suitable for determining candidate EoR fields. However, in this work, the RMS metric is intended as a simple diagnostic to compare different fields. Therefore, it is most important to ensure consistency in how these images are generated.

3.2 Amplitude solutions smoothness

Motivation: Spectral features in gain amplitudes propagate into final statistics of interest.

Expected behaviour: Smoothness metric should be minimised.

In the calibration solutions case, we investigate the amplitude and phase of the calibration solutions separately. For the amplitude solutions, we inspect the smoothness of the solutions across frequency. The smoothness is especially important since Hyperdrive, our calibration tool, makes no assumptions on how the complex gains should behave; Hyperdrive solves for solutions per-frequency, hence capturing most of the spectral features. In methods that use regularisation, assumptions about how the gains behave are made, and a penalty function is used to constrain the solutions (Ikeda et al. Reference Ikeda, Nakazato, Tsukagoshi, Takeuchi and Yamaguchi2025; Yatawatta Reference Yatawatta2015) which may bias the solutions towards a certain class of solutions.

Amplitude solutions that are not smooth – which may be the result of the particular observing field or an error at the time of observation – can contain spectral structure due to a poor sky model which will propagate into the final power spectrum (Byrne et al. Reference Byrne2019; Ewall-Wice et al. Reference Ewall-Wice, Dillon, Liu and Hewitt2017; Barry et al. Reference Barry, Hazelton, Sullivan, Morales and Pober2016). This can heavily affect the ability to make a measurement; hence, we suggest that smooth amplitude solutions are ideal. We calculate the smoothness of the amplitude solutions with the following

(3) \begin{equation} S = \frac{1}{M} \displaystyle\sum_{i = 1}^{N_+} m_i, \end{equation}

where M is the value of the DC mode of the Fourier transform (FT) of the amplitude solutions, $m_i$ is the value of the i-th positive Fourier mode (after normalising to the median), and $N_+$ is the number of positive Fourier modes. Since the amplitude solutions are not a periodic signal, we also apply the Blackman–Harris window to reduce spectral leakage (Lessard Reference Lessard2006).

3.3 Phase solutions metrics

Motivation: Phase solutions should be linear and similar between the two polarisations.

Expected behaviour: Root-mean-square error, mean absolute deviation, KS test, and average Euclidean distance should all be minimised.

Hyperdrive applies direction-independent calibration, meaning it calibrates for time delays in the system which arise from factors such as cable lengths from antenna to receiver. A time delay in the signal corresponds to a linear phase shift in Fourier space. Hence, we investigate the linearity of the phase solutions with the root-mean-square-error (RMSE) and mean absolute deviation (MAD). The RMSE and MAD are calculated using the following

(4) \begin{equation} \text{RMSE} = \sqrt{\frac{1}{N}\sum_1^N{(x_i - \hat{x}_i)^2}} \end{equation}

and

(5) \begin{equation} \text{MAD} = \frac{1}{N} \sum_1^{N} |{r_i - R}|, \end{equation}

where N is the number of frequency channels, $x_i$ is the antenna’s phase solution at frequency channel i, $\hat{x_i}$ is the phase solution predicted by the line of best fit, $r_i$ are the residuals between the observed and predicted phase solutions, and R is the mean value of the residuals.

Another intriguing aspect of the phase solutions arise from the fact that the NS and EW dipoles have differing response patterns. Hence, investigating the similarity between the NS and EW polarisation phase solutions may provide interesting insight in the behaviour of a field. Here, we suggest that solutions which are similar to each other are better than those which are not similar. A dissimilarity can indicate bright sources are present in one polarisation and not the other, or suggest an error occurred during observation in one polarisation. We use two metrics to describe the similarity, the average Euclidean distance between the two sets of solutions, and the Kolgomorov–Smirnov (KS) metric. The Euclidean distance is simply the average of the distance between the NS and EW phase solutions at each channel. A smaller value indicates that the solutions are closer. The KS metric tests how similar the underlying empirical distribution functions of the two sets of solutions are, a value of 0 indicates the two sets are identical while a value of 1 indicates they are the most dissimilar.

4.1 Gain uncertainties

For the theoretical analysis of gain uncertainties, we have used the Crameŕ–Rao bound (CRB) statistic. The CRB provides a lower bound on the variance of an unbiased estimator; or, equivalently, it provides the precision of an unbiased estimator. It does not provide a method for estimating the unknown parameters. The CRB has been used previously to investigate the effects of imprecise calibrator parameters such as the position, spectral index, and brightness (Trott et al. Reference Trott, Wayth, Macquart and Tingay2011). In this work, the same methods were applied to investigate the best a calibration pipeline can estimate gain uncertainties.

We begin with the Fisher Information Matrix, an element of which can be written in the form (Kay Reference Kay1993):

(6) \begin{equation} I_{ab}(\boldsymbol{\theta}) = \frac{2}{\sigma^2} \int{ \frac{\partial s^{H}(t;\, \boldsymbol{\theta})}{\partial\theta_a}\frac{\partial s(t;\,\boldsymbol{\theta})}{\partial\theta_b} dt}, \end{equation}

where $\sigma^2$ is the variance, $s(t;\, \boldsymbol{\theta})$ is a signal, H denotes the Hermitian conjugate, and $\boldsymbol{\theta}$ is a vector of unknown parameters which the signal is conditioned on. The signal is given by

(7) \begin{equation} V_{ab} = \sum_{j=1}^{N_c} B_j g_a \overline{g_b} \text{exp}({-}2 \pi i(u_{ab}l_j + v_{ab}m_j)), \end{equation}

where $g_a$ and $g_b$ are the complex gains of antennas a and b respectively, $u_{ab}$ and $v_{ab}$ are the baseline coordinates between antennas a and b, $B_j$ is the brightness of calibrator j, and $l_j$ and $m_j$ are the calibrator’s coordinates (direction cosines). In this study, we have assumed $g_a = g_b = 1$ , that is, unity gains.

We are concerned with how well we can estimate the gains $g_a$ and $g_b$ ; hence for a telescope with N antennas, we have $\boldsymbol{\theta} = [g_1, g_2, g_3, \dots, g_{ N }]$ . This means $s(t;\, \boldsymbol{\theta})$ in Equation (6) is a vector of visibilities measured by every baseline. Each element of this vector of length $N_{\text{baselines}}$ is given by Equation (7), where $N_{\text{baselines}} = N^2$ is the number of baselines, including conjugates, and auto-correlations to complete the matrix. Taking the partial derivative of this vector with respect to $g_a$ , we have

(8) \begin{equation} \frac{\partial s_{ab}(t;\, \boldsymbol{\theta})}{\partial g_a} = \sum_{j=1}^{N_c} B_j \overline{g_b} \text{exp}({-}2 \pi i(u_{ab}l_j + v_{ab}m_j)). \end{equation}

This means most of the vector, $\frac{\delta \boldsymbol{s}(t;\boldsymbol{\theta})}{\delta g_a}$ , will possess a value of 0 and N elements will be populated, since $N_{\text{baselines}} - N$ elements do not contain $g_a$ and will be differentiated to 0.

When $a \neq b$ , the product of the two partial derivatives is then only non-zero when the product contains $g_a$ and $g_b$ . For example, for antennas 1 and 2, the product is non-zero when the terms $g_1$ and $g_2$ are both present in the product. When $a = b$ , the product of the partial derivatives is non-zero at N elements, $N -1$ of those elements are identical, and 1 element will have a multiplicative factor of 4 due to the partial derivative of a $g_a^2$ term. Hence, when $a = b$ , an additional multiplicative factor of $N + 3$ is needed, this condition can be expressed with a Kronecker-delta function. The final expression for an element of the Fisher Information Matrix is

(9) \begin{align} &I_{ab}(\boldsymbol{\theta}) = (N + 3)^{\delta_{a b}} \times \frac{2}{\sigma^2} \times \nonumber\\ &\int{\displaystyle\sum_{j = 1}^{N_c}\displaystyle\sum_{k = 1}^{N_c}B_j B_k g_a \overline{g_b} \exp \big\{-2\pi i(u_{ab}(l_j - l_k) + v_{ab}(m_j - m_k)\big\} dt}. \end{align}

The variance $\sigma^2$ of a single visibility signal measured by a single baseline is given by

(10) \begin{equation} \sigma = \displaystyle\frac{2k_b T_{\text{sys}}}{A_{\text{eff}} \sqrt{\Delta \nu \Delta \tau}}, \end{equation}

where $k_b$ is the Boltzmann constant, $T_{\text{sys}}$ is the system temperature, $A_{\text{eff}}$ is the effective area of an antenna, $\Delta \nu$ is the bandwidth of a single measurement, and $\Delta \tau$ is the integration time for a single measurement.

The CRB matrix is then the inverse of Equation (9), and the gain uncertainties can be found on the square root of the diagonal of this matrix. Previous work with the CRB (Trott et al. Reference Trott, Wayth, Macquart and Tingay2011) has shown one bright calibrator produces better gain and phase precision compared to many lower brightness sources.

4.2 Propagation of gain uncertainties

To propagate gain uncertainties into the visibilities, we use the standard covariance matrix method given by

(11) \begin{equation} \sigma_{V_{ab}}^2 = \left[\mathbf{J} \mathbf{C}_{\theta} \mathbf{J}^\dagger \right]_{ab}, \end{equation}

where $\mathbf{J}$ is the Jacobian of partial derivatives of the visibility function with respect to the parameters $\theta_{ab} = (g_{a}, g_{b})$ , $\mathbf{J}^\dagger$ is the complex conjugate, and $\mathbf{C}_{\theta}$ is the covariance matrix of parameter uncertainties. Expanding out Equation (11), we arrive at the following

(12) \begin{align} \sigma_{V_{ab}}^2 &= \left| \frac{\partial{V_{ab}}}{\partial{g_a}} \right|^2 \sigma_{g_a}^2 + \frac{\partial{V_{ab}}}{\partial{g_a}}\overline{\frac{\partial{V_{ab}}}{\partial{g_b}}} \text{cov}_{g_bg_a}\nonumber \\ & \quad + \overline{\frac{\partial{V_{ab}}}{\partial{g_a}}}\frac{\partial{V_{ab}}}{\partial{g_b}}\text{cov}_{g_ag_b} + \left| \frac{\partial{V_{ab}}}{\partial{g_b}} \right|^2 \sigma_{g_b}^2, \end{align}

and since we assume the gains are unity, meaning the coefficients of the variance and covariance terms are the same, this expression can be simplified further to

(13) \begin{equation} \sigma_{V_{ab}}^2 = \left| \frac{\partial{V_{ab}}}{\partial{g_a}} \right|^2 \big(\sigma_{g_a}^2 + cov_{g_ag_b} + cov_{g_bg_a} + \sigma_{g_b}^2\big). \end{equation}

The variance and covariance terms can be directly taken from the diagonal and off-diagonals of the CRB matrix, respectively.

5.1 Metrics

The evaluation of both the image and most of the calibration metrics is simple. Equations (1) and (2) were evaluated after reading in the image data. The smoothness of the amplitude solutions of the antenna were found by evaluating Equation (3), after applying a Blackman–Harris window to the solutions. The phase linearity metrics were evaluated for each set of phase solutions by first applying a linear fit to the data and then evaluating Equations (4) and (5). The Euclidean distance was simply the average of the absolute value of the difference between the NS and EW phase solutions at each channel, and the KS metric was obtained using the kstest method from scipy (Virtanen et al. Reference Virtanen2020).

5.2 Detecting metric outliers

It is inevitable that some observational data are corrupted by some source of error that is independent of the sky, such as problems with the telescope or strong Radio Frequency Interference (RFI). For this work, these effects need to be disentangled from sky-based effects. Processing these bad observations will result in outliers in the metrics. These outliers were detected by first treating each observation’s metrics as a 128-dimensional point and then calculating the distance to every other observation’s metrics (another 128-dimensional point). We then calculate the MAD for these distances to use in the modified z-score. In this paper, a threshold of $5\sigma$ was used to determine if a result is an outlier. The number of observations left after applying this method for each field is given in Table 3.

Table 3. Number of 2-min observations per field remaining after bad data have been removed with the method described in Section 5.2.

5.3 Gain uncertainties

In this work, all simulations are zenith pointed meaning simulations of some of the fields may not be representative of reality, as some fields do not transit the zenith of the MWA. However, this allows us to compare the effects of source positions within the beam between the different fields. Furthermore, it is common to calibrate only with sources within the main lobe of the primary beam, because these are the most reliable. Hence, during these calculations, only sources within the main lobe are preserved. The main lobe will change slightly depending on pointing, but overall it will behave similarly whether at zenith or off-zenith. Therefore, simulating at zenith will still provide valuable insight while allowing for much easier implementation.

Additionally, all components in the sky model are treated as point sources in the simulation, as Equation (9) requires only the brightness and positions of the sources. This means that Fornax A, which is composed of shapelets in our sky model, is modelled as a set of point sources with large brightness. Moreover, systematics like sky model uncertainties are unaccounted for; all other parameters are assumed to be accurate.

Equation (10) was first evaluated for the given input parameters. In this work, the values used are representative of EoR science with the MWA telescope $T_{\text{sys}} =200\,\text{K}$ (value at the centre of the bandwidth), $N = 128$ , $A_{\text{eff}} =21\,\textrm{m}^2$ , $\Delta\nu = 80\, \text{kHz}$ , and $\Delta\tau = 8\,\text{s}$ . The field of view of a telescope can be approximated with $\theta \approx \frac{\lambda}{D}$ where D is the diameter of a dish/antenna. For the MWA, $D = 4.4\,\text{m}$ . We have assumed that $A_{\text{eff}}$ does not change with frequency (Tingay et al. Reference Tingay2013).

Next, the calculation and propagation of the gain uncertainties takes place within a large loop over the frequency range. At the beginning of each iteration, sources in the sky model were vetoed by the current field of view and also their estimated flux density for that frequency. Following this, the beam pattern is calculated for the current frequency, using the array factor method (Warnick et al. Reference Warnick, Maaskant, Ivashina, Davidson and Jeffs2018), in a 1 024 by 1 024 grid where each pixel corresponds to a portion of the sky in the (l, m) plane. These pixel values are used to attenuate sources located within the pixel by simply multiplying the source brightness by the pixel value.

The attenuated sources are then used to calculate the CRB matrix by first generating the FIM matrix with Equation (9), where each element represents an antenna pair (a, b). The FIM is Hermitian; hence, only the top triangle of the matrix needs to be calculated and the bottom triangle can be filled in by taking the complex conjugate. The CRB are the diagonal elements of the inverse of the FIM matrix, and the gain uncertainties for each antenna are their square root.

5.4 Power spectra

Taking the previously calculated gain uncertainties, still within the large frequency loop, the uncertainties for visibilities were evaluated with Equation (13). Once uncertainties were calculated, the real and complex components for the visibilities for each baseline were randomly generated from a normal distribution with mean 0 and standard deviation $\frac{\sigma_{V_{ab}}}{\sqrt{2}}$ . The visibilities were then gridded onto a common (u,v) grid and appended onto a growing list of gridded visibilities.

Once the loop completed, we were left with a data cube of gridded visibilities at each frequency channel. To transform this data cube into the final power spectrum, a FT was applied to each (u,v) cell along the frequency axis. The result was then multiplied by its complex conjugate to yield the unnormalised power. We then cylindrically average the unnormalised power at each $\eta$ slice (the FT of frequency) to arrive at the 2D temperature power spectrum.