11institutetext: Kapteyn Astronomical Institute, University of Groningen, PO
Box 800, 9700AV Groningen, The Netherlands. 22institutetext: ASTRON,
Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD,
Dwingeloo, The Netherlands. 33institutetext: Department of Physics, Technion,
Haifa 32000, Israel 44institutetext: Institute for Computational Science,
University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland.
55institutetext: Anton Pannekoek Institute, University of Amsterdam, Postbus
94249 1090 GE, Amsterdam, The Netherlands. 66institutetext: LERMA,
Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université,
F-75014 Paris, France. 77institutetext: The Oskar Klein Centre, Department of
Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden.
88institutetext: Leiden Observatory, Leiden University, PO Box 9513, NL-2300
RA Leiden, The Netherlands.

# Degree-Scale Galactic Radio Emission at 122 MHz around the North Celestial
Pole with LOFAR-AARTFAAC

B. K. Gehlot<EMAIL_ADDRESS>L. V. E. Koopmans 11 A. R. Offringa
2211 H. Gan 11 R. Ghara S. K. Giri M. Kuiack 334455 F. G. Mertens 6611 M.
Mevius 22 R. Mondal V. N. Pandey 772211 A. Shulevski R. A. M. J. Wijers
88552255 and S. Yatawatta 22

(Received —; accepted —)

###### Abstract

Aims. Contamination from bright diffuse Galactic thermal and non-thermal radio
emission poses crucial challenges in experiments aiming to measure the 21-cm
signal of neutral hydrogen from the Cosmic Dawn (CD) and Epoch of Reionization
(EoR). If not included in calibration, this diffuse emission can severely
impact the analysis and signal extraction in 21-cm experiments. We examine
large-scale diffuse Galactic emission at 122 MHz, around the North Celestial
Pole, using the Amsterdam-ASTRON Radio Transient Facility and Analysis Centre
(AARTFAAC)- High Band Antenna (HBA) system.

Methods. In this pilot project, we present the first-ever wide-field image
produced with a single sub-band of the data recorded with the AARTFAAC-HBA
system. We demonstrate two methods: multiscale CLEAN and shapelet
decomposition, to model the diffuse emission revealed in the image. We use
angular power spectrum metrics to quantify different components of the
emission and compare the performance of the two diffuse structure modelling
approaches.

Results. We observe that the point sources dominate the angular power spectrum
($\ell(\ell+1)C_{\ell}/2\piup\equiv\Delta^{2}(\ell)$) of the emission in the
field on scales $\ell\gtrsim 60$ ($\lesssim 3$ degree). The angular power
spectrum after subtraction of compact sources is flat within
$20\lesssim\ell\lesssim 200$ range, suggesting that the residual power is
dominated by the diffuse emission on scales $\ell\lesssim 200$. The residual
diffuse emission has a brightness temperature variance of
$\Delta^{2}_{\ell=180}=(145.64\pm 13.61)\leavevmode\nobreak\ {\rm K}^{2}$ at
122 MHz on angular scales of 1 degree, and is consistent with a power-law
following $C_{\ell}\propto\ell^{-2.0}$ in $20\lesssim\ell\lesssim 200$ range.
We also find that, in the current setup, the multiscale CLEAN is suitable to
model the compact and diffuse structures on a wide range of angular scales,
whereas the shapelet decomposition method better models the large scales,
which are of the order of a few degrees and wider.

###### Key Words.:

dark ages, reionization, first stars – Methods: observational – Methods:
statistical – Techniques: interferometric – Radio continuum: general

##  1 Introduction

Observations of the redshifted 21-cm hyperfine transition line of neutral
hydrogen (Hi hereafter) are expected to unveil the properties of the first
luminous objects and the evolution of the large-scale structure during the
Cosmic Dawn (CD; $12\lesssim z\lesssim 30$) and Epoch of Reionization (EoR;
$6\lesssim z\lesssim 12$) (Madau et al., 1997; Shaver et al., 1999; Furlanetto
et al., 2006). Several current and next-generation experiments such as GMRT
(Paciga et al., 2013), LOFAR (van Haarlem et al., 2013), MWA (Tingay et al.,
2013; Bowman et al., 2013), NenuFAR (Zarka et al., 2012; Mertens et al.,
2021), OVRO-LWA (Eastwood et al., 2019), AARTFAAC (Gehlot et al., 2020), HERA
(DeBoer et al., 2017), SKA (Koopmans et al., 2015), EDGES (Bowman et al.,
2018), SARAS2 (Singh et al., 2017), LEDA (Bernardi et al., 2016), and REACH
(de Lera Acedo, 2019) are attempting to measure the fluctuations and global
evolution of the 21-cm brightness temperature as a function of redshift.

The faint 21-cm signal from high redshifts, however, is contaminated by
astrophysical foreground emission that is 3-4 orders of magnitude brighter
than the signal of interest. The astrophysical foregrounds consist of Galactic
diffuse emission (synchrotron and free-free emission), and extra-galactic
compact sources such as radio-galaxies, supernova remnants and other sources
(Di Matteo et al., 2002; Zaldarriaga et al., 2004; Bernardi et al., 2009,
2010; Ghosh et al., 2012). Subtraction or avoidance of these bright
foregrounds poses a significant challenge in all 21-cm cosmology experiments.
The interferometric experiments such as PAPER, MWA, and HERA follow the so-
called “foreground-avoidance” technique for the separation of the 21-cm signal
from the foregrounds. This technique takes advantage of the fact that the
spectrally smooth foregrounds occupy a small number of spectral modes and tend
to reside within a “wedge-shaped” region of the Fourier space, whereas the
“EoR-window”, in contrast to the “wedge” region, should be dominated by 21-cm
signal and the noise (Datta et al., 2010; Morales et al., 2012; Vedantham et
al., 2012; Trott et al., 2012; Parsons et al., 2012; Hazelton et al., 2013;
Dillon et al., 2014). On the other hand, 21-cm experiments with LOFAR,
AARTFAAC, and NenuFAR (and also the planned SKA) follow the foreground-removal
approach in their analysis strategy (Gehlot et al., 2019, 2020; Mertens et
al., 2020). Various foreground removal techniques used in the latter
experiments currently include image deconvolution and source extraction
algorithms such as CLEAN (Högbom, 1974; Clark, 1980; Cornwell, 2008; Offringa
& Smirnov, 2017), Duchamp (Whiting, 2012), and PyBDSF (Mohan & Rafferty, 2015)
to model the extra-galactic compact sources as delta functions and Gaussians
in image space. This is followed by subtraction of the model from the observed
signal using Direction Dependent (DD) calibration, e.g., with sagecal
(Yatawatta, 2016; Yatawatta et al., 2017). However, these traditional methods,
largely developed for longer-baseline high brightness-temperature data, are
sub-optimal for modelling large-scale diffuse emission, extended sources and
foreground emission below the confusion noise that dominates the observed
signal. The latter are typically mitigated in these experiments using blind
foreground removal techniques such as Wp Smoothing (Harker et al., 2009),
FastICA (Chapman et al., 2012), GMCA (Chapman et al., 2013), and GPR (Mertens
et al., 2020) that exploit the spectral smoothness of the foregrounds to model
and remove foregrounds from the observed data.

