# MeerKLASS L-band deep-field intensity maps: entering the HI dominated regime

Matilde Barberi-Squarotti, José L Bernal, Philip Bull, Stefano Camera, Isabella P Carucci, Zhaoting Chen, Steven Cunnington, Brandon N Engelbrecht, José Fonseca, Keith Grainge, et al.

► To cite this version:

Matilde Barberi-Squarotti, José L Bernal, Philip Bull, Stefano Camera, Isabella P Carucci, et al. MeerKLASS L-band deep-field intensity maps: entering the HI dominated regime. Mon.Not.Roy.Astron.Soc., 2025, 537 (4), pp.3632-3661. 10.1093/mnras/staf195. hal-04679935

HAL Id: hal-04679935

https://hal.science/hal-04679935v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MeerKLASS L-band deep-field intensity maps: entering the HI dominated regime

MeerKLASS Collaboration, Matilde Barberi-Squarotti,¹,²,³,⁴ José L. Bernal,⁵ Philip Bull,⁶,⁷

Stefano Camera,¹,⁸,⁹,⁷ Isabella P. Carucci,¹⁰,¹¹ Zhaoting Chen,¹² Steven Cunningham,¹³,⁶★

Brandon N. Engelbrecht,⁷ José Fonseca,⁷,¹³,¹⁴ Keith Grainge,⁶ Melis O. Irfan,¹⁵,⁷ Yichao Li,¹⁶,⁷

Aishrila Mazumder,⁶ Sourabh Paul,⁶ Alkistis Pourtsidou,¹²,¹⁷,⁷ Mario G. Santos,⁷,¹⁸ Marta Spinelli,¹⁹,⁷

Jingying Wang,²,⁰,⁷,★ Amadeus Witzemann,⁷,⁶ and Laura Wolz⁶

*Affiliations are listed at the end of the paper*

Accepted 2025 January 24. Received 2025 January 13; in original form 2024 August 6

## ABSTRACT

We present results from MeerKAT single-dish HI intensity maps, the final observations to be performed in L-band in the MeerKAT Large Area Synoptic Survey (MeerKLASS) campaign. The observations represent the deepest single-dish HI intensity maps to date, produced from 41 repeated scans over 236 deg², providing 62 h of observational data for each of the 64 dishes before flagging. By introducing an iterative self-calibration process, the estimated thermal noise of the reconstructed maps is limited to ~1.21 mK (1.2 × the theoretical noise level). This thermal noise will be subdominant relative to the HI fluctuations on large scales ($k \lesssim 0.15$ h Mpc⁻¹), which demands upgrades to power spectrum analysis techniques, particularly for covariance estimation. In this work, we present the improved MeerKLASS analysis pipeline, validating it on both a suite of mock simulations and a small sample of overlapping spectroscopic galaxies from the Galaxy And Mass Assembly (GAMA) survey. Despite only overlapping with ~25 per cent of the MeerKLASS deep field, and a conservative approach to covariance estimation, we still obtain a > 4 σ detection of the cross-power spectrum between the intensity maps and the 2269 galaxies at the narrow redshift range 0.39 < z < 0.46. We briefly discuss the HI autopower spectrum from these data, the detection of which will be the focus of follow-up work. For the first time with MeerKAT single-dish intensity maps, we also present evidence of HI emission from stacking the maps onto the positions of the GAMA galaxies.

**Key words:** methods: data analysis – methods: statistical – large scale structure of Universe – cosmology: observations – radio lines: general.

## 1 INTRODUCTION

Probing density fluctuations in large-scale cosmic structure is a cornerstone of precision cosmology. Constructing maps that span wide areas and reach deep redshifts provides a rich resource for testing and constraining the standard ΛCDM cosmological model. Surveys at optical and near-infrared wavelengths have made impressive progress in this pursuit (eBOSS Collaboration 2021; Heymans et al. 2021; DES Collaboration 2022; DESI Collaboration 2024). However, extending experiments to a broader range of wavelengths boosts the quantity of available data. Furthermore, multiwavelength observations provide robust checks of systematic influence due to the alternative instrumental approach required to survey other wavelengths. There can also be vastly different complex astrophysics for tracers at different wavelengths. Thus alternative tracers can validate the consistency of cosmological conclusions.

Of particular interest for surveying large-scale structure is mapping redshifted neutral hydrogen (HI) gas at 21 cm radio wavelengths. In the post-reionization Universe, at redshifts $z \lesssim 6$, the remaining HI is overwhelmingly contained within dark matter haloes (Villaescusa-Navarro et al. 2018), meaning the 21 cm radiation spontaneously emitted from HI acts as a tracer of the large-scale structure. However, at the redshifts of interest for large-scale structure cosmology ($z \gtrsim 0.1$), the 21 cm emission becomes faint and hard to detect with the high signal-to-noise required for a comprehensive mapping survey. This problem can be circumvented by abandoning the requirement to resolve individual galaxies and instead integrate the combined unresolved emission with low-aperture radio telescopes. This process is now commonly known as HI intensity mapping (Bharadwaj et al. 2001; Battye, Davies & Weller 2004; Chang et al. 2008; Wyithe, Loeb & Geil 2008).

The novel approach of intensity mapping presents some challenges. Since it records all radiation entering the telescope receiver within its frequency bandpass, the contribution to the signal will not come from cosmological HI alone; in fact, only a small percentage is. There are multiple sources of foregrounds in the frequency

<footnote>* E-mail: steven.cunnington@manchester.ac.uk (SC); jywang@shao.ac.cn (JW)</footnote>
ranges corresponding to redshifted 21 cm emission. These mostly originate from our Galaxy and come from diffuse synchrotron and free-free emission, but there can also be strong extragalactic point sources (Mozdzen et al. 2019; Lian et al. 2020; Irfan et al. 2021; Matshawule et al. 2021; Spinelli et al. 2021; Hurley-Walker et al. 2022). These sources dominate the cosmological HI brightness temperature by three to five orders of magnitude. Furthermore, instrumental effects can distort the foreground spectra, making their removal additionally challenging (Shaw et al. 2015; Carucci, Irfan & Bobin 2020; Matshawule et al. 2021; Spinelli et al. 2022). In addition to foregrounds, the telescopes will detect terrestrial, human-made radio frequency interference (RFI) (Harper & Dickinson 2018; Engelbrecht et al. 2024). This can dominate even more than the foregrounds, but vary with time and sky position so can therefore be flagged and deleted from the final maps. Residual RFI will also average down if repeated scans are performed over the same patch of sky, as is the case with MeerKLASS observations.

Despite foreground contamination and RFI challenges, the potential benefits from HI intensity mapping are too large to ignore, thus justifying the ongoing endeavour to develop the technique into a tool for precision cosmology. Experiments have already delivered promising results, demonstrating how detections are possible even without purpose-built instruments (Masui et al. 2013; Wolz et al. 2017; Anderson et al. 2018; Wolz et al. 2022). New experiments are starting to arrive (e.g. Nan et al. 2011; Battye et al. 2012; Chen 2012; Newburgh et al. 2016; PUMA Collaboration 2019; Vanderlinde et al. 2019; CHIME Collaboration 2022; Pal et al. 2022), many designed specifically with HI intensity mapping in mind, creating a bright future ahead.

The MeerKAT telescope,<fnref n="1" /> an array of 64 dishes in the Upper Karoo region of South Africa, has already showcased its capability to perform HI intensity mapping both as an *interferometer*, claiming the first detection of the HI power spectrum (Paul et al. 2023), and also in *single-dish mode*, where the field-of-view can reach cosmological scales. Using the single-dish pilot survey (Wang et al. 2021), a detection of the cross-correlation power spectrum was achieved with overlapping WiggleZ galaxies (Drinkwater et al. 2010; Cunningham et al. 2022).

The MeerKAT Large Area Synoptic Survey (MeerKLASS) (Santos et al. 2017) will continue adding to these observations until it is merged with the SKA Observatory (SKAO) (SKA Cosmology SWG et al. 2020) on the same site.<fnref n="2" /> MeerKAT can observe in two frequency bands relevant for 21 cm-cosmology: L-band (900 MHz < *v* < 1420 MHz, reaching redshifts *z* < 0.58) and UHF-band (580 MHz < *v* < 1000 MHz, 1.45 > *z* > 0.4).

In this work, we present results from MeerKAT single-dish observations in the MeerKLASS L-band deep field, a set of observations performed over a single 236 deg<fnref n="2" /> patch with 41 repeated scans for each of the 64 dishes. Integrating over the sky in this way will average down *time-varying* contributions from thermal noise and RFI, thus boosting the signal-to-noise of *time-constant* sky signals, which includes the foregrounds and the cosmological HI. Previous analysis (e.g. Cunningham et al. 2022; Wolz et al. 2022) exploited the fact that the foreground cleaned intensity maps were thermal noise dominated when estimating error bars on the power spectrum measurements. For this MeerKLASS deep-field and future MeerKLASS observations, it can no longer be assumed that the thermal noise will be dominant therefore demanding more emphasis

on covariance estimation. We therefore explore using the scatter in the foreground transfer function (Cunnington et al. 2023), primarily used to correct for signal loss from foreground cleaning, as a measure of the covariance in measurements. The transfer function is constructed using mock signal injection into the real data, hence its scatter should encapsulate cosmic variance, thermal noise fluctuations, residual RFI contributions, foreground residuals and signal loss correction uncertainties.

We validate the precision and accuracy of the analysis pipeline to ensure estimator accuracy is within the statistical error associated with these MeerKLASS observations. We formalize the main steps in our pipeline, such as beam reconstitution, foreground cleaning, and sampling to a Cartesian grid, presenting validation steps on a suite of simulations which emulate the MeerKLASS L-band deep-field observations. We also use a small sample of overlapping spectroscopic galaxies from the Galaxy And Mass Assembly (GAMA) (Drinkwater et al. 2010) survey to study the cross-power spectrum and its covariance under different analysis choices.

The paper is outlined as follows: Section 2 outlines the observations for the MeerKLASS deep field, our new calibration strategy (with an iterative self-calibration process) and overlapping GAMA galaxies. Section 3 presents the methods for generating our suite of mock data which emulate the observations. Section 4 provides details of our foreground cleaning technique. Section 5 focuses on the power spectra, formalizing our estimation and modelling process, and presents the results. Section 6 shows the results from stacking the MeerKLASS deep-field maps onto the locations of GAMA galaxies and discusses the process and complications that come with this. We then finally conclude in Section 7.

## 2 OBSERVATIONAL DATA

### 2.1 MeerKLASS L-band deep-field

Using the same scan strategy as the MeerKLASS L-band pilot survey in Wang et al. (2021), a series of observations were conducted using the MeerKAT telescope between 2021 September 1 and December 29. Observations with single dishes typically require a scanning strategy where the dishes are rapidly moved across the sky. The MeerKAT antennas were set to scan in azimuth at constant elevation to minimize fluctuations of ground spill and airmass. To make sure the projected speed on the sky is $5 \text{ arcmin s}^{-1}$, the telescope scan speed was set to $5/\cos(\ell) \text{ arcmin s}^{-1}$ (where $\ell$ is the elevation angle) along the azimuth. This speed also ensures that the telescope pointing does not move significantly compared to the width of the primary beam (~1 deg) during a single time dump (2 s).

With a time resolution of 2 s, this gives a scan speed of no more than 10 arcmin per time sample, which is well within the beam size. Noise diodes attached to each receiver were fired for 0.585 s once every 19.5 s during the observation to provide a relative time-ordered data (TOD) calibration reference. The dishes are moved back and forth with a slew of about 10 deg sky projection in each direction, corresponding to an observing time of about 120 s per stripe. The duration of each set of scans with 48 strips is about 100 min, and before and after each scan we spent ~8 min tracking a nearby celestial point source to use as a bandpass calibrator and absolute flux calibrator.

At fixed elevation, two scans can be performed per night, corresponding to when the field is rising and setting respectively. The two scans will cross each other to achieve good sky coverage in the region of overlap. There were 41 observations in total (each one referred as a *block*), the details of each are listed in the Appendix in Table A1. The

<footnote n="1">¹www.sarao.ac.za/science/meerkat</footnote>

<footnote n="2">²www.skao.int/en</footnote>
observations repeated scans covering a patch of sky approximately 236 deg² in the southern sky, spanning in R.A. (330°, 360°) and Dec. (−36°, −25°). The 62 h obtained for each of the 64 dishes before flagging, provide a cumulative 3968 h of observations,<fnref n="3" /> likely to be the deepest single-dish intensity maps, in terms of observational hours per square degree, in the MeerKLASS campaign.

## 2.2 Standard MeerKLASS calibration

We employ the MeerKAT single-dish calibration pipeline KATcali to calibrate the deep-field observations and create the combined 3D data cubes. KATcali includes several cycles of RFI flagging, calibration, and map-making. The details of the calibration are thoroughly described in Wang et al. (2021). Here, we provide a brief overview of the main points. We take a two-step approach to the calibration process. The first step is to use the flux model of a strong point source during tracking observations to calibrate the bandpass and absolute flux, as well as the injection power of the noise diode. The second step then uses the noise diode injections to calibrate the TOD during the constant-elevation scanning phase.

To calibrate the scanning data, we build the model as

$$T_{\text{model}}(t, v) = T_{\text{diffuse}}(t, v) + T_{\text{el}}(t, v) + T_{\text{diode}}(t, v) + T_{\text{rec}}(t, v), \quad (1)$$

where $T_{\text{el}}$ is the antenna temperature model of elevation-dependent terrestrial emission (atmosphere and ground spill), $T_{\text{diode}}$ is the noise diode contribution, $T_{\text{rec}}$ is the receiver temperature, while the diffuse celestial components $T_{\text{diffuse}}$ includes Galactic emission and the cosmic microwave background (CMB),

$$T_{\text{diffuse}}(t, v) = T_{\text{Gal}}(t, v) + T_{\text{CMB}}, \quad (2)$$

where the CMB emission is set to be a uniform $T_{\text{CMB}} = 2.725$ K (Fixsen 2009), and the Galactic emission $T_{\text{Gal}}(t, v)$ is derived from PySM (Thorne et al. 2017) that is based on the reprocessed Haslam 408 MHz total intensity map (Remazeilles et al. 2015).

To compare with the raw data from the correlator, we multiply the signal model by the gain, $g(t, v)$, so that

$$\hat{T}_{\text{model}}(t, v) = g(t, v) T_{\text{model}}(t, v), \quad (3)$$

where the hat indicates a temperature in the (arbitrary) correlator units.

The signal and gain models are then fitted (with a Bayesian framework; see more details in section 3.3 of Wang et al. 2021) to the TOD for each polarization, frequency channel, dish, and observation scan, all of which are treated independently. The function of gain over time is a fourth-order polynomial (equation 17 in Wang et al. 2021). We checked that the order of the polynomial was sufficient to represent the data; choosing higher orders produced noise-dominated coefficients. Using the gain solution, we can then obtain the calibrated temperature,

$$T_{\text{cal}}(t, v) \equiv \frac{\hat{T}_{\text{raw}}(t, v)}{g(t, v)}, \quad (4)$$

where $\hat{T}_{\text{raw}}$ is the raw data collected by MeerKAT. The residual temperature difference between the signal model and calibrated data is then

$$T_{\text{res}}(t, v) \equiv T_{\text{cal}}(t, v) - T_{\text{model}}(t, v). \quad (5)$$

<figure><img src="image_2.png" /><figcaption>Figure 1. Flowchart showing each step in the calibration pipeline. This is an extension beyond the Wang et al. (2021) pipeline since it includes the new self-calibration strategy.</figcaption></figure>

In standard calibration, we assumed that $T_{\text{res}}$ is noise-dominated during the fitting, ignoring the faint point sources in the target field. This was a fair assumption when the noise level was relatively high (Wang et al. 2021). In this work, we seek a better calibration approach by implementing a self-calibration method.

## 2.3 Self-calibration based on the MeerKLASS sky map

The higher signal-to-noise from the MeerKLASS deep-field data may allow us to deal more carefully with the point sources, which were ignored in the standard calibration, to achieve a better calibration result. Additionally, we aim to use a more accurate sky model to replace the $T_{\text{Gal}}$ model that is based on the Haslam map. The 408 MHz radio continuum all-sky map, which has been widely used to study diffuse synchrotron radiation from our Galaxy and works as a Galactic foreground model for cosmological research, combines data from four different surveys that happened more than 40 yr ago (Haslam et al. 1981; Haslam et al. 1982). The instrument performance in the 1980s limited the data quality, especially for the Galactic substructures and the zero level of the maps (Remazeilles et al. 2015). The selection of the spectral index map, used to generate Galactic maps at MeerKAT frequencies, is also a matter of debate (Spinelli et al. 2022).

We therefore extend the calibration process described in Section 2.2 and test the behaviour of the results by introducing an iterative self-calibration process with several loops. This includes (1) replacing $T_{\text{diffuse}}(t, v)$ with $T_{\text{sky}}(t, v)$, obtained from the previous Level-6 $T_{\text{sky}}$ result, when the corresponding Level-6 pixel has a hit-count > 40 (which ensures that $T_{\text{sky}}$ values are obtained from high signal-to-noise data). Here Level-6 $T_{\text{sky}}$ are the combined calibrated temperature maps for all scans and all dishes at all frequencies, after several rounds of RFI flagging, as shown in Fig. 1.

(2) setting the time-order weight $w(t_i) = 1$ for all the scan time samples; instead of setting a higher weight for the diode (main calibrator during the standard calibration) injection samples (see

<footnote n="3">³This number is reduced for the final useable data since some observation blocks were entirely flagged and some dishes failed in observation. We provide details of this later in this section.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The standard deviations of $T_{res}$ for different self-calibration cycles (Loops 1–5) in the 971–1023 MHz region. The standard calibration value is also shown for reference (first point).</figcaption></figure>

more details in section 3.3 of Wang et al. 2021). Thus equation (1) becomes

$$T_{\text{model}}^{\text{upd}}(t, \nu) = T_{\text{sky}}(t, \nu) + T_{\text{el}}(t, \nu) + T_{\text{diode}}(t, \nu) + T_{\text{rec}}(t, \nu). \quad (6)$$

In this way, we have all components (especially point sources and Galactic substructures) included in the input ‘model’. The difference between the data and our calibrated sky data is

$$T_{\text{res}}^{\text{upd}}(t, \nu) \equiv T_{\text{cal}}(t, \nu) - T_{\text{model}}^{\text{upd}}(t, \nu), \quad (7)$$

which should be noise-dominated compared to that in equation (5). The self-calibration process is iterative for five loops, as shown in Fig. 1, in which the $T_{sky}(t, \nu)$ always uses the Level-6 result of the previous calibration loop. Each iteration refines the calibration, reducing the discrepancies between the observed data and the model. We stopped the self-calibration iteration after five loops because the standard deviation of $T_{res}$ converged to a stable value, as shown in Fig. 2.

## 2.4 Characteristics of final MeerKLASS data cubes

### 2.4.1 Final data cubes

The combined calibrated temperature cubes (with a pixel size of 0.3 deg and frequency resolution of 0.2 MHz) are derived from all available scans and dishes. However, predominantly due to RFI, an identified repeated source of which we discuss in the following sub-section, only 27 of the 41 observation blocks yielded well-calibrated results. We mark them in the last column of Table A1. In Fig. 3 we present the number count of good (flagged time < 50 per cent) channels within the 971–1075 MHz band (channels 550–1050) for each block and each dish. In particular, we found that antenna m059 has a higher chance of surviving calibration even when all other dishes fail.<fnref n="4" /> This may be because m059 is the easternmost dish of MeerKAT, located far from the array core, providing it with a more isolated environment.

