Limits on the 21 cm power spectrum at $z=6.5-7.0$ from MWA observations

The study of primordial hydrogen in the early Universe provides a deeper insight into the properties of the first galaxies, the physics of mini-quasars, the formation of very metal-poor stars, and the understanding of cosmic structure formation. Hydrogen was produced during Big Bang nucleosynthesis (Kamionkowski, 2007). It remained neutral until ultraviolet radiation from early massive stars or active galactic nuclei (AGN) escaped their host galaxies. This radiation began ionizing the surrounding IGM, initiating a period known as the Epoch of Reionization (EoR; Zaroubi 2013).

The EoR represents a pivotal stage in the Universe’s evolution, and various probes are used to study this epoch. These include observations of the Gunn-Peterson trough in the spectra of distant quasars (e.g. Fan et al. 2006; Becker et al. 2015; Mortlock 2016; Zhu et al. 2022), the Lyman-$\alpha$ forest (e.g. Stark et al. 2011; Dijkstra 2014; Pentericci et al. 2014; Bosman et al. 2022), large-scale Cosmic Microwave Background (CMB) anisotropies (e.g. Planck Collaboration et al. 2020) and fluctuations in the infrared extragalactic background light (Doré et al., 2014). Observations of the 21 cm hyperfine line from neutral hydrogen are considered the most promising tool, as they directly probe the IGM (e.g. Barkana & Loeb 2001; Furlanetto et al. 2006).

Measurements of the weak 21 cm signal are hindered by bright foregrounds, which are several orders of magnitude brighter, yielding systematic contamination that includes errors due to inaccurate calibration and spectral structures from radio galaxies and Galactic emission (Jelić et al., 2008). Efforts have been invested for the past few years into enhancing the EoR data processing through characterization of the foregrounds (e.g. Datta et al. 2010; Trott et al. 2012; Vedantham et al. 2012; Chapman et al. 2016; Thyagarajan et al. 2015a, b; Eastwood et al. 2018; Mertens et al. 2018; Barry et al. 2023), the calibration model (e.g. Offringa et al. 2015; Barry et al. 2016; Trott & Wayth 2016; Patil et al. 2017; Procopio et al. 2017; Dillon et al. 2018; Kern et al. 2019; Orosz et al. 2019; Lynch et al. 2021), the instrument model (e.g. de Lera Acedo et al. 2017; Trott et al. 2017; Li et al. 2018; Nunhokee et al. 2020), the power spectrum estimation (e.g. Parsons et al. 2010, 2012; Liu et al. 2014; Choudhuri et al. 2017; Barry et al. 2019a; Kern et al. 2020), RFI (e.g. Offringa 2010; Wilensky et al. 2019; Kerrigan et al. 2019; Wilensky et al. 2023), polarization leakage (e.g. Moore et al. 2017; Nunhokee et al. 2017) and the effect of the ionosphere (e.g. Mevius et al. 2016; Jordan et al. 2017; Trott et al. 2018).

Current low-frequency experiments actively geared towards measuring the 21 cm hydrogen line through statistical fluctuations of the IGM include the Hydrogen Epoch of Reionization Array (HERA; DeBoer et al. 2017; Berkhout et al. 2024), the Low-Frequency Array (LOFAR; van Haarlem et al. 2013), and the Murchison Widefield Array (MWA; Tingay et al. 2013; Wayth et al. 2018a). These experiments operate between redshifts of 6 to 13, while other instruments, such as the Long Wavelength Array (LWA; Ellingson et al. 2009) and the New Extension in Nançay Upgrading LOFAR (NenuFAR; Munshi et al. 2024), focus on statistical measurements at higher redshifts.

All the aforementioned experiments have reported upper limits on the 21 cm power spectrum. The best 21 cm upper limits of $\Delta^{2}\leq(21.4)^{2}\,\textrm{mK}^{2}$ at $k=0.34\,h$ Mpc^-1 and $\Delta^{2}\leq(59.1)^{2}\,\textrm{mK}^{2}$ at $k=0.36\,h$ Mpc^-1 at redshifts of 7.9 and 10.4, respectively, have been reported by HERA (HERA Collaboration et al., 2023). The latest findings indicate that heating of the intergalactic medium above the adiabatic cooling limit must have occurred by at least $z=10.4$, effectively ruling out a wide range of “cold–reionization” models (Abdurashidova et al., 2022; HERA Collaboration et al., 2023).