Besides the intrinsic interest of studying properties of the diffuse Galactic
emission, it is essential to include the diffuse (and compact) emission
component in the instrumental gain calibration step regardless of the signal
separation choices adopted by various experiments. Various investigations by
Patil et al. (2016), Barry et al. (2016), Ewall-Wice et al. (2017), and Mouri
Sardarabadi & Koopmans (2019) showed that using an incomplete sky-model in the
calibration step leads to significant suppression of the diffuse emission as
well as the 21-cm signal of interest on shorter baselines that are dominated
by diffuse emission. Since shorter baselines provide most of the sensitivity
towards the 21-cm signal, it is crucial to calibrate these baselines
accurately to prevent any signal suppression. One option to mitigate the
effect of sky-incompleteness is to remove the short baselines from calibration
and include only longer baselines where the diffuse emission is resolved out
(Patil et al., 2017). However, the exclusion of short baselines from
calibration leads to a so-called excess variance on the excluded baselines in
calibration (Patil et al., 2016; Barry et al., 2016; Ewall-Wice et al., 2017;
Mouri Sardarabadi & Koopmans, 2019), impacting the power spectrum estimation.
Data from several 21-cm experiments such as MWA (Byrne et al., 2021), OVRO-LWA
(Eastwood et al., 2018) and LWA New Mexico station (Dowell et al., 2017) are
already being used to develop diffuse sky maps to facilitate precise
calibration of 21-cm observations. As a part of the ACE project (Gehlot et
al., 2020), the data from AARTFAAC-LBA is also being utilized to produce a
wide-band map of northern sky ($\sim$10 arcmin resolution) at low frequencies
with high fidelity due to its superb $uv$-coverage.

In this work, we present the first-ever wide-field images produced using the
High-Band Antenna (HBA) observations of the LOFAR Amsterdam-ASTRON Radio
Transients Facility And Analysis Centre (AARTFAAC) wide-field imager (Prasad
et al., 2016) at 122 MHz. We explore multiscale CLEAN and shapelet
decomposition methods to model the diffuse emission in an extended region of
around $20^{\circ}$ radius centred at the North Celestial Pole (NCP) and
quantify this diffuse emission using the angular power spectrum technique.
Line et al. (2020) also presented a similar study comparing multiscale CLEAN
and shapelet decomposition to model the bright extended source Fornax A using
simulations and data from MWA. However, that study was focused on modelling
sources with extended emission from a few arcmins up to a degree.

The paper is organized as follows: section 2 and 3 briefly describe the
observational setup, preprocessing steps, calibration, and imaging scheme. The
two diffuse foreground modelling methods and their comparison are presented in
section 4. Section 5 describes the angular power spectrum to quantify the
properties of the diffuse emission. Finally, section 6 provides a summary of
the work and lays out the future steps and improvements.

Table 1: Observational details of the data. Parameter | value
---|---
Telescope | LOFAR-AARTFAAC
Antenna configuration | A12
Number of receivers | 576 (HBA tiles)
Number of pointings | 1
Pointing centre (RA, Dec; J2000): | 0h0m0s, $+90^{\circ}00^{\prime}00^{\prime\prime}$
Duration of observation | 11 hours
Observing frequency (and width) | 122 MHz (195.3 kHz)
Primary Beam FWHM |  $28^{\circ}$ at 122 MHz
Field of View |  $625\,\text{deg}^{2}$ at 122 MHz
Polarization | Linear X-Y
Time, frequency resolution: |
Raw Data | 2 s, 195.3 kHz
After flagging | 12 s, 195.3 kHz

##  2 Observations and preprocessing

We use the AARTFAAC-HBA wide-field imager to observe an extended region of
$\sim$20-degree radius around the NCP in the frequency range of $114-126$ MHz,
which is the primary observation window of the LOFAR-EoR KSP (Patil et al.,
2017; Mertens et al., 2020). The observational setup and the preprocessing
steps are discussed briefly in the following subsections.

###  2.1 The AARTFAAC Wide-Field Imager

AARTFAAC is a LOFAR-based all-sky radio transient monitor (Kuiack et al.,
2019). It piggybacks on ongoing LOFAR observations and taps the digital signal
streams of individual antenna elements from six (A6 mode that uses “superterp”
stations) or twelve core stations (A12 mode that uses 12 central stations).
The A6 mode consists of 288 dual-polarization receivers (e.g. Low Band Antenna
(LBA) dipoles or HBA tiles) within a 300 m diameter circle, and the A12 mode
consists of 576 such receivers spread across 1.2 km diameter. The digitized
signals from these receivers are tapped and transported to the AARTFAAC
correlator prior to beam-forming. Due to current network capacity limitations,
only 16 sub-bands can be correlated in the 16-bit mode. Each sub-band is 195.3
kHz wide and consists of up to 64 channels providing a maximum frequency
resolution of 3 kHz, with an instantaneous system bandwidth of 3.1 MHz. The
correlator subsystem, located at the Centre for Information Technology (CIT)
at the University of Groningen (the Netherlands), is a GPU based correlator
that maximally handles 576 input signal streams per polarization (1152 total
streams), and produces XX, XY, YX, YY correlations for all receiver pairs for
every frequency channel with 1-second integration. The output correlations can
either be dumped as raw correlations on disks or routed to the AARTFAAC real-
time calibration and imaging pipeline for transient detection (Prasad et al.,
2016; Kuiack et al., 2019). Readers are referred to van Haarlem et al. (2013)
and Prasad et al. (2016) for detailed information about LOFAR and AARTFAAC
system design and observing capabilities.

###  2.2 Observational setup