We show the final sky temperature map from the combination of the remaining 27 scans after 5 loops of self-calibration at 1023 MHz in the middle panel of Fig. 4. For comparison, the top panel shows the

PySM diffuse Galactic emission model with the same sky coverage. We also show in the bottom panel the total number of time samples (across all dishes and observations) that were combined in each pixel. The temperature offset between the top and middle panels is mainly attributed to the CMB, with additional contributions from point sources.

To demonstrate the improvement in data quality from the deep-field intensity maps, Fig. 5 plots the amplitude of each pixel along the line-of-sight as a function of frequency. We compare this to the previous MeerKLASS L-band pilot survey (top-panel) relative to the deep field used for this work. The maps should be dominated by the continuum synchrotron spectra, thus their smoothness is a rough proxy for how dominant any systematics are. There are numerous perturbations in the pilot survey spectra that have consequences for the effectiveness of foreground cleaning, required to access the cosmological HI. However, for the deep field, the spectra appear smoother and well-calibrated. In addition to the improvements shown in Fig. 2, we also find a marginal performance improvement on the line-of-sight rms when using the self-calibration technique relative to the previous standard calibration in the deep-field. The average variance along the line-of-sight of each pixel decreases to 977.2 mK² for the self-calibration, down from 984.3 mK² for the standard calibration.

### 2.4.2 RFI contaminations

From the 41 observation blocks, we had to remove 14 blocks due to excessive, low-level RFI contaminations. In affected blocks, these contaminations were sometimes recognizable on the map level of individual dish data as horizontal stripes. After careful investigations, we have now identified the source of the RFI contamination as the downlink of the mobile communication tower (~54 km away) that enters through the sidelobes and causes a non-linear gain effect in the entire band. Our aggressive flagging approach in this study, with the entire removal of 14 blocks, was to ensure high data quality for our analysis. We are currently developing a method to correct for this effect, which should allow us to use the currently flagged data for future projects.

### 2.4.3 Noise level estimation

Here we estimate the noise levels from the final maps to check if they are consistent with theoretical expectations. We use the final data cube, $T_{sky}$, to estimate the noise level using the difference between four neighbouring channels (the so-called ABBA technique). By adopting the ABBA pattern, this difference statistic removes both constant and linearly varying signals to cancel signals detected by the receiver, leaving only the noise. The noise level is estimated by

$$\Delta T_i(\nu^*) = \frac{1}{2} \left( T_{\text{sky}}^i(\nu) + T_{\text{sky}}^i(\nu + \delta\nu) \right) - \frac{1}{2} \left( T_{\text{sky}}^i(\nu - \delta\nu) + T_{\text{sky}}^i(\nu + 2\delta\nu) \right), \quad (8)$$

where the index i goes over all pixels in the map and * means the combined result from the four channels [$\nu - \delta\nu$, $\nu$, $\nu + \delta\nu$, $\nu + 2\delta\nu$]. From the observation data, we take a weighted rms,

$$\Delta T_{\text{RMS}}^2(\nu^*) = \frac{N_{\text{pix}}}{N_{\text{pix}} - 1} \left( \frac{\sum_i w_i \Delta T_i^2}{\sum_i w_i} - \left[ \frac{\sum_i w_i \Delta T_i}{\sum_i w_i} \right]^2 \right), \quad (9)$$

where the sum is over the number of pixels in the map, $N_{pix}$, and we have suppressed the dependence of $w_i$ and $\Delta T_i$ on $\nu^*$. For the

<footnote n="4">⁴For consistency with other MeerKAT work, we indicate MeerKAT’s 64 antennae with this notation spanning m000-m063.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Number count of good (flagged time < 50 per cent) channels within the 971–1075 MHz band (L-band channels 550–1050) for each dish (x-axis) and observation block (y-axis). We only show the 27 observation blocks used, out of the possible 41. The other blocks were entirely removed mainly due to RFI contamination, and one failed observation run. The details of all 41 blocks are listed in Table A1.</figcaption></figure>

weight, $w_i$, we used the theoretical expected variance itself, $w_i(v^*) = 1/\sigma_{th,i}^2(v^*)$. The corresponding variance in each pixel would be,

$$ \sigma_{th,i}^2(v^*) = \frac{1}{4} \left( \sigma_{th,i}^2(v) + \sigma_{th,i}^2(v+\delta v) + \sigma_{th,i}^2(v-\delta v) + \sigma_{th,i}^2(v+2\delta v) \right). \quad (10) $$

Here the noise in the initial polarized TOD is given by the radiometer equation

$$ \sigma_{\text{th,pol}}^2(t, v) = \frac{T_{\text{sys}}^2(t, v)}{\delta v \delta t}, \quad (11) $$

where $\delta v = 0.209$ MHz is the frequency width of the MeerKAT L-band channels and $\delta t = 2$ s is the length of each time stamp recording. For the system temperature we use the calibrated data itself, $T_{\text{sys}} = T_{\text{cal,HH}}(t, v)$ or $T_{\text{cal,VV}}(t, v)$ for horizontal and vertical polarization data, respectively. Thus the noise in the initial intensity TOD is

$$ \sigma_{\text{th}}^2(t, v) = \frac{T_{\text{cal,HH}}^2(t, v) + T_{\text{cal,VV}}^2(t, v)}{4 \delta v \delta t}. \quad (12) $$

In order to get the final data cube, $T_{\text{sky}}^i(v)$, the data goes through several stages of averaging. We propagate the variance taking into account this averaging in order to get to the final $\sigma_{th,i}^2(v)$.

The theoretical rms which can be compared to $\Delta T_{\text{RMS}}(v^*)$ in equation (9), is then simply

$$ \sigma_{\text{th}}^2(v^*) = N_{\text{pix}} \left( \sum_i \sigma_{th,i}^{-2}(v^*) \right)^{-1}. \quad (13) $$

More details can be found in section 5.6 of Wang et al. (2021).

In the upper panel of Fig. 6, we show the noise levels (from data and theory) in the 971–1075 MHz frequency band (channels 550–1050). Note here both the real and theory values are related to the real $T_{\text{sys}}$ and hit counts in corresponding map pixels. The count reduction due to RFI flagging, and unflagged low-level RFIs that are hiding in the data would cause wiggles. Four groups of spikes are present at higher frequencies, while below 1023 MHz, there is much less scatter compared to the previous MeerKLASS

pilot survey result. We therefore chose to restrain our cosmological analysis to the 971–1023 MHz region, for which we show the noise level in the middle panel, finding a median of $\Delta T_{\text{RMS}} = 1.21$ mK. In the lower panel of Fig. 6, we also show the histogram of the ratio $R_{\text{RMS}} = \Delta T_{\text{RMS}}(v^*)/\sigma_{\text{th}}(v^*)$ in the 971–1023 MHz band. We find a median value of $R_{\text{RMS}} = 1.2$, while some channels are noisy and $R_{\text{RMS}}$ can reach 1.5.

## 2.5 GAMA spectroscopic galaxies

Overlapping with the MeerKLASS deep field is a catalogue of spectroscopic galaxies recorded in the Galaxy and Mass Assembly (GAMA) survey. GAMA began observing in February 2008 (Driver et al. 2011; Liske et al. 2015). The survey aimed to construct catalogues of galaxies to test the cold dark matter paradigm and probe galaxy evolution (Driver et al. 2009). For our purposes, we exploit the fact that these galaxies will trace dark matter and should therefore correlate with the overlapping MeerKLASS intensity maps after foreground cleaning.

We use the fourth and final data release from the Galaxy and Mass Assembly survey, GAMA DR4 (Driver et al. 2022). This includes galaxies from the GAMA II survey phase which enlarged the GAMA I equatorial survey, by including extra survey regions. One of these extra regions was the 23 h (G23) field, which completely overlaps in area with the MeerKLASS L-band deep field, shown by the blue scatter points in the bottom panel of Fig. 4, which represent the GAMA G23 galaxies relative to the observational time footprint map for the L-band deep field. The G23 field covered 339 deg $< R.A. < 351$ deg and $-35$ deg $< \text{Dec.} < -30$ deg. The magnitude limit of the field was reduced from the planned 19.8 in r-band to 19.2 in i-band due to limited observational time, but this allowed the GAMA II survey to obtain a high 94.2 per cent redshift completion in the G23 region (Liske et al. 2015).

The GAMA DR2 only provides random catalogues that emulate the survey selection function of the real galaxies for the G09, G12, and G15 regions (not the G23 used in this work). Ideally, these would be used to construct a survey completion window function, which would be used in the cross-power spectrum estimation. In the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** MeerKLASS L-band deep-field map at 1023 MHz (middle panel) in comparison with the Galactic emission (top panel) derived from PySM. Blue circles mark the positions of known point sources with flux > 1 Jy at 1.4 GHz. The bottom panel shows the number of observational recordings, i.e. * timestamps* in each pixel at 1023 MHz, from the surviving 27 blocks. Each timestamp is 2 s and we use this for noise estimation and weighting in the power spectrum analysis. Blue scatter points are the positions of all overlapping GAMA galaxies.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Amplitude of each map pixel along the line of sight for the previous MeerKLASS pilot survey (top panel) and the deep field used in this work (bottom panel). We only show the ~250 channels most clean of RFI used in later cosmological analysis. This has the frequency range 971.2 < ν < 1023.6 MHz corresponding to 0.46 < z < 0.39 for redshifted 21 cm emission.</figcaption></figure>

absence of this, we assume a uniform survey completion across the G23 area and a constant distribution with redshift, given the very flat distribution for the redshift range used in this analysis (green region of Fig. 7). Given our signal-to-noise levels, survey masks will be a subdominant effect in the cross-power, thus this is an adequate approximation to make.

Due to our aggressive flagging strategy on the MeerKLASS L-band data, the overlap region in redshift is restricted to 0.39 < z < 0.46. As shown by Fig. 7, this catches the GAMA survey galaxy redshift distribution, $N_g(z)$, near the end of its optimal range, but still provides 2269 galaxies to cross-correlate with. At this redshift range, the comoving (volumetric) number density of the GAMA galaxies is $\bar{n}_g = 4.8 \times 10^{-4} h^3 \text{ Mpc}^{-3}$.

## 3 SIMULATION OF MOCK DATA

To accompany the calibrated observations outlined in Section 2, we produce a suite of 500 mock simulations, used for validating the analysis pipeline, calculating signal loss from foreground cleaning and estimating covariance in the power spectrum measurements. To create this large number of mocks we generate lognormal density fields (Coles & Jones 1991) from which a map of HI brightness temperatures can be cut, and galaxies Poisson sampled onto for the spectroscopic galaxy counterpart. The lognormal distribution has been shown to be a more accurate representation of the PDF for the distribution of cosmic matter, relative to a more simple Gaussian approximation (Agrawal et al. 2017; Shin et al. 2017). The consequences of only using a Gaussian approximation have not been fully investigated but could potentially impact the accuracy of covariance and signal loss estimation. We therefore adopt the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between noise levels derived from calibrated observational data [$\Delta T_{RMS}(v^*)$] and from theory equation [$\sigma_{th}(v^*)$] in 971–1075 MHz (upper panel), a zoom-in for above in 971–1023 MHz (middle panel), and the histogram of the ratio $R_{RMS} = \Delta T_{RMS}(v^*)/\sigma_{th}(v^*)$ in 971–1023 MHz (lower panel).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Redshift distribution of GAMA G23 galaxies in the 23 h field. The accessible redshifted 21 cm emission for the full L-band range is shown by the light-pink shading. The region used in this work between 0.39 < z < 0.46 (1023.6 < v < 971.2 MHz) is shown by the darker green shading, in which 2269 GAMA galaxies exist.</figcaption></figure>

lognormal approach since this adds no extra computational demand yet will provide a more realistic suite of mocks.

An important stage of the analysis pipeline is regridding the intensity mapping voxels with celestial sky coordinates (R.A., Dec., v) onto a regular 3D Cartesian grid with lengths $l_x, l_y, l_z$ in h⁻¹Mpc, on which a Fast Fourier Transform (FFT) can be applied for the power spectrum estimation (Cunnington & Wolz 2024). We therefore

require the mocks to begin in celestial coordinates so we can validate the regridding process. To achieve this, we begin by fitting a Cartesian cube around the full MeerKLASS deep-field footprint, creating a grid with dimensions $l_x, l_y, l_z = 774.9, 335.7, 252.2$ h⁻¹Mpc over $n_{cell} = 512^3$ cells. We use a large number of cells for this initial grid to approximately emulate a continuous sky, aiming to avoid the simulated maps resolving the finite input grid structure. Onto this grid, we simulate a lognormal density field. This is done by interpolating k-modes from a model power spectrum onto the FFT grid. The HI model power spectrum is given by

$$P_{\mathrm{HI}}(\mathbf{k}) \equiv P_{\mathrm{HI}}(\mathbf{k}, \mu) = \frac{\bar{T}_{\mathrm{HI}}^2 b_{\mathrm{HI}}^2 (1 + f \mu^2 / b_{\mathrm{HI}})^2 P_{\mathrm{m}}(\mathbf{k})}{1 + (\mathbf{k} \mu \sigma_v / H_0)^2} \quad (14)$$

The anisotropy in this model makes it convenient to parametrize k using $\mu$, defined as the cosine of the angle between the line of sight and k. The HI parameters are the mean HI gas temperature ($\bar{T}_{\mathrm{HI}}$) and the HI bias ($b_{\mathrm{HI}}$). $f$ is the growth rate of structure and is included to model linear redshift space distortions (RSD) which enhance clustering across all scales. The denominator is included to approximately model the non-linear RSD; the so-called Fingers-of-God (FoG), which suppresses clustering at small scales. The FoG are modulated by $\sigma_v$, which is the velocity dispersion of the tracer, in this case, the HI gas. All of these quantities have a redshift dependence, which we omit for brevity; they will also have little variation given our narrow redshift range. Lastly, $H_0$ is the Hubble parameter at z=0, known as the Hubble constant. An alternative process for generating brightness temperatures is pursued for our stacking analysis discussed later in the paper. More realistic lognormal line intensity mapping mocks are available. For example, in future work, we intend to adopt the approach in Niemeyer, Bernal & Komatsu (2023) where luminosities from the spectral lines are painted given a luminosity function. This has been shown to reproduce better RSD for galaxy distributions and models the whole temperature distribution expected in the maps.

To obtain the HI correlation function, $\xi_{\mathrm{HI}}$, needed for the lognormal transform we first inverse Fourier transform $P_{\mathrm{HI}}/\bar{T}_{\mathrm{HI}}^2$, where $P_{\mathrm{HI}}$ is the model in equation (14) and we factor out the background average temperature $\bar{T}_{\mathrm{HI}}$ to remove the temperature units. The correlation function is then log-transformed via $\xi'_{\mathrm{HI}} = \ln(1 + \xi_{\mathrm{HI}})$. Fourier transforming $\xi'_{\mathrm{HI}}$ back gives the modified lognormal power spectrum $P'_{\mathrm{HI}}$, without temperature units. Gaussian random numbers are then sampled into each cell of the Fourier grid $\delta_G(\mathbf{k})$ with variance $\sigma^2(\mathbf{k}) \propto P'_{\mathrm{HI}}(\mathbf{k})$. This Gaussian random field is then inverse Fourier transformed to real space giving $\delta_G(\mathbf{x})$ from which the lognormal HI density field is given by

$$\delta_{\mathrm{HI}}(\mathbf{x}) = \exp\left[\delta_{\mathrm{G}}(\mathbf{x}) - \frac{\sigma_{\mathrm{G}}^2}{2}\right] - 1, \quad (15)$$

where $\sigma_G^2$ is the measured variance of the Gaussian random number field $\delta_G$. This is converted into the HI brightness temperature field given by $\delta T_{\mathrm{HI}} = \bar{T}_{\mathrm{HI}}\delta_{\mathrm{HI}}$. Further details on this process can be found in the appendix A of Beutler et al. (2011). For the galaxy overdensity field $\delta_g(\mathbf{x})$, we follow the same procedure except the input model is for the galaxy power spectrum, given by

$$P_g(\mathbf{k}) \equiv P_g(k, \mu) = \frac{b_g^2 (1 + f \mu^2 / b_g)^2 P_m(k)}{1 + (\mathbf{k} \mu \sigma_v / H_0)^2}, \quad (16)$$

where $b_g$ is the bias for the galaxies. We assume that the velocity dispersion $\sigma_v$ for the HI gas and galaxies is the same for simplicity. This will be sufficient for this work where signal to noise is low at the small scales where $\sigma_v$ has the main effect.
The galaxy overdensity $\delta_g$ is unit-less whereas the intensity maps $\delta T_{HI}$ are brightness temperature fluctuations measured in mK. We use the same seed in the random number generation to ensure that both the $\delta T_{HI}$ and $\delta_g$ correlate. The lognormal code we adopt utilizes some functions from the publicly available package powerbox (Murray 2018).<fnref n="5" />

## 3.1 Transforming between Cartesian grids and sky maps

As already mentioned, the galaxies and the intensity mapping voxels from observations are in celestial coordinates at a certain frequency (or redshift) i.e. $(\theta, \nu) \equiv (\theta_{R.A.}, \theta_{Dec.}, \nu)$. However, for the final power spectrum estimation, we require the fields be in a Cartesian space with dimensions $(l_x, l_y, l_z)$ in $h^{-1}$Mpc. For this we follow the recent work in Cunnington & Wolz (2024) which we discuss further in Section 5.

Creating lognormal mocks in celestial coordinates that emulate observations requires the inverse of this step i.e. $\delta T_{HI}(x) \rightarrow \delta T_{HI}(\theta, \nu)$, since the lognormal fields are generated onto Cartesian grids measured in $h^{-1}$Mpc. Transforming mock galaxy coordinates (the Poisson sampling of which is discussed in Section 3.3) from $(l_x, l_y, l_z)[h^{-1}Mpc] \rightarrow (\theta_{R.A.}, \theta_{Dec.}, \nu)$ is trivial. We can use as- tropy<fnref n="6" /> (Astropy Collaboration 2022) transformation routines on each galaxy's coordinate. However, transforming the tessellating field of Cartesian cells, each with an associated brightness temperature, into a gap-less map of voxels with identical celestial coordinates and frequency binning as the real MeerKAT maps, requires a resampling of the field. For this we also utilize the techniques and public code<fnref n="7" /> outlined Cunnington & Wolz (2024).

This process involves a Monte-Carlo particle sampling technique where each Cartesian cell is filled with $N_p$ particles, randomly distributed within the cell. The particles are assigned the brightness temperature of the cell they fall within. The particles are then transformed to celestial coordinates $(\theta_{R.A.}, \theta_{Dec.}, z)$ and added into the map voxels associated with their celestial coordinate and the frequency channel corresponding to their redshift. Each voxel is normalized by the number of particles falling within it. Providing enough sampling particles per cell are used (in this work we use $N_p = 10$) and the Cartesian grid is of sufficient resolution (hence our choice of $n_{cell} = 512^3$), it will provide an intensity map in celestial sky coordinates, that traces the fluctuations generated by the lognormal field.

Throughout this work, we do not require the necessary extensions of higher order particle assignment or field interlacing to suppress the effects of aliasing (as trialled in Cunnington & Wolz 2024), since these effects are below the current statistical noise from this data. We are however, still required to implement a correction to the power spectrum estimation from the nearest-grid-point assignment which smooths the field, which we formalize in Section 5.

## 3.2 HI intensity map mocks

Using the regridding process described in the previous section we are provided with HI intensity mapping mocks $\delta T_{HI}(\theta, \nu)$, with identical coordinates and dimensions as the observed MeerKAT intensity maps. From here we can add in observational effects and simulate various systematics. The most dominant effect comes from