Given that MWA EoR observations are more sensitive at lower redshifts, we focus on deriving the deepest upper limits on 21 cm power spectrum at $z=6.5$, $6.8$ and $7.0$, using a decade’s worth of data. This work aims to place new constraints on, or further extend the exclusion of “cold–reionization” models to lower redshifts.

We present the 21 cm experimental set up in Section 2 followed by the breakdown of the various stages involved in data processing in Section 3. The results are presented in Section 4 and validated in Section 5. The astrophysics inference from our results is conducted in Section 6. We discuss the findings and future work in Section 7 and conclude in Section 8.

In this work, we used observations from the MWA, a low-frequency interferometric radio telescope located at Inyarrimanha Ilgari Bundara, the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Murchison Radio-astronomy Observatory. The instrument is situated in the mid-west of Western Australia, approximately 300 kilometers inland from the coastal town of Geraldton, due to its radio-quiet nature (Offringa et al., 2015).

Data were captured using the Phase I and Phase II compact configurations across four polarizations: East–West (EW), North–South (NS), East–West/North–South (EW–NS) and North–South/East–West (NS–EW) (for detailed information, see Tingay et al. 2013; Wayth et al. 2018b). They spanned a frequency range of 167 to 197 MHz ($7.56<z<6.25$) from the target field EoR0, centered at (RA$=0h$, DEC$=-27^{\circ}$). The telescope observes using coarse pointing centers, spaced at intervals of $6.8^{\circ}$, to maintain a consistent primary beam response by aligning with fixed analogue delay settings (Tingay et al., 2013). Phase tracking remains constant on the observation field, EoR0, for all pointings.

Most observations were recorded at a resolution of 40 kHz with a 2 second cadence, although some of the observations used a finer cadence of 0.5 seconds. This resulted in raw file sizes ranging from 5 to 11 GB per 2 minute observation, with variations also attributed to later observations incorporating baselines from 145 tiles. Approximately $10\%$ of the observations were recorded at higher resolutions, producing raw file sizes of up to $247$ GB per observation. In total, 260.9 TB of raw data were processed.

This paper adopts a similar data processing and systematic mitigation strategies as Nunhokee (2024). We made a few improvements to the existing methodology which will be discussed in the following section. A schematic illustration of the processes and relevant tools used throughout the processing pipeline is presented in Figure 1. The pipeline is implemented in Nextflow.¹¹1https://github.com/MWATelescope/MWAEoR-Pipeline

This section explores the results obtained from the implementation of various statistical metrics and presents the power spectrum results from the final integration set.

Validation of the reported results is crucial to demonstrate their robustness. We take two approaches to the validation, both of which support different components of the work.

Firstly, we employ a full end-to-end simulation of the data, instrument and analysis pipeline to show that the pipeline produces 21 cm power levels that are input from the beginning, in the presence of noise, diffuse emission and compact foregrounds. This validation pipeline and its results are reported in Line et al. (2024), and define an overall data normalization correction that we need to include to restore full cosmological power (Liu et al., 2014; Barry et al., 2019a). This validation work uses a realistic 21 cm signal model (Greig et al., 2022), and matches the observational strategy used by the MWA for the data used in this work. The validation calculation is performed again with the final pipeline used in this work, to ensure that the recovered 21 cm power matched that in the test dataset.

The second approach employs a comparison between an independent analysis pipeline and our pipeline on a subset of the data. This approach has been previously used by the international MWA EoR Collaboration to validate results (Jacobs et al., 2016; Beardsley et al., 2016; Barry et al., 2019b; Li et al., 2019). Here, we use the Fast Holographic Deconvolution (FHD; Sullivan et al. 2012) and $\epsilon$ppsilon pipeline (Barry et al., 2019a) for independent analysis.

The FHD/$\epsilon$ppsilon pipeline, whilst similar to the Hyperdrive/CHIPS pipeline, uses a variety of different techniques in analysis.

• •

  Calibration in FHD/$\epsilon$ppsilon uses auto-correlations to constrain a normalised spectral amplitude, and all catalog sources down to the horizon to constrain the phase and absolute amplitude. The spectral variation of the calibration solutions are then averaged over each season and pointing (see Barry et al. 2019a).

• •

  Calibrated visibilities are gridded with a Gaussian-decomposition of the beam to avoid aliasing (Barry & Chokshi, 2022) – for simplicity, only the Gaussian recreating the FWHM of the FEE pattern was used for this validation.