We use the AARTFAAC-HBA system in A12 mode (A12-HBA hereafter) to observe the
field around the NCP. For our observations, we target the 114-126 MHz
frequency range corresponding to the redshift range of $z=10.2-11.4$. Due to
the currently limited bandwidth of 3.1 MHz, we adopt an observation strategy
where we combine four different observations recorded within a span of a few
weeks to cover the 12 MHz band. For each observation, the sub-bands are
sparsely spread over the 12 MHz band in a frequency-comb configuration (every
4th sub-band correspond to the same observation) but covering interlacing and
non-overlapping sets of frequency channels for each observation. This
observing strategy allows each observation to independently span the 12 MHz
bandwidth. Although higher time and frequency resolution is preferred for RFI-
flagging, we restricted the data resolution to lower values of 2 s and 195 kHz
(to observe the NCP for longer duration) due to the limited capacity of the
AARTFAAC storage system. The raw correlation data is later transported to the
LOFAR-EoR KSP processing system (Dawn HPC cluster, Pandey et al. 2020) where
all of the data processing and analysis is performed. The raw data is
converted to a standard Measurement Set (MS) format using
aartfaac2ms111https://github.com/aroffringa/aartfaac2ms(Offringa et al., 2015)
and the raw visibilities are phased to the NCP before preprocessing. Note that
the data we used throughout the analysis is only a subset of one such
observation, i.e. a single sub-band at 122.06 MHz. The observational details
of this dataset are summarised in table 1. The analysis of the entire dataset,
including a power-spectrum analysis, will be presented in a future
publication.

###  2.3 Data Preprocessing

The preprocessing steps include RFI-flagging and averaging of the raw
visibilities. We use aoflagger (Offringa et al., 2010, 2012) to flag RFI-
corrupted data. We also flagged all the visibilities corresponding to non-
operational tiles in the datasets. The remaining data were averaged to a lower
time resolution of 12 seconds. The data volume of an 11-hour observation of
3.1 MHz bandwidth after averaging is $\sim 1.7$ TB. We noted that 5 out of 6
stations outside the “superterp” had significantly lower visibility amplitudes
in all observations recorded during this first pilot program. Furthermore, the
visibilities corresponding to these stations showed extremely erratic
behaviour even after calibration. Therefore, we also flagged the tiles
corresponding to these five stations before proceeding with the calibration
step.

##  3 Calibration and imaging

Visibilities measured by a radio interferometer are corrupted by errors due to
instrumental imperfections and environmental effects. These corruptions are
broadly classified into two categories, viz. Direction Independent (DI)
effects such as complex receiver gains, global band-pass and a global
ionospheric phase, and Direction Dependent (DD) effects that change with the
apparent direction of the incoming electromagnetic signals due to the antenna
voltage patterns (and consequently the tile beam), ionospheric phase
fluctuations and Faraday rotation. Calibration refers to estimating the gains
corresponding to these effects to obtain a reliable estimate of the true
visibilities after the gain correction. We perform gain calibration, in a
self-calibration manner, using the following steps:

  1. 1.

The first step is to remove the bright sources Cas A and Cyg A (with the
apparent flux of several hundred Jansky) that superpose significant PSF side-
lobes over the field. We use the DD-calibration option within sagecal-co
(Yatawatta, 2015, 2016; Yatawatta et al., 2017) to subtract these sources from
the raw visibilities. sagecal-co regularizes the gain solutions with a
spectrally smooth function in the gain calibration to improve the calibration
quality. We use the Cas A and Cyg A shapelet model adopted from the LOFAR-EoR
calibration sky-model to subtract these sources using their inferred
directional gains. We use a solution interval of 2 minutes and 5 ADMM
iterations with a regularization factor $\rho=5$ (Yatawatta, 2016). We exclude
$|\boldsymbol{u}|<10$ (in units of wavelength) baselines from the calibration
to avoid most of the large-scale diffuse emission biasing the calibration (but
note that diffuse emission extends over a much wider range of baselines as
shown in later sections).

Figure 1: The left panel shows the Stokes $I$ dirty image of a $\sim$20 degree
radius field around the NCP using a single sub-band (122 MHz) of A12-HBA data
and the baselines in the range $\boldsymbol{u}\leq 500$). The right panel
shows the Stokes $I$ image for the same field with the same imaging parameters
and baseline range, but using LOFAR-HBA station beam-formed data that lacks
the shortest baselines. The dotted circles and spokes represent different
declinations and right ascension, respectively.

  2. 2.

We perform DI-calibration on the resulting visibilities after step 1. The
initial calibration model consists of six sources, viz. 3C61.1 (intrinsic flux
$F_{\text{int}}\sim$ 38 Jy at 151.5 MHz; Baldwin et al. 1985), 3C220.3
($F_{\text{int}}\sim$38 Jy at 150 MHz; Cohen et al. 2007), LQAC 244+085 001
($F_{\text{int}}\sim$6.2 Jy at 151.5 MHz; Baldwin et al. 1985), NVSS
J011045+873822 ($F_{\text{int}}\sim$5.1 Jy at 151.5 MHz; Baldwin et al. 1985),
NVSS J190350+853648 ($F_{\text{int}}\sim$4.9 Jy at 151.5 MHz; Baldwin et al.
1985) and NVSS J062205+871948 ($F_{\text{int}}\sim$4.9 Jy at 151.5 MHz;
Baldwin et al. 1985). All sources, except for 3C61.1, are represented by delta
functions, and a power-law representing their corresponding source spectra
with an assumed, but typical, spectral index of $-0.8$. The radio galaxy
3C61.1 has a more complicated morphology, and its model is adapted from the
intrinsic sky-model used in the LOFAR-EoR calibration pipeline (Mertens et
al., 2020). The 3C61.1 model uses delta functions and shapelets to represent
the source, and a third-order log-polynomial represents the spectrum. We use a
calibration solution interval of 2 minutes to maintain a significant Signal-
to-Noise ratio per solution and carry out 5 ADMM iterations with a
regularization parameter of $\rho=5$ using sagecal-co (Yatawatta, 2016). We
also remove the $|\boldsymbol{u}|<10$ baselines during this step. We do not
use a beam model during calibration, and the flux scale of the visibilities
post-calibration is on an apparent scale (i.e. uncorrected for the average
beam). Currently, the primary beam model for AARTFAAC-HBA tiles is not
implemented in sagecal-co, and efforts are ongoing to utilize the present
LOFAR-HBA tile beam model for AARTFAAC-HBA calibration. Therefore, to avoid
notable differential primary beam effects, all six sources chosen for this
intermediate gain calibration step reside within a $\sim 7^{\circ}$ radius
around the NCP. Since the primary beam does not change substantially within
this radius, the corrected visibilities are close to the correct flux-scale.

  3. 3.

The calibrated visibilities are imaged and deconvolved using the multiscale
CLEAN feature of wsclean (Offringa et al., 2014; Offringa & Smirnov, 2017),
using a cleaning threshold of $0.7\sigma$, with a ‘Briggs –0.1’ weighting-
scheme and $|\boldsymbol{u}|>50$ baselines in order to avoid bright degree-
scale diffuse emission severely affecting the deconvolution process.

  4. 4.

The steps 2 and 3 were repeated three times in a ”self-cal” loop using the
improved calibration model (updated with the clean components) obtained after
each iteration. For each iteration, we use an auto-masking threshold (Offringa
& Smirnov, 2017) of $8\sigma$, $5\sigma$, and $3\sigma$ respectively, where
$\sigma$ is the image rms. The final CLEAN model contains $4924$ components
(delta functions and Gaussians) to represent the compact sources in the field.

  5. 5.