foreground contamination. Systematics will have a large influence on foreground contamination, and this can be extremely complex. Therefore, this is not something we attempt to directly simulate in this work since we risk underestimating the problem. Instead, we use the observational data itself which has the relevant systematics contained within, albeit only one realization of them. However, adding these directly onto different realizations of HI mocks, still recreates stochasticity in the HI signal's response to systematics. We discuss this further in Section 4.

Another important addition to the mocks will be a smoothing to the intensity maps corresponding to the angular pixelization. This is necessary to emulate the radio telescope's beam, which for single-dish intensity mapping is quite broad giving a low angular resolution to the maps and damping small-scale fluctuations perpendicular to the line of sight. In this work, we assume the MeerKAT beam profile to be Gaussian, which is a simplification but still valid given the signal to noise at small scales where the beam suppresses signal below the noise. We will investigate implementing a more realistic MeerKAT beam profile in future work (see e.g. Asad et al. 2021; Matshawule et al. 2021; Spinelli et al. 2022). The assumed Gaussian beam profile has a frequency ($\nu$) dependent standard deviation in radians given by

$$ \sigma_{\text{beam}}(\nu) = \frac{1}{2\sqrt{2\ln 2}} \frac{c}{\nu D_{\text{dish}}} \quad (17) $$

where we use $D_{\text{dish}} = 13.5$ m for the diameter of the MeerKAT dishes and $c$ is the speed of light. To emulate the beam effects we smooth the mock maps by creating a 2D Gaussian array kernel centred on the median angular pixel $\theta_0$ given by

$$ \mathcal{B}_{\text{beam}}(\theta, v_i) = \exp \left[ -\frac{1}{2} \left( \frac{\theta - \theta_0}{\sigma_{\text{beam}}(v_i)} \right)^2 \right] \quad (18) $$

This is normalized such that $\sum_{\theta} \mathcal{B}_{\text{beam}} = 1$, so that convolution with this beam kernel does not change the total brightness temperatures within the maps. The HI intensity map mocks are convolved with equation (18) at each frequency $v_i$ to give the beam-smoothed HI map;

$$ \delta T_{\oplus \text{HI}}(\theta, v_i) = \delta T_{\text{HI}}(\theta, v_i) \otimes \mathcal{B}_{\text{beam}}(\theta, v_i), \quad (19) $$

where the $\otimes$ symbol represents a convolution over the angular $\theta$ space.

To realistically emulate additional components such as foregrounds, noise, and their response to systematics, we add the real observations to each mock iteration. This is to avoid the requirement of simulating the complex and poorly understood systematics.

In some cases, mostly for studying the HI autocorrelation, it will be useful to simulate just the thermal noise originating from the thermal signature of the receivers. We do this by generating Gaussian random noise the power of which can then be measured with the same pipeline as the rest of the analysis, providing a robust estimate of the thermal noise floor. We begin with an analytical description of the thermal noise fluctuations per pixel using the radiometer equation (as in Section 2 but now at map-level)

$$ \sigma_{\text{N}}(\theta, \nu) = \frac{T_{\text{sys}}(\nu)}{\sqrt{2 \delta \nu t_{\text{p}}(\theta, \nu)}}, \quad (20) $$

where $t_p$ is the summed observation time contained in each voxel, which can be taken directly from the time stamp counts as shown in Fig. 4, with each count representing $\delta t = 2$ s. This emulates the real data by creating a non-uniform noise map, generating higher noise around the edges where scan coverage is lower. For the system

<footnote n="5">⁵github.com/steven-murray/powerbox</footnote>

<footnote n="6">⁶docs.astropy.org</footnote>

<footnote n="7">⁷github.com/stevecunnington/gridimp</footnote>
temperature, we follow Cunnington (2022) and use a model tuned to coincide with the estimated system temperature of the MeerKAT telescope (Wang et al. 2021),<fnref n="8" /> assuming

$$T_{\text{sys}}(\nu) = T_{\text{rx}}(\nu) + T_{\text{spl}} + T_{\text{CMB}} + T_{\text{gal}}(\nu), \quad (21)$$

where we assume the contribution from our own galaxy is $T_{\text{gal}} = 25\,\text{K} (408\,\text{MHz} \nu^{-1})^{2.75}$, the CMB temperature is $T_{\text{CMB}} = 2.725\,\text{K}$, as in Section 2. The spill-over temperature is assumed as $T_{\text{spl}} = 3\,\text{K}$ and the receiver temperature we model as

$$T_{\text{rx}}(\nu) = 7.5\,\text{K} + 10\,\text{K} \left( \frac{\nu}{\text{GHz}} - 0.75 \right)^2. \quad (22)$$

Due to the narrow frequency band we use, this gives a near-constant system temperature with the range $15.97 < T_{\text{sys}} < 16.02\,\text{K}$. Gaussian random numbers can then be generated across the MeerKLASS deep-field map with the voxel-dependent rms of $\sigma_N(\theta, \nu)$.

### 3.3 Galaxy mock catalogues

To accompany the mock intensity maps, we also require corresponding galaxy mocks with the same correlated clustering. Rather than transforming the Cartesian galaxy overdensity field $\delta_g(\mathbf{x})$ into sky coordinates, it is more realistic to instead Poisson sample the field to generate a catalogue of galaxies with the same number density as the GAMA G23 survey. This also allows additional validity checks such as testing the stacking pipeline which cannot be done without a catalogue of galaxy coordinates.

We generate the galaxy catalogue coordinates initially in the Cartesian domain to exploit the high $n_{\text{cell}} = 512^3$ resolution. We first calculate the average galaxy counts in cells $\langle n_g \rangle$ for the real GAMA G23 region in the redshift range where there is overlap with the intensity map frequency range. We then take the galaxy survey selection window function, which for the GAMA G23 region, we have to assume is constant between $339 < \text{R.A.} < 351$ deg and $-35 < \text{Dec.} < -30$ deg (see discussion in Section 2.5). We then transform this window function into the high resolution Cartesian space. A galaxy number counts field $n_g(\mathbf{x})$ is then given by

$$n_g(\mathbf{x}) = \langle n_g \rangle (\delta_g(\mathbf{x}) + 1) W_g(\mathbf{x}), \quad (23)$$

where we use the transformed window function $W_g(\mathbf{x})$ and the real GAMA average counts $\langle n_g \rangle$ to ensure the mock galaxy counts will share the same survey selection as the real galaxy survey, with the same number density. In this process we normalize the window function $W_g$ such that its maximum value is 1, so for the GAMA G23 survey, this becomes a binary mask of 0 and 1. We then Poisson sample the $n_g$ field to introduce some stochasticity which we expect to be present in the real galaxy catalogues. This provides an integer value in each Cartesian cell for the galaxy counts. The last step is to assign random Cartesian coordinates within the cell boundaries for the $n_g(\mathbf{x})$ galaxies within each cell. At this stage we rely on the Cartesian field having sufficiently high resolution such that the random assignment of a coordinate within the cell does not shift the galaxies away from the peaks of the density fluctuations. The catalogue of galaxy coordinates (x, y, z [$h^{-1}\text{Mpc}$]) can then be transformed into celestial coordinates ($\theta_{\text{R.A.}}$, $\theta_{\text{Dec.}}$, z).

In this work, we generate 500 mock intensity maps with corresponding mock GAMA galaxy catalogues. As some early validation of these mocks, we show in Fig. 8 how the power spectra for both the HI intensity map and galaxy autocorrelations, along with their

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Averaged power spectra (blue solid line) of the 500 regridded mocks and agreement with the fiducial model (black dashed line). Mocks include the beam but no other systematics. Green dotted lines in the bottom panels mark the edges of the 10 per cent residual region, while the shaded regions mark the diagonal variance inferred from the simulations.</figcaption></figure>

cross-correlation, agree with the fiducial model (black dashed line). This demonstrates that our estimator is unbiased and processes such as the regridding are performing to within an acceptable tolerance. The power spectrum estimation process used here on the mocks is identical to that used on the real data and is detailed later in Section 5, along with the modelling.

The mocks used in the test of Fig. 8 do not include any additional systematics beyond the telescope beam in the intensity maps. We will investigate the inclusion of foregrounds and systematics in Section 4. The large scatter shown by the shaded blue region, which represents the $1\sigma$ range across the 500 mocks, is therefore driven by cosmic variance due to the small volume covered by the data. We see that the errors for the galaxy power spectrum appear larger than the HI intensity maps. This is because we have not introduced any

<footnote n="8">⁸See skafrica.atlassian.net/wiki/x/AYCHE</footnote>
thermal noise (or other observational or systematic effects) into the intensity maps. However, the galaxy catalogues will have shot noise contributing to the uncertainties. This nicely showcases the potential for Hi intensity mapping if noise and systematics can be controlled, however, this is not the purpose of this exercise so is not a point we will labour. There are also additional caveats when considering high-density tracers such as Hi. As shown by Foreman, Obuljen & Simonović (2024), despite the low shot noise, stochastic noise and non-linear bias in the Hi field can limit the potential of constraints inferred from the power spectrum.

The green dashed lines in the residual panels mark the regions where there is > 90 per cent accuracy relative to the model which we can see is achieved at all scales for all three cases. This is also well within the statistical uncertainty (blue-shaded regions). The largest residual comes in the cross-correlation where discrepancies seem to increase with k nearly reaching the −10 per cent boundary. We speculate that this may be caused by the Poisson sampling of mock galaxies followed by random assignment of coordinates within the cell. This extra stochasticity is not repeated in the Hi intensity map mocks and could therefore be generating a pseudo-cross-correlation coefficient, which would depart from r = 1 at smaller scales (higher-k). However, this appears to be a minor effect. Furthermore, this could be seen as a feature, since we expect r < 1 in the real observations, which will also have a scale dependence. For these reasons, we deem this to be a suitable result for the purposes of this work, and defer a higher precision investigation and validation to future work.

## 4 FOREGROUND CLEANING

Foreground cleaning Hi intensity maps is a well-studied problem (Santos, Cooray & Knox 2005; Ansari et al. 2012; Wolz et al. 2014; Alonso et al. 2015; Shaw et al. 2015; Olivari, Remazeilles & Dickinson 2016; Carucci et al. 2020; Cunnington et al. 2021; Irfan & Bull 2021; Soares et al. 2022; Spinelli et al. 2022), and generally blind cleaning techniques are the currently preferred method for single-dish experiments. As their name suggests, blind methods do not rely on prior knowledge of the foreground structure. This is advantageous for current experiments because the foregrounds at these frequencies lack observational constraints. Furthermore, blind techniques avoid placing overwhelmingly stringent calibration requirements on the data, since the dominant foregrounds require beyond sub-percent precision modelling for a successful subtraction. Blind techniques instead work by exploiting the foreground dominance and their frequency-correlated structure to assume that they can be reduced to a distinct set of correlated modes separated from the oscillating, subdominant fluctuations which will contain most of the Hi information.

We outline the formalism for our blind foreground cleaning approach in Section 4.3, but first we discuss two pre-processing steps before foreground cleaning. These are designed to improve the foreground clean by reducing the impact from residual RFI or systematics, plus any known frequency-dependent observational effects.

### 4.1 Reconvolution

All radio intensity mapping experiments will have a chromatic beam whose size evolves with frequency. This can add to the number of degrees of freedom in the radial foreground structure, inhibiting foreground removal since more modes need to be removed to perform a suitable clean. Smoothing the maps differently in each frequency channel can homogenize the angular resolution, giving them a

consistent effective beam size at all frequencies. This reconvolution should also suppress residual contamination manifesting as localized amplitude spikes in the temperature maps, smoothing them out so there are fewer discontinuous fluctuations which will impact the PCA and Fourier transforms for power spectrum estimation.

The reconvolution step involves convolving the intensity maps, one frequency channel $v_i$ at a time with the kernel

$$K(\theta, v_i) = \exp \left[ -\frac{1}{2} \frac{(\theta - \theta_0)^2}{\gamma \sigma_{\text{max}}^2 - \sigma_{\text{beam}}^2(v_i)} \right], \quad (24)$$

where $\sigma_{\text{beam}}$ is the beam size as previously defined in equation (17), $\sigma_{\text{max}}$ is the maximum $\sigma_{\text{beam}}(v_{\text{min}})$ value at the lowest frequency, and $\gamma$ is a scaling factor which governs how much the final effective resolution decreases. The reconvolution kernel is normalized such that the sum over all pixels is unity i.e. $\sum_{\theta} K(\theta, v_i) = 1$.

As in previous work (Cunnington et al. 2022), we use a weighted reconvolution to mitigate the impact of edge pixels with low scan coverage. From hereon, the weights for the intensity maps derive directly from the time stamp count in each pixel, as shown in Fig. 4. The weighted reconvolved intensity maps are given by

$$\delta T_{\oplus}(\theta, v_i) = \frac{[\delta T(\theta, v_i)w(\theta, v_i)] \otimes K(\theta, v_i)}{w(\theta, v_i) \otimes K(\theta, v_i)}. \quad (25)$$

As in equation (19), the ⊗ symbol represents a convolution. The weights themselves must also be reconvolved to keep them consistent and optimal for the new $\delta T_{\oplus}$ maps. The reconvolved weights are given by

$$w_{\oplus}(\theta, v_i) = \frac{[w(\theta, v_i) \otimes K(\theta, v_i)]^2}{w(\theta, v_i) \otimes K^2(\theta, v_i)}. \quad (26)$$

In this work we use a value of $\gamma = 1.4$ for the reconvolution kernel in equation (24). We find that neglecting this resmoothing step entirely increases the influence of systematics in the results which is likely due to residual RFI or contaminants not being smoothed and depressing their dominance. We discuss this further in Appendix B1 but also plan a dedicated study into this aspect of the analysis where a more robust model of the MeerKAT beam will also be implemented.

The resmoothed, uncleaned map is shown in Fig. 9 (top panel). Careful comparison with Fig. 4 reveals how some of the foreground structure has been smoothed, albeit subtly, which is the aim. Too much reconvolution will suppress power in high $k_{\perp}$ modes in the final power spectrum estimate. The reconvolution is more noticeable in the time stamp count field, which we use as the weights $w_{\oplus}$, shown by the bottom panel of Fig. 9.

### 4.2 Map trimming

Shown by the time stamp count field of Fig. 4 (and the resmoothed version in Fig. 9), the MeerKLASS scanning strategy of combining numerous rising and setting scans does not generate a uniform coverage. The edges of the patch result in relatively low recorded time stamps due to a lack of overlapping scans in these regions.

We found that implementing an edge trim to the field improved the quality of the foreground clean and final power spectrum estimates. This is motivated and explored in greater detail in Carucci et al. (2024). The improvement from map trimming is because pixels are removed where residual time-varying contamination has been averaged down less, due to fewer scans.

In Fig. 9, the black contour shows the cut to the patch we make, removing everything outside the trimmed lines. This amounts to all pixels at 334° < R.A. < 357° being set to zero in both the weights and intensity map field.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Reconvolved MeerKLASS deep field (top-panel) using $\gamma = 1.4$ (see equation 25). The bottom panel shows the reconvolved time-stamp field which we use for weighting, thus also requires reconvolution (using equation 26) for consistency with the intensity maps. Both panels are at 1023 MHz to match Fig. 4.</figcaption></figure>

## 4.3 Blind foreground cleaning with PCA

We represent the intensity maps with matrices $X$ which are equivalent to the $T(\theta, \nu)$ temperature fields, discretized into arrays on which matrix algebra can be performed.

Blind foreground cleaning assumes the matrix of the observational data cube $X_{\text{obs}}$, which contains the intensity maps with $n_{\theta}$ pixels at each of the $n_{\nu}$ frequency channels, can be generalized by the linear equation

$$X_{\text{obs}} = AS + R. \qquad (27)$$

A is referred to as the *mixing matrix* and is the set of $N_{\text{fg}}$ spectral functions estimated to project out the foregrounds. S are the $N_{\text{fg}}$ source maps obtained by projecting A along the data;

$$S = A^T X_{\text{obs}}. \qquad (28)$$

The dot product AS thus has shape $[n_{\nu}, n_{\theta}]$, along with the other terms in equation (27), and represents a subset of the observations estimated to contain the bulk of the foreground contamination. The residual term R will contain the remaining HI signal and thermal noise. In reality, the situation is more complex and there is cross-term leakage. Both HI signal will leak into the subtracted AS term causing *signal loss*, and foreground will leak into the remaining R

causing *foreground residuals*. The control of both these is crucial and we discuss them further throughout the remaining sections.

Different blind foreground techniques only vary in their approach to estimating the mixing matrix $A$. In this work, we adopt the straightforward principal component analysis (PCA) approach, since it is well-tested and understood. Implementing more sophisticated blind techniques is planned in future work. PCA foreground cleaning begins by computing the frequency covariance matrix of the weighted observations

$$C = \frac{(w \circ X_{\text{obs}}) (w \circ X_{\text{obs}})^T}{n_{\theta} - 1}, \qquad (29)$$

where we have used the intensity mapping weights to optimize the estimation and $w \circ X_{\text{obs}}$ represents the element-wise product between the weights and the observational data matrix. Since the weights are the stacked average of the counts along the line of sight, they are constant in frequency by design to avoid increasing the rank of $X_{\text{obs}}$ i.e. not inject additional structure into the frequency spectra. Weighting in this way means $C$ is not a formally normalized covariance estimation, however, this is not a concern for us since we are only extracting unitless eigenvectors from the covariance and using these to project out modes in the real observed data so sensible amplitudes are recovered in the component separated pieces. We discuss this further, with some examples in Appendix B2. $C$ is then eigen-decomposed by

$$C = U \Lambda U^T, \qquad (30)$$

where $U$ is the matrix of eigenvectors with the corresponding eigenvalues contained in $\Lambda$. The mixing matrix is then formed by taking the first $N_{\text{fg}}$ column vectors $u_i$ from $U$ i.e. $A = [u_1, u_2, ..., u_{N_{\text{fg}}}]$. The first nine eigenmodes are shown in the Appendix plot Fig. B2. It is at this stage a choice is required for $N_{\text{fg}}$, which determines how many modes are removed and deemed as foregrounds. If the $N_{\text{fg}}$ is too high, the HI signal loss increases; on the contrary, if $N_{\text{fg}}$ is too low, residual contamination will dominate more. The estimate of the foregrounds is provided by $X_{\text{fg}} \equiv AS$ which is subtracted to provide the residuals, revealing the final cleaned data used in the analysis i.e. $X_{\text{clean}} \equiv R$ where

$$X_{\text{clean}} = X_{\text{obs}} - X_{\text{fg}}. \qquad (31)$$

Fig. 10 shows the resulting cleaned maps by removing $N_{\text{fg}} = 10$ PCA modes. We show the maps split into six frequency sub-bands (given in the bottom-left of each panel) to illustrate the changing intensity structures. Each panel is the frequency average of the 42 channels in each sub-band showing a reasonably consistent amplitude of temperature fluctuations across the sub-bands. The amplitudes in the maps are still likely driven by non-cosmological contributions, mostly thermal noise, but they are not far from the 0.1 mK fluctuations we expect from HI. As we will show statistically in Section 5, over large-scale fluctuations, the predominant contribution comes from the cosmological HI.

## 4.4 Foreground cleaning tests using mocks