• •

  Images are made at 160 kHz resolution in HEALPix projection (Górski et al., 2005) extending nearly to the horizon to avoid aliasing. These images are then coherently integrated in this basis.

• •

  The coherent image is transformed back to the Fourier basis. Maximum-likelihood estimates for the mean, observed noise, and uncertainties are propagated from the initial variance of the visibilities assuming fully independent measurements. These are then compared to an expected noise calculated from a fully independent noise simulation and propagated through the pipeline. The comparison yields near-unity, suggesting a lack of significant cross-correlations in measurements and robustness in the calculated uncertainties.

A subset of $335$ observations, chosen to be representative of the full sample, is processed through both pipelines. To keep the comparison solely between the analysis methodologies, some aspects are constrained to be the same. This includes RFI flagging, individual tile flagging, input sky catalogues, binning schemes, and Fourier transform window functions. The thermal noise from the data subset for each analysis is thus constrained to be comparable.

Figure 9 shows the measured power for redshifts z = 6.5, z = 6.8, and z = 7.0 and for each polarization from each pipeline. We set $0.015<k_{\perp}<0.06$ $h$Mpc^-1, and $k_{||}>0.11\,h$ Mpc^-1, and $k_{||}>2\times k_{\perp}+0.08\,h$Mpc^-1. This avoids the majority of foreground-contaminated modes and poorly sampled $uv$-modes while retaining as much data as possible.

The main difference between the validation results from Hyperdrive/CHIPS and FHD/$\epsilon$ppsilon is the strength of the coarse band flagging contamination (approximately $k$ = 0.043, 0.09, and 1.3 $h$Mpc^-1). The polyphase filter bank of the MWA requires flagging 40 kHz at the beginning and end of each 1.28 MHz band. Hyperdrive/CHIPS averages to 80 kHz frequency resolution throughout the post-calibration analysis. FHD/$\epsilon$ppsilon averages up to 160 kHz when creating images. Aliasing from the polyphase filter bank flagging depends on the final frequency resolution, and thus aliasing contamination is different between the two analyses. Other differences include contamination from the 150 m and 230 m cable reflection (Beardsley et al., 2016) being more prevalent in FHD/$\epsilon$ppsilon, which may reflect that the per-field averaging is accidentally propagating cable residuals.

Hyperdrive/CHIPS has significantly lower power in the East–West polarization of z = 6.8 at sensitive $k$-modes. We investigate this discrepancy by conducting several tests. FHD constructed power spectra in sequential groups of 50 observations and compared power levels generated from the best and worst sets, yielding comparable results in both Hyperdrive/CHIPS and FHD/$\epsilon$ppsilon. This result reduces the likelihood that the observed discrepancy stems from the specific $50$ faulty observations that may be treated differently by the pipelines. We then assess whether calibration influences this behavior by calibrating a subset of the observations using Hyperdrive and processing them through the power spectrum pipeline in both FHD/$\epsilon$ppsilon and CHIPS. The resulting power levels remain consistent, ruling out calibration errors as the source of the discrepancy.

Notably, this behavior is only prominent at $z=6.8$, corresponding to a center frequency of 184.5 MHz, which aligns with the broadband digital TV Channel 7 (DTV7). This suggests that low-level residual radio frequency interference (RFI), undetected by SSINS, may be responsible. To investigate further, we process the data using EAVILS, a complementary flagging algorithm currently in development. EAVILS identifies contamination within the DTV7 band, confirming the presence of residual RFI. Building on this RFI contamination analysis, we conduct simulations using WODEN (Line, 2022) to evaluate how both pipelines manage RFI contamination and to verify that CHIPS does not introduce signal decoherence when integrating datasets. Specifically, we simulate the following scenarios:

1. 1.

   We simulate five point sources within the main lobe, with flux densities ranging from 0.5 Jy to several Janskys, across 69 observations spanning two hours.

2. 2.

   We introduce a 7 Jy DTV7-like RFI, at the southern horizon during 17 randomly distributed observations throughout the night. If no mitigation is in place, this should lead to power spectrum contamination similar to that simulated in Wilensky et al. (2020).