We finally image the DI-calibrated visibilities obtained after step 4 with
wsclean using a ‘Briggs +0.5’ weighting scheme in all subsequent analyses.

Figure 1 shows single sub-band Stokes $I$ dirty images of the NCP field
produced using A12-HBA and LOFAR-HBA data at 122.06 MHz (195.3 kHz bandwidth).
The A12-HBA system has a $\sim 25$ times larger field of view compared to the
LOFAR-HBA system. Galactic large-scale diffuse emission is clearly visible
around NCP in the A12-HBA image (see also Bernardi et al. 2010). Although the
diffuse emission is also partly present in LOFAR-HBA, it is not visible
because LOFAR-HBA has a much lower short-baseline density and a total absence
of baselines shorter than 40 metres ($|\boldsymbol{u}|<15$ or $\ell<100$ at
122 MHz) compared to A12-HBA.

##  4 Modelling the Diffuse Galactic Emission

In this section, we present two approaches to model the diffuse structure
around the NCP as seen in figure 1. In the first approach, we use multiscale
CLEAN deconvolution with wsclean (Offringa & Smirnov, 2017) to model the
diffuse structure. In the second method, we use an orthonormal set of Hermite
polynomial basis functions called “shapelets” (Refregier, 2003; Yatawatta,
2011) to model the diffuse structure.

Figure 2: The left panel shows the deconvolved Stokes $I$ image of the NCP at
122 MHz, produced using a ’Briggs –0.1’ weighting scheme and with
$\boldsymbol{u}>50$ baselines. The right panel shows the Stokes $I$ residual
image after subtraction of the model with 4924 CLEAN components from the data.
The residuals are consistent with the confusion noise. Figure 3: The Stokes
$I$ dirty image of the residual visibilities after the subtraction of 4924
CLEAN model components. The image was produced using a ’Briggs +0.5’ weighting
scheme (which gives more weight to short baselines to better visualize the
diffuse emission), $\boldsymbol{u}\leq 120$ baseline range, and a Gaussian
taper of $30\arcmin$.

###  4.1 Removal of compact sources

Before proceeding with the diffuse foreground modelling, the bright compact
sources need to be removed from the map. We use the CLEAN component model
obtained during the gain calibration (see step 4 in Section 3) to subtract
4924 components from the data. Figure 2 shows Stokes $I$ image of the NCP
field before and after the subtraction of the CLEAN model. The images are
produced using ’Briggs –0.1’ weighting scheme with $\boldsymbol{u}>50$
baselines. We observe that subtracting the CLEAN model removes most of the
compact sources, leaving only faint residuals of the order of a few hundred
mJy rms. The images do not visually reveal signatures of diffuse emission
because of the baseline cut used during the imaging process. We also see that
the subtraction of some of the brightest sources leaves negative or ring-like
artifacts near their source location. Factors such as imperfect calibration
and imperfect source modelling are the likely cause. These artifacts can be
further mitigated by using DD-calibration to subtract the bright sources. We
plan to employ the DD-calibration for the subtraction of bright sources using
several directions in future analyses. For the purpose of this work, i.e.
analysing the properties of the diffuse foregrounds, these minor artifacts are
not important.

We derive an estimate of the confusion noise for AARTFAAC from the classical
confusion limit ($\sigma_{c}$) using the relation from van Haarlem et al.
(2013):

$\sigma_{c}=30\left(\dfrac{\theta}{1^{\arcsec}}\right)^{1.54}\left(\dfrac{\nu}{74\,\textrm{MHz}}\right)^{-0.7}\
[\mu\textrm{Jy}\,\textrm{beam}^{-1}],$ (1)

where $\theta$ is the angular resolution (FWHM), and $\nu$ is the observation
frequency. A12-HBA has an angular resolution of $\sim 7\arcmin$ at 122 MHz,
yielding a confusion limit of $\sigma_{c}\sim 230$ mJy beam-1. The standard
deviation of the residuals, after cleaning, for the inner $10^{\circ}$ region
of the field is $\sim 220$ mJy. These two values are consistent with each
other, assuming that the primary beam correction is small for this inner
$10^{\circ}$ region. In the future, we plan a deeper multi-frequency CLEAN by
combining subbands spanning a wide frequency range along with an improved
direction-dependent calibration step that includes an HBA-tile beam model.
Figure 3 shows the Stokes $I$ image of the residual visibilities after
subtraction of the 4924 component CLEAN model, produced using
$\boldsymbol{u}\leq 120$ baselines. We now model this residual diffuse
emission using the two methods mentioned above.

Figure 4: Top row: Stokes $I$ dirty image before subtraction of the diffuse
structure (left panel), the diffuse model obtained with multiscale CLEAN
(middle panel), and the residual Stokes $I$ image after subtraction of the
CLEAN model (right panel). Bottom row: same as the top row but with shapelet
decomposition technique to model the diffuse structure. These images were
produced using ’Briggs +0.5’ and $\boldsymbol{u}\leq 120$ baselines with a
Gaussian taper of $30\arcmin$.

###  4.2 Modeling with multiscale CLEAN

We use the multiscale deconvolution algorithm implemented in wsclean (Offringa
& Smirnov, 2017) to model the diffuse structure observed in figure 3. We
perform the deconvolution with the following settings: ‘Briggs +0.5’ weighting
scheme, auto-mask of $3\sigma$, major iteration cleaning gain (mgain) of
$0.8$, and a clean stopping threshold of $0.6\sigma$. The ‘Briggs’ weighting
scheme with a threshold greater than zero gives more weight to short baselines
and reveals the diffuse emission better. Figure 4 (top row) shows the Stokes
$I$ image of the field before and after subtracting the diffuse foreground
model obtained from multiscale CLEAN. We observe that multiscale CLEAN is able
to model the extended emission at different scales using Gaussians and point
sources, adequately capturing the diffuse flux. The algorithm also models the
compact residuals (both positive and negative) initially missed by the compact
model subtraction step. After the subtraction of the diffuse model, the
residual appears to be devoid of diffuse structure, with a standard deviation
of $\sim 500$ mJy inside the inner region of the image. However, we notice
that there is a lack of intermediate scales in the diffuse model. Note that
the CLEAN algorithm is sensitive to image weighting and sampling function
applied during the deconvolution process. A possible improvement in the model
may be achieved by tuning the image weighting, scales used for deconvolution,
and adding more data in $uv$-plane to cover the (unmodeled) spatial scales
with lower sensitivity. For the specific case presented in the paper, we think
that multiscale CLEAN, in its current implementation, is a good choice for
modelling large-scale diffuse emission.

