# The foreground transfer function for Hi intensity mapping signal
reconstruction: MeerKLASS and precision cosmology applications

Steven Cunnington1, Laura Wolz1, Philip Bull1,2, Isabella P. Carucci3,4, Keith
Grainge1, Melis O. Irfan2,5, Yichao Li6,2, Alkistis Pourtsidou7,8,2, Mario G.
Santos2,9, Marta Spinelli10,2, Jingying Wang11,2
1Jodrell Bank Centre for Astrophysics, Department of Physics & Astronomy, The
University of Manchester, Manchester M13 9PL, UK
2Department of Physics and Astronomy, University of the Western Cape, Robert
Sobukwe Road, Cape Town, 7535, South Africa
3Dipartimento di Fisica, Università degli Studi di Torino, via P. Giuria 1,
10125, Torino, Italy
4INFN – Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P.
Giuria 1, 10125, Torino, Italy
5Department of Physics and Astronomy, Queen Mary University of London, London,
E1 4NS, UK
6Department of Physics, College of Sciences, Northeastern University, Wenhua
Road, Shenyang, 11089, China
7Institute for Astronomy, The University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, UK
8Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK
9South African Radio Astronomy Observatory (SARAO), 2 Fir Street, Cape Town,
7925, South Africa
10Institute of Particle Physics & Astrophysics, Department of Physics, ETH
Zurich, Switzerland
11Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai, 200030, China
E-mail<EMAIL_ADDRESS>behalf of the MeerKLASS
Collaboration

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Blind cleaning methods are currently the preferred strategy for handling
foreground contamination in single-dish Hi intensity mapping surveys. Despite
the increasing sophistication of blind techniques, some signal loss will be
inevitable across all scales. Constructing a corrective transfer function
using mock signal injection into the contaminated data has been a practice
relied on for Hi intensity mapping experiments. However, assessing whether
this approach is viable for future intensity mapping surveys where precision
cosmology is the aim, remains unexplored. In this work, using simulations, we
validate for the first time the use of a foreground transfer function to
reconstruct power spectra of foreground-cleaned low-redshift intensity maps
and look to expose any limitations. We reveal that even when aggressive
foreground cleaning is required, which causes ${>}\,50\%$ negative bias on the
largest scales, the power spectrum can be reconstructed using a transfer
function to within sub-percent accuracy. We specifically outline the recipe
for constructing an unbiased transfer function, highlighting the pitfalls if
one deviates from this recipe, and also correctly identify how a transfer
function should be applied in an auto-correlation power spectrum. We validate
a method that utilises the transfer function variance for error estimation in
foreground-cleaned power spectra. Finally, we demonstrate how incorrect
fiducial parameter assumptions (up to ${\pm}100\%$ bias) in the generation of
mocks, used in the construction of the transfer function, do not significantly
bias signal reconstruction or parameter inference (inducing ${<}\,5\%$ bias in
recovered values).

###### keywords:

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

††pubyear: 2023††pagerange: The foreground transfer function for Hi intensity
mapping signal reconstruction: MeerKLASS and precision cosmology
applications–D

##  1 Introduction

Probing fluctuations in the Universe’s density field is an excellent tool for
furthering precision cosmology. A number of large sky surveys have been
commissioned with this aim and have contributed towards constraining
parameters in the standard cosmological model (eBOSS Collaboration, 2021;
Heymans et al., 2021; DES Collaboration, 2022). Despite cosmic microwave
background (CMB) experiments leading the way with constraints (Planck
Collaboration, 2020), it is expected that large-scale structure maps will soon
be the leading resource given the three-dimensional information they provide
(Slosar et al., 2008; Giannantonio et al., 2012; Alonso et al., 2015b).
Increases in survey volume will allow fluctuations across the largest scales
to be probed, improving constraints. It is within the relatively unexplored
ultra-large scales where novel tests of general relativity will be possible
and where there will be increased sensitivity to new physics such as non-
Gaussian fluctuations in the Universe’s primordial density field (Camera et
al., 2013; Camera et al., 2015; Alvarez et al., 2014; Fonseca et al., 2015;
Baker & Bull, 2015; Bull, 2016; Cunnington, 2022). Ultra-large scales are also
highly linear, avoiding the complex modelling challenges facing surveys
attempting to exploit non-linear regimes (D’Amico et al., 2020; Martinelli et
al., 2021; Pourtsidou, 2023).

An efficient method for surveying large volumes is using radio telescopes to
map the redshifted 21cm emission from neutral hydrogen (Hi). The Hi, which
mostly resides in galaxies in the post-reionization Universe, traces the
underlying dark matter, thus allowing the Universe’s large-scale structure to
be probed. By rapidly scanning the sky and recording all radiation as
unresolved diffuse emission, the faint 21cm signals are integrated allowing
for a comprehensive survey of Hi density in 3D to be obtained. This process is
known as Hi intensity mapping (Bharadwaj et al., 2001; Battye et al., 2004;
Wyithe et al., 2008; Chang et al., 2008).

The Hi power spectrum has been recently detected on small Mpc scales (Paul et
al., 2023) with intensity maps from the 64-dish MeerKAT array, a pathfinder
telescope for the Square Kilometre Array Observatory (SKAO) (SKA Cosmology
SWG, 2020). This detection used MeerKAT as an interferometer which has higher
sensitivity on small scales. However, within the next few years the MeerKAT
Large Area Synoptic Survey (MeerKLASS) plans to conduct a wide
(${\gtrsim}\,4{,}000\,\text{deg}^{2}$) Hi intensity mapping survey,
potentially spanning $0.4\,{<}\,z\,{<}\,1.45$ in redshift if performed using
the UHF band (Santos et al., 2017). Since the MeerKAT interferometer does not
have sufficiently short baselines to achieve such a field of view, the
observations will be gathered using the single-dish data from each element of
the array. This auto-correlation mode of observation, often referred to as
single-dish mode, is also planned for the full SKAO in order to probe large-
scale cosmology, which has been identified as a top priority science objective
(Weltman et al., 2020). In the pre-SKA era however, MeerKAT will pursue low-
redshift Hi intensity mapping and has already demonstrated calibration and
map-making from single-dish mode observations with a small pilot survey (Wang
et al., 2021). This same pilot survey was also used to achieve the first
single-dish mode cosmological detection with a multi-dish array in cross-
correlation with an overlapping galaxy survey (Cunnington et al., 2022).

Since Hi intensity mapping records all emission in the frequency range of the
instrument, the major challenge is removing any signals which are not
cosmological Hi. This can include radio frequency interference (RFI) and
astrophysical foregrounds, both of which can dominate by orders of magnitude
over the weak Hi signal111Additional contaminants come from atmospheric
emission and ground pickup which can be approximately modelled as constant
over time when a constant elevation scanning strategy is adopted.. In
principle, RFI should be time-varying and can be flagged when it enters the
observations. However, foregrounds will consistently enter the observations
due to their fixed sky coordinates, therefore a process for separating them
from the Hi is required. The dominant sources producing foregrounds in the low
redshift Hi frequency range (${\sim}\,300\,{<}\,\nu\,{<}\,1420$ MHz) are
synchrotron and free-free radiation from within our own galaxy, along with
extra-galactic concentrated emission from strong point sources such as active
galactic nuclei (Oh & Mack, 2003; Santos et al., 2005; Alonso et al., 2014).

To date, blind foreground cleaning techniques have been the only approach that
has led to cross-correlation detections of a cosmological power spectrum in
single-dish intensity mapping (Masui et al., 2013; Wolz et al., 2017; Anderson
et al., 2018; Wolz et al., 2022; Cunnington et al., 2022). Blind techniques
exploit the robust assumption that foregrounds are a dominant and correlated
contribution to the observations and can be statistically reduced into a few
components which are removed. This requires little prior knowledge of the
foregrounds which is a huge advantage since it is challenging to obtain a
detailed understanding of the foreground’s precise amplitude through
frequency, or how they respond to instrumental systematics. Blind foreground
removal performed at the map level has proven to be the most successful
approach and these have been validated and refined in simulations (Wolz et
al., 2014; Alonso et al., 2015a; Carucci et al., 2020; Cunnington et al.,
2021a; Spinelli et al., 2021). Interferometric intensity mapping can to some
extent adopt a foreground avoidance strategy, which assumes they are isolated
in a foreground wedge region in ($k_{\perp},k_{\parallel}$)-Fourier space (Liu
et al., 2014; Paul et al., 2023). The foreground avoidance technique has the
advantage of being immune to signal loss from foreground cleaning. However,
recent studies have shown some component separation improves foreground
mitigation relative to foreground avoidance alone in interferometric surveys
(Chen et al., 2022). Hence blind foreground removal is likely to be an adopted
technique beyond single-dish mode experiments.

Whilst blind foreground cleaning algorithms themselves have been well studied,
the precise effects they cause on the underlying Hi field have not been to the
same extent. Broadly speaking there are two unfavourable consequences from
blind foreground cleaning, and both will occur simultaneously to some extent.
The first is foreground residuals, i.e. foreground contamination not removed
from the data resulting from under-cleaning. The second is signal loss i.e.
the reduction in the Hi power spectrum amplitude resulting from the foreground
clean. It is this second issue that is the main focus of this paper. Whilst
foreground residuals are of course important, their influence on the data is
similar to RFI and instrumental noise, that is they cause an additive bias to
the estimated Hi power spectrum. However, for foregrounds, it is expected that
their residuals should be reducible to sub-dominant levels relative to the Hi
(Cunnington et al., 2021a). Furthermore, in cross-correlation with a
foreground-free tracer such as a galaxy survey, any additive bias from
foreground residuals will be absent and the only impact will be on the error
budget.

For signal loss, it has been shown that even in ideal simulations, some loss
is always inevitable across all scales when blind foreground removal methods
are applied, and this is not mitigated in cross-correlations (Cunnington et
al., 2021a). Ignoring or incorrectly estimating signal loss, unsurprisingly,
leads to a biased recovery of the Hi power spectrum. Thus signal loss is a
crucial concept to understand exhaustively for precision cosmology to be
possible with Hi intensity mapping. The necessary process of signal
reconstruction i.e. correcting for the signal loss, is where there is little
dedicated study. A foreground transfer function $\mathcal{T}$ can be simply
defined as the object which delivers a reconstructed signal power spectrum
that is unbiased to the underlying truth i.e. $\langle
P_{\text{rec}}(k)\rangle\,{=}\,P_{\text{true}}(k)$ where $\langle
P_{\text{rec}}(k)\rangle\,{\equiv}\,P_{\text{clean}}(k)/\mathcal{T}(k)$.
Previous intensity mapping detections (Masui et al., 2013; Anderson et al.,
2018; Wolz et al., 2022; Cunnington et al., 2022) have all relied on a process
of mock signal injection to estimate the foreground transfer function. By
subjecting the injected mocks to the same foreground cleaning process as the
observations, we can use the drop in the measured mock power spectra to
estimate the transfer function. This method was first extensively analysed in
Switzer et al. (2015) in the context of the Green Bank Telescope (GBT) Hi
intensity maps (Masui et al., 2013; Switzer et al., 2013). There have been
similar methods of signal injection to correct for signal loss implemented for
epoch of reionisation experiments where past analyses underestimated signal
loss leading to biased results (Cheng et al., 2018), highlighting the
importance of correctly understanding this issue. To date, there has been no
dedicated simulations-based investigation into the reliability of the transfer
function for low-redshift Hi intensity maps blindly cleaned at the map level,
despite its clear importance.

In this work, we use various Hi intensity mapping simulations to validate the
reliability of the transfer function for signal reconstruction in foreground-
cleaned maps. We explore how signal loss arises in a foreground clean and
illustrate the subtleties of this both analytically and empirically in
simulation tests, demonstrating how a transfer function can be correctly
estimated to account for all these subtleties. We focus exclusively on
Principal Component Analysis (PCA)-based foreground cleaning but much of the
formalism and results presented will be transferable to other foreground
cleaning techniques. Furthermore, whilst our focus is on single-dish intensity
mapping, the conclusions will also be applicable to interferometers. We
demonstrate how a foreground transfer function is a reliable tool for small
pilot surveys, validating it on simulations constructed using empirical
MeerKAT 2019 observations aiming to realistically emulate current MeerKAT
intensity maps. Lastly, we look to the future and pursue to what extent the
transfer function can be relied on for conducting precision cosmology with
intensity mapping where sub-percent accuracy on parameter estimation is the
aim.

This paper is structured as follows; in Section 2 we present an overview of
the formalism for a blind PCA-based foreground clean, explicitly highlighting
where signal loss arises. Section 3 presents how one should construct a
foreground transfer function to correct for signal loss. In Section 4 we test
the transfer function on a simulation of a MeerKAT-like intensity mapping
pilot survey, validating the transfer function in this low signal-to-noise
regime. Section 5 focuses on how suited the transfer function is for the
purposes of precision cosmology, showcasing the robustness of the transfer
function even where the fiducial cosmology assumed for its construction
disagrees with the truth in the observations. Finally, we conclude in Section
6.

##  2 Signal loss from foreground cleaning

We begin with a pedagogical introduction to the formalism describing blind
foreground cleaning and with the aid of simulations demonstrate some key
concepts of signal loss induced by the foreground clean. For consistency, we
largely follow the notation in Switzer et al. (2015). Whilst we focus on a
PCA-based method, the formalism we present is in principle transferable to
more sophisticated blind foreground removal techniques when they are used as
linear filters (e.g. Bobin et al., 2007; Chapman et al., 2012; Alonso et al.,
2015a; Carucci et al., 2020; Cunnington et al., 2021a; Irfan & Bull, 2021;
Spinelli et al., 2021). We blackuse a set of simulations with separable Hi
signal and foreground contributions, allowing us to provide examples of the
claims made in certain derivations. We begin by using some generic simulations
which are similar to that used in Cunnington et al. (2021a). The exact details
of these simulations are outlined in Appendix A.1 but we include a short
summary of points below.

  * •

The $1\,(\text{Gpc}/h)^{3}$ MultiDark (Klypin et al., 2016) $N$-body semi-
analytical simulation with approximate cold gas masses is used for the single
realisation of the underlying true Hi signal at a snapshot redshift of
$z\,{=}\,0.39$, gridded into
$n_{\text{x}},n_{\text{y}},n_{\text{z}}\,{=}\,256,256,256$ voxels. We include
redshift-space distortions (RSD) to provide the Hi signal with an anisotropic
signature. A frequency range of $900\,{<}\,\nu\,{<}\,1156\,\text{MHz}$ with
resolution $\delta\nu\,{=}\,1\,\text{MHz}$ is assumed which is consistent with
the snapshot redshift and is reasonably consistent with MeerKAT L-band
intensity mapping observations.

  * •

