Statistical estimation of full-sky radio maps from 21cm array visibility data
using Gaussian Constrained Realisations–A

# Statistical estimation of full-sky radio maps from 21cm array visibility
data
using Gaussian Constrained Realisations

Katrine A. Glasscock1 E-mail<EMAIL_ADDRESS>(KAG)
0000-0001-6894-0902 Philip Bull1,2 0000-0001-5668-3101 Jacob Burba1
0000-0002-8465-9341 Hugh Garsden1 0009-0001-3949-9342 Michael J. Wilensky1
1Jodrell Bank Centre for Astrophysics  0000-0001-7716-9312 University of
Manchester  Manchester  M13 9PL  UK
2Department of Physics and Astronomy  University of Western Cape  Cape Town
7535  South Africa

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

###### Abstract

An important application of next-generation wide-field radio interferometers
is making high dynamic range maps of radio emission. Traditional deconvolution
methods like CLEAN can give poor recovery of diffuse structure, prompting the
development of wide-field alternatives like Direct Optimal Mapping and
$m$-mode analysis. In this paper, we propose an alternative Bayesian method to
infer the coefficients of a full-sky spherical harmonic basis for a drift-scan
telescope with potentially thousands of baselines. The can precisely encode
the uncertainties and correlations between the parameters used to build the
recovered image. We use Gaussian Constrained Realisations (GCR) to efficiently
draw samples of the spherical harmonic coefficients, despite the very large
parameter space and extensive sky-regions of missing data. Each GCR solution
provides a complete, statistically-consistent gap-free realisation of a full-
sky map conditioned on the available data, even when the interferometer’s
field of view is small. Many realisations can be generated and used for
further analysis and robust propagation of statistical uncertainties. In this
paper, we present the mathematical formalism of the spherical harmonic GCR-
method for radio interferometers. We focus on the recovery of diffuse emission
as a use case, along with validation of the method against simulations with a
known diffuse emission component.

###### keywords:

reionisation, observations, large scale structure of universe – methods:
statistical – techniques: interferometric

##  1 Introduction

Through measurements of the Cosmic Microwave Background (CMB), significant
progress has been made in studying the very early Universe (Planck
Collaboration et al., 2020), but the subsequent epochs up to and including the
Epoch of Reionisation still leaves much to be discovered. A promising method
to explore these intermediate epochs of cosmic history is the redshifted 21 cm
line from neutral hydrogen. The signal can be used as a cosmological probe of
structure formation, as it traces the distribution and evolution of the
neutral hydrogen that fills the early Universe and forms the first stars and
galaxies during Cosmic Dawn and also provides an inverse tracer of the
reionisation of the intergalactic medium during the later Epoch of
Reionisation (EoR) (Pritchard & Loeb, 2012; Liu & Shaw, 2020). Importantly,
this tracer is redshift dependent, providing line-of-sight distance
information that adds an extra dimension compared to the CMB, (Furlanetto et
al., 2006; Morales & Wyithe, 2010; Liu & Shaw, 2020).

Studies of the 21 cm signal are often divided into either measuring the
spatially averaged _global signal_ , or measuring the statistical properties
of the spatial fluctuations in the brightness temperature field. The 21 cm
signal from the above-mentioned epochs is redshifted to the range $z\sim
27-6$, corresponding to frequencies of $50-$200\text{\,}\mathrm{MHz}$$. A
number of experiments have been built with the purpose of searching for this
signal, and many have already reported upper limits. A non-exhaustive list
includes the Murchison Widefield Array; MWA (Bowman et al., 2013; Tingay et
al., 2013; Wayth et al., 2018), the Donald C. Backer Precision Array to Probe
the Epoch of Reionisation; PAPER (Parsons et al., 2012; Ali et al., 2015), The
LOw Frequency Array; LOFAR (van Haarlem et al., 2013; Patil et al., 2017), and
the Hydrogen Epoch of Reionization Array; HERA (DeBoer et al., 2017; The HERA
Collaboration et al., 2022).

The biggest technical challenge to date is the proper handling of foreground
contaminants and how they are distorted by the receiving instrumentation. The
dominant contaminant is typically the diffuse emission that is present at all
parts of the sky, and is roughly a factor of around $10^{5}$ times brighter
than the underlying 21 cm signal (Liu & Shaw, 2020). The main component of
diffuse emission is the Galactic synchrotron emission, but there are also
contributions from extragalactic free-free emission, bright Galactic radio
point sources, and unresolved extragalactic point sources, (Santos et al.,
2005). In addition there may be as-yet unknown diffuse components, e.g.
associated with apparent excess background emission as claimed by ARCADE 2
(Fixsen et al., 2011; Singal et al., 2010) and OVRO-LWA (Dowell & Taylor,
2018), or potential dark matter annihilation signals (Evoli et al., 2014;
Lopez-Honorez et al., 2016). Being able to robustly separate known but
uncertain Galactic and extra-galactic components from novel sources of
emission would greatly aid the physical interpretation of all of these
effects, particularly as highly sensitive but complex next-generation arrays
such as SKAO come online.

The issue of modelling diffuse emission is especially important for close-
packed radio arrays. Sparse interferometers such as the VLA or SKA-MID are
essentially blind to large-scale emission, as the antenna configuration
resolves out large angular scales. Compact arrays with smaller antennas are
designed to make high dynamic range measurements on large angular scales
however, which means the bright diffuse emission becomes an unavoidable issue.
Imaging with interferometers is complicated, particularly for diffuse
emission, as the effect of the interferometer response (its point spread
function; PSF) must be deconvolved. For regular close-packed arrays, there is
often highly incomplete coverage of the $(u,v)$-plane, since there are many
baselines of the same length. This ‘redundancy’ of the baselines results in
artefacts like grating lobes in the PSF (Dillon & Parsons, 2016).

The typical lack of uniform $(u,v)$-coverage has given rise to the use of
specialised deconvolution algorithms, where the so-called _dirty image_ is
deconvolved to result in a more uniform image with imaging artifacts removed
or suppressed. Many traditional pipelines make use of the CLEAN algorithm,
which works by iteratively removing and replacing bright sources in the dirty
image with their calculated point spread functions (PSFs) to get rid of the
point source side lobes, (Högbom, 1974; Cornwell, T. J., 2009). CLEAN,
however, does not correct for diffuse emission and additionally the statistics
of the resulting image are not well known, particularly as CLEAN acts as a
non-linear transformation of the data. Some improvements have been made to the
CLEAN algorithm through Multi-Scale CLEAN (Cornwell, 2008) and WSCLEAN
(Offringa & Smirnov, 2017), which can be used in either (or both) multi-scale
and multi-frequency mode. Frequency dependence is introduced by dividing the
full bandwidth into several output channels and then treating them
individually. In multi-scale mode, further iterations are introduced using
sub-loops and going through the scales one-by-one by convolving the dirty
image with the corresponding space-kernel. Even so, both methods still build
upon the same CLEAN principle of using maxima to replace sources with an ideal
PSF response.

An alternative approach is Direct Optimal Mapping (DOM), which uses a maximum
likelihood method to estimate the beam-weighted map of the sky using a pixel-
basis, (Dillon et al., 2015; Xu et al., 2022). DOM is thus a one-shot process
that estimates both the mean and the covariance in a lossless way, i.e.
without losing information on model parameters. The model is built around a
measurement matrix $\boldsymbol{\mathbf{A}}$ that maps a pixelated sky map to
the visibilities. The choice of $\boldsymbol{\mathbf{A}}$ can be problematic
however, and regularisation is needed in order to overcome degeneracies and
ensure uniqueness of the map solutions that are obtained by ‘inverting’ the
measuring matrix. Furthermore, the noise covariance $\boldsymbol{\mathbf{N}}$
is also affected by the measurement matrix, requiring a choice for how to
beam-weight the recovered map. While in principle the DOM formulation accounts
for the whole sky, much smaller faceted maps are typically used to reduce the
computational costs.

Like DOM, $m$-mode analysis is also typically implemented as a maximum
likelihood estimator that will produce an estimate of the mean and covariance
of the sky, although now in a spherical harmonic basis instead (Shaw et al.,
2014, 2015; Eastwood et al., 2018). This method is developed specifically for
drift-scan telescopes, taking advantage of simplifications that arise after
applying a spherical transform to both the sky intensity and beam function, as
well as performing a Fourier transform of the visibilities along the time
direction, into $m$-mode space. As with DOM, proper care must be taken to
ensure uniqueness of the solutions, so the inverse of the beam transfer
function is replaced with the Moore-Penrose pseudo-inverse. Both DOM and
$m$-mode analysis are maximum likelihood estimators, and therefore only
provide moments of the underlying posterior distribution of the sky model
parameters (although see Chapter 6 of Kriele (2022), which recasts $m$-mode
analysis as a Bayesian linear model).

In this paper, we introduce a statistical method that, like $m$-mode analysis,
is based on a spherical harmonic description of the sky (although we do not
use an $m$-mode transform here). We use an explicitly Bayesian treatment of
the map reconstruction problem, introducing a prior to the system which
provides a regularisation, particularly in regions of missing data, and using
the Gaussian Constrained Realisation (GCR) method to make sampling of the
large number of spherical harmonic coefficients tractable. Each sample vector
drawn from the posterior is a valid realisation of the whole sky, without any
missing data regions, and we can treat it like we would treat ‘ideal’ data
without noise etc. Uncertainties can be propagated by passing multiple
realisations through subsequent processing steps and then inspecting their
distribution. By using an explicit forward model, we can also avoid ad-hoc
beam weighting steps, which makes interpretation of the results more
straightforward. Furthermore, this sampler can be incorporated in a larger
‘Gibbs sampling’ framework that also samples other sky and instrumental
parameters (e.g. Eriksen et al., 2008; Kennedy et al., 2023; CHIME
Collaboration, 2023).

Naturally, calculating multiple samples will require more computational power
than solving once for the maximum likelihood (or for the maximum _a
posteriori_) solution. We propose some approaches to reduce the computational
burden, such as specialising to the drift-scan telescope case and using Wigner
$\mathfrak{D}$-matrices to account for the sky rotating above the instrument.
This removes the need to explicitly calculate the visibility response function
(which maps the spherical harmonic coefficients to the visibilities) at
multiple times.

The paper is structured as follows: In Section 2 we build and present the
visibility response model and cover the Bayesian methods of Wiener filtering
and Gaussian Constrained realisations. Both methods are related to Gibbs
sampling, a future prospect of this work. Section 3 describes the diffuse
emission foreground model as well as the visibility simulation- and sampling
parameters. In Section 4 we present the performance for our reference
(standard) case as well as a comparative analysis of various noise- and prior
levels and ten other simulation scenarios including varying the field of view
of the array. Finally, in Section 5 we conclude.

##  2 Bayesian sampling of diffuse foregrounds

In this section we develop the mathematical formalism to statistically sample
spherical harmonic modes on the whole sky given radio interferometer
visibility data.

###  2.1 Data model

The data model in this work is based on a single diffuse component modelled as
a sum of spherical harmonic modes at each frequency. Ultimately, this will
only be one component of a more comprehensive sky model involving
contributions from point sources, the EoR signal, etc. (which do not need to
be modelled using a spherical harmonic basis). It would also be possible to
define a particular functional form for the frequency dependence and interpret
the spherical harmonic coefficients as amplitudes of this spectral template,
e.g. by defining the specific intensity of the sky as