###  4.3 Modelling with Shapelets

Yatawatta (2011) showed that the use of orthonormal basis functions e.g.
Cartesian shapelets (Refregier, 2003) or Prolate Spheroidal Wave Functions
(PWSF) (Slepian & Pollak, 1961; Landau & Pollak, 1961) in radio
interferometric image deconvolution, provides improved image fidelity and
larger dynamic range. In the image domain, $(l,m)$, shapelet basis functions
can be written as (Yatawatta, 2011):

$\displaystyle\phi_{n_{1},n_{2}}(l,m,\beta)=\dfrac{1}{2^{n_{1}+n_{2}}\,\piup\beta^{2}n_{1}!\,n_{2}!}\,H_{n_{1}}(l/\beta)\,H_{n_{2}}(m/\beta)$
$\displaystyle\times\exp{[-(l^{2}+m^{2})/2\beta^{2}]},$ (2)

where the functions $H_{n_{1}}$ and $H_{n_{2}}$ are the Hermite polynomials of
order $n_{1}$ and $n_{2}$ being integer values, the value of $\beta$ is the
model scale factor. We demonstrate the use of Cartesian shapelets to model the
diffuse structure observed in A12-HBA data using the shapelet_gui222developed
by Sarod Yatawatta. See LOFAR imaging cookbook for more details:
https://support.astron.nl/LOFARImagingCookbook/ tool. We use $25\times 25$
basis functions for the shapelet decomposition using L1-norm regularization.
We optimize the model scale during the process. Figure 4 (bottom row) presents
the Stokes $I$ images before and after shapelet model subtraction, as well as
the shapelet model of the diffuse structure obtained from the decomposition.
We note that the shapelet model captures the large-scale diffuse structure
within the primary beam. The two negative-flux holes in the model were masked
during the decomposition. The residuals within the primary beam have a
standard deviation of $\sim 540$ mJy, which is approximately the same as
residuals after CLEAN diffuse model subtraction. These residuals are most
likely the unmodeled compact sources and sources below the confusion noise.
However, the diffuse structure outside the primary beam is rather poorly
modelled. We suspect that this is mainly due to the finite support of the
basis functions and a sub-optimized scale factor ($\beta$) that corresponds to
the width of the basis functions, leading to a inadequate performance of the
shapelet model of capturing the structure outside the primary beam. We expect
that the performance of the shapelet decomposition modelling may be enhanced
by optimizing the number of basis functions used for the decomposition, and
their scale factor ($\beta$). Additionally, using multiple shapelet models to
capture the diffuse structure in different regions in the image space and
using primary-beam corrected images is expected to improve the shapelet model.

###  4.4 Comparing the two diffuse emission modelling approaches

We find that, in this particular setup, the multiscale CLEAN deconvolution
method captures compact sources and diffuse structures at a range of angular
scales. On the other hand, the shapelet decomposition captures the structure
on large scales significantly better, but it is less effective (without hugely
increasing the number of basis functions and shapelets) in modelling the
fluctuations on scales smaller than a degree. From the visual comparison of
the two methods in figure 4, we infer that in the current implementation, both
multiscale CLEAN and shapelet decomposition are able to model the large-scale
(order of several degrees) diffuse Galactic emission. However, the shapelet
model performs poorly outside the primary beam and is inefficient at capturing
the diffuse structure on small scales ($\lesssim 1$ degree). A possible
improvement in the CLEAN model may be achieved by tweaking the image weighting
scheme, and the scales used in deconvolution. Similarly, optimizing the number
of basis functions, their scale factor, and using multiple shapelets are
expected to enhance the shapelet model. However, the CLEAN deconvolution
provides the added benefit of fitting a spectrum to capture the spectral
behaviour of the diffuse structure, which is currently not the case with
shapelet decomposition. Inclusion of spectral indices (per basis function) in
shapelet model is currently under development. For future analyses, we plan to
improve the analysis and modelling process by adding more data over a wide
range of frequencies, and tweaking various parameters described above. We also
plan to explore a hybrid approach to build the sky model by combining the two
methods, where the CLEAN model is used to capture small-scale structures and
compact sources, and shapelets are used to represent the large-scale diffuse
structures. The more complete model of diffuse emission will consist of
information on a wide range of spatial scales along with spectral information.
This model can also be utilized to improve the calibration and foreground
removal (using DD-calibration) of the LOFAR-EoR and AARTFAAC-HBA observations.

Figure 5: Comparison of the angular power spectrum ($\Delta_{\ell}^{2}$) of
the A12-HBA Stokes $I$ image before and after subtraction of the 4924
component CLEAN model. The solid orange, blue, and grey curves correspond to
the emission in the field (data before point source CLEAN model subtraction),
the point source CLEAN model, and the residual emission, respectively. The
dashed blue and grey curves represent the best fit profiles of the point
source CLEAN model ($\Delta_{\textrm{CLEAN}}^{2}(\ell)$) and the residual
emission ($\Delta_{\textrm{residual}}^{2}(\ell)$), respectively. The error
bars on data correspond to the $2\sigma$ uncertainty due to sample variance.

##  5 The Angular power spectrum

The angular power spectrum is a commonly used metric to study the spatial
properties of the foreground emission. We use the angular power spectrum,
$C_{\ell}$, to describe the observed diffuse emission in the NCP field. The
angular power spectrum $C_{\ell}$ is defined as (see e.g., Seljak, 1997;
Bernardi et al., 2009):

$C_{\ell}=\dfrac{1}{N_{\ell}}\sum_{\boldsymbol{l}}\tilde{I}(\boldsymbol{l})\tilde{I}^{*}(\boldsymbol{l}),$
(3)

where $\tilde{I}(\boldsymbol{l})$ is the spatial Fourier transform of an image
$I$ (in units of Kelvins), $N_{\ell}$ is the number of samples in an azimuthal
bin, $\ell$ is the multipole moment and
$\ell=2\piup|\boldsymbol{u}|=180^{\circ}/\Theta$, and $\Theta$ is the angular
scale in degrees (Half-Width Half Maximum). We use the Stokes $I$ images to
estimate $C_{\ell}$. The images are Fourier transformed along the spatial axes
using Fast Fourier Transform (FFT). We multiply the images with a spatial
Tukey (tapered cosine) window before the FFT to restrict the emission outside
the $15^{\circ}$ radius, left after shapelet modelling, affecting the results.
We normalize the Fourier transform of the images with the PSF to remove the
effect of the weighting scheme applied during imaging. Note that we do not
apply a primary-beam correction during the estimation of $C_{\ell}$, nor
correct for curvature of the sky, which are both limited in their effect.