In order to validate our foreground removal and estimate the signal loss it induces, we run the same PCA cleaning technique on mock intensity maps. As discussed in Section 3.2, simulating foreground contamination that emulates the systematics and complex observational effects is beyond current capabilities. We therefore use the real observations themselves and inject mock HI intensity maps into these. We then apply the foreground removal pipeline on the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Foreground cleaned intensity maps for the MeerKLASS L-band deep field in frequency sub-bands (stated in bottom left). Maps have been cleaned by removing $N_{fg} = 10$ PCA modes.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Power spectra for the HI intensity map mocks under different scenarios, each averaged over 500 mock realizations. Including: HI-only mock (grey dotted line). HI plus thermal noise realizations (black dashed). The HI mock with added foreground and systematics from the real MeerKLASS data, then cleaned as in equation (32) (red solid). The cross-correlation between HI-only and cleaned mock (orange solid). Lastly, the blue solid line is the autocorrelation of the cleaned mock but with signal loss reconstructed.</figcaption></figure>

mock-injected data and estimate the signal loss from the injected mock whose original signal is known.

We follow the formalism outlined in Cunnington et al. (2022), beginning with a HI mock, whose matrices we notify by M to distinguish from the real observation X. The mock is injected into the real observations i.e. $X_{obs} + M_{HI}$. This is then cleaned by projecting out a subset of modes defined by the mixing matrix, as in the previous sections, such that

$$M_{clean} = (X_{obs} + M_{HI}) - AA^{T}(X_{obs} + M_{HI}). \quad (32)$$

The data in $M_{clean}$ will now contain signal loss emulating that of the real data clean, along with residual systematics and noise from the real data.

We show the results of mocks contaminated by the real data in Fig. 11. We show measured power spectra averaged over all 500 mocks under different scenarios (details of power spectrum estimation are provided in the following Section 5). First, the grey dotted line shows the HI-only mocks, $M_{HI}$, with no foreground cleaning effects, representing the true power. The red solid line shows the result of adding the real observations to the mock, then PCA cleaning as in equation (32). The power is higher than the grey line, despite the inevitable signal loss the clean will have caused. This is evidence that the $M_{clean}$ maps contain residual systematics and thermal noise, which are additively biasing the power.<fnref n="9" /> We can demonstrate this by considering the orange solid line which is the cross-correlation of the cleaned mocks $M_{clean}$ with the HI-only mocks $M_{HI}$, a test only possible with simulations. Since the residual systematics and noise will be uncorrelated with the HI-only mocks, the distortion in power comes purely from signal loss, hence why it is lower than the HI-only power at all scales.

To isolate the contribution from unknown systematics, we need to account for thermal noise and signal loss from foreground cleaning. The thermal noise contribution can be approximated using

<footnote n="9">⁹There will also be a small contribution from the true HI signal in the observed data.</footnote>
equation (20) to generate 500 noise maps and add these to the Hi-only mocks. The black-dashed line shows this and demonstrates the boost in power relative to the grey-dotted line, more so at small scales where, due to the beam, the HI signal is low relative to the scale-independent noise.<fnref n="10" /> The last result presented in Fig. 11 is the reconstructed power (i.e. correcting for foreground cleaning signal loss) shown by the blue solid line. This utilizes the signal-loss evaluated from the cross-power (orange line) to reconstruct the Mclean autopower. This is discussed more in the following subsection. The discrepancy between the blue-solid and black-dashed line must be caused by systematics, since signal-loss and thermal noise have now been accounted for. We find a similar excess power with the real observations, which we discuss in the following Section 5. The results in Fig. 11 are to validate the realism of the mock pipeline, and also to motivate the signal loss reconstruction process, which we formalize in the next subsection.

## 4.5 Signal loss correction with transfer functions

We can utilize the method of mock signal injection from the previous Section 4.4 for estimating the signal loss in the power spectrum of the real data. This is now a well-understood process (Switzer et al. 2015; Cunnington et al. 2023) and used in previous single-dish intensity mapping analyses (Masui et al. 2013; Anderson et al. 2018; Cunningham et al. 2022; Wolz et al. 2022).

We begin with the Mclean term in equation (32), which as demonstrated, will contain signal loss that should emulate the real data. In some cases, the cleaned real data Xclean can be subtracted from Mclean. This accelerates the convergence of the transfer function calculation by removing the additional noise caused by the real data but has no bearing on the signal loss calculation. Leaving this piece in the calculation of Mclean is however crucial when using the transfer function for estimating errors in the power spectrum. We also find that due to the large number of mocks used in this work, an equivalent accuracy of the transfer function is obtained whether subtracting Xclean or not. For simplicity, we therefore only use the case as defined by equation (32), without the Xclean subtraction.

The transfer function is defined by the power spectrum, signified by the operator $P()$, of the cleaned mock Mclean with the foreground-free HI-only mock MHI, and the ratio of this with the power spectrum of the HI-only mock

$$ T_{\text{HI}}(\mathbf{k}) = \left\langle \frac{P(\text{M}_{\text{clean}}, \text{M}_{\text{HI}})}{P(\text{M}_{\text{HI}}, \text{M}_{\text{HI}})} \right\rangle_{N_{\text{mock}}}. \quad (33) $$

The angled brackets indicate an averaging over a large number of mocks such that a converged result is obtained. In this work, we use all mocks in our suite, thus $N_{\text{mock}} = 500$. The operator $P()$ can represent any power spectrum estimator i.e. we use this same recipe for estimating signal loss in the spherically averaged power and also in the cylindrical $(k_{\perp}, k_{\parallel})$ power spectrum, which we later utilize in some analysis.

In this work, we will mostly focus on the MeerKLASS deep field's cross-correlation with the overlapping GAMA galaxies. The signal loss from PCA cleaning in the cross-power spectrum should be identical to that in the autopower. This was established empirically in Cunningham et al. (2023) which showed how removing a subset of PCA modes entirely but leaving their contribution in the other

half of the cross-power product (e.g. in the galaxy maps which have no cleaning, thus no signal loss) generates no more power than the auto-HI. We present further reasoning for why this is the case in Appendix B3. To summarize, we believe it is explained by the PCA modes being orthogonal, thus maximally uncorrelated. Therefore, the additional remaining modes in the uncleaned field are linearly uncorrelated with what survives in the cleaned HI, producing no extra power.

Despite the approximately equal signal loss between auto and cross-correlations, we still define a separate transfer function for the cross-power spectrum purely for covariance estimation purposes whereby the only difference is we switch in the GAMA galaxy mocks M$_g$ instead which will emulate the signal strength and shot-noise of the GAMA galaxies, thus providing a more robust error estimate in the cross-power.

$$ T_{\text{HI,g}}(\mathbf{k}) = \left\langle \frac{P(\text{M}_{\text{clean}}, \text{M}_{\text{g}})}{P(\text{M}_{\text{HI}}, \text{M}_{\text{g}})} \right\rangle_{N_{\text{mock}}}. \quad (34) $$

Any measured power spectra should be divided by $T(k)$, thus the transfer function represents a reconstruction of signal at the power spectrum level.

# 5 POWER SPECTRUM ESTIMATION

In this section we present power spectra estimations from the foreground cleaned MeerKLASS L-band deep field. We first outline the pipeline for performing these measurements then present the results. By default, we use maps that have been foreground cleaned by removing $N_{\text{fg}} = 10$ PCA modes, as outlined in the previous section. However, we also examine the results from varying $N_{\text{fg}}$.

## 5.1 Transforming to a Cartesian grid

The first task to perform in order to measure a power spectrum is to transform the foreground cleaned intensity maps onto a Cartesian grid with comoving $h^{-1}$Mpc coordinates. This allows for a Fast Fourier transform to be calculated on the Cartesian array, from which the power at different Fourier bandpasses $k$, can be averaged and presented.

We discussed the transformation between sky maps and Cartesian grids already in Section 3.1 when we performed the inverse routine of sampling Cartesian cells from lognormal mocks into celestial sky map voxels that emulate the MeerKLASS intensity maps. We again utilize the recent work and code in Cunningham & Wolz (2024) and first define a grid of size $l_x, l_y, l_z = 632.2, 341.5, 231.5 h^{-1}$Mpc with $n_x, n_y, n_z = 88, 48, 168$ cells along each dimension. The choice of cells is a factor of 1.5 smaller than the number of map voxels in R.A., Dec., and $v$. Making the Cartesian grid slightly lower resolution than the map and using 10 sampling particles per map voxel ensures each Cartesian cell is well-sampled, avoiding empty cells.

The foreground cleaned maps $\delta T_{\text{HI}}(\theta, v)$ are transformed to this Cartesian grid, giving $\delta T_{\text{HI}}(\mathbf{x})$. The intensity mapping weights also follow an identical transformation process $w_{\text{HI}}(\theta, v) \rightarrow w_{\text{HI}}(\mathbf{x})$. The catalogue of GAMA galaxies can be directly gridded by transforming their celestial coordinates into Cartesian, then assigning them onto the Cartesian grid $n_g(\mathbf{x})$. For all sampling particles and galaxies we use a straightforward nearest grid point (NGP) assignment. Cunningham & Wolz (2024) presented extensions to this method for avoiding biased power and obtaining sub-percent accuracy. However, given the current signal-to-noise of our data at small scales, aliasing effects will be negligible, thus higher order interpolations

<footnote n="10">¹⁰The noise is not completely scale-independent because it has gone through the same regridding and reconvolution processes as the data which causes some subtle damping.</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 12. Tukey window functions with differing shape parameters α. For comparison, we also show the Blackman window used in previous work.</figcaption></figure>

and interlacing are overkill and serve only to increase compute demand. However, we correct for the dominant effect which is the field's smoothing by the (NGP) particle assignment. This requires deconvolving the field, achieved by dividing through by the NGP mass assignment window function in Fourier space (Jing 2005)

$$ \delta \tilde{T}_{\text{HI}}(\mathbf{k}) \rightarrow \frac{\delta \tilde{T}_{\text{HI}}(\mathbf{k})}{W_{\text{ngp}}(\mathbf{k})}, \qquad (35) $$

where the window function $W_{\text{ngp}}$ accounts for the sampling particles discretization into cells of size $H$. These are separable functions along each Cartesian dimension, thus modelled by a Fourier transform of a top-hat function, which becomes a sinc function; $\text{sinc}(x) = \sin(x)/x$ and therefore

$$ \tilde{W}_{\text{ngp}}(\mathbf{k}) = \text{sinc}\left(\frac{k_x H_x}{2}\right) \text{sinc}\left(\frac{k_y H_y}{2}\right) \text{sinc}\left(\frac{k_z H_z}{2}\right). \qquad (36) $$

Previous work (e.g. Switzer et al. 2013; Cunnington et al. 2022) employed a tapering along the radial direction of the regridded field with the primary purpose of reducing edge effects in the Fourier transform. Our modelling approach should however account for these effects since the model is also convolved with the window function (Blake 2019). We still found results were worse when omitting a radial tapering though, which we suspect is instead caused by edge-effects from the foreground clean (Alonso et al. 2015) which are suppressed by the tapering. We experimented with a less aggressive taper by using a Tukey window, also known as a cosine window [this is also explored in Carucci et al. (2024)]. This has a shape parameter $\alpha$ which is simply a top-hat with no apodization when $\alpha = 0$ and equivalent to a Hann window when $\alpha = 1$. We show these window functions in Fig. 12. By opting for less aggressive Tukey windows we still effectively suppress the extreme edges where foreground cleaning effects may be present but avoid unnecessary tapering on the central data. We will explore this further with future data in MeerKAT's UHF band where the larger useable frequency ranges may render this less of an issue. In this work we found optimal results when using $\alpha = 0.8$ for the Tukey window which is applied to the intensity map weights $w_{\text{HI}}(\mathbf{x})$. We do not apply this to the galaxy data as we found this increased the correlations between k-bins in the power spectrum. A tapering in the angular directions could also be applied, but again as argued above, since we convolve our model with the survey window function, the bulk of these effects should be controlled. There is also no requirement for angular tapering

to suppress any radial-only foreground-cleaning edge effects. We therefore chose to omit angular tapering from our analysis.

## 5.2 Power spectrum formalism

With the fields assigned to a Cartesian space and appropriately tapered, the power spectra can then be measured. The Fourier transformed field of the HI temperature fluctuation maps $\delta T_{\text{HI}}$ is given by

$$ \tilde{F}_{\text{HI}}(\mathbf{k}) = \sum_{x} \delta T_{\text{HI}}(\mathbf{x}) w_{\text{HI}}(\mathbf{x}) \exp(ik \cdot \mathbf{x}), \qquad (37) $$

where the weights, now denoted by $w_{\text{HI}}$ to distinguish them from the galaxy weights, still derive from the field of time-stamp counts, which are regridded following the same process as outlined in the previous subsection. The Fourier transformed field of the galaxy count field $n_g$ is given by

$$ \tilde{F}_{\text{g}}(\mathbf{k}) = \sum_{x} n_{\text{g}}(\mathbf{x}) w_{\text{g}}(\mathbf{x}) \exp(ik \cdot \mathbf{x}) - N_{\text{g}} \tilde{W}_{\text{g}}(\mathbf{k}), \qquad (38) $$

where $N_g = \sum x$ is the total number of galaxies in the optical maps and $\tilde{W}_g$ is the weighted Fourier transform of the selection function

$$ \tilde{W}_{\text{g}}(\mathbf{k}) = \sum_{x} W_{\text{g}}(\mathbf{x}) w_{\text{g}}(\mathbf{x}) \exp(ik \cdot \mathbf{x}). \qquad (39) $$

The galaxy window function $W_g$ is normally constructed by stacking mock galaxies that carefully follow the survey selection characteristics of the experiment to account for its incompleteness. As discussed previously, there are no GAMA G23 mocks available for this purpose so we therefore construct a simple window function which is 1 within the survey boundaries and 0 everywhere else. Some inclusion of the redshift evolution of $N(z)$ could have been included, but as Fig. 7 showed, the variation over our narrow redshift range is minimal. The window function is normalized so that $\sum_x W_g(\mathbf{x}) = 1$. Since we lack this proper galaxy window function we do not implement the optimal weighting as per Feldman, Kaiser & Peacock (1994), but instead use uniform weighting, so equivalent of assuming $w_g(\mathbf{x}) \equiv W_g(\mathbf{x})$.

The cross-power spectrum estimator is then given by

$$ \hat{P}_{\text{HI,g}}(\mathbf{k}) = \frac{V_{\text{cell}}}{\sum_{x} w_{\text{HI}}(\mathbf{x}) w_{\text{g}}(\mathbf{x}) W_{\text{g}}(\mathbf{x})} \text{Re}\left\{\tilde{F}_{\text{HI}}(\mathbf{k}) \tilde{F}_{\text{g}}^{*}(\mathbf{k})\right\} \frac{1}{N_g}, \qquad (40) $$

where $V_{\text{cell}} = l_x l_y l_z / (n_x n_y n_z)$. The hat notation ($\hat{\cdot}$) indicates that this is the estimated power from the data, to distinguish from the modelled power, which we discuss shortly. The HI autopower spectrum estimator is more straightforwardly given by

$$ \hat{P}_{\text{HI}}(\mathbf{k}) = \frac{V_{\text{cell}}}{\sum_{x} w_{\text{HI}}^{2}(\mathbf{x})} |\tilde{F}_{\text{HI}}(\mathbf{k})|^2. \qquad (41) $$

These power spectra are either spherically averaged into bandpowers $|k| \equiv k$ to provide 1D power spectrum results or alternatively cylindrically averaged into 2D bandpowers defined by modes perpendicular ($k_{\perp}$) and parallel ($k_{\parallel}$) to the radial direction $z$, which we assume is the line of sight in a plane-parallel approximation, valid for the current MeerKLASS sky sizes. They are finally divided by their relevant foreground transfer functions [equation (33) and equation (34)] to reconstruct signal loss for each k-bin.

### 5.2.1 Modelling the power spectrum

Here we outline the model of the power spectrum estimations for which we mostly follow the work in Blake (2019), with some extensions. We model the HI power spectrum using the equation already
introduced in equation (14) when discussing mock simulations. We repeat the model here for convenience;

$$P_{\mathrm{H1}}(\mathbf{k}) \equiv P_{\mathrm{H1}}(k, \mu) = \frac{\overline{T}_{\mathrm{H1}}^2 (b_{\mathrm{H1}} + f \mu^2)^2 P_{\mathrm{m}}(k)}{1 + (k \mu \sigma_v / H_0)^2}. \quad (42)$$

Additional modelling of observational and mapping effects is also necessary since these will distort the power. To model the MeerKAT beam and the additional smoothing from the reconvolution process we assume the field has undergone a constant smoothing by a 2D Gaussian kernel whose standard deviation is given by $R_{\mathrm{beam}} = \sqrt{\gamma} R_{\mathrm{beam}}^{\mathrm{max}}$ in comoving units. $\gamma = 1.4$ accounts for the reconvolution to a common effective resolution and $R_{\mathrm{beam}}^{\mathrm{max}}$ is the Gaussian approximation of the MeerKAT beam size at the deepest redshift (and lowest frequency $\nu_{\mathrm{min}}$)

$$R_{\mathrm{beam}}^{\mathrm{max}} = \chi(z_{\mathrm{max}}) \sigma_{\mathrm{beam}}(\nu_{\mathrm{min}}), \quad (43)$$

where $\sigma_{\mathrm{beam}}$ is given as before from equation (17) and $\chi(z)$ is the comoving distance to redshift $z$, which in this case, is the maximum redshift for our frequency range. The perpendicular modes will then be damped according to

$$\tilde{B}_{\mathrm{beam}}(\mathbf{k}) = \exp\left[\frac{-k^2(1-\mu^2)R_{\mathrm{beam}}^2}{2}\right]. \quad (44)$$

For this work, we find $R_{\mathrm{beam}} = 14.1 \, h^{-1} \mathrm{Mpc}$.

Similar to the beam, there will also be distortions to the power caused by the finite resolution of the intensity maps. We approximate the pixel size in comoving units by $s_{\mathrm{pix}} = \chi(z_{\mathrm{eff}})\Delta\theta$ where $\Delta\theta = 0.3$ deg is the pixel size used for the MeerKLASS L-band mapping. The frequency channel depth is approximated by taking the mean across all spacings between redshifts corresponding to channel boundaries i.e. $s_{\mathrm{chan}} = \sum_i [\chi(z_i) - \chi(z_{i+1})]/N_{\mathrm{chan}}$. These quantities are used in damping functions to model the suppressed resolution from mapping into pixels and frequency channels, the former being the far more dominant effect for the wide-pixel MeerKLASS intensity maps. These functions are given as

$$\tilde{B}_{\mathrm{pix}}(\mathbf{k}) = \mathrm{sinc}\left(\frac{k\sqrt{1-\mu^2}s_{\mathrm{pix}}}{2}\right), \quad (45)$$

and

$$\tilde{B}_{\mathrm{chan}}(\mathbf{k}) = \mathrm{sinc}\left(\frac{k \mu s_{\mathrm{chan}}}{2}\right). \quad (46)$$

The model and damping functions are discretized onto the same 3D grid of k-modes as the FFT for the real data, then convolved with the survey window function. The model for the observed H1 autopower spectrum therefore becomes

$$P_{\mathrm{H1}}^{\mathrm{obs}}(\mathbf{k}) = \frac{W_{\mathrm{H1}}^2(\mathbf{k}) \circledast [P_{\mathrm{H1}}(\mathbf{k}) B_{\mathrm{beam}}^2(\mathbf{k}) B_{\mathrm{pix}}^2(\mathbf{k}) B_{\mathrm{chan}}^2(\mathbf{k})]}{\sum_x w_{\mathrm{H1}}^2(\mathbf{x}) W_{\mathrm{H1}}^2(\mathbf{x})}, \quad (47)$$

where we have normalized by the sum of the weights and window function, $W_{\mathrm{H1}}$, to obtain an amplitude that matches the weighted and windowed estimator. For the intensity mapping window function, we assume a simple binary function; 1 where we have scan coverage and 0 everywhere else. There is scope for improving this approximation and further work is needed to examine how intensity mapping surveys are impacted by any selection effects.