$\displaystyle I(\nu,\theta,\phi)$ $\displaystyle=\sum_{\ell,m}a_{\ell
m}Y_{\ell m}(\theta,\phi)f(\nu)$ (1)

with frequency dependence $f(\nu)$. A common choice of frequency dependence is
a power law with spectral index $\beta$. These scenarios involve only simple
modifications to the single component (and per-frequency) model that we
consider here however, and so for the sake of simplicity we will not consider
them further here.

We begin by writing out an expression for the complex visibilities observed by
an interferometer. A baseline separated by vector $\boldsymbol{\mathbf{b}}$
given by antennas $i,j$ roughly probes the sky brightness temperature at
Fourier mode $\boldsymbol{\mathbf{u}}=\boldsymbol{\mathbf{b}}/\lambda$, (Liu &
Shaw, 2020). The visibility that the interferometer measures by that baseline
$\boldsymbol{\mathbf{b}}$ is then given by,

$\displaystyle
V_{ij}(\nu,t)=\int\text{d}^{2}\Omega\,A_{i}(\nu,\boldsymbol{\mathbf{\theta}})A_{j}^{*}(\nu,\boldsymbol{\mathbf{\theta}})\,I(\nu,\boldsymbol{\mathbf{\theta}})\,e^{-2\pi\boldsymbol{\mathbf{u}}(\nu)\cdot\boldsymbol{\mathbf{\theta}}}+n_{ij},$
(2)

where $A_{i}$ and $A_{j}^{*}$ are the E-field beams for each antenna, $I$ is
the specific intensity of the sky, the exponential term describes the fringes
where $\boldsymbol{\mathbf{\theta}}$ is the topocentric coordinates of the
sources, and $n_{ij}$ is the baseline dependent noise.

The per-frequency spherical harmonic expansion of the sky model is given as,

$\displaystyle I(\nu,\theta,\phi)=\sum_{\ell,m}a_{\ell m}(\nu)Y_{\ell
m}(\theta,\phi),$ (3)

where $\theta$ and $\phi$ are the declination and right ascension in
equatorial coordinates respectively. For a typical drift scan array the motion
of the sky is especially simple in an equatorial coordinate system, as it is a
simple right ascension (or LST) rotation.

There are many different conventions on which $a_{\ell m}$ modes to include,
when working with spherical harmonics. For the sake of reducing computational
resources, it makes sense to manipulate the $a_{\ell m}$-vector to contain
only the minimal required modes to solve the system. Firstly, it should be
noted that even though the visibilities are complex entities, the actual sky
needs to be a real field. The following symmetry relation for the $a_{\ell m}$
modes must be satisfied,

$\displaystyle{a_{\ell m}}=(-1)^{m}a^{*}_{l,-m},$ (4)

or, split into real and imaginary parts,

$\displaystyle a_{\ell,+m}^{\rm re}=(-1)^{m}a_{\ell,-m}^{\rm re}\quad{\rm
and}\quad a_{\ell,+m}^{\rm im}=(-1)^{m+1}a_{\ell,-m}^{\rm im}\,.$ (5)

This means that only $m\geq 0$ modes are necessary and all negative $m$ modes
are naturally excluded from the $a_{\ell m}$-vector.

Generally, the spherical harmonic coefficients are complex-valued. This can
present complications when dealing with vectors of $a_{\ell m}$-values
numerically. We therefore split the $a_{\ell m}$-values into their real- and
imaginary parts. The $a_{\ell m}$-vector now consists first of all the real-
and then the imaginary modes. Moreover, since the $m=0$ imaginary parts will
always be zero in order to satisfy the (anti-)symmetry condition of Eq. 5,
these spherical harmonic modes are removed while sampling and injected back in
afterwards. For a given $\ell_{\text{max}}$ this leads to a total number of
$a_{\ell m}$ modes of $N_{\text{modes}}=({\ell_{\text{max}}}+1)^{2}$.

We next define the visibility response operator, $\delta V_{ij}^{\ell
m}(\nu,t)$, which gives the projection from a spherical harmonic vector to a
set of radio interferometer visibilities. The visibility response is dependent
on the antenna array configuration and location, which means the primary beam
function is also absorbed into this operator. The full visibility will then be
the visibility response function applied to the $a_{\ell m}$ modes summed over
all $(\ell,m)$-values,

$\displaystyle V_{ij}=\sum_{\ell m}\delta V_{ij}^{\ell m}(\nu,t)\,a_{\ell m}.$
(6)

where the visibility response is defined as,

$\displaystyle\delta V_{ij}^{\ell
m}(\nu,t)=\int\text{d}^{2}\Omega\,A_{i}(\nu,\boldsymbol{\mathbf{\theta}})A_{j}^{*}(\nu,\boldsymbol{\mathbf{\theta}})\,Y_{\ell
m}(\alpha,\delta)\,e^{-2\pi\boldsymbol{\mathbf{u}}_{ij}(\nu)\cdot\boldsymbol{\mathbf{\theta}}}.$
(7)

In this expression, $\boldsymbol{\mathbf{\theta}}$ are in topocentric
coordinates (e.g. local altitude/azimuth), and
$\alpha=\alpha(\boldsymbol{\mathbf{\theta}},t)$ and
$\delta=\delta(\boldsymbol{\mathbf{\theta}},t)$ are the RA and Dec coordinates
corresponding to a given topocentric pointing $\boldsymbol{\mathbf{\theta}}$
at local sidereal time $t$. This operator can be computed ahead of time for a
given array configuration, which determines the available baseline vectors
$\mathbf{u}$ and set of antenna E-field beams $A_{i}$.

###  2.2 Wigner $\mathfrak{D}$-matrix formalism

In the above section, the visibility response is defined as a function of both
time and frequency. Instead of simulating the visibility response operator for
each time separately, we can choose a single reference time and apply a
rotation to get the response at any other desired LST. This is made simpler by
choosing the spherical harmonic basis to align with the rotation axis of the
sky as seen by a drift-scanning telescope, i.e. by defining the spherical
harmonic basis in equatorial coordinates. In this case, a simple azimuthal
rotation around the celestial axis implements the mapping between RA and LST,
while the declination stays constant.

We begin by defining the rotational matrix $\mathcal{R}_{\ell
m\ell^{\prime}m}$ that transforms the spherical harmonic modes at a reference
sidereal time $t_{\text{ref}}$ to the updated sidereal time,

$\displaystyle a_{\ell m}(t)$ $\displaystyle=\mathcal{R}_{\ell
m\ell^{\prime}m^{\prime}}(t)\,a_{\ell^{\prime}m^{\prime}}(t_{\text{ref}}).$
(8)

For spherical harmonics, this is given by the Wigner $\mathfrak{D}$-matrix,
which is a unitary matrix in an irreducible representation of the SO(3) group
(Tung, 1985). The spherical harmonics transform as

$\displaystyle Y_{\ell}^{m}(\theta^{\prime},\phi^{\prime})$
$\displaystyle=\sum_{m^{\prime}=-\ell}^{\ell}Y_{\ell}^{m^{\prime}}(\theta,\phi)\,{\mathfrak{D}}_{m^{\prime}m}^{\ell}(\alpha,\beta,\gamma),$
(9)

where ${\mathfrak{D}}_{mm^{\prime}}^{\ell}(\alpha,\beta,\gamma)$ is the
$\mathfrak{D}$-matrix given by the three Euler angles (that describe
sequential rotations around three axes). Note, that for any rotation with
$\mathfrak{D}$-matrices, the spherical harmonic of degree $\ell$ and order $m$
transforms into a linear combination of spherical harmonics to the same degree
$\ell$. The $\mathfrak{D}$-matrix itself is given as

$\displaystyle{\mathfrak{D}}^{\ell}_{m^{\prime}m}(\alpha,\beta,\gamma)$
$\displaystyle=e^{-im^{\prime}\alpha}\,d^{\ell}_{m^{\prime}m}(\beta)\,e^{-im\gamma},$
(10)

where $d^{\ell}_{m^{\prime}m}(\beta)$ are the corresponding reduced Wigner
matrices. A table of the (small) $d$-matrices can be found in Dong (1994).
Note the lack of mixing between $\ell$ modes.

For a zenith-pointing drift-scan telescope, the sky rotation in a time
interval $t-t_{\rm ref}$ can be mapped directly to an azimuthal rotation
angle. The full visibility given in Eq. 6 can then be written as,

$\displaystyle V_{ij}(\nu,t)=\sum_{\ell m}\delta V_{ij}^{\ell
m}(\nu,t_{\text{ref}})\,\sum_{m^{\prime}}{\mathfrak{D}}^{\ell}_{m\,m^{\prime}}(t)\,a_{\ell
m^{\prime}}(t_{\rm ref}),$ (11)

where $\delta V_{ij}^{\ell m}(\nu,t_{\text{ref}})$ is the pre-computed
visibility response at a reference time $t_{\text{ref}}$,
${\mathfrak{D}}^{\ell}_{m\,m^{\prime}}(t)$ is the appropriate
$\mathfrak{D}$-matrix for a given LST, and $a_{\ell m^{\prime}}$ are the
spherical harmonic modes of the sky at $t_{\text{ref}}$. If we let $p$ and $q$
label the frequency and time samples, $\nu_{p}$ and $t_{q}$, we can rewrite
the above in index notation as

$\displaystyle V_{ijpq}={X}_{ijpq\ell m}\,\,a_{\ell m},$ (12)

where repeated indices denote summation, and ${X}_{ijpq\ell m}$ is the
combined operator for the two sums in Eq. 11. The combined visibility response
and $\mathfrak{D}$-matrices now form the full
$\boldsymbol{\mathbf{X}}$-operator, which is linear and equivalent to that of
Eq. 6.

###  2.3 Wiener filter – the maximum a posteriori solution

Now, having established our visibility model based on the linear operator
$\boldsymbol{\mathbf{X}}$ for the visibility response using a spherical
harmonics basis — we want to obtain the posterior distribution for the
spherical harmonic modes given the data, noise, and priors, in order to
generate samples of the diffuse emission sky. The steps outlined in this
section and the next are steps of increasing complexity in a Bayesian
hierarchy, ultimately preparing for implementing this model into a full Gibbs
sampling scheme. Here, we go from the Wiener filter solution; a simple maximum
a posteriori (MAP) solution, to drawing independent samples from the Gaussian
Constrained Realisation (GCR) equation.

Thereby, as a first step in the Bayesian hierarchy we find the Wiener filter
solution, which also acts as the central value for drawing the realisations.
From Bayes’s theorem we can get the conditional distribution on the $a_{\ell
m}$ modes,

$\displaystyle
P(\boldsymbol{\mathbf{a}}\,|\,\boldsymbol{\mathbf{d}},\boldsymbol{\mathbf{N}},\boldsymbol{\mathbf{a}}_{0},\boldsymbol{\mathbf{S}})\propto\,$
$\displaystyle
P(\boldsymbol{\mathbf{d}}|\,\boldsymbol{\mathbf{a}},\boldsymbol{\mathbf{N}})\,P(\boldsymbol{\mathbf{a}}|\,\boldsymbol{\mathbf{a}}_{0},\boldsymbol{\mathbf{S}}),$
(13)