Figure 5 presents the angular power spectrum
$\Delta^{2}_{\ell}\equiv\ell(\ell+1)C_{\ell}/2\piup$ of the Stokes $I$ image
before and after subtraction of the 4924 components of the point source CLEAN
model. We do not subtract the noise bias, however, it is expected to be a
small fraction of the observed power. We observe that the power on smaller
scales (large $\ell$ values) is dominated by compact sources and is mostly
removed after subtraction of the CLEAN model. However, power on $\ell<60$
remains intact after model subtraction. Moreover, the angular power spectrum
of the residual emission is flat over a wide range of multipole moments,
ranging from $\ell=20$ to $200$, which corresponds to angular scales of
$\Theta=0.9-10^{\circ}$. To further quantify the point source CLEAN model and
residual emission power spectra, we fit the corresponding angular power
spectra respectively with the following functions:

${\rm Point\leavevmode\nobreak\ source\leavevmode\nobreak\
CLEAN\leavevmode\nobreak\ model:\leavevmode\nobreak\ }\Delta^{2}_{\ell,{\rm
CLEAN}}=A\left(\dfrac{\ell}{\ell_{0}}\right)^{\alpha},$ (4) ${\rm
Residual\leavevmode\nobreak\ emission:\leavevmode\nobreak\
}\Delta^{2}_{\ell,{\rm
Residual}}=A^{\prime}+B^{\prime}\left(\dfrac{\ell}{\ell_{0}}\right)^{\alpha^{\prime}},$
(5)

where $A$, $\alpha$, $A^{\prime}$, $B^{\prime}$, and $\alpha^{\prime}$ are the
free parameters for the CLEAN model fit, and the residual emission fit,
respectively. We choose $\ell_{0}=180$ which correspond to an angular scale
$\Theta=1^{\circ}$. The best-fit parameter values for the two models are
listed in Table 2. The data and the point source CLEAN model follow power-law
like behaviour on $\ell\gtrsim 60$. Fitting the angular power spectrum of the
point source CLEAN Model with equation 4 yields a best fit power-law index
$\alpha=1.99\pm 0.09$, which is consistent with the angular power spectrum of
unresolved and unclustered sources ($\Delta_{\ell}^{2}\propto\ell^{2}$). This
suggests that the compact sources begin dominating the power on baselines
$\boldsymbol{u}\gtrsim 10$ or angular scales of $\Theta\lesssim 3^{\circ}$
($\ell\gtrsim 60$). The angular power spectrum ($\Delta_{\ell}^{2}$) of the
residual emission is flat within $20\lesssim\ell\lesssim 200$, suggesting that
$C_{\ell}\propto\ell^{-2}$. This is consistent with the angular power spectrum
of the Galactic diffuse emission, $C_{\ell}\propto\ell^{-2.2}$ observed by
Bernardi et al. (2009). However, the power on $\ell>200$ is possibly dominated
by thermal noise and low level emission below confusion noise, and starts to
follow a power-law like behaviour. As revealed from the fit, the residual
diffuse emission has a brightness temperature variance of
$\Delta^{2}_{\ell=180}=(145.64\pm 13.61)\,{\rm K}^{2}$ at 122 MHz on angular
scales of 1 degree, and is consistent with $C_{\ell}\propto\ell^{-2.0}$ in
$20\lesssim\ell\lesssim 200$ range. Using observations of the NCP with
Westerbork Synthesis Radio Telescope (WSRT), Bernardi et al. (2010) reported
$\Delta^{2}_{\ell=120}\sim 14\,{\rm K}^{2}$ at 150 MHz, which corresponds to
$\Delta^{2}_{\ell=120}\sim 40\,{\rm K}^{2}$ at 122 MHz when extrapolated with
the spectral index of $-2.55$. Additionally, Offringa et al. (2022) observed
$\Delta^{2}_{\ell=120}\sim 3\,{\rm K}^{2}$ for the field centred at “coldest
patch” of the northern Galactic hemisphere (Kogut et al., 2011), and
$\Delta^{2}_{\ell=120}\sim 10\,{\rm K}^{2}$ for another field $\sim
10^{\circ}$ away from the “coldest patch” at 140 MHz with LOFAR. Whereas, we
observe $\Delta^{2}_{\ell=120}\sim 130\,{\rm K}^{2}$ in our analysis, that is
around quadruple the power observed by Bernardi et al. (2010). We hypothesize
that the high power observed in our analysis is possibly due to the limited
bandwidth of AARTFAAC data ($195$ kHz) used in this study compared to the
analyses in Bernardi et al. (2010) ($40$ MHz) and Offringa et al. (2022) ($48$
MHz). Moreover, we do not apply any noise correction (instrumental and
confusion) contrary to Bernardi et al. (2010). We defer the latter for the
future analyses.

Table 2: Best-fit Parameters Parameter | Best-fit value
---|---
CLEAN point sources model
$A$ | $1438.75\pm 111.08$
$\alpha$ | $1.99\pm 0.09$
Residual emission model
$A^{\prime}$ | $128.16\pm 13.45$
$B^{\prime}$ | $17.48\pm 2.08$
$\alpha^{\prime}$ | $3.00\pm 0.47$
Figure 6: The angular power spectra ($\Delta_{\ell}^{2}$) of the residual
emission (grey curve) same as the grey curve in figure 5, CLEAN residual after
subtraction of the CLEAN diffuse model (blue curve), and shapelet residual
after subtraction of the shapelet diffuse model (orange curve) shown in figure
4.

###  5.1 Power spectra from diffuse foreground modelling methods

The angular power spectrum is also a useful metric to compare the two
foreground modelling methods, i.e., multiscale CLEAN and shapelet
decomposition, that we demonstrated in Section 4. We use the images shown in
figure 4 for the two methods to determine the corresponding angular power
spectra $\Delta_{\ell}^{2}$. We use the same methodology as mentioned in the
previous section to determine $C_{\ell}$. Figure 6 shows the angular power
spectra $\Delta_{\ell}^{2}$ of the residuals after subtraction of CLEAN and
shapelet diffuse models obtained from the two methods. We observe that the
multiscale CLEAN method leads to lower residuals than the shapelet model on
$\ell<70$ modes ($\Theta\lesssim 2.3^{\circ}$). However, the shapelet model
produces slightly lower residuals on $70\lesssim\ell<200$ compared to the
multiscale CLEAN. This is also evident from the visual inspection of the two
models in figure 4. The residuals on $\ell\gtrsim 200$ behave similarly to the
data, probably due to the noise in the data. Note that the two methods (and
the residuals obtained from these) compared in this analysis, use a very
specific setup, and a detailed investigation to study and quantify finer
differences between the two methods will be performed in future analyses.

##  6 Summary and future work