A similar approach is followed for modelling the cross-power with GAMA. The cross-correlation power spectrum is given by

$$P_{\mathrm{H1,g}}(\mathbf{k}) = \frac{\overline{T}_{\mathrm{H1}}[rb_{\mathrm{H1}}b_{\mathrm{g}} + b_{\mathrm{H1}}f\mu^2 + b_{\mathrm{g}}f\mu^2 + f^2\mu^4]P_{\mathrm{m}}(k)}{1 + (k\mu\sigma_v/H_0)^2}, \quad (48)$$

<table><caption>**Table 1.** Fiducial parameters assumed in model fitting.</caption><thead><tr><th>Parameter</th><th>b<sub>HI</sub></th><th>b<sub>g</sub></th><th>r</th><th>f</th><th>σ<sub>v</sub></th></tr></thead><tbody><tr><td>Fiducial value</td><td>1.13</td><td>1.90</td><td>0.90</td><td>0.74</td><td>400 km s<sup>-1</sup></td></tr></tbody></table>

where we introduce $r$, the correlation coefficient between the H1 fluctuations and the GAMA galaxy overdensities. This accounts for the different astrophysics in each tracer and non-linear clustering, which can mean the two fields are not exact multiples of each other and there is instead some stochasticity which will suppress their cross-power (see e.g. Moradinezhad Dizgah et al. 2022; Obuljen et al. 2023; Sato-Polito, Kokron & Bernal 2023, for applications in the context of line-intensity mapping). We assume $r$ to be constant, which is sufficient for the current signal to noise, but in reality it will exhibit some scale dependence. The correlation coefficient only enters on the single term in equation (48) because it only influences the biased tracer fields, not their *unbiased* velocity fields which are modulated by $f$. The observational cross-power is then similarly convolved with the survey window functions but only requires single factors of the damping functions $B$;

$$P_{\mathrm{H1,g}}^{\mathrm{obs}}(\mathbf{k}) = \frac{W_{\mathrm{H1}}(\mathbf{k}) W_{\mathrm{g}}^{*}(\mathbf{k}) \circledast [P_{\mathrm{H1,g}}(\mathbf{k}) B_{\mathrm{beam}}(\mathbf{k}) B_{\mathrm{pix}}(\mathbf{k}) B_{\mathrm{chan}}(\mathbf{k})]}{\sum_{\mathbf{x}} w_{\mathrm{H1}}(\mathbf{x}) w_{\mathrm{g}}(\mathbf{x}) W_{\mathrm{H1}}(\mathbf{x}) W_{\mathrm{g}}(\mathbf{x})}. \quad (49)$$

We do not include modelling of uncertainties in the GAMA galaxy positioning, such as redshift estimation. The GAMA radial kernel equivalent of equation (46) will likely have a larger impact due to spectroscopic redshift estimation errors. However, this will still be a subdominant contribution to the overall modelling relative to the beam and intensity map pixelization, thus it is reasonable to neglect observational effects from the GAMA galaxies.

We perform a simple least-squares fit for an overall amplitude to our cross-power, equivalent to assuming fixed values for all parameters except for a fitted $\overline{T}_{\mathrm{H1}}$. This is deliberately simplistic to keep the model and fitting in a low-dimensional form. The small overlap sample in the cross-correlation means the statistical noise from such measurements is not yet suited to high-dimensional modelling. We are not therefore extending much from previous work (e.g. Cunningham et al. 2022; Wolz et al. 2022) in this regard. Furthermore, the auto-H1 still has evidence of additive bias, as we will discuss, and therefore can only be trusted for upper limits in parameter inference.

For the power amplitude fitting, we assume the parameter values outlined by Table 1. The H1 bias comes from interpolated results from hydrodynamical simulations (Villaescusa-Navarro et al. 2018), which can be approximated by the polynomial

$$b_{\mathrm{H1}}(z) = 0.842 + 0.693z - 0.046z^2. \quad (50)$$

The galaxy bias is obtained by fitting the amplitude of the GAMA autopower spectrum. The growth rate is obtained from $f = \Omega_{\mathrm{m}}(z_{\mathrm{eff}})^{0.545}$ for our Planck 18 cosmology. The correlation coefficient $r$ and velocity dispersion $\sigma_v$ are poorly constrained nuisance parameters tuned by eye to provide an adequate fit. We discuss model and parameter fitting further in Section C where we present results from a simultaneous MCMC fit on several parameters, demonstrating the degeneracies between them and discuss plans for breaking such degeneracies, but this is not a central focus of this work.
## 5.3 Error estimation with the transfer function

In previous analysis (Cunnington et al. 2022), an analytical approach was primarily used for estimating errors on the power spectrum measurements. However, this approach makes various assumptions, mainly that thermal noise dominates the error budget and that Gaussian statistics are valid. Since we are now reaching regions in k-space where thermal noise no longer dominates with MeerKLASS, we are required to relax these assumptions. The assumption of Gaussianity is arguably not solid either, given the additive systematics we know to still be present within the data. Furthermore, an extra contribution from variance in the PCA foreground cleaning process is not accounted for in this analytical error estimation.

A more robust approach to estimating the errors is to use the transfer function scatter across the different mock realizations $T_i(k)$. This approach has been used in previous analysis (e.g. Anderson et al. 2018) and has also been validated in simulations (Cunnington et al. 2023). This error will estimate the contribution from sample variance, galaxy shot-noise, foreground and systematic residual, and variance from signal loss caused by eigenmode variation. The contribution from residual foreground and systematics to the error using this transfer function has not been robustly validated and is difficult to do so. Whilst they are present in each transfer function realization they are not varying; only the response to them from each iteration's mock signal varies. The ideal way to robustly investigate this would be to simulate all residual systematics and foregrounds, but as we have discussed elsewhere, this is not currently possible for any Hi intensity mapping experiments, given their complexity and novelty. We defer this detailed investigation to future work and, for now, are content with our approach. As we will show, our estimated errors, if anything, appear conservative based on some analysis and are therefore not a major concern for the confidence gleaned from our conclusions. We believe this is because it is the variance from signal loss caused by eigenmode variation that is the most dominant contribution to the error with our current data. We concluded this by comparing the variance between mock injection either before or after the foreground cleaning. When mocks are injected before foreground cleaning, the variance in the power spectrum is boosted by nearly an order of magnitude across most scales. This highlights the importance of optimizing contaminant subtraction which is the focus of Carucci et al. (2024). For it to perform effectively, the mocks in the transfer function calculation must emulate the real data effectively, which we have strived to do so, as outlined in Section 3.

The transfer function for one single realization is defined by

$$T_i^{\text{Hi,g}}(k) = \frac{\mathcal{P}(\text{M}_{\text{clean}}, \text{M}_g)}{\mathcal{P}(\text{M}_{\text{Hi}}, \text{M}_g)}, \qquad (51)$$

i.e. equivalent to equation (34) but without the averaging over $N_{\text{mock}}$. To reconstruct signal in a measured power spectrum, we take the ensemble average of equation (34) i.e. $T(k) = \langle T_i(k) \rangle_{N_{\text{mock}}}$, and divide the estimated power by this; $\hat{P}_{\text{rec}}(k) = \hat{P}(k)/T(k)$. Note we are dropping the subscript Hi,g for brevity and also because the formalism we present here is identically applicable to the Hi-auto. The assumption made when using the transfer function for error estimation is that the scatter in the reconstructed power comes from the transfer function distribution $T(k) + \delta T_i(k)$ contributed by each $i$th mock. We can isolate the scatter component from the array of calculated transfer functions by simply subtracting the ensemble mean for each $i$th mock

$$\delta T_i(k) = T_i(k) - T(k). \qquad (52)$$

Assuming any errors add in quadrature, error propagation suggests

$$\frac{\delta \hat{P}_i^{\text{rec}}(k)}{\hat{P}_{\text{rec}}(k)} = \sqrt{\left(\frac{\delta \hat{P}_i(k)}{\hat{P}(k)}\right)^2 + \left(\frac{\delta T_i(k)}{T(k)}\right)^2}, \qquad (53)$$

however, we assume all quantifiable error lies in the transfer function distribution, hence $\delta \hat{P}_i = 0$ and so

$$\delta \hat{P}_i^{\text{rec}}(k) = \hat{P}_{\text{rec}}(k) \frac{\delta T_i(k)}{T(k)}. \qquad (54)$$

Our setup thus provides $N_{\text{mock}}$ realizations of $\delta \hat{P}_i^{\text{rec}}$ and the covariance can be computed from this for any $k$-bin. We reiterate that this formalism is valid whichever way the $k$-modes are binned e.g. for a 2D cylindrical averaged or 1D spherically averaged power. In this work, our error bars in any 1D power spectrum estimate represent the limits of the 68th percentile regions of the distribution of $\delta \hat{P}_i^{\text{rec}}(k)$. If the uncertainties are Gaussian symmetric, this is equivalent to $1\sigma$ error bars, however, as we will later show, this is not perfectly the case for all $k$-modes.

## 5.4 Spherical averaging $k$-cuts

Before presenting spherically averaged power spectra from the MeerKLASS data, we first present results from the cross-correlation with GAMA galaxies in cylindrical-(k⊥, k‖) space. This is shown in Fig. 13 and motivates a discussion on where in k-space we have the most signal-to-noise so we can make optimal spherically averaged $\hat{P}(k)$ measurements. The dashed grey lines on all panels represent the |k|-bin boundaries we use for the spherically averaged power spectra, restricting this to $k < 0.3\ h\ \text{Mpc}^{-1}$ to avoid highly non-linear cosmology and increasingly poor signal-to-noise. The orange solid lines represent boundaries for potential $k$-cuts which we will discuss.

The first panel of Fig. 13 shows the cross-power between GAMA G23 galaxies and the MeerKLASS deep field, cleaned by removing $N_{\text{fg}} = 10$ PCA modes. The $k_{\perp}k_{\parallel}$ scaling is for demonstrative purposes so the high power at low-k does not saturate the plot. At $k \gtrsim 0.15\ h\ \text{Mpc}^{-1}$ the large MeerKAT beam begins to have an impact and will damp Hi fluctuations, which explains why the cross-power signal tends to zero at high-k⊥. Except for the extreme low-k regions, the rest of the k-space reveals a reasonably consistent positive cross-correlation signal, providing strong evidence of consistent clustering between the MeerKLASS maps and GAMA galaxies. As discussed in the previous Section 4.5, the signal loss in this power has been corrected with the foreground transfer function. The second panel shows the corresponding transfer function for the $N_{\text{fg}} = 10$ clean. Low values of $T_{\text{Hi,g}}$ signify high signal loss. This is evidence that nearly all signal is lost for the lowest-k‖ modes which is expected since these modes should be highly degenerate with the removed long radial modes of the continuum foregrounds. The high signal loss at low-k⊥ is interesting and not seen in previous MeerKLASS work (Cunnington et al. 2022) with a similar method, albeit with less robust mocks. This is likely caused by the small GAMA footprint (see Fig. 4) which means the large-ranging modes at these low-k are not adequately resolved and deliver little signal which is easily coupled to foreground modes and stripped away in the clean. This can be explored further in follow-up work with wider surveys where we expect this not to be a problem.

Due to the noisy power at low-k and damped power at high-k⊥ in Fig. 13, we can choose to cut these regions of k-space from our spherically averaged power to obtain a more robust result. The orange solid lines represent chosen boundaries, where modes beyond these boundaries are excluded in the averaging of each k-bin (shown by
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Cylindrical ($k_{\perp}, k_{\|}$) analysis of power, signal-loss, and signal-to-noise. The first panel from the left shows the cross-power between the MeerKLASS deep field and GAMA galaxies. The second panel is the transfer function used to reconstruct the power in $\hat{P}_{\text{HI,g}}$. The third panel shows the signal-to-noise by dividing the cross-power by the error estimation discussed in the text. The final fourth panel counts the number of discrete $k$ wavenumber in the Fourier grid, demonstrating the concentration of available modes. The grey contours mark the boundaries of the $k$-bins used in spherically averaged power $P(k)$. The orange lines indicate regions used for $k$-cuts discussed in the text.</figcaption></figure>

the dashed grey lines). Providing the modelling of any power spectra also incorporates the $k$-cuts, then it will not bias any results. The $k$-cuts approach is further motivated by considering the signal-to-noise of the cross-power which we show in the third panel of Fig. 13. Here the errors $\sigma(\hat{P}_{\text{HI,g}})$ are computed using the transfer function by taking the standard deviation from the $\delta\hat{P}_{\|}^{\text{rec}}(k_{\perp}, k_{\|})$ distribution (see previous Section 5.3). This result provides a strikingly clear window where signal-to-noise is optimal; the orange to red regions in the third panel. The orange $k$-cut boundary lines have been chosen to target this region. One limitation of this approach is that it excludes the regions of $k$-space with the most modes, as shown by the fourth and final panel. We will present comparisons between the spherically averaged $P(k)$ with and without these $k$-cuts in the upcoming Section 5.6.

## 5.5 HI autocorrelation power spectrum

We next examine results for the autocorrelation HI power spectrum from the MeerKLASS L-band deep field. Previous single-dish intensity mapping results have been thermal noise-dominated across all scales. We are now crossing an important threshold with MeerKLASS where the HI fluctuations should be dominating the noise at large scales. We demonstrate this in Fig. 14, where we show the estimated HI power spectrum from the MeerKLASS deep field (blue data points) and overlay an estimate of the thermal noise (red solid line). As outlined in Section 3.2, this noise is modelled by generating Gaussian random fields with the expected rms given by equation (20) which utilizes the observational hit map. Despite being the HI autopower, we still restrict the measurement to the GAMA-field footprint, for consistency with the cross-correlation results. We tested using the full field but found no improvement in the autopower. We leave exploration of the full-field to follow-up

<figure><img src="image_15.png" /><figcaption>**Figure 14.** The HI autopower spectrum for the MeerKLASS L-band deep-field foreground cleaned by removing $N_{\text{fg}} = 10$ PCA modes. The red solid line is the noise estimate from measuring the average power of 10 Gaussian random fields with rms given by $\sigma_{\text{N}}$ (equation 20). The grey dotted line is the observational HI power spectrum model (equation 47) and the black dashed line is the combination of this HI model plus the thermal noise estimate.</figcaption></figure>

work which focuses on the HI-auto power spectrum. The results in Fig. 14 include the scale $k$-cuts discussed in the previous section.

Fig. 14 also shows the expected HI power (grey dotted line) from equation (47). For this auto-HI model, we do not fit the amplitude to avoid additive bias from residual systematics over inflating the model estimate. Instead we use the scale-independent $\bar{T}_{\text{HI}}$ fit from
the upcoming cross-power measurement. This model shows how the cosmological HI power surpasses the thermal noise on large scales ($k \lesssim 0.15 h \text{ Mpc}^{-1}$). The thermal noise will be present even in the absence of clustering, thus the total expected power can be reasonably approximated as a simple addition of the HI power, plus the thermal noise (shown by the black dashed line).

The small discrepancy between the deep-field power estimate and the model represents the lowest upper limit from a HI autocorrelation on these large scales. We find the mean residual factor between model and data to be just 2.36 across all k. Whilst the model and data are close in agreement, inferred parameter values for the bias and mean HI temperature from fitting the amplitude of the autopower yields inconsistent values with the GAMA cross-correlation (discussed next). This is therefore evidence of low-level residual systematics that are additively biasing the power. In this evaluation, we inherently assume that the transfer function reconstructs any signal loss in the Gaussian thermal noise. If signal loss in the Gaussian noise is underestimated, the residual between the model and data may be higher than we report here. We therefore intend to robustly test this assumption when pursuing a more dedicated study of the HI autocorrelation.

In follow-up work, we will focus on improving this upper limit from the HI autopower which can be done in two ways: First, by implementing some corrections to remove low-level systematic contributions such as the residual RFI structures discussed in Section 2.4.2. Secondly, we can split the observations into subsets i.e. combinations of different maps from different scan times and/or dishes. This is useful for isolating dish- and time-dependent systematics and removing their biasing contributions to the power. Since each scan and dish is observing the same patch of sky, we would ideally expect all maps to agree, thus any differences will represent thermal noise, or systematics such as residual RFI, varying polarization leakage etc. We do not pursue these ideas in this paper since they warrant their own detailed investigation.

## 5.6 Cross-correlation power with GAMA

We now focus on the cross-correlation with the GAMA galaxies for the remainder of the paper. In Fig. 15 we present the cross-correlation power spectrum between the GAMA galaxies and the MeerKLASS deep field. All points are positive so none are lost to the log scaling, indicating a strong, positive correlation. Generally, a good agreement with the model is obtained. We performed a null test on the result by substituting the GAMA mock galaxies (from Section 3.3) in place of the real galaxies. The mocks should not correlate with the intensity maps, and indeed, the resulting cross-power converges to zero once averaged over a sufficient number of mocks. The fitted amplitude of the model is equivalent to an inferred parameter of $\bar{T}_{\text{HI}} = 0.166 \text{ mK}$, with all other values in Table 1 assumed. We discuss this further in Appendix C, but do not focus on parameter inference in detail due to the large statistical noise owing to the small overlap footprint between MeerKLASS and GAMA.

The error bars for the cross-power in Fig. 15 are obtained using the scatter from the foreground transfer function as discussed in Section 5.3. The error size appears to be overpredicted based on the scatter of the data points, which have a reasonably tight agreement with the model relative to what the error bars reflect. This is supported by examining the reduced $\chi^2$ which we find to be $\chi^2_{\text{dof}} \equiv \chi^2/N_{\text{dof}} = 0.42$. A $\chi^2_{\text{dof}} < 1$ is conventionally indicative of overinflated errors. Using the reduced $\chi^2$ however is not a reliable test where there are non-linear contributions, as is likely the case for these observations that are still prone to residual systematics. In

<figure><img src="image_16.png" /><figcaption>**Figure 15.** The cross-correlation power spectrum between GAMA galaxies and the MeerKLASS L-band deep field, foreground cleaned by removing $N_{\text{fg}} = 10$ PCA modes. The black dashed line represents a model given by equation (49) with a fitted scale-independent amplitude.</figcaption></figure>

these cases, the number of degrees of freedom cannot be robustly determined (Andrae, Schulze-Hartung & Melchior 2010). For our case we simply assume $N_{\text{dof}} = 13$, for our 14 k-bins minus a single fitted degree of freedom, representing the amplitude fit. Furthermore, noise within the data causes high uncertainty in $\chi^2$ estimates, which despite the relatively low noise levels we are now reaching, may still be a concern. Alternatives to $\chi^2$ have been proposed, and specifically for 21 cm applications (e.g. Tauscher, Rapetti & Burns 2018), but these have not been pursued in this work. We only use the $\chi^2$ as an approximate guide, mainly for relative comparisons between different cases.

For the estimation of $\chi^2$, we include the full contribution from off-diagonals in the k-bin covariance C, which can be trivially estimated for our error estimation method, simply the covariance over $\hat{P}_{i}^{\text{rec}}(k) = \hat{P}_{\text{rec}}(k) + \delta \hat{P}_{i}^{\text{rec}}(k)$ (see Section 5.3). The $\chi^2$ is then estimated from

$$ \chi^2 = (\hat{p} - p_{\text{mod}})^{\text{T}} C^{-1} (\hat{p} - p_{\text{mod}}), \quad (55) $$

where $\hat{p}$ is the data vector and $p_{\text{mod}}$ the corresponding model. We present the normalized covariance matrix (or correlation matrix) in Fig. 16, showing evidence that the k-bins are slightly correlated, especially at low-k. This could be contributing to the large error bars, hence why we include the full covariance in the $\chi^2$ estimation.

