# The Correlation Calibration of PAPER-64 data

Tamirat G. Gogo,1,2 Yin-Zhe Ma1,2, Piyanat Kittiwisit1,2, Jonathan L.
Sievers3,1, Aaron R. Parsons4,5, Jonathan C. Pober6, Daniel C. Jacobs7, Carina
Cheng 4, Matthew Kolopanis7, Adrian Liu8,9, Saul A. Kohn10, James E.
Aguirre10, Zaki S. Ali,5, Gianni Bernardi11,12,13, Richard F. Bradley14,15,16,
David R. DeBoer5, Matthew R. Dexter5, Joshua S. Dillon4, Pat Klima15, David H.
E. MacMahon5, David F. Moore10, Chuneeta D. Nunhokee4, William P. Walbrugh13,
and Andre Walker11
1 School of Chemistry and Physics, University of KwaZulu-Natal, Westville
Campus, Private Bag X54001, Durban, 4000, South Africa
2 NAOC-UKZN Computational Astrophysics Centre (NUCAC), University of KwaZulu-
Natal, Durban, 4000, South Africa
3 Department of Physics, McGill University, 3600 University Street Montreal,
QC H3A 2T8, Canada
4 Department of Astronomy, UC Berkeley, Berkeley, CA 94720, USA
5 Radio Astronomy Laboratory, UC Berkeley, Berkeley, CA 94720, USA
6 Department of Physics, Brown University, Providence, RI 02912, USA
7 School of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287, USA
8 Department of Physics and McGill Space Institute, McGill University,
Montreal, QC, Canada
9 CIFAR Azrieli Global Scholar, Gravity & the Extreme Universe Program,
Canadian Institute for Advanced Research, 661 University
Ave., Suite 505, Toronto, Ontario M5G 1M1, Canada
10 Department of Physics and Astronomy, U. Penn., Philadelphia PA, USA
11 INAF-Istituto di Radioastronomia, via Gobetti 101, 40129, Bologna, Italy
12 Department of Physics and Electronics, Rhodes University, PO Box 94,
Grahamstown, 6140, South Africa
13 South African Radio Astronomy Observatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa
14 Department of Electrical and Computer Engineering, U. Virginia,
Charlottesville VA, USA
15 National Radio Astronomy Obs., Charlottesville VA, USA
16 Department of Astronomy, U. Virginia, Charlottesville VA, USA
Corresponding author<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observation of redshifted 21-cm signal from the Epoch of Reionization (EoR) is
challenging due to contamination from the bright foreground sources that
exceed the signal by several orders of magnitude. The removal of this very
high foreground relies on accurate calibration to keep the intrinsic property
of the foreground with frequency. Commonly employed calibration techniques for
these experiments are the sky model-based and the redundant baseline-based
calibration approaches. However, the sky model-based and redundant baseline-
based calibration methods could suffer from sky-modeling error and array
redundancy imperfection issues, respectively. In this work, we introduce the
hybrid correlation calibration (CorrCal) scheme, which aims to bridge the gap
between redundant and sky-based calibration by relaxing redundancy of the
array and including sky information into the calibration formalisms. We
demonstrate the slight improvement of power spectra, about $-6\%$ deviation at
the bin right on the horizon limit of the foreground wedge-like structure,
relative to the power spectra before the implementation of CorrCal to the data
from the Precision Array for Probing the Epoch of Reionization (PAPER)
experiment, which was otherwise calibrated using redundant baseline
calibration. This small improvement of the foreground power spectra around the
wedge limit could be suggestive of reduced spectral structure in the data
after CorrCal calibration, which lays the foundation for future improvement of
the calibration algorithm and implementation method.

###### keywords:

telescopes – instrumentation: interferometers – methods: observational –
cosmology: observations

††pubyear: 2021††pagerange: The Correlation Calibration of PAPER-64 data–B

##  1 Introduction

The cosmic dark ages have ended at around 100 million years after the big
bang, marking the formation of the first structures such as galaxies and
stars. These objects started to emit radiation that influences the nearby
primordial intergalactic medium (IGM), which predominantly filled with the
neutral hydrogen atom. As the number of these very first sources increases,
they emitted sufficient radiation to convert the neutral IGM into the fully
ionized medium of the hydrogen atom eventually. The period that marks the
transition of neutral IGM into the fully ionized state is known as the Epoch
of Reionization (EoR) in the history of the Universe. Thus, investigating the
era of EoR will answer the fundamental questions of the formation of the first
structures and their properties.

The measurement of the radiation from the high-redshifted 21-cm emission line
is employed to probe the EoR era. This emission line is due to the hyperfine
transition of the neutral hydrogen atom in the IGM empowered by the first
sources. The information imprinted into the radiation associated with the
redshifted 21-cm line can be used as a powerful tool in the modern cosmology
to study the ionization state and large-scale temperature fluctuation of IGM,
owing to the formation of the first structures (see, e.g., Oh, 2001; Barkana &
Loeb, 2001; Furlanetto et al., 2006; Morales & Wyithe, 2010; Pritchard & Loeb,
2012).

A direct probe of the EoR is the tomographic mapping of redshifted 21-cm
signals from the neutral hydrogen atom (Madau et al., 1996; Barkana & Loeb,
2001; Furlanetto et al., 2006; Morales & Wyithe, 2010; Pritchard & Loeb,
2012). Significant progress has made in the past decades. Several experiments
dedicated to studying this signal have been built or are being upgraded. The
first-generation 21-cm experiments, including the Murchison Widefield Array
(MWA; Tingay et al., 2013), the Donald C. Backer Precision Array for Probing
the Epoch of Reionization (PAPER; Parsons et al., 2010), the LOw-Frequency
ARray (LOFAR; van Haarlem et al., 2013), and the GiantMeter wave Radio
Telescope EoR experiment (GMRT; Paciga et al., 2011) are already operating and
taking data. Observations from these first-generation experiments have set
upper limits to the statistical power spectrum of 21-cm brightness temperature
from cosmic reionization. In particular, results from MWA, LOFAR, and PAPER
observations (e.g., Paciga et al., 2011; Dillon et al., 2014; Patil et al.,
2017; Cheng et al., 2018; Kolopanis et al., 2019; Mertens et al., 2020) have
placed upper limits on the statistical power spectrum of the 21-cm emissions
over a broad range of redshifts, providing evidence for the heating of the IGM
before reionization. The most recent result from MWA constrains the EoR power
spectrum as $\Delta^{2}=(43\,{\rm mK})^{2}=1.8\times 10^{3}\,{\rm mK}^{2}$ for
$k=0.14\,h\,{\rm Mpc}^{-1}$ at $z=6.5$, which is the lowest measurement
currently (Trott et al., 2020; Rahimi et al., 2021). Lessons learned from
these experiments have led to the construction of the second-generation
instruments, which include the Hydrogen Epoch of Reionization Array (HERA;
DeBoer et al., 2017), the upgraded MWA (MWA-II; Wayth et al., 2018), and the
LOFAR 2.0 Survey (Edler et al., 2021). With more collecting areas and improved
electronics, high-significant measurements of the 21-cm power spectrum are
anticipated from these second-generation instruments (Pober et al., 2014; Liu
& Parsons, 2016; Wayth et al., 2018).

Detecting the 21-cm signal from EoR is challenging due to the weakness of
21-cm signal comparing to the brightness temperature of the foreground
originated from synchrotron radiation from our Galaxy, strong point sources
and thermal Bremsstrahlung from the HII region (Shaver et al., 1999; Santos et
al., 2005; Furlanetto et al., 2006; Bernardi et al., 2009; Parsons et al.,
2014). To disentangle strong foregrounds from weak 21-cm signal, several
mitigation techniques have been developed (See e.g., Wang et al., 2006; Liu et
al., 2009; Bowman et al., 2009; Liu & Tegmark, 2011; Parsons et al., 2012;
Dillon et al., 2013; Liu et al., 2014a; Beardsley et al., 2016; Pober et al.,
2016; Hothi et al., 2021). These techniques, in general, can be categorised
into two approaches – foreground subtraction and foreground avoidance (Chapman
& Jelić, 2019). Foreground subtraction involves modelling of the foreground
and subtracting them from the data. In contrast, foreground avoidance
completely discards foreground contaminated data from analysis and try to
reconstruct the power spectra of EoR in the Fourier regime which are not
affected by foreground.

In both foreground mitigation and subtraction, accurate calibration of the
array is critical because calibration errors can contaminate the Fourier mode
of the EoR power spectrum (Morales et al., 2012; Ewall-Wice et al., 2017).
Furthermore, simulations in Orosz et al. (2018) suggested that calibration
error introduced by quasi-redundancy in antenna positioning and slight
variations of the beam responses between different antennas across the array
are significant sources of contamination of the EoR signal. Even in the limit
of perfect antenna positioning and telescope response, the EoR window can
still be contaminated if the instrument calibration is susceptible to sky
modelling error (Byrne et al., 2019).

The array configuration of PAPER is designed to enhance sensitivity for a
first power-spectrum detection, by measuring same Fourier modes on the sky
with a group of redundant baselines. These redundant visibilities need to be
calibrated accurately by using the redundancy information of the array. The
calibration technique implemented for such kind of array layout is generally
known as redundant baseline calibration (Wieringa, 1992; Liu et al., 2010).
This technique intends to calibrate measurements from the array of antennas
configured on a regular grid spacing with maximum redundancy in baseline
distribution and tightly packed with the drift-scan mode of observations. Then
one can use statistical tool to calibrate the visibilities for independent
baselines.

In this work, we introduce a recently developed calibration scheme called
_correlation calibration_ (CorrCal111The prototype version of the software
located on https://github.com/sievers/corrcal2.git is being used for this
work.), first proposed in Sievers (2017). CorrCal aims to bridge the gap
between sky-based and redundant calibration by incorporating partial sky
information and information on the non-redundancies of the array into the
calibration solutions. We use this new scheme to re-calibrate data from the
64-element PAPER array (hereafter PAPER64), which is a 9-hour integration data
set. Ali et al. (2015) calibrated this data set with a redundant baseline
calibration scheme and published power spectrum results. This work used the
same initial calibration techniques employed in Ali et al. (2015), but have an
independent pipeline from that point on. We will show that CorrCal can
potentially constrain the wedge-like structure of foreground in the Fourier
space.

The rest of paper is organized as follows. In Section 2, we give a general
review of common calibration approaches for radio interferometric
measurements. In Section 3, we present CorrCal calibration scheme. In Section
4 we apply our method to the 9-hour PAPER-64 data set. In Section 5, we
present the results of the CorrCal calibration. The conclusion will be in the
last section.

##  2 A Review of Calibration Formalism

In radio interferometry, the correlated signal (visibility) from two
independent receiving elements (antennas) is the sum of the system noise and
the actual sky signal (or true sky visibility) multiplied by the complex gain
parameters of the antennas. The gain parameters depend on the electromagnetic
properties of the antennas, whereas the additive noise is mostly dominated by
the thermal noise of the receiver system and the sky components.

Performing a calibration is to solve for the complex gain parameters to
recover the true sky signal from the measured signal. The mathematical form of
this statement is the measurement equation,

$\displaystyle d_{pq}=g_{p}\,g^{\ast}_{q}\,v_{pq}+\,n_{pq}.$ (1)

Here, $g_{p}$ and $g_{q}$ are the unknown complex gain parameters of antenna
$p$ and antenna $q$, which depend on the frequency $f$. $v_{pq}$ is the true
sky visibility corresponding to the baseline between antenna $p$ and antenna
$q$, for which we ultimately want to recover from the measured visibility
(observed data) $d_{pq}$. $n_{pq}$ is the system noise on this baseline, which
is usually assumed to be Gaussian. The asterisk superscript ∗ indicates the
complex conjugate. All terms in Eq. (1) are also implicitly time-dependent.
The antenna gain $g$ can also be written as a complex exponential form with
amplitude $a$ and phase $\phi$ for each of the antenna $p$ (or $q$),

$\displaystyle g_{p}=a_{p}\exp\left(\mathrm{i}\mkern 1.0mu\phi_{p}\right).$
(2)

Note that we used “$\mathrm{i}\mkern 1.0mu$” to indicate the imaginary unit
$\sqrt{-1}$ throughout the paper.

For multiple baselines, the measurement equation can be written in a more
compact form by using vector and matrix notations

$\displaystyle\textbf{{d}}=\textbf{{G}}\textbf{{v}}+\textbf{{n}}.$ (3)

Here, d is a vector containing the measured visibilities of different
baselines. v is vector of the corresponding true sky visibility. G is a
diagonal gain matrix with $g_{p}g^{\ast}_{q}$ as its elements placed
diagonally, and n is a vector containing noise of the corresponding baseline
formed from antenna $p$ and $q$. For simplicity, we will omit the frequency
(and time) dependency in the vector and matrix notations through the rest of
the paper unless otherwise stated.