In this work, we presented the first-ever wide-field images obtained with the
LOFAR AARTFAAC-HBA system in A12 mode. In particular, we find strong degree-
scale diffuse Galactic Stokes $I$ radio emission at 122 MHz within a
$\sim$20-degree radius field around the North Celestial Pole, which is one of
the primary windows of the LOFAR EoR project (Yatawatta et al., 2013; Mertens
et al., 2020). We have compared two different methods for modelling this
diffuse emission namely, multiscale CLEAN deconvolution and shapelet
decomposition. We use angular power spectrum to quantify the behaviour of
different foreground components, viz. point sources and diffuse emission. The
main results of this work are summarised as follows:

  * •

Stokes $I$ radio emission, as seen by the LOFAR AARTFAAC-HBA system around the
NCP at 122 MHz on baselines that are shorter than 120 wavelengths, is
dominated by large-scale diffuse emission. The angular power spectrum of the
emission in the field is dominated by the point sources on scales smaller than
3 degrees ($\ell\gtrsim 60$). After subtraction of the point source CLEAN
model, the residuals are dominated by the diffuse emission on scales larger
than a degree ($\ell\lesssim 200$). This diffuse emission,
with$\leavevmode\nobreak\ 10^{4}$ times more power than the 21-cm signal, can
have a considerable impact on the calibration of any radio-interferometric
instrument (e.g., LOFAR, MWA, HERA, NenuFAR, LWA, and SKA) if this emission is
not part of the sky model during instrumental gain calibration (or filtered
out before calibration) or if their baselines are not excluded during
calibration.

  * •

We show that, with our particular setup, multiscale CLEAN can model the small
and intermediate scales well, but shows a slightly lower performance on
intermediate scales than the shapelet decomposition, which captures the large-
scale diffuse emission well but it is incapable of modelling the emission on
scales smaller than several degrees.

  * •

The diffuse emission has a brightness temperature variance of
$\Delta^{2}_{\ell=180}=(145.64\pm 13.61)\,{\rm K}^{2}$ at 122 MHz on angular
scale of 1 degree, and is consistent with $C_{\ell}\propto\ell^{-2.0}$ in
$20\lesssim\ell\lesssim 200$ range.

The analysis in this pilot work is admittedly based only on data from a single
sub-band at 122 MHz. Despite this restriction, for the first time using
AARTFAAC-HBA as a wide-field imaging instrument, we have convincingly shown
that large-scale diffuse emission dominates the emission in Stokes $I$ on
degree scales and larger. This was previously much harder to assess due to the
limited short-baseline coverage of LOFAR in beam-formed mode. To further
improve the imaging and modelling of this emission, we plan to include the
following future analyses:

  * •

Including an HBA-tile beam model in the gain calibration, which allows us to
(i) set an absolute flux scale, (ii) improve the calibration quality, (iii)
perform deeper deconvolution during imaging and (iv) improve the modelling of
compact sources and diffuse emission.

  * •

Expand this analysis to a broader frequency range to obtain a complete spatial
and spectral model of the diffuse emission.

  * •

Explore a hybrid approach to model the emission on a wide range of spatial
scales by combining multiscale CLEAN to model small-scale compact emission and
shapelet decomposition to model large-scale diffuse emission.

  * •

Extend the analysis to model diffuse polarized emission in Stokes $Q$ and $U$.
Wide-field polarization studies will enable us to understand the behaviour and
morphology of polarized foregrounds on very large scales and how it
instrumentally leaks to Stokes $I$ (Asad et al., 2015, 2016, 2018; Nunhokee et
al., 2017; Byrne et al., 2021).

  * •

Include the diffuse emission model in the sky-model as part of the LOFAR-EoR
and NenuFAR-CD calibration pipeline.

## Software