We simulate galactic synchrotron, galactic free-free, and bright point source
emission at these frequencies to provide a foreground sky. We use the Planck
Legacy Archive222pla.esac.esa.int/pla FFP10 simulations for the synchrotron
and free-free emission. The point sources catalogue is produced following the
same approach as in Battye et al. (2013).

  * •

We cut a patch of sky consistent with the $1\,(\text{Gpc}/h)^{2}$ Hi survey
size and chose this to be centered on the Stripe 82 region of the sky, where a
real survey could be targeted. The foreground component with
$n_{\theta}\,{=}\,n_{\text{x}}\,{\times}\,n_{\text{y}}$ angular pixels and
$n_{\nu}\,{\equiv}\,n_{\text{z}}$ frequency channels is added onto the Hi
signal simulation.

  * •

To increase the complexity of the foreground clean, we simulate instrumental
polarisation leakage (Carucci et al., 2020; Cunnington et al., 2021a) which
disrupts the smooth frequency spectra of the foreground simulations, requiring
a clean which is more aggressive and consistent with real data. For this we
used the CRIME333intensitymapping.physics.ox.ac.uk/CRIME.html software (Alonso
et al., 2014). This is used by default and we highlight any cases where this
has not been used.

  * •

By default we add no further instrumental effects, but in some cases we
introduce instrumental noise and smoothing perpendicular to the line-of-sight
to emulate the telescope beam. For the noise, we assume isotropic Gaussian
white noise with $\sigma_{\text{n}}\,{=}\,1\,\text{mK}$, approximately
corresponding to $30\,\text{hrs}$ of observation time on a MeerKAT-like survey
of ${\sim}\,3{,}000\,\text{deg}^{2}$ sky (see Equation 25 and 26 for more
details). This noise dominates over the Hi signal which has an rms of
$\sigma_{\textrm{H{i}}}\,{\sim}\,0.14\,\text{mK}$. The beam we approximate as
a Gaussian with comoving transverse length scale
$R_{\text{beam}}\,{=}\,10\,h^{-1}\text{Mpc}$ (see Equation 24 for a
definition). We clearly indicate cases where noise or a beam has been added.
We discuss the limitations of these approximations and also explore some more
realistic systematics in Section 4 based on real MeerKAT pathfinder data,
which we introduce there.

Throughout we will refer to these as the MD1GPC simulations. Later in the
paper, we use some more specific simulations to explore different scenarios
which we will introduce then, but for the majority of our results, we use the
MD1GPC by default unless otherwise clearly stated.

To begin a PCA-based clean of the Hi \+ foregrounds combination, we first
calculate the $\nu,\nu^{\prime}$ covariance of the foreground contaminated
data $\textbf{{X}}_{\text{obs}}$, where the data matrices X have dimensions
$[n_{\nu},\,n_{\theta}]$. The covariance is estimated by
$\textbf{{C}}\,{=}\,(n_{\theta}\,{-}\,1)^{-1}\textbf{{X}}_{\text{obs}}^{\text{T}}\textbf{{X}}_{\text{obs}}\,{=}\,\textbf{{U}}\boldsymbol{\Lambda}\textbf{{U}}^{\text{T}}$,
where the last equality is the eigen-decomposition (or diagonalisation) of the
covariance matrix, with U representing a matrix with the $n_{\nu}$ spectral
eigenvectors $\textbf{{U}}_{i}$ and $\boldsymbol{\Lambda}$ is the diagonal
matrix of eigenvalues. Neglecting noise contributions,444This is mainly done
for brevity and incorporating noise into this formalism is not overly
complicated. It acts as additional perturbations to the pure foreground modes
in a similar way to the Hi signal, as we later discuss. i.e.
$\textbf{{X}}_{\text{obs}}\,{=}\,\textbf{{X}}_{\text{f+s}}\,{\equiv}\,\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{s}}$,
we can write this as

$\textbf{{C}}=(n_{\theta}\,{-}\,1)^{-1}\left(\textbf{{X}}_{\text{f}}+\textbf{{X}}_{\text{s}}\right)^{\text{T}}\left(\textbf{{X}}_{\text{f}}+\textbf{{X}}_{\text{s}}\right)\,,$
(1)

which expands to

$\displaystyle\textbf{{C}}_{\mathrm{f}+\mathrm{s}}$
$\displaystyle=(n_{\theta}\,{-}\,1)^{-1}\left(\textbf{{X}}_{\mathrm{f}}\textbf{{X}}_{\mathrm{f}}^{\text{T}}+\textbf{{X}}_{\mathrm{f}}\textbf{{X}}_{\mathrm{s}}^{\text{T}}+\textbf{{X}}_{\mathrm{s}}\textbf{{X}}_{\mathrm{f}}^{\text{T}}+\textbf{{X}}_{\mathrm{s}}\textbf{{X}}_{\mathrm{s}}^{\text{T}}\right)$
(2) $\displaystyle=\textbf{{C}}_{\mathrm{f}}+\textbf{{C}}_{\Delta}\,,$

where
$\textbf{{C}}_{\Delta}\,{=}\,(n_{\theta}\,{-}\,1)^{-1}(\textbf{{X}}_{\mathrm{f}}\textbf{{X}}_{\mathrm{s}}^{\text{T}}+\textbf{{X}}_{\mathrm{s}}\textbf{{X}}_{\mathrm{f}}^{\text{T}}+\textbf{{X}}_{\mathrm{s}}\textbf{{X}}_{\mathrm{s}}^{\text{T}})$
are residual contributions to the estimate of the foreground covariance. The
estimate of the foregrounds, which is to be removed from the observations, is
then given by

$\hat{\textbf{{X}}}_{\text{f}}=\textbf{{U}}\textbf{{S}}\textbf{{U}}^{\text{T}}\textbf{{X}}_{\text{obs}}\,,$
(3)

where, following Switzer et al. (2015), we introduce the selection matrix S
which is zero everywhere except for the first $N_{\text{fg}}$ elements along
the diagonal which are set equal to one, assigning the number of contaminated
modes projected out from each line of sight.

Whilst we expect the foregrounds to be orders of magnitude larger than the Hi
signal, and $\textbf{{C}}_{\mathrm{f}}$ to be the dominant term in Equation 2,
it is important for the additional perturbations from the signal through
$\textbf{{C}}_{\Delta}$ to be considered. Due to this mix of foregrounds and
signal, the eigenvectors we obtain are

$\textbf{{U}}\equiv\textbf{{U}}_{\mathrm{f}+\mathrm{s}}=\textbf{{U}}_{\mathrm{f}}+\boldsymbol{\Delta}\,,$
(4)

where $\textbf{{U}}_{\mathrm{f}}$ are pure foreground modes and
$\boldsymbol{\Delta}$ are the perturbed contributions caused by the signal. We
show the difference between the unperturbed ($\textbf{{U}}_{\text{f}}$) and
perturbed ($\textbf{{U}}_{\text{f+s}}$) eigenvectors in Figure 1. This shows
the first 9 most dominant eigenvectors from the MD1GPC simulations. For the
unperturbed, pure foreground eigenvectors (dashed black lines), the modes are
smooth in frequency, only showing longer wavelength oscillations caused by the
polarisation leakage. However, for the eigenvectors estimated from the
foreground and signal mix (solid coloured lines), perturbations to the modes
caused by the signal start to arise. These perturbations are more severe the
higher the order of the eigenvector. This is because the eigenvectors become
increasingly less dominant and are more easily perturbed by the presence of
the signal whose contributions remain fairly consistent through all the
eigenvectors due to its high-rank properties.

Figure 1: The first 9 eigenvectors (first 5 in the top panel and next four in
the bottom) from the PCA on the MD1GPC foreground contaminated simulations.
The dashed-black lines show purely foreground modes taken from the PCA on
foreground-only sims. The solid colour lines show the eigenvectors perturbed
by the inclusion of the Hi signal in the data. These perturbations are the
origins of signal loss in foreground cleaning. The thin-faint solid color
lines show eigenvectors perturbed by the inclusion of the Hi signal + noise.

Figure 1 begins to demonstrate how signal loss can enter in a foreground
clean. By using the eigenvectors $\textbf{{U}}_{\text{f{+}s}}$ as the basis
functions which are projected out in the PCA clean, this will project out
modes that have some Hi structure shown by the perturbations on the lower
modes. It is tempting to try and address this issue at this early stage and
use a low-pass filter or smooth the perturbed modes to correct the
perturbations from the signal. However, we briefly experimented with various
smoothing routines with this aim and in all cases, the results were made
worse. Since the aim of this work is not to enhance foreground cleaning
efficiency, but instead to ensure we can control signal loss, we defer any
investigations into cleaning optimisation to future work.

We also show in Figure 1 the impact on the eigenmodes by including the
dominant instrumental noise ($\textbf{{U}}_{\text{f+s+n}}$ shown by faint
colored lines). These create much larger perturbations to the pure foreground
modes, hence large noise can impact foreground cleaning. This is an important
issue for early pilot surveys where observation time is low since in these
cases, the noise will dominate over the Hi and will be the main source of
perturbations to the eigenvectors. We will discuss this in more detail later
and demonstrate how intensity maps with a high level of noise, and other
additive systematics like residual RFI, can still undergo signal
reconstruction using a transfer function. For the remainder of this section,
we omit the instrumental noise for simplicity.

The perturbations in Figure 1 are dependent on the ratio between the
foreground amplitude and the other components e.g. the Hi. Whilst the Hi
signal amplitude will be consistently uniform due to the cosmological
principle, the foreground emission can vary with the choice of sky patch, e.g.
being orders of magnitude higher near the galactic plane relative to the South
Celestial Pole. This was explored in Cunnington et al. (2021a) where different
foreground regions were tested. For the remainder of this work, however, we
will stick to one region as most of the conclusions we draw are generic
regardless of how strong the foregrounds are, within a physically reasonable
range.

###  2.1 Toy model foreground cleaning

Here we investigate some idealised foreground cleaning scenarios to
demonstrate the nature of signal loss in blind foreground cleaning. We begin
with the most ideal toy case scenario where we project out pure foreground
modes from pure foreground-only data and subtract this from the observed
combination $\textbf{{X}}_{\text{f+s}}$. Of course, if we could access perfect
foreground-only data $\textbf{{X}}_{\text{f}}$ this could be simply subtracted
from the observed foreground and signal mix, so there would be no need for any
mode projection cleaning process. However, we proceed with this example since
it provides valuable insight from which we can add complication. This first
toy-case is given by

$\textbf{{X}}_{\text{toy:clean1}}=\textbf{{X}}_{\text{f+s}}-\textbf{{U}}_{\text{f}}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{f}}\,,$
(5)

In this ideal case, since we are projecting out perfect foreground modes from
pure foreground data, the optimal selection matrix S would have the diagonals
filled with ones ($N_{\text{fg}}\,{=}\,n_{\nu}$) i.e. the identity matrix, to
remove all foreground without the consequence of signal loss. In reality, this
is not possible and a balance is sought between projecting out enough modes to
remove foregrounds but not so many that large signal loss is sustained.

Figure 2: Auto power spectra for different foreground cleaning cases in
comparison with the true Hi-only simulations (black-dashed line). The first
three panels represent idealised scenarios. Far-left shows projecting pure
foreground modes ($\textbf{{U}}_{\text{f}}$) out from pure foreground data
($\textbf{{X}}_{\text{f}}$). Second panel shows projecting pure foreground
modes out from foreground data mixed with signal
($\textbf{{X}}_{\text{f+s}}$). Third panel shows projecting modes perturbed by
the signal ($\textbf{{U}}_{\text{f+s}}$) out from pure foreground data. Far-
right shows the realistic scenario where modes are perturbed by the signal and
these are then projected out from foreground data mixed with signal.

In Figure 2 we show the measured power spectrum for the idealised toy clean in
Equation 5 (far-left panel). We also show other cases of foreground clean in
the other panels which we will discuss shortly. Details on the power spectrum
estimation formalism used throughout the paper are presented in Appendix B. We
show the power spectra in comparison with the original Hi-only (foreground-
free) data, which we are aiming to agree with. In all cases we show two
cleaned examples with 8 and 12 modes projected out, i.e. the first
$N_{\text{fg}}\,{=}\,8$ and $N_{\text{fg}}\,{=}\,12$ elements in the diagonal
of S are set to 1. This is why we do not reach perfect agreement in the ideal
top-left panel, because we are not projecting out all pure foreground modes,
leaving some foreground residual in the remaining modes. Thus, there is slight
disagreement with the Hi-only, albeit at a sub-percent level for the
$N_{\text{fg}}\,{=}\,12$ case. The values of $N_{\text{fg}}\,{=}\,8,\,12$ are
chosen to sufficiently suppress the polarised foregrounds in the increasingly
more realistic cases shown by the other panels which we discuss next.

In reality, the situation in Equation 5 is not possible, because we can not
project out modes from foreground-only data $\textbf{{X}}_{\text{f}}$ because
the observed data we have, $\textbf{{X}}_{\text{f+s}}$, is inherently mixed
with signal. Signal loss begins to manifest in the case where we project out
foreground modes from the observed data mix

$\textbf{{X}}_{\text{toy:clean2}}=\textbf{{X}}_{\text{f+s}}-\textbf{{U}}_{\text{f}}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{f+s}}\,.$
(6)

This is still an idealised scenario since we are projecting out pure
foreground modes from the data. In reality, a further complication arises
since the modes identified in the PCA will be perturbed by the presence of
signal i.e.
$\textbf{{U}}_{\text{f}}\,{\rightarrow}\,\textbf{{U}}_{\text{f+s}}$, which we
will discuss shortly. Comparing the first panel with the second, where the
difference is that in the former case, pure foreground modes are now being
projected out from the mix of foreground and signal (Equation 6), evidence of
signal loss in the cleaned power spectrum begins to show. The signal loss is
clearly dominating over any foreground residuals remaining in the data from
only projecting out a finite number of modes. In other words, the small
${\sim}\,5\%$ additive bias in the far-left panel caused by foreground
residuals is not seen in the second panel, due to a more dominant impact from
signal loss. Of course if a smaller number for $N_{\text{fg}}$ were chosen,
foreground residuals would cause more of a problem.

The second panel of Figure 2 confirms that signal loss begins to manifest when
modes (even purely foreground ones) are projected out of the data, which is a
combination of foregrounds and signal. The reason for this is because signal
will unavoidably have degeneracies with some foreground structure. Thus when a
set of foreground functions are projected out of the data, some signal will
leak into this subtraction, mainly large-scale (small-$k$) modes since these
are most degenerate with the foreground structure which is highly correlated
through frequency.

The third panel of Figure 2 shows a final idealised toy case where we are only
projecting out modes from the pure foreground map, but the modes are now
perturbed by the presence of signal, $\textbf{{U}}_{\text{f+s}}$ (see Equation
4). This is something we have to deal with in reality where we can not form a
perfect estimation for the foreground-only eigenmodes
$\textbf{{U}}_{\text{f}}$, from the true observed data where foreground and
signal are mixed. In this case, the eigenvectors themselves are perturbed and
it is these perturbed modes that we project along the data;