To investigate the errors from the cross-power further, we plot the full distributions for the reconstructed power spectra $\hat{P}_{i}^{\text{rec}}$ from each i'th mock in the foreground transfer function computation. Fig. 17 shows this for each k-bin, including an additional low-k bin which we otherwise exclude due to the large error on these scales. We centre the distributions on the model and normalize by it, to aid comparison with the case of no k-cuts (red results), which includes the areas of k-space with low signal-to-noise (discussed in Section 5.4). The plot shows how the distributions have some non-Gaussian features, which could also impact error estimation, along with the correlated k-bins shown by the covariance. We find the average kurtosis across all k-bins for the k-cut case to be -0.476, indicating the distribution is flatter than a normal distribution with more contribution in the
<figure><img src="image_17.png" /><figcaption>**Figure 16.** k-bin correlation matrix for the cross-power between the MeerK-LASS deep field and GAMA galaxies, shown in Fig. 15.</figcaption></figure>

tails. This could also widen the central percentiles and thus increase the error bar size.

Fig. 17 shows a mild improvement in performance for the case with *k*-cuts. The distributions around the model are less broad at low-*k*, which is expected since this was one of the regions where signal-to-noise was poorest (see third panel of Fig. 13). In terms of the Gaussianity of the distributions, we find a better kurtosis score for the case with *k*-cuts, suggesting the profiles are more Gaussian. In contrast, we find the results *without* *k*-cuts to be less skewed, although this is marginal and in general, the skewness of the distributions is not concerning. Any improvement coming from the *k*-cuts is small and we also find little difference in detection significance, as we will explore shortly. As the volume of data grows beyond this small overlap sample, the *k*-cut results may prove more optimal since we will no longer be cutting modes from an already sample variance-limited data set. This is something that future data sets with MeerKLASS will reveal.

Perfect Gaussian distributions in Fig. 17 is an unrealistic expectation, even in the high-precision limit. This is explored in Wilensky, Brown & Hazelton (2023), albeit for 21 cm epoch of reionization analyses. They highlight how in cross-multiplication estimators such as our galaxy cross-correlation, implicit assumptions made in the central limit theorem have consequences for the Gaussianity of error distributions. However, our distributions in Fig. 17 serve mainly as

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Estimated power distribution for each *k*-bin from the cross-correlation between MeerKLASS intensity maps and GAMA. This is computed using the power distribution in the 500 mocks used in the foreground transfer function. *k*-bin value is given in the top-left corner of each panel. The power distributions are centred on the fitted model at each *k* so a zero median in the distribution represents perfect accuracy. Shaded regions show the 68th percentiles for the distributions which are used as the limits for the 1 σ error bars in the power spectra plots.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Variation of cross-correlation results with changing $N_{fg}$ i.e. the number of PCA modes removed in the foreground clean. The left (black) axis shows the $\chi^2$, the lower values generally represent a better fit between cross-power data and model. The right (red) axis shows the cross-correlation power spectrum detection significance. Solid and thin-dashed lines represent the case with and without k-cuts respectively.</figcaption></figure>

a guide where approximate Gaussianity is sought to demonstrate reasonable robustness in our data and error estimation procedure. This is supported by Wilensky et al. (2023), which concludes that extreme violations of the central limit theorem are required to break the approximate Gaussianity in error distributions. This will require careful consideration in future precision cosmology involving parameter inference. However, for current pathfinding applications, the Gaussianity and noise within the distributions of Fig. 17 provide deeper context into error estimation, without overly stretching central limit theorem assumptions.

### 5.6.1 Varying the foreground clean

So far we have only considered the chosen case of removing $N_{fg} = 10$ PCA modes for the foreground clean. Fig. 18 shows how the $\chi^2$ (which we continue to use as a rough goodness-of-fit metric) varies with $N_{fg}$. As seen in previous work (Cunnington et al. 2022), too low an $N_{fg}$ does not sufficiently clean the foregrounds, causing a poor fit. Going too high with $N_{fg}$ also diminishes the fit. Ideally, we would expect high $N_{fg}$ to cause more signal loss which should be reconstructed by the transfer function, giving unbiased results and thus not diminishing the fit. A higher uncertainty could be expected with greater signal loss, but this should be accounted for by the increasing error bars i.e. the $\chi^2$ should stabilize and not start to rise again at high $N_{fg}$. On further scrutiny, we found that it is a single data point at $k \sim 0.2 \, h \, \text{Mpc}^{-1}$ with very little power distorting the fits at higher $N_{fg}$. This is evidence of low-level residual systematics that become more influential as signal is stripped away in the more aggressive foreground clean. Therefore, with future results where systematics are more controlled, this is expected to be less of a problem. For the case with no k-cuts (thin dashed lines), there is a plateau in $\chi^2$ at higher $N_{fg}$, as we argued should be the expected situation. This is because the no k-cut power spectra were less impacted by the $k \sim 0.2 \, h \, \text{Mpc}^{-1}$ anomaly at high $N_{fg}$ which is potentially suggesting that the k-cuts are overly aggressive for this scale. With larger overlapping fields, the more accessible k-bins will naturally mitigate this discrepancy. The no k-cut results in Fig. 18 are therefore encouraging that a model and data fit should reach a plateau for high $N_{fg}$ as expected with the current analysis pipeline.

The right hand side (red) y-axis in Fig. 18, also shows the detection significance defined by $\sqrt{\Delta\chi^2} \equiv \sqrt{\chi^2 - \chi_{\text{null}}^2}$, where $\chi_{\text{null}}^2$

is the result's $\chi^2$ with a null model with zero power. The cross-correlation with GAMA is returning a modest $> 4\sigma$ detection for the most optimal choices of foreground clean, lower than previous MeerKLASS detections (Cunnington et al. 2022). This is however still impressive, given that galaxies only overlap with approximately a quarter of the full MeerKLASS deep field. The detection significance is also restricted by the more robust approach to error estimation, which as we have already discussed, provides conservative estimates potentially inflated by residual systematics and correlated k-bins. To demonstrate this point, we artificially construct a scenario where a $\chi^2_{\text{dof}} \sim 1$ is obtained for the $N_{fg} = 10$ case. This is achieved by assuming a perfectly diagonal covariance and decreasing the diagonal errors by 40 per cent. This boosts the detection significance to $9.8 \, \sigma$, verifying that it is the conservative approach to errors lowering the detection significance, rather than the data-model residual. This is of course only a toy demonstration but it validates the argument that if large errors, which are suspected to be driven by systematics and k-bin correlations, can be contained, then detection robustness should improve. This is expected in future work where we will pursue dedicated studies into systematic mitigation. Some mode coupling is always likely, driven at least by non-linear clustering. However, we are confident that improvements can be made in covariance. MeerKLASS observations ongoing now are covering wider areas with more overlap with Stage-IV optical galaxy surveys. These observations are also using MeerKAT's UHF-band, which will allow a deeper usable redshift range, extending the $0.39 < z < 0.46$ used in this work to approximately $0.4 < z < 1.4$. This will naturally improve the large statistical noise due to the small volume overlap in this work.

## 6 STACKING

In this section, we explore the possibility of a stacking detection of the HI signal in the MeerKLASS L-band deep-field intensity map. The technique of spectral line stacking (Zwaan, van Dokkum & Verheijen 2001) can be used to detect the HI content of galaxies below noise limit. In combination with spectroscopic galaxy surveys, spectral line stacking probes the average emission line profile of selected galaxy samples. Complementary to the intensity mapping technique, the stacking measurement of the emission line profile probes the HI mass function of galaxy samples (e.g. Bera et al. 2022), the relation of HI content with its host environment (e.g. Guo et al. 2020; Bera et al. 2023; Dev et al. 2023; Rhee et al. 2023), and the baryonic Tully-Fisher relation (e.g. Meyer et al. 2016). Interferometric observations of MeerKAT using the MeerKAT International GigaHertz Tiered Extragalactic Exploration survey (MIGHTEE; Jarvis et al. 2016; Maddox et al. 2021) have been used to produce HI measurements with regard to the scaling relation (Sinigaglia et al. 2022) with galaxy properties and the dependence on the large-scale structure environment (Sinigaglia et al. 2024). See also Ponomareva et al. (2021), Tudorache et al. (2022), Pan et al. (2023), and Ponomareva et al. (2023) for direct detections of HI galaxies in the MIGHTEE survey.

For intensity mapping experiments, stacking analysis has been performed on the HI intensity maps using the Parkes telescope (Anderson et al. 2018; Tramonte et al. 2019; Tramonte & Ma 2020) and the CHIME telescope (CHIME Collaboration 2023), as well as on CO intensity maps from the CO Mapping Array Project (Dunne et al. 2024) and Ly $\alpha$ emission of [O III]-selected galaxies from the Hobby-Eberly Telescope Dark Energy Experiment (Lujan Niemeyer et al. 2022a, b).
MeerKLASS, on the other hand, faces unique challenges. The ~1 deg instrument power beam introduces heavy contamination from nearby sources. For the GAMA galaxy sample used, the number density of the sources is ~35 deg⁻² when integrated along the frequency subband. Therefore, when stacking each source, the stacking procedure will inevitably include emissions from more than one GAMA galaxy. The velocity width of the HI sources further scatter the 21 cm emission into nearby channels. The double-counting of sources from the beam creates a plateau of excess signal as shown in Fig. D2, and can be resolved by PCA cleaning, which we discuss in Appendix D. However, the residual after PCA cleaning also incorporates systematic effects. The chromatic structure of the systematic effects produces oscillations in the stacked profile, with their amplitudes comparable to the signal. For this work, we are only interested in using stacking as part of the validation test. As discussed later, we obtain a detection of a stacked signal in both angular and spectral space, which can be validated by performing a null test with random galaxy positions. The existence of a stacked signal further demonstrates that we have reached the depth needed for the possible detection of the HI power spectrum using the MeerKLASS L-band intensity maps. We leave a more detailed study into refining the stacking technique for probing HI science to follow-up work.

## 6.1 The stacking procedure of the intensity maps

We use the PCA-cleaned map described in Section 4 and the positions of the GAMA galaxies to perform a 3D stacking in $(\delta\alpha, \delta\phi, \delta\nu)$ space so that

$$I(\delta\alpha, \delta\phi, \delta\nu) = \frac{\sum_i I(\alpha_i + \delta\alpha, \phi_i + \delta\phi, \nu_i + \delta\nu)w_i}{\sum_i w_i}, \quad (56)$$

where $(\alpha_i, \phi_i)$ are the right ascension and declination of the map pixel in which the i-th galaxy resides, $I$ is the flux density, and $w_i = w_{HI}(\alpha_i + \delta\alpha, \phi_i + \delta\phi, \nu_i + \delta\nu)$ is the reconvolved weight of each pixel defined in Section 4.1. We choose the angular and the spectral bins to be consistent with the map resolution with $\Delta\alpha = \Delta\phi = 0.3$ deg. The frequency resolution corresponds to an effective velocity resolution of $\Delta\nu = 62.83$ km s⁻¹ at the centre frequency of the subband 997.37 MHz. Note that no Doppler correction has been performed on the data. For our analysis, the intensity maps are coherently averaged in frequencies for PCA cleaning and therefore cannot be averaged in velocity channels.

In order to calculate the statistical noise for the stacked cube, we perform the same stacking routine on random galaxy positions, similar to the shuffling null test in Section 5.6. The positions are randomly sampled uniformly within the region of the GAMA galaxy sample. The redshifts are randomly sampled following the redshift distribution of the survey. The random process is repeated for 100 realizations to calculate the mean and the standard deviation of the reference cube. The reference cube is used for performing the null test, i.e. checking that the stacked emission signal in the reference cube is consistent with zero.

After obtaining the stacked cubes, we create the angular and spectral stacking by summing over the cubes along the radial and transverse directions. To control the systematics generated by the foreground cleaning, we limit the voxels that go into the summation. For angular stacking, only voxels that are within 150 km s⁻¹ of the central channel are considered, corresponding to the typical width of the HI emission line profile. For spectral stacking, only the voxels whose centre positions are within the 1 deg area of the centre voxel in each channel are included, corresponding to the size of the primary

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Left panel: The angular stacking of the GAMA galaxies on the HI intensity map. The dashed circle denotes the FWHM of the beam. Central panel: The reference null detection image from 100 realizations of random galaxy positions. Right panel: The standard deviation of the reference image multiplied by 3 to visualize the 99.7 per cent confidence interval.</figcaption></figure>

beam. The same summation is performed on the random realization to calculate the mean and the standard deviation of the reference image and spectrum. We note that there is signal loss due to an insufficient number of voxels included. However, as the random realizations follow the same procedure, it does not affect the validation of the detection which is our focus in this work.

## 6.2 The detection of stacked HI signal

In Fig. 19, we show the results for the angular stacking of the GAMA galaxies on the HI intensity map. The stacked map shows a clear detection at the centre of the image, with the centre pixel achieving a detection with a statistical significance of 8.35 σ. The stacked HI emission follows the structure of the beam. The extended structure of emission outside the beam full-width half-maximum (FWHM) is a result of the relatively large thermal noise fluctuation and contamination from nearby sources in the stacking procedure as discussed above. We note that, as shown in Fig. 4, the lower part of the GAMA galaxy sample region sits near the edge of the observed intensity maps. As a result, the lower half of the stacked image has a higher level of statistical fluctuations as shown in the right panel of Fig. 19. Furthermore, the reference image is consistent with a null detection, demonstrating that the stacked signal does not arise from systematic effects in the map.

In Fig. 20, we show the stacked spectrum along the line of sight. The peak of the stacked spectrum shows a detection with a significance of 6.30 σ. The detection is validated by the reference spectrum which is consistent with a null detection. However, as seen in the reference spectrum, there is a systematic noise component with fluctuations on large frequency intervals. The systematics can be most visibly seen at the peaks and troughs of the stacked signal at $v \sim 2500, 4500, \text{ and } 7000$ km s⁻¹. We have identified the source of this fluctuation as the residual systematics in the map after PCA cleaning, which is demonstrated in Appendix D. As stated previously, we are only interested in validating the detection against the null test and leave the suppressing of the systematic effects for future work.

The stacking results show a clear detection of HI signal in the positions of the GAMA galaxies. The detection of an emission line signal is statistically significant against the null test. At the same time, the residual systematics in the data, while being subdominant comparing to the HI signal, induce structures in the stacked spectrum. We intend to follow up on the findings in this work and develop analysis techniques to extract the underlying emission line signal
<figure><img src="image_21.png" /><figcaption>**Figure 20.** The stacked spectrum of the GAMA galaxies on the Hi intensity map. The blue solid line marks the measured stacked spectrum. The orange dashed line denotes the reference spectrum from 100 random realizations, and the shaded region denotes the standard deviation of the realizations. The black dotted line denotes the 5 σ limit of the reference spectrum rms and is smoothed with a Gaussian kernel of 300 km s⁻¹ for better illustration.</figcaption></figure>

in the future, as well as using the stacked signal to understand the residual systematics.

## 7 CONCLUSIONS

In this work we have presented L-band observations from the MeerKLASS intensity mapping campaign. The survey covered a single patch of approximately 200 deg², repeatedly scanned 41 times (27 scans remained after aggressive RFI flagging), which we have referred to as the MeerKLASS L-band deep field. The observations used in our analysis came from only a small section of the full L-band (971.2 < *v* < 1023.6 at redshift 0.39 < *z* < 0.46) to avoid RFI which we find is more present in this band relative to the lower frequency UHF band (580 < *v* < 1000) which all future MeerKLASS surveys will be conducted in. We detailed the calibration strategy used on these observations which included a new approach of self-calibration (Section 2.3), exploiting the data's higher signal-to-noise to devise a less model-dependent calibration. The improvement in data quality relative to previous MeerKLASS pilot surveys is evident (see e.g. Fig. 5) and encouraging for future MeerKLASS surveys gathering ever-larger data sets.

We demonstrated the enhanced signal-to-noise from the calibrated deep-field intensity maps by presenting the autocorrelation power spectrum (Fig. 14) which revealed the lowest upper-limit from a Hi power spectrum at these cosmological scales *k* < 0.3 *h* Mpc⁻¹. This remains a factor of ~2 above where we would expect the Hi power to lie. We will pursue this limit further in follow-up work where we will explore a systematic removal strategy to reduce the quantity of flagged data plus the cross-correlation of divided data sets to mitigate the additive bias from time-dependent residual systematics.

We also presented results from the cross-correlation power spectrum between the MeerKLASS deep-field and 2269 overlapping GAMA galaxies with spectroscopic redshifts. Despite the small field of overlap (~25 per cent of the deep field patch), we still obtained a clear > 4 σ detection of the cross-power spectrum (Fig. 15). Furthermore, we presented a detection of Hi emission from stacking the MeerKLASS maps onto the positions of the GAMA galaxies. This yielded a 8.35 σ and 6.30 σ detection from angular and spectral stacking, respectively, the first stacking detections with MeerKAT in single-dish mode (see Section 6).

We implemented several novel strategies into the power spectrum analysis pipeline, including a selective procedure to target the region of *k*-space with the highest signal-to-noise. The third panel of

Fig. 13 demonstrated the clear window of Fourier modes where power dominates over the error. By implementing some minima and maxima cuts on *k*⊥ and *k*∥ we were able to exclude regions of *k*-space dominated by error. This provides a more robust result since error-dominated modes will not be included in final spherical averaging. Comparison of the error profiles for each *k*-bin (Fig. 17) demonstrated a mild improvement in terms of Gaussianity when using the *k*-cuts, however, the sacrifice of many modes contained at high-*k* could be restricting this improvement. We plan to investigate this strategy in future work where larger data sets will reveal the optimal approach.

We formalized our approach to covariance estimation for the power spectra measurements which utilizes the foreground transfer function distribution from each mock realization. This distribution will be determined by all foreseeable sources of error; residual systematics, thermal noise fluctuations, galaxy shot-noise, cosmic variance, and any uncertainty in signal-loss correction. This should therefore be a robust method for error estimation. However, we found the errors in the spherically averaged cross-power spectrum appear to be large relative to their scatter about the model. This was supported by a χ² analysis. Further inspection revealed that the *k*-bins are mildly correlated and the error profiles for each *k*-bin are not perfectly Gaussian, both of which could be inflating the errors. We hypothesize that this should be naturally mitigated with future MeerKLASS analysis where systematics can be more reliably controlled.

This work also presented validation tests of our improved analysis pipeline. This was achieved by generating a suite of simulations that emulate the real data products. The mock data were then put through the same pipeline as the real data where we could validate its performance. Results in Fig. 8 showed good agreement with the expected model based on the known inputs to generate the mocks, thus validating processes such as Cartesian field regridding and power spectrum estimation. Measured power from the mocks in Fig. 11 which included the foreground cleaning also showed good consistency with real data.

The MeerKLASS campaign is ongoing and will continue to be so in the pre-SKAO era. Further observations have already begun in the telescope's UHF band. These lower frequencies should be less prone to RFI, meaning far fewer channels will be flagged. Furthermore, the lower frequency and wider wave-band mean a much larger portion of redshift space will be covered, providing intensity maps of unprecedented volume.
ACKNOWLEDGEMENTS

The authors would like to acknowledge the ongoing contributions to the MeerKLASS collaboration from members Suman Chatterjee, Karin Fornazier, Tamirat Gogo, Wenkai Hu, Piyanat Kittiwisit, and Sifiso Mahlalela which provided valuable insights for this paper. We would also like to thank Jacob Burba, whose UK SKA Regional Centre project involved a full run-through of our post-calibration analysis pipeline, providing useful validation.