when conditioning on the known covariances of the signal
$\boldsymbol{\mathbf{S}}$ and noise $\boldsymbol{\mathbf{N}}$, the data-vector
$\boldsymbol{\mathbf{d}}$, and the mean of the prior on the $a_{\ell m}$ modes
$\boldsymbol{\mathbf{a}}_{0}$. For simplification of the notation we have
dropped the $(\ell,m)$-indices and simply denote the $a_{\ell m}$-vector as
$\boldsymbol{\mathbf{a}}$ as well as dropping the baseline indices ($i,j$). It
is assumed that generally the signal prior and covariance will be independent
of the data and noise covariance. For simplicity we are keeping
$\boldsymbol{\mathbf{S}}$ fixed, however, in a full Gibbs sampling scheme it
would be possible to sample $\boldsymbol{\mathbf{S}}$ as well. Furthermore, we
assume that the noise is Gaussian thus rewriting Eq. 13 as:

$\displaystyle
P(\boldsymbol{\mathbf{a}}\,|\,\boldsymbol{\mathbf{d}},\boldsymbol{\mathbf{N}},\boldsymbol{\mathbf{a}}_{0},\boldsymbol{\mathbf{S}})\propto
e^{-(\boldsymbol{\mathbf{d}}-\boldsymbol{\mathbf{Xa}})^{\dagger}\boldsymbol{\mathbf{N}}^{-1}(\boldsymbol{\mathbf{d}}-\boldsymbol{\mathbf{Xa}})}e^{-(\boldsymbol{\mathbf{a}}-\boldsymbol{\mathbf{a}}_{0})^{\dagger}\boldsymbol{\mathbf{S}}^{-1}(\boldsymbol{\mathbf{a}}-\boldsymbol{\mathbf{a}}_{0})}.$
(14)

For a Gaussian distribution the maximum of the posterior distribution is the
same as the maximum of the log-posterior, hence we can find the Wiener filter
solution by maximising the partial derivative of the exponent of Eq. 14 with
respect to the signal, $\boldsymbol{\mathbf{a}}$:

$\displaystyle\left.\frac{\partial}{\partial\boldsymbol{\mathbf{a}}}\right|_{\boldsymbol{\mathbf{a}}=\boldsymbol{\mathbf{\hat{a}}}}\left(-(\boldsymbol{\mathbf{d}}-\boldsymbol{\mathbf{Xa}})^{\dagger}\boldsymbol{\mathbf{N}}^{-1}(\boldsymbol{\mathbf{d}}-\boldsymbol{\mathbf{Xa}})-(\boldsymbol{\mathbf{a}}-\boldsymbol{\mathbf{a}}_{0})^{\dagger}\boldsymbol{\mathbf{S}}^{-1}(\boldsymbol{\mathbf{a}}-\boldsymbol{\mathbf{a}}_{0})\right)=0,$
(15)

resulting in,

$\displaystyle\boldsymbol{\mathbf{d}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}+\boldsymbol{\mathbf{a_{0}}}^{\dagger}\boldsymbol{\mathbf{S}}^{-1}=\boldsymbol{\mathbf{a}}^{\dagger}\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}+\boldsymbol{\mathbf{a}}^{\dagger}\boldsymbol{\mathbf{S}}^{-1}.$
(16)

Now, taking advantage of the fact that covariance matrices are Hermitian, we
can complex conjugate the entire expression and rearrange to obtain the Wiener
filter solution,

$\displaystyle\left[\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}+\boldsymbol{\mathbf{S}}^{-1}\right]\boldsymbol{\mathbf{a}}_{\textup{wf}}=\left(\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}+\boldsymbol{\mathbf{S}}^{-1}\boldsymbol{\mathbf{a}}_{0}\right).$
(17)

Returning to Eq.14; multiplying two multivariate Gaussians will result in a
new multivariate Gaussian. By completing the square it can be shown that this
can be written as proportional to a multivariate Gaussian with inverse
covariance matrix given as

$\displaystyle\boldsymbol{\mathbf{\Sigma}}^{-1}=\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}+\boldsymbol{\mathbf{S}}^{-1},$
(18)

and where the Wiener filter solution acts as the mean,

$\displaystyle\boldsymbol{\mathbf{\hat{a}}}=\boldsymbol{\mathbf{\Sigma}}\left(\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}+\boldsymbol{\mathbf{S}}^{-1}\boldsymbol{\mathbf{a}}_{0}\right).$
(19)

Even though the Wiener filter is the maximum a posteriori solution, it is
generally a biased estimator as described in more detail in Kennedy et al.
(2023). Furthermore, the Wiener filter is a summary statistic but we are
interested in generating actual samples that are complete and statistically-
consistent realisations of the full-sky map. For this, the Wiener filter can
instead act as the central value for drawing the samples.

Figure 1: The absolute values of the visibility response operator
$X_{\text{re}}$ on a log-scale and with dimensions (Nmodes) $\times$ (NLSTs
$\times$ Nfreq $\times$ Nbl). Here it is shown for $\ell_{\text{max}}=20$, LST
$=$ 0–8 $\mathrm{h}$, $\nu=$ 100–101 $\mathrm{MHz}$ and 10 close packed
antennas that form 45 baselines as shown in Fig. 2. The darker is is, the less
that frequency/LST/baseline contributes to measuring that $(m,\ell)$-value.
Note that this is the _real-part only_ of $X$ (the imaginary part showing
similar structure) but it covers both the real- and imaginary parts of the
$a_{\ell m}$ modes.

###  2.4 Gaussian Constrained Realisations

In order to draw samples from the full conditional distribution
$P(\boldsymbol{\mathbf{a}}\,|\,\boldsymbol{\mathbf{d}},\boldsymbol{\mathbf{N}},\boldsymbol{\mathbf{a}}_{0},\boldsymbol{\mathbf{S}})$,
we can take advantage of the fact that we can describe our model as a
multivariate Gaussian distribution. The Wiener filter solution described by
its mean and covariance given in eqs. 18 and 19 acts as a first step in the
hierarchy to which we can add random normal realisation terms of the signal
$\boldsymbol{\mathbf{\omega}}_{a}$ and noise terms
$\boldsymbol{\mathbf{\omega}}_{d}$ correctly scaled by their respective
covariances to draw the constrained realistations
$\boldsymbol{\mathbf{a}}_{\textup{cr}}$, (Eriksen et al., 2008). This leads to
the Gaussian Constrained Realisation (GCR) equation,

$\displaystyle\left[\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}+\boldsymbol{\mathbf{S}}^{-1}\right]\boldsymbol{\mathbf{a}}_{\textup{cr}}=\left(\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}+\boldsymbol{\mathbf{S}}^{-1}\boldsymbol{\mathbf{a}}_{0}+\boldsymbol{\mathbf{S}}^{-\frac{1}{2}}\boldsymbol{\omega}_{a}+\boldsymbol{\mathbf{X}}^{\dagger}\boldsymbol{\mathbf{N}}^{-\frac{1}{2}}\boldsymbol{\mathbf{\omega}}_{d}\right).$
(20)

In turn Eq. 20 reduces back to the mean
$\langle\boldsymbol{\mathbf{a}}_{\textup{cr}}\rangle=\boldsymbol{\mathbf{\hat{a}}}$,
since the unit variance Gaussian vectors $\boldsymbol{\mathbf{\omega}}_{a}$
and $\boldsymbol{\mathbf{\omega}}_{d}$ have zero mean. By solving Eq. 20
repeatedly we can draw independent samples of the $a_{\ell m}$ modes on the
sky consistent with the given data vector $\boldsymbol{\mathbf{d}}$, chosen
signal prior $\boldsymbol{\mathbf{a}}_{0}$, the visibility response operator
$\boldsymbol{\mathbf{X}}$, and the covariances $\boldsymbol{\mathbf{S}}$ and
$\boldsymbol{\mathbf{N}}$. For the sake of keeping computational costs low, we
have kept to 100 samples per set of parameters for this paper.

In the case studied in this paper, there is only missing data outside of the
FOV defined by the primary beam and the horizon; both of which is dealt with
through the visibility response function. In case of actual gaps in the data,
for instance due to RFI flagging, one can define a set of weights
$\boldsymbol{\mathbf{w}}$ to redefine the inverse noise covariance,

$\displaystyle\boldsymbol{\mathbf{N}}_{\boldsymbol{\mathbf{w}}}^{-1}=\boldsymbol{\mathbf{ww}}^{T}\circ\boldsymbol{\mathbf{N}}^{-1},$
(21)

where $\circ$ denotes element-wise multiplication. Using this weighted noise
covariance in Eq. 20 means we automatically zero the contribution from the
data inside the mask. The signal covariance, however, does not necessarily go
to zero inside the masked regions and thereby _takes over_ the signal
estimation in the absence of information from the likelihood function. Note
that this is not equivalent to drawing samples from the prior distribution and
simply filling-in the masked regions. The samples generated within the flagged
regions will be informed by both the prior and data from the un-flagged
regions.

In the future this process will be done within the context of Gibbs sampling
(Geman & Geman, 1984), which is a Bayesian method to sample directly from the
joint posterior of more complicated high-dimensional distributions. In this
instance, we focus only on the distribution for the spherical harmonic
coefficients, and use a conjugate gradient solver to obtain
$\boldsymbol{\mathbf{a}}_{\textup{cr}}$ from the GCR equation. Most solvers do
not handle complex quantities well so to overcome this, we saw it necessary to
‘realify’ the full system as described in App. A.

When implemented into a Gibbs sampling scheme, it will enable us to also
sample the signal covariance $\boldsymbol{\mathbf{S}}$. Gibbs sampling is an
iterative method, that samples from each conditional distribution in turn and
thereby effectively samples from the joint posterior. For each iteration it
updates the conditional variables with the samples obtained from the previous
step,

$\displaystyle\boldsymbol{\mathbf{a}}_{i+1}$ $\displaystyle\leftarrow
P\\!\left(\boldsymbol{\mathbf{a}}_{i}\,|\,\boldsymbol{\mathbf{d}},\boldsymbol{\mathbf{N}},\boldsymbol{\mathbf{a_{0}}},\boldsymbol{\mathbf{S}}_{i}\right),$
(22a) $\displaystyle\boldsymbol{\mathbf{S}}_{i+1}$ $\displaystyle\leftarrow
P\\!\left(\boldsymbol{\mathbf{S}}_{i}\,|\,\boldsymbol{\mathbf{a}}_{i+1}\right),$
(22b)

where $\leftarrow$ represents generating samples through evaluating the
conditional distribution and $i$ indexes the Gibbs iteration number. First,
the ${a_{\ell m}}$ modes are sampled and it is noticed that the conditional
distribution takes the same form as in Eq. 13. We can therefore use the same
arguments as before; we are still dealing with a multivariate Gaussian
distribution and therefore end up at the GCR equation of Eq. 20. The sample
obtained from Eq. 22a already contains the conditioning on the noise
covariance $\boldsymbol{\mathbf{N}}$, data $\boldsymbol{\mathbf{d}}$, and
prior mean $\boldsymbol{\mathbf{a_{0}}}$, so there is no need to condition
again on these parameters when sampling $\boldsymbol{\mathbf{S}}$ in Eq. 22b.
We leave a full exploration of sampling the signal covariance (which encodes
the spherical harmonic angular power spectra) to later work.

Figure 2: The standard array layout used in the majority of the simulations,
consisting of 10 closed packed antennas with a diameter of
$14\text{\,}\mathrm{m}$ and a separation of $14.6\text{\,}\mathrm{m}$.

##  3 Simulations and sky model