We then generate power spectra for three cases: (i) the 69 RFI-free observations, (ii) the 17 RFI-contaminated observations, and (iii) the full dataset containing both RFI-free and RFI-contaminated observations. Our results show that when processing RFI-contaminated visibilities, FHD/$\epsilon$ppsilon produces power levels approximately an order of magnitude higher than CHIPS, whereas power levels remain consistent between the two pipelines for RFI-free observations. This discrepancy persists even after incorporating RFI-free observations, indicating that CHIPS more effectively suppresses horizon contamination without decohering the signal. The primary difference in power levels arises from the gridding kernel, which influences horizon power suppression in each pipeline. CHIPS uses a 2D Blackman-Harris taper as a gridding kernel; the kernel size is matched to the instrumental Fourier beam, with extensive testing undertaken to ensure no signal loss in the main lobe and correct power recovery. The taper’s response in image space is to retain power within the main lobe, but heavily suppress power outside. This is in contrast to the single Gaussian reconstructing the FWHM of the main lobe used in this implementation of FHD/$\epsilon$ppsilon (recent published limits from FHD/$\epsilon$ppsilon have used a Blackman Harris to suppress horizon emission, e.g., Barry et al. 2019b; Wilensky et al. 2023). While this difference is minor, it is enough to indicate that the discrepancy in the data could be caused by horizon RFI suppression from differing gridding kernels.

Other polarizations and redshifts have equivalent measured power levels at sensitive $k$-modes, while FHD generally performs better at high $k$. The overall agreement indicates that there are no outstanding biases present in the Hyperdrive/CHIPS pipeline that are not present in the FHD/$\epsilon$ppsilon pipeline. Given the significant analysis differences, including calibration, gridding kernels, integration space, and noise calculation methods, this is an indication of robustness.

Some of the processing for Hyperdrive/CHIPS was carried out in collaboration with industry partner Downunder Geosolutions, as their resources were better suited for the task.

To understand the implications of the 21-cm power spectrum upper limits from the MWA, we use a Bayesian framework based on 21cmFAST (Mesinger et al., 2011; Greig & Mesinger, 2018; Park et al., 2019; Murray et al., 2020) to infer the properties of the high-redshift IGM and the galaxies that drive the 21-cm signal. 21cmFAST forward-models the 21-cm power spectra based on bulk galaxy properties parameterized as functions of the host halo mass. These include the stellar-to-halo mass ratio, UV ionizing escape fraction, star formation timescales, duty cycles, and X-ray spectral energy distribution (SED), controlled by 9 free parameters. Among these parameters, the soft-band X-ray luminosity per unit star formation rate (SFR) is most constrained by current 21-cm limits (Greig et al., 2021a, b; Abdurashidova et al., 2022; HERA Collaboration et al., 2023), as it determines the strength of IGM heating during the EoR. For computational efficiency, we employ the 21cmFAST emulator, 21cmEMU (Breitman et al., 2024), that has sub-percent median accuracy among the prior range in all forward-modeled observables explored in this work. For additional details on 21cmFAST or 21cmEMU, see the referenced works.

We perform Bayesian inference using the Dynesty sampler (Speagle, 2020) under two scenarios: one that incorporates the MWA upper limits into the likelihood as an additional error functional form proposed by Abdurashidova et al. (2022) to penalize models with power higher than the observed limits, and one that does not. We consider the East–West and North–South polarizations as independent constraints. Both analyses account for additional multi-wavelength EoR probes, including:

1. 1.

   The CMB optical depth measured by Planck (Planck Collaboration et al., 2020), using a double-sided Gaussian likelihood (Qin et al., 2020);

2. 2.

   The galaxy UV luminosity function established by HST (Bouwens et al., 2015, 2016; Oesch et al., 2018), following the likelihood formalism in Park et al. (2019); and

3. 3.

   The neutral fraction ($x_{\rm HI}$) inferred from Lyman-$\alpha$ forests observed by VLT (Bosman et al., 2022; D’Odorico et al., 2023) at $z=6$, 7, and 8, with values of $x_{\rm HI}$ taken from Table 2 of Qin et al. (2024) and incorporated as Gaussian likelihoods.

These complementary probes place robust constraints on the star formation efficiency and the UV ionizing escape fraction of high-redshift galaxies, reinforcing that reionization was primarily driven by faint, low-mass galaxies. However, these multi-wavelength observations are largely insensitive to the heating properties of the IGM, highlighting the unique role of 21-cm cosmology in probing IGM thermal evolution during reionization.