Most calibration software packages developed so far to solve for the gain
parameter G can be categorized into either one of the two conventional
approaches _sky-based_ calibration or _redundant-baseline_ calibration222It is
worth noticing that Liu et al. (2010) investigated the Taylor expansion of the
sky to deal with direction-dependence or non-redundancy.. In the sub-sections
that follow, we present the underlying mathematical formalism of these two
calibration approaches and discuss their inherent advantages and
disadvantages.

###  2.1 Sky-based Calibration

Sky-based calibration uses a model visibility that approximates the true sky
visibility, $m_{pq}\approx v_{pq}$ to solve for the gain solutions. Eq. (1)
can be written equivalently in the matrix form,

$\displaystyle\textbf{{d}}=\textbf{{G}}\textbf{{m}}+\textbf{{n}}.$ (4)

The model visibility is usually constructed from known sources and the
telescope beam response. Solving for the gains G is then essentially
equivalent to minimizing the $\chi^{2}$ function of the difference between the
measured visibility d and the model visibility m,

$\displaystyle\chi^{2}=(\textbf{{d}}-\textbf{{G}}\textbf{{m}})^{\dagger}\textbf{{N}}^{-1}(\textbf{{d}}-\textbf{{G}}\textbf{{m}}).$
(5)

Here, N is an ensemble average noise matrix in which its diagonal elements are
the noise variance $\sigma_{pq}^{2}$ on the baseline formed from antenna $p$
and antenna $q$,

N $\displaystyle=$
$\displaystyle\langle\textbf{{n}}\textbf{{n}}^{\dagger}\rangle$ (6)
$\displaystyle=$ $\displaystyle\sigma_{pq}^{2}\textbf{{I}},$

where the angle brackets $\langle...\rangle$ denote ensemble average, and in
the second equality in Eq. (6) assumes that the noise matrix is diagonal,
i.e., the noise on visibilities of different baselines are uncorrelated. I is
the identity matrix, and $\dagger$ denotes a complex conjugate transpose. The
measured visibility can be corrected (or calibrated) by using the gain
solutions estimated from optimization of Eq. (5). However, the precision of
sky-based calibration is highly dependent on the sky models, which are
imperfect due to missing sources or an inaccurate beam model. This
imperfection can lead to calibration errors that introduce foreground
contamination in the cosmological Fourier space for 21 cm signal (Beardsley et
al., 2016; Barry et al., 2016; Ewall-Wice et al., 2017; Byrne et al., 2019).

###  2.2 Redundant Calibration

Redundant calibration relies on multiple copies of identical baselines in
arrays with the regular layout, to relatively solve for both the gain
parameter and the true sky visibility from an over-determined system of linear
equations without the need for the sky model (Wieringa, 1992). The measurement
equation for redundant calibration replaces the true sky visibilities $v_{pq}$
with visibility terms that are constrained to be equal across redundant
baselines $u_{\alpha}$, where $\alpha$ indexes the redundant baseline sets
(redundant groups),

$\displaystyle d_{pq}=g_{p}g^{\ast}_{q}u_{\alpha}+n_{pq}.$ (7)

Since the number of measured visibility $d_{pq}$ are greater than the number
of $u_{\alpha}$ that averaged over redundant group of the array and the gain
parameter $g$, this equation becomes over-determined and can be solved,
provided there are sufficient redundant baselines. The solutions for Eq. (7)
can be then determined by minimizing the $\chi^{2}$ function of the form,

$\displaystyle\chi^{2}=(\textbf{{d}}-\textbf{{G}}\textbf{{u}}_{\alpha})^{\dagger}\textbf{{N}}^{-1}(\textbf{{d}}-\textbf{{G}}\textbf{{u}}_{\alpha}).$
(8)

In contrast to Eq. (5), the redundant visibility parameter
$\textbf{{u}}_{\alpha}$ here is not pre-defined and will be solved iteratively
alongside with the antenna gain factor $\textbf{{G}}\equiv g_{p}g_{q}$. To
minimise the $\chi^{2}$ function, a least square estimator is usually
utilised. We direct the interested reader to the existing literature (e.g.,
Liu et al., 2010; Zheng et al., 2014; Dillon et al., 2018, 2020) for detailed
discussion of least square solver implementation. Minimizing Eq. (8) gives us
a power spectrum estimation

$\hat{\textbf{{x}}}=[\textbf{{A}}^{\intercal}\textbf{{N}}^{-1}\textbf{{A}}]^{-1}\textbf{{A}}^{\intercal}\textbf{{N}}^{-1}\textbf{{d}}.$
(9)

Here, the matrix A, which depends on the array configuration, maps the
measured visibility d to the calibration parameters G and v for which we want
to solve. The vector $\hat{\textbf{{x}}}$ contains the least-square estimates
power spectrum estimates for per antenna gain parameter, and per-redundant
group averaged visibility. Eq. (9) is solvable iteratively, using the fitted
calibration parameter as the fiducial guesses for the next cycle of the
fitting, until the solutions converge. The method is implemented in the
OmniCal333https://github.com/jeffzhen/omnical.git software package (Zheng et
al., 2014), which was used to calibrate PAPER and the MWA Phase II
observations (Parsons et al., 2014; Ali et al., 2015; Li et al., 2018;
Kolopanis et al., 2019).

Due to the _degeneracies_ in the measurement equation, Eq. (7), the least
square estimator of Eq. (8) will not converge to a unique solution. Here, the
term degeneracies refer to the linear combinations of gains and visibilities
that redundant calibration cannot solve for, and which results in the
$\chi^{2}$ function is invariant across different choices of the degenerate
parameters. In-depth discussions of degeneracy issues in redundant baseline
calibration can be found in (Liu et al., 2010; Zheng et al., 2014; Dillon et
al., 2018; Byrne et al., 2019).

In general, redundant baseline calibration suffers from four types of
degeneracies for array elements configured on the same plane. First, if one
multiplies the amplitude of the antenna gain by a constant factor $A$ and
divides the sky parameter by $A^{2}$, the $\chi^{2}$ remains unchanged so that
the overall amplitude, $A$, of the instrument can not be solved. Second,
moving the phase of the gain parameter by certain amount $\Delta$, i.e.,
$g_{p}=|g_{p}|e^{j\phi_{p}}\rightarrow|g_{p}|e^{j(\phi_{p}+\Delta)}$ does not
affect the $\chi^{2}$ due to cancellation effect of $g_{p}\times g^{*}_{q}$ in
Eq. (1) so that the redundant baseline calibration still can not be set the
overall phase degeneracy $\Delta$. The third and fourth degeneracies are
caused by phase gradients $\Delta_{x}$ and $\Delta_{y}$ along $x$ and $y$
direction, respectively, for a planar array. For example, if one shifts
$g_{p}=|g_{p}|e^{j\phi_{p}}\rightarrow|g_{q}|e^{j(\phi_{q}+\Delta_{x}x_{q}+\Delta_{y}y_{q})}$
it does not affect the $\chi^{2}$ if it is complemented by shift in sky
parameter
$V_{pq}^{\rm{true}}=|V_{pq}^{\rm{true}}|e^{j\phi_{pq}}\rightarrow|V_{pq}^{\rm{true}}|e^{j(\phi_{pq}-\Delta_{x}b_{x}-\Delta_{y}b_{y})}$,
where $x_{q}$ and $y_{q}$ are the coordinates for antenna $j$, and
$b_{x}=x_{q}-x_{p}$ and $b_{y}=y_{q}-y_{p}$ are $x$ and $y$ coordinates of a
baseline vector b connecting antenna $p$ and $q$, respectively. So, the phase
gradients $\Delta_{x}$ and $\Delta_{y}$ are still degenerate in the relative
calibration stage. In a sense, this means that inclining the whole array in
either direction (x or y) is exactly equivalent to moving the phase centre of
the sources on the sky in opposite direction. Hence, these two effects can
cancel each other and result in the phase gradient degeneracies, which in turn
results in the sources appearing to be offset from the phase centre.

To fix the degeneracy issues, the redundant calibration can be divided into
two calibration steps. These are relative calibration and absolute
calibration. First, relative calibration is performed using the algorithm that
we have just described to solve the antenna gains and the true sky. However,
these solutions are degenerate, and an overall amplitude, an overall phase,
and phase gradient calibration parameters are remain unknown. Then, absolute
calibration is performed to set these degeneracies to a reference and obtain
the final calibration solutions. One particular approach for absolution
calibration, as discussed in Byrne et al. (2019), is to map a pre-defined sky
model of bright point sources to the degenerate parameters. However, setting
degeneracies using the sky information in redundant calibration is difficult
because incompleteness of the sky model can introduce frequency-dependent gain
errors, which may lead to spectral structures in the otherwise smooth
observations (Barry et al., 2016; Byrne et al., 2019).

Thus, the fundamental shortcomings of redundant calibration are imperfection
of array redundancy and sky-modeling error. During the absolute calibration
stage, the sky-modeling error enters into the redundant calibration. Redundant
calibration assumes that the same sky is seen by all baselines in the same
redundant group, requiring these baselines to have the same physical length
and orientation. This assumption also dictates that all antenna elements in
the array have the same primary beam responses. In the real situation, both
are impossible to achieve for the level of accuracy required. These non-
redundancies will result in calibration errors that lead to more foreground
contamination in the cosmological Fourier space for 21 cm signal (Ewall-Wice
et al., 2017; Orosz et al., 2018). The powerful alternative is to include the
relaxed-redundancy information and position of sources into the calibration
formalism (Sievers, 2017).

##  3 Correlation Calibration (CorrCal)

Sievers (2017) proposed a hybrid calibration scheme called CorrCal, which aims
to further bridge the gap between sky-based and redundant calibrations by
taking into account sky and array information into the calibration algorithm.
The CorrCal framework relies on the assumption that sky information is
statistically Gaussian. While the real sky is never exactly Gaussian, Sievers
(2017) argued that the assumption still provide a good approximation of the
reality. This assumption allows CorrCal to relax the redundancy of the array
and including the known sky information in its formalism through covariance
based calculation of $\chi^{2}$. The covariance matrix thus forms a
_statistical model_ of the sky power spectrum, weighted by the instrumental
response. We describe the underlying mathematical formalism of CorrCal in this
section.

Given a vector of the measured visibilities d as is defined in Eq. (3), its
covariance $\boldsymbol{\Sigma}$ follows,

$\displaystyle\boldsymbol{\Sigma}=\langle\textbf{{d}}\textbf{{d}}^{\dagger}\rangle.$
(10)

Using the definition of d from Eq. (3), $\boldsymbol{\Sigma}$ can be expanded
as,

$\displaystyle\boldsymbol{\Sigma}$ $\displaystyle=$
$\displaystyle\langle(\textbf{{G}}\textbf{{v}}+\textbf{{n}})(\textbf{{G}}\textbf{{v}}+\textbf{{n}})^{\dagger}\rangle$
(11) $\displaystyle=$
$\displaystyle\textbf{{G}}\langle\textbf{{v}}\textbf{{v}}^{\dagger}\rangle\textbf{{G}}^{\dagger}+\textbf{{G}}\langle\textbf{{v}}\textbf{{n}}^{\dagger}\rangle+\langle\textbf{{n}}\textbf{{v}}^{\dagger}\rangle\textbf{{G}}^{\dagger}+\langle\textbf{{n}}\textbf{{n}}^{\dagger}\rangle.$

By assuming that the instrument noise vector n does not correlate with the
true sky visibility vector v, the middle two terms drop out. The last term
describes the instrument noise, which is generally assumed to be Gaussian and
uncorrelated between baselines; therefore, it can be written as a noise matrix
N as defined in Eq. (6). The first term describes the correlation of the sky
multiplied by the gains. If we define the variance-covariance matrix of the
true sky as

$\displaystyle\textbf{{C}}=\langle\textbf{{v}}\textbf{{v}}^{\dagger}\rangle,$
(12)

Eq. (11) takes

$\displaystyle\boldsymbol{\Sigma}=\textbf{{G}}\textbf{{C}}\textbf{{G}}^{\dagger}+\textbf{{N}}.$
(13)

Under an assumption that the sky is random Gaussian distributed, we note that
by virtue of the central limit theorem the multitude of random phenomena that
produce the random character of the observed data d, implies their
distributions are nearly Gaussian in general. Thus, the likelihood
$\mathcal{L}$ of data d given their covariance matrix $\boldsymbol{\Sigma}$
takes the standard form of a multi-dimensional Gaussian distribution,