The structure of the simulated visibility response function has already been
laid out in Sec. 2.1 where it is also made clear that the simulation is
dependent on which LSTs are chosen, $\ell_{\text{max}}$, frequency $\nu$,
baseline orientation and length $b$, and dish size $\theta_{D}$. This section
will describe how the simulations were obtained as well as the standard
parameters used for the reference results, which from here on will be referred
to as the _standard_ case. More realisations have been drawn using a variation
of both the simulation/observation parameters as well as samples from using a
variation of noise and prior levels. The details of these parameter variations
are explained in Sections 4.2–4.4.

###  3.1 Precomputing the visibility response operator

To simulate the visibility response function we used the
hydra111https://github.com/HydraRadio/ code, which builds upon matvis
(Kittiwisit et al., 2023) to simulate the visibility per $a_{\ell m}$-mode as
well as per baseline, frequency, and LST. The absolute values of the real-part
of the visibility response operator, $\boldsymbol{\mathbf{X}}_{\rm re}$, is
shown in Fig. 1 also showing how the operator’s dimensions are packed. The
visibility response is only simulated for Stokes-I. We have assumed the most
general case where we want to solve for a different set of ${a_{\ell m}}$
modes for each frequency channel. The simulations cover the narrow bandwidth
of $100$–$102\text{\,}\mathrm{MHz}$ in just two channels, in order to keep the
problem small for the time being. The LST range spans
$0$–$8\text{\,}\mathrm{h}$ in 10 steps. The primary beam model is calculated
using pyuvsim222https://github.com/RadioAstronomySoftwareGroup/pyuvsim which
provides an analytic Gaussian beam model assumed to be identical for each
antenna, with the width $\theta\simeq\lambda/b$ where $b$ is length of the
baseline $\boldsymbol{\mathbf{b}}$. As standard the distance between each
antenna is $14.6\text{\,}\mathrm{m}$ (HERA-like) and the dishes themselves are
modelled as hyperbolic dishes of diameter
$\theta_{\text{D}}=$14\text{\,}\mathrm{m}$$ with no side-lobes. The model is
simple, but has been specifically chosen so that we can clearly understand and
validate the results. The impact of different beams is an interesting study in
its own right, which we defer to future work.

The standard array layout used in the simulations is a small closed-packed
redundant array to represent a subsection of the full HERA array. A redundant
layout benefits from sampling the same modes many times, thus increasing
sensitivity especially to the relatively large angular scales relevant for EoR
experiments. In this configuration we consider only 10 receivers as shown in
Fig. 2. The subset is sufficiently small to both keep computational costs low
and with its resulting 45 baselines it still offers multiple redundant
baseline groups and variation in both direction and baseline length. As with
HERA the array used in the simulations is a drift-scan instrument pointing
towards zenith placed at latitude
$\sim$-31\text{\,}\mathrm{\SIUnitSymbolDegree}$$. From Fig. 1 it is already
clear that some baselines contribute more than others. For instance the
baseline of antennas $(3,6)$ is not very responsive, which is the long E-W
$43.8\text{\,}\mathrm{m}$ baseline (Fig. 2) whereas all the short
$14.6\text{\,}\mathrm{m}$ baselines contribute significantly more.

###  3.2 Diffuse emission sky model

As we aim to sample the diffuse emission foregrounds specifically, no other
foregrounds are included in the foreground model. For the purely diffuse
emission foreground sky model we use the Global Sky Model (GSM2016; Zheng et
al., 2017) as implemented by pyGDSM (Price, 2016). GSM2016 is based on
principal component analysis of a large set of multi-frequency datasets
(spanning $10\text{\,}\mathrm{MHz}$ to $5\text{\,}\mathrm{THz}$) as well as
performing a blind separation into physical components revealing five
components identified as synchrotron emission, free-free emission, cold dust,
warm dust, and the CMB anisotropy. The sky model is provided as frequency
dependent HEALpix maps where we use only the GSM2016 map at reference
frequency $100\text{\,}\mathrm{MHz}$. The foreground map is then passed to
HEALpy to get the true $a_{\ell m}$ modes,
$\boldsymbol{\mathbf{a}}_{\text{true}}$. As standard, the resolution is set to
nside $=128$, which corresponds to a HEALPix angular pixel size of
$0.46\text{\,}\mathrm{\SIUnitSymbolDegree}$. The maximum mode of the operator
is set to ${\ell_{\text{max}}}=20$, however, which corresponds to
$\sim$9\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The true sky is used as input
for the full data model, which is obtained by applying the visibility response
to the true $a_{\ell m}$ modes.

###  3.3 Data and noise model

The data input for the analysis in this work is based on the visibility
response simulations. In the future, this will be done with real visibility
data instead, but for now simulations with a known diffuse sky component serve
as a means of validation. Here, we use the visibility response operator
$\boldsymbol{\mathbf{X}}$ to map the true-sky spherical harmonic coefficient
vector $\boldsymbol{\mathbf{a}}_{\rm true}$ into visibilities and add a noise
vector $\boldsymbol{\mathbf{n}}$ to the simulated visibilities,

$\displaystyle\boldsymbol{\mathbf{d}}=\boldsymbol{\mathbf{X}}\boldsymbol{\mathbf{a}}_{\text{true}}+\boldsymbol{\mathbf{n}}.$
(23)

The noise on the data model is added as uncorrelated Gaussian random noise,

$\displaystyle\boldsymbol{\mathbf{n}}\sim\mathcal{N}(0,\boldsymbol{\mathbf{N}}),$
(24)

given by the noise covariance $\boldsymbol{\mathbf{N}}$ with components
$N_{ij}$ modelled by the simulated auto-correlation visibilities $V_{ii}$ and
$V_{jj}$ as given by the radiometer equation,

$\displaystyle N_{ij}=\sigma_{ij}^{2}=\frac{V_{ii}V_{jj}}{\Delta t\Delta\nu}.$
(25)

Since the foregrounds are spectrally smooth, there is no need for very high
spectral resolution and we choose $\Delta\nu=$1\text{\,}\mathrm{MHz}$$. The
time resolution is set to $\Delta t=$60\text{\,}\mathrm{s}$$, both to avoid
issues with sky rotation (smearing) and to still ensure good signal-to-noise
ratio.

Lastly, we have applied a $10\%$ prior on the ${a_{\ell m}}$ values. To
implement this, we have defined a prior covariance $\boldsymbol{\mathbf{S}}$
that is diagonal, with values
${S}_{nm}=(0.1\,a_{\text{true},n})^{2}\delta_{nm}$. The prior mean,
$\boldsymbol{\mathbf{a}}_{0}$, is not set equal to
$\boldsymbol{\mathbf{a}}_{\text{true}}$ however, as this would effectively be
inputting information about the correct answer into the inference ahead of
time, which is not realistic. Instead, we draw values for
$\boldsymbol{\mathbf{a}}_{0}$ from a Gaussian distribution centered on the
true $a_{\ell m}$ values, with the same prior covariance, i.e.
$\boldsymbol{\mathbf{a}}_{0}\sim\mathcal{N}(\boldsymbol{\mathbf{a}}_{\text{true}},\boldsymbol{\mathbf{S}})$.
This ensures that the prior mean that we input in the GCR equation (Eq. 20) is
set consistently with the chosen prior model (GSM2016), while deviating from
the true values of the parameters (as would be the case in a real analysis).

Figure 3: The real- and imaginary parts of the difference of the mean of 100
samples of $a_{\ell m}$ modes from the GCR solver and the true sky using the
_standard_ configuration and parameters. Any _outliers_ from the central
region is marked with $\blacktriangle$ and it is notable that these occur more
frequently in the low-$\ell$ region. As described in Sec. 2.1 the $m=0$
imaginary modes ($\times$) should always be zero and the GCR solver therefore
does not solve for this.

Figure 4: Maps generated from the estimated spherical harmonic modes on the
sky using the _standard_ configuration and parameters. _Upper left:_ The true
sky given by pyGSM with ${\ell_{\text{max}}}=20$ and nside$=128$. Note that
the rippled structure comes from the true sky and not from the GCR solver.
_Upper right:_ The estimated sky based on the mean of 100 samples from the GCR
solver. The spherical harmonic modes can be seen in Fig. 3. _Bottom left:_
Fractional difference between the mean and the true sky adjusted to show
differences $<10\%$, which coincides with the beam region, see also Fig. 5 for
a closeup of the region. _Bottom right:_ The standard deviation of the
difference of the estimated and true sky.

##  4 Results

In the following we show results of the samples obtained under various
different parameter configurations of the sampler itself and by investigating
the effects of different observation scenarios. The purpose of this analysis
is to demonstrate the basic behaviours of the method and to validate that it
can recover the true sky to a reasonable level. The comparative analysis uses
simulations in order to understand the behaviour we would expect from real
data as a function of different LST coverage, the effects of effectively
truncating the visibility response operator at different levels of
$\ell_{\text{max}}$, and the impact of increasing the FoV of the array. We
also dive into a stress-test of the prior- vs likelihood levels, making sure
that neither is too broad nor too narrow.

Before going into the comparative analysis, we show in more detail in Sec. 4.1
the results of our chosen reference – or, _standard_ – case, using the general
parameters described in Sec. 3, considering not only the fractional difference
between the posterior mean and true sky map but also looking at the recovered
$a_{\ell m}$ modes, followed by the full comparative analysis in Secs.
4.2–4.3. Unless otherwise stated, we draw 100 samples per scenario, as we
deemed this sufficient for calculating the relevant statistics of the
recovered spherical harmonic modes. As a precautionary check of this
assumption, we also consider a high-sample size case of
$N_{\text{samples}}\approx 3500$. Finally, Sec. 4.4 is dedicated to a closer
study of the special case of increasing the field of view by decreasing the
diameter of the antennas to $\theta_{\text{D}}=$3\text{\,}\mathrm{m}$$, which
results in a FWHM at $100\text{\,}\mathrm{MHz}$ of
$\sim$60\text{\,}\mathrm{\SIUnitSymbolDegree}$$. Having a greater sky coverage
should improve on the sensitivity on lower $\ell$ modes. A large FoV can be
achieved with instruments like the MWA.

Figure 5: Cartesian projection of the fractional difference from the true sky
using the _standard_ configuration and parameters. The RA and dec ranges have
been narrowed to focus on the primary beam region. The FWHM of
$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ is illustrated in the corner of
the figure. The white horizontal lines indicate a FWHM distance from the
centre of the primary beam at $-31.7\text{\,}\mathrm{\SIUnitSymbolDegree}$.
The simulated LST range is $0$-$8\text{\,}\mathrm{h}$ corresponding to an RA
range of
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$-$120\text{\,}\mathrm{\SIUnitSymbolDegree}$.
The centre of the beam-region has fractional differences to the true sky of
$<5\%$ and at the bounds $<10\%$. Figure 6: Approximate range of $\ell$-value
sensitivity for several baselines, as a function of frequency. The baselines
displayed here are the five unique baseline lengths of the 10-dish standard
configuration. The width (shaded regions) is given by the FWHM of the beam
$\delta\ell\simeq\pi/{\rm FWHM}$, shown here for the HERA-like case with FWHM
$=$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The two horizontal lines
indicate an ${\ell_{\text{max}}}=30$ (top, dark purple) and
${\ell_{\text{max}}}=20$ (bottom, magenta) and the black vertical line is set
at $\nu=$100\text{\,}\mathrm{MHz}$$, the frequency used for the simulations in
this work.