Fig. 10 shows the inferred 21-cm power spectra at $z=6.5$ to 10, alongside the MWA results (this work and Trott et al. 2020) and other experiments including HERA (HERA Collaboration et al., 2023), LOFAR (Acharya et al., 2024), GMRT (Paciga et al., 2013) and PAPER (Kolopanis et al., 2019) for comparison. We see a range of models with excessive power at $z=6.5$ to 7.0 are now ruled out by the MWA limits. These models, often referred to as “cold reionization”, predict high 21-cm power spectra due to strong spatial fluctuations in the IGM ionized fraction and density which are amplified by minimal IGM heating. The exclusion of these models implies a more thermally evolved IGM than simple adiabatic cooling, with sufficient X-ray heating to suppress power at the observed redshifts. Extrapolating our results to higher redshifts where MWA constraints are weaker (Trott et al., 2020), we see the predicted power remains lower than (i.e. consistent with) most existing measurements⁶⁶6This includes the recent update from LOFAR where they further reduced the systemic noises using a variational autoencoder (see more in Acharya et al. 2024)., except for HERA⁷⁷7In other words, HERA still holds the deepest limit among all existing interferometers and, based on galaxy-evolutionary models such as 21cmFAST, provides the strongest constraints on the properties of high-redshift IGM and galaxies that drive the 21-cm signal(HERA Collaboration et al., 2023) at $z=8$.

In the inset of each panel in Fig. 10, we also present the lower limits of the predicted gas spin temperature ($T_{\rm S}$). Without the MWA constraints, our parameter space is dominated by models with little or no heating, where the gas temperature reaches as low as the adiabatic cooling limit. However, incorporating the MWA limits eliminates these low-temperature models, yielding the first evidence for a heated IGM at $z=6.5$ to 7.0. The highest posterior density (HPD) 95% confidence intervals (C.Is.) for the spin temperature are:

1. 1.

   $T_{\rm S}=2.5$ – 797.8 K at $z=6.5$;

2. 2.

   $T_{\rm S}=2.3$ – 758.6 K at $z=6.8$; and

3. 3.

   $T_{\rm S}=2.2$ – 727.3 K at $z=7.0$.

Ruling out “cold reionization” models also has implications for the sources of IGM heating, which in 21cmFAST are assumed to be high mass X-ray binaries (HMXBs) in star-forming galaxies whose soft-band X-ray luminosity per SFR ($L_{\rm X<2keV}/{\rm SFR}$) governs heating in the early universe. In Fig. 11, we present our updated constraints on $L_{\rm X<2keV}/{\rm SFR}$ derived from the MWA limits, together with previous results obtained from deep Chandra observations targeting low-redshift star forming galaxies (Mineo et al., 2012; Brorby et al., 2016; Lehmer et al., 2016, 2019; Fornasini et al., 2019, 2020; Saxena et al., 2021) and HERA Collaboration et al. (2023). For consistency, the Chandra results have been converted to the same soft-band energy interval and initial mass function as the default Salpeter (1955) in 21cmFAST. Overall, constraints from 21-cm power spectra follow the trend of an enhanced X-ray luminosity towards higher redshits, which are typically expected for the lower-metallicity environment of early galaxies (Fragos et al., 2013). However, unlike HERA Collaboration et al. (2023), the current MWA limits provide only a mild constraint, yielding $L_{\rm X<2keV}/{\rm SFR}>10^{39}{\rm erg\ s^{-1}\ M_{\odot}^{-1}yr}$ and still leaving local values permissible. This suggests that even modest X-ray emission can produce sufficient heating to explain the observed limits of the 21-cm power spectra.

The MWA results mark an important step in constraining the thermal state of the early IGM and the nature of heating during reionization. Future observations with increased sensitivity, such as those anticipated from the Square Kilometre Array (SKA), will tighten constraints on $L_{\rm X<2keV}/{\rm SFR}$ as well as other properties of the first galaxies (Park et al., 2019; Qin et al., 2021), enabling a more detailed understanding of the interplay between star formation, heating and the 21-cm signal. Combined with James Webb Space Telescope observations of early galaxies, these studies will refine our picture of galaxy formation and IGM evolution in the early Universe.

Additionally, MWA results can be improved by refining our current data processing and power spectrum estimation tools. Below are some potential identified avenues:

• •

  Van Vleck Correction: Addressing non-linearities in the correlator responses caused by signal loss during quantization (see Benkevitch et al. 2016; Star 2024 for details).