$\displaystyle\mathcal{L}(\textbf{{d}},\boldsymbol{\Sigma})\propto\exp\left[-\frac{1}{2}\textbf{{d}}^{\dagger}\boldsymbol{\Sigma}^{-1}\textbf{{d}}\right].$
(14)

Using Eq. (13) in this equation, we obtain the likelihood

$\displaystyle\mathcal{L}(\textbf{{d}},\textbf{{C}}|\textbf{{G}})\propto\exp\left[-\frac{1}{2}\textbf{{d}}^{\dagger}\left(\textbf{{G}}\textbf{{C}}\textbf{{G}}^{\dagger}+\textbf{{N}}\right)^{-1}\textbf{{d}}\right].$
(15)

Taking the logarithm of both sides of Eq. (15), we obtain the negative log-
likelihood function $\log\mathcal{L}$,

$\displaystyle\log\mathcal{L}(\textbf{{d}},\textbf{{C}}|\textbf{{G}})\propto-\frac{1}{2}\textbf{{d}}^{\dagger}\left(\textbf{{G}}\textbf{{C}}\textbf{{G}}^{\dagger}+\textbf{{N}}\right)^{-1}\textbf{{d}}.$
(16)

Maximising the likelihood function $\mathcal{L}$ is equivalent to minimising
the negative log-likelihood $\log\mathcal{L}$. The term on the right-hand side
of Eq. (16) measures the square of the distances from predicted (or true)
visibility to the measured visibility in the standard deviation unit.
Minimising this term against G is equal to maximising the likelihood of true
visibility being observed. The exponential term is thus simply equal to the
$\chi^{2}$, given by

$\displaystyle\chi^{2}=\textbf{{d}}^{\dagger}\left(\textbf{{G}}\textbf{{C}}\textbf{{G}}^{\dagger}+\textbf{{N}}\right)^{-1}\textbf{{d}}.$
(17)

The covariance based calculation of $\chi^{2}$ in this equation is the
underlying mathematical formalism of CorrCal. We further add an additional
term $\chi^{2}\rightarrow\chi^{2}+(\sum_{i}{\rm Im}(g_{i}))^{2}+(\sum_{i}{\rm
Re}(g_{i})-N_{\rm ant})^{2}$ as a regularisation factor to bound the gain
solution. Here, the $\chi^{2}$ function does not explicitly depends on the sky
parameter or written in terms of the redundant sets of visibility, but, still,
the array and sky information are integrated into the covariance matrix C of
the true sky visibility.

Since we are modeling the data as being Gaussian distributed, a form of the
$\chi^{2}$ function in Eq. (17) is different from the one defined in Eq. (5)
and Eq. (8) which explicitly depends on the sky parameters m and
$\textbf{{u}}_{\alpha}$, respectively. However, the CorrCal approach of
keeping the array and sky information in the C matrix allows one to reproduce
an exact copy of the redundant calibration algorithm in a limiting case. For
instance, all baselines within the perfectly redundant group see the same sky,
so the effective noise $\mathbf{N}_{\rm eff}$ for this group is written as the
sum of per-visibility diagonal noise variance matrix N and covariance between
visibilities. That is,

$\displaystyle\mathbf{N}_{\rm
eff}=\textbf{{N}}+(a\mathbf{1})\otimes(a\mathbf{1})^{\dagger},$ (18)

where $a$ is a parameter controlling the variance of the sky signal,
$\mathbf{1}=[1,1,1,\dots]$ is a vector of ones, and an operator $\otimes$
denotes the outer product of vectors. Applying the Woodbury identity (see Eq.
(24)) to $\mathbf{N}_{\rm eff}^{-1}$, and taking $a\to\infty$, the $\chi^{2}$
becomes the same to one for redundant calibration algorithm.

Although complete information of the true sky can never be obtained, we can
incorporate partial sky information into the covariance matrix by using known
sky models. Thus, we further split the covariance matrix in Eq. (17) into two
components as per suggestion in Sievers (2017),

$\displaystyle\textbf{{C}}=\textbf{{S}}+\textbf{{R}},$ (19)

where the matrix S contains the covariance of the expected visibilities of
known sky sources, whereas R is the covariance matrix of the visibilities
within the redundant group that contains everything else that is not in S,
i.e. ‘the rest,’ which is mostly the diffuse sky signal from our Galaxy’s
emission, weighted by the instrumental response. Note that Eq. (19) assumes
that the sky components in S and R are uncorrelated since most point sources
are extragalactic. With this definition, Eq. (17) becomes

$\displaystyle\chi^{2}=\textbf{{d}}^{\dagger}\left(\textbf{{G}}\big{(}\textbf{{S}}+\textbf{{R}})\textbf{{G}}^{\dagger}+\textbf{{N}}\right)^{-1}\textbf{{d}}.$
(20)

Note that the explicit formalisms of S and R will be presented in Section 3.1
in details.

With the redundant group diffuse sky R and sources S expressed as vector outer
products, Eq. (19) can be rewritten as

$\displaystyle\textbf{{C}}=\textbf{{s}}\textbf{{s}}^{\dagger}+\textbf{{r}}\textbf{{r}}^{\dagger},$
(21)

where r and s are $N$-dimensional vectors for the diffuse sky and source
components, respectively. We substitute Eq. (21) for C in Eq. (20) to obtain a
$\chi^{2}$ function

$\displaystyle\chi^{2}=\textbf{{d}}^{\dagger}\left(\textbf{{N}}+\textbf{{G}}\big{(}\textbf{{s}}\textbf{{s}}^{\dagger}+\textbf{{r}}\textbf{{r}}^{\dagger}\big{)}\textbf{{G}}^{\dagger}\right)^{-1}\textbf{{d}}.$
(22)

Applying the gain matrix to the source and redundant vectors in this equation,
we have

$\displaystyle\chi^{2}=\textbf{{d}}^{\dagger}\Big{(}\textbf{{N}}+\hat{\textbf{{s}}}\hat{\textbf{{s}}}^{\dagger}+\hat{\textbf{{r}}}\hat{\textbf{{r}}}^{\dagger}\Big{)}^{-1}\textbf{{d}},$
(23)

where
$\hat{\textbf{{s}}}=\textbf{{G}}\textbf{{s}},\hat{\textbf{{r}}}=\textbf{{G}}\textbf{{r}}$,
so $\hat{\textbf{{s}}}^{\dagger}=\textbf{{s}}^{\dagger}\textbf{{G}}^{\dagger}$
and
$\hat{\textbf{{r}}}^{\dagger}=\textbf{{r}}^{\dagger}\textbf{{G}}^{\dagger}$.

Inversion of the matrices inside the parenthesis term in Eq. (23) can be done
with Woodbury inversion formula (Woodbury, 1950), which we write here in the
special case of a Hermitian matrix if both $\textbf{{B}}^{-1}$ and
$\textbf{{I}}+\textbf{{v}}^{\dagger}\textbf{{B}}^{-1}\textbf{{v}}$ are
invertible

$\displaystyle\Big{(}\textbf{{B}}+\textbf{{v}}\otimes\textbf{{v}}^{\dagger}\Big{)}^{-1}=\textbf{{B}}^{-1}-\textbf{{B}}^{-1}\textbf{{v}}\Big{(}\textbf{{I}}+\textbf{{v}}^{\dagger}\textbf{{B}}^{-1}\textbf{{v}}\Big{)}^{-1}\textbf{{v}}^{\dagger}\textbf{{B}}.$
(24)

Where B and the identity matrix I are square matrices with dimension $N\times
N$, while v is $N$-dimensional vector.

For the sake of computational tractability, the implementation of the
Woodburry inversion to Eq. (23) in CorrCal follows the following steps:

  * •

First, the Woodburry inversion can be applied to $\Gamma^{-1}$ where
$\Gamma=\big{(}\textbf{{N}}+\hat{\textbf{{r}}}\hat{\textbf{{r}}}^{\dagger}\big{)}$
for each redundant group separately, ignoring the source vectors s. To apply
Woodburry identity to $\Gamma^{-1}$, let $\textbf{{B}}=\textbf{{N}}$ and
$\textbf{{v}}=\textbf{{r}}$. Then, the inverted form of $\Gamma$ is kept
separate for the each redundant group, and saved in a factored form by
applying the Cholesky factorization of the matrix
$\big{[}\textbf{{I}}+\textbf{{v}}^{\dagger}\textbf{{B}}^{-1}\textbf{{v}}\big{]}^{-1}$
in Eq. (24) for the redundant vectors r.

  * •

Finally, using the source vectors s, and the already factored form of
$\Gamma^{-1}$ from the first step the Woodburry inversion can be computed
using
$\big{(}\Gamma+\hat{\textbf{{s}}}\hat{\textbf{{s}}}^{\dagger}\big{)}^{-1}.$

###  3.1 Determining the Sky and Noise Variance-Covariance Matrices

To perform the CorrCal, we first need to estimate the matrices R, S, and N.
Thus, in this section, we will derive these variance-covariance matrices
explicitly.

####  3.1.1 Diffuse Sky Component Matrix (R)

For an instrument with a small field of view, the covariance matrices can be
directly constructed from the visibility in the 2D UV plane. However, most
21-cm experiments, including PAPER, use receiving elements that have wide
field-of-views. Thus, the curvature of the sky becomes important. To account
for this effect and simplify the calculation, we will first project the
visibility function into the spherical harmonic space before calculating the
covariance.

Given a vector $\textbf{{u}}=\textbf{{b}}/\lambda$ that expresses the baseline
vector b between antennas $p$ and $q$ in wavelength $\lambda$, the measured
cross-correlation signal from these antennas, or the visibility, ignoring
noise term, is given by,

$\displaystyle
v_{pq}(\textbf{{u}}|f)=\int\textrm{d}\Omega\,I(\hat{\textbf{{r}}},f)B(\hat{\textbf{{r}}},f),$
(25)

where

$\displaystyle B(\hat{\textbf{{r}}},f)$ $\displaystyle=$ $\displaystyle
A_{p}(\hat{\textbf{{r}}},f)A^{\ast}_{q}(\hat{\textbf{{r}}},f)\exp\left[-2\mathrm{i}\mkern
1.0mu\pi\textbf{{u}}\cdot\hat{\textbf{{r}}}\right]$ (26) $\displaystyle=$
$\displaystyle|A_{p}(\hat{\textbf{{r}}},f)|^{2}\exp\left[-2\mathrm{i}\mkern
1.0mu\pi\textbf{{u}}\cdot\hat{\textbf{{r}}}\right].$

In this equation, $I(\hat{\textbf{{r}}},f)$ is the flux density measured over
a solid angle $\Omega$ subtended by the observed region of the sky in the
direction of a unit vector $\hat{\textbf{{r}}}$. $A_{p}(\hat{\textbf{{r}}},f)$
and $A_{q}(\hat{\textbf{{r}}},f)$ are the primary beam of antenna $p$ and $q$,
repectively. Here, the squaring of the antenna beam term
$A_{p}(\hat{\textbf{{r}}},f)$ in the second equality of Eq. (26) comes from
the fact that we assume all antennas have the same primary beam response. The
antenna position information is integrated into the exponential factor. Thus,
the second line of the equation carries array information (position and beam
response).

Expanding the terms for the sky ($I(\hat{\textbf{{r}}},f)$) and array
($B(\hat{\textbf{{r}}},f)$) over the spherical harmonics, we have

$\displaystyle I(\hat{\textbf{{r}}},f)=\sum_{\ell m}a_{\ell m}(f)Y_{\ell
m}(\hat{\textbf{{r}}}),$ (27) $\displaystyle
B(\hat{\textbf{{r}}},f)=\sum_{\ell m}b_{\ell m}(f)Y_{\ell
m}(\hat{\textbf{{r}}}),$ (28)

where $a_{\ell m}$ and $b_{\ell m}$ are the amplitudes of the flux density and
the beam function for the spherical harmonics $Y_{\ell m}$ respectively. The
angular number $\ell=0,1,2,3,\dots$, and the azimuthal number
$m=-\ell,-\ell+1,\cdots,\ell-1,\ell$ follow standard convention, allowing
$2\ell+1$ values of $m$ for each value of $\ell$. Using Eq. (27) and Eq. (28),
the visibility in Eq. (25) can be expressed in spherical harmonics as,

$\displaystyle v_{pq}(f)$ $\displaystyle=$
$\displaystyle\int\textrm{d}\Omega\,\sum_{\ell m}a^{\ast}_{\ell
m}(f)Y^{\ast}_{\ell
m}(\hat{\textbf{{r}}})\sum_{\ell^{\prime}m^{\prime}}b_{\ell^{\prime}m^{\prime}}(f)Y_{\ell^{\prime}m^{\prime}}(\hat{\textbf{{r}}})$
(29) $\displaystyle=$
$\displaystyle\sum_{\ell\ell^{\prime}mm^{\prime}}a^{\ast}_{\ell
m}(f)b_{\ell^{\prime}m^{\prime}}(f)\int\textrm{d}\Omega\,Y^{\ast}_{\ell
m}(\hat{\textbf{{r}}})Y_{\ell^{\prime}m^{\prime}}(\hat{\textbf{{r}}})$
$\displaystyle=$ $\displaystyle\sum_{\ell m}a^{\ast}_{\ell m}(f)b_{\ell
m}(f),$