$\textbf{{X}}_{\text{toy:clean3}}=\textbf{{X}}_{\text{f+s}}-\textbf{{U}}_{\text{f+s}}\textbf{{S}}\textbf{{U}}_{\text{f+s}}^{\text{T}}\textbf{{X}}_{\text{f}}\,.$
(7)

This provides an interesting result with signal loss again appearing to be the
more dominant effect with little evidence of additive bias from foreground
residuals. However we are only projecting out modes from pure foreground data,
so it seems counter-intuitive that there is signal loss. As we will explicitly
show in the following section, this is caused by the perturbation to the modes
from the presence of signal ($\textbf{{U}}_{\text{f+s}}$), which creates a
complicated mix of subtracted terms that can have signal correlating and anti-
correlating contributions, as identified in Switzer et al. (2015).

###  2.2 The origins and subtleties of signal loss

Despite seeing signal loss in the second and third panels of Figure 2, both
cases still represent unrealistic scenarios. In reality, we see a combination
of both where the presence of signal perturbs the eigenmodes
($\textbf{{U}}_{\text{f+s}}$) as well as complicating the clean since
information is projected out from data which contain not just foregrounds, but
signal too
($\textbf{{U}}_{\text{f+s}}\textbf{{S}}\textbf{{U}}_{\text{f+s}}^{\text{T}}\textbf{{X}}_{\text{f+s}}$).
Thus, the true resulting cleaned data is given by

$\textbf{{X}}_{\text{clean}}=\textbf{{X}}_{\text{f+s}}-\textbf{{U}}_{\text{f+s}}\textbf{{S}}\textbf{{U}}_{\text{f+s}}^{\text{T}}\textbf{{X}}_{\text{f+s}}\,.$
(8)

The result from this foreground clean is shown in the final far-right panel of
Figure 2. Results appear similar to $\textbf{{X}}_{\text{toy:clean2}}$ but
some differences can be seen on large scales caused by the increased
complexity of having perturbed eigenmodes ($\textbf{{U}}_{\text{f+s}}$).

To understand the complexity of foreground cleaning, we expand the above
Equation 8 into all its terms, giving

$\displaystyle\textbf{{X}}_{\text{clean }}$
$\displaystyle=\left[1-\left(\textbf{{U}}_{\mathrm{f}}+\boldsymbol{\Delta}\right)\textbf{{S}}\left(\textbf{{U}}_{\mathrm{f}}+\boldsymbol{\Delta}\right)^{\text{T}}\right]\left(\textbf{{X}}_{\mathrm{f}}+\textbf{{X}}_{\mathrm{s}}\right)$
(9)
$\displaystyle=\textbf{{X}}_{\text{f}}+\textbf{{X}}_{\text{s}}-\textbf{{U}}_{\mathrm{f}}\textbf{{S}}\textbf{{U}}_{\mathrm{f}}^{\text{T}}\textbf{{X}}_{\text{f}}-\textbf{{U}}_{\mathrm{f}}\textbf{{S}}\textbf{{U}}_{\mathrm{f}}^{\text{T}}\textbf{{X}}_{\text{s}}-\Delta\textbf{{S}}\textbf{{U}}_{\mathrm{f}}\textbf{{X}}_{\mathrm{f}}-\textbf{{U}}_{\mathrm{f}}\textbf{{S}}\Delta^{\text{T}}\textbf{{X}}_{\mathrm{f}}$
$\displaystyle\quad-\boldsymbol{\Delta}\textbf{{S}}\boldsymbol{\Delta}^{\text{T}}\textbf{{X}}_{\text{s}}-\boldsymbol{\Delta}\textbf{{S}}\boldsymbol{\Delta}^{\text{T}}\textbf{{X}}_{\text{f}}-\boldsymbol{\Delta}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{s}}-\textbf{{U}}_{\text{f}}\textbf{{S}}\boldsymbol{\Delta}^{\text{T}}\textbf{{X}}_{\text{s}}\,.$

In Figure 3 we show power spectra for the subtracted decomposed terms in
Equation 9, plotting their cross-correlation with the pure Hi signal to
demonstrate where signal loss originates. Since these are the subtracted
terms, the higher their cross-power with pure-Hi, the more they are
contributing to signal loss. For reference, we also show the pure-Hi (i.e. the
Hi auto-correlation) as the black-dotted line, and the fully cleaned result
(all terms from Equation 9 combined) as the grey dotted line.

Figure 3: Contributions to signal loss from the decomposed terms in Equation
9. Each power spectra shows the cross-correlation between the original-Hi and
a subtracted term in the foreground clean. The different panels represent
different numbers of PCA modes removed ($N_{\text{fg}}$). Solid lines indicate
residual foreground contributions whereas dashed lines indicate Hi signal
contributions. Thin lines indicate perturbed contributions from
$\mathbf{\Delta}$ caused by the presence of signal in the eigenmode
estimation. For reference, we also include the pure Hi-power (black-dotted)
along with the total cleaned result $\textbf{{X}}_{\text{clean}}$ (grey-
dotted).

As expected, a large bulk of signal loss is caused by the subtraction of the
$\textbf{{U}}_{\text{f}}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{s}}$
term (blue dashed line). This direct signal loss will clearly increase for
higher $N_{\text{fg}}$ i.e. more ones along the diagonal of S, and this is
demonstrated by the growing amplitude of the blue dashed line mostly at
small-$k$, going left to right in the panels. The projected out foregrounds
$\textbf{{U}}_{\text{f}}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{f}}$
are entirely uncorrelated with the Hi signal as shown by the consistent with
zero power spectrum (orange line). However, we still decrease the amplitude of
$\textbf{{U}}_{\text{f}}\textbf{{S}}\textbf{{U}}_{\text{f}}^{\text{T}}\textbf{{X}}_{\text{f}}$
by three orders of magnitude (indicated in the legend) since this dominant
term still has large purely statistical fluctuations in the Hi cross-
correlation dependent only on the foreground realisation. The thinner lines
represent terms including perturbative contributions from the Hi signal
$\mathbf{\Delta}$. This is where the issue of signal loss begins to
complicate. As $N_{\text{fg}}$ increases, the contribution to signal loss from
$\mathbf{\Delta}\textbf{{S}}\mathbf{\Delta}^{\text{T}}\textbf{{X}}_{\text{s}}$
(brown dashed line) becomes non-negligible. Complicating matters further is
the removal of the anti-correlating contribution in
$\mathbf{\Delta}\textbf{{S}}\textbf{{U}}^{\text{T}}_{\text{f}}\textbf{{X}}_{\text{s}}$.
Lastly, there are also noticeable correlations in the perturbed foreground
removed terms. The thin red line shows how the removed term
$\mathbf{\Delta}\textbf{{S}}\textbf{{U}}^{\text{T}}_{\text{f}}\textbf{{X}}_{\text{f}}$
will also introduce a contribution to signal loss. This explains the previous
result presented in the bottom-left panel of Figure 2 where, despite only
projecting out modes from pure foreground data $\textbf{{X}}_{\text{f}}$,
signal loss still arose in the cleaned power spectrum, albeit at a
${\sim}\,5\%$ level in the largest scales. This will be caused by the signal
perturbations $\mathbf{\Delta}$ which introduce a small correlation with the
Hi signal, shown by the red line in Figure 3.

The complex mix of signal correlation and anti-correlation caused by the
perturbations $\mathbf{\Delta}$ from non-foreground modes can clearly affect
the overall signal loss. The impact of signal perturbations on the foreground
modes becomes increasingly more important the more aggressive the foreground
clean due to $\mathbf{\Delta}$ having more influence over less dominant
foreground modes, as seen in Figure 1. It is therefore crucial to model or
emulate all these contributions in any signal reconstruction to avoid unbiased
results as highlighted in previous literature (Switzer et al., 2015; Cheng et
al., 2018). In Section 3 we will explore how signal injection can be used to
construct a foreground transfer function and validate with simulations how it
is able to successfully emulate all the subtle contributions to signal loss.
We will also explicitly highlight cases where a transfer function can be
incorrectly assembled such that some of the contributions demonstrated by
Figure 3 are not accounted for, leading to incorrect estimations of signal
loss.

##  3 Validating the transfer function with simulations

As demonstrated in the previous section, signal loss from blind foreground
cleaning is complicated by the subtraction of sub-sets of data that have
spurious correlations and anti-correlations with the Hi signal. The spurious
correlations arise because the estimated modes projected out in the blind
foreground clean are perturbed by the presence of the Hi signal itself. Thus
the signal loss is dependent on the specific realisation of the signal,
foregrounds, and their combination. Modelling the signal loss, or measuring it
in pure simulations, is therefore potentially problematic and could lead to
biased results. In this section, we explore how we can utilise the observed
contaminated data itself to emulate the complex spurious correlations in
injected mock data and use the signal loss experienced in the mocks to
construct a foreground transfer function. This data-driven approach has been
extensively used in single-dish intensity mapping detections (Masui et al.,
2013; Anderson et al., 2018; Wolz et al., 2022; Cunnington et al., 2022). Here
we present the formalism for an unbiased application of the transfer function
and for the first time validate its performance on simulations whilst also
trying to expose any limitations.

Throughout this section, we treat the MD1GPC simulation as the observed
intensity mapping data, and use lognormal mocks as completely separate
simulations in the construction of the transfer functions. This maintains a
certain independence between the injected mocks and the simulated signal that
we are trying to recover, as would be the case in real observations. There is
also the option of generating more complex mocks to inject into the data, for
example using field level forward modelling (Obuljen et al., 2022) or a Hi-
halo prescription as trialed in Wolz et al. (2022), however we found lognormal
mocks sufficient for our purposes.

In this section, in cases where we are investigating the cleaned or
reconstructed power spectrum, unless stated otherwise, we will use the cross-
correlation power spectrum between the foreground cleaned MD1GPC map and the
Hi-only (foreground-free) MD1GPC map. This is so foreground residuals will be
less of an issue and their additive bias does not confuse the investigation of
signal loss and reconstruction accuracy. In cross-correlation with the Hi-only
map, any difference relative to the original Hi-only auto-power spectrum will
be caused solely by signal loss.

###  3.1 Summarised recipe for the transfer function and its unbiased results

We begin by providing a summary of how a transfer function can be used for
correcting signal loss from foreground cleaning, and validate the performance
of the process. We then go into more detail in the remainder of this section,
clarifying exactly how the transfer function can be constructed and used for
various scenarios. In short, the foreground transfer function is constructed
by injecting mock data into the observed maps. Then, by running the same
foreground removal routine, one will subject the mock data to a similar signal
loss that is experienced by the true underlying Hi signal, thus giving an
estimate of the true signal loss.

Below is the step-by-step recipe for how to construct and apply an unbiased
transfer function;

  1. (i)

Firstly, foreground clean the observed data $\textbf{{X}}_{\text{obs}}$ by
projecting out $N_{\text{fg}}$ PCA modes i.e.
$\textbf{{X}}_{\text{clean}}\,{=}\,\textbf{{X}}_{\text{obs}}\,{-}\,\textbf{{U}}\textbf{{S}}\textbf{{U}}^{\text{T}}\textbf{{X}}_{\text{obs}}$,
where U is a matrix of eigenvectors from the diagonalisation of the
$\nu,\nu^{\prime}$ covariance matrix estimated empirically from the data, and
S is the selection matrix with ones along the first $N_{\text{fg}}$ diagonal
elements and zeros elsewhere.

  2. (ii)

Compute the power spectrum for the foreground-cleaned data, which will be
negatively biased due to the signal loss from the foreground clean. The power
spectrum is given here by
$P_{\text{clean}}(k)\,{=}\,\mathcal{P}(\textbf{{X}}_{\text{clean}},\textbf{{X}}_{\text{tr}})$,
where $\mathcal{P}(\textbf{{X}}_{\text{clean}},\textbf{{X}}_{\text{tr}})$ is
an operator which measures the cross-power spectrum between data sets
$\textbf{{X}}_{\text{clean}}$ and $\textbf{{X}}_{\text{tr}}$, then reduces the
power into the spherically-averaged $k$-bins. Here, $\textbf{{X}}_{\text{tr}}$
is any overlapping tracer, which can be the intensity map itself for an auto-
correlation survey, or a galaxy survey as a common example of cross-
correlation.

  3. (iii)

Generate mock Hi signal maps $\textbf{{X}}_{\text{m}}$ with the same
dimensions as the observed data. In this work we use a fast lognormal
transform process to generate mocks from the Hi power spectrum model given in
Equation 30. We investigate the consequences of variation in the input mocks
in Section 5.

  4. (iv)

Emulate signal loss in the mock by injecting it into the real data and
foreground cleaning the observed data and mock combination,

$\textbf{{X}}^{\text{m}}_{\text{clean}}\,{=}\,(\textbf{{X}}_{\text{obs}}\,{+}\,\textbf{{X}}_{\text{m}})\,{-}\,\textbf{{U}}_{\text{obs+m}}\textbf{{S}}\textbf{{U}}_{\text{obs+m}}^{\text{T}}(\textbf{{X}}_{\text{obs}}\,{+}\,\textbf{{X}}_{\text{m}})\,{-}\,[\textbf{{X}}_{\text{clean}}]\,.$
(10)

The term in the square brackets is subtracting the cleaned observed data
without mocks to remove contributions in the map uncorrelated to the mock
signal thus reducing the variance of the transfer function, as we will
explicitly demonstrate.

  5. (v)

The foreground transfer function is then given by

$\mathcal{T}(k)=\left\langle\frac{\mathcal{P}(\textbf{{X}}^{\text{m}}_{\text{clean}}\,,\,\textbf{{X}}_{\text{m}})}{\mathcal{P}(\textbf{{X}}_{\text{m}}\,,\,\textbf{{X}}_{\text{m}})}\right\rangle_{N_{\text{mock}}}\,,$
(11)

where the angled brackets denote an averaging over iterations of a suitably
large number of mocks ($N_{\text{mock}}$) until a converged transfer function
is achieved. We use $N_{\text{mock}}\,{=}\,100$ by default unless otherwise
mentioned.

  6. (vi)

De-bias the cleaned power spectrum using the transfer function to reconstruct
the signal loss with
$P_{\text{rec}}(k)\,{=}\,P_{\text{clean}}(k)[\mathcal{T}(k)]^{-1}$. Note the
index of $-1$ should also be used in auto-correlation i.e. an auto-correlation
of an intensity map should not have signal loss corrected for twice, as we
will demonstrate in Section 3.4.

  7. (vii)