JW acknowledges support from the National SKA Program of China (No. 2020SKA0110100). SCu is supported by a UK Research and Innovation Future Leaders Fellowship grant [MR/V026437/1]. This result is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 948764; PB). JLB acknowledges funding from the Ramón y Cajal Grant RYC2021-033191-I, financed by MCIN/AEI/10.13039/501100011033 and by the European Union 'NextGenerationEU'/PRTR, as well as the project UC-LIME (PID2022-140670NA-I00), financed by MCIN/AEI/10.13039/501100011033/FEDER, UE. SCa acknowledges support from the Italian Ministry of University and Research, PRIN 2022 'EXSKALIBUR—Euclid-Cross-SKA: Likelihood Inference Building for Universe Research', from the Italian Ministry of Foreign Affairs and International Cooperation (grant no. ZA23GR03), and from the European Union-Next Generation EU. IPC is supported by the European Union within the Next Generation EU programme [PNRR-4-2-1.2 project No. SOE.0000136, RadioGaGa]. JF acknowledges support of Fundação para a Ciência e a Tecnologia through the Investigador FCT Contract No. 2020.02633.CEECIND/CP1631/CT0002, the FCT project PTDC/FIS-AST/0054/2021, and the research grants UIDB/04434/2020 and UIDP/04434/2020. AP is a UK Research and Innovation Future Leaders Fellow [grant MR/X005399/1]. ZC's research is supported by a UK Research and Innovation Future Leaders Fellowship [grant MR/X005399/1]. MGS acknowledges support from the South African Radio Astronomy Observatory and National Research Foundation (Grant No. 84156). SCa acknowledges support from the Italian Ministry of University and Research (MUR), PRIN 2022 'EXSKALIBUR—Euclid-Cross-SKA: Likelihood Inference Building for Universe's Research', Grant No. 20222BBYB9, CUP D53D2300252 0006, from the Italian Ministry of Foreign Affairs and International Cooperation (MAECI), Grant No. ZA23GR03, and from the European Union-Next Generation EU. LW is a UK Research and Innovation Future Leaders Fellow [grant MR/V026437/1].

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. We acknowledge the use of the Ilifu cloud computing facility, through the Inter-University Institute for Data Intensive Astronomy (IDIA).

GAMA is a joint European–Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel–ATLAS, GMRT, and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is https://www.gama-survey.org/.

For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

DATA AVAILABILITY

The data underlying this article will be shared upon reasonable request to the corresponding author. Access to the raw data used in the analysis is public (for access information please contact archive@ska.ac.za).

REFERENCES

Agrawal A., Makiya R., Chiang C.-T., Jeong D., Saito S., Komatsu E., 2017, J. Cosmol. Astropart. Phys., 10, 003

Alonso D., Bull P., Ferreira P. G., Santos M. G., 2015, MNRAS, 447, 400

Anderson C. J. et al., 2018, MNRAS, 476, 3382

Andrae R., Schulze-Hartung T., Melchior P., 2010, preprint (arXiv:1012.3754)

Ansari R. et al., 2012, A&A, 540, A129

Asad K. M. B. et al., 2021, MNRAS, 502, 2970

Astropy Collaboration 2022, ApJ, 935, 167

Battye R. A., Davies R. D., Weller J., 2004, MNRAS, 355, 1339

Battye R. A. et al., 2012, preprint (arXiv:1209.1041)

Bera A., Kanekar N., Chengalur J. N., Bagla J. S., 2022, ApJ, 940, L10

Bera A., Kanekar N., Chengalur J. N., Bagla J. S., 2023, ApJ, 950, L18

Bernal J. L., Breysse P. C., Gil-Marín H., Kovetz E. D., 2019, Phys. Rev. D, 100, 123522

Beutler F. et al., 2011, MNRAS, 416, 3017

Bharadwaj S., Nath B., Nath B. B., Sethi S. K., 2001, J. Astrophys. Astron., 22, 21

Blake C., 2019, MNRAS, 489, 153

CHIME Collaboration 2022, ApJS, 261, 29

CHIME Collaboration 2023, ApJ, 947, 16

Carucci I. P., Irfan M. O., Bobin J., 2020, MNRAS, 499, 304

Carucci I. P. et al., 2024, preprint (arXiv:2412.06750)

Chang T.-C., Pen U.-L., Peterson J. B., McDonald P., 2008, Phys. Rev. Lett., 100, 091303

Chen X., 2012, Int. J. Mod. Phys. Conf. Ser., 12, 256

Coles P., Jones B., 1991, MNRAS, 248, 1

Cunnington S., 2022, MNRAS, 512, 2408

Cunnington S., Wolz L., 2024, MNRAS, 528, 5586

Cunnington S., Pourtsidou A., Soares P. S., Blake C., Bacon D., 2020, MNRAS, 496, 415

Cunnington S., Irfan M. O., Carucci I. P., Pourtsidou A., Bobin J., 2021, MNRAS, 504, 208

Cunnington S. et al., 2022, MNRAS, 518, 6262

Cunnington S. et al., 2023, MNRAS, 523, 2453

DES Collaboration, 2022, Phys. Rev. D, 105, 023520

DESI Collaboration 2024, preprint (arXiv:2404.03002)

Dev A. et al., 2023, MNRAS, 523, 2693

Drinkwater M. J. et al., 2010, MNRAS, 401, 1429

Driver S. P. et al., 2009, Astron. Geophys., 50, 5

Driver S. P. et al., 2011, MNRAS, 413, 971

Driver S. P. et al., 2022, MNRAS, 513, 439

Dunne D. A. et al., 2024, ApJ, 965, 7

eBOSS Collaboration, 2021, Phys. Rev. D, 103, 083533

Engelbrecht B. et al., 2024, preprint (arXiv:2404.17908)

Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23

Fixsen D. J., 2009, ApJ, 707, 916

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Foreman S., Obuljen A., Simonović M., 2024, preprint (arXiv:2405.18559)

Guo H., Jones M. G., Haynes M. P., Fu J., 2020, ApJ, 894, 92

Harper S., Dickinson C., 2018, MNRAS, 479, 2024

Haslam C. G. T., Klein U., Salter C. J., Stoffel H., Wilson W. E., Cleary M. N., Cooke D. J., Thomasson P., 1981, A&A, 100, 209

Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

Heymans C., et al., 2021, A&A, 646, A140

Hurley-Walker N. et al., 2022, Publ. Astron. Soc. Aust., 39, e035

Irfan M. O., Bull P., 2021, MNRAS, 508, 3551
Irfan M. O. et al., 2021, MNRAS, 509, 4923

Jarvis M. et al., 2016, in MeerKAT Science: On the Pathway to the SKA. p. 6 preprint (arXiv:1709.01901)

Jing Y. P., 2005, ApJ, 620, 559

Jones M. G., Haynes M. P., Giovanelli R., Moorman C., 2018, MNRAS, 477, 2

Lian X., Xu H., Zhu Z., Hu D., 2020, MNRAS, 496, 1232

Liske J. et al., 2015, MNRAS, 452, 2087

Lujan Niemeyer M. et al., 2022a, ApJ, 929, 90

Lujan Niemeyer M. et al., 2022b, ApJ, 934, L26

Maddox N. et al., 2021, A&A, 646, A35

Masui K. W., et al., 2013, ApJ, 763, L20

Matshawule S. D., Spinelli M., Santos M. G., Ngobese S., 2021, MNRAS, 506, 5075

Meyer S. A., Meyer M., Obreschkow D., Staveley-Smith L., 2016, MNRAS, 455, 3136

Moradinezhad Dizgah A., Nikakhtar F., Keating G. K., Castorina E., 2022, J. Cosmol. Astropart. Phys., 02, 026

Mozdzen T. J., Mahesh N., Monsalve R. A., Rogers A. E. E., Bowman J. D., 2019, MNRAS, 483, 4411

Murray S. G., 2018, J. Open Source Softw., 3, 850

Nan R. et al., 2011, Int. J. Mod. Phys. D, 20, 989

Newburgh L. B. et al., 2016, in Hall H. J., Gilmozzi R., Marshall H. K., eds, Proc. SPIE Conf. Ser. Vol. 9906, Ground-based and Airborne Telescopes VI. SPIE, Bellingham, p. 99065X

Niemeyer M. L., Bernal J. L., Komatsu E., 2023, ApJ, 958, 4

Obuljen A., Simonović M., Schneider A., Feldmann R., 2023, Phys. Rev. D, 108, 083528

Olivari L. C., Remazeilles M., Dickinson C., 2016, MNRAS, 456, 2749

PUMA Collaboration 2019, preprint (arXiv:1907.12559)

Pal S., et al., 2022, MNRAS, 516, 2851

Pan H. et al., 2023, MNRAS, 525, 256

Paul S., Santos M. G., Chen Z., Wolz L., 2023, preprint (arXiv:2301.11943)

Ponomareva A. A. et al., 2021, MNRAS, 508, 1195

Ponomareva A. A. et al., 2023, MNRAS, 522, 5308

Remazeilles M., Dickinson C., Banday A. J., Bigot-Sazy M. A., Ghosh T., 2015, MNRAS, 451, 4311

Rhee J. et al., 2023, MNRAS, 518, 4646

SKA Cosmology SWG 2020, Publ. Astron. Soc. Aust., 37, e007

Santos M. G., Cooray A., Knox L., 2005, ApJ, 625, 575

Santos M. G. et al., 2017, in MeerKAT Science: On the Pathway to the SKA. preprint (arXiv:1709.06099)

Sato-Polito G., Kokron N., Bernal J. L., 2023, MNRAS, 526, 5883

Shaw J. R., Sigurdson K., Sitwell M., Stebbins A., Pen U.-L., 2015, Phys. Rev. D, 91, 083514

Shin J., Kim J., Pichon C., Jeong D., Park C., 2017, ApJ, 843, 73

Sinigaglia F. et al., 2022, ApJ, 935, L13

Sinigaglia F. et al., 2024, MNRAS, 529, 4192

Soares P. S., Cunningham S., Pourtsidou A., Blake C., 2021, MNRAS, 502, 2549

Soares P. S., Watkinson C. A., Cunningham S., Pourtsidou A., 2022, MNRAS, 510, 5872

Spinelli M., Bernardi G., Garsden H., Greenhill L. J., Fialkov A., Dowell J., Price D. C., 2021, MNRAS, 505, 1575

Spinelli M., Carucci I. P., Cunningham S., Harper S. E., Irfan M. O., Fonseca J., Pourtsidou A., Wolz L., 2022, MNRAS, 509, 2048

Switzer E. R. et al., 2013, MNRAS, 434, L46

Switzer E. R., Chang T.-C., Masui K. W., Pen U.-L., Voytek T. C., 2015, ApJ, 815, 51

Tauscher K., Rapetti D., Burns J. O., 2018, J. Cosmol. Astropart. Phys., 12, 015

Thorne B., Dunkley J., Alonso D., Næss S., 2017, MNRAS, 469, 2821

Tramonte D., Ma Y.-Z., 2020, MNRAS, 498, 5916

Tramonte D., Ma Y.-Z., Li Y.-C., Staveley-Smith L., 2019, MNRAS, 489, 385

Tudorache M. N. et al., 2022, MNRAS, 513, 2168

Vanderlinde K. et al., 2019, preprint (arXiv:1911.01777)

Villaescusa-Navarro F., et al., 2018, ApJ, 866, 135

Wang J., et al., 2021, MNRAS, 505, 3698

Westmeier T., Jurek R., Obreschkow D., Koribalski B. S., Staveley-Smith L., 2014, MNRAS, 438, 1176

Wilensky M. J., Brown J., Hazelton B. J., 2023, MNRAS, 521, 5191

Wolz L., Abdalla F. B., Blake C., Shaw J. R., Chapman E., Rawlings S., 2014, MNRAS, 441, 3271

Wolz L., et al., 2017, MNRAS, 464, 4938

Wolz L., et al., 2022, MNRAS, 510, 3495

Wyithe S., Loeb A., Geil P., 2008, MNRAS, 383, 1195

Zwaan M. A., van Dokkum P. G., Verheijen M. A. W., 2001, Science, 293, 1800

## APPENDIX A: OBSERVATION DETAILS

Table A1 provides the details for the MeerKLASS L-band deep field observations, for all 41 attempted observation blocks.
<table><caption>**Table A1.** Information for all 41 observation blocks i.e. each of the scans for the MeerKLASS L-band deep field patch. The last column marks if the block was selected for the final data set. An absent tick means the block was flagged, primarily due to RFI contamination.</caption><thead><tr><th>Block ID<br>(Timestamp)</th><th>Short name<br>(in this paper)</th><th>Observation start time<br>(UTC time)</th><th>Sunset<br>(UTC time)</th><th>az range<br>(°)</th><th>el<br>(°)</th><th>Calibrator</th><th>Motion of field</th><th>Calibrated</th></tr></thead><tbody><tr><td>obs1630519596</td><td>obs210901</td><td>2021-09-01 18:06:53</td><td>09-01 16:14:50</td><td>[98.0, 111.9]</td><td>42.3</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631379874</td><td>obs210911a</td><td>2021-09-11 17:04:51</td><td>09-11 16:20:28</td><td>[100.1, 113.2]</td><td>37.8</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1631387336</td><td>obs210911b</td><td>2021-09-11 19:09:13</td><td>09-11 16:20:28</td><td>[86.6, 108.8]</td><td>63.0</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631552188</td><td>obs210913a</td><td>2021-09-13 16:56:49</td><td>09-13 16:21:35</td><td>[100.1, 113.2]</td><td>37.7</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1631559762</td><td>obs210913b</td><td>2021-09-13 19:02:58</td><td>09-13 16:21:35</td><td>[86.4, 108.8]</td><td>63.3</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631659886</td><td>obs210914</td><td>2021-09-14 22:52:47</td><td>09-14 16:22:08</td><td>[-108.8, -85.8]</td><td>64.1</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1631667564</td><td>obs210915a</td><td>2021-09-15 00:59:40</td><td>09-14 16:22:08</td><td>[-113.1, -99.9]</td><td>38.1</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1631724508</td><td>obs210915b</td><td>2021-09-15 16:48:49</td><td>09-15 16:22:42</td><td>[100.3, 113.4]</td><td>37.2</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1631732038</td><td>obs210915c</td><td>2021-09-15 18:54:17</td><td>09-15 16:22:42</td><td>[86.9, 108.7]</td><td>62.6</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631810671</td><td>obs210916a</td><td>2021-09-16 16:44:51</td><td>09-16 16:23:16</td><td>[100.2, 113.3]</td><td>37.5</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631818149</td><td>obs210916b</td><td>2021-09-16 18:49:27</td><td>09-16 16:23:16</td><td>[86.8, 108.8]</td><td>62.7</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631982988</td><td>obs210918a</td><td>2021-09-18 16:36:47</td><td>09-18 16:24:23</td><td>[100.2, 113.3]</td><td>37.4</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1631990463</td><td>obs210918b</td><td>2021-09-18 18:41:19</td><td>09-18 16:24:23</td><td>[86.9, 108.8]</td><td>62.6</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1632069690</td><td>obs210919a</td><td>2021-09-19 16:41:48</td><td>09-19 16:24:57</td><td>[99.5, 112.9]</td><td>38.9</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1632077222</td><td>obs210919b</td><td>2021-09-19 18:47:18</td><td>09-19 16:24:57</td><td>[85.6, 108.8]</td><td>64.4</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1632184922</td><td>obs210921</td><td>2021-09-21 00:44:34</td><td>09-20 16:25:31</td><td>[-113.7, -100.7]</td><td>36.3</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1632505883</td><td>obs210924</td><td>2021-09-24 17:51:44</td><td>09-24 16:27:48</td><td>[90.4, 109.0]</td><td>57.0</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1632760885</td><td>obs210927</td><td>2021-09-27 16:41:43</td><td>09-27 16:29:33</td><td>[96.6, 111.1]</td><td>45.2</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1633365980</td><td>obs211004</td><td>2021-10-04 16:46:39</td><td>10-04 16:33:46</td><td>[93.1, 109.7]</td><td>52.1</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs16333970780</td><td>obs211011</td><td>2021-10-11 16:46:38</td><td>10-11 16:38:15</td><td>[89.9, 108.9]</td><td>58.0</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1634252028</td><td>obs211014</td><td>2021-10-14 22:55:17</td><td>10-14 16:40:15</td><td>[-113.0, -99.8]</td><td>38.5</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1634402485</td><td>obs211016</td><td>2021-10-16 16:41:45</td><td>10-16 16:41:37</td><td>[88.0, 108.7]</td><td>61.0</td><td>PKS 1934-638</td><td>Rising</td><td>✓</td></tr><tr><td>obs1634748682</td><td>obs211020</td><td>2021-10-20 16:51:50</td><td>10-20 16:44:25</td><td>[84.1, 109.1]</td><td>66.2</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1634835083</td><td>obs211021</td><td>2021-10-21 16:51:43</td><td>10-21 16:45:08</td><td>[83.4, 109.2]</td><td>67.0</td><td>PKS 1934-638</td><td>Rising</td><td>-</td></tr><tr><td>obs1637346562</td><td>obs211119a</td><td>2021-11-19 18:31:59</td><td>11-19 17:08:34</td><td>[-109.0, -84.6]</td><td>65.7</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1637354605</td><td>obs211119b</td><td>2021-11-19 20:44:44</td><td>11-19 17:08:34</td><td>[-113.4, -100.4]</td><td>37.2</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1637691677</td><td>obs211123a</td><td>2021-11-23 18:22:58</td><td>11-23 17:11:59</td><td>[-108.8, -86.4]</td><td>63.2</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1637699408</td><td>obs211123b</td><td>2021-11-23 20:31:33</td><td>11-23 17:11:59</td><td>[-113.8, -101.0]</td><td>35.9</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638130295</td><td>obs211128</td><td>2021-11-28 20:13:37</td><td>11-28 17:16:12</td><td>[-113.9, -101.1]</td><td>35.5</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638294319</td><td>obs211130a</td><td>2021-11-30 17:47:13</td><td>11-30 17:17:51</td><td>[-108.9, -84.9]</td><td>65.3</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638301944</td><td>obs211130b</td><td>2021-11-30 19:53:44</td><td>11-30 17:17:51</td><td>[-113.1, -100.0]</td><td>38.1</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638386189</td><td>obs211201</td><td>2021-12-01 19:18:31</td><td>12-01 17:18:39</td><td>[-111.2, -96.8]</td><td>44.7</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638639082</td><td>obs211204a</td><td>2021-12-04 17:33:15</td><td>12-04 17:21:03</td><td>[-108.9, -85.4]</td><td>64.6</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1638647186</td><td>obs211204b</td><td>2021-12-04 19:47:48</td><td>12-04 17:21:03</td><td>[-113.9, -101.1]</td><td>35.7</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1638898468</td><td>obs211207</td><td>2021-12-07 17:36:20</td><td>12-07 17:23:20</td><td>[-108.7, -87.7]</td><td>61.4</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1639157507</td><td>obs211210</td><td>2021-12-10 17:33:41</td><td>12-10 17:25:31</td><td>[-108.8, -88.8]</td><td>59.7</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr><tr><td>obs1639331184</td><td>obs211212</td><td>2021-12-12 17:48:20</td><td>12-12 17:26:53</td><td>[-109.3, -91.8]</td><td>54.6</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1639935088</td><td>obs211219</td><td>2021-12-19 17:34:00</td><td>12-19 17:31:08</td><td>[-109.5, -92.8]</td><td>52.7</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1640540184</td><td>obs211226</td><td>2021-12-26 17:39:00</td><td>12-26 17:34:21</td><td>[-111.0, -96.3]</td><td>45.8</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1640712986</td><td>obs211228</td><td>2021-12-28 17:38:01</td><td>12-28 17:35:03</td><td>[-111.6, -97.5]</td><td>43.2</td><td>Pictor A</td><td>Setting</td><td>-</td></tr><tr><td>obs1640799689</td><td>obs211229</td><td>2021-12-29 17:43:29</td><td>12-29 17:35:22</td><td>[-112.2, -98.5]</td><td>41.2</td><td>Pictor A</td><td>Setting</td><td>✓</td></tr></tbody></table>

