The correlation calibration of PAPER-64 data

Gogo, Tamirat G; Ma, Yin-Zhe; Kittiwisit, Piyanat; Sievers, Jonathan L; Parsons, Aaron R; Pober, Jonathan C; Jacobs, Daniel C; Cheng, Carina; Kolopanis, Matthew; Liu, Adrian; Kohn, Saul A; Aguirre, James E; Ali, Zaki S; Bernardi, Gianni; Bradley, Richard F; DeBoer, David R; Dexter, Matthew R; Dillon, Joshua S; Klima, Pat; MacMahon, David H E; Moore, David F; Nunhokee, Chuneeta D; Walbrugh, William P; Walker, Andre

doi:10.1093/mnras/stab3516

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. 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, which can suffer from sky-modelling error and array redundancy imperfection 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 apply the CorrCal to the data of Precision Array for Probing the Epoch of Reionization (PAPER) experiment, which was pre-calibrated using redundant baseline calibration. We show about |$6\%$| suppression at the bin right on the horizon limit of the foreground wedge-like structure, relative to the pre-calibrated power spectra. 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.

instrumentation: interferometers, methods: observational, telescopes, cosmology: dark ages, reionization, first stars

1 INTRODUCTION

The cosmic dark ages 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. Barkana & Loeb 2001; Oh 2001; Furlanetto, Oh & Briggs 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, because 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 constrained 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 tightest 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, de Gasperin & Rafferty 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 H ii region (Shaver et al. 1999; Santos, Cooray & Knox 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; Bowman, Morales & Hewitt 2009; Liu et al. 2009; Liu & Tegmark 2011; Parsons et al. 2012; Dillon, Liu & Tegmark 2013; Liu, Parsons & Trott 2014a; Beardsley et al. 2016; Pober et al. 2016; Hothi et al. 2021). These techniques, in general, can be categorized 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 tries 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 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 (CorrCal¹), 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-h 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 has 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-h 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,

$$\begin{eqnarray*} d_{pq} = g_{p}\, g^{\ast }_{q}\, v_{pq}+\, n_{pq}. \end{eqnarray*}$$