###  4.1 Realisations of the standard configuration

Figure 7: Cartesian projection of the fractional difference between the
posterior mean and true sky for multiple sets of samples with varied levels of
noise covariance $\boldsymbol{\mathbf{N}}$, prior covariance
$\boldsymbol{\mathbf{S}}$, prior mean $\boldsymbol{\mathbf{a}}_{0}$, and
increased sample size compared to the _standard_ case. The fractional
difference is defined as defined as
$(\mu({\boldsymbol{\mathbf{a}}}_{\textup{cr}})-\boldsymbol{\mathbf{a}}_{\textup{true}})/\boldsymbol{\mathbf{a}}_{\textup{true}}$,
where $\mu$ indicates the sample mean. The two horizontal lines indicate a
FWHM distance from the centre of the beam.

The results from the Gaussian constrained realisations of the $a_{\ell m}$
modes for the standard case are shown in Figs. 3, 4, and 5. First, we define
the true-sky subtracted mean as the mean of the spherical harmonic
coefficients sampled by the GCR method subtracted by the true sky (input
model) spherical harmonic modes from pyGDSM. In Fig. 3 the true-sky subtracted
mean of 100 realisations of the $a_{\ell m}$ modes is shown with corresponding
error bars defined by the standard deviations of the samples. To be able to
show the details of the region around a difference of zero, the figure is
cropped to $\pm 100$ with any modes outside of this region are marked with a
$\blacktriangle$ and their values. Since the $m=0$ imaginary modes are
artificially injected back in with a value of 0 (cf. Eq. 5), they are without
error bars but have been included and marked with $\times$ in the figure for
clarity. It is clear that the higher order modes are easiest to constrain and
that the lower order modes are far from the true sky.

In Fig. 6 the $\ell$-mode sensitivity given the baseline length and frequency
is shown. The shortest, and most frequent, baseline length in the standard
configuration is $b=$14.6\text{\,}\mathrm{m}$$ corresponding to an
$\ell_{\textup{min}}\sim(\pi b)/\lambda\sim 15$ at a frequency of
$100\text{\,}\mathrm{MHz}$. Note, that at this frequency the baseline length
(in metres) roughly corresponds to the $\ell$-value. Thus we are not expecting
to be sensitive to $\ell\sim$ few, since modes with
$\Delta\theta\gtrsim\lambda/b_{\rm min}$ are resolved out by the
interferometer, and so we should not expect the data to constrain them. The
recovered values of these low-$\ell$ modes are instead prior-driven.

In fact, Fig. 3 is showing the spherical harmonic coefficients, each of which
integrates information across the whole sky. Since the observed region is
small (only $\sim$ 8.6% of the sky), the recovered $a_{\ell m}$ values are
generally mostly determined by the prior even when the noise level is very low
in the observed region, we suspect this is because the prior-dominated area is
much larger. Hence, plotting the statistics of the recovered $a_{\ell m}$
values directly is not a very sensitive test of this method.

Instead, it is much more useful to look at the recovered sky in map-space, as
the difference between the observed and unobserved region is much clearer laid
out. In Fig. 4 we show the maps of the true input sky (top left) given as
described in Sec. 3.2 compared to the recovered map defined by the mean of the
posterior (top right). It is noted that the recovered map shows lower values
for the large scale modes (just as the ${a_{\ell m}}$-values in Fig. 3)
compared to the true sky map, which we expect is due to a random fluctuation
rather than a bias. Fig. 4 also shows the fractional difference between these
two maps (bottom left) and the standard deviation of the recovered sky (bottom
right). The fractional difference between the posterior mean and true sky is
lowest within the primary beam region and the standard deviation is much
smaller in this area too, thus demonstrating the specific patch of sky that is
directly observed versus that constrained by the prior. This explains how we
are still able to recover most of the galactic structure, despite the majority
of it being outside the observed region of the sky. Note, that between
RA$=$150\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and
RA$=$180\text{\,}\mathrm{\SIUnitSymbolDegree}$$ there is another small low
residual region in the fractional difference map, however, this does not
appear in the standard deviation map. This low residual is simply due to a
random fluctuation rather than the true sky being well recovered in this
region, and it is therefore ignored.

Now focusing on the primary beam region alone, we show a Cartesian projection
of the fractional difference between the posterior mean and true sky in Fig.
5. The region is defined to be within a FWHM on either side of the central
latitude of $-31.7\text{\,}\mathrm{\SIUnitSymbolDegree}$. For the standard
case with the antenna diameter of $\theta_{\text{D}}=$14\text{\,}\mathrm{m}$$
the FWHM is $12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$. For the standard case
the sky is recovered to within 5% close to the centre of the beam region and
to within 10% within the FWHM bounds.

###  4.2 Noise and prior levels

Figure 8: The absolute value of the fractional difference between the posterior mean and true sky taken as a slice through the centre of the primary beam region at dec=$-31.7\text{\,}\mathrm{\SIUnitSymbolDegree}$. Artefacts due to the pixelisation have been smoothed away by applying a pair of moving average filters to each of these curves. For all runs in this plot the LST range is $0-$8\text{\,}\mathrm{h}$$ corresponding to RA values of $0-$120\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The light shaded region covers the FWHM of $12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ and the dark shaded region is outside of the primary beam. Within the LST range (white region) the standard and A1 case have the lowest fractional difference between the posterior mean and true sky. Figure 9: Cartesian projection of the fractional difference between the posterior mean and true sky for runs with variations to the standard observation/simulation parameters such as increasing the ${\ell_{\text{max}}}$, adding antennas for longer baselines, increasing the LST range, and increasing the cadence of observations. The fractional difference is defined as defined as $(\mu({\boldsymbol{\mathbf{a}}}_{\textup{cr}})-\boldsymbol{\mathbf{a}}_{\textup{true}})/\boldsymbol{\mathbf{a}}_{\textup{true}}$, where $\mu$ indicates the sample mean. The two horizontal lines indicate a FWHM distance from the centre of the beam. Figure 10: The absolute value of the fractional difference between the posterior mean and true sky taken as a slice through the centre of the primary beam region at dec=$-31.7\text{\,}\mathrm{\SIUnitSymbolDegree}$. Artefacts due to the pixelisation have been smoothed away by applying a pair of moving average filters to each of these curves. Results B4-B5 have different LST/RA ranges as indicated by the additional x-axes. Note that for B2, only the first half of the LST/RA range is shown, to match the $8\text{\,}\mathrm{h}$ span of the other runs. The light shaded region covers the FWHM of $12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ and the dark shaded region is outside of the primary beam. It is noticeable that all three ${\ell_{\text{max}}}=30$ results have higher fractional differences to the true than the ${\ell_{\text{max}}}=20$ results. Table 1: Overview of all the cases analysed in this paper and their labels along with which parameters are varied. A “—” indicates that the parameter is set as in the standard case. The labels are numbered after the order of appearance in Figs. 7 and 9.  | $N_{\rm samples}$ | Noise covariance | Prior covariance | Prior mean
---|---|---|---|---
A1 | $\sim 3500$ | — | — | —
A2 | — | — | — | $\boldsymbol{\mathbf{a}}_{0}=0$
A3 | — | $\boldsymbol{\mathbf{N}}\times 10^{4}$ | — | $\boldsymbol{\mathbf{a}}_{0}=0$
A4 | — | $\boldsymbol{\mathbf{N}}\times 10^{2}$ | — | —
A5 | — | $\boldsymbol{\mathbf{N}}\times 10^{2}$ | $\boldsymbol{\mathbf{S}}^{-1}=0$ | —
A6 | — | — | $\boldsymbol{\mathbf{S}}\times 10^{2}$ | —
A7 | — | $\boldsymbol{\mathbf{N}}\times 10^{4}$ | — | —
A8 | — | $\boldsymbol{\mathbf{N}}\times 10^{4}$ | $\boldsymbol{\mathbf{S}}^{-1}=0$ | —
A9 | — | — | $\boldsymbol{\mathbf{S}}\times 10^{4}$ | —
| ${\ell_{\text{max}}}$ | Number of LSTs | LST range |
B1 | ${\ell_{\text{max}}}=30$ | — | — |
B2 | — | $N_{\rm LST}=20$ | $0-$16\text{\,}\mathrm{h}$$ |
B3 | — | $N_{\rm LST}=20$ | $0-8\,\,\,\,$\mathrm{h}$$ |
B4 | — | $N_{\rm LST}=10$ | $16-$24\text{\,}\mathrm{h}$\,\,\,$ |
B5 | — | $N_{\rm LST}=10$ | $8-$16\text{\,}\mathrm{h}$$ |
| $N_{\rm ants}$ | Dish diameter | FWHM |
C1 | — | $3\text{\,}\mathrm{m}$ | $57.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ |

In this section we explore how varying the noise covariance $\mathbf{N}$,
prior covariance $\mathbf{S}$, and prior mean $\boldsymbol{\mathbf{a}}_{0}$
will affect the results from the GCR-sampler compared to our chosen standard
case. The motivation for this is to assure that a sufficient balance is
reached between the data- and prior-information, and that the prior is not
dominating the system. The results of the various noise- and prior level runs
can be seen in Figs. 7 and 8. All the various runs have been labelled in order
of there appearance in Fig. 7 and these labels will be used throughout the
paper. A full overview of all labels can also be seen in Table 1.

Before varying the noise and prior-levels, we examine whether using many more
samples ($\approx 3500$) would significantly change the results (case A1), in
order to check that only 100 samples is enough to produce statistically valid
results. The fractional difference between the posterior mean and true sky for
the A1 result is shown in Fig. 7 along with re-displaying the standard case
for comparison. It is clear that the fractional difference of the two cases
look almost identical. For even closer comparison we take a slice of the
centre of the beam of the absolute values of the fractional difference to the
true value, as shown in Fig. 8. As the figure shows, even when comparing the
fractional difference slices, the two cases only have a slight difference at
${\rm RA}\sim$10\text{\,}\mathrm{\SIUnitSymbolDegree}$$, thus emphasising that
these two results are almost identical, which (together with the full beam
region fractional difference) justifies the use of 100 samples as adequate for
our purposes.

Next, we set the inverse prior covariance to zero, which may be interpreted as
the case where only the likelihood term is used and the prior is neglected. We
then multiply the noise covariance with a factor of $10^{2}$ and $10^{4}$,
case A5 and A8 respectively. This is done in order to confirm, that removing
the information provided by the likelihood and/or prior does indeed lead to no
recovery of information of the sky, i.e. we demonstrate that the good recovery
in other cases is not simply due to the structure of the sky model or similar.
The fractional difference between the posterior mean and true sky for the A5
case ($\boldsymbol{\mathbf{S}}^{-1}=0,\mathbf{N}\times 10^{2}$) shows a slight
hint of the primary beam structure, but the high noise levels washes out the
data and there is no prior-information to assist the recovery of the sky. When
increasing the noise further in the A8 case
($\boldsymbol{\mathbf{S}}^{-1}=0,\mathbf{N}\times 10^{4}$), the recovered sky
is all gone and there is only noise left. In order to see how well the sampler
can recover the sky with no prior information, an extra run (not shown) was
made with standard noise level (but keeping the inverse prior covariance at
zero), which resulted in a fractional difference between the posterior mean
and true sky very similar to that of case A8 (i.e. where $\mathbf{S}\times
10^{4}$). Essentially the very high value of the prior covariance $\mathbf{S}$
is the same as saying the inverse prior covariance
$\boldsymbol{\mathbf{S}}^{-1}$ is closer to zero. It is therefore not that
surprising that these two results would be similar. For both these cases the
prior is now so broad and uncertain that the only well-recovered part is the
observed region.