The covariance of the reconstructed power spectrum can also be extracted from
the mocks used in the transfer function calculation. Whilst the mean of
$P_{\text{clean}}\mathcal{T}^{-1}_{i}$ over all $N_{\text{mock}}$ iterations
provides the reconstructed power spectrum, the covariance estimates the errors
inclusive of signal loss uncertainty. However, as we will show, it is crucial
not to subtract the square-bracket $\textbf{{X}}_{\text{clean}}$ term in
Equation 10 when estimating the covariance, as this will include foreground
residuals, instrumental noise, etc. all of which should contribute to the
error budget. We discuss this in detail in Section 3.3.

The numerator in Equation 11 is taking the cross-correlation between the
cleaned mock $\textbf{{X}}^{\text{m}}_{\text{clean}}$ and original mock
$\textbf{{X}}_{\text{m}}$ with no cleaning effects. This should therefore not
be overly influenced by foreground residuals and differences between this
cross-correlation and the auto-correlation in the denominator should only be
caused by signal loss from the foreground clean, thus their ratio provides the
level of the original signal remaining in the power spectrum of the cleaned
mock $\textbf{{X}}^{\text{m}}_{\text{clean}}$. The crucial part for Equation
11 is having a process for obtaining $\textbf{{X}}^{\text{m}}_{\text{clean}}$
such that the signal loss it experiences across all scales, is the same as the
signal loss in the actual data $\textbf{{X}}_{\text{s}}$. To achieve this we
inject mock signal $\textbf{{X}}_{\text{m}}$ into the observed data with
foregrounds and true Hi signal
($\textbf{{X}}_{\text{f+s+m}}\,{\equiv}\,\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{s}}\,{+}\,\textbf{{X}}_{\text{m}}$)
then project out the same number of modes as in the original foreground clean
of the observations i.e.;

$\textbf{{X}}^{\text{m}}_{\text{clean}}=\textbf{{X}}_{\text{f+s+m}}-\textbf{{U}}_{\text{f+s+m}}\textbf{{S}}\textbf{{U}}_{\text{f+s+m}}^{\text{T}}\textbf{{X}}_{\text{f+s+m}}-[\textbf{{X}}_{\text{clean}}]\,.$
(12)

This is equivalent to what we presented in the summarised recipe in Equation
10. As discussed, the term in the square bracket is subtracting the cleaned
observed data (with no mock injection) to reduce the transfer function
variance, which we discuss in more detail later. The presence of mock signal
will cause perturbations to the eigenmodes, and will emulate the signal loss
coming from both projecting out the modes with signal perturbations and the
complex correlations between all the cross terms discussed in Section 2.2 and
Equation 9.

The presence of the true observed Hi signal $\textbf{{X}}_{\text{s}}$ in
$\textbf{{X}}_{\text{f+s+m}}\,{\equiv}\,\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{s}}\,{+}\,\textbf{{X}}_{\text{m}}$
in Equation 12 creates unwanted complications to the transfer function
construction which should ideally only be concerned with the mock signal
$\textbf{{X}}_{\text{m}}$ and its relationship with $\textbf{{X}}_{\text{f}}$.
The presence of $\textbf{{X}}_{\text{s}}$ will not matter from a direct signal
loss perspective, since this will not affect the cross-correlation with
$\textbf{{X}}_{\text{m}}$ in Equation 11. However, $\textbf{{X}}_{\text{s}}$
will perturb the estimation of the eigenvectors. This is unwanted because we
ideally only want the mock signal to perturb the eigenvectors and just produce
$\textbf{{U}}_{\text{f+m}}$, but by injecting mocks into the true signal we
will actually measure

$\textbf{{U}}_{\text{f+s+m}}=\textbf{{U}}_{\text{f}}+\boldsymbol{\Delta}_{\text{s}}+\boldsymbol{\Delta}_{\text{m}}\,,$
(13)

where we have introduced the subscripts m and s to the perturbations
$\mathbf{\Delta}$ to distinguish perturbations from mock signal and true Hi
signal respectively. The two sources of perturbation is not seen in the
foreground clean of just the observations
($\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{s}}$). In other words the
eigenvectors are now being perturbed twice. As we will show from our results
shortly, this appears to have little impact and we still obtain an unbiased
transfer function. We tested the transfer function in an idealised case where
mock signal was injected into just pure foreground
($\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{m}}$) and found little
difference in performance compared to the realistic case where true signal is
present
($\textbf{{X}}_{\text{f}}\,{+}\,\textbf{{X}}_{\text{s}}\,{+}\,\textbf{{X}}_{\text{m}}$).

The accuracy of the reconstructed power spectra is demonstrated in Figure 4.
The simulated observations are cleaned by removing either
$N_{\text{fg}}\,{=}\,8$ or $12$ PCA modes, then the transfer function is used
to correct for the signal loss, with the reconstructed result being divided by
the original foreground-free simulation $P_{\textrm{H{i}}}(k)$. Thus, a
perfect reconstruction would give unity across all scales. We see excellent
performance with sub-percent accuracy achieved across most scales above
$k\,{>}\,0.1\,h\,\text{Mpc}^{-1}$ for the $N_{\text{fg}}\,{=}\,8$ case.
Performance is still good for the $N_{\text{fg}}\,{=}\,12$ case, albeit with a
noticeable drop in accuracy relative to $N_{\text{fg}}\,{=}\,8$ mostly at
large scales (small-$k$). This will be caused by the increased effect from
spurious correlations between foregrounds and mock signal, which, as we
demonstrated in Figure 3, increases for more aggressive (higher
$N_{\text{fg}}$) foreground cleans. This will not necessarily bias the results
since the variance in the transfer function also increases for higher
$N_{\text{fg}}$, as shown by the shaded regions, thus can be reflected in the
error estimations (discussed in a later section).

Figure 4: Accuracy of reconstructed foreground cleaned power spectra relative
to the foreground-free Hi-only data ($P_{\textrm{H{i}}}$) from simulated
intensity maps. The foreground cleaned power spectrum has been reconstructed
using the transfer function to correct for the signal loss from foreground
cleaning. The transfer functions are calculated using Equation 11 and 12 and
averaging over 100 lognormal mocks. The shaded bands show the rms over these
100 mocks. Results for a mild ($N_{\text{fg}}\,{=}\,8$, blue lines) and more
aggressive ($N_{\text{fg}}\,{=}\,12$, red lines) foreground cleans are shown.
blackDark-thick (light-thin) green horizontal lines indicate sub 1% (5%)
accuracy regions of the reconstructed power spectrum.

In general on large ($k\,{<}\,0.1\,h\,\text{Mpc}^{-1}$) scales, we see a less
reliable result in terms of pure accuracy, but this is also accounted for by
the transfer function variance, which can reach ${\gtrsim}\,5\%$ on these
scales. The performance at large scales is however dependent on the size of
the intensity mapping survey. The depth of the $1\,(\text{Gpc}/h)^{3}$ MD1GPC
simulation at $z\,{\sim}\,0.39$ is reasonably consistent with a MeerKAT L-band
survey, assuming it uses the complete band range ($0.2\,{<}\,z\,{<}\,0.58$).
However, future surveys in UHF band, and then eventually the SKAO, will cover
much wider frequency ranges. This will mean reconstructed modes at
$k\,{<}\,0.1\,h\,\text{Mpc}^{-1}$ become more reliable due to a suppression of
sample variance and less signal loss which will now be contained to even
larger scales. We will demonstrate this point later in Section 5 with some
additional simulations which cover a larger volume.

The presence of the true observed Hi signal in the transfer function
calculation will increase its variance because there will be residual true
signal after the foreground clean. This will be uncorrelated from the mocks
and act like noise and increase the variance across all of the mocks being
averaged over. This is why we subtract the cleaned data (the
$\textbf{{X}}_{\text{clean}}$ term in the square brackets of Equation 12),
since this is only contributing variance to the result. We will revisit this
point in the next section where we will demonstrate that the increased
variance coming from $\textbf{{X}}_{\text{clean}}$ can be utilised for error
estimation. With the $\textbf{{X}}_{\text{clean}}$ subtraction, this version
of the transfer function is not only achieving a good accuracy but also a good
uncertainty on most scales, shown by the shaded region.

The validation of the transfer function demonstrated by Figure 4 is an
important result. This is a method for reconstructing signal loss which is
applicable on real data and delivers unbiased results across all scales and
within sub-percent precision across smaller scales where the particular survey
volume allows those modes to be well sampled
($k\,{>}\,0.1\,h\,\text{Mpc}^{-1}$ for the case of the MD1GPC simulations).
The compromise of having to inject mock signal into a combination of both
foreground and true Hi signal is an unavoidable complication, however, there
is no evidence that this causes any bias in the reconstructed power spectrum.
Furthermore, by subtracting the cleaned data ($\textbf{{X}}_{\text{clean}}$
term in the square brackets of Equation 10), we found the increase in variance
relative to an ideal case where no true signal ($\textbf{{X}}_{\text{s}}$) is
present in the transfer function calculation was only ${\sim}\,20\%$. It is
crucial that the form of the transfer function as defined by Equation 11 and
Equation 10 be followed and in Appendix C we explicitly highlight the
consequences of deviating from this prescription, demonstrating the
significant biases caused when different definitions of the transfer function
are used.

In Figure 6 and Figure 6 we demonstrate the shape of the transfer function in
$k$-space and in doing so, analyse where signal loss is most severe. Figure 6
shows transfer functions for different PCA modes removed (given by
$N_{\text{fg}}$). This confirms that signal loss increases with
$N_{\text{fg}}$ and is higher at smaller-$k$, both as expected. We also show
the impact from adding the dominant instrumental noise. Perhaps counter-
intuitively, this causes less signal loss. This is because the noise is the
dominant source of perturbations to the pure foreground modes (as shown by
Figure 1), hence these noise-dominant modes will have less contribution from
the Hi signal and removing them causes less signal loss. However, this will
result in a poorer overall foreground clean. Thus foreground residuals and the
high-level noise already present will cause problems for additive biases in
auto-correlation and would also lead to higher errors in cross-correlation. So
the presence of high noise is not beneficial as Figure 6 alone may appear to
suggest. We explore the high noise scenario further in Section 4 where we use
simulations which emulate an early pathfinder-like intensity mapping survey.

Figure 5: Foreground transfer functions $\mathcal{T}(k)$ constructed using
Equation 11 and Equation 10 for the MD1GPC simulations for different numbers
of PCA modes removed (given by $N_{\text{fg}}$). Solid lines indicate noise-
free simulations, thin dashed-lines are for cases where white noise with rms
$\sigma_{\text{n}}\,{=}\,1\,\text{mK}$, which dominates over the Hi, is added
to the simulated observations.

Figure 6: Similar to Figure 6, this shows the foreground transfer function but
now in cylindrical $k_{\perp},k_{\parallel}$ space where
$N_{\text{fg}}\,{=}\,12$ PCA modes have been removed. This is for the noise-
free case.

Similarly to Figure 6, an example transfer function is also shown in Figure 6
but now decomposed into cylindrical contributions in
$k_{\perp},k_{\parallel}$. As already well established, signal loss is
overwhelmingly a function of small-$k_{\parallel}$. However, there is also
some slight $k_{\perp}$ dependence with signal loss being slightly higher at
small-$k_{\perp}$ caused by the large angular structures in the foregrounds.
It is important to highlight that the nature of signal loss will vary
depending on not just the foreground’s strength and spectral smoothness, but
also on the depth of the survey in frequency. This means that the signal loss
presented in Figure 6 and Figure 6 is specific to the MD1GPC simulation.
However, the conclusions we have drawn from this are still mostly generic. For
example, signal loss is still contained in small-$k_{\parallel}$ modes in more
systematic dominated intensity maps as shown in recent MeerKAT analysis
(Cunnington et al., 2022) and as we will show later in the simulations
designed to emulate a small MeerKAT pilot intensity mapping survey. Whilst
signal loss appears widespread throughout all modes in the spherically
averaged Figure 6, where ${>}\,5\%$ signal loss is evident even on small
scales, it is clear from Figure 6 that signal loss does tend to zero (where
$\mathcal{T}(k_{\perp},k_{\parallel})\,{\sim}\,1$) above some $k_{\parallel}$
cut. This raises an intriguing possibility of adopting a hybrid foreground
cleaning/avoidance strategy where a blind foreground clean is run on the full
data, but then an avoidance strategy is used where only modes above some
$k_{\parallel}$ are kept for further analysis. At the very least this would
limit the dependence on the transfer function but would unlikely be reliable
enough to completely avoid any use of signal reconstruction. Furthermore, the
methodology of the transfer function would still be a required tool for
robustly assessing where an optimum cut in $k_{\parallel}$ should be made.
Scale cuts will also limit the scope and constraints possible with the
experiment. We defer any investigations of $k_{\parallel}$ cuts to future work
and continue to test the transfer function on small-$k_{\parallel}$,
especially since these are the scales that will test the performance of the
transfer function most stringently.

###  3.2 1D vs 2D bandpowers in the transfer function

Our results so far have reduced all power directly into 1D spherically
averaged $k$-bins and the results in Figure 4 suggest this can be sufficient.
Providing the same $k$-bins are used in the 1D spherical averaging for the
transfer function construction and cleaned power spectrum, then the transfer
function should encapsulate the same anisotropic signal loss in each $k$-bin
as inflicted on the cleaned data. Furthermore, going straight to 1D $k$-bins
avoids extra compression steps which could potentially lead to results being
lossier. However, it has been proposed in the literature that because the
signal loss is anisotropic (demonstrated by Figure 6) that the transfer
function should be estimated and applied in 2D cylindrical
$k_{\perp},k_{\parallel}$ space, with these bandpowers then re-binned to
provide the final spherically averaged 1D power spectrum (Masui et al., 2013;
Switzer et al., 2015).

We tested the 2D-cylindrical transfer function approach and found evidence of
higher variance in the standard case where complicated polarised foregrounds
are present in the observations. We demonstrate this in Figure 7. Here, the 1D
reconstruction (blue shading) shows our default setup used everywhere else in
the paper, averaging straight into 1D spherical $k$-bins. The 2D
reconstruction refers to a case where we average all power into
$100\,k_{\perp}\times 100\,k_{\parallel}$ cylindrical linear-spaced bins, with
$0\,{<}\,k_{\perp},k_{\parallel}\,{<}\,0.4\,h\,\text{Mpc}^{-1}$, in the
transfer function construction. The measured power for the cleaned
observations is also reduced into the same 2D bins and a reconstructed power
spectrum for a single $i$th mock iteration is given by
$P_{\text{rec},i}(k_{\perp},k_{\parallel})\,{=}\,P_{\text{clean}}(k_{\perp},k_{\parallel})/\mathcal{T}_{i}(k_{\perp},k_{\parallel})$.
The 2D powers then undergo a weighted average into 1D $k$-bins to give the
rebinned 1D power spectrum, defined by