where we have used the orthogonality of spherical harmonics to turn the
integral in the second line into a Kronecker delta function. Eq. (29) is the
spherical harmonics version of the visibility equation in the absence of
noise.

We calculated the array function that is defined in Eq. (26) using the
exponential argument ($\textbf{{b}}\cdot\hat{\textbf{{r}}}{}/c$) in Eq. (25)
and the PAPER-64 beam model depicted in Fig. 2. The simulated exponential
argument is shown in Fig. 1 for the PAPER-64 baselines oriented along the
East-West and the North-South directions. To generate the beam model of the
PAPER-64 array within the frequency range of 120 MHz to 168 MHz, we use the
Astronomical Interferometry in Python (AIPY)444https://github.com/HERA-
Team/aipy.git, a software package primarily developed for PAPER/HERA data
analysis. Fig. 2 shows the AIPY simulated PAPER-64 beam model projected on the
HEALPIx555https://github.com/healpy/healpy.git coordinates at three arbitrary
right accessions for frequency of 150 MHz. We will use this beam model for all
calculations of the covariance matrices.

Figure 1: A delay term ($\textbf{{b}}\cdot\hat{\mathbf{r}}/c$) pattern in
nanoseconds (ns) over PAPER latitude for baselines pointed along the East-West
(left) and the North-South (right) directions. These baselines are calculated
from the actual PAPER-64 antenna coordinates that shown in Fig. 8. Figure 2:
The AIPY simulated antenna beam pattern of PAPER-64 array at 150 MHz at
different longitudes, plotted for the PAPER latitude ($30.7^{\circ}$S) at
Karoo desert, South Africa.

Using the expected visibility function defined in Eq. (29), in this section we
construct the diffuse sky covariance matrix R. While determining this matrix,
all degrees of freedom including frequencies and baselines must be taken into
accounts. We estimate this matrix from the relation

$\displaystyle\textbf{{R}}_{({\textbf{{u}}}f);({\textbf{{u}}}^{\prime}f^{\prime})}$
$\displaystyle=$ $\displaystyle\left\langle v^{\rm
diff}(\textbf{{u}},f)v^{{\rm
diff},*}(\textbf{{u}}^{\prime},f^{\prime})\right\rangle$ (30) $\displaystyle=$
$\displaystyle\left\langle\sum_{\ell m}a^{{\rm diff},\ast}_{\ell m}(f)b^{{\rm
diff}}_{\ell m}({\textbf{{u}}},f)\right.$ $\displaystyle\times$
$\displaystyle\left.\sum_{\ell^{\prime}\mathrm{m}^{\prime}}a^{{\rm
diff}}_{\ell^{\prime}\mathrm{m}^{\prime}}(f^{\prime})b^{{\rm
diff},\ast}_{\ell^{\prime}\mathrm{m}^{\prime}}({\textbf{{u}}}^{\prime},f^{\prime}))\right\rangle$
$\displaystyle=$ $\displaystyle\sum_{\ell
m}\sum_{\ell^{\prime}\mathrm{m}^{\prime}}b^{{\rm diff}}_{\ell
m}({\textbf{{u}}},f)\left\langle a^{{\rm diff},\ast}_{\ell m}(f)a^{{\rm
diff}}_{\ell^{\prime}\mathrm{m}^{\prime}}(f^{\prime})\right\rangle$
$\displaystyle\times$ $\displaystyle b^{{\rm
diff},\ast}_{\ell^{\prime}\mathrm{m}^{\prime}}({\textbf{{u}}}^{\prime},f^{\prime}),$

where "diff" is to denote $\textbf{{v}}_{pq}(f)$ is coming from the diffuse
sky component and the angle bracket $\langle...\rangle$ shows the averaging
over an ensemble of realizations of the fluctuation. With the assumption that
the sky is statistically isotropic Gaussian random field, the covariance for
sky flux density spherical harmonics coefficients becomes

$\displaystyle\left\langle a^{{\rm diff}}_{\ell
m}(f)a_{\ell^{\prime}\mathrm{m}^{\prime}}^{{\rm
diff},*}(f^{\prime})\right\rangle$ $\displaystyle=\left\langle|a_{\ell
m}(f,f^{\prime})|^{2}\right\rangle\delta_{\ell\ell^{\prime}}\delta_{\mathrm{mm}^{\prime}}$
(31) $\displaystyle=C^{{\rm
diff}}_{\ell}(f,f^{\prime})\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}},$

where $C^{{\rm diff}}_{\ell}(f,f^{\prime})=\langle|a^{{\rm diff}}_{\ell
m}(f,f^{\prime})|^{2}\rangle$ is the the diffuse sky angular power spectrum,
which only depends on multipole number $\ell$ that corresponds to the angular
scale $\theta=180^{\circ}/\ell$, and $\delta$ is the Kronecker delta-function.
This is because the assumption of isotropy ensures that $\langle|a^{{\rm
diff}}_{\ell m}(f,f^{\prime})|^{2}\rangle$ is a function of only $\ell$, not
$m$. Therefore, there is no correlation in the $m$ modes. Thus, the covariance
matrix (Eq. (30)) can be simplified to

$\displaystyle\textbf{{R}}=\sum_{\ell m}b^{{\rm diff}}_{\ell
m}({\textbf{{u}}},f)b^{{\rm diff},*}_{\ell
m}({\textbf{{u}}}^{\prime},f^{\prime})C^{{\rm diff}}_{\ell}(f,f^{\prime}).$
(32)

Equation (32) is a general expression for the covariance matrix corresponding
to the diffuse sky. It depends on the baseline and antenna beam model through
a function defined in Eq. (26). The measured antenna positions (instead of
assuming they are perfectly redundant) are integrated into this equation.
Hence, through this expression, the covariance-based calculation of $\chi^{2}$
in Eq. (20) carries array information during optimization processes. Fig 3
displays the covariance matrix between visibilities at 150 MHz in a given
redundant group of baselines whose length about $30$ m.

Figure 3: Covariance matrix between baselines obtained from the diffuse sky
component at $150\,{\rm MHz}$ for a redundant group containing about $30$ m
length baselines. The axes of this figure spanning the true measured variation
in the baseline lengths for the baselines in the 30m group. The entries of the
covariance are not precisely equal as displayed in the figure. The covariance
between baselines drops off as baselines become less redundant.

In Eq. (32), however, the term $C^{{\rm diff}}_{\ell}(f,f^{\prime})$ is not
known. The measured data covariance matrix should be constructed from the same
redundant group that R had formed to estimate $C^{{\rm
diff}}_{\ell}(f,f^{\prime})$. The data covariance matrix can be determined
from

$\displaystyle\textbf{{D}}^{\alpha}=\langle\textbf{{d}}^{\alpha}({\textbf{{u}}},f)\textbf{{d}}^{\alpha\dagger}({\textbf{{u}}}^{\prime},f^{\prime})\rangle,$
(33)

where $\textbf{{D}}^{\alpha}$ is a covariance matrix of measured visibility
data $\textbf{{d}}^{\alpha}$ from the redundant group $\alpha$. To determine
the $C^{\rm diff}_{\ell}(f,f^{\prime})$, thus, we equate the covariance matrix
$\textbf{{D}}^{\alpha}$ with the diffuse sky matrix R from the same redundant
group $\alpha$ in Eq. (32). Note that the diffuse galactic emissions dominate
visibilities from shorter baselines, while visibilities measured by longer
baselines are unlikely to be dominated by diffuse emissions, so we only equate
D with R (i.e., $\textbf{{R}}^{\alpha}\approx\textbf{{D}}^{\alpha}$) if both
are computed from the same redundant group $\alpha$. This assumption may avoid
mode-mixing between baselines with different lengths. It follows that,

$\displaystyle C^{{\rm
diff}}_{\ell}(f,f^{\prime})\approx\frac{\textbf{{D}}^{\alpha}}{\sum_{\ell
m}b^{{\rm diff}}_{\ell m}({\textbf{{u}}},f)b^{{\rm diff}}_{\ell
m}({\textbf{{u}}}^{\prime},f^{\prime})}.$ (34)

Figure 4: Left–Cross-frequency sky power spectrum $C_{\ell}^{{\rm
diff}}(f,f^{\prime})$ obtained from Eq. (34). Middle–the fitted power spectrum
$C_{\ell,{\rm fits}}^{{\rm diff}}(f,f^{\prime})$. Right–the residuals $C^{\rm
diff}_{\ell,\rm resd}(f,f^{\prime})=C_{\ell}^{{\rm
diff}}(f,f^{\prime})-C_{\ell,{\rm fits}}^{{\rm diff}}(f,f^{\prime})$. In all
plot we assume the $\ell$-dependence of $C^{\rm diff}_{\ell}(f,f^{\prime})$
are the same. The whiteout regions are corresponding to the flagged frequency
channels from the data.

To pull $C^{\rm diff}_{\ell}(f,f^{\prime})$ out of the summation sign in Eq.
(32), we assume that for a given redundant block $\alpha$ the sky is slowly
changing function of $\ell\approx 2\pi|\textbf{{u}}|$. In other words, this
means that the angular power spectrum of the sky does not evolve significantly
over the baseline length within the redundant group.

Since sky power spectrum is real we discard the complex part of $C^{{\rm
diff}}_{\ell}(f,f^{\prime})$ in Eq. (34). Then, we implement the power-law
curve fitting into Eq. (34) to model the angular cross-frequency power
spectrum statistically. We make explicit assumptions about the functional form
of the fitting function based on the fact that the power spectrum of the
diffuse emission follows a simple power-law with frequency (Santos et al.,
2005). We implement the Scipy optimization routine for curve fitting to fit
data with the pre-specified function. With the fitting parameters $A$
(amplitude in ${\rm mK^{2}}$) and the spectral index $\alpha$, the fitted
power spectrum takes

$C^{\rm diff}_{{\rm
fits}}(f,f^{\prime})=A\left(\frac{ff^{\prime}}{f^{2}_{0}}\right)^{-\alpha},$
(35)

where the fitted parameters $A=249.69\,{\rm mK^{2}}$ and $\alpha=1.5$, and the
$f_{0}=150\,{\rm MHz}$ is the reference frequency. Since the power spectrum of
the diffuse emission does not evolve significantly over baselines in the
quasi-redundant group, we assume that the fitted power spectrum is not a
function of $\ell$.

Then the fitted angular power spectrum $C^{{\rm diff}}_{{\rm
fits}}(f,f^{\prime})$ (see Fig. 4) for that particular quasi-redundant group
will be inserted in place of $C^{{\rm diff}}_{\ell}(f,f^{\prime})$ in Eq. (32)
to estimate the model covariance matrix of the group. Likewise, one needs to
derive the best-fit cross-frequency power spectrum statistically from the
respective quasi-redundant sets to estimate the diffuse sky covariance matrix
in each redundant group. In terms of the fitted $C^{{\rm diff}}_{{\rm
fits}}(f,f^{\prime})$ parameter, the diffuse sky covariance matrix R finally
takes

$\displaystyle\textbf{{R}}=C^{\rm diff}_{{\rm
fits}}(f,f^{\prime})\sum_{\ell}b^{{\rm diff}}_{\ell}({\textbf{{u}}},f)b^{{\rm
diff},*}_{\ell}({\textbf{{u}}}^{\prime},f^{\prime}).$ (36)

We use Eq. (36) to represent the diffuse sky component for $\chi^{2}$
optimization process in Eq. (17). The dimension of the covariance matrix R
depends on the number of baseline $N_{\mathbf{\mathrm{u}}}$ and number of
frequency channels $N_{f}$, i.e. its dimension is equal to
$N_{\mathbf{\mathrm{u}}}N_{f}\times N_{\mathbf{\mathrm{u}}}N_{f}$.

####  3.1.2 Point Source Component Matrix (S)

The partial information about the point sources on the sky is one of the
components the CorrCal calibration scheme employs to estimates the best
likelihood antenna gain parameters. Simulation in Sievers (2017) showed that
the inclusion of the bright sources information in CorrCal has a significant
effect in reducing both amplitude and phase calibration error. However, the
phase error has reduced remarkably (about a factor of 5 better than non-source
informed correlation calibration) when the correlation information of point
sources has implemented in Sievers (2017), which substantially improves the
quality of phase calibration.

