2.1 AusEoRPipe overview

For a full overview of the AusEoRPipe, including various software packages, see Line et al. (Reference Line, Trott, Cook, Greig, Barry and Jordan2024). Here we include the most pertinent details.

MWA observations are undertaken in 2-min snapshots, sampled at 2 s and 40 kHz resolution. Direction-independent calibration is applied at native resolution, before the data are averaged to 8 s and 80 kHz for further analysis. During the 2-min snapshot, the phase centre is fixed in celestial coordinates, and the beam pointing centre is fixed, meaning the sky drifts a small amount through the primary beam. The calibration is performed on a per-channel basis. This choice for the vanilla version of the hyperdrive software was made in response to the spectral polynomial fitting that was performed by the earlier RTS software (Mitchell et al. Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008), which imparted systematic power into the power spectrum (PS), and could not handle cable reflections. In this work, we maintain the per-channel calibration strategy and present the results within that framework.

Direction-independent calibration uses a sky model constructed from unresolved (point) and extended sources that are above the horizon at the observation time. The sky model is constructed from a combination of GLEAM (Wayth et al. Reference Wayth2015; Hurley-Walker et al. Reference Hurley-Walker2017) and LoBES (Lynch et al. Reference Lynch2021) for point and multi-component sources, and additional custom models for A-Team radio galaxies and supernova remnants (Cook, Trott, & Line Reference Cook, Trott and Line2022; Line et al. Reference Line2020). Within the EoR fields, this catalogue is 90% complete to 32 mJy (Lynch et al. Reference Lynch2021). Notably, the calibration sky model does not include diffuse emission. This omission is handled by only calibrating with baselines longer than 30 $\lambda$ , where the diffuse power is sub-dominant. Nonetheless, this is a shortcoming of the calibration sky model. The same sky model can then be used for foreground subtraction, whereby model visibilities formed from the catalogue are directly subtracted from the calibrated visibilities. This is a good path for removing foreground power, but assumes that there are no direction-dependent (DD) effects such that the model deviates from the measurements. DD calibration, such as peeling (Mitchell et al. Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008), is not implemented for this work.

2.2 Simulation strategy

We make one further concession to computational and development time and only focus on instrumental effects that are independent of direction upon the sky. These kinds of effects, such as visibility flagging and internal cable reflections, can be added to visibilities post simulation. In this way, a base set of simulations can be used to test a number of instrumental effects in isolation, for little extra compute. Sky directional effects, such as the ionosphere and directional RFI, must be simulated during the calculations of the visibilities themselves. This is outside the scope of this paper and left for future work.

2.3 Data selection

In all AusEoRPipe testing, we aim to reproduce as many real effects as possible. As such, we simulate with observational settings matching real MWA data for comparison. We select a set of 15 EoR0 zenith pointings from 2015, spanning observation GPS times 1125938488 to 1125940192 (UTC 2015-09-10 16:41:11 - UTC 2015-09-10 17:09:35). Note that GPS start time is used for the observation identifier within the MWA. These observations are known to produce good calibration solutions and have little ionospheric activity (Jordan et al. Reference Jordan2017), allowing for a reasonable comparison to simulation. We only simulate a single pointing (zenith) to reduce computational load, particularly with the diffuse and 21-cm emission. In Trott et al. (Reference Trott2020), the behaviour of individual versus combined pointings in the PS did not show significant differences, justifying use of a single pointing for this work. We simulate all data at a 2 s, 40 kHz to match the resolution of real data, ensuring any effects of averaging to 8 s, 80 kHz when applying calibration solutions or creating power spectra are captured. All power spectra shown come from the north-south aligned dipoles, because these are shown to consistently produce cleaner power spectra in previous MWA publications with this pipeline (Trott et al. Reference Trott2016). A weighted FFT is used throughout for spectral analysis.

3.1 21-cm model

We use the 21-cm model described in Paper I (Line et al. Reference Line, Trott, Cook, Greig, Barry and Jordan2024), which was designed to match EoR0 high band. It covers a $\sim 50\times50$ squared degree field, at a spectral resolution of 80kHz, and angular resolution of $\sim$ $27\,$ arcsec. The model is a TAN FITS projection with each pixel represented as a point source in the sky model, with the $\sim$ $27\,$ arcsec resolution over-sampling the $\sim$ $2\,$ arcmin MWA resolution. The visibilities are analytically predicted from a point source via assuming a point source is a dirac-delta function and applying that to the measurement equation. For the exact mechanics of the WODEN simulator, please see Line (Reference Line2022). For details of how the 21-cm model was projected and interpolated onto the sky, along with the ability of the AusEoRPipe to recover the expected PS, please see Line et al. (Reference Line, Trott, Cook, Greig, Barry and Jordan2024).