$P_{\text{rec},i}(k)=\frac{\sum_{\alpha}N_{\alpha}P_{\text{rec},i}(k_{\perp,\alpha},k_{\parallel,\alpha})}{\sum_{\alpha}N_{\alpha}}\,,$
(14)

where all unique 2D $k_{\perp},k_{\parallel}$ bandpowers are indexed by
$\alpha$ and the summation is over all 2D powers contained in the 1D bin
$k\,{\equiv}\,\sqrt{k_{\perp}^{2}\,{+}\,k_{\parallel}^{2}}\,{\in}\,(k\,{-}\,\Delta
k/2,k\,{+}\,\Delta k/2)$. $N_{\alpha}$ is the number of 3D Fourier modes
contained in the particular 2D $k_{\perp,\alpha},k_{\parallel,\alpha}$
bandpower. Similar to the direct 1D reconstruction, the mean over all $i$th
mocks in $P_{\text{rec},i}(k)$ gives the final estimated reconstructed power
spectrum, and the variance provides an estimate of the expected errors. The
results for the 2D rebinned power are shown by the red-hatched shading in
Figure 7 where the large increase in variance is clear. We found by switching
off the complexity caused to the foregrounds by the simulated polarisation
leakage made the 2D transfer function more reliable (purple hatched results),
although still a higher variance is returned in these results, particularly at
small-$k$, relative to the 1D case.

Figure 7: Demonstration of the large increase in the variance of the
reconstructed power spectrum, when the transfer function is constructed and
applied in 2D $k_{\perp},k_{\parallel}$ space. Blue results show the standard
1D cases (same as Figure 4) with $N_{\text{fg}}\,{=}\,8$. Red-hatched results
show the 2D construction following steps in Section 3.2. Purple-hatched
results also show results with a 2D construction, but with simplified
foregrounds excluding the polarisation leakage.

It is likely that outliers in isolated iterations are causing the large
variance in Figure 7. These outliers would then be suppressed in the simpler
case where there is no polarisation leakage, thus less complex residuals or
signal loss to cause extreme spurious correlations in the mocks. blackIn an
attempt to suppress the variance, we increased the number of mocks used in the
2D polarisation leakage construction to 500, but this yielded no improvement.
An extension that may be necessary to avoid the blow-up in variance, is a more
detailed weighting to the rebinning procedure we perform in Equation 14. In
Switzer et al. (2015) they discuss applying an inverse covariance weighting to
maximize the 1D signal-to-noise ratio. This could down-weight some of the
outliers in our iterations and suppress the large variance in the 2D polarised
results. We defer this extension to future work, ideally with even more
realistic sims where a conclusive study can be performed into whether the 2D
transfer function construction can be more optimal. Since our results suggest
that direct 1D construction is better performing, at least in the case of the
spherically averaged power spectrum, we blackuse this approach for the rest of
the paper, unless presenting a 2D cylindrical transfer functions which we now
only do for demonstration purposes.

###  3.3 Error estimation for power spectra reconstructed with a foreground
transfer function

Evaluating how to correctly estimate the contributions to the error from
foreground contamination and signal loss uncertainty will be crucial for
future precision cosmology with Hi intensity mapping. This is the focus of
this section. An approach taken in some previous intensity mapping detections
(Anderson et al., 2018) has been to blackuse the variance over the mock
simulations used in the transfer function construction for the error estimate
on cross-correlation measurements. blackIt is possible to capture this
uncertainty from the variance in the transfer function i.e. the errors can be
estimated using

$\hat{\sigma}_{P_{\text{c}}}=\sigma\\{P_{\text{clean}}\mathcal{T}^{-1}_{i}\\}\,,$
(15)

where $\mathcal{T}_{i}$ is the transfer function from the $i$th mock in the
construction, and $\sigma\\{\\}$ is taking the rms over all $i$ iterations.
The rms over all transfer function iterations which include the injected true
observations should provide an error estimate which incorporates thermal
noise, foreground residuals, residual RFI, sample variance and signal loss
from foreground cleaning. blackHowever, crucially this approach relies on
modifying the transfer function definition so that the
$\textbf{{X}}_{\text{clean}}$ term in Equation 10 is not subtracted.

We begin by demonstrating the impact subtracting the cleaned observed data
$\textbf{{X}}_{\text{clean}}$ has on the transfer function. Since we are
investigating error estimation, for this section it is helpful to use data
with noticeable error-bar size, so we therefore add the dominant
$\sigma_{\text{n}}\,{=}\,1\,\text{mK}$ noise to the MD1GPC simulations. Figure
8 shows the performance of the transfer function for the high-noise
simulations, calculated using Equation 11 and Equation 10, both with and
without the $\textbf{{X}}_{\text{clean}}$ subtraction. It is encouraging to
see that the addition of noise is not majorly affecting the performance of the
transfer function. For the case where $\textbf{{X}}_{\text{clean}}$ has been
subtracted (blue results), the accuracy is only mildly affected relative to
the noise-free results in Figure 4. For the noise-inclusive results of Figure
8, we divide by $P_{\text{FGfree}}$ in the $y$-axis which contains the same
noise as the reconstructed power. This is to divide out the fluctuations
caused by the presence of the dominant noise, allowing analysis into the
performance of the reconstruction alone. For the case where
$\textbf{{X}}_{\text{clean}}$ has not been subtracted (orange results), there
is a slight drop in accuracy at small scales. This small bias is absent when
we use a Hi signal without RSD thus it appears to be caused by the addition of
uniform noise in the presence of an anisotropic Hi signal. We also found this
small bias is decreased when we use the 2D transfer function construction
outlined in Section 3.2, although this relied on there being no polarisation
leakage which otherwise causes the variance of the result to blow up, as we
showed. We discuss the performance of the transfer function in the presence of
anisotropic phenomena later, but given this is a small bias and is absent in
the subtracted $\textbf{{X}}_{\text{clean}}$ case, it is not overly important
for this discussion on error covariance estimation.

Figure 8: Impact from subtracting the cleaned observational data
$\textbf{{X}}_{\text{clean}}$ in Equation 10. The shaded regions represent the
rms over the 100 mocks used to construct the transfer function. The rms is
significantly reduced when the cleaned data is subtracted. This test was done
using a foreground clean with $N_{\text{fg}}\,{=}\,8$ on simulations with high
instrumental noise. Here the reconstructed result is divided by the
foreground-free power spectrum $P_{\text{FGfree}}$ which contains identical
instrumental noise.

Figure 8 suggests that subtracting $\textbf{{X}}_{\text{clean}}$ is the
favorable strategy if one is purely pursuing the most accurate transfer
function possible. However, the main point of this plot is the difference in
variance between the two cases. If one is using the variance on the transfer
function as a basis for the error estimation, the reduced variance caused by
subtracting $\textbf{{X}}_{\text{clean}}$ has implications and can lead to
under-estimated errors on the reconstructed power spectrum, as we will now
demonstrate.

To quantitatively evaluate error estimation performance, it is helpful to
analyse how errors for a power spectrum measurement are analytically derived
in a foreground-free case. The variance on a cross-correlation power spectrum
$P_{\text{c}}$ between two tracers, 1 and 2, can be estimated as (Feldman et
al., 1994)

$\sigma_{\text{theory}}^{2}=\frac{1}{2N_{\text{m}}}\left[P_{\text{c}}^{2}(k)+\left(P_{1}+\frac{V\sigma_{1}^{2}}{\langle
f_{1}\rangle^{2}}\right)\left(P_{2}+\frac{V\sigma_{2}^{2}}{\langle
f_{2}\rangle^{2}}\right)\right]\,,$ (16)

where $N_{\text{m}}$ is the number of modes spherically averaged in each
$k$-bin, $V$ is the volume of a single voxel on the Fourier grid, uncorrelated
noise in the field is represented by the variance $\sigma^{2}$, which for an
ideal intensity map would be the variance of the instrumental noise, and
$\langle f\rangle$ is the background mean for the field e.g. the mean
brightness temperature for intensity mapping. For galaxy surveys, the noise
component given by the second terms in the curved brackets will reduce to
shot-noise i.e. $V\sigma^{2}/\langle f\rangle^{2}\,{=}\,1/\bar{n}_{\text{g}}$,
where $\bar{n}_{\text{g}}$ is the galaxy number density. We refer the reader
to Blake (2019) where a detailed derivation of the above is provided with
applications to intensity mapping and its cross-correlations with galaxy
surveys.

Extending Equation 16 to incorporate contributions from foreground
contamination is challenging. Foreground residuals could presumably be
estimated and would provide an additive variance, or be assumed sub-dominant
enough not to warrant inclusion. However, the uncertainty from the transfer
function, which for high signal loss is non-negligible at large scales,
requires careful inclusion. The uncertainty in the transfer function can be
estimated from the variance over the mocks used to construct it as we have
shown, however, analytically adding this into the error budget of Equation 16
is non-trivial since this will not necessarily be a contribution entirely
independent of the noise and cosmic variance already being factored for in
Equation 16. This is why some previous analyses have used the transfer
function variance as a basis for overall error estimation. To evaluate whether
this is a robust method, we can use the analytical errors as a benchmark.
Errors estimated based on the variance in the transfer function should
approximately agree on small scales with the analytical ones where foreground
contamination and signal loss are minimal, but the large noise still
dominates.

We validated that the analytical errors (Equation 16) are a good estimate for
the foreground-free case. Using the MD1GPC simulations with the large
$\sigma_{\text{n}}\,{=}\,1\,\text{mK}$ Gaussian white-noise but with no
foregrounds, we measure the cross-correlation power spectrum with a noise-free
equivalent, then estimate the errors using Equation 16 and evaluate the
$\chi^{2}_{\text{dof}}$ given by

$\chi^{2}_{\text{dof}}=\sum_{k}\frac{P_{\text{data}}(k)-P_{\text{mod}}(k)}{\sigma_{P}(k)}\Big{/}(N_{\text{k}}-1)\,,$
(17)

where $N_{\text{k}}$ is the number of $k$-bins. $P_{\text{mod}}$ is the model
defined in Appendix B.1 with parameters matched to the Multi-Dark inputs used
in the MD1GPC simulation. The analytical errors, return a
$\chi^{2}_{\text{dof}}\,{\sim}\,1$ as expected, evidence that the errors are a
reasonable size, given the reliable model.

Figure 9: Comparison between methods of error estimation for simulations with
signal loss reconstructed by a transfer function. The black dashed-line shows
an analytical error estimate given by Equation 16 which we have shown to be
reliable for foreground-free data. The coloured lines show error estimation
based on the variance of the transfer function, under two difference
scenarios, with $\textbf{{X}}_{\text{clean}}$ subtraction (see Equation 10)
and without. This is with the MD1GPC simulation where $N_{\text{fg}}\,{=}\,8$
PCA modes removed.

Figure 10: Covariance and correlation matrices for foreground-free (left
column) and foreground cleaned with transfer function reconstructed signal
loss (right column), without the $\textbf{{X}}_{\text{clean}}$ subtraction in
the transfer function. Produced using the MD1GPC simulations with
$N_{\text{fg}}\,{=}\,8$ PCA modes removed, as in Figure 10. The covariance
matrices have been multiplied by the product of the two-$k$-bins cubed times
$10^{3}$ for demonstration purposes so high-$k$ covariance can be seen.

Using the analytical errors $\sigma_{\text{theory}}$ as a validated benchmark,
Figure 10 shows how the error estimation based on the variance in the transfer
function compares for the same simulations but with cleaned foregrounds. The
blue line shows the case using a transfer function defined by Equation 11 and
Equation 10, where the cleaned observed data ($\textbf{{X}}_{\text{clean}}$)
has been subtracted to reduce the variance. This is underestimating the errors
relative to the analytical ones, likely caused by subtracting the
$\textbf{{X}}_{\text{clean}}$ term which will contain noise and foreground
residuals, which contribute to the error budget. Figure 10 therefore shows
that subtracting the cleaned data $\textbf{{X}}_{\text{clean}}$ in the
transfer function construction results in the transfer function variance no
longer being a reliable means for estimating the errors. However, it is clear
from the orange line, which is equivalent to the blue but without the
$\textbf{{X}}_{\text{clean}}$ subtraction, that the variance from this version
of the transfer function leads to an excellent agreement with the analytical
errors at high-$k$ as required. Furthermore, it also shows an increase in
error at small-$k$ as one would expect where signal loss is highest, thus it
is incorporating the increased uncertainty from signal reconstruction, not
accounted for in the analytical errors.

Using the variance in the transfer function mocks for analysing uncertainties
also has the advantage of being able to examine off-diagonal covariance of the
data, something not trivially possible with the analytical approach. In Figure
10 we show the $k_{i}k_{j}$ covariance matrix $C_{ij}$ (top row) as well as
the normalised correlation matrix defined by
$R_{ij}\,{=}\,C_{ij}/\sqrt{C_{ii}C_{jj}}$. The left column shows the
foreground-free scenario where we inject mocks into the Hi \+ high-noise
simulations to get an estimate of the covariance without a foreground cleaning
step. The right column is equivalent but with cleaned foregrounds and a
transfer function constructed without the subtraction of
$\textbf{{X}}_{\text{clean}}$. The covariance over all iterations in the
reconstructed power spectra (Equation 15, but now including off-diagonal
elements $i\,{\neq}\,j$) estimates the covariance of the observed data. As
expected, and consistent with the orange line of Figure 10, the cleaned
foregrounds are increasing covariance on large scales, but encouragingly they
do not appear to increase off-diagonal correlations between $k$-bins.

We conclude from this investigation that the variance in the transfer function
is a reliable tool for error estimation in the final reconstructed power
spectrum, provided the cleaned data $\textbf{{X}}_{\text{clean}}$ term is not
subtracted in its construction. This becomes similar in approach to galaxy
surveys which use vast suites of mocks as their primary method for estimating
the covariance in their data (e.g. Zhao et al., 2021). If opting to use the
transfer function variance for error estimation, it becomes important to
ensure that all aspects of the error budget are emulated in the mocks used in
the transfer function construction. An example of this would be in a galaxy
survey cross-correlation where the galaxy shot noise would need to be captured
by using galaxy mocks with the correct number densities and survey coverage.
With the observational data injected into the mocks, we are also including
variance from signal loss, foreground residuals, residual RFI, instrumental
noise, etc. all of which are currently not well enough understood to reliably
emulate in mock intensity maps555Jackknife resampling will also be a useful
tool when unknown systematics are present..

###  3.4 Auto- and cross-correlation applications