In this study, we use the published astronomical catalogues in Jacobs et al.
(2013) to model the known sources. Using sources information in this
catalogue, we form the matrix S that carries the sources information.

The source position is implicit in time and usually expressed in terms of
right ascension (RA) and declination (Dec) of catalogue source position. The
time dependence of sources position vector $\hat{\textbf{{r}}}_{i}$ thus comes
from source motion relative to the observer meridian, which best visualised
using the Hour-Angle (HA) and the observer Local Sidereal Time (LST). Because
HA defines the amount of time since the source transited the observer
meridian, hence, it tells how distant a source is from the meridian. Based on
the HA values (HA=LST-RA) of each source from the catalogue, we use sources
that are above the horizon during observations time to create S. We then
project each source position onto the sphere. Then, we predict the
visibilities in the spherical sky from a set of point sources with known
positions. Using the source statistics in this catalogue, we calculate the
predicted visibility with the help of Eq. (26) for known sources. That is,
suppose the point source profile is a Dirac-delta function on the sky
$\hat{\textbf{{r}}}_{i}$ ($i=1,2,...,N_{\rm ps}$, $N_{\rm ps}$ is the total
number of point sources), then the specific intensity (Eq. (25)) can be
written as

$\displaystyle I(\hat{\textbf{{r}}},f)=\sum^{N_{\rm
ps}}_{i=1}I(\hat{\textbf{{r}}}_{i},f)\delta_{\rm
2D}\left(\hat{\textbf{{r}}},\hat{\textbf{{r}}}_{i}\right),$ (37)

where $\delta_{\rm 2D}\left(\hat{\textbf{{r}}},\hat{\textbf{{r}}}_{i}\right)$
is a two-dimensional Dirac delta function. Substituting this function into Eq.
(26), we obtain the visibilities vector for point sources

$\displaystyle\textbf{{v}}^{\rm ps}(\textbf{{u}},f)=\sum^{N_{\rm
ps}}_{i=1}I(\hat{\textbf{{r}}}_{i},f)\left|A_{pq}(\hat{\textbf{{r}}}_{i},f)\right|^{2}{\rm
e}^{-{\rm i}2\pi\textbf{{u}}\cdot\hat{\textbf{{r}}}_{i}}.$ (38)

To simulate beam function for sources, we use the known position on the sky
from the established catalogue in Jacobs et al. (2013) and the PAPER-64
simulated beam model that displayed in Fig. 2. Taking the outer product
betweenthe visibility vector $\textbf{{v}}^{\rm ps}(\textbf{{u}},f)$ and its
complex conjugate transpose $\textbf{{v}}^{\rm ps,\dagger}(\textbf{{u}},f)$,
we form the covariance information for the discrete sources,

$\displaystyle\textbf{{S}}=\textbf{{v}}^{\rm
ps}(\textbf{{u}},f)\textbf{{v}}^{\rm ps,\dagger}(\textbf{{u}},f).$ (39)

Thus, we will use this source information covariance matrix to calculate the
$\chi^{2}$ function in Eq. (20).

####  3.1.3 Noise Matrix (N)

In this work, we estimate the per-diagonal noise variance matrix N from the
radiometer equation

$\displaystyle\sigma_{\rm rms}=\frac{T_{\rm sys}}{\sqrt{2\Delta f\tau}},$ (40)

where $\Delta f=97\,{\rm KHz}$ is the frequency band width in KHz, $\tau=49.2$
s is the integration time. The system temperature $T_{\mathrm{sys}}$ from the
relation defined in eq. (23) of Cheng et al. (2018) is,

$T_{\rm sys}=180\,{\rm K}\left(\frac{f}{0.18\,{\rm GHz}}\right)^{-2.55}+T_{\rm
rcvr}$ (41)

where $f$ is frequency, and $T_{\rm rcvr}$ is a receiver temperature which
equal to 144 K.

Figure 5: Noise variance estimated from the system temperature of PAPER.

###  3.2 Eigendecomposition of sky signal covariance matrix

The amplitude of R in Eq. (36) for a perfectly redundant group containing two
baselines is equal to a parameter controlling the variance of the sky,
$C_{\ell}$, times a square matrix of ones. That is,

$\displaystyle\textbf{{R}}=C_{\ell}^{2}\begin{bmatrix}1&1\\\ 1&1\\\
\end{bmatrix}\qquad,$ (42)

However, for a quasi-redundant group, all the off-diagonal elements of the
correlation matrix R are slightly different from one. As shown in Fig. 3, all
the elements of the covariance matrix are not precisely equal. It shows that
the correlation between baselines decreases as they become less redundant. In
that case, we approximate R from its eigendecomposition in our analysis
following Sievers (2017).

Figure 6: Left–Eigenvalues of cross-frequency covariance matrix described in
Eq. (36). The vertical axis represents the $\log_{10}$ eigenvalues, while the
horizontal axis labels the corresponding principal component index. Right–The
first three eigenmodes for R corresponding to the first three eigenvalues for
a frequency range from 120 to 168 MHz.

If each redundant group in the array is not represented by few eigenmodes that
are much less than the visibilities in the group, it shows that there is no
sufficient redundancy in the array and therefore any calibration scheme that
relies on the array redundancy may not give the accurate results. In the
eigen-decomposition analysis, few modes of the principal components accounted
for most of the variances in the original data. Thus, those eigenmodes of R
corresponding to sufficiently large eigenvalues could replace its original
data with very little loss of information. In the left panel of Fig. 6, we
show the eigenvalues of cross-frequency covariance matrix described in Eq.
(36). The vertical axis represents the $\log_{10}$ eigenvalues, while the
horizontal axis labels the corresponding principal component numbers. The
sharp fall-off in the eigenvalues against the eigenmode index suggests that,
by measuring a few modes of an eigenvalue one can account for most of the
variation in the original data that forms R. In addition, all the eigenvalues
are greater than zero is because R is a positive-definite covariance matrix.
In the right panel of Fig. 6, we depict the first three eigenmodes for R
corresponding to the first three eigenvalues for a frequency range from 120 to
168 MHz. Specially, the first and the second eigenmodes are more important to
explain the diffuse sky covariance between frequency channel because they are
a relatively smooth function of frequency.

We keep eigenmodes with amplitude more than $10^{-5}$ times the largest. This
allows the keeping of $\sim$2-3 complex eigenmodes per quasi-redundant block.
These few modes could represent most of the variations in a given redundant
group $\alpha$. For instance, the proportion of variance related to the
eigenvalues $\lambda_{1}$ and $\lambda_{2}$ can be calculated from

$\displaystyle\left(\frac{\lambda_{1}+\lambda_{2}}{\lambda_{1}+\lambda_{2}+\cdots+\lambda_{p}}\right)\times
100,$ (43)

where $p$ is the total eigenmodes. If the value of Eq. (43) is greater than
the threshold percentage $99\%$, we use eigenmodes corresponding to
eigenvalues $\lambda_{1}$ and $\lambda_{2}$ to create a matrix that carries
information about weights contributed by the variables in the original data to
these eigenmodes. Variables with the highest correlation (in absolute value)
with a principal component receive the highest weight in that component. This
vector can be calculated from the statistical relation

$\displaystyle\textbf{{p}}=\mathbf{e}_{i}\sqrt{\lambda_{i}},$ (44)

where $(\mathbf{e}_{i},\lambda_{i})$ is the eigenvector-eigenvalue pairs for
R, and the subscript $i$ is the principal component index. Using the vector
outer product of p thus we estimate a matrix that represents the diffuse sky
covariance, given by

$\displaystyle\textbf{{R}}_{\rm
PCA}=\textbf{{p}}\otimes\textbf{{p}}^{\dagger},$ (45)

where $\textbf{{R}}_{\rm PCA}$ is the diffuse sky covariance matrix calculated
from the principal component analysis (PCA) of R for a given redundant block
$\alpha$. The dimension of $\textbf{{R}}_{\rm PCA}$ would significantly be
reduced if very few eigenmodes could represent most of the variation in the
original data. For instance, if we keep only two eigenmodes of R for a given
redundant group containing four baselines, a complex vector p could be

$\displaystyle\textbf{{p}}=[p_{1},p_{2},0,0],$ (46)

where $p_{1}$ and $p_{2}$ are complex elements of a vector p that retained
from eigenmode analysis. The outer product of this vector by its complex
conjugate transpose takes

$\displaystyle\textbf{{p}}\otimes\textbf{{p}}^{\dagger}=\begin{bmatrix}p_{1}p^{*}_{1}&p_{1}p^{*}_{2}&0&0\\\
p_{2}p^{*}_{1}&p_{2}p^{*}_{2}&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{bmatrix},\quad$ (47)

here, the operator ∗ denotes a complex conjugate. Adding the noise to the
above matrix, we have

$\displaystyle\Gamma=\textbf{{N}}+\textbf{{p}}\otimes\textbf{{p}}^{\dagger}.$
(48)

All off-diagonal elements of the noise variance matrix N, which has the same
size as $\textbf{{p}}\otimes\textbf{{p}}^{\dagger}$, are zero. Thus, $\Gamma$
has only a relatively small number of nonzero elements and it can be
considered to be sparse because most of its elements are zero. Therefore, it
is wasteful to use the general dense matrix algebraic methods to inverting
$\Gamma$ since most of the $O(N^{3})$ operations devoted to inverting the
matrix involve zero operands. Hence, in terms of memory storage and
computational efficiency, performing the sparse matrices operation that
traversing only nonzero elements have significant advantages over their dense
matrix counterparts. Taking advantage of the sparse structure of the
correlation matrix $\Gamma$, CorrCal utilize the sparse matrix operations to
increases its computational performance in a large dataset.

To implement the sparse matrix computations in CorrCal, per-redundant group
correlation matrices $\Gamma$ are stored as a sub-matrices on the diagonal of
the square block-diagonal matrix. These sub-matrices will be square, and their
dimensions vary depending on the number of baselines in the redundant group.
Since the off-diagonal elements of the block-diagonal matrix are zero, and
there can be many sparse sub-matrices stored as its diagonal elements, the
overall matrix must be sparse. It can be easily compressed and thus requires
significantly less storage which also enhances the computational speed.

###  3.3 Implementation of CorrCal

The crucial steps to perform CorrCal are grouping data according to their
corresponding redundant baseline and determining the matrices R, S, and N to
each of these redundant groups. Then, one needs to set up the sparse matrix
representations for these matrices, determine the $\chi^{2}$ function and its
gradient information using the CorrCal routines. Utilizing these ingredients
as inputs, CorrCal finally adopts the
Scipy666https://github.com/scipy/scipy.gitbuilt-in optimization algorithm to
fits antenna gains. We summarise the procedures of implementation of CorrCal
as follows:

  * •

First, we need to create the sparse matrices representing R, S, and N set up
to each redundant group using the CorrCal routine, sparse_2level(). The
following arguments must be supplied to this routine to create sparse
matrices:

    1. 1.

The real-imaginary separated per-redundant group diffuse sky vector p. In this
vector, we make all its entries equal to zero except those entries
corresponding to eigenmodes that we kept per redundant block (see Eq. (46)).

    2. 2.

The real-imaginary separated vector s that carries information about the known
point sources.

    3. 3.

Per-diagonal noise variance matrix N.

    4. 4.

A vector that contains indices setting off the redundant group.

  * •

Second, the $\chi^{2}$ function and its gradient information should be
determined using the CorrCal subroutines. The arguments passed to these
subroutines in order to calculate the $\chi^{2}$ function and the gradient
information of $\chi^{2}$ are:

    1. 1.

A vector containing the real-imaginary separated initial guess for antenna
gains.

    2. 2.

Per-redundant group real-imaginary separated measured visibility data vector
with imaginaries follow their respective reals (e.g., see vectorised
visibility in Eq. (3.3.1)).

    3. 3.

The sparse matrix representations for R, S and N from the first step.

    4. 4.

A vector containing per-visibility antenna indices.

  * •

Finally, using the $\chi^{2}$, its gradient information, and a vector of the
initial guess for antenna gains that scaled by a large number as an input,
CorrCal adopts the Scipy built-in nonlinear conjugate gradient algorithm,
Scipy.optimize.fmin_cg(), to solve for antenna gains.

####  3.3.1 Antenna gain solutions as a function of frequency-bandpass
solution

For the EoR studies, highly precise frequency-dependent bandpass calibration
plays a crucial role. If one’s calibration solution has a noise-like structure
along the frequency axis, it may introduce fine-frequency structure to the
otherwise smooth sky observations even if the instrument’s response is the
slow function of observing frequencies. It poses the serious challenges toward
efforts to isolate spectrally structured weak 21-cm signal from spectrally-
smooth strong foreground signal. Several strategies have proposed for ensuring
smooth calibration solutions as a function of frequency. For instance,
simulation in Ewall-Wice et al. (2017) and Orosz et al. (2018) has shown that
down-weighting measurements from longer baselines in a redundant baseline
calibration solution would eliminate the migration of fine spectral structure
from longer baselines to the shorter baselines regimes. This effect is because
the longer baselines have more chromatic nature than the shorter ones as shown
in Morales et al. (2012). Other strategies implement fitting the calibration
solutions with smooth function in frequency for a bandpass calibration.