3.2 Discrete foreground model

For this model we use the current calibration catalogue as used by the AusEoRPipe. This hybrid catalogue uses GLEAM as a base (Hurley-Walker et al. Reference Hurley-Walker2017), with LoBES (Lynch et al. Reference Lynch2021) covering the EoR0 field, and the work of Procopio et al. (Reference Procopio2017) covering the EoR1 field. LoBES is used to supplant the EoR0 field as it was created with the MWA Phase II array layout, offering a resolution of $\sim$ $80\,$ arcsec, compared to the $\sim$ $2\,$ arcmin resolution of GLEAM, which was created with the MWA phase I layout. As a result, LoBES is 90% complete at 32 mJy whereas GLEAM is 90% complete at 170 mJy. This difference is visualised in the denser central fields in Fig. 1b.

A number of extended sources are modelled including Fornax A (Line Reference Line2022), and Galactic supernova remnants (Cook et al. Reference Cook, Trott and Line2022). All source spectral energy distributions are modelled by either a power-law or curved power-law model, making all discrete sources spectrally smooth. Most sources are represented as point sources in the sky model, with a number Gaussian components used for sources with more structure on the sky. The most complicated sources are represented by Shapelets, including Fornax A. Again, all visibilities are analytically calculated from the point source, Gaussian, or Shapelet parameters. The visibility calculations for a Gaussian and Shapelet component are again detailed in (Line Reference Line2022). Fig. 1b shows the sky coverage of this model; the missing area in the top right is due to the underlying GLEAM catalogue coverage. In total, the discrete sky model has 338 797 sources.

3.3 Diffuse foreground model

The diffuse sky model is based on an m-mode analysis map as described in Kriele et al. (Reference Kriele, Wayth, Bentum, Juswardy and Trott2022), made using the Engineering Development Array 2 (EDA2 Wayth et al. Reference Wayth2022). The EDA2 is a low angular resolution telescope based at the same site as the MWA, and using the same MWA dipoles, making the sky coverage perfectly matched for our purposes. The original Stokes I map was imaged at 159 MHz into a HEALPix (Górski et al. Reference Górski, Hivon, Banday, Wandelt, Hansen, Reinecke and Bartelmann2005) $N_{\mathrm{side}}=64$ projection (a HEALPixel resolution of $\sim$ $0.9^\circ$ ). Kriele et al. (Reference Kriele, Wayth, Bentum, Juswardy and Trott2022) also produced an accompanying spectral index map. To match the resolution of the MWA, we take the original Stokes I and spectral index maps and upgrade to $N_{\mathrm{side}}=2\,048$ ( $\sim$ $1.7\,$ arcmin) using the healpy (Zonca et al. Reference Zonca, Singer, Lenz, Reinecke, Rosset, Hivon and Gorski2019) ud_grade function. We then smooth the map using the healpy.smoothing function with a $0.9^\circ$ kernel, to remove the edges of the original HEALPixels. It should be noted that this diffuse map includes all sky emission, meaning all discrete sources as described in 3.2 are present. The consequence being any simulation containing both sky models will result in some double counting of flux. However, for the purposes of this paper, both models combined need not perfectly match the sky; as long as they produce comparable power to that seen in real data, they can be used in pipeline validation.

Figure 2. Comparison of a real zenith pointed 2 minsnapshot to simulated data. Both real and simulated data were averaged to 8 s, 80 kHz, and the full bandwidth imaged via WSClean. The top row were both imaged with Briggs 0 weighting, where (a) shows the real data after calibration and (b) shows the simulated data with no calibration. The bottom row were both imaged with natural weighting, where (c) shows the real data after subtracting 8 000 sources and (d) shows the simulated data with the same 8 000 sources subtracted.

4.1 Direct subtraction

In this subSection we test how much discrete sky model flux must be subtracted before we can recover the 21-cm signal. We start by creating a two minute simulation containing the full discrete and 21-cm sky models as a visibility test bed. With the EoR0 field centre at zenith, there are 227 585 discrete sources above the horizon, all of which are present in the test bed visibilities. We then generate three sub sky models to subtract from the test bed, by cutting the sky model at different flux thresholds ( $10^{-1},\,10^{-2},\,10^{-3}\,$ Jy). These flux thresholds are in reference to the apparent flux, where we have weighted the discrete sky model by the primary beam to calculate the apparent flux of all sources at the central frequency of 182 MHz. See Table 1 for the resultant numbers after these cuts were applied.