## APPENDIX B: PCA FOREGROUND CLEANING

### B1 Impact of reconvolution

As discussed in the main text, the reconvolution of the maps before foreground cleaning (outlined in Section 4.1), slightly improves the foreground clean performance. This is not necessarily just because it homogenizes the effective resolution of the beam, as is the primary aim of this technique. The smoothing will also suppress systematics in the data. Residual contributions, for example from RFI, will distort the true amplitude of the sky map. If these contributions vary with frequency, which will be the case for most systematics, then this will create small perturbations to the spectra along each pixel. Despite smoothing the maps in the *angular* direction only, this will still suppress perturbations along the line of sight, hence providing

slightly smoother spectra, which may be linked to the improved foreground clean.

We demonstrate this in Fig. B1 which plots the amplitudes of each pixel along frequency, as previously shown in Fig. 5 (here we also trim the edges of the map, as outlined in Section 4.2). The top panel shows the spectra without reconvolution. A clear improvement in the smoothness of the spectra is evident when reconvolution is used (bottom panel).

Future work is planned where a more realistic MeerKAT beam will be implemented into our modelling. In this work, we will also investigate this technique of reconvolution and aim to conclude if the smoothing to a common resolution is the source of any improvement or if it is just entirely from the suppression of systematics, reducing the perturbations they cause in the spectra which aids the PCA process.
<figure><img src="image_22.png" /><figcaption>**Figure B1.** The impact of reconvolution on the smoothness of the sky map spectra. Each plotted line is the amplitude along frequency, which, for these pre-foreground clean sky maps, should trace the dominant synchrotron. The top panel is with no reconvolution. The bottom panel shows the smoothing effect from the reconvolution.</figcaption></figure>

## B2 Weighted frequency covariance estimation

Here we explore the small impact that opting to use weights in the covariance estimation has on the PCA cleaning process. The weights used are constant along the frequency direction, so not to add structure to the modes and diminish the PCA performance. Furthermore, after trimming the map edges, the time stamps have a more uniform coverage (see Fig. 9). This means the weights do not vary considerably across our whole field, so we would expect only minor changes by including them.

Fig. B2 shows the first nine eigenmodes from U, extracted from the diagonalization of the covariance matrix computed using equation (29). We show the difference between the weighted case (blue solid) and the unweighted (orange dashed), equivalent to omitting the w factor in the covariance calculation. As expected, the differences are marginal, especially when considering that in some modes the discrepancy is mainly an arbitrary change of sign. However, some small differences are still noticeable, particularly for the higher eigenmodes where noise will become more of a component which is effectively what we are weighting for.

We found little evidence of improved performance in the final results whether weights were used or not in the foreground clean. The difference between the cleaned maps was at the sub-1µK level, and testing the process on the mocks showed no clear performance gain. However, we included them in our final results for completeness. Furthermore, the weighted formalism may be useful in future work focusing on H1 autocorrelations, where cross-correlating subsets of maps from different dishes or time blocks will be pursued. Here, different map weights can be used in a combined cross-covariance computation; then asymmetric singular vectors can be extracted from

<figure><img src="image_23.png" /><figcaption>**Figure B2.** The first nine eigenvectors (descending order of dominance from top to bottom) extracted from the frequency covariance of the MeerKLASS deep field intensity maps. Eigenmode 1 shows the overwhelming synchrotron signal, a smooth slope with a spectral index power law. The different lines show the contrast between an unweighted and weighted PCA clean.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure B3.** The first seven weighted eigenmodes from Fig. B2, projected onto the map space, thus showing the map content removed in each PCA mode in the foreground cleaning. We do not show Eigenmode 1 as this is dominated by the synchrotron and is indistinguishable from the sky map of Fig. 4. Each map is averaged along frequency.</figcaption></figure>

this, as in the singular value decomposition (SVD) approach of Switzer et al. (2013).

To accompany the eigenmode plot of Fig. B2, we also show the projected map from the first seven eigenmodes in Fig. B3 for the weighted PCA case. We omit the first eigenmode since it is indistinguishable from the synchrotron sky maps we have shown

previously. These show at map level, each component that is being removed in the PCA clean. For example, a $N_{fg}=7$ clean would remove all these components. No obvious systematic structure is apparent which could otherwise be targeted for removal by other means, rather than relying on an aggressive PCA clean. However, as we discussed in Section 2.4.2, when cleans are performed individually on smaller subsets of time blocks, constant declination stripes become apparent in some of the components which we will aim to correct for in calibration in future work.

## B3 Power index of the transfer function for signal reconstruction in auto- and cross-correlations

The transfer function defined in equation (33) is applied to the measured power spectrum to reconstruct the loss of information due to foreground cleaning. A priori, intuition would lead to applying the transfer function as defined for the cross-power spectrum of HI and galaxies, and the square of the transfer function for the HI autopower spectrum, as signal loss occurs twice in that calculation (Switzer et al. 2013). However, Cunnington et al. (2023) found empirically that squaring the transfer function overcorrected the signal loss and that the transfer function should be applied equally to cross- and autopower spectra for correct modelling. In this appendix, we provide simplified analytic arguments to support such a claim.

Let us assume the idealized scenario where we can decompose the cosmological fluctuations $\delta T_{HI}(\mathbf{k}) \equiv r(\mathbf{k}) + l(\mathbf{k})$, where $r$ and $l$ correspond to the residual and loss modes of the cosmological signal, respectively. Under this separation, the transfer function can be understood as

$$ \mathcal{T}_{HI}(\mathbf{k}) = \frac{\langle r_m(\mathbf{k}) [r_m(\mathbf{k}') + l_m(\mathbf{k}')]^* \rangle}{\langle [r_m(\mathbf{k}) + l_m(\mathbf{k})] [r_m(\mathbf{k}') + l_m(\mathbf{k}')]^* \rangle}, \quad (B1) $$

where the subscript ‘m’ refers to quantities obtained from mock injection. Due to the component separation process in PCA, $r$ and $l$ are a combination of independent components, and are therefore uncorrelated. The power spectrum of their sum is

$$ \begin{aligned} \langle [r_m(\mathbf{k}) + l_m(\mathbf{k})] [r_m(\mathbf{k}') + l_m(\mathbf{k}')]^* \rangle &= \frac{1}{(2\pi)^3} \delta_D(\mathbf{k} - \mathbf{k}') [P_m^{(r)}(k) + P_m^{(l)}(k)] \\ &= \frac{1}{(2\pi)^3} \delta_D(\mathbf{k} - \mathbf{k}') P_m^{HI}(k). \end{aligned} \quad (B2) $$

Note that the equation above applies equally to the actual observed cosmological signal. Therefore, the transfer function can be now understood as

$$ \mathcal{T}_{HI}(\mathbf{k}) \sim \frac{P_m^{(r)}(k)}{P_m^{(r)}(k) + P_m^{(l)}(k)}. \quad (B3) $$

Let us now consider the application of this conceptual transfer function to the power spectra measured from the clean maps. In this case, we have

$$ \frac{\langle r(\mathbf{k}) r^*(\mathbf{k}') \rangle}{\mathcal{T}_{HI}(\mathbf{k})} \sim P^{(r)}(k) \left[ \frac{P_m^{(r)}(k) + P_m^{(l)}(k)}{P_m^{(r)}(k)} \right] \quad (B4) $$

for the autopower spectrum, and

$$ \frac{\langle r(\mathbf{k}) g^*(\mathbf{k}') \rangle}{\mathcal{T}_{HI}(\mathbf{k})} \sim P^{(rg)}(k) \left[ \frac{P_m^{(r)}(k) + P_m^{(l)}(k)}{P_m^{(r)}(k)} \right] \quad (B5) $$

for the cross-power spectrum. Assuming $P^{(r)} \sim P_m^{(r)}$, it is clear that a single power of the transfer function (i.e. $\mathcal{T}^{1}$) correctly reconstructs the autopower to the desired $P^{(r)} + P^{(l)}$ level. For the cross-power
spectrum, let us simplify the situation and assume linear theory, for which $\delta_g \propto \delta T_{HI} \sim r+l$. Then, $P^{(rg)} \propto P^{(r)}$; if we denote the proportionality factor generically as $\mathcal{A}$, we have

$$ \frac{\langle r(\mathbf{k})g^*(\mathbf{k}') \rangle}{\mathcal{T}_{\mathrm{HI}}(\mathbf{k})} \sim \mathcal{A} P^{(r)}(k) \left[ 1 + \frac{P_m^{(l)}(k)}{P_m^{(r)}(k)} \right], \quad (B6) $$

which similarly shows that a single power of the transfer function must also be applied to the cross-power spectrum.

## APPENDIX C: PARAMETER INFERENCE FROM THE CROSS-POWER

We have not dedicated a lot of attention to the modelling of the power spectra in this work. This is mostly due to the large errors in the cross-power due to the small sample overlap and the small additive bias still present in the auto-HI power spectrum. This makes high dimensional models unnecessary and we opted instead for a simple amplitude fit of a fiducial model with some assumptions on parameters (see Table 1). Some of these parameters are poorly constrained however, so here we explore what happens if we extend the model and include an MCMC fit of four parameters; $\varphi = \{\Omega_{HI}, b_{HI}, b_g, \sigma_v\}$. $\Omega_{HI}$ is the dimension-less HI density parameter and can be directly linked to the mean HI temperature through

$$ \bar{T}_{\mathrm{HI}}(z) = 180 \Omega_{\mathrm{HI}}(z) h \frac{(1+z)^2}{\sqrt{\Omega_m(1+z)^3 + \Omega_{\Lambda}}} \mathrm{mK}, \quad (C1) $$

where $\Omega_m$ and $\Omega_\Lambda$ are the density parameters for matter and the cosmological constant, respectively.

We assume wide and flat priors for the $\varphi$ parameters given by

(i) $0 < 10^3 \Omega_{\mathrm{HI}} < 2$

(ii) $0.5 < b_{\mathrm{HI}} < 2.0$

(iii) $0 < \sigma_v < 600 \, [\mathrm{km s}^{-1}]$

(iv) $1.8 < b_g < 2.0$

A tighter prior is assumed for the GAMA galaxy bias since we fit the galaxy autocorrelation power spectrum, finding a best-fitting bias value of $b_g \sim 1.9$. The remaining priors are sufficiently wide, even for these poorly constrained parameters. We stress that this is still an overly simplistic approach, but our aim here is a demonstration of current modelling capabilities, not rigorous parameter inference.

Fig. C1 shows the inferred parameter constraints from the MCMC fit. We use $\mathrm{emc} \equiv \{\mathrm{Foreman-Mackey et al. 2013}\}$ with 500 iterations of 200 walkers. We checked the chains had converged and discarded a burn-in of 40 iterations. Since it is entirely prior dominated in the cross-power, we marginalize over $b_g$ and do not show its parameter contour. The nuisance parameter $\sigma_v$ is also completely prior dominated and suggests that a more refined model is needed on smaller scales, something explored in CHIME Collaboration (2023). There is a constraint on $\Omega_{HI}$ but this is only due to the restricted prior on $b_{HI}$. Relaxing the HI bias prior would inflate the $\Omega_{HI}$ constraint indefinitely.

Fig. C1 thus demonstrates the degeneracy between $\Omega_{HI}$ and $b_{HI}$ which has led previous efforts to merge this parameter fit in different ways (Anderson et al. 2018; Cunnington et al. 2022; Wolz et al. 2022; CHIME Collaboration 2023). This degeneracy can be broken by utilizing RSD which introduces an anisotropic dependence on $b_{HI}$ but not on $\Omega_{HI}$. We tested a measurement of the power spectrum quadrupole (a smoking-gun of RSD) (Cunnington et al. 2020; Soares et al. 2021) and found a clear non-zero measurement. However, we found that this was overwhelmingly driven by observational effects and higher signal-to-noise is required to make a distinction between this and a true RSD-induced anisotropy. Further discussion

<figure><img src="image_25.png" /><figcaption>**Figure C1.** Parameter constraints from MCMC for $\varphi = \{\Omega_{HI}, b_{HI}, b_g, \sigma_v\}$, showing the degeneracies that exist from a monopole-only measurement. Flat priors are given in the text. We do not show $b_g$ since this is entirely prior dominated in the cross-power.</figcaption></figure>

on the degeneracies between astrophysics and cosmology and the capabilities of the LIM power spectrum multipoles to constrain them can be found in Bernal et al. (2019), where optimal parametrizations are proposed.

## APPENDIX D: EFFECTS OF PCA IN THE STACKED SIGNAL

In this section, we use mock HI galaxy signals to generate signal-only maps for stacking, and illustrate the effects of PCA in creating additional systematic noise structure.

As we are only interested in simulating the effects for stacking in this section, we use generate Poisson realizations of HI sources with no clustering information as shown in Fig. D1. For galaxies inside the survey coverage of GAMA, we use the exact same positions as the GAMA galaxies to replicate the double-counting effects. The rest

<figure><img src="image_26.png" /><figcaption>**Figure D1.** The positions of the mock galaxies on the sky. The blue points indicate the galaxies outside the GAMA coverage, and are not used for the mock stacking. The red points follow the same positions as the GAMA galaxies with random HI mass assigned.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure D2.** Top panel: The stacked spectra of the signal-only mock cube with and without PCA cleaning against the input mock signal. Bottom panel: The stacked spectrum of the signal-only mock cube with PCA cleaning, compared against the reference spectrum stacked on random positions. The dashed line shows the mean of the reference spectrum from 100 realizations and the shaded region shows the standard deviation among the realizations.</figcaption></figure>

of the map is then filled with galaxies at random positions assuming the same number density as the GAMA field. Random redshifts are assigned to the galaxies following the redshift kernel shown in Fig. 7. Each galaxy is then assigned a H1 mass following the H1 mass function reported in Jones et al. (2018). Furthermore, we assign the velocity FWHM of the emission line profile $w_{50}$ to each galaxy using the Tully-Fisher relation reported in Ponomareva et al. (2021). The emission line profile of each galaxy is then calculated and added to the mock signal cube following the simplified busy function (Westmeier et al. 2014)

$$B_2(v) = \frac{a}{2} \times (\text{erf}[b(w_{50}^2 - v^2)] + 1) \times (cv^2 + 1), \quad (\text{D1})$$

where erf is the error function and a, b, c are parameters to be set. For each galaxy, we randomly choose b logarithmically uniformly in $[10^{-3}, 10^{-2}]$ km$^{-2}$s$^2$ and c in $[10^{-2}, 10^0]$ km$^{-2}$s$^2$. The amplitude parameter a is then set by the H1 mass of the galaxy. We find that varying the choice of the busy function parameters does not significantly impact the results due to the relatively low velocity resolution. The signal is then convolved with the resmoothed beam described in Section 4.1.

As discussed in Section 6, the low angular resolution of the map induces leakage from nearby sources in the stacked profile. To illustrate this, we perform the stacking routine on the signal-only mock cube and compare against the input mock signal in Fig. D2. The input signal is calculated by averaging the emission line profile of each source according to equation (D1). The double-counting introduces a large plateau of contamination across the frequencies. We note that the increasing slope of the plateau is due to the redshift distribution of sources. We then apply the PCA cleaning to the mock signal cube and repeat the stacking procedure. Instead of calculating the mixing matrix from the mock cube, we use the exact same mixing matrix obtained from the map for PCA and show the results in Fig. D2 (orange dashed line). The PCA removes the plateau to recover an emission line profile. However, it introduces several effects. There is signal loss due to the cleaning as well as the fact that we did not sum over an area much larger than the beam. Moreover, the PCA extends the width of the profile as seen in the lower panel of Fig. D2.

<figure><img src="image_28.png" /><figcaption>**Figure D3.** The stacked spectra of the signal-only mock cube with and without simulated systematics. The reference spectrum and its standard deviation, which are consistent with a null detection, is also shown for reference.</figcaption></figure>

Apart from the effects on the recovered signal, PCA also creates a systematic noise component as seen in the stacked spectrum. To further demonstrate this, we repeat the reference spectrum calculation described in Section 6 on the mock signal and present the results in Fig. D2. Compared to the reference spectrum, the stacked mock signal outside the centre area of the peak shows deviations. The deviation from the reference spectrum appears to be coherent over large velocity intervals $\Delta v \sim 5000$ km s$^{-1}$. The scale of the interval is consistent with the systematic effects seen in the stacked spectrum in Fig. 20. We emphasize that this is purely the combination of the H1 signal and the PCA cleaning since the mock map does not contain any other component. We further find that, when increasing the angular area for summation from $(1 \text{ deg})^2$ to $(3 \text{ deg})^2$, the amplitude of the coherent structure increases.

The structure of the PCA cleaning is caused by the residual systematics in the data. Following the structure of the stacked spectrum in Fig. 20, we generate random Gaussian realizations of systematic noise components and convolve the systematics to the mock signal. We then repeat the stacking procedure as well as reference spectrum calculation and show the results in Fig. D3. As one can see, the mock signal affected by the systematics shows similar oscillating structure as seen in the data. The reference spectrum, on the other hand, shows consistency with a null detection, and demonstrates the fact that the detection reported in Section 6 is a result of the detection H1 signal with residual systematics, and does not arise purely from contamination in the data.

We intend to simulate the mock data cube with foregrounds, noise, and systematics to fully understand the impact of residual systematics on the stacked spectrum in the follow-up work.

<footnote n="1">¹Dipartimento di Fisica, Università degli Studi di Torino, via P. Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="2">²Dipartimento di Fisica, Università degli Studi di Milano, via G. Celoria 16, I-20133 Milano, Italy</footnote>

<footnote n="3">³INFN-Istituto Nazionale di Fisica Nucleare, Sezione di Milano, via G. Celoria 16, I-20133 Milano, Italy</footnote>

<footnote n="4">⁴INAF-Istituto Nazionale di Astrofisica, Osservatorio Astrofisico di Brera-Merate, via Brera 28 I-20121 Milano, Italy</footnote>

<footnote n="5">⁵Instituto de Física de Cantabria (IFCA), CSIC-Univ. de Cantabria, Avda. de los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="6">⁶Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="7">⁷Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, Cape Town 7535, South Africa</footnote>

<footnote n="8">⁸INFN-Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="9">⁹INAF-Istituto Nazionale di Astrofisica, Osservatorio Astrofisico di Torino, I-10025 Pino Torinese, Italy</footnote>

<footnote n="10">¹⁰INAF-Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="11">¹¹IFPU-Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34151 Trieste, Italy</footnote>
This paper has been typeset from a *TEX/LATEX* file prepared by the author.

<footnote n="12">¹²Institute for Astronomy, The University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="13">¹³Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto CAUP, P-4150-762 Porto, Portugal</footnote>

<footnote n="14">¹⁴Departamento de Física e Astronomia, Faculdade de Ciências, Universi-
dade do Porto, Rua do Campo Alegre 687, P-4169-007 Porto, Portugal</footnote>

<footnote n="15">¹⁵Institute of Astrophysics, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="16">¹⁶Department of Physics, College of Sciences, Northeastern University, Wenhua Road, Shenyang 11089, China</footnote>

<footnote n="17">¹⁷Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK</footnote>

<footnote n="18">¹⁸South African Radio Astronomy Observatory (SARAO), 2 Fir Street, Cape Town 7925, South Africa</footnote>

<footnote n="19">¹⁹Observatoire de la Côte d'Azur, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, F-06304 Nice cedex 4, France</footnote>

<footnote n="20">²⁰Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China</footnote>