The prior drives the good recovery outside of the observed region, which is
very clear when comparing for instance case A7 ($\mathbf{S}=$ standard,
$\mathbf{N}\times 10^{4}$) to case A8 ($\mathbf{S}^{-1}=0,\mathbf{N}\times
10^{4}$). Without any prior information nothing is recovered – but with the
prior most of the sky is recovered well to within $\sim 10\%$. However, A7
lacks the primary beam structure that we saw in the standard case, so the
recovery here is due to the prior driving the solution since the broad noise
covariance is washing out the data. When the noise is lowered again in case A4
to $\mathbf{N}\times 10^{2}$ (while keeping $\mathbf{S}=$ standard), the
central beam structure starts to reappear. The good recovery of the sky
outside of the observed region is without a doubt driven by the prior – as the
data cannot say anything about this region anyway – but the _very_ good
recovery of the observed region is due to the data. In the standard case the
well-recovered sky is further limited to be close to the observed region. The
fact that the spherical harmonics are continuous means the sky will still be
constrained by the data a little outside of the observed region. This explains
why the standard deviation in Fig. 4 smoothly goes from low to a high value.

An important test is to make sure that our choice of prior mean,
$\boldsymbol{\mathbf{a}}_{0}$, is not over-informing the sampler. Comparing
the A2 ($\mathbf{a}_{0}=0,\mathbf{N}=$ standard) and A3
($\mathbf{a}_{0}=0,\mathbf{N}\times 10^{4}$) cases to the standard case in
Fig. 7, it initially looks like setting the prior mean to zero results in a
poorly recovered sky. Indeed, this is the case if the noise is also increased
(as in A3), although generally increasing the noise with a factor of $10^{4}$
has yielded poor results (A7, A8, A9) and does not seem to be due to the
change of the prior mean alone. Instead, when the noise level is kept at
standard level (case A2), the true sky is still well-recovered near the centre
of the beam (i.e. recover is not biased), but the residual is larger beyond
the beam FWHM than in the standard case. Again, taking a slice of the absolute
fractional difference at the central declination only (around
$\textup{dec}=$-31.7\text{\,}\mathrm{\SIUnitSymbolDegree}$$), Fig. 8), shows
that the two $\mathbf{a}_{0}=0$ cases are not performing as good as the
standard case. Although, especially for the A2 case, at
RA$=60-$100\text{\,}\mathrm{\SIUnitSymbolDegree}$$ it is performing similarly
or (very) slightly better than the standard case. In most of our runs we have
a 10% prior that has prior means which are scattered around the true value.
This is somewhat conservative (we are not assuming a ’correct’ prior mean and
thereby not enforcing the true solution), but there is still the question of
how sensitive the results are to prior choices. What the
$\boldsymbol{\mathbf{a}}_{0}=0$ runs (A2, A3) show is what happens when our
prior is typically wrong by $10\sigma$.

###  4.3 Effect of changing array and observing configurations

The performance of the sampler is affected by which spherical harmonic modes
are available to it. This can be tested by changing the ${\ell_{\text{max}}}$
or for instance by including longer baselines to increase sensitivity to
higher $\ell$ modes, as can also be gauged from Fig. 6. Since most of the
diffuse foreground power comes from the Galaxy, it is also interesting to
examine how changing the LST range will affect how well the sky is recovered.
All the test cases in the following section has been labelled with a similar
style to the previous section, now based on the order of appearance in Fig. 9.
A full overview of all labels can be seen in Table 1.

To capture regions of the entire HERA-strip we examine LST ranges of
$8-$16\text{\,}\mathrm{h}$$ and $16-$24\text{\,}\mathrm{h}$$ on top of the
$0-$8\text{\,}\mathrm{h}$$ LST range; case B5, B4, and the standard case,
respectively. The B4 case covers the brightest regions of the galactic diffuse
emission, which can be seen in the RA range of
$$240\text{\,}\mathrm{\SIUnitSymbolDegree}$-$360\text{\,}\mathrm{\SIUnitSymbolDegree}$$
(the readers right hand side) on the top-left plot of Fig. 4. We also examine
the effects of increasing the data set to $N_{\textup{LST}}=20$ either by
doubling the LST range and keeping the same cadence (case B2) – or, vice versa
(case B3). The former results in a greater sky coverage and the latter results
in more information of the same sky region. The standard case cadence is one
integration of duration 60 sec per 48 minutes equivalent of 10 samples evenly
spread over an 8 hour LST range. Note, since the sky rotates
$15\text{\,}\deg\mathrm{/}\mathrm{h}$, the beam-crossing time of the HERA-like
array is $\sim$49\text{\,}\mathrm{min}$$ given the FWHM of
$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$. In this section we go over the
results from each of these cases and how they affect the resulting recovered
sky.

For the standard case we use ${\ell_{\text{max}}}=20$ because it is a
reasonable minimum for testing purposes that includes the shortest (most
sensitive) HERA baseline, see Fig. 6. In case B1, we see from Fig. 10 that
increasing ${\ell_{\text{max}}}$ to 30 results in a significantly degraded
fractional residual compared to ${\ell_{\text{max}}}$ of 20. While increasing
${\ell_{\text{max}}}$ brings the $b=$25.3\text{\,}\mathrm{m}$$ baseline group
into the _directly_ constrained $\ell$ range (see Fig. 6), each baseline is
actually sensitive to a broad range of $\ell$-values, as for example the
primary beam contributes to smearing of the $\ell$-mode sensitivity. For
instance the $b=$29.2\text{\,}\mathrm{m}$$ baseline (antennas $(0,2)$) in Fig.
1 is clearly sensitive to $\ell$-values down to around $\ell\sim 8$ and the
$(0,5)$-baseline, which is a shorter $25.3\text{\,}\mathrm{m}$ baseline,
contributes mostly at $\ell\gtrsim 5$ but still has some visibility response
below this level. These baselines will therefore already contribute to some
degree at the ${\ell_{\text{max}}}=20$ level. However, when we increase the
${\ell_{\text{max}}}$ to 30 the total number of modes more than doubles, from
441 to 961. This stretches the available signal-to-noise of the data much
farther than in the ${\ell_{\text{max}}}=20$ case.

Next, we vary the LST ranges. First, we keep to 10 time steps and 8 hr ranges
(as in the standard case), but shift the ranges to cover different parts of
the sky. Fig. 9 shows that the range of $8-$16\text{\,}\mathrm{h}$$ (B5)
performs better than both the standard $0-$8\text{\,}\mathrm{h}$$ and the B4
case ($16-$24\text{\,}\mathrm{h}$$), which covers the brightest part of the
galactic emission. As before, we inspect the slice at the central declination
for a closer comparison, see Fig. 10. It is clear from the figure, that within
their respective observed regions all three sets of LST ranges actually
perform comparatively well with only small differences. As one moves closer to
the edge of the observed region, however, the B5 case shows some improvement
to the standard case – and, especially, when compared to the B4 case. The
difference in recovery level can simply be due to the shifting of the sampling
points. For all the runs the simulated data is the same and have one specific
realisation of the noise. Areas that have less noise will therefore always
have a higher signal-to-noise and should perform better.