• •

  Kernel Density Estimation: Reducing contamination within the EoR window by isolating non-Gaussian components introduced by foregrounds, using a kernel density estimator.

• •

  Refinement of RFI Flagging: While our current RFI mitigation techniques have successfully detected many instances of bright, or rapidly changing RFI, there is significant room for improvement in this area. Other flagging algorithms currently under development, such as EAVILS, show great potential for detecting faint and slowly moving RFI.

• •

  Using auto-correlations for calibration: Applying bandpass solutions derived from auto-correlated visibilities can mitigate spurious spectral structures introduced into the calibration solutions by imperfections in sky and instrument models (Barry et al., 2019a, b; Li et al., 2019). This method was first implemented in Ewall-Wice et al. (2016) to reduce residual spectral structures imprinted on foregrounds by reflections in the signal chain.

• •

  Re-performing direction-independent calibration: Yoshiura et al. (2021) demonstrated that ionospheric correction is essential even during direction-independent calibration. This process involves recalibrating for direction-independent gains using an updated calibration model that incorporates ionospherically corrected sources.

In this study, we presented the data processing and power spectrum estimation pipeline implemented using the Hyperdrive/CHIPS pipeline. Our approach builds upon the systematic mitigation strategy described in Nunhokee (2024), with several modifications designed to address the challenges posed by incorporating data from Phase II MWA configurations.

Starting with an initial dataset of 19,698 observations, our methodology identified 52$\%$ of the data as high quality and suitable for power spectrum analysis. Observations from pointings at -3 were excluded due to the prominent aliasing effects of the Galactic plane.

We constructed power spectra at three redshifts: $z=6.5$, $z=6.8$ and $z=7.0$. At $z=6.5$, the most stringent upper limit was determined to be $(30.2)^{2}$ mK² at $k=0.18,h$ Mpc^-1. At $z=6.8$, the upper limit was measured as $(31.3)^{2}$ mK² at $k=0.18,h$ Mpc^-1, while at $z=7.0$, it was evaluated as $(39.1)^{2}$ mK² at $k=0.21,h$ Mpc^-1.

Our results provide the first evidence of a heated IGM at redshfits $z=6.5$ to $z=7.0$. By ruling out models with excessive power, these limits establish a lower bound of $L_{\rm X<2keV}/{\rm SFR}>10^{39}{\rm erg\ s^{-1}\ M_{\odot}^{-1}yr}$ for star-forming galaxies, where the soft-band X-ray luminosity per unit star formation rate governs the heating of the early Universe.

We validated the results produced from Hyperdrive/CHIPS against an end-to-end simulation and an independent pipeline, FHD/$\epsilon$ppsilon. The simulation, performed using WODEN⁸⁸8https://github.com/JLBLine/WODEN, derived a normalization factor that matched the observed power levels. Additionally, FHD/$\epsilon$ppsilon produced power spectra consistent with Hyperdrive/CHIPS at low $k$-modes when applied to a subset of the final integration. However, at high $k$-modes, Hyperdrive/CHIPS exhibited significantly higher power levels. These discrepancies are attributed to the 160 kHz averaging performed by FHD/$\epsilon$ppsilon. Beyond this, Hyperdrive/CHIPS does not display any unique bias.

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. This scientific work uses data obtained from Inyarrimanha Ilgari Bundara / the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. Establishment of CSIRO’s Murchison Radio-astronomy Observatory is an initiative of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund. Support for the operation of the MWA is provided by the Australian Government (NCRIS), under a contract to Curtin University administered by Astronomy Australia Limited. This work was supported by resources provided by the Pawsey Supercomputing Research Centre with funding from the Australian Government and the Government of Western Australia. This work was supported by resources awarded under Astronomy Australia Ltd’s merit allocation scheme on the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). The International Centre for Radio Astronomy Research (ICRAR) is a Joint Venture of Curtin University and The University of Western Australia, funded by the Western Australian State government. This research received technical support from the Australian SKA Regional Centre (AusSRC). The inferences were performed on the Gadi supercomputer. YQ acknowledges HPC resources from the ANUMAS and support from the ARC Discovery Early Career Researcher Award (DECRA) through fellowship DE240101129. NB acknowledges support by the ARC Discovery Early Career Researcher Award (DECRA) through project number DE240101377. JCP acknowledges support from the U.S. National Science Foundation grant #2106510. MFM, BH, and EL acknowledge support from the U.S. National Science Foundation grants #2107538 and #2228990.