So far in this section, we have considered the cross-power spectrum between a
foreground cleaned Hi intensity map and the original foreground-free, Hi-only
map. This is so that any foreground residuals or noise in the cleaned maps do
not complicate the analysis of Hi signal loss in the power spectra. Since
foreground residuals and noise will not correlate with the Hi-only maps, the
additive biases they cause are avoided in cross-correlation thus the only
departure from the true-Hi power should be just from signal loss. However, Hi
intensity mapping surveys will also aim to conduct analysis in auto-
correlation and we need to consider how signal loss behaves in this scenario.

It has been previously suggested that there would be twice as much signal lost
in the auto-correlation power spectrum because its effects are present twice
in the map product
$P_{\textrm{H{i}}}(\boldsymbol{k})\,{\propto}\,|\tilde{\textbf{{X}}}(\boldsymbol{k})|^{2}$.
This would mean a $\mathcal{T}(k)^{-2}$ correction factor is needed to
reconstruct the power spectrum. However, we found from our simulations that
this is not the case, and the same degree of signal loss is also present in an
auto-correlation as is in cross-correlation. In other words, the same
correction of $\mathcal{T}(k)^{-1}$ is also needed in the auto-correlation as
well as in cross-correlation. Figure 11 demonstrates this finding showing how
the cross-correlation (blue-dotted line) and auto-correlation (orange-solid)
appear to have approximately equivalent levels of signal loss. For this test,
we return to the noise-free simulations and use an aggressive
$N_{\text{fg}}\,{=}\,12$ PCA clean which will suppress foreground residuals
significantly, making it reasonable to ignore their influence on the results.
The green line shows what appears to be the correct application of the
$\mathcal{T}(k)^{-1}$ transfer function whereas the red line shows the
consequential over-correction from applying the $\mathcal{T}(k)^{-2}$ to the
auto-correlation.

Figure 11: Demonstrating the correct application of the transfer function in
an auto-correlation analysis. Black-dashed line shows the original foreground
free power spectrum. The blue-dotted line shows the cross-correlation with the
foreground-free. The orange line shows the auto-correlation. For this
aggressive ($N_{\text{fg}}\,{=}\,12$) foreground clean, foreground residuals
should be minimal and the similar amplitude between $P_{\text{auto}}$ and
$P_{\text{cross}}$ suggests the signal loss is similar in both. The green line
shows the correct application of the transfer function and the red line shows
the over-correction where $\mathcal{T}^{-2}$ is used.

To provide a deeper understanding for why signal loss to the power spectrum is
the same for auto- and cross-correlations, we present in Figure 12 the
amplitude of all 3D-Fourier mode products for a randomly chosen $k$-bin. The
spherically averaged $P(k)$ value for the chosen
$0.0590\,{<}\,k\,{<}\,0.0687\,\,h\,\text{Mpc}^{-1}$ bin is then simply the
average of all these amplitudes, which is stated in the legend for each
scenario. The top-panel shows the comparison between an auto-correlation with
foreground cleaning and the cross-correlation between foreground cleaned and
foreground free (Hi-only) data. As can be seen, the average of the modes is
approximately the same in both cases and thus consistent with Figure 11,
demonstrating that signal loss is equivalent in auto- and cross-correlation.
The reason for this is related to the fact that the same modes are projected
out of the analysis in both cases and signal loss does not compound when two
maps with the same removed modes are combined in an auto-correlation. We
confirm this to be the case in the bottom panel of Figure 12 where we blackuse
simulations with a foreground from a different region of sky so that we
produce a cleaned map $\textbf{{X}}^{\text{clean}}_{\text{FG2}}$ which will
have different foreground modes removed compared to the original
$\textbf{{X}}^{\text{clean}}_{\text{FG1}}$ used in the rest of the paper for
the MD1GPC simulation. The exact regions are not overly important just the
fact that they will generate a different set of modes which are projected out
in the foreground clean. When these two foreground cleaned maps are cross-
correlated (red results) we now get a drop in power relative to the cross-
correlation between $\textbf{{X}}^{\text{clean}}_{\text{FG1}}$ and the Hi-only
map (see the mean power in the legend) showing that the signal loss is related
to the modes being projected out in the foreground clean, and it is only a
difference in these which will create further signal loss in a Hi auto-
correlation.

The results demonstrated by Figure 12 have consequences for auto-correlation
analyses with Hi intensity mapping. Not just because it further confirms that
signal loss should be the same in auto- and cross-correlation where the same
foreground modes are projected out, but also because the results in the bottom
panel show when two differently cleaned maps are cross-correlated, the signal
loss becomes more complex to estimate. This is relatable to a method that is
likely to be pursued when attempting an auto-correlation detection whereby
cross-correlations are measured between different sub-sets of the
observations, created either by splitting data into different time-blocks
("sub-seasons") as done in GBT analysis (Masui et al., 2013; Wolz et al.,
2022), or by splitting different dishes as is possible with a multi-dish
telescope such as MeerKAT. This is pursued in order to avoid the additive
biases from noise and time- or dish-dependent systematics. Whilst this method
would still observe the same foreground, the response to systematics may be
different in each sub-set creating a scenario similar to the red results in
Figure 12. We leave further investigation into this specific form of auto-
correlation method for future dedicated studies.

Figure 12: Amplitude of all Fourier modes
$P(\boldsymbol{k})\,{\propto}\,\operatorname{Re}\left\\{\tilde{\textbf{{X}}}_{\text{A}}(\boldsymbol{k}){\cdot}\tilde{\textbf{{X}}}_{\text{B}}^{*}(\boldsymbol{k})\right\\}$
in the range $0.0590\,{<}\,k\,{<}\,0.0687\,\,h\,\text{Mpc}^{-1}$, which
corresponds to one chosen $k$-bin in the spherically averaged power spectrum.
The average of these modes, stated in the legend for each scenario, will
represent the spherically averaged power spectrum value for the particular
$k$-bin. Top-panel shows the equivalence in signal-loss between auto- and
cross-correlation. Lower panel shows how signal loss can be larger where a
different set of modes has been projected out in the foreground clean shown by
the red results which is the cross-correlation between two different simulated
foreground regions.

##  4 Applications to pathfinder intensity mapping with MeerKLASS

In the previous sections, the MD1GPC simulations used have been deliberately
kept free of further observational effects (except for a couple of identified
cases) besides the foreground contamination which included simulated
polarisation leakage. However, a relevant question for current pathfinder
single-dish experiments is whether the early pilot surveys, which typically
have low signal-to-noise and additional systematic observational effects, can
also rely on the foreground transfer function to correct for signal loss. In
these pathfinder surveys (e.g. Wolz et al., 2022; Cunnington et al., 2022)
foreground cleaning is typically aggressive and signal loss can reach high
levels, thus one could argue that we are more reliant on signal reconstruction
in these early surveys, compared with future surveys where foreground cleaning
and systematics will be more controlled.

The cosmological detection in Cunnington et al. (2022) (like all other
intensity mapping detections preceding it) relied on a foreground transfer
function to reconstruct the signal loss from foreground cleaning. The
$7.7\sigma$ cross-correlation detection significance fell to ${\sim}\,4\sigma$
without signal reconstruction. The survey covered just $200\,\text{deg}^{2}$
and only the frequency channels spanning $1015-973\,\text{MHz}$ (0.4 < z <
0.46) were used to avoid the worst RFI. Furthermore, the observation gathered
just $10.5\,\text{hours}$ of data per dish. This means sky coverage and
signal-to-noise were low and foregrounds could be easily impacted by
systematics, rendering signal loss more complex and widespread than the
examples we have investigated so far.

The reliability of the transfer function was validated for the results in
Cunnington et al. (2022) and here we demonstrate these validation tests by
utilising a different set of simulations which we refer to as the MDMK
simulations, which aim to emulate the MeerKAT 2019 pilot intensity mapping
survey (Wang et al., 2021). We use the survey’s non-uniform mask and use
fluctuations in the uncleaned foreground sky to generate systematic
perturbations in the MDMK foreground simulations creating a demanding cleaning
requirement. We outline the details of how this is achieved in the following
section.

###  4.1 MDMK MeerKLASS simulations

To emulate current MeerKAT pathfinder data and investigate the performance of
the transfer function when signal loss is spread across a wide range of
scales, we utilise the MeerKAT 2019 pilot survey data (Wang et al., 2021;
Cunnington et al., 2022). The survey targeted a single patch of
${\sim}\,200\,\text{deg}^{2}$ in the WiggleZ 11hr field, covering
$153^{\circ}\,{<}\,\text{R.A.}\,{<}\,172^{\circ}$ and
$-1^{\circ}\,{<}\,\text{Dec.}\,{<}\,8^{\circ}$. The telescope observed at
constant elevation, scanning back and forth through azimuth taking
$1.5\,\text{hours}$ to complete one time-block. 7 time-blocks were obtained
with a mix of rising and setting scans, creating offset coverage providing the
footprint which can be seen in the MDMK simulation maps in Figure 14. Holes in
the footprint are evident and are caused from the multi-stage RFI flagging
which can leave gaps in the scanning. For the MeerKAT Hi simulation, (top-
panel of Figure 14), we use the same Multi-Dark simulation as in MK1GPC
(Section A.1) but calculate the physical volume covered by the MeerKAT 2019
data and cut a volume of this size from the $1\,(\text{Gpc}/h)^{3}$ cube. We
use a similar pixelisation as the 2019 data
($n_{\text{x}},n_{\text{y}},n_{\text{z}}\,{=}\,133,41,250$), and then apply
the exact same footprint mask. The MeerKAT 2019 observations were performed in
L-band but to avoid dominant RFI, only 199 channels with a
$973.2-1014.6\,\text{MHz}$ frequency range ($0.400\,{<}\,z\,{<}\,0.459$) were
used. Additionally, of the 199 channels selected to use, a further 32 are
removed due to evidence of RFI contribution in their eigenmodes. We also
replicated this exact channel flagging in the MDMK simulations.

####  4.1.1 Frequency perturbed foregrounds

Evidence of systematics was seen in the MeerKAT 2019 data, which was expected
given the low amount of observational time and it being a first of its kind
pilot survey. One way systematics were evident was in the perturbations to
what should be smooth spectra in the raw foreground sky. The exact cause of
these systematic perturbations is beyond the aim of this paper but we can
still use the distorted spectra from the real data to perturb an idealised
foreground simulation, emulating the main impact from these systematics on the
foreground clean and signal reconstruction.

Figure 13: MDMK simulated maps aiming to emulate MeerKAT 2019 intensity
mapping data, averaged along the $973.2{-}1014.6\,\text{MHz}$ frequency range.
Top panel shows HI only, produced using the Multi-Dark $N$-body semi-
analytical simulation. Middle-panel shows the foregrounds simulated using a
perturbed version from the Planck Sky Model. The points on the foreground map
indicate the positions of the example spectra plotted in Figure 14, with green
and red points representing pixels near the centre or edge of the MeerKAT
footprint respectively. The bottom panel shows an estimate for some MeerKAT
noise and residual systematics obtained by subtracting data observed at
different times (see Section 4.1.2).

Figure 14: Example spectra from the MDMK foreground simulations aiming to
emulate MeerKAT pilot survey data. The perturbations to the spectra have been
produced using the MeerKAT 2019 data to create realistic foreground
simulations with systematic effects (see Section 4.1.1). The green-dashed
lines represent pixels taken from a central position in the MeerKAT footprint
(corresponding to green points in Figure 14). The red-solid lines represent
pixels near from the edge of the footprint (corresponding to red points in
Figure 14) and are more vulnerable to systematics, hence the more noticeable
perturbations. The grey-shaded regions represent channels which were flagged
in the cross-correlation analysis (Cunnington et al., 2022), which we also
flag in this simulation for consistency.

To create foregrounds for the MDMK simulations we begin by using the Planck
Sky Model to generate synchrotron and free-free emission at the relevant
frequencies and sky position, as with the MD1GPC simulation. The bottom panel
of Figure 14 shows the foreground simulation for the MeerKAT 2019 footprint.
Unlike the MD1GPC simulation, we do not include any simulated polarisation
leakage and instead use the MeerKAT 2019 data itself, aiming to create more
realistic systematic perturbations to the foregrounds. This is done by fitting
a smooth polynomial to each line-of-sight in the 2019 data, then the
systematic perturbations to the foreground spectra can be approximated by the
ratio between the data and the polynomial i.e.

$T_{\text{perturbation}}(\boldsymbol{x},\nu)=\frac{T_{\text{2019-data}}(\boldsymbol{x},\nu)}{T_{\text{smooth-
poly}}(\boldsymbol{x},\nu)}\,.$ (18)

We found on average these perturbations were small sub-percent values however
they could be as high as $3.7\%$. We multiply these perturbations with the PSM
foreground and Figure 14 shows some example perturbed spectra from the final
simulation. The perturbations are worse near the edge of the map (shown as
red-solid lines) where due to the lower coverage, systematics can have more
impact. This is consistent with what was found in the actual data and in the
cross-correlation analysis these edge pixels are down-weighted (Cunnington et
al., 2022). Figure 14 also shows the flagged channels (grey regions) used in
the cross-correlation analysis which we also adopt in the MDMK simulation.

####  4.1.2 Anisotropic systematics and RFI residuals

The analysis of the auto-correlation power spectrum for the 2019 MeerKAT
survey in Cunnington et al. (2022) showed evidence of additive biases most
likely from instrumental noise, residual RFI, or other systematics. Their
contribution appears to dominate over the Hi signal because the auto-power
spectra amplitude was larger than one would expect from Hi only power. We also
include a contribution to the MDMK simulation which attempts to emulate these
types of additional components. Again, we use the real MeerKAT 2019 survey
itself to produce a map of time-varying anisotropic contributions and add
these onto the MDMK Hi and perturbed foreground maps. This is achieved by
taking the residuals from different time blocks in the MeerKAT 2019 survey. We
take the difference between the first four time-blocks and the last three
where these residuals will represent components that vary in time. Therefore,
in principle, this should not include the Hi signal or the foregrounds since
these would be consistent in time, but instead only include time-varying
systematic contributions, which is what we are aiming to emulate.

The map of the MeerKAT time-block residuals is shown in the bottom panel of
Figure 14. In some instances, there are no shared pixels in both time-block
groups so the residual is undefined. For these pixels we instead resort to
adding a large level of Gaussian random noise whose variance dominates the Hi
signal by one order of magnitude. However, when these are plotted in Figure
14, which is the average along the line-of-sight, their contribution is
averaged down which is why the amplitude appears relatively low around the
edges where most of the missing pixels between time-blocks are. The pixels
from the actual MeerKAT residuals, which are most concentrated in the centre
where shared coverage is better, appear higher relative to the Gaussian noise.
This will be due to the residuals being more correlated in frequency, thus do
not average down in the plot. The frequency correlation and apparent
anisotropies of the residuals as evident in Figure 14, suggest they are
contributions beyond simple instrumental noise. Whilst this is a complication
for the pilot survey analysis, it is useful for our purposes, providing
additional complications to foreground cleaning and signal loss in the MDMK
simulations.