(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 equation (1) are implicitly time-dependent. The antenna gain g can also be written as a complex exponential form with amplitude a and phase ϕ for each of the antenna p (or q),

$$\begin{eqnarray*} g_{p} = a_{p} \exp (\mathrm{i}\phi _{p}). \end{eqnarray*}$$

(2)

Note that we used ‘|$\mathrm{i}$|’ 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

$$\begin{eqnarray*} {\boldsymbol {d}}= {\bf \textsf {G}}{\bf \textsf {v}}+ {\boldsymbol {n}}. \end{eqnarray*}$$

(3)

Here, d is a vector containing the measured visibilities of different baselines. |${\bf \textsf {v}}$| is vector of the corresponding true sky visibility. |${\bf \textsf {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 |${\bf \textsf {G}}$| can be categorized into either one of the two conventional approaches, i.e. sky-based calibration or redundant-baseline calibration.² 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 ≈ v_pq to solve for the gain solutions. equation (1) can be written equivalently in the matrix form,

$$\begin{eqnarray*} {\boldsymbol {d}}= {\bf \textsf {G}}{\boldsymbol {m}}+ {\boldsymbol {n}}. \end{eqnarray*}$$

(4)

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

$$\begin{eqnarray*} \chi ^{2} = ({\boldsymbol {d}}-{\bf \textsf {G}}{\boldsymbol {m}})^{\dagger } {\bf \textsf {N}}^{-1} ({\boldsymbol {d}}- {\bf \textsf {G}}{\boldsymbol {m}}). \end{eqnarray*}$$

(5)

Here, |${\bf \textsf {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,

$$\begin{eqnarray*} {\bf \textsf {N}}&=& \langle {\boldsymbol {n}}{\boldsymbol {n}}^{\dagger } \rangle \nonumber \\ &=& \sigma _{pq}^{2} {\bf \textsf {I}}, \end{eqnarray*}$$

(6)

where the angle brackets 〈...〉 denote ensemble average, and in the second equality in equation (6) assumes that the noise matrix is diagonal, i.e. the noise on visibilities of different baselines is uncorrelated. |${\bf \textsf {I}}$| is the identity matrix, and † denotes a complex conjugate transpose. The measured visibility can be corrected (or calibrated) by using the gain solutions estimated from optimization of equation (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 (Barry et al. 2016; Beardsley 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_α, where α indexes the redundant baseline sets (redundant groups),

$$\begin{eqnarray*} d_{pq}= g_{p}g^{\ast }_{q} u_{\alpha } + n_{pq}. \end{eqnarray*}$$

(7)

Since the number of measured visibility d_pq is greater than the number of u_α 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 equation (7) can be then determined by minimizing the χ² function of the form,

$$\begin{eqnarray*} \chi ^{2} = ({\boldsymbol {d}}- {\bf \textsf {G}}{\boldsymbol {u}}_{\alpha })^{\dagger } {\bf \textsf {N}}^{-1} ({\boldsymbol {d}}- {\bf \textsf {G}}{\boldsymbol {u}}_{\alpha }). \end{eqnarray*}$$

(8)

In contrast to equation (5), the redundant visibility parameter u_α here is not pre-defined and will be solved iteratively alongside with the antenna gain factor |${\bf \textsf {G}}\equiv g_{p}g_{q}$|⁠. To minimize the χ² function, a least-square estimator is usually utilized. 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 the least-square solver implementation. Minimizing equation (8) gives us a power spectrum estimation

$$\begin{eqnarray*} \hat{{\boldsymbol {x}}}= [{\bf \textsf {A}}^{\intercal }{\bf \textsf {N}}^{-1}{\bf \textsf {A}}]^{-1} {\bf \textsf {A}}^{\intercal } {\bf \textsf {N}}^{-1} {\boldsymbol {d}}. \end{eqnarray*}$$

(9)

Here, the matrix |${\bf \textsf {A}}$|⁠, which depends on the array configuration, maps the measured visibility d to the calibration parameters |${\bf \textsf {G}}$| and |${\bf \textsf {v}}$| for which we want to solve. The vector |$\hat{{\boldsymbol {x}}}$| contains the least-square estimates for per antenna gain parameter and per-redundant group averaged visibility. Equation (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 OmniCal³ 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 (equation (7)), the least-square estimator of equation (8) will not converge to a unique solution. Here, the term degeneracies refers to the linear combinations of gains and visibilities that redundant calibration cannot solve for, which results in the χ² function 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), and 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², the χ² remains unchanged so that the overall amplitude, A, of the instrument cannot be solved. Secondly, moving the phase of the gain parameter by certain amount Δ, i.e. |$g_{p}=|g_{p}|e^{j\phi _{p}}\rightarrow |g_{p}|e^{j(\phi _{p}+\Delta)}$| does not affect the χ² due to the cancellation of |$g_{p}\times g^{*}_{q}$| in equation (1). The third and fourth degeneracies are caused by phase gradients Δ_x and Δ_y along x and y directions, 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 χ² 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 Δ_x and Δ_y are still degenerate in the relative calibration stage. In a sense, 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 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-modelling error. During the absolute calibration stage, the sky-modelling 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, which is a good approximation (Sievers (2017)). 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 χ². 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 defined in equation (3), its covariance |$\boldsymbol{\Sigma }$| follows,

$$\begin{eqnarray*} \boldsymbol{\Sigma } = \langle {\boldsymbol {d}}{\boldsymbol {d}}^{\dagger } \rangle . \end{eqnarray*}$$

(10)

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

$$\begin{eqnarray*} \boldsymbol{\Sigma } &=& \langle ({\bf \textsf {G}}{\bf \textsf {v}}+ {\boldsymbol {n}}) ({\bf \textsf {G}}{\bf \textsf {v}}+ {\boldsymbol {n}})^{\dagger } \rangle \nonumber \\ &=& {\bf \textsf {G}}\langle {\bf \textsf {v}}{\bf \textsf {v}}^{\dagger }\rangle {\bf \textsf {G}}^{\dagger } + {\bf \textsf {G}}\langle {\bf \textsf {v}}{\boldsymbol {n}}^{\dagger } \rangle + \langle {\boldsymbol {n}}{\bf \textsf {v}}^{\dagger }\rangle {\bf \textsf {G}}^{\dagger } + \langle {\boldsymbol {n}}{\boldsymbol {n}}^{\dagger } \rangle . \end{eqnarray*}$$

(11)

By assuming that the instrumental noise vector n does not correlate with the true sky visibility vector |${\bf \textsf {v}}$|⁠, the middle two terms drop out. The last term describes the noise covariance, which is generally assumed to be Gaussian and uncorrelated between baselines; therefore, it can be written as a noise matrix |${\bf \textsf {N}}$| as defined in equation (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

$$\begin{eqnarray*} {\bf \textsf {C}}= \langle {\bf \textsf {v}}{\bf \textsf {v}}^{\dagger } \rangle , \end{eqnarray*}$$

(12)

equation (11) takes

$$\begin{eqnarray*} \boldsymbol{\Sigma } = {\bf \textsf {G}}{\bf \textsf {C}}{\bf \textsf {G}}^{\dagger } + {\bf \textsf {N}}. \end{eqnarray*}$$

(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 that 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 multidimensional Gaussian distribution,

$$\begin{eqnarray*} \mathcal {L}({\boldsymbol {d}},\boldsymbol{\Sigma }) \propto \exp \left[-\frac{1}{2}{\boldsymbol {d}}^{\dagger }\boldsymbol{\Sigma }^{-1}{\boldsymbol {d}}\right]. \end{eqnarray*}$$

(14)

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

$$\begin{eqnarray*} \mathcal {L}({\boldsymbol {d}},{\bf \textsf {C}}|{\bf \textsf {G}}) \propto \exp \left[-\frac{1}{2}{\boldsymbol {d}}^{\dagger }({\bf \textsf {G}}{\bf \textsf {C}}{\bf \textsf {G}}^{\dagger }+{\bf \textsf {N}})^{-1}{\boldsymbol {d}}\right]. \end{eqnarray*}$$

(15)

Taking the logarithm of both sides of equation (15), we obtain the log-likelihood function |$\log \mathcal {L}$|⁠,

$$\begin{eqnarray*} \log \mathcal {L}({\boldsymbol {d}},{\bf \textsf {C}}|{\bf \textsf {G}}) \propto -\frac{1}{2}{\boldsymbol {d}}^{\dagger}({\bf \textsf {G}}{\bf \textsf {C}}{\bf \textsf {G}}^{\dagger }+{\bf \textsf {N}})^{-1}{\boldsymbol {d}}. \end{eqnarray*}$$

(16)

Maximizing the likelihood function |$\mathcal {L}$| is equivalent to minimizing the (negative) log-likelihood |$\log \mathcal {L}$|⁠. The term on the right-hand side of equation (16) measures the square of the distances from predicted (or true) visibility to the measured visibility in the standard deviation unit. Minimizing this term against |${\bf \textsf {G}}$| is equal to maximizing the likelihood of true visibility being observed. The exponential term is thus simply equal to the χ², given by

$$\begin{eqnarray*} \chi ^{2}={\boldsymbol {d}}^{\dagger }({\bf \textsf {G}}{\bf \textsf {C}}{\bf \textsf {G}}^{\dagger }+{\bf \textsf {N}})^{-1}{\boldsymbol {d}}. \end{eqnarray*}$$

(17)

The covariance-based calculation of χ² in this equation is the underlying mathematical formalism of CorrCal. We further add an additional term χ² → χ² + (∑_iIm(g_i))² + (∑_iRe(g_i) − N_ant)² as a regularization factor to bound the gain solution. Here, the χ² function does not explicitly depend 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 |${\bf \textsf {C}}$| of the true sky visibility.

Since we are modelling the data as Gaussian distributed, a form of the χ² function in equation (17) is different from the one defined in equation (5) and equation (8), which explicitly depends on the sky parameters m and u_α, respectively. However, the CorrCal approach of keeping the array and sky information in the |${\bf \textsf {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 |${\bf \textsf {N}}$| and covariance between visibilities. That is,

$$\begin{eqnarray*} \mathbf {N}_{\rm eff}={\bf \textsf {N}}+(a\mathbf {1})\otimes (a\mathbf {1})^{\dagger }, \end{eqnarray*}$$

(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 ⊗ denotes the outer product of vectors. Applying the Woodbury identity [see equation (24)] to |$\mathbf {N}_{\rm eff}^{-1}$|⁠, and taking a → ∞, the χ² 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 equation (17) into two components as per suggestion in Sievers (2017),

$$\begin{eqnarray*} {\bf \textsf {C}}= {\bf \textsf {S}}+ {\bf \textsf {R}}, \end{eqnarray*}$$

(19)

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

$$\begin{eqnarray*} \chi ^{2}={\boldsymbol {d}}^{\dagger }({\bf \textsf {G}} ({\bf \textsf {S}}+{\bf \textsf {R}}){\bf \textsf {G}}^{\dagger }+{\bf \textsf {N}})^{-1}{\boldsymbol {d}}. \end{eqnarray*}$$

(20)

Note that the explicit formalisms of |${\bf \textsf {S}}$| and |${\bf \textsf {R}}$| will be presented in Section 3.1 in details.

With the redundant group diffuse sky |${\bf \textsf {R}}$| and sources |${\bf \textsf {S}}$| expressed as vector outer products, equation (19) can be rewritten as

$$\begin{eqnarray*} {\bf \textsf {C}}={\bf \textsf {s}}{\bf \textsf {s}}^{\dagger }+{\bf \textsf {r}}{\bf \textsf {r}}^{\dagger }, \end{eqnarray*}$$

(21)

where |${\bf \textsf {r}}$| and |${\bf \textsf {s}}$| are N-dimensional vectors for the diffuse sky and source components, respectively. We substitute equation (21) for |${\bf \textsf {C}}$| in equation (20) to obtain a χ² function

$$\begin{eqnarray*} \chi ^{2}={\boldsymbol {d}}^{\dagger }({\bf \textsf {N}}+{\bf \textsf {G}} ({\bf \textsf {s}}{\bf \textsf {s}}^{\dagger }+{\bf \textsf {r}}{\bf \textsf {r}}^{\dagger }){\bf \textsf {G}}^{\dagger })^{-1}{\boldsymbol {d}}. \end{eqnarray*}$$

(22)

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

$$\begin{eqnarray*} \chi ^{2}={\boldsymbol {d}}^{\dagger } ({\bf \textsf {N}}+\hat{{\bf \textsf {s}}}\hat{{\bf \textsf {s}}}^{\dagger }+\hat{{\bf \textsf {r}}}\hat{{\bf \textsf {r}}}^{\dagger })^{-1}{\boldsymbol {d}}, \end{eqnarray*}$$

(23)

where |$\hat{{\bf \textsf {s}}}={\bf \textsf {G}}{\bf \textsf {s}},\hat{{\bf \textsf {r}}}={\bf \textsf {G}}{\bf \textsf {r}}$|⁠, so |$\hat{{\bf \textsf {s}}}^{\dagger }={\bf \textsf {s}}^{\dagger }{\bf \textsf {G}}^{\dagger }$| and |$\hat{{\bf \textsf {r}}}^{\dagger }={\bf \textsf {r}}^{\dagger }{\bf \textsf {G}}^{\dagger }$|⁠.

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

$$\begin{eqnarray*} ({\bf \textsf {B}}+{\bf \textsf {v}}\otimes {\bf \textsf {v}}^{\dagger})^{-1}={\bf \textsf {B}}^{-1} -{\bf \textsf {B}}^{-1}{\bf \textsf {v}} ({\bf \textsf {I}}+{\bf \textsf {v}}^{\dagger }{\bf \textsf {B}}^{-1}{\bf \textsf {v}})^{-1} {\bf \textsf {v}}^{\dagger }{\bf \textsf {B}}, \end{eqnarray*}$$

(24)

where |${\bf \textsf {B}}$| and the identity matrix |${\bf \textsf {I}}$| are square matrices with dimension N × N, while |${\bf \textsf {v}}$| is N-dimensional vector.

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

First, the Woodburry inversion can be applied to Γ⁻¹, where |$\Gamma = ({\bf \textsf {N}}+\hat{{\bf \textsf {r}}}\hat{{\bf \textsf {r}}}^{\dagger})$| for each redundant group separately, ignoring the source vectors |${\bf \textsf {s}}$|⁠. To apply Woodburry identity to Γ⁻¹, let |${\bf \textsf {B}}={\bf \textsf {N}}$| and |${\bf \textsf {v}}={\bf \textsf {r}}$|⁠. Then, the inverted form of Γ is kept separate for each redundant group and saved in a factored form by applying the Cholesky factorization of the matrix |$[{\bf \textsf {I}}+{\bf \textsf {v}}^{\dagger }{\bf \textsf {B}}^{-1}{\bf \textsf {v}}]^{-1}$| in equation (24) for the redundant vectors |${\bf \textsf {r}}$|⁠.
Finally, using the source vectors |${\bf \textsf {s}}$|⁠, and the already factored form of Γ⁻¹ from the first step, the Woodburry inversion can be computed using |$(\Gamma +\hat{{\bf \textsf {s}}}\hat{{\bf \textsf {s}}}^{\dagger })^{-1}.$|

3.1 Determining the sky and noise variance–covariance matrices

To perform the CorrCal, we first need to estimate the matrices |${\bf \textsf {R}}$|⁠, |${\bf \textsf {S}}$|⁠, and |${\bf \textsf {N}}$|⁠. Thus, in this section, we will derive these variance–covariance matrices explicitly.

3.1.1 Diffuse sky component matrix (⁠|${\bf \textsf {R}}$|⁠)

For an instrument with a small field of view, the covariance matrices can be directly constructed from the visibility in the two-dimensional (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 u = b/λ that expresses the baseline vector b between antennas p and q in wavelength λ, the measured cross-correlation signal from these antennas, or the visibility, ignoring noise term, is given by,

$$\begin{eqnarray*} v_{pq}({\boldsymbol {u}}| f) = \int \textrm {d}\Omega \, I(\hat{{\boldsymbol {r}}}, f) B(\hat{{\boldsymbol {r}}}, f), \end{eqnarray*}$$

(25)

where

$$\begin{eqnarray*} B(\hat{{\boldsymbol {r}}},f) &=& A_{p}(\hat{{\boldsymbol {r}}},f)A^{\ast }_{q}(\hat{{\boldsymbol {r}}}, f) \exp \left[ -2\mathrm{i}\pi {\boldsymbol {u}}\cdot \hat{{\boldsymbol {r}}}\right]\nonumber \\ &=&|A_{p}(\hat{{\boldsymbol {r}}},f)|^{2}\exp \left[ -2\mathrm{i}\pi {\boldsymbol {u}}\cdot \hat{{\boldsymbol {r}}}\right]. \end{eqnarray*}$$

(26)

In this equation, |$I(\hat{{\boldsymbol {r}}}, f)$| is the flux density measured over a solid angle Ω subtended by the observed region of the sky in the direction of a unit vector |$\hat{{\boldsymbol {r}}}$|⁠. |$A_{p}(\hat{{\boldsymbol {r}}}, f)$| and |$A_{q}(\hat{{\boldsymbol {r}}}, f)$| are the primary beam of antenna p and q, respectively. Here, the squaring of the antenna beam term |$A_{p}(\hat{{\boldsymbol {r}}},f)$| in the second equality of equation (26) comes from the fact that we assume that 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{{\boldsymbol {r}}}, f)$|⁠) and array (⁠|$B(\hat{{\boldsymbol {r}}}, f)$|⁠) over the spherical harmonics, we have

$$\begin{eqnarray*} I(\hat{{\boldsymbol {r}}}, f)=\sum _{\ell m}a_{\ell m}(f)Y_{\ell m}(\hat{{\boldsymbol {r}}}), \end{eqnarray*}$$

(27)

$$\begin{eqnarray*} B(\hat{{\boldsymbol {r}}}, f)=\sum _{\ell m}b_{\ell m}(f)Y_{\ell m}(\hat{{\boldsymbol {r}}}), \end{eqnarray*}$$

(28)

where a_ℓm and b_ℓm are the amplitudes of the flux density and the beam function for the spherical harmonics Y_ℓm, respectively. The angular number ℓ = 0, 1, 2, 3, …, and the azimuthal number m = −ℓ, −ℓ + 1, ⋅⋅⋅, ℓ − 1, ℓ follow standard convention, allowing 2ℓ + 1 values of m for each value of ℓ. Using equation (27) and equation (28), the visibility in equation (25) can be expressed in spherical harmonics as,

$$\begin{eqnarray*} v_{pq}(f)&=& \int \textrm {d}\Omega \, \sum _{\ell m}a^{\ast }_{\ell m}(f)Y^{\ast }_{\ell m}(\hat{{\boldsymbol {r}}}) \sum _{\ell ^{\prime } m^{\prime }}b_{\ell ^{\prime } m^{\prime }}(f)Y_{\ell ^{\prime } m^{\prime }}(\hat{{\boldsymbol {r}}}) \nonumber \\ &=& \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{{\boldsymbol {r}}})Y_{\ell ^{\prime } m^{\prime }}(\hat{{\boldsymbol {r}}}) \nonumber \\ & =& \sum _{\ell m}a^{\ast }_{\ell m}(f)b_{\ell m}(f), \end{eqnarray*}$$

(29)

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

We calculated the array function that is defined in equation (26) using the exponential argument (⁠|${\boldsymbol {b}}\cdot \hat{{\boldsymbol {r}}}{}/c$|⁠) in equation (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 –168 MHz, we use the Astronomical Interferometry in Python (AIPY),⁴ a software package primarily developed for PAPER/HERA data analysis. Fig. 2 shows that the AIPY-simulated PAPER-64 beam model projected on the HEALPIx⁵ coordinates at three arbitrary right accessions for frequency of 150 MHz. We will use this beam model for all calculations of the covariance matrices.

$A delay term (${\boldsymbol {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 are shown in Fig. 8.$

Figure 1.

A delay term (⁠|${\boldsymbol {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 are shown in Fig. 8.

Open in new tab Download slide

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^○S) at Karoo desert, South Africa.

Open in new tab Download slide

Using the expected visibility function defined in equation (29), in this section, we construct the diffuse sky covariance matrix |${\bf \textsf {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

$$\begin{eqnarray*} {\bf \textsf {R}}_{({{\boldsymbol {u}}}f);({{\boldsymbol {u}}}^{\prime }f^{\prime })} &=& \langle v^{\rm diff}({\boldsymbol {u}},f)v^{{\rm diff},*}({\boldsymbol {u}}^{\prime },f^{\prime })\rangle \nonumber \\ &= & \left\langle \sum _{\ell m}a^{{\rm diff},\ast }_{\ell m}(f)b^{{\rm diff}}_{\ell m}({{\boldsymbol {u}}},f) \right. \nonumber \\ &&\times \, \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 }}({{\boldsymbol {u}}}^{\prime },f^{\prime }))\right\rangle \nonumber \\ &= & \sum _{\ell m}\sum _{\ell ^{\prime }\mathrm{m}^{\prime }} b^{{\rm diff}}_{\ell m}({{\boldsymbol {u}}},f)\big\langle a^{{\rm diff},\ast }_{\ell m}(f)a^{{\rm diff}}_{\ell ^{\prime }\mathrm{m}^{\prime }}(f^{\prime })\big\rangle \nonumber \\ &&\times \, b^{{\rm diff},\ast }_{\ell ^{\prime }\mathrm{m}^{\prime }}({{\boldsymbol {u}}}^{\prime },f^{\prime }), \end{eqnarray*}$$

(30)

where |$^{\rm `diff^{\prime }}$| is to denote that |${\bf \textsf {v}}_{pq}(f)$| is coming from the diffuse sky component and the angle bracket 〈...〉 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

$$\begin{eqnarray*} \big\langle a^{{\rm diff}}_{\ell m}(f)a_{\ell ^{\prime }\mathrm{m}^{\prime }}^{{\rm diff},*}(f^{\prime })\big\rangle &=& \langle |a_{\ell m}(f,f^{\prime })|^{2}\rangle \delta _{\ell \ell ^{\prime }}\delta _{\mathrm{mm}^{\prime }} \nonumber \\ &=& C^{{\rm diff}}_{\ell }(f,f^{\prime })\delta _{\ell \ell ^{\prime }}\delta _{mm^{\prime }}, \end{eqnarray*}$$

(31)

where |$C^{{\rm diff}}_{\ell }(f,f^{\prime })=\langle |a^{{\rm diff}}_{\ell m}(f,f^{\prime })|^{2}\rangle$| is the diffuse sky angular power spectrum, which depends only on multipole number ℓ that corresponds to the angular scale θ = 180^○/ℓ, and δ 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 ℓ, not m. Therefore, there is no correlation in the m modes. Thus, the covariance matrix [equation (30)] can be simplified to

$$\begin{eqnarray*} {\bf \textsf {R}}=\sum _{\ell m}b^{{\rm diff}}_{\ell m}({{\boldsymbol {u}}},f)b^{{\rm diff},*}_{\ell m}({{\boldsymbol {u}}}^{\prime },f^{\prime })C^{{\rm diff}}_{\ell }(f,f^{\prime }). \end{eqnarray*}$$

(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 equation (26). The measured antenna positions (instead of assuming that they are perfectly redundant) are integrated into this equation. Hence, through this expression, the covariance-based calculation of χ² in equation (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 is about 30 m.

$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 30-m 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.$

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 30-m 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.

Open in new tab Download slide

In equation (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 |${\bf \textsf {R}}$| had formed to estimate |$C^{{\rm diff}}_{\ell }(f,f^{\prime })$|⁠. The data covariance matrix can be determined from

$$\begin{eqnarray*} {\bf \textsf {D}}^{\alpha } =\langle {\boldsymbol {d}}^{\alpha }({{\boldsymbol {u}}},f){\boldsymbol {d}}^{\alpha \dagger }({{\boldsymbol {u}}}^{\prime },f^{\prime })\rangle , \end{eqnarray*}$$

(33)

where |${\bf \textsf {D}}^{\alpha }$| is a covariance matrix of measured visibility data d^α from the redundant group α. To determine the |$C^{\rm diff}_{\ell }(f,f^{\prime })$|⁠, thus, we equate the covariance matrix |${\bf \textsf {D}}^{\alpha }$| with the diffuse sky matrix |${\bf \textsf {R}}$| from the same redundant group α in equation (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 |${\bf \textsf {D}}$| with |${\bf \textsf {R}}$| (i.e. |${\bf \textsf {R}}^{\alpha }\approx {\bf \textsf {D}}^{\alpha }$|⁠) if both are computed from the same redundant group α. This assumption may avoid mode-mixing between baselines with different lengths. It follows that,

$$\begin{eqnarray*} C^{{\rm diff}}_{\ell }(f,f^{\prime })\approx \frac{{\bf \textsf {D}}^{\alpha }}{\sum _{\ell m}b^{{\rm diff}}_{\ell m}({{\boldsymbol {u}}},f)b^{{\rm diff}}_{\ell m}({{\boldsymbol {u}}}^{\prime },f^{\prime })}. \end{eqnarray*}$$

(34)

To pull |$C^{\rm diff}_\ell (f,f^{\prime })$| out of the summation sign in equation (32), we assume that for a given redundant block α, the sky is slowly changing function of ℓ ≈ 2π|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 equation (34). Then, we implement the power-law curve fitting into equation (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 mK²) and the spectral index α, the fitted power spectrum takes

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

(35)

where the fitted parameters |$A=249.69\, {\rm mK^{2}}$| and α = 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 ℓ.

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 equation (32) to estimate the model covariance matrix of the group. Likewise, one needs to derive the best-fitting 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 |${\bf \textsf {R}}$| finally takes

$$\begin{eqnarray*} {\bf \textsf {R}}=C^{\rm diff}_{{\rm fits}}(f,f^{\prime })\sum _{\ell }b^{{\rm diff}}_{\ell }({{\boldsymbol {u}}},f)b^{{\rm diff},*}_{\ell }({{\boldsymbol {u}}}^{\prime },f^{\prime }). \end{eqnarray*}$$

(36)

$Left: Cross-frequency sky power spectrum $C_{\ell }^{{\rm diff}}(f,f^{\prime })$ obtained from equation (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 that the ℓ-dependence of $C^{\rm diff}_{\ell }(f,f^{\prime })$ is the same. The whiteout regions are corresponding to the flagged frequency channels from the data.$

Figure 4.

Left: Cross-frequency sky power spectrum |$C_{\ell }^{{\rm diff}}(f,f^{\prime })$| obtained from equation (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 that the ℓ-dependence of |$C^{\rm diff}_{\ell }(f,f^{\prime })$| is the same. The whiteout regions are corresponding to the flagged frequency channels from the data.

Open in new tab Download slide

Figure 5.

Noise variance estimated from the system temperature of PAPER.

Open in new tab Download slide

We use equation (36) to represent the diffuse sky component for χ² optimization process in equation (17). The dimension of the covariance matrix |${\bf \textsf {R}}$| depends on the number of baseline N_u and number of frequency channels N_f, i.e. its dimension is equal to N_uN_f × N_uN_f.

3.1.2 Point source component matrix (⁠|${\bf \textsf {S}}$|⁠)

The partial information about the point sources on the sky is one of the components the CorrCal calibration scheme employs to estimate 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 |${\bf \textsf {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{{\boldsymbol {r}}}_{i}$| thus comes from source motion relative to the observer meridian, which best visualized 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 |${\bf \textsf {S}}$|⁠. We then project each source position on to 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 equation (26) for known sources. That is, suppose the point source profile is a Dirac-delta function on the sky |$\hat{{\boldsymbol {r}}}_{i}$| (i = 1, 2,..., N_ps, N_ps is the total number of point sources), then the specific intensity [equation (25)] can be written as

$$\begin{eqnarray*} I(\hat{{\boldsymbol {r}}}, f) =\sum ^{N_{\rm ps}}_{i=1}I(\hat{{\boldsymbol {r}}}_{i},f)\delta _{\rm 2D}\left(\hat{{\boldsymbol {r}}},\hat{{\boldsymbol {r}}}_{i} \right), \end{eqnarray*}$$

(37)

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

$$\begin{eqnarray*} {\bf \textsf {v}}^{\rm ps}({\boldsymbol {u}},f)=\sum ^{N_{\rm ps}}_{i=1}I(\hat{{\boldsymbol {r}}}_{i},f)\left|A_{pq}(\hat{{\boldsymbol {r}}}_{i},f) \right|^{2}{\rm e}^{-{\rm i}2\pi {\boldsymbol {u}}\cdot \hat{{\boldsymbol {r}}}_{i}}. \end{eqnarray*}$$

(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 is displayed in Fig. 2. Taking the outer product between the visibility vector |${\bf \textsf {v}}^{\rm ps}({\boldsymbol {u}},f)$| and its complex conjugate transpose |${\bf \textsf {v}}^{\rm ps,\dagger }({\boldsymbol {u}},f)$|⁠, we form the covariance information for the discrete sources,

$$\begin{eqnarray*} {\bf \textsf {S}}={\bf \textsf {v}}^{\rm ps}({\boldsymbol {u}},f){\bf \textsf {v}}^{\rm ps,\dagger }({\boldsymbol {u}},f). \end{eqnarray*}$$

(39)

Thus, we will use this source information covariance matrix to calculate the χ² function in equation (20).

3.1.3 Noise matrix (⁠|${\bf \textsf {N}}$|⁠)

In this work, we estimate the per-diagonal noise variance matrix |${\bf \textsf {N}}$| from the radiometer equation

$$\begin{eqnarray*} \sigma _{\rm rms}=\frac{T_{\rm sys}}{\sqrt{2\Delta f \tau }}, \end{eqnarray*}$$

(40)

where |$\Delta f=97\, {\rm kHz}$| is the frequency bandwidth in kHz, τ = 49.2 s is the integration time. The system temperature T_sys from the relation defined in equation (23) of Cheng et al. (2018) is,

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

(41)

where f is frequency, and T_rcvr is a receiver temperature which equals to 144 K. The noise variance is shown in Fig. 5.

3.2 Eigendecomposition of sky signal covariance matrix

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

$$\begin{eqnarray*} {\bf \textsf {R}}=C_{\ell }^{2} {\begin{bmatrix}1 &\quad 1 \\ 1 &\quad 1 \\ \end{bmatrix}} , \end{eqnarray*}$$

(42)

However, for a quasi-redundant group, all the off-diagonal elements of the correlation matrix |${\bf \textsf {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 |${\bf \textsf {R}}$| from its eigendecomposition in our analysis following Sievers (2017). 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 eigendecomposition analysis, few modes of the principal components accounted for most of the variances in the original data. Thus, those eigenmodes of |${\bf \textsf {R}}$| corresponding to sufficiently large eigenvalues could replace its original data with very little loss of information. In the left-hand panel of Fig. 6, we show the eigenvalues of cross-frequency covariance matrix described in equation (36). The vertical axis represents the log₁₀ eigenvalues, while the horizontal axis labels the corresponding principal component numbers. The sharp falloff 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 |${\bf \textsf {R}}$|⁠. In addition, all the eigenvalues are greater than zero because |${\bf \textsf {R}}$| is a positive-definite covariance matrix. In the right-hand panel of Fig. 6, we depict the first three eigenmodes for |${\bf \textsf {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.

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

Figure 6.

Left: Eigenvalues of cross-frequency covariance matrix described in equation (36). The vertical axis represents the log₁₀ eigenvalues, while the horizontal axis labels the corresponding principal component index. Right: The first three eigenmodes for |${\bf \textsf {R}}$| corresponding to the first three eigenvalues for the frequency range from 120 to 168 MHz.

Open in new tab Download slide

We keep eigenmodes with amplitude more than 10⁻⁵ times the largest. This allows the keeping of ∼2–3 complex eigenmodes per quasi-redundant block. These few modes could represent most of the variations in a given redundant group α. For instance, the proportion of variance related to the eigenvalues λ₁ and λ₂ can be calculated from

$$\begin{eqnarray*} \left(\frac{\lambda _{1}+\lambda _{2}}{\lambda _{1}+\lambda _{2}+\cdots +\lambda _{p}}\right)\times 100, \end{eqnarray*}$$

(43)

where p is the total eigenmodes. If the value of equation (43) is greater than the threshold percentage |$99{{\ \rm per\ cent}}$|⁠, we use eigenmodes corresponding to eigenvalues λ₁ and λ₂ 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

$$\begin{eqnarray*} {\bf \textsf {p}}=\mathbf {e}_{i}\sqrt{\lambda _{i}}, \end{eqnarray*}$$

(44)

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

$$\begin{eqnarray*} {\bf \textsf {R}}_{\rm PCA}={\bf \textsf {p}}\otimes {\bf \textsf {p}}^{\dagger }, \end{eqnarray*}$$

(45)

where |${\bf \textsf {R}}_{\rm PCA}$| is the diffuse sky covariance matrix calculated from the principal component analysis of |${\bf \textsf {R}}$| for a given redundant block α. The dimension of |${\bf \textsf {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 |${\bf \textsf {R}}$| for a given redundant group containing four baselines, a complex vector |${\bf \textsf {p}}$| could be

$$\begin{eqnarray*} {\bf \textsf {p}}=[p_{1},p_{2},0,0], \end{eqnarray*}$$

(46)

where p₁ and p₂ are complex elements of a vector |${\bf \textsf {p}}$| that retained from eigenmode analysis. The outer product of this vector by its complex conjugate transpose takes

$$\begin{eqnarray*} {\bf \textsf {p}}\otimes {\bf \textsf {p}}^{\dagger }={\begin{bmatrix} p_{1}p^{*}_{1} &\quad p_{1}p^{*}_{2} &\quad 0 & \quad 0 \\ p_{2}p^{*}_{1} &\quad p_{2}p^{*}_{2} &\quad 0 & \quad 0\\ 0 &\quad 0 &\quad 0 & \quad 0 \\ 0 &\quad 0 &\quad 0 & \quad 0\\ \end{bmatrix}}. \end{eqnarray*}$$

(47)

Here, the operator * denotes a complex conjugate. Adding the noise to the above matrix, we have

$$\begin{eqnarray*} \Gamma ={\bf \textsf {N}}+{\bf \textsf {p}}\otimes {\bf \textsf {p}}^{\dagger }. \end{eqnarray*}$$

(48)

All off-diagonal elements of the noise variance matrix |${\bf \textsf {N}}$|⁠, which has the same size as |${\bf \textsf {p}}\otimes {\bf \textsf {p}}^{\dagger }$|⁠, are zero. Thus, Γ has only a relatively small number of non-zero 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 Γ since most of the O(N³) operations devoted to inverting the matrix involve zero operands. Hence, in terms of memory storage and computational efficiency, performing the sparse matrices operation traversing only non-zero elements has significant advantages over their dense matrix counterparts. Taking advantage of the sparse structure of the correlation matrix Γ, CorrCal utilizes the sparse matrix operations to increase its computational performance in a large data set.

To implement the sparse matrix computations in CorrCal, per-redundant group correlation matrices Γ are stored as 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 |${\bf \textsf {R}}$|⁠, |${\bf \textsf {S}}$|⁠, and |${\bf \textsf {N}}$| to each of these redundant groups. Then, one needs to set up the sparse matrix representations for these matrices and determine the χ² function and its gradient information using the CorrCal routines. Utilizing these ingredients as inputs, CorrCal finally adopts the Scipy⁶ built-in optimization algorithm to fit antenna gains. We summarize the procedures of implementation of CorrCal as follows:

First, we need to create the sparse matrices representing |${\bf \textsf {R}}$|⁠, |${\bf \textsf {S}}$|⁠, and |${\bf \textsf {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:
- The real-imaginary separated per-redundant group diffuses sky vector |${\bf \textsf {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 equation (46)].
- The real-imaginary separated vector |${\bf \textsf {s}}$| that carries information about the known point sources.
- Per-diagonal noise variance matrix |${\bf \textsf {N}}$|⁠.
- A vector that contains indices setting off the redundant group.
Secondly, the χ² function and its gradient information should be determined using the CorrCal subroutines. The arguments passed to these sub-routines in order to calculate the χ² function and the gradient information of χ² are:
- A vector containing the real-imaginary separated initial guess for antenna gains.
- Per-redundant group measured visibility data vector. The real-imaginary parts should be separated such that imaginaries follow their respective reals [e.g. see vectorized visibility in equation (49)].
- The sparse matrix representations for |${\bf \textsf {R}}$|⁠, |${\bf \textsf {S}}$|⁠, and |${\bf \textsf {N}}$| from the first step.
- A vector containing per-visibility antenna indices.
Finally, using the χ², 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 non-linear 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 towards 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 downweighting 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 fit 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, visibility data d from different frequency channels put together in a longer vector as,

$$\begin{eqnarray*} {\boldsymbol {d}}= {\begin{bmatrix}\operatorname{\mathbb {R}e}\lbrace \mathrm{v({f_{1})}}\rbrace ,\ \operatorname{\mathbb {I}m}\lbrace \mathrm{v(f_{1})}\rbrace ,\ \operatorname{\mathbb {R}e}\lbrace \mathrm{v(f_{2}})\rbrace ,\ \operatorname{\mathbb {I}m}\lbrace \mathrm{v(f_{2}})\rbrace ,\ \cdots \end{bmatrix}}, \end{eqnarray*}$$

(49)

where |$\operatorname{\mathbb {R}e}\lbrace \mathrm{v}(f_{1})\rbrace$|⁠, |$\operatorname{\mathbb {I}m}\lbrace \mathrm{v}(f_{1})\rbrace$|⁠, |$\operatorname{\mathbb {R}e}\lbrace \mathrm{v(f_{2}})\rbrace$|⁠, and |$\operatorname{\mathbb {I}m}\lbrace \mathrm{v(f_{2})}\rbrace$| are the real (⁠|$\operatorname{\mathbb {R}e}$|⁠) and imaginary (⁠|$\operatorname{\mathbb {I}m}$|⁠) parts of the visibility measurements corresponding to frequency channel f₁ and f₂, respectively. Concatenating visibility from all frequency channels in a longer vector like equation (49) enforces 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 10-min 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 affected by RFI contamination.

Open in new tab Download slide

The correlation calibration suffers from the same degeneracy problem as with redundant baseline calibration 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 1. Because we are applying CorrCal to data that have already been absolutely calibrated, we do not want to rescale the overall amplitude. To set the overall phase degenerate parameter Δ, 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,

$$\begin{eqnarray*} \chi^{2}_{\text{phase}} = \sum _{p=1}^{N}(\text{Arg}[\hat{g}_{p}(f)]-\Delta_{x}b_{x,p}-\Delta _{y}b_{y}, p), \end{eqnarray*}$$

(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-modelled sky catalogue 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 that the source catalogue is complete and well modelled with precisely known source positions. The method works well if the sky-modelling error could be mostly caused by uncertain/missed source fluxes in the sky-based calibration approach.

The covariance-based formulation of χ² 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. equation (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 and 200 MHz, corresponding to the 21-cm signal from redshifts 6–13, with a frequency resolution of ∼97 kHz. Data from the 32-elements PAPER array have been analysed 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) and Kolopanis et al. (2019) have placed significant upper limits on the power spectrum amplitude of 21-cm signal.

We will perform a re-calibration 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 these data has been fixed. The detailed discussions about the absolute calibration stage for these data are presented in Ali et al. (2015). The data set has dual polarization products. We use one night of observations taken between 2012 November 10 (JD2456242.17382) and 2012 November 11 (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 × (64 − 1)/2 = 2016 and thus measure 2016 visibilities. For our calibration, we use all baselines except those associated with three antennas that 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).

Open in new tab Download slide

Table 1.

Open in new tab

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	2012 November 10 to 2012 November 11
Observing mode	Drift-scan
Time resolution	∼42.9 s
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	4 m
Longest baseline	212 m

Parameter	Value
PAPER array location	\|$\mathrm{30.7^{o} S, 21.4^{o}\, E}$\|
Observation dates	2012 November 10 to 2012 November 11
Observing mode	Drift-scan
Time resolution	∼42.9 s
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	4 m
Longest baseline	212 m

Table 1.

Open in new tab

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	2012 November 10 to 2012 November 11
Observing mode	Drift-scan
Time resolution	∼42.9 s
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	4 m
Longest baseline	212 m

Parameter	Value
PAPER array location	\|$\mathrm{30.7^{o} S, 21.4^{o}\, E}$\|
Observation dates	2012 November 10 to 2012 November 11
Observing mode	Drift-scan
Time resolution	∼42.9 s
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	4 m
Longest baseline	212 m

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. Fig. 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).

Open in new tab Download slide

5.2 Delay-space power spectra

The statistical analysis of the Fourier mode of the power spectrum (k_⊥, k_∥) for low-frequency observations is a powerful tool to probe the EoR (Morales & Hewitt 2004). This method shows the existence of distinctive regions dominated by 21-cm signal and foreground contaminants in a 2D line of sight (k_∥) and transverse co-moving |$\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_∥ values, while the vast majority of k_∥ 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_∥ values because of chromatic interaction of interferometry with intrinsically smooth foreground (Morales et al. 2012; Parsons et al. 2012; Vedantham, Shankar & Subrahmanyan 2012; Thyagarajan et al. 2013; Dillon et al. 2014; Liu et al. 2014a; Liu, Parsons & Trott 2014b). This effect leads to the distinctive wedge-like structure in k_⊥ and k_∥ Fourier space, leaving the foreground-free EoR window beyond the wedge (see results in Figs. 10 and 11 and discussions there).

To convert the measured visibility to 2D power spectrum in k_∥ and k_⊥ 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 visibilities⁷ 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-h observation over frequencies 140–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 one-dimensional (1D) CLEAN algorithm on the delay-transformed data to reduce the effect introduced by RFI flagging and to increase 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 equation (A5).

Fig. 10 shows the delay power spectra where the positive and negative values of k_∥ 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_∥, respectively. The delay power spectra from emissions near-zero delays (where k_∥ ∼ 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.

Figure 10.

Delay power spectra on either side of k_∥ 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 re-calibration. 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 re-calibration.

Open in new tab Download slide

However, from Fig. 10, we see that the strength of power spectra around zero delay decreases as one moves from the lower k_⊥ to higher k_⊥ regions. Furthermore, for both calibration approaches, more emission goes beyond the horizon limit for shortest baselines (smallest k_⊥ 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_∥ (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 equation (A7). The spectral-smooth emissions are supposed to be confined within that limit, and any emissions that 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; Li et al. 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 phenomenon suggests a reduction in the spectral structure of the data after re-calibration with CorrCal.

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_∥ directions and folding it over the centre of the zero-delay line (where k_∥ = 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 Fig. 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, Figs 11(a) and (b), respectively, display the wedge-shaped power spectra to data already calibrated using OmniCal and to data after CorrCal re-calibration. 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_∥ at a longer baseline (k_⊥) regime, demonstrating how the theoretical wedge limit increases as a function of baseline length. In addition, sources that are located far away from the pointing centre of the primary beam will also appear at higher delays, corresponding to higher k_∥. The relative deviation of power in percent between CorrCal and OmniCal calibration approaches is 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{{\ \rm per\ cent}}}$| and |$\mathrm{[(C-A)/A]\times 100{{\ \rm per\ cent}}}$| of Fig. 11 indicates that the power spectra around the wedge limit have become slightly tighter after CorrCal re-calibration. This effect is not considerably saturated if we include the sources’ information in CorrCal formalism as shown in the right-hand panel of the figure.

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 re-calibrated using without source information in CorrCal; (C) power spectra after re-calibration 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 (τ_H).

Open in new tab Download slide

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_⊥ axis over the range of |$0.0733\, h\, {\rm Mpc}^{-1} \lt k_{\perp } \lt 0.0736\, h\, {\rm Mpc}^{-1}$| to generate the corresponding 1D power spectra shown on the L.H.S of Fig. 12. From this figure, we observe that there is a steep falloff 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 falloff follows the same trend for both calibration schemes. However, the inclusion of sky sources statistics as a prior during CorrCal re-calibration step does not show the significant change in the power spectra, as shown on the L.H.S of Fig. 12.

$Right: Two-dimensional power spectra after (A) OmniCal, (B) CorrCal re-calibration without sources info, and (C) CorrCal re-calibration 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 (τ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∥ for a portion of baselines ($\sim 150\rm m$) to both calibration strategies. The bottom panel is the percentage deviation.$

Figure 12.

Right: Two-dimensional power spectra after (A) OmniCal, (B) CorrCal re-calibration without sources info, and (C) CorrCal re-calibration 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 (τ_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_∥ for a portion of baselines (⁠|$\sim 150\rm m$|⁠) to both calibration strategies. The bottom panel is the percentage deviation.

Open in new tab Download slide

6 CONCLUSIONS

In this work, we used a new calibration scheme (CorrCal) to re-calibrate the PAPER64 observations. Our new calibration scheme boils down to the covariance-based calculation of χ², 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{{\ \rm per\ cent}}$| 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 21-cm signal.

ACKNOWLEDGEMENTS

We would like to acknowledge the PAPER collaboration for providing us data that have been used in this work. YZM acknowledges the support of NRF-120385, NRF-120378, and NSFC-11828301. PK 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.

Footnotes

1

The prototype version of the software located on https://github.com/sievers/corrcal2.git is being used for this work.

2

It is worth noticing that Liu et al. (2010) investigated the Taylor expansion of the sky to deal with direction dependence or non-redundancy.

3

https://github.com/jeffzhen/omnical.git

4

https://github.com/HERA-Team/aipy.git

5

https://github.com/healpy/healpy.git

6

https://github.com/scipy/scipy.git

7

Note that these are not true Stokes I would rather the average of XX and YY visibilities (see Section A).

8

https://github.com/HERA-Team/hera_sim.git

REFERENCES

Ali

Z. S.

et al. ,

2015

,

ApJ

,

809

,

61

10.1088/0004-637X/809/1/61

Crossref

Search ADS

Barkana

R.

,

Loeb

A.

,

2001

,

Phys. Rep.

,

349

,

125

Barry

N.

,

Hazelton

B.

,

Sullivan

I.

,

Morales

M. F.

,

Pober

J. C.

,

2016

,

MNRAS

,

461

,

3135

10.1093/mnras/stw1380

Crossref

Search ADS

Beardsley

A. P.

et al. ,

2016

,

ApJ

,

833

,

102

10.3847/1538-4357/833/1/102

Crossref

Search ADS

Bernardi

G.

et al. ,

2009

,

A&A

,

500

,

965

10.1051/0004-6361/200911627

Crossref

Search ADS

Bowman

J. D.

,

Morales

M. F.

,

Hewitt

J. N.

,

2009

,

ApJ

,

695

,

183

10.1088/0004-637X/695/1/183

Crossref

Search ADS

Byrne

R.

et al. ,

2019

,

ApJ

,

875

,

70

10.3847/1538-4357/ab107d

Crossref

Search ADS

Chapman

E.

,

Jelić

V.

,

2019

,

preprint (arXiv:1909.12369)

Cheng

C.

et al. ,

2018

,

ApJ

,

868

,

26

10.3847/1538-4357/aae833

Crossref

Search ADS

Datta

A.

,

Bowman

J. D.

,

Carilli

C. L.

,

2010

,

ApJ

,

724

,

526

10.1088/0004-637X/724/1/526

Crossref

Search ADS

DeBoer

D. R.

et al. ,

2017

,

PASP

,

129

,

045001

10.1088/1538-3873/129/974/045001

Crossref

Search ADS

Dillon

J. S.

et al. ,

2014

,

Phys. Rev. D

,

89

,

023002

10.1103/PhysRevD.89.023002

Crossref

Search ADS

Dillon

J. S.

et al. ,

2018

,

MNRAS

,

477

,

5670

10.1093/mnras/sty1060

Crossref

Search ADS

Dillon

J. S.

et al. ,

2020

,

MNRAS

,

499

,

5840

Dillon

J. S.

,

Liu

A.

,

Tegmark

M.

,

2013

,

Phys. Rev. D

,

87

,

043005

10.1103/PhysRevD.87.043005

Crossref

Search ADS

Edler

H. W.

,

de Gasperin

F.

,

Rafferty

D.

,

2021

,

A&A

,

652

,

A37

10.1051/0004-6361/202140465

Crossref

Search ADS

Ewall-Wice

A.

,

Dillon

S. J.

,

Liu

A.

,

Hewitt

J.

,

2017

,

MNRAS

,

470

,

1849

10.1093/mnras/stx1221

Crossref

Search ADS

Furlanetto

S. R.

,

Oh

S. P.

,

Briggs

F. H.

,

2006

,

Phys. Rep.

,

433

,

181

Hothi

I.

et al. ,

2021

,

MNRAS

,

500

,

2264

10.1093/mnras/staa3446

Crossref

Search ADS

Hurley-Walker

N.

et al. ,

2017

,

MNRAS

,

464

,

1146

10.1093/mnras/stw2337

Crossref

Search ADS

Jacobs

D. C.

et al. ,

2013

,

ApJ

,

776

,

108

10.1088/0004-637X/776/2/108

Crossref

Search ADS

Kohn

S. A.

et al. ,

2016

,

ApJ

,

823

,

88

10.3847/0004-637X/823/2/88

Crossref

Search ADS

Kolopanis

M.

et al. ,

2019

,

ApJ

,

883

,

133

10.3847/1538-4357/ab3e3a

Crossref

Search ADS

Li

W.

et al. ,

2018

,

ApJ

,

863

,

170

10.3847/1538-4357/aad3c3

Crossref

Search ADS

Li

W.

et al. ,

2019

,

ApJ

,

887

,

141

10.3847/1538-4357/ab55e4

Crossref

Search ADS

Liu

A.

,

Parsons

A. R.

,

2016

,

MNRAS

,

457

,

1864

10.1093/mnras/stw071

Crossref

Search ADS

Liu

A.

,

Tegmark

M.

,

2011

,

Phys. Rev. D

,

83

,

103006

10.1103/PhysRevD.83.103006

Crossref

Search ADS

Liu

A.

,

Tegmark

M.

,

Bowman

J.

,

Hewitt

J.

,

Zaldarriaga

M.

,

2009

,

MNRAS

,

398

,

401

10.1111/j.1365-2966.2009.15156.x

Crossref

Search ADS

Liu

A.

,

Tegmark

M.

,

Scott

M.

,

Lutomirski

A.

,

Zaldarriaga

M.

,

2010

,

MNRAS

,

408

,

1029

10.1111/j.1365-2966.2010.17174.x

Crossref

Search ADS

Liu

A.

,

Parsons

A. R.

,

Trott

C. M.

,

2014a

,

Phys. Rev. D

,

90

,

023018

10.1103/PhysRevD.90.023018

Crossref

Search ADS

Liu

A.

,

Parsons

A. R.

,

Trott

C. M.

,

2014b

,

Phys. Rev. D

,

90

,

023019

10.1103/PhysRevD.90.023019

Crossref

Search ADS

Madau

P.

,

Ferguson

H. C.

,

Dickinson

M. E.

,

Giavalisco

M.

,

Steidel

C. C.

,

Fruchter

A.

,

1996

,

MNRAS

,

283

,

1388

10.1093/mnras/283.4.1388

Crossref

Search ADS

Mertens

F. G.

et al. ,

2020

,

MNRAS

,

493

,

1662

10.1093/mnras/staa327

Crossref

Search ADS

Moore

D. F.

,

Aguirre

J. E.

,

Parsons

A. R.

,

Jacobs

D. C.

,

Pober

J. C.

,

2013

,

ApJ

,

769

,

154

10.1088/0004-637X/769/2/154

Crossref

Search ADS

Morales

M. F.

,

Hewitt

J.

,

2004

,

ApJ

,

615

,

7

10.1086/424437

Crossref

Search ADS

Morales

M. F.

,

Wyithe

J. S. B.

,

2010

,

ARA&A

,

48

,

127

10.1146/annurev-astro-081309-130936

Crossref

Search ADS

Morales

M. F.

,

Hazelton

B.

,

Sullivan

I.

,

Beardsley

A.

,

2012

,

ApJ

,

752

,

137

10.1088/0004-637X/752/2/137

Crossref

Search ADS

Oh

S. P.

,

2001

,

ApJ

,

553

,

499

10.1086/320957

Crossref

Search ADS

Orosz

N.

,

Dillon

J. S.

,

Ewall-Wice

A.

,

Parsons

A. R.

,

Thyagarajan

N.

,

2018

,

MNRAS

,

487

,

537

Paciga

G.

et al. ,

2011

,

MNRAS

,

413

,

1174

10.1111/j.1365-2966.2011.18208.x

Crossref

Search ADS

Parsons

A. R.

et al. ,

2010

,

AJ

,

139

,

1468

10.1088/0004-6256/139/4/1468

Crossref

Search ADS

Parsons

A. R.

et al. ,

2014

,

ApJ

,

788

,

106

10.1088/0004-637X/788/2/106

Crossref

Search ADS

Parsons

A. R.

,

Backer

D. C.

,

2009

,

AJ

,

138

,

219

10.1088/0004-6256/138/1/219

Crossref

Search ADS

Parsons

A. R.

,

Pober

J. C.

,

Aguirre

J. E.

,

Carilli

C. L.

,

Jacobs

D. C.

,

Moore

D. F.

,

2012

,

ApJ

,

756

,

165

10.1088/0004-637X/756/2/165

Crossref

Search ADS

Patil

A. H.

et al. ,

2017

,

ApJ

,

838

,

65

10.3847/1538-4357/aa63e7

Crossref

Search ADS

Pober

J. C.

et al. ,

2013

,

ApJ

,

768

,

L36

10.1088/2041-8205/768/2/L36

Crossref

Search ADS

Pober

J. C.

et al. ,

2014

,

ApJ

,

782

,

66

10.1088/0004-637X/782/2/66

Crossref

Search ADS

Pober

J. C.

et al. ,

2016

,

ApJ

,

819

,

8

10.3847/0004-637X/819/1/8

Crossref

Search ADS

Pritchard

J. R.

,

Loeb

A.

,

2012

,

Rep. Prog. Phys.

,

75

,

086901

10.1088/0034-4885/75/8/086901

Crossref

Search ADS

PubMed

Rahimi

M.

et al. ,

2021

,

MNRAS

,

508

,

5954

10.1093/mnras/stab2918

Crossref

Search ADS

Santos

M. G.

,

Cooray

A.

,

Knox

L.

,

2005

,

ApJ

,

625

,

575

10.1086/429857

Crossref

Search ADS

Shaver

P. A.

,

Windhorst

R. A.

,

Madau

P.

,

de Bruyn

A. G.

,

1999

,

A&A

,

345

,

380

Sievers

J. L.

,

2017

,

preprint (arXiv:1701.01860)

Thyagarajan

N.

et al. ,

2013

,

ApJ

,

776

,

6

10.1088/0004-637X/776/1/6

Crossref

Search ADS

Tingay

S. J.

et al. ,

2013

,

PASA

,

30

,

e007

10.1017/pasa.2012.007

Crossref

Search ADS

Trott

C. M.

et al. ,

2020

,

MNRAS

,

493

,

4711

10.1093/mnras/staa414

Crossref

Search ADS

van Haarlem

M. P.

et al. ,

2013

,

A&A

,

556

,

A2

10.1051/0004-6361/201220873

Crossref

Search ADS

Vedantham

H.

,

Shankar

N. U.

,

Subrahmanyan

R.

,

2012

,

ApJ

,

745

,

176

10.1088/0004-637X/745/2/176

Crossref

Search ADS

Wang

X.

,

Tegmark

M.

,

Santos

M. G.

,

Knox

L.

,

2006

,

ApJ

,

650

,

529

10.1086/506597

Crossref

Search ADS

Wayth

R. B.

et al. ,

2018

,

PASA

,

35

,

33

10.1017/pasa.2018.37

Crossref

Search ADS

Wieringa

M. H.

,

1992

,

Exp. Astron.

,

2

,

203

10.1007/BF00420576

Crossref

Search ADS

Woodbury

M. A.

,

1950

,

Memorandum Rep.

,

42

,

336

Zhang

Z.

et al. ,

2020

,

PASA

,

37

,

e045

10.1017/pasa.2020.37

Crossref

Search ADS

Zheng

H.

et al. ,

2014

,

MNRAS

,

445

,

1084

10.1093/mnras/stu1773

Crossref

Search ADS

APPENDIX A: DELAY SPECTRUM

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

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

where V_XX and V_YY, respectively, are visibilities from XX and YY linear polarization. We then use V_I (Stokes I polarization 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 equation (25) over the flat-sky model, we have

$$\begin{eqnarray*} {\boldsymbol {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)}, \end{eqnarray*}$$

(A1)

where (l, m) are direction cosines of a unit vector |$\hat{{\boldsymbol {r}}}$|⁠, 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 equation (A1) in terms of geometric time delay τ_g for a single baseline b as

$$\begin{eqnarray*} {\boldsymbol {d}}_{{\boldsymbol {b}}}(f)=\int \textrm {d}l\, \textrm {d}m\, I(l,m,f)A(l,m,f)e^{-2i\pi f\tau _{\rm g}}, \end{eqnarray*}$$

(A2)

where

$$\begin{eqnarray*} \tau _{\mathrm{g}}=\frac{1}{c}{\boldsymbol {b}}\cdot \hat{{\boldsymbol {r}}}=\frac{1}{c}(b_{x}l+b_{y}m), \end{eqnarray*}$$

(A3)

with b = (b_x, b_y), |$(u,v)=f\mathbf {b}/c$| and |$\hat{{\boldsymbol {r}}}=(l,m)$|⁠. Here, b_x and b_y are projected baseline length in metre along east and north directions in the antenna plane, respectively. Note that τ_g in equation (A2) 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 equation (A2) and then Fourier transforming it along its frequency axis, the delay transformed visibility takes a form

$$\begin{eqnarray*} \tilde{V}_{\mathrm{b}}(\tau) &=& \int _{B} \textrm {d}f\, W(f)v^{\mathrm{mea}}_{\rm b}(f)e^{2\pi if\tau }\nonumber \\ &=& \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] \nonumber \\ &&\times \, e^{2\pi if\tau } \nonumber \\ &=& \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}), \nonumber \\ \end{eqnarray*}$$

(A4)

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 τ 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, and * stands for convolution. Hence, from equation (A4), one can clearly see that there is mapping between point sources on the sky in |$\tilde{V}_{\mathrm{b}}$| space and the Dirac-delta function δ_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₂₁(k_∥, k_⊥) of 21-cm signal from the EoR fluctuations to the delay-transformed visibility |$\tilde{V}$| has been derived in Parsons et al. (2012) as

$$\begin{eqnarray*} 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}, \end{eqnarray*}$$

(A5)

where λ is the observing wavelength, k_B is the Boltzmann constant, B is the observing bandwidth, Ω is the solid angle of the primary beam of the antenna, and X and Y are cosmological scalars that 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 timestamp 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

$$\begin{eqnarray*} \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}. \end{eqnarray*}$$

(A6)

Here, Δt ≈ 42.9 s (a time resolution of PAPER-64 data) and θ_zen(Δt) is the zenith re-phasing factor. The line-of-sight k_∥, and transverse co-moving k_⊥ components of k in terms of the instrumental and cosmological parameters have been shown in Morales & Hewitt (2004),

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

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

where |$|\mathbf {u}|=\sqrt{u^{2}+v^{2}}$|⁠, η 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₁₀ = 1420 MHz is the rest frequency of 21-cm emission from cosmic re-ionization, c is the speed of light in free space, D_m(z) is the transverse co-moving distance at the reference redshift z, the present-day Hubble constant |$H_{0}=70\, {\rm km}\, {\rm s}^{-1}\, {\rm Mpc}^{-1}$|⁠, and the function |$E(z)=\sqrt{\Omega _{\mathrm{m}}(1+z)^{3}+\Omega _{k}(1+z)^{2}+\Omega _{\Lambda }}$|⁠. Here, the cosmological density parameters Ω_m, Ω_Λ, and Ω_k are total matter density, the dark energy density, and curvature, respectively. The geometric delay τ_g in equation (A3) 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{{\boldsymbol {r}}}$| is set to be zero. Such an alignment to both vectors is achieved, the maximum delay is commonly referred to ‘horizon limit (τ_H)’. Theoretically, over the (k_⊥, k_∥) cylindrical space, it can be determined from

$$\begin{eqnarray*} k_{\parallel }=\frac{H_{0}E(z)D_{\rm m}(z)}{c(1+z)}k_{\perp }. \end{eqnarray*}$$

(A7)

As 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 −τ_H ≤ τ ≤ τ_H. However, intrinsic sky components and instrumental effects that are generally non-smooth can scatter beyond the horizon limit. This is due to the sidelobes of the convolving kernel of sky and beam, i.e. |$\tilde{I}*\tilde{A}$|⁠, which broadens the footprint of k_∥ 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_∥ 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_⊥ regions (Datta, Bowman & Carilli 2010; Parsons et al. 2012; Vedantham 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_sim⁸). 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), and 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. B1 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.

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

Figure B1.

Top: (A) Expected delay power spectra, (B) recovered power spectra after CorrCal calibration, and (C) the power spectra ratio of B to A. The solid black dashed line is the horizon limit, τ_H. Bottom: Averaged 1D power spectra as a function of k_∥ for a portion of baselines in the |$\sim 30~\rm m$| group.

Open in new tab Download slide

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
December 2021	24
January 2022	30
February 2022	7
March 2022	12
April 2022	8
May 2022	10
June 2022	3
July 2022	5
August 2022	13
September 2022	4
October 2022	8
November 2022	3
December 2022	2
January 2023	2
February 2023	2
March 2023	6
April 2023	2
May 2023	2
June 2023	3
October 2023	5
November 2023	2
December 2023	1
January 2024	10
February 2024	15
March 2024	5
April 2024	5
May 2024	12
June 2024	6
July 2024	8
August 2024	20
September 2024	11
October 2024	3
November 2024	10
December 2024	6
January 2025	2
February 2025	3
March 2025	3
April 2025	6
May 2025	4

Article Contents

The correlation calibration of PAPER-64 data

ABSTRACT

1 INTRODUCTION

2 A REVIEW OF CALIBRATION FORMALISM

2.1 Sky-based calibration

2.2 Redundant calibration

3 CORRELATION CALIBRATION (CORRCAL)

3.1 Determining the sky and noise variance–covariance matrices

3.1.1 Diffuse sky component matrix (⁠|${\bf \textsf {R}}$|⁠)

3.1.2 Point source component matrix (⁠|${\bf \textsf {S}}$|⁠)

3.1.3 Noise matrix (⁠|${\bf \textsf {N}}$|⁠)

3.2 Eigendecomposition of sky signal covariance matrix

3.3 Implementation of CorrCal

3.3.1 Antenna gain solutions as a function of frequency–bandpass solution

3.4 Comparisons of CorrCal to other methods

4 THE PAPER-64 DATA SET

5 RESULTS

5.1 Post-CorrCal spectral structure of data

5.2 Delay-space power spectra

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: DELAY SPECTRUM

APPENDIX B: DETERMINING THE EOR SIGNAL LOSS IN CorrCal ANALYSIS

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

The correlation calibration of PAPER-64 data

ABSTRACT

1 INTRODUCTION

2 A REVIEW OF CALIBRATION FORMALISM

2.1 Sky-based calibration

2.2 Redundant calibration

3 CORRELATION CALIBRATION (CORRCAL)

3.1 Determining the sky and noise variance–covariance matrices

3.1.1 Diffuse sky component matrix (⁠|${\bf \textsf {R}}$|⁠)

3.1.2 Point source component matrix (⁠|${\bf \textsf {S}}$|⁠)

3.1.3 Noise matrix (⁠|${\bf \textsf {N}}$|⁠)

3.2 Eigendecomposition of sky signal covariance matrix

3.3 Implementation of CorrCal

3.3.1 Antenna gain solutions as a function of frequency–bandpass solution

3.4 Comparisons of CorrCal to other methods

4 THE PAPER-64 DATA SET

5 RESULTS

5.1 Post-CorrCal spectral structure of data

5.2 Delay-space power spectra

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: DELAY SPECTRUM

APPENDIX B: DETERMINING THE EOR SIGNAL LOSS IN CorrCal ANALYSIS

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only