We simulate these three sub sky models through WODEN and then subtract them from the test bed directly in visibility space. Resultant 1D PS are shown in Fig. 3. Before averaging from a 2D to a 1D PS, the foreground wedge is avoided by excluding any modes where $k_\parallel \lt 0.08 \, h \mathrm{Mpc}^{-1}$ , $k_\perp \gt 0.06 \,h \mathrm{Mpc}^{-1}$ , and performing a horizon line cut (discarding any modes in the wedge as marked by the solid black line as shown in Fig. 10. Fig. 3 shows that the apparent flux cut at $10^{-3}$ Jy is necessary to reliably recover the 21-cm signal.

Table 1. Cuts placed on discrete sky model and resultant number of sources and their contribution to the apparent flux. We include a column detailing the total number of sources above the horizon without a flux cut for reference.

Figure 3. 1D PS depicting data from a single zenith EoR0 observation where only the discrete foregrounds were included are shown. All PS were obtained after performing a wedge cut.

Note this is a flux cut on this particular discrete sky model; this is not saying that leaving all sources below $10^{-3}$ Jy in real data will allow a detection. This result is limited by the completeness of this discrete sky model, and these simulations do not include confusion noise. It should also be noted that as we cut by apparent flux, at the lower the flux cut thresholds, we are adding more sources outside the main lobe of the MWA primary beam. These sources that are further from field centre contribute to higher $k_\parallel$ values due to projection effects. So although they contribute less overall power than sources in the main primary beam lobe, they contribute power closer to the so-called 2D PS ‘window’, an area of the PS expected to have the least contamination from foreground sources.

We repeated this experiment after adding the diffuse sky model to the test bed. We call this combination of 21-cm, discrete, and diffuse sky models as the full sky model. We find that subtracting the three sub discrete sky models makes little difference to the results and instead subtract the entire discrete sky model. Results are shown in Fig. 4. As mentioned in Section 3.3, the diffuse sky model contains the entire visible radio sky and so may contain duplicate flux at higher k-modes (small angular scales). Residual power left at high k-modes after subtraction therefore may be false power. However, the large residual power seen at lower k-modes comes exclusively from the diffuse model. These results show for this single observation, the 21-cm signal is unrecoverable at low k-modes without some removal of power from the diffuse foregrounds.

Figure 4. 1D PS depicting data from a single zenith EoR0 observation where both discrete and diffuse foregrounds were included.

4.2 Calibration and subtraction

We repeat the subtraction experiment in Section 4.1 but now include calibration. We use hyperdrive di-calibrate to calibrate the full sky model visibilities, using a calibration catalogue of the apparent brightest 8 000 sources in the discrete sky model. We then use hyperdrive solutions-apply to apply those calibration solutions to the full sky model visibilities. We then subtract the uncalibrated discrete sky model visibilities from the calibrated visibilities.Footnote ^a To summarise, we calibrate using an incomplete discrete sky model and then subtract the entire discrete sky model from the visibilities, thereby including a systematic error in the calibration model. We also repeat this procedure using a combination of just the discrete and 21-cm sky models. If we had perfect calibration, we would expect this to result in a perfect recovery of the 21-cm signal. Results are shown in Fig. 5. These results show that calibration induces up to three orders of magnitude of power into the window, when using a single observation. This comes from percent-level amplitude fluctuations in the calibration solutions, that vary quickly as a function of frequency. This couples power at low k-modes from the spectrally smooth foregrounds, into higher k-modes. This leakage completely masks the 21-cm signal. We note here that applying the calibration solutions to the 21-cm simulation alone does not bias the recovery of the signal. Only a small fraction of power in the sky is perturbed, and so the 21-cm alone is greater than the leakage caused by calibration solutions from itself.

Figure 5. 1D PS depicting data from a single zenith EoR0 observation, showing the effects of calibration on leakage into the window. Note that the dark purple and brown lines showing uncalibrated simulations are colour-coded to lines showing the same simulations in Figs. 3 and 4, for easy comparison.

5.1 Edge and centre channel flagging

The legacy MWA correlator used a polyphase filterbank (Tingay et al. Reference Tingay2013) to frequency channelise data. The entire bandwidth was split across 24 ‘coarse bands’, each of 1.28 MHz bandwidth, with a frequency-dependent bandpass imparted by the filterbank. These coarse bands suffered from low SNR at the edges, and aliasing causes the central channel to also degrade. For this reason, calibration is more effective when flagging these channels. Given we are testing 40 kHz resolution data, the first systematic instrumental effect is to simply flag the first two, central, and final two channels in each coarse band.Footnote ^c These regular spectral flags cause harmonic modes, linking foreground power from low to high $k_\parallel$ . The current version of the CHIPS PS processing software can optionally include a non-uniform FFT calculator to handle these missing channels. Besides flagging for bandpass effects, real data are also flagged for malfunctioning receiving elements and RFI. We leave investigating these flagging effects for future work.

5.2 Tile based gain error

Each receiving element (often called a ‘tile’ for MWA,) in a dual polarisation interferometer can have a gain ( $g_x, g_y$ ) and leakage ( $D_x,D_y$ ) term for each polarisation. Each visibility is a correlation of two tiles. WODEN implements any tile gains and leakages by multiplying the visibility between tiles 1 and 2 through the following operation:

(1) \begin{align} \begin{bmatrix} V^{\prime}_{12\,XX} & V^{\prime}_{12\,XY} \\ V^{\prime}_{12\,YX} & V^{\prime}_{12\,YY} \end{bmatrix} = \begin{bmatrix} g_{x1} & g_{x1}D_{x1} \\ g_{y1}D_{y1} & g_{y1} \end{bmatrix} \begin{bmatrix} V_{12\,XX} & V_{12\,XY} \\ V_{12\,YX} & V_{12\,YY} \end{bmatrix} \begin{bmatrix} g_{x2}^{\ast} & g_{x2}D_{x2}^{\ast} \\ g_{y2}D_{y2}^{\ast} & g_{y2}^{\ast} \end{bmatrix}^T\end{align}

where V is a visibility, $\ast$ means the complex conjugate, T is the transpose, and $V^{\prime}$ is the visibility after the tile gains and leakages have been applied. The leakage terms are calculated via Equation A4.5 from Thompson, Moran, & Swenson (Reference Thompson and Moran2017):

(2) \begin{align}D_x = \Psi - i \chi \\[-30pt] \nonumber\end{align}

(3) \begin{align}D_y = -\Psi + i \chi,\end{align}

where $\Psi, \chi$ are alignment errors of the dipoles. This equation is really designed for single antennas, but in the MWA case, one could assume all dipoles in a tile are aligned perfectly to the mesh, and the mesh is slightly offset. This would mean the alignment errors for all dipoles are the same.

In this work, we apply a different random gain error to each tile and polarisation. We draw the amplitudes of $g_x, g_y$ from a uniform distribution between 0.7 and 1.3, $U(0.7,1.3)$ . We add a phase slope to $g_x, g_y$ as a function of frequency with a maximum phase offset drawn from $U(\!-60^{\circ},+60^{\circ})$ . We draw $\Psi$ from $U(0,0.02^{\circ})$ and $\chi$ from $U(0,0.05^{\circ})$ . The random gain errors are kept constant across all 15 observations and are therefore a systematic error. The manifestation of the gain errors can be seen in Fig. 11.

5.3 Cable reflections

Mismatched impedance at coaxial cable ends can cause internal reflections that setup standing waves, adding frequency-dependent ripples to the visibilities. We follow the formalism from Beardsley et al. (Reference Beardsley2016) to define the cable reflection gain seen by tile i for polarisation pol as

(4) \begin{equation} R_{\mathrm{pol},i}(\nu) = R_{0,i} \exp(\!-2\pi i \nu \tau_i),\end{equation}

where $R_{0,i}$ is the complex reflection coefficient, and $\tau_i$ is the time delay caused by the cable length connected to tile i. The time delay is given by

(5) \begin{equation}\tau_i = \frac{2l_i}{0.81c},\end{equation}

where $l_i$ is the length of the cable connected to tile i, and c is the speed of light. The factor 0.81 comes from the velocity factor of the cable, which we again take from Beardsley et al. (Reference Beardsley2016).

The amplitudes of these cable reflections have been measured to average between 0.02 and 0.1 using Fast Holographic Deconvolution - FHD (Barry, private communication). We therefore draw the amplitude of $R_{0,i}$ from $U(0.02, 0.1)$ and add a random phase offset of $U(0.0, 180^{\circ})$ for each tile.

Figure 6. Comparison of an integration over 15 real two-minute zenith observations to simulated data. All real and simulated data have been calibrated using 10 000 sources, with the real data calibrated using 8 000 sources. In the ratios, blue means more power in the real data, red means less power. In the differences, purple means more power in the real data, and orange less. Negative pixels in a), b), e) are shown in grey. These pixels have also been masked in the ratio and differences.