In these simulations, whilst we do not explicitly include the effects from a
realistic telescope beam, some of the impacts it has on the foregrounds will
be included in the perturbations we add to the simulations from Equation 18. A
simple Gaussian beam is trivial to include and assuming it is approximately
matched in the transfer function construction, it makes no difference to the
performance of the transfer function as we will explicitly show in the
following section Section 5.1. However, in reality, the MeerKAT beam will be
more complex with wide-reaching frequency-dependent side lobes which could
complicate foreground cleaning (Matshawule et al., 2021; Spinelli et al.,
2021). Trying to replicate this in the mocks in the transfer function
construction may be difficult and an investigation into what level this needs
to be considered should be pursued. This requires a more detailed simulation
which we will pursue in follow-up work.

###  4.2 Correcting signal loss in pathfinder data

Figure 15 shows power spectra for the MDMK (MeerKAT-based) simulations
presented in the previous sub-section. The black dashed line shows the
foreground-free Hi-only result, and the red solid line shows the result from
adding the perturbed foregrounds and residual time-varying systematics based
on real MeerKAT data, then performing a $N_{\text{fg}}\,{=}\,10$ PCA clean. We
find that using the perturbed foregrounds (described in Section 4.1.1) is the
main cause for requiring an aggressive foreground clean, highlighting the
importance of instrument calibration so that smooth spectra are maintained in
the observed data. Similar to MeerKAT data (Cunnington et al., 2022), signal
loss in the foreground cleaned data is widespread throughout all scales, with
noticeable signal loss occurring even in the highest-$k$. The main reason for
this is due to the decreased depth of the frequency/redshift range. We tested
this with the main MD1GPC simulation. Reducing the number of pixels along the
LoS by a factor of 4 to 64 pixels with $0.25\,\text{Gpc}/h$ of depth, produced
${>}\,75\%$ signal loss in the smallest 4 $k$-bins, and still $13\%$ in the
highest $k$-bins. This was even without the polarisation leakage and just
removing 4 PCA modes. Whereas for an equivalent scenario but using the full
$1\,\text{Gpc}/h$ depth, the signal loss is never greater than $30\%$ and only
$3\%$ at the highest-$k$. This can be understood by considering that modes
projected out of narrow frequency-range data will be confined to a
higher-$k_{\parallel}$ space. Thus the signal loss will also spread into
higher $k$-modes. Encouragingly this means that signal loss should be
naturally mitigated in future observations by using a larger frequency range.
This will be possible with MeerKAT UHF-band observations which will probe
lower frequencies where RFI is expected to be less dominant, thus a more
complete frequency range can be used.

Figure 15: Power spectra for MDMK simulations which are small in volume, have
a more complex non-continuous survey footprint, high anisotropic noise and
residual systematics, as well as systematically perturbed foregrounds, all to
emulate actual MeerKAT pathfinder data. We show the cross-correlation with the
Hi-only maps. Signal loss from the $N_{\text{fg}}\,{=}\,10$ PCA clean (red
line) is larger and more widespread into high-$k$ compared with previous
results in the more idealised MD1GPC simulation. Despite this the transfer
function is encouragingly still able to reconstruct a reasonably unbiased
result shown by the blue data points. Error bars are given by the limits of
the central 68th percentile region from the distribution of reconstructed
power spectra using the transfer function mocks.

Despite the more complex and widespread signal loss in Figure 15 (red line),
when we construct a transfer function using the process summarised in Section
3.1, we are able to reconstruct the correct Hi power spectrum. As in previous
tests, the clean and corrected power spectra are cross-correlations with the
original-Hi to avoid any issues with residual foreground contamination
confusing the assessment of signal loss. The power spectra are naturally more
noisy than the previous MD1GPC simulations due to the decreased volume, the
systematically perturbed foregrounds (see Figure 14) and the large time-
varying systematics (Figure 14 bottom panel) inserted into the simulation. The
instrumental noise and additive systematics is an important consideration that
we did not include in the default MD1GPC simulations of previous sections.
This additional noise will introduce extra perturbations to the foreground
modes, in the same way the signal introduces perturbations (see Figure 1). If
the noise is large, as is the case for pathfinder observations with low
observational time, these perturbations will be large. Encouragingly, this
does not appear to cause noticeable problems for the transfer function,
evidenced by the corrected result in Figure 15, which includes large
additional contributions that dominate over the Hi signal.

Given the more complex nature of the pilot survey simulation, some
reconstructed power spectra from the transfer function mocks produced outliers
and returned non-Gaussian distributions for each $k$-mode. In this scenario,
using the rms over the mocks for the errors would be a poor estimation and
would be overly distorted by the outliers. We therefore instead use the 68th
percentiles limits to provide the asymmetric error bars, which is what are
presented in Figure 15. To obtain the converged distribution, we used 1000
mocks in the transfer function calculation. This presents a further advantage
of using the transfer function mocks for error estimation, providing more
options to handle non-Gaussian uncertainties. This of course would have
complex implications for further analysis and parameter estimation, which we
do not investigate here. However, errors should naturally become more Gaussian
as intensity map quality improves and noise, systematics, etc. are reduced.

Figure 16: Foreground transfer function in cylindrical
$k_{\perp},k_{\parallel}$ space for the MDMK simulations emulating MeerKAT
pathfinder data. $N_{\text{fg}}\,{=}\,10$ PCA modes have been removed for the
foreground clean. The 1D spherically averaged version of this is used for the
corrected results in Figure 15.

Figure 16 shows the computed transfer function in cylindrical
$k_{\perp},k_{\parallel}$ space for the MDMK simulations. It is interesting to
analyse the differences between this more realistic case and that from the
more idealised MD1GPC simulations in Figure 6. This again reveals that signal
loss is more widespread into larger $k_{\parallel}$ modes, which is consistent
with the more widespread signal loss evident in Figure 15. We still see large
regions where $\mathcal{T}\,{\sim}\,1$ suggesting that the approach of
discarding some regions of $k$-space to massively reduce the dependence on the
transfer function (as we discussed in Section 3) could still be pursued even
with small intensity mapping pilot surveys.

The results from the MDMK simulations provide validation for using the
foreground transfer function with pathfinder survey intensity maps. We have
used MeerKAT data itself to attempt to complicate the foreground clean and
signal loss to stress-test the current signal reconstruction process. The
simulations produced make no attempt to understand the source of the
perturbations to foreground spectra or additive time-varying systematics,
simply emulating them as realistically as possible to mimic the challenge they
pose. Whilst the success is encouraging, the investigation would be completed
further by including specific simulations of known observational effects such
as directly simulating RFI, 1/$f$ correlated noise, non-uniform noise, a non-
Gaussian beam etc. This would allow more analysis into exactly what
observational effects are most troublesome for foreground cleaning and signal
loss. The development of a robust simulation pipeline for MeerKAT single-dish
intensity maps including a realistic beam and all these observational effects
is outside the aims of this paper but is being pursued in other MeerKLASS
collaboration projects (e.g. Irfan et al., 2023).

##  5 Precision cosmology suitability

In this section we look to future Hi intensity mapping observations and test
how reliable a transfer function would be where sub-percent accuracy on
parameter estimates is required. We return to the more generic simulations of
MD1GPC to avoid the investigation being confused by the large statistical
noise present in the previous section’s realistic simulations of a MeerKAT
pilot survey.

###  5.1 Mock parameter dependence

Up until now, it has not been investigated how robust the accuracy of the
transfer function is when there are discrepancies between parameters used in
the transfer function construction and their true fiducial values in the real
observed data. If the parameter assumptions used for the generation of 100
mocks strongly influence the final accuracy of the reconstructed power
spectrum then this is a large concern for precision cosmology since this would
lead to biased cosmological parameter estimates.

In this section, we demonstrate how mild the mock parameter dependency is and
show how large $+/-\,100\%$ discrepancies between the assumed parameters in
the mocks used for transfer construction and the underlying truth in the data,
mostly only yield small ${\lesssim}\,1\%$ inaccuracies in the recovered
parameter estimations. We test this by treating the MD1GPC as the observed
data with the underlying "true" parameters, then vary some of the values
$\\{\Omega_{\textrm{H{i}}}\,{\propto}\,\overline{T}_{\textrm{H{i}}},f,\sigma_{\text{v}},R_{\text{beam}}\\}$
in the lognormal mocks which are used to construct the transfer functions. The
model power spectrum we use to generate the lognormal mocks is described in
more detail in Appendix B.1 but we repeat it here for convenience

$P_{\text{mod}}(k,\mu)=\overline{T}^{2}_{\textrm{H{i}}}\frac{\left(b^{2}_{\textrm{H{i}}}+f\mu^{2}\right)^{2}}{1+(k\mu\sigma_{\text{v}}/H_{0})^{2}}P_{\text{m}}(k)\\\
\times\,\exp\left[-(1-\mu^{2})k^{2}R_{\text{beam}}^{2}\right]\,.$ (19)

For this test we used the MD1GPC simulations with a
$\sigma_{\text{n}}\,{=}\,1\,\text{mK}$ dominant white noise and a default
Gaussian beam where $R_{\text{beam}}\,{=}\,10\,h^{-1}\text{Mpc}$.

Figure 17 shows how subtle the impact on parameter estimation is when
parameters used in the mock’s model power spectra, given by the panel titles,
are biased by an amount indicated by the $x$-axis. The $y$-axis shows the
percentage bias relative to the foreground-free parameter estimation. In all
cases, we sample the parameter posterior distribution in a Bayesian MCMC only
varying one parameter at a time, fixing all other parameters in the model
(Equation 19) to fiducial values fitted to the foreground-free MD1GPC
simulation. The error bars represent the 68% confidence regions in the
posterior distributions and they are plotted relative to the median of the
foreground-free posterior. To avoid non-linear complications in the modeling
we only use modes where $k\,{<}\,0.3\,\,h\,\text{Mpc}^{-1}$. Whilst these
scales are still quite non-linear, our model worked reasonably well on these
scales and since we are testing its performance relative to the foreground-
free case, any shortcomings caused by non-linear effects will be present in
both.

Figure 17: Robustness of the transfer function in response to biased parameter
assumptions. Each panel shows a different parameter used in the mocks which
construct the transfer function. The $x$-axis shows different %-biases
relative to the correct fiducial values. The $y$-axis shows the estimated
parameter posterior from an MCMC on the reconstructed data using a transfer
function. The error bars represent the 68% confidence region in the estimated
posteriors. These are plotted relative to the foreground-free parameter
estimates to demonstrate that only small biases are induced from incorrect
parameter assumptions in the transfer function construction. These results are
produced with the MD1GPC simulation with
$\sigma_{\text{n}}\,{=}\,1\,\text{mK}$ dominant noise and a Gaussian beam with
$R_{\text{beam}}\,{=}\,10\,h^{-1}\text{Mpc}$. All results are for a
$N_{\text{fg}}\,{=}\,12$ PCA clean.

$\Omega_{\textrm{H{i}}}$, which is proportional to
$\overline{T}_{\textrm{H{i}}}$ (see Equation 31), will only change the
amplitude of the Hi power spectrum in our simple linear model, and the
transfer function appears extremely robust to these scale-independent
amplitude changes, with no ${>}\,1\%$ bias being induced. We even tested going
to $4\times$ the fiducial truth on $\Omega_{\textrm{H{i}}}$, since this is a
very unconstrained parameter, but this still yielded a sub-percent bias. Note
that the red-cross indicates an undefined result for the
$\Omega_{\textrm{H{i}}}\,{-}\,100\%$ case since the intensity maps are zero
when
$\overline{T}_{\textrm{H{i}}}\,{\propto}\,\Omega_{\textrm{H{i}}}\,{=}\,0$. The
increase in uncertainties with a decreasing $\Omega_{\textrm{H{i}}}$ mock
parameter input is caused by the white noise having a more dominant impact
relative to these lower amplitude mocks in the transfer function construction.

The growth rate $f$ is included to model the redshift-space distortions (RSD)
which are present in the MD1GPC simulations. $f$ will be an anisotropic
parameter since RSD are a line-of-sight-only effect. This appeared to induce
more of a bias in the reconstructed power spectra, but it still only caused
${<}\,2\%$ in most cases tested. For certain cosmological parameters such as
the growth rate, much tighter priors will be applicable such that the biased
parameter values we have used for $f$ will be unrealistic. It is interesting
to see that failing to include a linear RSD model in the mocks at all, as
shown by the ${-}100\%$ result (i.e. $f\,{=}\,0$), induces a ${\sim}\,5\%$
bias which suggests that it is important to include some basic anisotropic RSD
in the mocks for reconstruction accuracy.

$\sigma_{\text{v}}$ is again used to model RSD but as a phenomenological
attempt to model fingers-of-god (FoG) on mildly non-linear scales. This will
also be an anisotropic parameter but is also directly scale-dependent too,
having a greater influence at high-$k$. This shows ${\sim}\,3\%$ bias if the
parameter is set to zero. Even though this is a phenomenological parameter
without a physically defined fiducial value, it is still unlikely that the
parameter will be completely omitted without a suitable method for replacing
its modelling effects. We highlight that the increase in error bars for the
$f$ and $\sigma_{\text{v}}$ parameters is not necessarily foreground-related.
Some parameters will be constrained better than others and $\sigma_{\text{v}}$
relies more on small scales which are damped by the beam. Furthermore, both
parameters would be more suitably constrained by including the quadrupole as
opposed to the spherically averaged monopole we are using here.

Lastly, we introduce the size of the beam $R_{\text{beam}}$ as a varying
nuisance parameter. This is defined as the rms of the beam profile in comoving
units at the probed redshift. Similarly to $\sigma_{\text{v}}$ this is another
scale-dependent anisotropic parameter, although this time affecting
high-$k_{\perp}$ modes. This can reach a ${\sim}\,5\%$ negative bias if
completely unaccounted for, shown by the ${-}100\%$ result. For nuisance
parameters linked to the instrument such as $R_{\text{beam}}$, we should have
a much tighter prior on its value effectively ruling out a ${>}20\%$ incorrect
assumption. The small ${\sim}\,1\%$ bias relative to the truth in the case
where the correct fiducial beam has been used (see $0\%$ input mock parameter
bias result), suggests there could be some discrepancy in high-$k_{\perp}$
modes between foreground-free and reconstructed foreground-cleaned data, where
$R_{\text{beam}}$ has the most impact, leading to inaccuracies in its
estimation. Given this is only a very small bias and is only in a nuisance
parameter, we defer this to future work, where we will investigate the impact
on signal reconstruction in the presence of more complex beams.