This research made use of publicly available software developed for LOFAR and
AARTFAAC telescopes. Listed below are software packages used in the analysis:
aartfaac2ms (https://github.com/aroffringa/aartfaac2ms), aoflagger
(https://gitlab.com/aroffringa/aoflagger), sagecal-co
(http://sagecal.sourceforge.net/), and wsclean
(https://gitlab.com/aroffringa/wsclean). The analysis also relies on the
Python programming language (https://www.python.org) and several publicly
available python software modules: astropy (https://www.astropy.org/; The
Astropy Collaboration et al. 2013), matplotlib (https://matplotlib.org/),
scipy (https://www.scipy.org/) and numpy (https://numpy.org/).

###### Acknowledgements.

BKG and LVEK acknowledge the financial support from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 884760, ”CoDEX”). FGM and MK acknowledge
support from a SKA-NL roadmap grant from the Dutch Ministry of OCW. This work
was supported in part by ERC grant 247295 “AARTFAAC” to RAMJW. LOFAR, the Low
Frequency Array designed and constructed by ASTRON, has facilities in several
countries, that are owned by various parties (each with their own funding
sources), and that are collectively operated by the International LOFAR
Telescope (ILT) foundation under a joint scientific policy.

## References

  * Asad et al. (2018) Asad, K. M. B., Koopmans, L. V. E., Jelić, V., et al. 2018, MNRAS, 476, 3051
  * Asad et al. (2016) Asad, K. M. B., Koopmans, L. V. E., Jelić, V., et al. 2016, MNRAS, 462, 4482
  * Asad et al. (2015) Asad, K. M. B., Koopmans, L. V. E., Jelić, V., et al. 2015, MNRAS, 451, 3709
  * Baldwin et al. (1985) Baldwin, J. E., Boysen, R. C., Hales, S. E. G., et al. 1985, MNRAS, 217, 717
  * Barry et al. (2016) Barry, N., Hazelton, B., Sullivan, I., Morales, M. F., & Pober, J. C. 2016, MNRAS, 461, 3135
  * Bernardi et al. (2009) Bernardi, G., de Bruyn, A. G., Brentjens, M. A., et al. 2009, A&A, 500, 965
  * Bernardi et al. (2010) Bernardi, G., de Bruyn, A. G., Harker, G., et al. 2010, A&A, 522, A67
  * Bernardi et al. (2016) Bernardi, G., Zwart, J. T. L., Price, D., et al. 2016, MNRAS, 461, 2847
  * Bowman et al. (2013) Bowman, J. D., Cairns, I., Kaplan, D. L., et al. 2013, PASA, 30, e031
  * Bowman et al. (2018) Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67
  * Byrne et al. (2021) Byrne, R., Morales, M. F., Hazelton, B., et al. 2021, arXiv e-prints, arXiv:2107.11487
  * Chapman et al. (2013) Chapman, E., Abdalla, F. B., Bobin, J., et al. 2013, MNRAS, 429, 165
  * Chapman et al. (2012) Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, MNRAS, 423, 2518
  * Clark (1980) Clark, B. G. 1980, A&A, 89, 377
  * Cohen et al. (2007) Cohen, A. S., Lane, W. M., Cotton, W. D., et al. 2007, AJ, 134, 1245
  * Cornwell (2008) Cornwell, T. J. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793
  * Datta et al. (2010) Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526
  * de Lera Acedo (2019) de Lera Acedo, E. 2019, in 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA), 0626–0629
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001
  * Di Matteo et al. (2002) Di Matteo, T., Perna, R., Abel, T., & Rees, M. J. 2002, ApJ, 564, 576
  * Dillon et al. (2014) Dillon, J. S., Liu, A., Williams, C. L., et al. 2014, Phys. Rev. D, 89, 023002
  * Dowell et al. (2017) Dowell, J., Taylor, G. B., Schinzel, F. K., Kassim, N. E., & Stovall, K. 2017, MNRAS, 469, 4537
  * Eastwood et al. (2018) Eastwood, M. W., Anderson, M. M., Monroe, R. M., et al. 2018, AJ, 156, 32
  * Eastwood et al. (2019) Eastwood, M. W., Anderson, M. M., Monroe, R. M., et al. 2019, AJ, 158, 84
  * Ewall-Wice et al. (2017) Ewall-Wice, A., Dillon, J. S., Liu, A., & Hewitt, J. 2017, MNRAS, 470, 1849
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep, 433, 181
  * Gehlot et al. (2019) Gehlot, B. K., Mertens, F. G., Koopmans, L. V. E., et al. 2019, MNRAS, 488, 4271
  * Gehlot et al. (2020) Gehlot, B. K., Mertens, F. G., Koopmans, L. V. E., et al. 2020, MNRAS, 499, 4158
  * Ghosh et al. (2012) Ghosh, A., Prasad, J., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2012, MNRAS, 426, 3295
  * Harker et al. (2009) Harker, G., Zaroubi, S., Bernardi, G., et al. 2009, MNRAS, 397, 1138
  * Hazelton et al. (2013) Hazelton, B. J., Morales, M. F., & Sullivan, I. S. 2013, ApJ, 770, 156
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Kogut et al. (2011) Kogut, A., Fixsen, D. J., Levin, S. M., et al. 2011, ApJ, 734, 4
  * Koopmans et al. (2015) Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 1
  * Kuiack et al. (2019) Kuiack, M., Huizinga, F., Molenaar, G., et al. 2019, MNRAS, 482, 2502
  * Landau & Pollak (1961) Landau, H. J. & Pollak, H. O. 1961, Bell System Technical Journal, 40, 65
  * Line et al. (2020) Line, J. L. B., Mitchell, D. A., Pindor, B., et al. 2020, PASA, 37, e027
  * Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429
  * Mertens et al. (2020) Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662
  * Mertens et al. (2021) Mertens, F. G., Semelin, B., & Koopmans, L. V. E. 2021, arXiv e-prints, arXiv:2109.10055
  * Mohan & Rafferty (2015) Mohan, N. & Rafferty, D. 2015, PyBDSF: Python Blob Detection and Source Finder, Astrophysics Source Code Library
  * Morales et al. (2012) Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, ApJ, 752, 137
  * Mouri Sardarabadi & Koopmans (2019) Mouri Sardarabadi, A. & Koopmans, L. V. E. 2019, MNRAS, 483, 5480
  * Nunhokee et al. (2017) Nunhokee, C. D., Bernardi, G., Kohn, S. A., et al. 2017, ApJ, 848, 47
  * Offringa et al. (2010) Offringa, A. R., de Bruyn, A. G., Biehl, M., et al. 2010, MNRAS, 405, 155
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606
  * Offringa et al. (2022) Offringa, A. R., Singal, J., Heston, S., Horiuchi, S., & Lucero, D. M. 2022, MNRAS, 509, 114
  * Offringa & Smirnov (2017) Offringa, A. R. & Smirnov, O. 2017, MNRAS, 471, 301
  * Offringa et al. (2012) Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, A95
  * Offringa et al. (2015) Offringa, A. R., Wayth, R. B., Hurley-Walker, N., et al. 2015, PASA, 32, e008
  * Paciga et al. (2013) Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639
  * Pandey et al. (2020) Pandey, V. N., Koopmans, L. V. E., Tiesinga, E., Albers, W., & Koers, H. U. A. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, 473
  * Parsons et al. (2012) Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, ApJ, 756, 165
  * Patil et al. (2017) Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, ApJ, 838, 65
  * Patil et al. (2016) Patil, A. H., Yatawatta, S., Zaroubi, S., et al. 2016, MNRAS, 463, 4317
  * Prasad et al. (2016) Prasad, P., Huizinga, F., Kooistra, E., et al. 2016, Journal of Astronomical Instrumentation, 5, 1641008
  * Refregier (2003) Refregier, A. 2003, MNRAS, 338, 35
  * Seljak (1997) Seljak, U. 1997, ApJ, 482, 6
  * Shaver et al. (1999) Shaver, P. A., Windhorst, R. A., Madau, P., & de Bruyn, A. G. 1999, A&A, 345, 380
  * Singh et al. (2017) Singh, S., Subrahmanyan, R., Udaya Shankar, N., et al. 2017, ApJ, 845, L12
  * Slepian & Pollak (1961) Slepian, D. & Pollak, H. O. 1961, Bell System Technical Journal, 40, 43
  * The Astropy Collaboration et al. (2013) The Astropy Collaboration, Robitaille, Thomas P., Tollerud, Erik J., et al. 2013, A&A, 558, A33
  * Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007
  * Trott et al. (2012) Trott, C. M., Wayth, R. B., & Tingay, S. J. 2012, ApJ, 757, 101
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * Vedantham et al. (2012) Vedantham, H., Udaya Shankar, N., & Subrahmanyan, R. 2012, ApJ, 745, 176
  * Whiting (2012) Whiting, M. T. 2012, MNRAS, 421, 3242
  * Yatawatta (2011) Yatawatta, S. 2011, in 2011 XXXth URSI General Assembly and Scientific Symposium, 1–4
  * Yatawatta (2015) Yatawatta, S. 2015, MNRAS, 449, 4506
  * Yatawatta (2016) Yatawatta, S. 2016, in 2016 24th European Signal Processing Conference (EUSIPCO), 265–269
  * Yatawatta et al. (2013) Yatawatta, S., de Bruyn, A. G., Brentjens, M. A., et al. 2013, A&A, 550, A136
  * Yatawatta et al. (2017) Yatawatta, S., Diblen, F., & Spreeuw, H. 2017, in 2017 IEEE 7th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 1–5
  * Zaldarriaga et al. (2004) Zaldarriaga, M., Furlanetto, S. R., & Hernquist, L. 2004, ApJ, 608, 622
  * Zarka et al. (2012) Zarka, P., Girard, J. N., Tagger, M., & Denis, L. 2012, in SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. S. Boissier, P. de Laverny, N. Nardetto, R. Samadi, D. Valls-Gabaud, & H. Wozniak, 687–694