However, the natural formalism of CorrCal would allow one to fit for gain
solutions as a function of all frequency channels at once. It means that,
instead of solving an independent set of equations for every single frequency,
a visibility data d from different frequency channels put together in a longer
vector as,

$\displaystyle\textbf{{d}}=\begin{bmatrix}\operatorname{\mathbb{R}e}\\{\mathrm{v({f_{1})}}\\},\
\operatorname{\mathbb{I}m}\\{\mathrm{v(f_{1})}\\},\
\operatorname{\mathbb{R}e}\\{\mathrm{v(f_{2}})\\},\
\operatorname{\mathbb{I}m}\\{\mathrm{v(f_{2}})\\},\ \cdots\end{bmatrix},$ (49)

where $\operatorname{\mathbb{R}e}\\{\mathrm{v}(f_{1})\\}$,
$\operatorname{\mathbb{I}m}\\{\mathrm{v}(f_{1})\\}$,
$\operatorname{\mathbb{R}e}\\{\mathrm{v(f_{2}})\\}$, and
$\operatorname{\mathbb{I}m}\\{\mathrm{v(f_{2})}\\}$ are the real,
$\operatorname{\mathbb{R}e}$, and imaginary, $\operatorname{\mathbb{I}m}$,
parts of the visibility measurements corresponding to frequency channel
$f_{1}$ and $f_{2}$, respectively. Putting together visibility from all
frequency channels in a longer vector like that enforce CorrCal to execute
covariance information between frequencies. This imposes one to solve for
smooth bandpass gain solutions for all frequencies as a single calibration
step as depicted in Fig. 7. To some extent, the frequency structure of the
gain solutions in Fig. 7 is comparable to the autocorrelation result in Li et
al. (2019) and Byrne et al. (2019). Once these bandpass calibration solutions
are found, we divide these solutions to the raw data for CorrCal calibration.

Figure 7: The bandpass gain solution derived from a snapshot of 10minutes
observation (JD2456242.25733) of PAPER-64 elements. Amplitude (left) and phase
(right) of the gain solutions against frequency. Each line is a different
antenna. Discontinuity in the gain solution shows the flagged frequency
channel that affected by RFI contamination.

The correlation calibration suffers from the same degeneracy problem as with
redundant baseline calibration that presented in section 2.2. Therefore, the
gain solutions in CorrCal are degenerate, and these degeneracy issues will
have to be dealt with. In this work, the overall amplitude degeneracy in
CorrCal gain solutions is constrained by setting the average absolute value of
the gains equal to one. Because we are applying CorrCal to data that has
already been absolutely calibrated, we do not want to rescale the overall
amplitude. To set the overall phase degenerate parameter $\Delta$, we define
the reference antenna. We then set the overall phase by bringing the argument
of the gain of the reference antenna equal to zero. The phase gradient
degeneracy issues are resolved using the known point sources information on
the sky. Then, we solve for phase gradient parameters by fitting them with the
phase solutions of gain that have been obtained using sources information.
That means,

$\displaystyle\chi^{2}_{\rm phase}=\sum_{p=1}^{N}\left({\rm
Arg}[\hat{g}_{p}(f)]-\Delta_{x}b_{x,p}-\Delta_{y}b_{y,p}\right),$ (50)

where $b_{x,p}$ and $b_{y,p}$ are the $x$ and $y$ components of the baseline
$b$ for antenna $p$. Then, we subtract the fitted gradient parameters off from
the phase calibration solutions across the array to set for phase gradient
issues.

###  3.4 Comparisons of CorrCal to other methods

The CorrCal approach has a distinct feature to alleviate issues like
imperfections of array redundancy and miss-modeled sky catalog that redundant
baseline and sky-based calibration techniques suffer from. In CorrCal, sources
with known positions can be included as a prior even if their fluxes are not
precisely known, because the position of the sources is better known than
their intensities during observation time and frequency. This concept allows
CorrCal to assume the source catalog is complete and well-modeled with
precisely known source positions. The method works well if the sky-modeling
error could be mostly caused by uncertain/missed source fluxes in the sky-
based calibration approach.

The covariance-based formulation of $\chi^{2}$ in CorrCal provides significant
flexibility by relaxing the assumption of explicit redundancy to accounts for
imperfections of array redundancy. For instance, if baselines in a quasi-
redundant group are assumed to be alike but they are not perfectly alike, then
the CorrCal method incorporates this offset from perfection in its formalism
by suppressing the off-diagonal elements of covariance matrix within the
group. However, the redundant baseline method assumes that the baselines are
perfectly redundant in the group, and therefore all the off-diagonal elements
of the covariance matrix in the group are exactly equal(e.g. Eq. (42)).

##  4 The PAPER-64 Data Set

PAPER is a low-frequency radio interferometer experiment dedicated to probing
the EoR through the measurements of the 21-cm power spectrum (Parsons et al.,
2010). It was the first 21-cm array experiment with a redundant array
configuration, where antennas were arranged in a regular grid to generate
multiple copies of the same baselines. The redundant baselines probe the same
modes on the sky, resulting in significant improvement in the sensitivity of
the 21-cm power spectrum, thus, it is reasonable to use redundant calibration.
The array was operating at the South African SKA site in the Karoo desert
until its decommissioning in 2015. It was first deployed with 16 antennas and
continued to grow, eventually increasing to 128 antennas in 2015. PAPER has
1024 frequency channels across the 100 MHz band. It observes between 100-200
MHz, corresponding to the 21 cm signal from redshifts 6-13, with a frequency
resolution of $\sim$97 KHz. Data from the 32-elements PAPER array has been
analyzed in Parsons et al. (2014), which results in one of the very first
upper limits on the 21-cm power spectrum at $z=7.7$. Analyses of the data set
from the 64-elements PAPER array (hereafter PAPER-64) in (Cheng et al., 2018;
Kolopanis et al., 2019) have placed significant upper limits on the power
spectrum amplitude of 21 cm signal.

We will perform a recalibration of the PAPER-64 data set taken from Ali et al.
(2015) with CorrCal. This data set has been calibrated with redundant
calibration using the OmniCal software, and then absolute calibrated. The
absolute calibration has been made by fitting to the bright sources such as
Fornax A, Pictor A, and the Crab Nebula to set the overall phase. Using the
Pictor A calibrator, the overall amplitude in this data has been fixed. The
detailed discussions about the absolute calibration stage for this data are
presented in Ali et al. (2015). The data set has dual polarization products.
We use one night of observations taken between 10 November 2012
(JD2456242.17382) and 11 November 2012 (JD 2456242.65402). As depicted in Fig.
8, the PAPER-64 array system comprises 64 dipole antennas arranged on a
regular grid spacing with the total number of baselines equal to
$64\times(64-1)/2=2016$ and thus measure 2016 visibilities. For our
calibration, we use all baselines except those associated with three antennas
which were flagged during the observation due to known spectral instability.
The observation parameters of this data set are summarized in Table 1. We
direct readers to Ali et al. (2015) for more details of observation
strategies, flagging of radio frequency interference (RFI), channelization of
data, antenna cross-talk minimization, and other technical information of the
PAPER array. Our method proposes a new possibility to fully gauge the
systematics and noise of radio interferometry.

Figure 8: Antenna position (left) and UV coordinate of the PAPER-64 array. The North-South and East-West antenna coordinates are in meter, and the UV coordinates are in wavelengths at 150 MHz (right). Table 1: Summary of observational parameters and their values used in this work Parameter |  Value
---|---
PAPER array location |  $\mathrm{30.7^{o}S,21.4^{o}\,E}$
Observation dates |  10 Nov 2012-11 Nov 2012
Observing mode |  drift-scan
Time resolution |  $\sim$42.9 seconds
Frequency range |  120 MHz-168 MHz
Frequency resolution |  0.49 MHz
Number of antenna |  61
Field-of-view |  60∘
Visibility polarization |  $XX$,$YY$
Shortest baseline |  4m
Longest baseline |  212m

##  5 Results

###  5.1 Post-CorrCal spectral structure of data

After CorrCal implementation, we have tried to compare the structure of data
as a function of frequency. Figure 9 depicts the real part of visibility
against frequency for a single integration time of a given baseline. The
careful examination of the plot shows that the real-visibility gets spectrally
smoother (red curve) after CorrCal re-calibration. It signifies the reduction
of the spectral structure of the visibility measurements after re-calibration.
Considering that, in the following section, we present the effect of re-
calibration on the delay transformed power spectra (Parsons & Backer, 2009).

Figure 9: The real part of visibility data (JD2456242.25733) as a function of
frequency at a single timestamp for an EW baseline of length about $30$ m for
both calibrations approaches. Line discontinuity corresponds to the flagged
frequency channels. The spectral structure looks slightly reduced after
CorrCal (red curve).

###  5.2 Delay-space power spectra

The statistical analysis of the Fourier mode of the power spectrum
($k_{\perp},k_{\parallel}$) for low-frequency observations is a powerful tool
to probe the EoR (Morales & Hewitt, 2004). This method shows the existence of
distinctive regions that dominated by 21-cm signal and foreground contaminants
in a 2D line-of-sight ($k_{\parallel}$) and transverse comoving
$\left(k_{\perp}=\sqrt{k^{2}_{x}+k^{2}_{y}}\right)$ Fourier plane. The
spectrally smooth foreground is expected to contaminate the region of low
$k_{\parallel}$ values, while the vast majority of $k_{\parallel}$ values are
free from these contaminants (Liu et al., 2014a). Nevertheless, several
studies have shown that foreground contamination may go further to the higher
$k_{\parallel}$ values because of chromatic interaction of interferometry with
intrinsically smooth foreground (Parsons et al., 2012; Morales et al., 2012;
Dillon et al., 2014; Vedantham et al., 2012; Thyagarajan et al., 2013; Liu et
al., 2014a, b). This effect leads to the distinctive wedge-like structure in
$k_{\perp}$ and $k_{\parallel}$ Fourier space, leaving the foreground-free EoR
window beyond the wedge (see results in Fig. 11 and discussions there).

Figure 10: Delay power spectra on either side of $k_{\parallel}$ direction.
Foreground wedges are confined under the region bounded by the black line or
the horizon line. (A) the delay power spectra for data already calibrated with
OmniCal, (B) power spectra after CorrCal recalibration, and the lower panel
displays linear scale percentage deviation. The negative percentage deviation
in power at the bin right on the horizon shows that the wedge has become
tighter after CorrCal recalibration.

To convert the measured visibility to 2D power spectrum in $k_{\parallel}$ and
$k_{\perp}$ space, the per-baseline delay transform method is used. The
technique was first introduced by Parsons & Backer (2009). As thoroughly
presented in Parsons et al. (2012), the delay transformed visibility is
obtained by inverse Fourier transforming of each baseline along the frequency
axis. The method allows the smooth foreground to localize at the central
region between the negative and positive geometric delay limits, and any
spectrally unsmooth structure such as 21-cm signal may spread across delay
axis (Parsons et al., 2012). We will review the mathematical formalism of the
per baseline delay-transform method in Appendix A. In our analysis, we use
Stokes $I$ visibilities777Note that these are not true Stokes $I$ would rather
the average of $XX$ and $YY$ visibilities.(see Sec. A) taken by PAPER-64 array
from Julian date 2456242.17382 to 2456242.65402 (2012 November 10-11) to form
the 2D power spectrum using delay transform method, a total of about 9 hours
observation over frequencies 140 to 160MHz.

Following Parsons et al. (2012) and Pober et al. (2013), we apply the
Blackman-Harris window function during the delay-transform stage to minimize
the foreground leakage caused by the finite bandwidth of frequency. This
technique allows the EoR window to be mostly dominated by 21-cm brightness
temperature by reducing the foreground spectral leakage to the smallest
possible amount (Parsons et al., 2012; Pober et al., 2013). We employ the 1D
CLEAN algorithm on the delay transformed data to reduce the effect introduced
by RFI flagging and to increases the resolution in delay space (Pober et al.,
2013). We then cross-multiply the consecutive integration time to avoid the
noise bias. After these stages, we form the cylindrical delay power spectrum
in $k$ space using Eq. (55).