Given that the small biases appear to be caused by anisotropic divergences
between mocks and observations, we investigated whether the 2D transfer
function (discussed in Figure 7) could yield improvements. We found that when
there is no polarisation leakage in the simulations, the average bias on the
reconstructed power spectrum from using $f\,{=}\,0$ in the mocks is $3.8\%$
relative to the foreground-free power spectrum. Interestingly though, when we
construct and apply the transfer function in 2D, then rebin into 1D following
the same procedure in Equation 14, the bias is only $1.9\%$. However, as we
found from the results in Figure 7, the variance blows up when we reintroduce
the more complex foreground with polarisation leakage. We thus leave further
investigation to future work with more realistic simulations where we will
definitely test if the large variance in the 2D transfer function can be
reliably reduced, and if it then still decreases the small biases from
anisotropic inconsistencies we see in Figure 17.

There are some important conclusions to draw from the results in Figure 17.
Firstly, the results demonstrate how the transfer function works. It is not a
process where an exact replica of the real data is required to measure the
precise impact of signal loss. Rather the mocks injected act as a test field
to construct a response function caused by the foreground cleaning. Broadly
speaking, it appears that the parameters used in the mocks for the transfer
function do not have a strong influence on the final accuracy of the
reconstructed power spectrum. However, where sub-percent accuracy is the aim,
it is clearly beneficial to have the mocks attempt emulation of some of the
features inherent in the observational data. For example, completely
neglecting the telescope beam or linear RSD in the mocks, can have a
noticeable impact on the reconstruction accuracy. This lays the foundation for
many further inquiries into this topic. For example; will a more realistic
frequency-dependent beam with a side-lobe structure be sufficiently emulated
by a Gaussian beam in the mocks for the purposes of the foreground transfer
function? Do any of the other multitude of parameters that we want to probe or
are forced to include in our model as nuisance parameters, have a stronger
influence on reconstruction accuracy? What happens when we allow multiple
parameters to vary simultaneously as apposed to varying one parameter at a
time as done for the results in Figure 17? Whilst these results are strong
evidence that a transfer function will not bias parameter inference, the most
robust way to confirm this would be with more realistic simulated
observational effects e.g. a MeerKAT model of the beam, and to perform
detailed modelling with a multi-parameter MCMC fit to the reconstructed power
spectrum, including shape parameteres e.g. $h$, $n_{\text{s}}$,
$\omega_{\text{c}}$, $\omega_{\text{b}}$, under a range of different
scenarios. This is beyond the scope of this work but is something we will aim
to showcase in a follow-up study.

Since it seems the fiducial mock parameters have a sub-dominant influence over
the reconstruction accuracy, one prospect to consider is a process whereby the
mock parameters are updated based on the parameter inference from the real
data. In this way an iterative transfer function could be developed whereby as
the parameter posteriors for the observed data are estimated and converge on a
final parameter estimate, these values can be used to update the transfer
function calculation and avoid any possibility of it biasing the parameter
inference. We discuss some of the early investigations for this possibility in
Appendix D, but also largely leave this to further dedicated investigation.

###  5.2 Probing exotic physics on ultra-large scales

As a final test of the transfer function’s performance, we examine its ability
to reconstruct signal on the largest scales, even when the underlying true
signal has some unknown "non-standard" properties. For this test, we focus on
primordial non-Gaussianity (PNG) which can be probed on the largest scales in
galaxy surveys (Mueller et al., 2022) and soon in Hi intensity maps (Li & Ma,
2017; Witzemann et al., 2019; Karagiannis et al., 2021).

The nature of the fluctuations in the primordial Universe which arise during
inflation carry a wealth of information regarding the physical mechanisms that
shaped the early Universe. The parameter $f_{\text{NL}}$ quantifies the
departure from Gaussianity in the primordial Universe (Komatsu & Spergel,
2001) and for the so-called local-type of PNG, $f_{\text{NL}}\,{\neq}\,0$
would be evidence of non-Gaussian fluctuations, ruling out slow roll, single-
field inflation in favour of more exotic multi-field models (Creminelli &
Zaldarriaga, 2004). Constraints on PNG so far come from CMB anisotropies and
results are consistent with Gaussian fluctuations with
$f_{\text{NL}}\,{=}\,0.9\,{\pm}\,5.1$ (Planck Collaboration, 2020). However,
large-scale structure surveys, in particular intensity mapping, are expected
to soon lead the way in improving PNG precision. Evidence for PNG in large-
scale structure surveys will manifest as a scale-dependent correction to the
linear bias (Dalal et al., 2008). This correction scales as $k^{-2}$ thus it
is at ultra-large scales where sensitivity to $f_{\text{NL}}$ becomes most
prominent.

In this section, we generate a new underlying Hi intensity map simulation, no
longer using the MD1GPC simulation. The reason for this is first because we
wish to add a clear signature of PNG into the field, and secondly, because we
need to cover much larger scales where we will be able to probe the scales
that are sensitive to the $f_{\text{NL}}$ parameter. We use the $N$-body
COmoving Lagrangian Acceleration (COLA)666pypi.org/project/pycola3 (Tassev et
al., 2013; Tassev et al., 2015), code to generate a fast $N$-body simulation
on a $(8{,}000\,h^{-1}\text{Mpc})^{3}$ grid with $256^{3}$ pixels which
approximately represents a wide and deep Hi intensity mapping survey with
something like SKAO. We seed the simulation with a Hi power spectrum given by

$\displaystyle P_{\textrm{H{i}}}(k,$ $\displaystyle\mu,z)=$ (20)
$\displaystyle\overline{T}_{\textrm{H{i}}}^{2}(z)\left[b_{\textrm{H{i}}}(z)+\Delta
b_{\textrm{H{i}}}(k,z)f_{\mathrm{NL}}+f(z)\mu^{2}\right]^{2}P_{\mathrm{m}}(k,z)\,,$

where

$\Delta
b_{\textrm{H{i}}}(k,z)=\left[b_{\textrm{H{i}}}(z)-1\right]\frac{3\Omega_{\mathrm{m}}H_{0}^{2}\delta_{\mathrm{c}}}{c^{2}k^{2}T(k)D(z)}$
(21)

and $\delta_{\text{c}}\,{=}\,1.686$ is the critical matter density contrast
for spherical collapse, $T(k)$ is the matter (not foreground) transfer
function, and lastly the growth function can be defined by

$D(z)=\frac{5}{2}\Omega_{\mathrm{m}}H_{0}^{2}H(z)\int_{z}^{\infty}\frac{1+z^{\prime}}{H^{3}\left(z^{\prime}\right)}\mathrm{d}z^{\prime}\,.$
(22)

We set $f_{\text{NL}}\,{=}\,100$ in Equation 20 to provide a clear scale-
dependent bias on large scales in the new observed data simulation. We add
foreground contamination to this following the same steps as the MD1GPC
foreground model, run the PCA foreground clean, then construct the foreground
transfer function following the same methods in the rest of this section. The
lognormal mocks which construct the transfer function will assume
$f_{\text{NL}}\,{=}\,0$ and we will be testing if this transfer function can
still recover the true underlying Hi power spectrum with a scale-dependent
bias induced by the $f_{\text{NL}}\,{=}\,100$ PNG signature. Despite
$f_{\text{NL}}\,{=}\,100$ being confidently ruled out by Planck18
observations, we still use this to set up a highly diverged underlying
cosmology from that assumed in the transfer function construction, similar to
the extreme biases we tested in Figure 17.

Figure 18: Capability to detect "non-standard" cosmology on ultra-large scales
when the fiducial cosmology used in the transfer function construction assumed
standard. For this test, we use the PNG parameter $f_{\text{NL}}$. We let the
underlying truth have an extreme $f_{\text{NL}}\,{=}\,100$, produced using
COLA simulations (black-dashed line). Despite an incorrect
$f_{\text{NL}}\,{=}\,0$ assumption used in the mocks for the transfer
function, the correct $f_{\text{NL}}$ signature is still recovered in the
reconstructed power spectrum (blue data points) across most scales.

Figure 18 shows the results for the large-scale PNG COLA simulation. The
black-dashed line shows the foreground-free result where the impact from the
$f_{\text{NL}}\,{=}\,100$ input is clearly evident on small-$k$. For reference
we also plot the $f_{\text{NL}}\,{=}\,0$ equivalent case shown by the purple-
dotted line. Adding foreground contamination and cleaning is enough to greatly
distort the PNG signature (shown by the red-solid line). For this simulation,
we do not add polarised foregrounds, therefore only a $N_{\text{fg}}\,{=}\,4$
PCA clean is required. This is more representative of a future large sky
survey where it is safer to assume an enhanced level of calibration has been
achieved. The blue data points show the result where we have used the transfer
function to reconstruct the signal loss from the foreground clean. The errors
are estimated from the variance of the transfer function as outlined in
Section 3.3. There is good agreement between the reconstructed power spectrum
and the true underlying power spectrum (black-dashed) where
$f_{\text{NL}}\,{=}\,100$, shown in more detail by the residuals in the bottom
panel. We emphasise that this is an extreme case where the underlying true
cosmology is very different from that assumed by the mocks in the transfer
function, and despite this only three of the measured modes in the
reconstructed power spectrum are more than $1\sigma$ beyond sub-percent
accuracy.

Previous work in Cunnington et al. (2020) investigated signal loss from
foreground cleaning in the context of PNG and demonstrated how a
phenomenological model could be implemented to account for the signal loss.
Whilst this delivered successful simulation-based results, there was a
worrying degeneracy between the parameters in the signal loss model and
$f_{\text{NL}}$, thus tight priors would be needed for such a method to allow
good constraints on $f_{\text{NL}}$. These tight priors would only come from a
very good understanding of foreground contamination and signal loss. The
advantage behind the foreground transfer function approach is that no
phenomenological model is required and the signal loss is reconstructed using
the observed data itself without the need for additional nuisance parameters.
We emphasise that this is a preliminary investigation into PNG with a
foreground transfer function and a more complete study with robust large-scale
simulations is required, including more detailed modelling and MCMC parameter
estimation, which we defer to future work.

The previous results in Figure 17 supported the claim that a discrepancy
between parameters assumed for the transfer function and the underlying truth
in the data, only has a mild influence on results. Figure 18 extends this
conclusion to the largest scales where foreground contamination is most
troublesome and thus is a very encouraging result. A slight negative bias in
the reconstructed power spectrum is apparent, likely caused by the lower
$f_{\text{NL}}$ in the mocks, thus we cannot claim complete independence from
the mocks. Sensible assumptions should therefore still be made when fixing a
fiducial cosmology for the transfer function construction. However, there is
potential for this small bias to be eradicated by adopting the iterative
approach to transfer function calculation as discussed in Appendix D. In this
case, the recovered $f_{\text{NL}}$ inferred from the reconstructed power
spectrum, could be used again to seed new mocks used in a new transfer
function calculation, thus achieving enhanced accuracy. This is another item
we leave for follow-up work, extending beyond the initial discussion and tests
in Appendix D.

##  6 Conclusions

Hi intensity mapping has the potential to be a leading resource for precision
cosmology. A major challenge involves removing astrophysical foregrounds that
dominate the underlying Hi cosmological signal by several orders of magnitude.
Simulations and real observations are providing evidence that foregrounds can
be sufficiently cleaned using blind separation techniques. However,
quantifying the precise impact foreground cleaning has on the Hi power
spectrum is crucial for avoiding biased analyses. In this work, we validate a
method involving mock signal injection into the observed data as a means for
accurately estimating the signal loss induced in the Hi as a consequence of
the foreground clean. This method, referred to as the foreground transfer
function, has been used in real data analysis before and its accuracy was
validated for the results in Cunnington et al. (2022), but its reliability for
the purposes of precision cosmology has not been studied until now. For the
first time, we present simulation-based tests demonstrating the foreground
transfer function’s accuracy as a tool for reconstructing estimated Hi power
spectra to sub-percent accuracies.

This work blackused a selection of simulations to enable tests on a range of
scenarios. In all cases we used an underlying Hi intensity map generated from
an $N$-body simulation, from which we could measure the "true" signal.
Simulated foreground maps were then added and a PCA-based foreground clean was
performed, the consequences of which were analysed relative to the truth. We
varied the complexity of the foregrounds and observational effects providing
scenarios with differing demands from the foreground clean, hence presenting a
range in signal loss. Foreground transfer functions were constructed using
lognormal Hi intensity mapping mocks and their ability to recover the true
signal could be scrutinised. In this work, our focus was on signal loss from
over-cleaning, as opposed to the other undesirable consequence of foreground
residuals caused by under-cleaning. Foreground residuals should be reducible
to sub-dominant contributions or circumvented in cross-correlation with e.g.
galaxy surveys, hence are less problematic than signal loss for precision
cosmology. In the majority of tests, we therefore inspected the reconstructed
cross-correlation of the cleaned intensity map with the original foreground-
free (Hi-only) map. This way, only the effects from signal loss would manifest
and the impact from foreground residuals would be mitigated, which cause a
positive bias in the Hi auto-correlation.

We summarise our main conclusions below;

  * •

We summarised the recipe for estimating a foreground transfer function
(Section 3.1) using mock signal injection into the observed data, which
delivers an unbiased reconstruced power spectrum (Figure 4). These results
included simulated polarisation leakage which demanded an aggressive
foreground clean, resulting in ${>}\,50\%$ signal loss on the largest scales.
In Figure 19 we demonstrated the potential consequences of deviating from this
unbiased recipe resulting in under-estimated signal loss of up to
${\sim}\,30\%$.

  * •

We validated a technique for estimating the covariance in reconstructed power
spectra which involves calculating the covariance across the reconstructed
power spectra from all injected mock realisations (Figure 10). Crucially, when
adopting this approach, the cleaned observed data
$\textbf{{X}}_{\text{clean}}$ must not be subtracted when calculating the
transfer function, which is otherwise removed to reduce the variance in the
transfer function and optimise its accuracy.

  * •

It has been previously assumed that the transfer function should be applied
twice to correct for auto-correlation signal-loss i.e.
$P_{\text{rec}}^{\text{auto}}\,{=}\,P_{\text{clean}}^{\text{auto}}\mathcal{T}^{-2}$.
However, our simulations tests show that this is not the case (Figure 11). At
the fundamental 3D Fourier transform level, the average suppression in signal
is the same in auto-correlation and cross-correlations with a foreground-free
tracer (Figure 12). Thus a $\mathcal{T}^{-1}$ transfer function is the correct
reconstruction in both cross- and auto-correlation power spectra.

  * •

When calculating the transfer function, we estimated power directly in 1D
spherically averaged $k$-bins, finding this performed well on all our
simulations. This deviates from some previous approaches which instead
calculate and apply a transfer function in 2D cylindrically averaged
$k_{\perp},k_{\parallel}$-bins, then re-projects these bandpowers into 1D
$k$-bins. However, we found evidence that this made results prone to a higher
variance (Figure 7), which may require a more detailed weighted average in the
2D to 1D projection to mitigate the issue.

  * •