Increasing the data volume should also help constrain the sky better. One way
to do this is to double the number of observation times from
$N_{\textup{LST}}=10$ to $20$, either by doubling the range and keeping the
cadence the same (case B2) or by increasing the cadence within the same LST
range (case B3). At first glance both results seem to perform similar or
perhaps slightly better than the standard case, see Fig. 9. The recovered
signal will be worse at the edges of the observed region. The B2 case is twice
the size and encompasses the far edge of the B3 case at LST
$=$8\text{\,}\mathrm{h}$$, and will naturally be better recovered at this LST
than the B3 case. We will therefore limit the comparison to the
$0-$8\text{\,}\mathrm{h}$$ range for the fractional residual slice in Fig. 10.
The B3 result mostly shows improvement to the standard case with the exception
of RA
$$20\text{\,}\mathrm{\SIUnitSymbolDegree}$-$40\text{\,}\mathrm{\SIUnitSymbolDegree}$$
although exceeded by the performance of the B2 case, that not only performs
better than the standard case but also remains closer to true at the lower
edge ($0\text{\,}\mathrm{h}$ of the observed region. This hints that
increasing the observation time range (and thus increasing the size of the
observable sky) has a larger impact on the recovery of the sky than doubling
the data we have inside the same region. This can also be explained with a
larger observable sky being sensitive to larger scales (lower spherical
harmonic modes).

###  4.4 Improving sensitivity to larger scales

Figure 11: The real- and imaginary parts of the difference of the mean of 100
samples of $a_{\ell m}$ modes from the GCR solver and the true sky with the
diameter of the receivers set to $\theta_{D}=$3\text{\,}\mathrm{m}$$ and
thereby increasing the field of view. As described in Sec. 2.1 the $m=0$
imaginary modes ($\times$) should always be zero and the GCR solver therefore
does not solve for this. All residual ${a_{\ell m}}$ modes now lie within the
$\pm 100$ boundaries and including the uncertainties are consistent with zero,
except the zeroth $\ell$-mode that still has some loss of power.

Figure 12: Resulting maps from the estimated spherical harmonic modes on the
sky with the diameter of the receivers set to
$\theta_{D}=$3\text{\,}\mathrm{m}$$ and thereby increasing the field of view.
_Upper left:_ The true sky given by pyGSM with ${\ell_{\text{max}}}=20$ and
nside$=128$. Note that the rippled structure comes from the true sky and not
from the GCR solver. _Upper right:_ The estimated sky based on the mean of 100
samples from the GCR solver. The spherical harmonic modes can be seen in Fig.
11. _Bottom left:_ Fractional difference between the mean and the true sky
adjusted to show differences $<10\%$ (i.e. highlighting improvements over the
prior). The beam region is no longer clear cut as in the standard case.
_Bottom right:_ The standard deviation of the difference of the estimated and
true sky. The lowest noise levels are now larger than with the standard case,
but overall noise is lower.

One reason that the smaller scales are more difficult to constrain comes from
the array layout itself, as discussed above. The prior helps breaking the
degeneracy of the low-$\ell$ modes, that can otherwise only be determined as a
linear combination (i.e. they are degenerate, as there are no sufficiently
short baselines to resolve them individually). It is clear, however, from the
$a_{\ell m}$ modes presented in Fig. 3 that these modes are still difficult to
constrain up to $\ell\sim 5$. Moving on the same logic as above, where the
observable sky is increased to increase sensitivity to the larger scales, we
can decrease the size of the dishes to increase the field of view (FoV). With
the HERA-like dishes the FWHM $=$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$$.
Altering the array to have a FoV more similar to that of OVRO-LWA (Eastwood et
al., 2018) or MWA (Yoshiura et al., 2021), the diameter of the receivers is
changed to $3\text{\,}\mathrm{m}$, now with a ${\rm
FWHM}=$57.3\text{\,}\mathrm{\SIUnitSymbolDegree}$$ (case C1). A broader FoV
now means that the smearing around the $\ell$ modes in Fig. 6 is much narrower
($\delta\ell\simeq\pi/{\rm FWHM}$), thus more clearly picking out the specific
$\ell$ modes. Unlike when we increase the LST range, we expand the observable
sky not only in the E-W direction but also N-S. As increasing the FoV with
more than $4\times{\rm FWHM}$ of the standard case is a much more radical
change than those made in Sec. 4.3, the C1 results are presented in more
detail in Figs. 11 and 12, although now leaving out the Cartesian projections
of the fractional residual, since the beam is now so wide that a close up of
the primary beam region is no longer relevant. Starting with the true-
subtracted mean of the recovered ${a_{\ell m}}$ modes in Fig. 11, it is
noticeable that there are no longer any _outliers_ from the displayed region,
as was the case in Fig. 3 with the standard configuration. The zeroth-mode is
still the most difficult to constrain and is again underestimated, however,
not as severely underestimated as earlier. Earlier it was argued that the
${a_{\ell m}}$ modes are not the most representative performance indicator of
the code, since they are all-sky quantities and the observed region is only a
narrow strip, but now that we have much larger coverage it clearly shows that
full-sky features as the ${a_{\ell m}}$ modes can be well-constrained.

By-eye, it is now almost impossible to tell the difference between the true
sky map and the mean of the recovered sky in Fig. 12. Taking the fractional
difference between the two maps reveals a very different structure than
earlier. Now, almost all parts of the sky are recovered within $10\%$ with the
Galaxy showing the smallest residuals. It is striking that the primary beam
structure no longer shows up in either the fractional difference map or the
standard deviation. Since the noise covariance $\mathbf{N}$ and subsequently
the noise $n_{ij}$ on the data is defined by the auto-visibilities as per Eq.
25, the noise level will be affected by whether the Galactic emission is
within the beam. Since the FoV is now so large, the noise has increased with
the extra power from the Galaxy as most of this falls within the observed
area.

Additional runs have been made where the noise was reduced with a factor of
$10^{-1}$ and $10^{-2}$ respectively (not shown), and not only does this
improve the estimate of the ${a_{\ell m}}$ modes, it is also clear that the
area of the sky with the largest fractional residual is the furthest region
from the observed region, corresponding to the same area in Fig. 12 where the
standard deviation is higher. Decreasing the noise also reduced the standard
deviation inside the observed region to the same level as the standard case.
Earlier, in Sec. 4.2, we showed that reducing the noise level could lead to
worse recovery of the sky. Since we are dealing with constrained realisations,
the solution in the prior-dominated region can still depend on the data-
dominated region (and vice-versa), as the spherical harmonic basis functions
connect the two. The influence of the data on the prior-dominated region is
minimal if the noise covariance is very narrow. For the $3\text{\,}\mathrm{m}$
beam case it is not as big an issue, since the large FoV makes up for it.
Before the beam was also very narrow, so when the noise covariance was made
narrow, it did not allow the prior to contribute as much, which made it
difficult to recover the sky even at the edges of the observed region. Now,
that the FoV is larger, the prior is not as crucial to determine full-sky
features and we can reduce the noise without losing the constraining power.

##  5 Conclusions

For next-generation wide-field radio interferometer arrays, particularly those
targeting the 21 cm fluctuations at low and high redshift, the biggest
challenge remains proper handling of foreground contamination following its
distortion/modulation by the instrument. To improve on this issue, it is
crucial to be able to make accurate maps of the radio emission on both small
and large scales from the measured visibilities. Traditional deconvolution
algorithms such as CLEAN (and its extensions) can struggle to properly recover
diffuse emission, which is the dominant foreground contaminant on large
angular scales. To address this, methods such as Direct Optimal Mapping and
$m$-mode analysis have been developed, both of which are maximum likelihood
estimators that focus on accurately recovering wide field maps of the
emission, i.e. on scales where the curvature of the sky and the shape of the
primary beam becomes important.

In this work we have presented an alternative approach to recovering diffuse
foregrounds from visibilities – by constructing a parametric model of the
emission, represented by spherical harmonics on the sky, and drawing samples
from the joint posterior of the spherical harmonic coefficients given the
visibility data and a choice of prior. A linear model can be constructed by
writing down an operator that encodes the response of the interferometer (and
thus the measured visibilities) to each spherical harmonic mode. This operator
includes all the relevant instrumental degrees of freedom. We can then
estimate the joint posterior distribution of all the coefficients by
repeatedly solving the Gaussian constrained realisation (GCR) equation, an
extension of the Wiener filter that includes random realisation terms. This
can be solved efficiently even for a very large number of coefficients, making
this large inference problem computationally tractable. Furthermore, each
sample of coefficients is a complete realisation of the spherical harmonic
coefficients, and therefore the full sky (i.e. with no gaps or mask) that is
statistically-consistent with the data. This ultimately allows us to recover
the diffuse emission signal in a statistically robust manner while also
avoiding difficulties with missing data in subsequent steps of the data
analysis pipeline.

After presenting the mathematical formalism for this sampler, the primary aim
of this paper was to validate the method by applying it to a suite of
simulations. We tested the performance of the spherical harmonic sampler by
comparing the recovered sky, defined as the mean of the sky maps of the
samples, to the (known) true sky modelled using pyGSM. For the analysis we
chose a standard set of parameters to act as a reference for the various tests
of noise and prior uncertainty levels, sample size, ${\ell_{\text{max}}}$, the
specific LST range of the observations, and the influence of the field of
view. The _standard case_ is based on 100 realisations with a $10\%$ prior
centred on the true ${a_{\ell m}}$-modes, and the noise is given by a Gaussian
distribution with zero mean and covariance defined by the auto-visibilities
through the radiometer equation. The standard case uses a HERA-like closed-
packed redundant array of 10 antennas and covers the LST range
$0-$8\text{\,}\mathrm{h}$$. The dishes are 14 m in diameter (FWHM
$=$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$$), which means the directly
observed sky is a narrow $24.6\text{\,}\mathrm{\SIUnitSymbolDegree}$ (dec) by
$120\text{\,}\mathrm{\SIUnitSymbolDegree}$ (RA) strip covering only $\sim
8.6\%$ of the full sky.

The standard case performs well within the observed region and recovers the
true sky to within $10\%$ within a distance of 1 FWHM from the centre and to
within $5\%$ at the centre of the beam. When increasing the number of samples
to $N_{\textup{samples}}\approx 3500$, we found that there were minimal
changes to the result allowing us to keep the computational costs lower by
continuing with 100 samples per scenario only. Next, the impact of the noise
level on recovery was studied. When the noise variance is increased by a
factor of $>10^{2}$, the noise will dominate over the sky signal and downgrade
the recovery. At this noise level, the samples will only show diffuse sky
structure if there is a prior as well to drive the solution. For the standard
case we chose a $10\%$ prior, as the pyGSM model should describe the diffuse
emission sky to roughly within this margin. If the prior covariance is
increased, however, the prior quickly becomes too broad and uncertain and the
only well-recovered part is the very centre of the observed region. Including
a prior helps break the degeneracy of the low-$\ell$ modes. It was also found
that even if the prior mean is set to zero (a highly pessimistic assumption),
the sky can still be recovered, although to a slightly worse degree than the
standard case. This shows the benefit of including a prior; the boundary
region between observed and unobserved sky is constrained by both the prior
and the data, while regions of missing data are still filled with a
statistically plausible realisation. We conclude that, for the standard case,
the prior is not overly-informative, and contributes about the same as the
data to the diffuse emission recovery around the observed region.

For close-packed arrays like HERA, the most abundant baselines are the
shortest ones, which have lengths $<$30\text{\,}\mathrm{m}$$. This justifies
choosing a relatively low ${\ell_{\text{max}}}$ for our tests. When increasing
the ${\ell_{\text{max}}}$ from 20 to 30, the number of modes in the ${a_{\ell
m}}$-vector increases significantly from 441 to 961, more than doubling the
number of parameters we are trying to constrain, and thus stretching the
available signal-to-noise across more degrees of freedom. A more informative
prior would help to balance this increase in the size of the parameter space.

Most of the diffuse emission on the sky originates from the Galaxy and is not
uniformly distributed. This results in field-dependent effects on the recovery
of the signal. For instance, the recovery was best for a
$8-$16\text{\,}\mathrm{h}$$ region, whereas the region covering the brightest
part of the Galaxy ($16-24$\mathrm{h}$$) was recovered worst, due largely to
the extra noise given by the auto-visibilities. When increasing the data
volume we expect the results to improve, for instance if we increase the
cadence and take more data samples within the same LST range, or, if we keep
the cadence but increase the range. When doubling the number of data points
within the range, there was only a slight improvement to the recovered sky,
suggesting that there is not much information to be gained by making more
measurements of nearby (correlated) pointings, e.g. within a beam width of one
another. In contrast, doubling the range itself not only improved the results
in the centre, but also at the boundaries of the LST range, as additional sky
coverage led to improved constraints on the spherical harmonic modes across
the board.

Ultimately, the biggest impact on our simulated results occurred by increasing
the field of view. This was done by decreasing the effective antenna diameter
to $3\text{\,}\mathrm{m}$, thus increasing the FWHM from
$12.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ to
$57.3\text{\,}\mathrm{\SIUnitSymbolDegree}$. For the large-FoV case, the
spherical harmonic coefficients now all coincide with the true sky within
their given uncertainties except the zeroth mode. The standard deviation map
of the large-FoV results no longer shows the clear strip of the observed sky,
as almost all of the modes are now well-constrained. Since the noise is
defined by the auto-visibilities, and this setup now directly observes most of
the Galaxy, the noise level is larger than the standard case. However, as the
FoV is now so large, we are also less dependent on the prior as there are
fewer gaps of information to fill in, and the boundary region (between
directly observed and totally unobserved regions) is wider.

The results presented in this paper are based on a handful of simplified
example cases, using a reduced-size array and only two frequencies for the
visibility response operator. A simple extension of this work would be to
rerun the analysis while simulating the visibility response for more
frequencies or, alternatively, the frequency dependence could be included via
a parametric model, e.g. power law with spectral index $\beta$. Likewise, we
have only considered cases of ${\ell_{\text{max}}}<30$, but in reality a
higher angular resolution would be desirable if a larger portion of the array
is considered (i.e. including more longer baselines). Finally, the Gaussian
constrained realisation method has been developed with inclusion into a Gibbs
sampling framework in mind, whereby beam, 21cm signal, and point source
foreground model parameters would also be sampled. A simpler and more focused
Gibbs sampling scheme would also enable us to sample the signal covariance
$\boldsymbol{\mathbf{S}}$ (e.g. as parametrised by the angular power
spectrum), as well as the spherical harmonic coefficients themselves.

## Acknowledgements

This result is part of a project that has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 948764; KAG, PB, JB, MJW). We
acknowledge use of the following software: matplotlib (Hunter, 2007), numpy
(van der Walt et al., 2011), and scipy (Virtanen et al., 2020). This work used
the DiRAC@Durham facility managed by the Institute for Computational Cosmology
on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was
funded by BEIS capital funding via STFC capital grants ST/P002293/1,
ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant
ST/R000832/1. DiRAC is part of the National e-Infrastructure.

## Data Availability

The Python code used to produce the results in this paper is available from
https://github.com/HydraRadio/Hydra. The Jupyter notebook used for the data
analysis and to produce plots is available from
https://github.com/katrinealice/sph_harm_GCR_analysis. Simulated data are
available on request.