Figure 10 shows the delay power spectra where the positive and negative values
of $k_{\parallel}$ have not been averaged together. This form of power spectra
is expected because emissions from the negative and positive delay directions
are proportional to the negative and positive $k_{\parallel}$ respectively.
The delay power spectra from emissions near-zero delays (where
$k_{\parallel}\sim 0$) are confined within the central region of the 2D delay
power spectra. It corresponds to foreground emissions from the peak primary
beam response of the instrument.

However, from Fig. 10, we see that the strength of power spectra around zero
delay decreases as one moves from the lower $k_{\perp}$ to higher $k_{\perp}$
regions. Furthermore, for both calibration approaches, more emission goes
beyond the horizon limit for shortest baselines (smallest $k_{\perp}$ values).
This effect is because shorter baselines resolve out less of the Galactic
synchrotron so that emission will be brighter and sidelobes extend further in
$k_{\parallel}$ (Pober et al., 2013).

The black diagonal line represents the horizon limit determined by baseline
length and the source location relative to the main field of view of the
array, marking the boundaries of the foreground wedge. Its analytic expression
has been defined in Eq. (57). The spectral-smooth emissions are supposed to be
confined within that limit, and any emissions which are intrinsically unsmooth
with frequency are moving beyond the horizon limit. The calibration error
plays a non-negligible role in this regard in imparting the spectral structure
to emissions that were originally spectrally smooth as discussed earlier. This
unnatural spectral structure, imparted due to calibration error or other
effects, will affect the EoR window by scattering power beyond the wedge
limits. To see CorrCal’s impact on the delay power spectra, we compare
foreground power in a wedge region of the power spectrum, and outside of it in
which the 21-cm signal is dominating by taking the percentage deviation
between power after CorrCal and OmniCal, shown in the bottom panel of Fig. 10
. It is apparent that the power spectra inside and outside the wedge region
are comparable for both calibration approaches. However, we observe a
consistent drop in power at the bin right on the horizon after CorrCal. Most
importantly, this is consistent with some of the analyses done with MWA in (Li
et al., 2018, 2019; Zhang et al., 2020) that demonstrate redundant
calibration’s ability to make small improvements in the modes nearest (but not
inside) the wedge. This could be more suggestive of a reduction in the
spectral structure of the data after re-calibration with CorrCal.

Figure 11: Same as Fig. 10 but averaged and folded over the delay centre to
obtain the “wedge”-shaped delay power spectra. From left to right: (A) The
power spectra generated from data already calibrated with OmniCal package; (B)
the power spectra obtained after data have recalibrated using without source
information in CorrCal; (C) power spectra after recalibration using source
information in CorrCal. The right two panels show the linear scale percentage
deviation without and with sources information in CorrCal formalism
respectively. The black diagonal line marks the horizon limit ($\tau_{\rm
H}$).

The well-known wedge-shaped foreground power spectra from our analysis are
shown in Fig. 11. We generate this form of power spectra by averaging power
from both positive and negative $k_{\parallel}$ directions and folding it over
the centre of the zero-delay line (where $k_{\parallel}=0$) of Fig. 10. A
similar feature of power spectra has been observed from PAPER-64 data in Pober
et al. (2013), as displayed in Figure 3 of their work. Moreover, Kohn et al.
(2016) have reported a wedge-like structure from the analysis of PAPER-32
data. In our analysis, Fig. 11(A) and (B), respectively, displaying the wedge-
shaped power spectra to data already calibrated using OmniCal and to data
after CorrCal recalibration. The middle panel in Fig. 11 (C) shows power
spectra after data have been calibrated using point sources information in
CorrCal. For each case, the wedge extends to higher values of $k_{\parallel}$
at a longer baseline ($k_{\perp}$) regime, demonstrating how the theoretical
wedge limit increases as a function of baseline length. This effect is also
due to sources located far away from the pointing centre of the primary beam
of the instrument appears at higher delays, which corresponds to higher
$k_{\parallel}$. The relative deviation of power in percent between CorrCal
and OmniCal calibration approaches are shown in the third (CorrCal without sky
sources information ) and fourth (CorrCal with sky sources information) panels
of Fig. 11 from left to right , respectively. The negative percentage
deviation at the bin right on the horizon in $\mathrm{[(B-A)/A]\times 100\%}$
and $\mathrm{[(C-A)/A]\times 100\%}$ of Fig. 11 indicating that the power
spectra around the wedge limit have become slightly tighter after CorrCal
recalibration. This effect is not considerably saturated if we include the
sources’ information in CorrCal formalism as shown in the right panel of the
figure.

Figure 12: Right–Two-dimensional power spectra after (A) OmniCal, (B) CorrCal
recalibration without sources info, (C) CorrCal recalibration with sources
information, the lower panels are the percentage deviation of B to A and C to
A, respectively. The black dashed vertical line is representing the horizon
($\tau_{\rm H}$) limit. Left–The averaged slice of power spectra from
$k_{\perp}=0.0733~{}h\,{\rm Mpc}^{-1}$ to $k_{\perp}=0.0736~{}h\,{\rm
Mpc}^{-1}$ to both calibration schemes. The top panel is the one-dimensional
average power spectrum as a function of $k_{\parallel}$ for a portion of
baselines ($\sim 150\rm m$) to both calibration strategies. The bottom panel
is the percentage deviation.

The R.H.S panel of Fig. 12 depicts a portion of 2D power spectra for a given
group of baselines whose lengths are about $150$ m. We average over
$k_{\perp}$ axis over the range of $0.0733\,h\,{\rm
Mpc}^{-1}<k_{\perp}<0.0736\,h\,{\rm Mpc}^{-1}$ to generate the corresponding
1D power spectra that shown on the L.H.S of Fig. 12. From this figure, we
observe that there is a steep fall-off of the power spectrum right away the
horizon limit as expected, showing that the foreground power confined within
the slice of the wedge is brighter than the 21-cm signal by 4-5 orders of
magnitude. This sharp fall-off follows the same trend for both calibration
schemes. However, the inclusion of sky sources statistics as a prior during
CorrCal recalibration step doesn’t show the significant change in the power
spectra, as shown on the L.H.S of Fig. 12.

##  6 Conclusions

In this work, we used a new calibration scheme (CorrCal) to recalibrate the
PAPER64 observations. Our new calibration scheme boils down to the covariance-
based calculation of $\chi^{2}$, allowing one to include the relaxed array
redundancy and sky information into its formalism. This formalism also
provides space for cross-frequency bandpass calibration. Furthermore, the
inclusion of known sources information in CorrCal breaks the statistical
isotropy of the sky, which sets the overall phase gradient degenerate
parameter across the array. We have compared the spectral smoothness of the
observed data before and after CorrCal’s application. We observed that the
spectral structure of the visibility data after CorrCal becomes slightly
reduced, for instance, as shown in Fig. 9. The significant factor for this
improvement could be the CorrCal natural formalism to carry out the cross-
frequency bandpass calibration.

We implemented the delay-transform method (Parsons & Backer, 2009; Parsons et
al., 2012) to filter out smooth foregrounds from the spectrally structured
weak 21-cm signal in 2D Fourier $k-$space. The comparison results in the delay
space show that the power spectra around the foreground wedge limit have been
reduced by about $6\%$ after the CorrCal re-calibration. Nevertheless, this
effect does not significantly improve the power spectrum in any bin around the
horizon limit of the foreground wedge and further investigation will be needed
to test the algorithm effectiveness and its implementation methods. The
current numerical tests lay the foundation of calibration radio instruments
using correlated signals on the sky. Our future work will include possible
beam shape errors, antenna pointing errors, and position errors, and improve
the CorrCal method by using sky and antenna simulation (e.g., pyuvsim) and
compare it with sky-based calibration and redundant baseline calibration. In
conclusion, the CorrCal may be an alternative calibration method for radio
interferometry, which potentially helps to probe the EoR using the 21cm
signal.

## Acknowledgements

We would like to acknowledge the PAPER collaboration for providing us data
that have been used in this work. Y.Z.M. acknowledges the support of
NRF-120385, NRF-120378, and NSFC-11828301. P.K. acknowledges funding support
from the South African–China Collaboration program in Astronomy at the
University of KwaZulu-Natal.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ali et al. (2015) Ali Z. S., et al., 2015, ApJ, 809, 61
  * Barkana & Loeb (2001) Barkana R., Loeb A., 2001, Phys. Rep., 349, 125
  * Barry et al. (2016) Barry N., Hazelton B., Sullivan I., Morales M. F., Pober J. C., 2016, MNRAS, 461, 3135
  * Beardsley et al. (2016) Beardsley A. P., et al., 2016, The Astrophysical Journal, 833, 102
  * Bernardi et al. (2009) Bernardi G., et al., 2009, A&A, 500, 965
  * Bowman et al. (2009) Bowman J. D., Morales M. F., Hewitt J. N., 2009, ApJ, 695, 183
  * Byrne et al. (2019) Byrne R., et al., 2019, ApJ, 875, 70
  * Chapman & Jelić (2019) Chapman E., Jelić V., 2019, arXiv e-prints, p. arXiv:1909.12369
  * Cheng et al. (2018) Cheng C., et al., 2018, ApJ, 868, 26
  * Datta et al. (2010) Datta A., Bowman J. D., Carilli C. L., 2010, The Astrophysical Journal, 724, 526
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Dillon et al. (2013) Dillon J. S., Liu A., Tegmark M., 2013, Phys. Rev. D, 87, 043005
  * Dillon et al. (2014) Dillon J. S., et al., 2014, Phys. Rev. D, 89, 023002
  * Dillon et al. (2018) Dillon J. S., et al., 2018, MNRAS, 477, 5670
  * Dillon et al. (2020) Dillon J. S., et al., 2020, arXiv e-prints, p. arXiv:2003.08399
  * Edler et al. (2021) Edler H. W., de Gasperin F., Rafferty D., 2021, A&A, 652, A37
  * Ewall-Wice et al. (2017) Ewall-Wice A., Dillon S. J., Liu A., Hewitt J., 2017, Monthly Notices of the Royal Astronomical Society, 470, 1849
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Hothi et al. (2021) Hothi I., et al., 2021, MNRAS, 500, 2264
  * Hurley-Walker et al. (2017) Hurley-Walker N., et al., 2017, MNRAS, 464, 1146
  * Jacobs et al. (2013) Jacobs D. C., et al., 2013, The Astrophysical Journal, 776, 108
  * Kohn et al. (2016) Kohn S. A., et al., 2016, The Astrophysical Journal, 823, 88
  * Kolopanis et al. (2019) Kolopanis M., et al., 2019, ApJ, 883, 133
  * Li et al. (2018) Li W., et al., 2018, The Astrophysical Journal, 863, 170
  * Li et al. (2019) Li W., et al., 2019, ApJ, 887, 141
  * Liu & Parsons (2016) Liu A., Parsons A. R., 2016, MNRAS, 457, 1864
  * Liu & Tegmark (2011) Liu A., Tegmark M., 2011, Phys. Rev. D, 83, 103006
  * Liu et al. (2009) Liu A., Tegmark M., Bowman J., Hewitt J., Zaldarriaga M., 2009, Monthly Notices of the Royal Astronomical Society, 398, 401
  * Liu et al. (2010) Liu A., Tegmark M., Scott M., Lutomirski A., Zaldarriaga M., 2010, Monthly Notices of the Royal Astronomical Society, 408, 1029
  * Liu et al. (2014a) Liu A., Parsons A. R., Trott C. M., 2014a, Phys. Rev. D, 90, 023018
  * Liu et al. (2014b) Liu A., Parsons A. R., Trott C. M., 2014b, Phys. Rev. D, 90, 023019
  * Madau et al. (1996) Madau P., Ferguson H. C., Dickinson M. E., Giavalisco M., Steidel C. C., Fruchter A., 1996, MNRAS, 283, 1388
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Moore et al. (2013) Moore D. F., Aguirre J. E., Parsons A. R., Jacobs D. C., Pober J. C., 2013, ApJ, 769, 154
  * Morales & Hewitt (2004) Morales M. F., Hewitt J., 2004, The Astrophysical Journal, 615, 7
  * Morales & Wyithe (2010) Morales M. F., Wyithe J. S. B., 2010, Annual Review of Astronomy and Astrophysics, 48, 127
  * Morales et al. (2012) Morales M. F., Hazelton B., Sullivan I., Beardsley A., 2012, The Astrophysical Journal, 752, 137
  * Oh (2001) Oh S. P., 2001, The Astrophysical Journal, 553, 499
  * Orosz et al. (2018) Orosz N., Dillon J. S., Ewall-Wice A., Parsons A. R., Thyagarajan N., 2018, arXiv e-prints, p. arXiv:1809.09728
  * Paciga et al. (2011) Paciga G., et al., 2011, Monthly Notices of the Royal Astronomical Society, 413, 1174
  * Parsons & Backer (2009) Parsons A. R., Backer D. C., 2009, The Astronomical Journal, 138, 219
  * Parsons et al. (2010) Parsons A. R., et al., 2010, AJ, 139, 1468
  * Parsons et al. (2012) Parsons A. R., Pober J. C., Aguirre J. E., Carilli C. L., Jacobs D. C., Moore D. F., 2012, The Astrophysical Journal, 756, 165
  * Parsons et al. (2014) Parsons A. R., et al., 2014, The Astrophysical Journal, 788, 106
  * Patil et al. (2017) Patil A. H., et al., 2017, ApJ, 838, 65
  * Pober et al. (2013) Pober J. C., et al., 2013, The Astrophysical Journal, 768, L36
  * Pober et al. (2014) Pober J. C., et al., 2014, The Astrophysical Journal, 782, 66
  * Pober et al. (2016) Pober J. C., et al., 2016, ApJ, 819, 8
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Rahimi et al. (2021) Rahimi M., et al., 2021, MNRAS, 508, 5954
  * Santos et al. (2005) Santos M. G., Cooray A., Knox L., 2005, The Astrophysical Journal, 625, 575
  * Shaver et al. (1999) Shaver P. A., Windhorst R. A., Madau P., de Bruyn A. G., 1999, A&A, 345, 380
  * Sievers (2017) Sievers J. L., 2017, arXiv e-prints, p. arXiv:1701.01860
  * Thyagarajan et al. (2013) Thyagarajan N., et al., 2013, The Astrophysical Journal, 776, 6
  * Tingay et al. (2013) Tingay S. J., et al., 2013, Publications of the Astronomical Society of Australia, 30, e007
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Vedantham et al. (2012) Vedantham H., Shankar N. U., Subrahmanyan R., 2012, The Astrophysical Journal, 745, 176
  * Wang et al. (2006) Wang X., Tegmark M., Santos M. G., Knox L., 2006, ApJ, 650, 529
  * Wayth et al. (2018) Wayth R. B., et al., 2018, Publ. Astron. Soc. Australia, 35, 33
  * Wieringa (1992) Wieringa M. H., 1992, Experimental Astronomy, 2, 203
  * Woodbury (1950) Woodbury M. A., 1950, Memorandum report, 42, 336
  * Zhang et al. (2020) Zhang Z., et al., 2020, Publ. Astron. Soc. Australia, 37, e045
  * Zheng et al. (2014) Zheng H., et al., 2014, Monthly Notices of the Royal Astronomical Society, 445, 1084
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, A&A, 556, A2