5.4 Noise

Thermal noise on cross-correlations from an interferometer are estimated from both internal receiver temperature ( $T_{\mathrm{rec}}$ ) and the sky temperature ( $T_{\mathrm{sky}}$ ), via Equation 6.50 in Thompson et al. (Reference Thompson and Moran2017):

(6) \begin{equation}\sigma_{\mathrm{cross}} = \frac{\sqrt{2}k_b(T_{\mathrm{sky}} + T_{\mathrm{rec}})}{A_{eff}\sqrt{\Delta\nu\Delta t}},\end{equation}

where $A_{eff}$ is the effective area of the tile, $\Delta\nu$ is the channel width, and $\Delta t$ is the integration time. $\sigma_{\mathrm{cross}}$ describes the standard deviation of a zero mean Gaussian noise distribution. Throughout this paper, we set $T_{\mathrm{sky}} = 228\,\mathrm{K}, T_{\mathrm{rec}} = 50\,\mathrm{K}, A_{eff} = 20.35\,\mathrm{m}^2$ . A different realisation of the noise is added to each simulated observation, inducing a random error.

5.5 Applying instrumental effects

We take each effect detailed in Section 5 and apply them to the 15 zenith observations in isolation, and all in combination. We calibrate all simulations using the apparent 10 000 brightest sources in the sky for the given LST of that simulation. We integrate all 15 observations containing all instrumental effects into 2D PS and compare to real data in Fig. 6. We compare simulations containing both the discrete and diffuse sky models, as well as only the discrete sky model. Qualitatively it can be seen from the top row that the added instrumental errors match the real data, given the matching power in coarse band harmonics, and leakage at $k_\parallel \sim 0.8$ . However, there is clearly less leakage into the window from the simulations, indicating further unmodelled systematics. As expected, Fig. 6 shows missing power at the large spatial scales and along the horizon when only simulating the discrete sky model. hyperdrive calibrates cable reflections well due to its per-channel calibration strategy.