## References

  * Ali et al. (2015) Ali Z. S., et al., 2015, ApJ, 809, 61
  * Bowman et al. (2013) Bowman J. D., et al., 2013, PASA, 30, e031
  * CHIME Collaboration (2023) CHIME Collaboration 2023, The Astrophysical Journal, 947, 16
  * Cornwell (2008) Cornwell T. J., 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793–801
  * Cornwell, T. J. (2009) Cornwell, T. J. 2009, A&A, 500, 65
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Dillon & Parsons (2016) Dillon J. S., Parsons A. R., 2016, The Astrophysical Journal, 826, 181
  * Dillon et al. (2015) Dillon J. S., et al., 2015, Phys. Rev. D, 91, 023002
  * Dong (1994) Dong R., 1994, Nuclear Magnetic Resonance of Liquid Crystals. Partially ordered systems, Springer-Verlag, https://books.google.co.uk/books?id=tAbwAAAAMAAJ
  * Dowell & Taylor (2018) Dowell J., Taylor G. B., 2018, The Astrophysical Journal Letters, 858, L9
  * Eastwood et al. (2018) Eastwood M. W., et al., 2018, The Astronomical Journal, 156, 32
  * Eriksen et al. (2008) Eriksen H. K., Jewell J. B., Dickinson C., Banday A. J., Górski K. M., Lawrence C. R., 2008, ApJ, 676, 10
  * Evoli et al. (2014) Evoli C., Mesinger A., Ferrara A., 2014, J. Cosmology Astropart. Phys, 2014, 024
  * Fixsen et al. (2011) Fixsen D. J., et al., 2011, The Astrophysical Journal, 734, 5
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Physics Reports, 433, 181
  * Geman & Geman (1984) Geman S., Geman D., 1984, IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6, 721
  * Högbom (1974) Högbom J. A., 1974, A&AS, 15, 417
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Kennedy et al. (2023) Kennedy F., Bull P., Wilensky M. J., Burba J., Choudhuri S., 2023, The Astrophysical Journal Supplement Series, 266, 23
  * Kittiwisit et al. (2023) Kittiwisit P., et al., 2023, arXiv e-prints, p. arXiv:2312.09763
  * Kriele (2022) Kriele M. A., 2022, PhD thesis, Eindhoven University of Technology, Curtin University, ASTRO3D
  * Liu & Shaw (2020) Liu A., Shaw J. R., 2020, Publications of the Astronomical Society of the Pacific, 132, 062001
  * Lopez-Honorez et al. (2016) Lopez-Honorez L., Mena O., Moliné Á., Palomares-Ruiz S., Vincent A. C., 2016, J. Cosmology Astropart. Phys, 2016, 004
  * Morales & Wyithe (2010) Morales M. F., Wyithe J. S. B., 2010, ARA&A, 48, 127
  * Offringa & Smirnov (2017) Offringa A. R., Smirnov O., 2017, Monthly Notices of the Royal Astronomical Society, 471, 301–316
  * Parsons et al. (2012) Parsons A., Pober J., McQuinn M., Jacobs D., Aguirre J., 2012, The Astrophysical Journal, 753, 81
  * Patil et al. (2017) Patil A. H., et al., 2017, The Astrophysical Journal, 838, 65
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A1
  * Price (2016) Price D. C., 2016, 2.0.0, Astrophysics Source Code Library, 1603.013
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Santos et al. (2005) Santos M. G., Cooray A., Knox L., 2005, ApJ, 625, 575
  * Shaw et al. (2014) Shaw J. R., Sigurdson K., Pen U.-L., Stebbins A., Sitwell M., 2014, ApJ, 781, 57
  * Shaw et al. (2015) Shaw J. R., Sigurdson K., Sitwell M., Stebbins A., Pen U.-L., 2015, Phys. Rev. D, 91, 083514
  * Singal et al. (2010) Singal J., Stawarz L., Lawrence A., Petrosian V., 2010, Monthly Notices of the Royal Astronomical Society, 409, 1172–1182
  * The HERA Collaboration et al. (2022) The HERA Collaboration et al., 2022, The Astrophysical Journal, 924, 51
  * Tingay et al. (2013) Tingay S. J., et al., 2013, Publications of the Astronomical Society of Australia, 30, e007
  * Tung (1985) Tung W., 1985, Group Theory in Physics. World Scientific, https://books.google.co.uk/books?id=O89tgpOBO04C
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wayth et al. (2018) Wayth R. B., et al., 2018, Publications of the Astronomical Society of Australia, 35, e033
  * Xu et al. (2022) Xu Z., et al., 2022, The Astrophysical Journal, 938, 128
  * Yoshiura et al. (2021) Yoshiura S., et al., 2021, Monthly Notices of the Royal Astronomical Society, 505, 4775–4790
  * Zheng et al. (2017) Zheng H., et al., 2017, MNRAS, 464, 3486
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, A&A, 556, A2
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22

##  Appendix A Further “realification”

Inherently visibilities are complex quantities and the spherical harmonic
coefficients are as well. Many standard linear solvers are not set up to
handle complex numbers however. The first step to make the system real (i.e.
_realify_) has already been done in Sec. 2.1 bt reordering the spherical
harmonic coefficient vector, $\boldsymbol{\mathbf{a}}_{\ell m}$, to contain
the values of first the real- and then the imaginary part of the spherical
harmonic coefficients. This leaves us with a fully real-valued vector to solve
for, denoted $\boldsymbol{\mathbf{a}}_{\rm cr}$ (constrained realisations). A
real-valued spherical harmonic coefficient vector does not automatically
result in real-valued visibilities however. Instead, the visibility response
is still complex valued, now containing a complex visibility response per
real- and per imaginary part of the $\boldsymbol{\mathbf{a}}_{\rm cr}$-vector.
Further _realification_ of the system is therefore needed, while making sure
that all the mixing of real and imaginary parts that you get when multiplying
complex numbers together is handled correctly.

First, we define vectors (bold, lower case) and matrices (bold, upper case) in
this real-valued system as

$\displaystyle\widetilde{\boldsymbol{\mathbf{v}}}=\begin{pmatrix}\boldsymbol{\mathbf{v}}_{\rm
re}\\\ \boldsymbol{\mathbf{v}}_{\rm
im}\end{pmatrix},\quad\widetilde{\boldsymbol{\mathbf{M}}}=\begin{pmatrix}\boldsymbol{\mathbf{M}}_{\rm
re}&-\boldsymbol{\mathbf{M}}_{\rm im}\\\\[3.0pt] \boldsymbol{\mathbf{M}}_{\rm
im}&\boldsymbol{\mathbf{M}}_{\rm re}\end{pmatrix}$ (26)

with the Hermitian conjugate of the matrix
$\widetilde{\boldsymbol{\mathbf{M}}}$ given as

$\displaystyle\widetilde{\boldsymbol{\mathbf{M}}}^{{\dagger}}=\begin{pmatrix}\boldsymbol{\mathbf{M}}_{\rm
re}^{T}&\boldsymbol{\mathbf{M}}_{\rm im}^{T}\\\\[3.0pt]
-\boldsymbol{\mathbf{M}}_{\rm im}^{T}&\boldsymbol{\mathbf{M}}_{\rm
re}^{T}\end{pmatrix},$ (27)

where the superscript $T$ denotes the transpose and the minus sign has been
swapped due to complex conjugation.

With this we can redefine the complex-valued visibility response
$\boldsymbol{\mathbf{X}}$ into a purely real-valued version,

$\displaystyle\widetilde{\boldsymbol{\mathbf{X}}}=\begin{pmatrix}\boldsymbol{\mathbf{X}}_{\textup{re}}&-\boldsymbol{\mathbf{X}}_{\textup{im}}\\\\[3.0pt]
\boldsymbol{\mathbf{X}}_{\textup{im}}&\boldsymbol{\mathbf{X}}_{\textup{re}}\end{pmatrix},\quad\widetilde{\boldsymbol{\mathbf{X}}}^{\dagger}=\begin{pmatrix}\boldsymbol{\mathbf{X}}_{\textup{re}}^{T}&\boldsymbol{\mathbf{X}}_{\textup{im}}^{T}\\\\[3.0pt]
-\boldsymbol{\mathbf{X}}_{\textup{im}}^{T}&\boldsymbol{\mathbf{X}}_{\textup{re}}^{T}\end{pmatrix}\,,$
(28)

and we define the new noise covariance as a diagonal matrix

$\displaystyle\widetilde{\boldsymbol{\mathbf{N}}}=\begin{pmatrix}\boldsymbol{\mathbf{N}}/2&0\\\
0&\boldsymbol{\mathbf{N}}/2\end{pmatrix}\,.$ (29)

Since the spherical harmonic coefficient vector is already constructed so that
it is real-valued, as described in Sec. 2.1, so is the prior mean
$\boldsymbol{\mathbf{a}}_{0}$ and $\boldsymbol{\mathbf{\omega}}_{a}$,

$\displaystyle\widetilde{\boldsymbol{\mathbf{a}}}=\begin{pmatrix}\boldsymbol{\mathbf{a}}_{\rm
cr}\\\
0\end{pmatrix},\quad\widetilde{\boldsymbol{\mathbf{a}}}_{0}=\begin{pmatrix}\boldsymbol{\mathbf{a}}_{0}\\\
0\end{pmatrix},\quad\widetilde{\boldsymbol{\mathbf{\omega}}}_{a}=\begin{pmatrix}\boldsymbol{\mathbf{\omega}}_{a}\\\
0\end{pmatrix}\,.$ (30)

Lastly, the signal covariance is defined in Sec. 3.3 by the real-valued
$\boldsymbol{\mathbf{a}}_{\text{true}}$-vector and so, already, is real-valued
itself. We therefore define the new signal covariance as,

$\displaystyle\widetilde{\boldsymbol{\mathbf{S}}}=\begin{pmatrix}{\color[rgb]{0,0,0}\boldsymbol{\mathbf{S}}}&0\\\
0&0\end{pmatrix}\,.$ (31)

Using these definitions, we can derive the Gaussian constrained realisation
equation once again, and we are left with the final _realified_ GCR equation,

$\displaystyle\left[\boldsymbol{\mathbf{S}}^{-1}+2\boldsymbol{\mathbf{X}}_{\textup{re}}^{T}\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}_{\textup{re}}+2\boldsymbol{\mathbf{X}}_{\textup{im}}^{T}\right.$
$\displaystyle\left.\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{X}}_{\textup{im}}\right]\boldsymbol{\mathbf{a}}_{\textup{cr}}=$
$\displaystyle\quad\boldsymbol{\mathbf{X}}_{\textup{re}}^{T}\left(2\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}_{\textup{re}}+\sqrt{2}\boldsymbol{\mathbf{N}}^{-\frac{1}{2}}(\boldsymbol{\mathbf{\omega}}_{d})_{\textup{re}}\right)$
$\displaystyle\quad+\boldsymbol{\mathbf{X}}_{\textup{im}}\left(2\boldsymbol{\mathbf{N}}^{-1}\boldsymbol{\mathbf{d}}_{\textup{im}}+\sqrt{2}\boldsymbol{\mathbf{N}}^{-\frac{1}{2}}(\boldsymbol{\mathbf{\omega}}_{d})_{\textup{im}}\right)$
$\displaystyle\quad+\boldsymbol{\mathbf{S}}^{-1}\boldsymbol{\mathbf{a}}_{0}+\boldsymbol{\mathbf{S}}^{-\frac{1}{2}}\boldsymbol{\mathbf{\omega}}_{a}.$
(32)