##  Appendix A Delay Spectrum

After data with $XX$ and $YY$ linear polarisation being recalibrated using
CorrCal separately, we form the Stokes $I$ parameter by combining each linear
polarisation following Moore et al. (2013):

$V_{\mathrm{I}}=\frac{1}{2}(V_{XX}+V_{YY}),$

where $V_{XX}$ and $V_{YY}$, respectively, are visibilities from $XX$ and $YY$
linear polarisation. We then use $V_{\mathrm{I}}$ (Stokes $I$ polarisation
visibility) to form the power spectra using a delay-transform approach
following that of Ali et al. (2015).

Approximating the interferometric visibility equation that is defined in Eq.
(25) over the flat-sky model, we have

$\displaystyle\textbf{{d}}(u,v,f)=\int\textrm{d}l\,\textrm{d}m\,I(l,m,f)A(l,m,f)e^{-2\pi
i(ul+vm)},$ (51)

where $(l,m)$ are direction cosines of a unit vector $\hat{\textbf{{r}}}$,
which pointing to the source on the sky. These sky-plane direction cosines
have the corresponding antenna-plane Fourier dual $u$ and $v$, respectively,
that is defined in Section 3.1. Following the representation of Parsons &
Backer (2009), we can rewrite Eq. (51) in terms of geometric time delay
$\tau_{\mathrm{g}}$ for a single baseline b as

$\displaystyle\textbf{{d}}_{\textbf{{b}}}(f)=\int\textrm{d}l\,\textrm{d}m\,I(l,m,f)A(l,m,f)e^{-2i\pi
f\tau_{\rm g}},$ (52)

where

$\displaystyle\tau_{\mathrm{g}}=\frac{1}{c}\textbf{{b}}\cdot\hat{\textbf{{r}}}=\frac{1}{c}(b_{x}l+b_{y}m),$
(53)

with $\textbf{{b}}=(b_{x},b_{y})$, $(u,v)=f\mathbf{b}/c$ and
$\hat{\textbf{{r}}}=(l,m)$. Here, $b_{x}$ and $b_{y}$ are projected baseline
length in meter along east and north direction in the antenna plane,
respectively. Note that $\tau_{\rm g}$ in Eq. (52) shows the time-shifting of
signal arriving one antenna relative to other due to source location on the
sky with respect to the baseline vector orientation in the antenna plane.

Delay transform technique begins by converting the frequency spectrum of the
visibility into the delay spectrum through application of Fourier transform on
the visibility along its frequency direction. Applying the windowing function
$W(f)$ over Eq. (52) and then Fourier transforming it along its frequency
axis, the delay transformed visibility takes a form

$\displaystyle\tilde{V}_{\mathrm{b}}(\tau)$ $\displaystyle=$
$\displaystyle\int_{B}\textrm{d}f\,W(f)v^{\mathrm{mea}}_{\rm b}(f)e^{2\pi
if\tau}$ $\displaystyle=$
$\displaystyle\int_{B}\textrm{d}f\,W(f)\left[\int\textrm{d}l\,\textrm{d}m\,I(l,m,f)A(l,m,f)e^{-2\pi
if\tau_{\rm g}}\right]$ $\displaystyle\times$ $\displaystyle e^{2\pi if\tau}$
$\displaystyle=$
$\displaystyle\int_{B}\textrm{d}l\,\textrm{d}m\,\tilde{W}(\tau)*\tilde{I}(l,m,\tau)*\tilde{A}(l,m,\tau)*\delta_{\rm
D}(\tau-\tau_{\rm g}),$

where the window function $W(f)$ chosen by the data analyst is used to control
the quality of delay spectrum as suggested in Vedantham et al. (2012) and
Thyagarajan et al. (2013), $B$ is the bandwidth of observation over which the
integration has to be carried out, the delay parameter $\tau$ is the Fourier
dual of $f$ in unit of time, $\tilde{V_{\mathrm{b}}}(\tau)$ is the delay
transformed visibility observed by a baseline $b$, * stands for convolution.
Hence, from Eq. (A) one can clearly sees that there is mapping between point
sources on the sky in $\tilde{V}_{\mathrm{b}}$ space and the Dirac-delta
function $\delta_{D}$, convolved by Fourier transform of window function
$\tilde{W}(\tau)$, sky flux density $\tilde{I}(l,m,\tau)$, and the antenna
beam response $\tilde{A}(l,m,\tau)$, as shown in Parsons et al. (2012).
Analytically, the delay formalism that relates the spatial power-spectrum
$P_{21}(k_{\parallel},k_{\perp})$ of 21-cm signal from the EoR fluctuations to
the delay-transformed visibility $\tilde{V}$ has been derived in Parsons et
al. (2012) as

$\displaystyle
P_{21}(k_{\parallel},k_{\perp})\approx\tilde{V}_{21}^{2}\left(\frac{\lambda^{2}}{2k_{\mathrm{B}}}\right)^{2}\frac{X^{2}Y}{\Omega
B},$ (55)

where $\lambda$ is the observing wavelength, $k_{\mathrm{B}}$ is the Boltzmann
constant, $B$ is the observing bandwidth, $\Omega$ is the solid angle of the
primary-beam of the antenna, $X$ and $Y$ are cosmological scalars which
convert observed angles and frequencies into $h\,\mathrm{Mpc^{-1}}$ as
calculated in Parsons et al. (2012).

Following the method in Pober et al. (2013), the successive time stamp in
delay-transformed visibilities has been cross-multiplied to avoid the noise-
bias. Using the formalism in Kohn et al. (2016), it can be written as

$\displaystyle\tilde{V}_{21}^{2}\approx|\tilde{V}_{21}(\tau,t)\times\tilde{V}_{21}(\tau,t+\Delta
t)e^{i\theta_{\rm zen}(\Delta t)}|^{2},$ (56)

here $\Delta t\approx$ 42.9 seconds (a time resolution of PAPER-64 data) and
$\theta_{\rm zen}(\Delta t)$ is the zenith re-phasing factor. The line-of-
sight $k_{\parallel}$, and transverse co-moving $k_{\perp}$ components of $k$
in terms of the instrumental and cosmological parameters have shown in Morales
& Hewitt (2004),

$k_{\perp}=\frac{2\pi|\mathbf{u}|}{D_{\mathrm{m}}(z)},$

$k_{\parallel}\approx\eta\frac{2\pi H_{0}f_{10}E(z)}{c(1+z)^{2}},$

where $|\mathbf{u}|=\sqrt{u^{2}+v^{2}}$, $\eta$ is the Fourier conjugate of
$f$ which used to denote the spatial frequency along the line-of-sight and it
has the unit of time, $f_{10}=1420$ MHz is the rest frequency of 21-cm
emission from cosmic reionization, $c$ is the speed of light in free space,
$D_{\mathrm{m}}(z)$ is the transverse comoving distance at the reference
redshift $z$, the present day Hubble constant $H_{0}=70\,{\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$, the function
$E(z)=\sqrt{\Omega_{\mathrm{m}}(1+z)^{3}+\Omega_{k}(1+z)^{2}+\Omega_{\Lambda}}$.
Here the cosmological density parameters $\Omega_{\mathrm{m}}$,
$\Omega_{\Lambda}$ and $\Omega_{k}$ are total matter density, the dark energy
density and curvature, respectively. The geometric delay $\tau_{\mathrm{g}}$
in Eq. (53) sets the maximum delay limit beyond which no emissions enter the
interferometry involving that baseline. This limit will be attained when the
angle between b and $\hat{\textbf{{r}}}$ is set to be zero. Such an alignment
to both vectors is achieved, the maximum delay is commonly referred to
"horizon limit ($\tau_{\rm H}$)". Theoretically, over the
$(k_{\perp},k_{\parallel})$ cylindrical space, it can be determined from

$\displaystyle k_{\parallel}=\frac{H_{0}E(z)D_{\rm m}(z)}{c(1+z)}k_{\perp}.$
(57)

As an exhaustively discussed in Parsons et al. (2012), any strictly smooth
foreground emissions with power-law spectra are confined within the narrow
range of horizon limits such that $-\tau_{\rm H}\leq\tau\leq\tau_{\rm H}$.
Non-smooth power either intrinsic to the sky or instrumental effects scattered
beyond the horizon limit as sidelobes of the convolving kernel of sky and
beam, $\tilde{I}*\tilde{A}$, which broadens the footprint of $k_{\parallel}$
modes in cosmological $k$ space. This kernel is narrow for smooth spectrum
foreground sources and confined within the wedge. However, flux beyond the
horizon limit is representing the power-spectrum of unsmooth emissions such as
21-cm signal. Therefore, regions in a wide range of $k_{\parallel}$ modes are
considerably affected by the 21-cm signal if one keeps the frequency response
of the instrument smoother. More of these modes are free from foreground
contamination on shortest baseline, i.e., at smallest $k_{\perp}$ regions
(Datta et al., 2010; Vedantham et al., 2012; Parsons et al., 2012; Pober et
al., 2013).

##  Appendix B Determining the EoR signal loss in CorrCal analysis

To determine whether our analysis methods are lossless or not in the EoR power
spectrum estimate, first, we simulate visibility data using HERA simulation
pipeline (hera_sim888https://github.com/HERA-Team/hera_sim.git). For this
simulation, we use the analytic airy beam model (the same for all antenna
elements), the simulated antenna gains, the GaLactic and Extragalactic All-sky
MWA (GLEAM) sky sources(Hurley-Walker et al., 2017), the simulated EoR signal
using hera_sim. Then, to quantify our new calibration approach does not lead
to the EoR power spectrum signal suppression, we propagate simulated data into
the CorrCal for calibration. We then compute the 2D power spectrum for both
simulated and recovered data after CorrCal calibration using the delay
spectrum approach described in Appendix A. Fig. 13 shows the delay space power
spectra for both simulated and recovered data. From the figure, it is apparent
that the recovered power spectra after CorrCal are comparable to the expected
one, confirming that CorrCal calibration method does not lead to significant
signal loss.

Figure 13:  Top–(A) Expected delay power spectra, (B) recovered power spectra
after CorrCal calibration , (C) the power spectra ratio of B to A. The solid
black dashed line is the horizon limit, $\tau_{\rm H}$. Bottom–Averaged 1D
power spectra as a function of $k_{\parallel}$ for a portion of baselines in
the $\sim 30~{}\rm m$ group.