The manifestations of each instrumental effect in isolation are shown in 1D PS in Fig. 7. In general, all instrumental effects add some systematic power, but the application of calibration exacerbates this. Source subtraction removes overall power in foreground-dominated modes, but does not help to correct instrumental errors. In most cases, random errors add overall power, whereas systematic errors (e.g. flags and cable reflections) impart structured power.

Figure 7. Comparison of different instrumental effects and their manifestation in the window. These 1D PS were made from 15 zenith observations, the 2D PS of which are shown in Fig. 6. Wedge cuts have been applied before averaging into 1D. In general, all instrumental effects add some systematic power, but the application of calibration exacerbates this. Source subtraction removes overall power in foreground-dominated modes, but does not help to correct instrumental errors.

5.6 Incomplete sky model

Figure 8. Comparison of calibrating and subtracting with various numbers of sources, with simulated data on left, real data on the right, when instrumental errors are included. Power spectra are made from integrating 15 EoR0 zenith observations. The simulations are noiseless as the effects of changing the calibration are below the noise threshold. Calibration was run on simulations including noise, so noise effects are carried into calibration solutions and then applied to a noiseless simulation.

In Section 4.1, completeness of the sky model used for subtraction revealed that discrete sources above an apparent flux density of 10 $^{-3}$ Jy needed to be removed to detect the 21-cm signal. In addition, Section 4.2, which used a set calibration sky model, showed that calibration alone can introduce significant power. Here we test whether, in the presence of instrumental effects, does increasing the number of calibrator sources improve leakage In this test, all dipoles are assumed to be functional (hyperdrive only has to calculate one primary beam, and so it can do 15 000 calibration sources in about 15 min. With real data, we flag dead dipoles and have to calculate numerous primary beam patterns, which is slower. So in these tests on real data, I have assumed all dipoles are alive. This might introduce a different calibration systematic.)

A full discrete + diffuse sky model is simulated, with 5 000, 10 000, or 15 000 sources used for calibration and subtraction. A matched set of 15 zenith EoR0 observations are then calibrated/subtracted with the same parameter sets. Fig. 8 displays the results as 1D power spectra. The simulations are noiseless as the effects of changing the calibration are below the noise threshold. Calibration was run on simulations including noise, so noise effects are carried into calibration solutions and then applied to a noiseless simulation.

As observed previously, the simulated data outperform the real data, with significantly more systematic power in the latter. In the simulated datasets, inclusion of more calibration (and subtraction) sources, does reduce the leakage into the EoR window. Similar results are observed for the real data, but the relative amplitude of improvement is reduced compared with the simulations.