# The Tracking Tapered Gridded Estimator for the 21-cm power spectrum from MWA
drift scan observations I: Validation and preliminary results

Suman Chatterjee  0000-0001-8852-5888 Department of Physics and Astronomy,
University of the Western Cape, 7535 Bellvill, Cape Town, South Africa  [
Khandakar Md Asif Elahi  0000-0003-1206-8689 Department of Physics, Indian
Institute of Technology Kharagpur, Kharagpur - 721 302, India.  Somnath
Bharadwaj  0000-0002-2350-3669 Department of Physics, Indian Institute of
Technology Kharagpur, Kharagpur - 721 302, India.  [ Shouvik Sarkar
0009-0003-3096-7028 Centre for Strings, Gravitation and Cosmology, Department
of Physics, Indian Institute of Technology Madras, Chennai 600036, India
Samir Choudhuri  Centre for Strings, Gravitation and Cosmology, Department of
Physics, Indian Institute of Technology Madras, Chennai 600036, India  Shiv
Sethi  Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bengaluru
560080, India.  Akash Kumar Patwa  0000-0002-6216-2430 Raman Research
Institute, C. V. Raman Avenue, Sadashivanagar, Bengaluru 560080, India.

###### Abstract

Drift scan observations provide the broad sky coverage and instrumental
stability needed to measure the Epoch of Reionization (EoR) 21-cm signal. In
such observations, the telescope’s pointing center (PC) moves continuously on
the sky. The Tracking Tapered Gridded Estimator (TTGE) combines observations
from different PC to estimate $P(k_{\perp},k_{\parallel})$ the 21-cm power
spectrum, centered on a tracking center (TC) which remains fixed on the sky.
The tapering further restricts the sky response to a small angular region
around TC, thereby mitigating wide-field foregrounds. Here we consider
$154.2\,{\rm MHz}$ ($z=8.2$) Murchison Widefield Array (MWA) drift scan
observations. The periodic pattern of flagged channels, present in MWA data,
is known to introduce artefacts which pose a challenge for estimating
$P(k_{\perp},k_{\parallel})$. Here we have validated the TTGE using simulated
MWA drift scan observations which incorporate the flagged channels same as the
data. We demonstrate that the TTGE is able to recover
$P(k_{\perp},k_{\parallel{}})$ without any artefacts, and estimate $P(k)$
within $5\%$ accuracy over a large $k$-range. We also present preliminary
results for a single PC, combining 9 nights of observation $(17\,{\rm min}$
total). We find that $P(k_{\perp},k_{\parallel})$ exhibits streaks at a fixed
interval of $k_{\parallel{}}=0.29\,{\rm Mpc}^{-1}$, which matches
$\Delta\nu_{\rm per}=1.28\,{\rm MHz}$ that is the period of the flagged
channels. Since the simulations demonstrate that the TTGE is impervious to the
flagged channels, the streaks seen for the actual data are possibly caused by
some systematic that has the same period as the flagged channels. These
streaks are more than $3-4$ orders of magnitude smaller than the peak
foreground power $\mid P(k_{\perp},k_{\parallel{}})\mid\approx 10^{16}\,{\rm
mK^{2}}\,{\rm Mpc^{3}}$ at $k_{\parallel{}}=0$. The streaks are not as
pronounced at larger $k_{\parallel}$, and in some cases they do not appear to
extend across the entire $k_{\perp}$ range. The rectangular region $0.05\leq
k_{\perp}\leq 0.16\,{\rm Mpc^{-1}}$ and $0.9\leq k_{\parallel{}}\leq 4.6\,{\rm
Mpc^{-1}}$ is found to be relatively free of foreground contamination and
artefacts, and we have used this to place the $2\sigma$ upper limit
$\Delta^{2}(k)<(1.85\times 10^{4})^{2}\,{\rm mK^{2}}$ on the EoR 21-cm mean
squared brightness temperature fluctuations at $k=1\,{\rm Mpc}^{-1}$.

###### keywords:

large-scale structure of universe–first stars–cosmology:reionization–diffuse
radiation, methods: statistical, technique–interferometric

National Centre for Radio Astrophysics, Tata Institute of Fundamental
Research, Post Bag 3, Ganeshkhind, Pune - 411 007, India. Suman
<EMAIL_ADDRESS><EMAIL_ADDRESS>Somnath
<EMAIL_ADDRESS>
##  1 Introduction

The epoch of reionization (EoR) when the diffuse neutral hydrogen (H I) in the
inter-galactic medium (IGM) underwent a transition to the ionized state is one
of the least understood phases in the evolution of our Universe. The
redshifted 21-cm radiation from H I is a promising observational probe of EoR
Madau1997; Bharadwaj2005; McQuinn2006; Morales2010; Pritchard2012. Several
radio interferometers, such as the Murchison Widefield Array (MWA;
Tingay2013), LOw Frequency ARray (LOFAR; vanHarlem2013), Hydrogen Epoch of
Reionization Array (HERA; Deboer2017), Giant Metrewave Radio Telescope (GMRT;
Swarup1991; Gupta2017) and the upcoming SKA-low Mellema2013; Koopmans2015 aim
to detect the power spectrum (PS) of the EoR 21-cm brightness temperature
fluctuations. Despite substantial observational efforts, it has not been
possible to detect the EoR 21-cm PS to date, and we only have upper limits
Paciga2013; Kolopanis2019; Mertens2020; Trott2020; Pal2020; Abdurashidova2022;
Kolopanis2023. At present we have the best upper limit of
$\Delta^{2}(k)<(30.76)^{2}\,{\rm mK}^{2}$ at $k=0.192\,h\,{\rm Mpc}^{-1}$ for
$z=7.9$ from HERA Abdurashidova2022.

The main challenge is that the faint H I 21-cm signal is buried in foregrounds
which are observed to be three to four orders of magnitude brighter Ali2008;
Bernardi2009; Ghosh2012; Paciga2013; Patil2017. In recent years, there have
been significant developments to mitigate the foregrounds, relying upon the
fact that foregrounds are spectrally smooth compared to the 21-cm signal. Two
approaches are generally taken to deal with the foregrounds. In the
‘foreground removal’ approach, one tries to entirely remove the foregrounds
(e.g. Chapman2012; Mertens2018; Elahi2023b). An alternative approach,
‘foreground avoidance’, only uses the $(k_{\perp},k_{\parallel{}})$ region
outside the ‘foreground wedge’ Datta2010; Morales2012; Vedantham2012;
Trott2012; Pober2016 to estimate $P(k)$ the spherical PS of the EoR 21-cm
signal (e.g. Dillon2014; Dillon2015b; Trott2020; Pal2020; Abdurashidova2022).

The Tapered Gridded Estimator (TGE; Choudhuri2014; Choudhuri2016b) is a
visibility-based 21-cm PS estimator that has been extensively used for
measuring the 21-cm PS Pal2020; Pal2022; Elahi2023; Elahi2023b; Elahi2024, and
for characterizing the diffused Galactic foregrounds Choudhuri2017a and
magnetohydrodynamic turbulence Saha2019; Saha2021. The main attribute of TGE
is that it suppresses the sidelobe responses of the telescope to mitigate the
effects of extra-galactic point source foregrounds Ghosh2011a; Ghosh2011b.
Additionally, the TGE is computationally efficient as it deals with gridded
visibilities. It is also an unbiased estimator as it internally estimates the
noise bias from the self-correlation of the visibilities to yield an unbiased
estimate of the PS. In a recent work Chatterjee2022 (hereafter,
Chatterjee2022) have introduced the Tracking Tapered Gridded Estimator (TTGE),
which generalises the TGE for estimating the 21-cm PS from drift scan
observations. Chatterjee2022 has also validated the TTGE considering simulated
MWA drift scan observations at a single frequency, where it was demonstrated
that the estimated $C_{\ell}$ matches the input model $C^{M}_{\ell}$.

Missing frequency channels, fagged to remove Radio Frequency Interferences
(RFI) or for other reasons, pose a serious problem for visibility-based PS
estimation. The Fourier transform from frequency to delay space (e.g.
Morales2004; Parsons2009) introduces artefacts in the estimated PS, and there
has been substantial work to address this problem Parsons2009; Parsons2014;
Trott2016; Kern2021; Ewall-Wice2021; Kennedy2023. We, however, note that this
problem does not arise if we correlate the visibilities directly in the
frequency domain Bharadwaj2001a; Bharadwaj2005 and use this to estimate the
21-cm signal. The TGE and the TTGE first correlate the visibilities to
estimate $C_{\ell}(\Delta\nu)$ the multi-frequency angular power spectrum
(MAPS; Datta2007; Mondal2018), and then Fourier transforms
$C_{\ell}(\Delta\nu)$ along the $\Delta\nu$ to estimate
$P(k_{\perp},k_{\parallel{}})$ the cylindrical PS. Using simulations,
Bharadwaj2018 have shown that $C_{\ell}(\Delta\nu)$ does not exhibit any
missing $\Delta\nu$ values even when $80\%$ of randomly chosen frequency
channels are flagged in the visibility data, and it is possible to estimate
$P(k_{\perp},k_{\parallel{}})$ without any artefacts due to the missing
channels. This has been further borne out in the analysis of actual data using
the TGE Pal2020; Pal2022; Elahi2023; Elahi2023b; Elahi2024.

MWA has a periodic pattern of flagged channels in the visibility data, which
introduce horizontal streaks in $P(k_{\perp},k_{\parallel{}})$ Paul2016;
Li2019; Trott2020; Patwa2021. In this paper, we consider the MWA drift scan
observation used in Patwa2021. Here we investigate if the TTGE can overcome
this issue. For this, we first simulated the MWA drift scan observation with
exactly the same flagging as in the actual data and used this to verify if the
TTGE can faithfully recover the input model power spectrum used for the
simulations. We have subsequently applied the TTGE to actual observations and
present preliminary results here.

We have arranged the paper in the following manner. First, we describe the MWA
drift scan observation in Section 2. In Section 3, we present the TTGE and the
formalism for the power spectrum. Section 4 describes the simulations and the
validation of the TTGE. We have shown the preliminary results from actual MWA
observation in Section 5. We summarize and discuss our findings in Section 6.

##  2 MWA Drift-scan Observation

Figure 1: This shows the $408\,{\rm MHz}$ Haslam map Haslam1982 scaled to 154
MHz assuming the brightness temperature spectral index $\alpha=-2.52$
Rogers2008. The iso-contours in green, magenta, red and black show the MWA
primary beam at values 0.9, 0.5, 0.05, and 0.005, respectively, for a pointing
center at ($6.1^{\circ},-26.7^{\circ}$) which corresponds to the data analysed
here. The scan starts roughly at the location of the ‘$\star$’ on the right
(RA=$349^{\circ}$) and lasts until the ‘$\star$’ on the left
(RA=$70.3^{\circ}$). Blue filled circles mark the fields EoR
0($0^{\circ},-26.7^{\circ}$) and EoR 1($60^{\circ},-26.7^{\circ}$). The red
circle shows the position of Fornax A.

The data analysed here is from Phase II (compact configuration) of the
Murchison Widefield Array (MWA, Lonsdale2009, Wayth2018) which is a radio
interferometer with 128 tiles (antennas in rest of the paper) located in
Western Australia (latitude $-26.7^{\circ}$, longitude $116.7^{\circ}$). Each
antenna consists of 16 crossed dipoles placed in a 4x4 arrangement on a square
mesh of side $\sim 4\leavevmode\nobreak\ {\rm m}$. MWA operates in several
frequency bands ranging from 80 MHz to 300 MHz. The observed visibilities are
recorded with the time resolution of 0.5 s, and they are written out at an
interval of 2 minutes (one snapshot) in which the actual duration of
observation is 112 s.

For this work we consider a particular drift scan observation (project ID
G0031) that is described in Patwa2021. In short, the observation is at a fixed
declination (DEC) $-26.7^{\circ}$ which corresponds to the zenith and it
covers a region of the sky from right ascension (RA) $349^{\circ}$ to
$70.3^{\circ}$ which spans $81.3^{\circ}$ in RA over a time duration of 5 hr
24 min. The beginning and end of the drift scan observation are marked with
$\star$ in Figure 1. Visibilities are recorded every 2 min which leads to
$162$ different pointing centers (PCs, labelled PC=1,2,…,162), located at an
interval of $0.5^{\circ}$ along RA. This region includes two well observed MWA
fields namely EoR 0 ($0^{\circ},-26.7^{\circ}$) and EoR 1
($60^{\circ},-26.7^{\circ}$) which are also shown in Figure 1. The same drift
scan observation was carried out on 10 consecutive nights. The present
observation has been performed at a nominal frequency of $\nu_{c}=154.2\,{\rm
MHz}$ (nominal redshift $z_{c}=8.2$) with $N_{c}=768$ channels of resolution
of $\Delta\nu_{c}=40\,{\rm kHz}$ covering the observing bandwidth of $B_{\rm
bw}=30.72$ MHz. This is further divided into 24 coarse bands each containing
32 channels or $1.28\,{\rm MHz}$.

Figure 2: This shows the periodic channel flagging in the observed MWA
visibility data. The entire frequency bandwidth is divided into $24$ coarse
bands of $32$ channels or $1.28\,{\rm MHz}$ width each. The colors here shows
arbitrarily normalised visibility amplitudes, and the black vertical lines
indicate the flagged channels.

The data has been pre-processed with COTTER (Offringa2015) which flags RFI and
non-working antennas. For each coarse band, COTTER also flags four channels at
both ends and one channel at the center resulting in channels (1-4,17,29-32)
to be flagged. The produces a period pattern of flagged channels as shown in
Figure 2. We apply COTTER individually to the 0.5 s time resolution visibility
data, and then average them to 10 s time resolution. This pre-processed data
is written in CASA111https://casa.nrao.edu/ readable Measurement Sets (MS).
The MS for 10 nights are calibrated separately using the steps mentioned in
Patwa2021. We note that the first $2\,{\rm hr}$ of data is missing from the
$6$-th night, and as a consequence the nights of observations is $N_{\rm
nights}=10$ for some PCs whereas it is $N_{\rm nights}=9$ for others. Since
the observations covers the same region of the sky, we perform Local Sidereal
Time (LST) stacking (e.g. Bandura2014; Amiri2022) and obtain the equivalent
one night drift scan data. We finally have $162$ MS, each corresponding to a
different pointing direction on the sky. Each MS contains visibility data with
$11$ different time stamps each with $t_{\rm int}=N_{\rm nights}\times
10\,{\rm s}$ effective integration time.

In addition to the sky signal, each measured visibility also has a system
noise contribution. For the present data, $\sigma_{\rm N}$ the r.m.s. system
noise for the real (and also imaginary) part of the measured visibilities is
estimated using

$\displaystyle\sigma_{\rm N}=\frac{1}{\eta_{s}}\frac{2K_{\rm B}T_{\rm
sys}}{A_{\rm eff}\sqrt{\Delta\nu t_{\rm int}}}$ (1)

where the $A_{\rm eff}=21.4\,{\rm m}^{2}$ and $T_{\rm sys}/\eta_{s}=290{\rm
K}$ Patwa2021. This yields $\sigma_{\rm N}=[N_{\rm nights}]^{-0.5}\times
60\,{\rm Jy}$ for a single visibility.

Considering the baseline distribution, $99$ percent of the baselines are found
to be smaller than 500 m, i.e., $257.8\lambda$ (Figure 2 of Chatterjee2022),
and we have restricted the baseline range in our analysis to $\mid U\mid\leq
245\lambda$.

##  3 The Tracking TGE

$T(\hat{\mathbf{n}},\nu)$ the brightness temperature distribution on the sky
is decomposed into spherical harmonics $Y_{\ell}^{m}(\hat{\mathbf{n}})$ using

$T(\hat{\mathbf{n}},\nu)=\sum_{\ell,m}a_{\ell}^{m{}}(\nu)Y_{\ell}^{m}(\hat{\mathbf{n}})$
(2)

where $a_{\ell}^{m{}}(\nu)$ are the expansion coefficients. The MAPS which is
defined as

$C_{\ell}(\nu_{a},\nu_{b})=\langle
a_{\ell}^{m{}}(\nu_{a})a_{\ell}^{m{*}}(\nu_{b})\rangle\,,$ (3)

jointly characterizes the angular and frequency dependence of the two-point
statistics of $T(\hat{\mathbf{n}},\nu)$. Further, the MAPS is a function of
the frequency separation

$C_{\ell}(\nu_{a},\nu_{b})=C_{\ell}(\Delta\nu)\,\,{\rm
where}\,\,\Delta\nu=\mid\nu_{b}-\nu_{a}\mid$ (4)

if $T(\hat{\mathbf{n}},\nu)$ is assumed to be statistically homogeneous
(ergodic) along the line of sight direction.

Chatterjee2022 presents the tracking tapered gridded estimator (TTGE) to
determine the MAPS $C_{\ell}(\nu_{a},\nu_{b})$ using the measured visibility
data from drift scan radio-interferometric observations. Chatterjee2022 also
explained how $C_{\ell}(\nu_{a},\nu_{b})$ can be used to determine the
cylindrical power spectrum $P(k_{\perp},k_{\parallel})$. However, the
validation in that paper is restricted to a single frequency $\nu$. There, we
have simulated observations where the sky signal $T(\hat{\mathbf{n}})$
corresponds to an input model $C^{M}_{\ell}$. We have applied TTGE on the
simulated visibilities to estimate $C_{\ell}$ and demonstrated that the
estimated $C_{\ell}$ matches the input model. In this section of the present
paper, we briefly summarize the formalism of the TTGE, and in Section 4, we
validate it considering multi-frequency observations and the full three-
dimensional power spectrum $P(k)$.

The visibility $\mathcal{V}\left(\mathbf{U},\nu\right)$ measured at a baseline
$\mathbf{U}$ and frequency $\nu$ is given by,

$\mathcal{V}\left(\mathbf{U},\nu\right)=Q_{\nu}\int
d\Omega_{\hat{\mathbf{n}}}T\left(\hat{\mathbf{n}},\nu\right)A\left(\Delta\hat{\mathbf{n}},\nu\right)e^{2\pi
i\mathbf{U}\cdot\Delta\mathbf{n}},$ (5)

where $Q_{\nu}=2k_{B}/\lambda^{2}$ is the conversion factor from brightness
temperature to specific intensity in the Raleigh-Jeans limit,
$T\left(\hat{\mathbf{n}},\nu\right)$ is the brightness temperature
distribution on the sky, $d\Omega_{\hat{\mathbf{n}}}$ is the elemental solid
angle in the direction $\hat{\mathbf{n}}$,
$A\left(\Delta\hat{\mathbf{n}},\nu\right)$ is the antenna primary beam (PB)
pattern, and $\Delta\hat{\mathbf{n}}=\hat{\mathbf{n}}-\hat{\mathbf{p}}$ where
$\hat{\mathbf{p}}$ is the telescope’s pointing direction or pointing center
(PC). In this work, we consider drift scan observations where the telescope is
held fixed (on earth) to point towards the zenith. Here $\hat{\mathbf{p}}$ has
a fixed declination $\delta_{0}$ (latitude of the array), whereas the right-
ascension $\alpha_{p}$ varies due to the earth’s rotation. In such a
situation, it is convenient to use
$\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})$ where we have included
$\alpha_{p}$ as a parameter which tells us the pointing direction
$\hat{\mathbf{p}}$ for which the visibilities were recorded. As noted earlier,
for the observations considered here, we have $\delta_{0}=-26.7^{\circ}$ and
$349^{\circ}\leq\alpha_{p}\leq 70.3^{\circ}$ at an interval
$\Delta\alpha_{p}=0.5^{\circ}$.

The question is “How do we combine the visibility data from different
$\hat{\mathbf{p}}$?”. To address this, we choose a tracking center (TC)
$\hat{\mathbf{c}}$ which remains fixed on the sky and consider the convolved
visibilities defined in the $uv$ plane as

$\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)=\int
d^{2}U^{{}^{\prime}}\,\tilde{w}(\mathbf{U}-\mathbf{U}^{{}^{\prime}})\,e^{2\pi
i\mathbf{U}^{{}^{\prime}}\cdot\bm{\chi}_{p}}\,\mathcal{V}(\alpha_{p},\mathbf{U}^{{}^{\prime}},\nu)\,$
(6)

where $\bm{\chi}_{p}=\hat{\mathbf{p}}-\hat{\mathbf{c}}$ and
$\tilde{w}(\mathbf{U})$ is the Fourier transform of a tapering function
$W\left(\bm{\theta}\right)$ which is typically chosen to be considerably
narrower than the antenna’s primary beam pattern
$A\left(\bm{\theta},\nu\right)$. Here we have used a Gaussian
$W\left(\bm{\theta}\right)=e^{-\theta^{2}/\theta_{w}^{2}}$. Adopting the flat-
sky approximation, eq. (6) can be expressed as

$\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)=Q_{\nu}\int
d^{2}\theta\,T\left(\bm{\theta},\nu\right)\,W(\bm{\theta})\,A\left(\bm{\theta}-\bm{\chi}_{p},\nu\right)e^{2\pi
i\mathbf{U}\cdot\bm{\theta}}$ (7)

where the phase center of $\mathcal{V}_{c}(\mathbf{U},\nu)$ is shifted to
$\hat{\mathbf{c}}$, and the function $W(\bm{\theta})$ restricts the sky
response to a small region around $\hat{\mathbf{c}}$ . For the purpose of this
discussion, we may assume that the convolved visibility
$\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)$ only records the signal from a
small region of the sky centered around $\hat{\mathbf{c}}$. We can now
coherently combine the data from different pointings

$\mathcal{V}_{c}(\mathbf{U},\nu)=\sum_{p}s_{p}\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)$
(8)

to estimate the sky signal from a small region around $\hat{\mathbf{c}}$. The
contribution to $\mathcal{V}_{c}(\mathbf{U},\nu)$ from
$\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)$ at a particular PC
$\hat{\mathbf{p}}$ goes down as $\sim A\left(-\bm{\chi}_{p},\nu\right)$ (eq.
7) which declines rapidly as the separation
$\bm{\chi}_{p}=\hat{\mathbf{p}}-\hat{\mathbf{c}}$ increases. In addition to
the sky signal, each $\mathcal{V}_{c}(\alpha_{p},\mathbf{U},\nu)$ also
contains a system noise contribution. This implies that the PC
$\hat{\mathbf{p}}$, which are close to $\hat{\mathbf{c}}$ contribute to
$\mathcal{V}_{c}(\mathbf{U},\nu)$ with a relatively high signal-to-noise ratio
(SNR) as compared to those which are at a large angular distance from
$\hat{\mathbf{c}}$. We account for this by suitably choosing the factor
$s_{p}$, which assigns different weights to every PC.

Considering the observational data, we evaluate the convolved visibilities on
a rectangular grid (labelled using $g$) on a $uv$ plane. Note that our
analysis does not account for the fact that the baseline
$\mathbf{U}=\mathbf{d}\nu/c$ corresponding to a fixed antenna separation
$\mathbf{d}$ changes with frequency, and it is held fixed at the value
corresponding to $\nu_{c}$. The convolved gridded visibilities are evaluated
using

$\displaystyle\mathcal{V}_{cg}(\nu_{a})=\sum_{p}s_{p}\sum_{n}$
$\displaystyle\tilde{w}(\mathbf{U}_{g}-\mathbf{U}_{i})e^{2\pi
i\mathbf{U}_{i}\cdot\bm{\chi}_{p}}\times$
$\displaystyle\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})F_{p,n}(\nu_{a})\,$
(9)

where $F_{p,n}(\nu_{a})$ has a value $0$ if the particular visibility is
flagged and is $1$ otherwise. Here the sum over different baselines
$\mathbf{U}_{i}$ is evaluated for a fixed PC (labelled $p$), whereas the outer
sum combines all the different PCs covered in the drift scan observation. The
final $\mathcal{V}_{cg}(\nu_{a})$ refers to the sky signal centered at a
particular TC $\hat{\mathbf{c}}$.

Following eq. (6) of Pal2022, we define the tracking tapered gridded estimator
(TTGE) as

$\displaystyle\hat{E}_{g}(\nu_{a},\nu_{b})=$ $\displaystyle
M^{-1}_{g}(\nu_{a},\nu_{b}){\mathcal{R}e}\,\Big{[}\mathcal{V}_{cg}(\nu_{a})\mathcal{V}^{*}_{cg}(\nu_{b})$
$\displaystyle-\sum_{p,i}F_{p,i}(\nu_{a})F_{p,i}(\nu_{b})\mid
s_{p}\tilde{w}(\mathbf{U}_{g}-\mathbf{U}_{i})\mid^{2}$
$\displaystyle\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})\mathcal{V}^{*}(\alpha_{p},\mathbf{U}_{i},\nu_{b})\Big{]}$
(10)

where $M_{g}(\nu_{a},\nu_{b})$ is a normalisation factor and
${\mathcal{R}e}[]$ denotes the real part. The second term in the r.h.s. of eq.
(10) subtract out the contribution from the correlation between the
visibilities measured at the same baseline and pointing direction. This is
primarily introduced to cancel out the noise bias which arises when we
correlate a visibility with itself (i.e. $a=b$; Chatterjee2022).

We use all-sky simulations (Section 4.1) to estimate $M_{g}(\nu_{a},\nu_{b})$.
The simulated sky signal $T(\hat{\mathbf{n}},\nu)$ is a realisation of a
Gaussian random field corresponding to an unit multi-frequency angular power
spectrum (UMAPS) where $C_{\ell}(\nu_{a},\nu_{b})=1$. The simulations consider
identical drift scan observations as the actual data, with exactly the same
flagging, and frequency and baseline coverage. The simulated visibilities
$[\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})]_{\rm UMAPS}$ are analysed
exactly the same as the actual data. We have estimated
$M_{g}(\nu_{a},\nu_{b})$ using

$\displaystyle
M_{g}(\nu_{a},\nu_{b})=\Big{\langle}{\mathcal{R}e}\,\Big{[}\mathcal{V}_{cg}(\nu_{a})\mathcal{V}^{*}_{cg}(\nu_{b})-\sum_{p,n}F_{p,n}(\nu_{a})F_{p,n}(\nu_{b})\times$
$\displaystyle\mid
s_{p}\tilde{w}(\mathbf{U}_{g}-\mathbf{U}_{i})\mid^{2}\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})\mathcal{V}^{*}(\alpha_{p},\mathbf{U}_{i},\nu_{b})\Big{]}\Big{\rangle}_{\rm
UMAPS}$ (11)

where the angular brackets $\langle...\rangle$ denote an average over multiple
realisations of the UMAPS. Here we have used $100$ random realisations of the
UMAPS to reduce the statistical uncertainties in the estimated
$M_{g}(\nu_{a},\nu_{b})$.

The MAPS TTGE defined in eq. (10) provides an unbiased estimate of
$C_{\ell}{{}_{g}}(\nu_{a},\nu_{b})$ at the angular multipole $\ell_{g}=2\pi
U_{g}$ i.e.

$\langle{\hat{E}}_{g}(\nu_{a},\nu_{b})\rangle=C_{\ell}{{}_{g}}(\nu_{a},\nu_{b})$
(12)

To enhance the SNR and also reduce the data volume, we have divided the $uv$
plane into annular bins. We use this to define the binned TTGE

$\hat{E}_{G}[q](\nu_{a},\nu_{b})=\frac{\sum_{g}w_{g}{\hat{E}}_{g}(\nu_{a},\nu_{b})}{\sum_{g}w_{g}}\,.$
(13)

where $w_{g}$ refers to the weight assigned to the contribution from any
particular grid point $g$. For the analysis presented in this paper we have
used the weight $w_{g}=M_{g}(\nu_{a},\nu_{b})$ which roughly averages the
visibility correlation
$\mathcal{V}_{cg}(\nu_{a})\,\mathcal{V}_{cg}^{*}(\nu_{b})$ across all the grid
points which are sampled by the baseline distribution and lie within the
boundaries of bin $q$.

The binned estimator has an expectation value

$\bar{C}_{\bar{\ell}_{q}}(\nu_{a},\nu_{b})=\frac{\sum_{g}w_{g}C_{\ell}{{}_{g}}(\nu_{a},\nu_{b})}{\sum_{g}w_{g}}$
(14)

where $\bar{C}_{\bar{\ell}_{q}}(\nu_{a},\nu_{b})$ is the bin averaged MAPS at

$\bar{\ell}_{q}=\frac{\sum_{g}w_{g}\ell_{g}}{\sum_{g}w_{g}}$ (15)

which is the effective angular multipole for bin $q$. For brevity of notation,
we use $C_{\ell}(\nu_{a},\nu_{b})$ instead of
$\bar{C}_{\bar{\ell}_{q}}(\nu_{a},\nu_{b})$ in the subsequent discussion.

For the subsequent analysis, we assume that the 21-cm signal is statistically
homogeneous (ergodic) along the line of sight whereby
$C_{\ell}(\nu_{a},\nu_{b})=C_{\ell}(\Delta\nu)$ (eq. 4). Such an assumption is
valid for the redshifted 21-cm signal if the observing bandwidth $B_{w}$ is
sufficiently small ($\approx 8\,{\rm MHz}$; Mondal2018) such that the mean
neutral hydrogen fraction does not evolve significantly across the
corresponding redshift interval. Although $B_{w}=30.72\,{\rm MHz}$ used here
is quite a bit larger and we may expect some signal loss at the small $k$, we
still adopt this assumption as it significantly simplifies the analysis as we
can quantify the signal using $P(k_{\perp},\,k_{\parallel})$ the cylindrical
PS of the 21-cm brightness temperature fluctuations which is related to a
$C_{\ell}(\Delta\nu)$ through a Fourier transform. We have Datta2007

$P(k_{\perp},\,k_{\parallel})=r^{2}\,r^{\prime}\int_{-\infty}^{\infty}d(\Delta\nu)\,e^{-ik_{\parallel}r^{\prime}\Delta\nu}\,C_{\ell}(\Delta\nu)$
(16)

where $k_{\parallel}$ (the Fourier conjugate of $\Delta\nu$) and and
$k_{\perp}=\ell/r$ are component of $\mathbf{k}$ respectively parallel and
perpendicular to the line of sight. Here $r=9210\,{\rm Mpc}$ is the comoving
distance corresponding to $\nu_{c}$ and $r^{\prime}=dr/d\nu=16.99\,{\rm
Mpc\,MHz}^{-1}$ evaluated at $\nu_{c}$. The reader is referred to
Bharadwaj2018 for further details.

To proceed further with the TTGE, we have collapsed the measured
$C_{\ell}(\nu_{a},\nu_{b})$ to obtain $C(\Delta\nu_{n})$ where
$\Delta\nu_{n}=n\,\Delta\nu_{c}$ with $n=0,1,2,...,N_{c}-1$. We then use a
maximum likelihood estimator (MLE; Pal2022) to obtain the cylindrical PS

$P(k_{\perp},k_{\parallel
q})=\sum_{n}\left[\left(\textbf{A}^{\dagger}\textbf{N}^{-1}\textbf{A}\right)^{-1}\textbf{A}^{\dagger}\textbf{N}^{-1}\right]_{qn}\mathcal{W}(\Delta\nu_{n})\,C_{\ell}(\Delta\nu_{n})$
(17)

where the Hermitian matrix A contains the Fourier coefficients and N is the
noise covariance matrix. We have assumed N to be diagonal, and estimated it
using noise only simulations. In these simulations the real and imaginary
parts of the measured visibilites were both replaced with Gaussian random
noise of variance $\sigma_{\rm N}^{2}$. These simulated data were run through
the TTGE pipeline to estimate $C_{\ell}(\Delta\nu_{n})$, and we have used $20$
realisations of the simulated visibilities to estimate $[\delta
C_{\ell}(\Delta\nu_{n})]^{2}$ the variance of $C_{\ell}(\Delta\nu_{n})$ which
are the diagonal elements of N. Note that TTGE avoids the noise bias, and the
mean $C_{\ell}(\Delta\nu_{n})$ is expected to be zero for noise only
simulations.

The window function $\mathcal{W}(\Delta\nu_{n})$, which is normalized to unity
at $\Delta\nu=0$, is used to avoid a discontinuity in the measured
$C_{\ell}(\Delta\nu_{n})$ at the band edges. For the present work we have used
a Blackman-Nuttall (BN; Nuttall1981) window function.

##  4 Validating the TTGE

###  4.1 All-sky simulation

The aim here is to simulate the visibilities that would be measured in MWA
drift scan observations of $T(\hat{\mathbf{n}},\nu)$ the redshifted 21-cm
brightness temperature distribution on the sky. We simulate
$T(\hat{\mathbf{n}},\nu)$ using the package HEALPix (Gorski2005, Hierarchical
Equal Area isoLatitude Pixelization of a sphere;). For the simulations
presented here, we have set $N_{\rm side}=512$ in HEALPix, which results in
$\ell_{max}=1535$ and a pixel size of $6.87^{{}^{\prime}}$.

We assume $T(\hat{\mathbf{n}},\nu)$ to be a Gaussian random field (GRF) with a
given input model power spectrum $P^{m}(k_{\perp},k_{\parallel{}})$. To
simulate the signal, we consider

$C_{\ell}(\Delta\nu_{n})=\frac{1}{N_{c}\Delta\nu_{c}\,r^{2}\,r^{\prime}}\sum_{q=0}^{N_{c}-1}P(k_{\perp},k_{\parallel
q})\,e^{2\pi inq/N_{c}}\,$ (18)

which is the discrete representation of the inverse of eq. (16). We generate
the expansion coefficients $a_{\ell}^{m{}}(\nu_{n})$ at the $n$-th frequency
channel (eq. 3) using

$a_{\ell}^{m{}}(\nu_{n})=\sum_{q=0}^{N_{c}-1}\left[\sqrt{\frac{P^{m}(k_{\perp},k_{\parallel{q}})}{N_{c}\Delta\nu_{c}\,r^{2}\,r^{\prime}}}\left(\frac{\hat{\mathbf{x}}_{q}+i\hat{\mathbf{y}}_{q}}{\sqrt{2}}\right)\right]\,e^{\frac{2\pi
inq}{N_{c}}}\,$ (19)

where $\hat{\mathbf{x}}_{q},\hat{\mathbf{y}}_{q}$ are independent Gaussian
random variables of unit variance, and use eq. (2) to calculate
$T(\hat{\mathbf{n}},\nu_{n})$. We now use eq. (4) of Chatterjee2022 to
calculate the simulated visibilities
$\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})$ for all the PCs covered in
the actual data. Note that we have simulated the visibilities for the same
baseline distribution, frequency channels, and flagging of the actual data.
For the present work we have assumed the model PS to be

$P^{m}(k)=(k/k_{0})^{s}\,{\rm K}^{2}\,{\rm Mpc}^{3}\,$ (20)

with $k_{0}=1\,{\rm Mpc}^{-1}$ and $s=-1$. We have used $20$ independent
relaisations of the simulated signal to estimate the mean and $1\sigma$ errors
for all the results presented below.

In order to simulate the sky signal corresponding to UMAPS, we first consider
a fixed frequency channel (say $n=0$) and generate a GRF
$T(\hat{\mathbf{n}},\nu_{0})$ corresponding to a unit angular power spectrum
$C_{\ell}=1$. We then assign the same brightness temperature map to all the
other frequency channels within our observing bandwidth i.e.
$T(\hat{\mathbf{n}},\nu_{n})=T(\hat{\mathbf{n}},\nu_{0})$ for
$n=1,2,3,...,N_{c}-1$. This ensures that the simulated sky signal
$T(\hat{\mathbf{n}},\nu_{n})$ corresponds to $C_{\ell}(\nu_{a},\nu_{b})=1$
(UMAPS). The UMAPS visibilities
$[\mathcal{V}(\alpha_{p},\mathbf{U}_{i},\nu_{a})]_{\rm UMAPS}$ were simulated
using HEALPix exactly the same way as for the model PS, with the exception
that the sky signal is different. As mentioned earlier, we have used $100$
realisations of UMAPS to estimate $M_{g}(\nu_{a},\nu_{b})$.

We have introduced noise in the simulations as a GRF with zero mean and
standard deviation $\sigma_{\rm N}=10\,{\rm Jy}$. This is equivalent to an
observation where $N_{\rm nights}=36$ (eq. 1) and sets the ${\rm SNR}=3$ at
the smallest $k$-bin for the input model considered here. The actual EoR
signal is significantly fainter and requires much longer observations for a
detection.

###  4.2 Single pointing

Figure 3: This shows $C_{\ell}(\Delta\nu)$ as a function of $\Delta\nu$ for
four values of $\ell$. The data points with $1\sigma$ error bars are estimated
from 20 realizations of the all-sky simulations. The lines show the
theoretical predictions calculated using the input model power spectrum
$P^{m}(k)=(1{\rm Mpc}^{-1}/k)\,{\rm K^{2}}\,{\rm Mpc^{3}}$ in equation (16).
The $\Delta\nu=0$ points have been slightly shifted for the convenience of
plotting on a logarithmic scale.

In this subsection, and also the subsequent subsection, we use the simulated
visibilities to validate the TTGE considering a single tracking center (TC) at
$({\rm RA,DEC)}=(6.1^{\circ},-26.7^{\circ})$. We have used a window function
(eq. 7) $W\left(\bm{\theta}\right)=e^{-\theta^{2}/\theta_{w}^{2}}$ where
$\theta_{w}=0.6\,\theta_{\rm FWHM}$ with $\theta_{\rm FWHM}=15^{\circ}$ . This
effectively restricts the sky signal to an angular region of extent $\sim
15^{\circ}$ centered around the TC, and this remians fixed even as the
pointing center PC drifts across the sky. In this section, we have considered
a single PC (=34) which exactly coincides with the TC. Note that the
simulations used here do not contain any system noise.

Figure 3 shows the estimated MAPS $C_{\ell}(\Delta\nu)$ along with the model
predictions $C_{\ell}^{\rm M}(\Delta\nu)$. Note that the range $50<\ell<1493$
has been divided into 20 $\ell$ bins, and the results are shown for $4$
representative bins. We see that model predictions are within the $1\sigma$
error bars of the estimated $C_{\ell}(\Delta\nu)$, indicating that the two are
in good agreement. It is important to note that the estimated
$C_{\ell}(\Delta\nu)$ shows a smooth $\Delta\nu$ dependence with no missing
frequency separations $\Delta\nu$. In particular, we do not see any artefacts
in the estimated $C_{\ell}(\Delta\nu)$ due to the periodic pattern of flagged
channels present in the simulations (Figure 2).

Figure 4: Left panel shows the cylindrical power spectrum
$P(k_{\perp},k_{\parallel{}})$ estimated from simulations with MWA coarse
channel flagging. For comparison, the right panel shows the
$P(k_{\perp},k_{\parallel{}})$ estimated from simulations without coarse
channel flagging. We do not notice any significant difference.

The left panel of Figure 4 shows the cylindrical PS
$P(k_{\perp},k_{\parallel{}})$ estimated from the $C_{\ell}(\Delta\nu)$ shown
in Figure 3 using eq. (17). For comparison, the right panel shows
$P(k_{\perp},k_{\parallel{}})$ estimated using identical simulations which do
not incorporate the channel flagging. We note that the
$P(k_{\perp},k_{\parallel{}})$ shown in the two panels are visually
indistinguishable. Several earlier works which have first transformed from
frequency to delay space and then estimated the PS (e.g. Paul2016; Li2019;
Trott2020; Patwa2021) have reported horizontal streaks in the estimated
$P(k_{\perp},k_{\parallel{}})$ due to the periodic pattern of flagged channels
present in the MWA data. In our approach Bharadwaj2018, we first estimate
$C_{\ell}(\Delta\nu)$, which does not have any missing $\Delta\nu$ even when
the visibility data contains flagged channels. We then Fourier transform
$C_{\ell}(\Delta\nu)$ to obtain a clean estimate of
$P(k_{\perp},k_{\parallel{}})$. We see that the missing frequency channels do
not introduce any artefacts in the estimated $P(k_{\perp},k_{\parallel{}})$.

Figure 5:  The upper panel shows the estimated spherically binned power
spectrum $P(k)$ and $1-\sigma$ error bars for simulations for PC=34 with no
noise and coarse channel flagging. For comparison, the input model
$P^{m}(k)=(1{\rm Mpc}^{-1}/k)\,{\rm K^{2}}\,{\rm Mpc^{3}}$ is also shown by
the solid line. The lower panels show the percentage error
$\Delta=[P(k)-P^{m}(k)]/P^{m}(k)$ (data points) and the relative statistical
fluctuation $\sigma/P^{m}(k)\times 100\%$ (between the solid lines). The four
lower panels consider situations for combining different PCs mentioned in the
figure legends.

The top panel of figure 5 shows $P(k)$ the spherical PS estimated directly
from the $C_{\ell}(\Delta\nu)$ shown in Figure 3 using a maximum likelihood
estimator (MLE) Elahi2023, for comparison we also show $P^{m}(k)$ the input
model PS. We see that the model predictions $P^{m}(k)$ are within the
$1\sigma$ error bars of the estimated $P(k)$, indicating that the two are in
good agreement through the entire $k$ range 0.01 to 4 ${\rm Mpc}^{-1}$ probed
here. We have quantified the percentage deviation between $P(k)$ and
$P^{m}(k)$ using $\Delta=[P(k)-P^{m}(k)]/P^{m}(k)\times 100$ shown in the
second panel (from top) of figure 5. We see that for nearly all $k$ the values
of $\Delta$ are within the predicted $1\sigma$ error bars. The values of
$\mid\Delta\mid$ are within $1\%$ for $k\geq 0.13\,{\rm Mpc^{-1}}$ and within
$2\%$ for $k\geq 0.035\,{\rm Mpc^{-1}}$. We see that $\mid\Delta\mid$
increases at smaller $k$ due to the convolution with the window function and
the primary beam pattern. We have the largest deviation $|\Delta|=6.9\%$ at
$k=0.013\,{\rm Mpc^{-1}}$. Overall, we find very good agreement between $P(k)$
and $P^{m}(k)$. Further, we do not find any artefacts due to the periodic
pattern of flagged channels present in the MWA data.

Figure 6: The upper panel shows a comparison of SNR achievable for a single
pointing with (circles) and without (squares) the periodic pattern of flagged
channels. The triangles show the expected SNR values in the presence of system
noise with $\sigma_{\rm N}=10\,{\rm Jy}$ (eq. 1). The lower panel shows the
ratio of the SNR values without and with flagging. The shaded region indicates
the $k$-range where the SNR values remain mostly unaffected due to flagging.

The flagged channels present in the MWA data cause $\sim 28\%$ data loss, and
we expect this to degrade the SNR of the 21-cm PS relative to the situation
when there are no flagged channels. The upper panel of Figure 6 shows a
comparison of the SNR for the 21-cm PS estimated from the simulations with
(squares) and without (circles) the periodic pattern of flagged channels. In
the absence of flagging, the SNR has value $\sim 4$ at the smallest $k$ bin
and rises monotonically to $\sim 300$ at the largest $k$ bin. The results with
flagging are similar, but the SNR values are somewhat smaller. The lower panel
shows the ratio of the SNR without flagging to with flagging. We find that the
ratio is very close to unity at $k<0.09\,{\rm Mpc}^{-1}$, and the ratio
increases only at large $k$ where it varies in the range $1$ to $2$. We find
that The ratio has a maximum value $(\sim 2)$ at $k\sim 0.5\,{\rm Mpc}^{-1}$.
Note that the discussion, till now, has not considered the system noise. The
upper panel of Figure 6 also shows the SNR for the simulations which include
the system noise. We see that the SNR at the two smallest $k$ bins are not
much affected even if we introduce the system noise. The total noise budget in
these two bins is dominated by the cosmic variance (CV), and the system noise
make only a small contribution. We find that there is a very substantial drop
in the SNR at the larger $k$ bins when we introduce the system noise. There is
very little difference in the SNR between with and without flagging once we
introduce the system noise contribution.

###  4.3 Multiple pointings

Figure 7: This shows a measure of visibility correlation between two different
PCs for a fixed TC. The correlation = $C_{\ell}(\Delta\nu=0,\Delta{\rm
PC})/C_{\ell}(\Delta\nu=0,\Delta{\rm PC}=0)$. Here we have fixed the TC at
PC=34 and estimated the MAPS $C_{\ell}(\Delta\nu)$ by correlating visibilities
separated by $\Delta{\rm PC}$ pointing centers. Here we only show the results
for $\Delta\nu=0$, however, $\Delta\nu>0$ results are quite similar.

In this subsection we coherently combine the visibilities measured at multiple
PCs to estimate the signal in a small angular region centered at the fixed TC
$\hat{\mathbf{c}}$ whose sky coordinates we have mentioned in the previous
subsection. As mentioned earlier, the contribution from a PC declines as $\sim
A\left(-\bm{\chi}_{p},\nu\right)$ (eq. 7) as the separation
$\bm{\chi}_{p}=\hat{\mathbf{p}}-\hat{\mathbf{c}}$ increases. The PC which are
close to the TC contribute to $\mathcal{V}_{cg}(\nu_{a})$ with the higher SNR
as compared to the PC which are at a large angular separation, and we account
for this by suitably choosing the factor $s_{p}$ (eq. 9).

In the previous subsection we have considered a single pointing center PC=34
which exactly coincides with the TC i.e. $\hat{\mathbf{p}}=\hat{\mathbf{c}}$
and $\bm{\chi}_{p}=0$. In this subsection we consider other pointing
directions PC=34 + $\Delta$ PC. Before combining the signal from multiple
pointing directions, we analyse how the correlation between the signal from
two different pointing directions PC=34 and PC=34 + $\Delta$ PC changes as
$\Delta\,{\rm PC}$ is varied. To evaluate this we consider the dimensionless
correlation coefficient

$c_{\ell}(\Delta\,{\rm PC})=[C_{\ell}(0)]_{\Delta\,{\rm PC}}/C_{\ell}(0)$ (21)

where we have evaluated $[C_{\ell}(0)]_{\Delta\,{\rm PC}}$ using eq (10) with
the difference that $\mathcal{V}_{cg}(\nu_{a})$ and
$\mathcal{V}^{*}_{cg}(\nu_{b})$ refer to PC=34 and PC=34 + $\Delta$ PC
respectively. Note that these simulations do not contain system noise, and it
is not necessary to subtract out the self-correlation in eq (10).

Figure 7 shows $c_{\ell}(\Delta\,{\rm PC})$ as a function of $\ell$ for
different values of $\Delta\,{\rm PC}$ in the range $1$ to $5$. We see that in
all cases, the correlation drops as $\Delta\,{\rm PC}$, the offset between TC
and PC, is increased. This is roughly consistent with the expected $\sim
A\left(-\bm{\chi}_{p},\nu\right)$ (eq. 7) decline, and also the finding of
Patwa2019. This analysis sets the choice for combining different PC later in
the analysis.

Considering $\Delta\,{\rm PC}=1$, we see that $c_{\ell}(\Delta\,{\rm
PC})\approx 1$ for $\ell\leq 300$, it is $>0.9$ for $\ell\leq 800$ and it
drops at larger $\ell$ to $\sim 0.8$ at $\ell=2,000$. We find a similar
behaviour for larger $\Delta\,{\rm PC}$, but the values of
$c_{\ell}(\Delta\,{\rm PC})$ are smaller. Overall, the signal at small
$\ell\,(\leq 300)$ remains correlated for large $\Delta\,{\rm PC}$, even
beyond $5$. However, at larger $\ell$ $(>800)$ the correlation falls below
$\sim 0.7$ for $\Delta\,{\rm PC}\geq 3$. We have also investigated
$c_{\ell}(\Delta\,{\rm PC})$ for other values of $\Delta\nu$ (eq. 21), and we
find that the behaviour is very similar to that for $\Delta\nu=0$ shown here.
We note that the variation of $c_{\ell}(\Delta\,{\rm PC})$ with $\Delta\,{\rm
PC}$ depends on the width of the window function where we have used
$\theta_{\rm FWHM}=15^{\circ}$, and we expect the correlation to decrease
faster with $\Delta\,{\rm PC}$ if $\theta_{\rm FWHM}$ is reduced.

The results shown in Figure 7 indicate that it is most optimal to consider the
weights $s_{p}$ in eq. (9) as functions of $\ell$. However, we have not
considered this here. We have considered three simple schemes where we use
uniform wights $s_{p}=1$ to combine multiple pointings. In the first scheme we
combine PC=33, 34 and 35 (33-35) for which the maximum $\mid\Delta\,{\rm
PC}\mid$ has value 1. We similarly also consider PC= 32-36 and 31-37 for which
the maximum $\mid\Delta\,{\rm PC}\mid$ has values $2$ and $3$ respectively.
The correlation at large $\ell$ falls considerably for $\mid\Delta\,{\rm
PC}\mid\geq 4$, and we have not considered combining such pointings.
Considering PC=31-37, the RA difference between PC=31 and PC=37 is $6^{\circ}$
and the total time duration is $14\,{\rm min}$.

The $C_{\ell}(\Delta\nu)$, $P(k_{\perp},k_{\parallel})$ and $P(k)$ estimated
after combining multiple pointings are very similar to those for a single
pointing, and we have not shown these separately. Considering the three
schemes for combining multiple pointings, the three lower panels of Figure 5
show $\Delta$ the percentage deviation between the estimated $P(k)$ and the
input model $P^{m}(k)$. We see that the results from the three schemes are
practically indistinguishable from those for a single pointing. Note that this
is expected as the simulations considered here do not have any system noise.
We expect the system noise to drop when multiple pointings are combined,
however, we do not expect the sky signal or the cosmic variance to change.

Figure 8:  Left upper panel shows the estimated spherical PS $P(k)$ after
combining multiple pointing centers coherently. The solid line shows the input
model ($P^{m}(k)=(1{\rm Mpc}^{-1}/k)\,{\rm K^{2}}\,{\rm Mpc^{3}}$) used for
the simulations. 20 realizations of the simulations are used to estimate the
mean and $1\sigma$ errors shown here. The red triangles and the error-bars
show the results for PC=34, whereas the black circles and error-bars
correspond to PC=31-37 i.e. all the pointing centers between PC=31 and 37 have
been combined. Considering the $P(k)$ estimates in the left upper panel, the
left lower panel shows the corresponding percentage deviation relative to the
input model, the hatched regions (in red) and the region between black solid
line show the $1\sigma$ statistical fluctuations for PC=34 and PC=31-37
respectively. Blue dotted lines show $|\Delta|=10\%$. The right panel shows
the SNR after combining multiple PCs. Here, PC=31-37(CV) refers to the SNR in
the absence of system noise, where we only have cosmic variance (CV). This is
very close to the SNR for a single pointing (PC=34 - Flagging) shown in the
right panel of Figure 6.

Figure 8 shows the estimated $P(k)$ and the SNR when the simulations include
system noise. Considering the left panel, we see that the estimated $P(k)$ is
consistent with the input model $P^{m}(k)$ within the $1\sigma$ error bars
through the entire $k$ range. However, the predicted $1\sigma$ error bars and
the percentage deviation from $P^{m}(k)$ depend on the number of PCs that we
combine. For a single PC we find that $\mid\Delta\mid<10\%$ for $k<0.2\,{\rm
Mpc}^{-1}$ whereas it increases up to $50\%$ at larger $k$. Considering
multiple pointings where we combine PC=31-37, we find that the $1\sigma$ error
bars reduce considerably, and $\mid\Delta\mid<10\%$ for the entire $k$ range.
The right panel of Figure 8 shows how the SNR improves as we combine multiple
PCs. Note that this is very similar to the right panel of Figure 6 which shows
the SNR for a single pointing PC=34. For reference, PC=31-37(CV) shows the SNR
in the absence of system noise, where we only have cosmic variance. We do not
expect the cosmic variance to change when we combine multiple pointings, and
the SNR for PC=31-37(CV) is very close to that for ‘PC=34-Flagging’ shown in
the right panel of Figure 6. Considering the results with system noise, we
note that the first $2$ $k$ bins are cosmic variance (CV) dominated, and these
SNR values do not change when we introduce system noise or combine PCs. The
SNR in all the other $k$ bins degrades substantially when we introduce system
noise. We see that we have the lowest SNR when we consider a single PC as
compared to multiple PCs. For a single PC, the SNR peaks in the second and
third $k$ bins where ${\rm SNR}=6$, the SNR falls at larger $k$ and it is
below unity at $k>0.09\,{\rm Mpc}^{-1}$. The SNR improves significantly as we
combine multiple PCs. Combining PC = 33-35, we can achieve ${\rm SNR}\sim 10$
in the range $0.02<k<0.07\,{\rm Mpc}^{-1}$ and the SNR is greater than unity
up to $k\sim 0.6\,{\rm Mpc}^{-1}$. Considering PC = 32-36, we find ${\rm
SNR}\sim 10$ in the range $0.02<k<0.09\,{\rm Mpc}^{-1}$ and the SNR is greater
than unity up to $k\sim 1\,{\rm Mpc}^{-1}$. We find that the SNR improves
further for PC=31-37. However, this improvement is very small, which indicates
that $\Delta{\rm PC}=3$ is indeed an optimal choice for combining multiple
PCs.

##  5 Preliminary results from the observed data

In this section, we present the results from the actual MWA observation
described in Section 2. The analysis is restricted to a single pointing PC=34
which we have considered in section 4.2. For this PC, we have $N_{\rm
Nights}=9$ whereby $\sigma_{N}=20\,{\rm Jy}$ (eq. 1). Considering the measured
visibilities, we find that the amplitudes at small baselines are much larger
than those at the longer baselines. The estimated $C_{\ell}(\Delta\nu)$ is
expected to be correlated across an extent $\delta\ell\sim 40$ (and possibly
larger) due to the tapering (Figure 5, Chatterjee2022), and there is a risk of
the high foreground level at the small baselines leaking out to the larger
baselines. To avoid this, we have discarded the baselines shorter than
$6\lambda$ for the subsequent analysis. We also note an earlier work Li2019,
which has discarded the baselines shorter than $12\lambda$ to avoid foreground
leakage into the EoR window. We have applied the TTGE on the data to estimate
$C_{\ell}(\Delta\nu)$, for which the results are presented in A. We have used
these $C_{\ell}(\Delta\nu)$ values to estimate the cylindrical power spectrum
$P(k_{\perp},k_{\parallel{}})$ which we now discuss in some detail.

Figure 9: The left panel shows the cylindrical PS $\mid
P(k_{\perp},k_{\parallel{}})\mid$. The grey and black dashed lines show the
theoretically predicted boundary of foreground contamination expected from a
monochromatic source located at the horizon and the FWHM of the telescope’s
PB, respectively. The region inside the black rectangle is used to constrain
the 21-cm signal. The middle panel shows the histogram of the quantity
$X=P(k_{\perp},k_{\parallel{}})/\delta P_{N}(k_{\perp},k_{\parallel{}})$
considering the modes inside the rectangle. The right panel shows
$\mid\Delta^{2}(k)\mid$ the absolute values of the mean squared brightness
temperature fluctuations and the corresponding $2\sigma$ error bars. The
negative values of $\Delta^{2}(k)$ are marked with a cross.

The left panel of Figure 9 shows the estimated $\mid
P(k_{\perp},k_{\parallel{}})\mid$ that is found to have values in the range
$10^{7}-10^{15}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$. The grey dashed line shows the
boundary of the foreground wedge, which corresponds to the foreground
contamination expected from a monochromatic source located at the horizon. The
black dashed line shows the same for a source located at $\theta=23^{\circ}$,
which corresponds to the FWHM of the telescope’s PB. We find that the large
values $\mid P(k_{\perp},k_{\parallel{}})\mid>10^{15}\,{\rm mK^{2}}\,{\rm
Mpc^{3}}$, which are due to the foregrounds, are localised within the FWHM
line. We also notice a relatively smaller foreground contribution $(\mid
P(k_{\perp},k_{\parallel{}})\mid\sim 10^{13}\,{\rm mK^{2}}\,{\rm Mpc^{3}})$
between the FWHM line and the wedge boundary. In addition to this, we also
find a significant amount of foreground leakage beyond the wedge boundary, and
the values of $\mid P(k_{\perp},k_{\parallel{}})\mid$ gradually decrease to
$\leq 10^{11}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$ at $k_{\parallel{}}\sim 0.3\,{\rm
Mpc}^{-1}$ and beyond.

Figure 10:  This figure shows $\mid P(k_{\perp},k_{\parallel{}})\mid$ as a
function of $k_{\parallel{}}$ for fixed values of $k_{\perp}$. The blue and
orange curves are from the observed data. The green and red curves are from
simulations with $P(k)\propto k^{-1}$ and $P(k)\propto k^{-2}$, respectively.

The curve corresponding to $k_{\perp}=0.14\,{\rm Mpc}^{-1}$ (for the data) has
been divided by a factor of $10^{3}$, while the curves corresponding to the
simulations have been scaled arbitrarily for better visualization. The grey
shaded region shows $2\sigma$ uncertainties for the simulation with
$P(k)\propto k^{-2}$.

We notice a few horizontal streaks in the power spectrum that extend across
several $k_{\perp}$ bins at some fixed values of $k_{\parallel{}}$. In
particular, we find two very closely spaced horizontal streaks at
$k_{\parallel{}}\approx 0.29\,{\rm Mpc^{-1}}$ where $\mid
P(k_{\perp},k_{\parallel{}})\mid\sim 10^{12}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$
that is in excess of the values ($\mid P(k_{\perp},k_{\parallel{}})\mid\sim
10^{11}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$) in the neighbouring $k_{\parallel}$
bins. Note that this is the lowest $k_{\parallel{}}$ value where we have a
horizontal streak. The horizontal streaks at larger $k_{\parallel}$ are not as
pronounced, and in some cases they do not appear to extend across the entire
$k_{\perp}$ range. In order to understand these better, we consider Figure 10
which shows $\mid P(k_{\perp},k_{\parallel{}})\mid$ as a function of
$k_{\parallel{}}$ for two fixed values of $k_{\perp}$ ($=0.01$ and $0.14\,{\rm
Mpc}^{-1}$). Considering $k_{\perp}=0.01\,{\rm Mpc}^{-1}$, which is the lowest
$k_{\perp}$ bin, we find spikes in the power spectrum at a regular interval of
$\Delta k_{\parallel{}}\sim 0.29\,{\rm Mpc}^{-1}$. The power in these peaks is
roughly of the order of $10^{11}-10^{12}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$, which
is nearly an order of magnitude larger than the adjacent values of $\mid
P(k_{\perp},k_{\parallel{}})\mid$. However note that these spikes are nearly
$3-4$ orders of magnitude smaller than the peak foreground power $\mid
P(k_{\perp},k_{\parallel{}})\mid\approx 10^{16}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$
that occurs at $k_{\parallel{}}=0$. We find that $k_{\parallel{}}=0.29\,{\rm
Mpc^{-1}}$, which is the position of the first spike and also the spacing
between the successive spikes, corresponds to a period of $\Delta\nu_{\rm
per}=1.28\,{\rm MHz}$ in the measured $C_{\ell}(\Delta\nu)$. Note that this is
exactly the spacing between the coarse bands of MWA (Figure 2), which is also
the period of the flagged channels. We find that for several $\ell$ values,
the measured $C_{\ell}(\Delta\nu)$ (Figure 11) exhibits tiny oscillatory
features that become clearly visible once we fit and subtract out the dominant
smoothly varying component. These oscillatory features, which are shown in the
inset of the panels of Figure 11, turn out to have a period of $\Delta\nu_{\rm
per}=1.28\,{\rm MHz}$. Note that the amplitude of these oscillatory components
is around two orders of magnitude smaller than the total measured
$C_{\ell}(\Delta\nu)$. We can attribute the spikes seen in $\mid
P(k_{\perp},k_{\parallel{}})\mid$ for $k_{\perp}=0.01\,{\rm Mpc}^{-1}$ (
Figure 10) to the oscillatory features in $C_{\ell}(\Delta\nu)$ for $\ell=51$.
We next consider $k_{\perp}=0.14\,{\rm Mpc}^{-1}$ which is the $17$-th
$k_{\perp}$ (and $\ell$) bin. In this case, we do not see such prominent
spikes as those seen for $k_{\perp}=0.01\,{\rm Mpc}^{-1}$. Considering
$C_{\ell}(\Delta\nu)$ for $\ell=1273$, here also we see that the oscillatory
pattern is not so well defined as for $\ell=51$. For comparison, Figure 10
shows $\mid P(k_{\perp},k_{\parallel{}})\mid$ for $k_{\perp}=0.01\,{\rm
Mpc}^{-1}$ considering the simulations where $P(k)\propto k^{-1}$ ($s=-1$ in
eq. 20), for which the results are shown in Figure 4. The results from these
simulations, which have exactly the same periodic pattern of flagged channels
as the actual MWA data, do not show any of the spikes and other artefacts seen
in the corresponding results for the actual data. However, the results for
these simulations only exhibits a 2 order of magnitude dynamic range, whereas
the periodic structures in the data exhibits a 4 order of magnitude dynamic
range with respect to $k_{\parallel}=0$. In order to achieve a higher dynamic
range in the simulations, we have also considered $P(k)\propto k^{-2}$ for
which the results are also shown in Figure 10. The simulation now exhibits a 4
order of magnitude dynamic range, for which we still do not see the spikes. We
however notice enhanced statistical fluctuations at some of the
$k_{\parallel}$ where we have spikes in the actual data, which arise because
these $k_{\parallel}$ values are poorly sampled due to the periodic pattern of
flagged channels. We further note that the large dynamic range in the actual
data is a direct consequence of the fact that foregrounds that dominate these
observations are essentially two dimensional (in the plane of the sky), and
their smooth spectral behaviour leads to a large value at $k_{\parallel}=0$.
It is difficult to replicate this in the three dimensional simulations that we
have considered here. We plan to consider foreground simulations in future
work.

The results from our simulations suggests that the spikes and other artefacts
seen in the measured $\mid P(k_{\perp},k_{\parallel{}})\mid$ presented here
are not just a straightforward consequence of the periodic pattern of flagged
channels. Our results seem to indicate that the spikes are possibly caused by
some systematic effect in the data, which has the same period as the missing
channels, periodic systematic errors in the calibration being one such
possibility. We further note that it may be possible to model and remove these
artefacts by subtracting out small period components from the measured
$C_{\ell}(\Delta\nu)$ using Gaussian process regression or some similar
technique Mertens2018; Elahi2023b.

We next attempt to constrain the $z_{c}=8.2$ EoR 21-cm signal using the
measured $P(k_{\perp},k_{\parallel{}})$. To avoid significant foreground
contamination, it is necessary to restrict the region of the
$(k_{\perp},k_{\parallel{}})$ plane which is used to estimate the 21-cm
signal. Here, we have limited the $k_{\parallel}$ range to choose a region
that is well above the foreground wedge. We further wish to avoid the large
spikes seen in the measured $\mid P(k_{\perp},k_{\parallel{}})\mid$. We find
that these spikes are relatively more prominent in the first six $k_{\perp}$
bins where they extend over the entire $k_{\parallel{}}$ range. In contrast,
the spikes are not so prominent in the subsequent $k_{\perp}$ bins where they
also do not extend out to the large $k_{\parallel{}}$ range. This feature is
visible when we compare the results for $k_{\perp}=0.01$ and $0.14\,{\rm
Mpc}^{-1}$ in Figure 10. Based on the above considerations we have discarded
the first six $k_{\perp}$ bins and chosen the rectangular region $0.05\leq
k_{\perp}\leq 0.16\,{\rm Mpc^{-1}}$ and $0.9\leq k_{\parallel{}}\leq 4.6\,{\rm
Mpc^{-1}}$, which is shown in the top right corner of the left panel of Figure
9. The subsequent analysis is restricted to this rectangular region where
$\mid P(k_{\perp},k_{\parallel{}})\mid$ is found to have values in the range
$10^{7}-10^{11}\,{\rm mK^{2}}\,{\rm Mpc^{3}}$.

We now consider the statistics of the measured $P(k_{\perp},k_{\parallel{}})$,
particularly to assess whether these values can be utilized to meaningfully
constrain the EoR 21-cm signal. Following Pal2020, we consider the quantity

$X=\frac{P(k_{\perp},k_{\parallel{}})}{\delta
P_{N}(k_{\perp},k_{\parallel{}})}\,$ (22)

which is the ratio of the measured $P(k_{\perp},k_{\parallel{}})$ to $\delta
P_{N}(k_{\perp},k_{\parallel{}})$ the statistical uncertainty expected from
the system noise contribution only. We expect $X$ to have a symmetric
distribution with mean $\mu=0$ and standard deviation $\sigma_{\rm Est}=1$ if
the estimated values of $P(k_{\perp},k_{\parallel{}})$ are consistent with the
uncertainties predicted due to the system noise contribution only. Residual
foreground contamination will manifest itself as a positive mean, whereas a
negative mean would indicate some sort of negative systematics in the power
spectrum estimates. A value $\sigma_{\rm Est}>1$ would indicate the presence
of additional sources of uncertainty, beyond the predicted system noise
contribution.

The middle panel of Figure 9 shows the probability density function (PDF) of
$X$, which is found to be largely symmetric around $0$ with $\mu=5.88$ and
$\sigma_{\rm Est}=95.81$. Here, the substantially large value of $\sigma_{\rm
Est}$ indicates that the measured $P(k_{\perp},k_{\parallel{}})$ has
fluctuations which are considerably larger than those predicted from the
system noise contribution alone. We note that several earlier works have
reported such ‘excess variance’ for the estimated 21-cm PS at $\sim 150\,{\rm
MHz}$ Mertens2020; Pal2020 and also at higher frequencies $\sim 430\,{\rm
MHz}$ Pal2022; Elahi2023; Elahi2023b; Elahi2024. The exact cause of this
excess variance is not known at present. We note that the earlier works which
have applied the TGE Pal2020; Pal2022; Elahi2023; Elahi2023b; Elahi2024 have
obtained values of $\sigma_{\rm Est}$ in the range $2$ to $5$, whereas the
value $\sigma_{\rm Est}=95.81$ obtained here is more than an order of
magnitude larger. Although the exact cause is not known, it is plausible that
the systematics responsible for the spikes in the PS at the lower values of
$k_{\perp}$ and $k_{\parallel{}}$ also causes excess fluctuations in
$P(k_{\perp},k_{\parallel{}})$ within the rectangular region used to constrain
the 21-cm signal. We find that the mean $\mu=5.88$ is not consistent with zero
within the expected statistical fluctuations predicted using $\sigma_{\rm
Est}$. This indicates that the measured $P(k_{\perp},k_{\parallel{}})$ has
some residual foreground contamination which introduces a positive bias in the
distribution if $X$. We further note that the PDF of $X$ is not Gaussian, and
it nicely fit by a Lorentzian distribution with a peak location $x_{0}=1.68$
and spread $\gamma=2.73$, consistent with the earlier studies Elahi2023;
Elahi2023b; Elahi2024 all of which find that $X$ follows a Lorenztian
distribution. For the subsequent analysis, we have scaled the system noise
only error predictions with $\sigma_{\rm Est}$ to account for the excess
variance i.e. we have used $\delta P(k_{\perp},k_{\parallel{}})=\sigma_{\rm
Est}\times\delta P_{N}(k_{\perp},k_{\parallel{}})$ to predict the statistical
fluctuations expected in the measured $P(k_{\perp},k_{\parallel{}})$.

We use the $(k_{\perp},k_{\parallel{}})$ modes in the rectangular region
mentioned earlier to estimate the spherical power spectrum $P(k)$. The right
panel of Figure 9 shows the mean squared brightness temperature fluctuations
$\Delta^{2}(k)=k^{3}P(k)/(2\pi^{2})$ as a function of $k$ along with the
$2\sigma$ error bars corresponding to the predicted statistical fluctuations
$\delta P(k_{\perp},k_{\parallel{}})$. The expected statistical fluctuations
are much larger than the EoR 21-cm signal which is predicted to typically have
values of the order of $\Delta^{2}(k)\sim 10^{1}-10^{2}\,{\rm mK}^{2}$ over
the range $k=1-4\,{\rm Mpc}^{-1}$ considered here Mondal2017. In the absence
of foreground contamination, we expect $\Delta^{2}(k)$ to fluctuate around $0$
with both positive and negative values that are consistent with the predicted
statistical fluctuations. We find that $\Delta^{2}(k)$ in most of the $k$-bins
are consistent with the statistical fluctuations, and the values are within
$0\pm 2\sigma$. However, the values of $\Delta^{2}(k)$ in the first, second
and sixth $k$-bins are somewhat larger, and these are consistent with the
statistical fluctuations only at $0\pm 5\sigma$. It is possible that these
$k$-bin may have some residual foreground contamination, however, it is also
possible that these are more extreme statistical fluctuations in the extended
tail of the Lorentzian distribution. Note that the two $k$-bins with negative
$\Delta^{2}(k)$ values are both consistent with $0\pm 2\sigma$ indicating that
the estimated PS is free from any negative systematics. Since the measured PS
is largely consistent with the statistical fluctuations and there are no
negative systematics, we use the measured values to place $2\sigma$ upper
limits on $\Delta^{2}(k)$. The tightest upper limit is found to be
$\Delta^{2}(k)<(1.85\times 10^{4})^{2}\,{\rm mK^{2}}$ at the first $k$-bin
$k=1\,{\rm Mpc}^{-1}$.

##  6 Summary and Discussion

Drift scan observations provide an economic and stable option with the broad
sky coverage required for 21-cm EoR experiments. In this paper, we consider
the Tracking Tapered Gridded Estimator (TTGE; Chatterjee2022), which aims to
measure the power spectrum directly from the visibilities recorded in radio
interferometric observations in the drift scan mode. The estimator is based on
the TGE Choudhuri2016b that has been widely used for analysing GMRT/uGMRT
radio interferometric observations Choudhuri2020; Pal2020; Pal2022; Elahi2023;
Elahi2023b; Elahi2024. The TGE works with gridded visibilities, a feature that
makes it computationally efficient and also allows it to taper the sky
response so as to suppress foreground contamination from the sources far away
from the telescope’s pointing direction. While the original TGE was developed
for observations where the telescope tracks a single field. The TTGE has been
designed to work with drift scan observation. This estimator is tailored for
observations with radio-interferometric arrays where the individual elements
(antennas) are fixed on the ground (e.g. LOFAR, MWA, HERA, CHIME, upcoming
SKA-low), and in the present work we have considered data from drift scan
observations Patwa2021 with the compact configuration of MWA-II Tingay2013. In
drift scan observations, the sky moves across the telescope’s field of view
with the passage of time. The TTGE allows us to follow a fixed position on the
sky, namely the tracking center (TC), as it moves across the telescope’s field
of view. It further allows us to taper the sky response to a small angular
region around the TC so as to suppress the foreground contamination from
bright sources located at large angular separations from the TC.

In Chatterjee2022, we have presented the mathematical framework for the TTGE
and validated this for a single frequency channel. The analysis there
considers an input model angular power spectrum $C_{\ell}^{M}$, for which the
simulated sky signal and resulting drift scan visibilities are processed
through the TTGE pipeline. The estimated $C_{\ell}$ is found to be in good
agreement with $C_{\ell}^{M}$, thereby validating the TTGE. The present work
validates the TTGE considering observations that cover multiple frequency
channels spanning a finite bandwidth. Here, we have considered an input model
power spectrum $P^{m}(k)=(1{\rm Mpc}^{-1}/k)\,{\rm K^{2}}\,{\rm Mpc^{3}}$, for
which simulated all-sky maps were used to calculate drift scan visibilities.
These were processed through the TTGE pipeline to estimate the MAPS
$C_{\ell}(\Delta\nu)$ which was found to be in good agreement with the
predicted $C_{\ell}^{M}(\Delta\nu)$ corresponding to $P^{\rm m}(k)$ (Figure
3). We have used a Maximum Likelihood estimator Pal2022; Elahi2023 to directly
compute both $P(k_{\perp},k_{\parallel{}})$ and $P(k)$ from the estimated
$C_{\ell}(\Delta\nu)$. The estimated $P(k)$ is found to be in good agreement
with the input model $P^{\rm m}(k)$ (Figure 5), thereby validating the TTGE
for 3D multi-frequency data.

The MWA has a periodic pattern of flagged channels, the period being
$\Delta\nu_{\rm per}=1.28\,{\rm MHz}$ (Figure 2). In addition, there is the
possibility of other frequency channels also being flagged in order to avoid
man-made radio frequency interference (RFI), etc. Several earlier works which
have first transformed from frequency to delay space, and then estimated the
PS (e.g. Paul2016; Li2019; Trott2020; Patwa2021) have reported a very
prominent pattern of horizontal streaks in the estimated
$P(k_{\perp},k_{\parallel{}})$ due to the periodic pattern of flagged channels
present in the MWA data. In our approach Bharadwaj2018, we first estimate
$C_{\ell}(\Delta\nu)$, which does not have any missing $\Delta\nu$ even when
the visibility data contains flagged channels. We then Fourier transform
$C_{\ell}(\Delta\nu)$ to obtain a clean estimate of
$P(k_{\perp},k_{\parallel{}})$. Considering the simulated data, we see that
the missing frequency channels do not introduce any artefacts in the estimated
$P(k_{\perp},k_{\parallel{}})$ (Figure 4). We find that the estimated
$P(k_{\perp},k_{\parallel{}})$ are visually indistinguishable irrespective of
whether there are flagged channels or not. We conclude that our estimator is
impervious to the periodic pattern of flagged channels, and these do not
introduce any artefacts in the estimated $P(k_{\perp},k_{\parallel{}})$.

We have applied the TTGE to estimate the PS for a single pointing of the
actual MWA data, which effectively corresponds to using the TGE. The pointing
center PC=34 corresponds to $({\rm RA,DEC)}=(6.1^{\circ},-26.7^{\circ})$
(Figure 1) with an observing time of approximately $17\,{\rm min}$. For the
tapering, we have used a Gaussian window function with $\theta_{\rm
FWHM}=15^{\circ}$. The estimated MAPS $C_{\ell}(\Delta\nu)$, arising mainly
from foregrounds, is shown in Figure 11. This does not appear to exhibit any
prominent features which can be associated with the periodic pattern of
flagged channels. However, on fitting and subtracting out a smooth polynomial
in $\Delta\nu$, for several $\ell$ values the residual $C_{\ell}(\Delta\nu)$
exhibit oscillations of period $\Delta\nu_{\rm per}=1.28\,{\rm MHz}$. The
amplitude of these oscillating residuals is approximately two orders of
magnitude smaller than the total $C_{\ell}(\Delta\nu)$. The left panel of
Figure 9 shows the cylindrical power spectrum $P(k_{\perp},k_{\parallel{}})$
estimated from the measured $C_{\ell}(\Delta\nu)$. We find that the
foregrounds are mainly localized within the expected foreground wedge
boundary, with some leakage beyond. We notice a few horizontal streaks in the
power spectrum that extend across several $k_{\perp}$ bins at some fixed
$k_{\parallel{}}$ values.

To understand the origin of these streaks, we consider Figure 10 where we find
the power spectrum at $k_{\perp}=0.01\,{\rm Mpc}^{-1}$ shows spikes at a
regular interval of $\Delta k_{\parallel{}}\sim 0.29\,{\rm Mpc}^{-1}$ which
matches $\Delta\nu_{\rm per}=1.28\,{\rm MHz}$ in $C_{\ell}(\Delta\nu)$. We
attribute the spike in $P(k_{\perp},k_{\parallel{}})$ to the periodic
oscillation in $C_{\ell}(\Delta\nu)$. We also find that
$P(k_{\perp},k_{\parallel{}})$ does not exhibit such prominent spikes at the
larger $k_{\perp}$ bins. Further, such spikes are also not present in the
$P(k_{\perp},k_{\parallel{}})$ estimated from simulations, which incorporate
the periodic pattern of flagged channels. Note, however, that the dynamic
range of the present simulations is an order of magnitude lower than the
actual data. In future work, we plan to carry out foreground simulations where
we expect to achieve a higher dynamic range.

The origin of the small periodic oscillation in $C_{\ell}(\Delta\nu)$, which
is responsible for the spikes in $P(k_{\perp},k_{\parallel{}})$, is not known
at present. The fact that these features are not present in the simulations
suggests that our estimator is impervious to the periodic pattern of flagged
channels, and the spikes are not just a straightforward consequence of the
missing channels. Our results seem to indicate that the spikes are possibly
caused by some systematic effect in the data, which has the same period as the
missing channels, periodic systematic errors in the calibration being one such
possibility. We note that the streaks found here are an order of magnitude
smaller than those found in earlier analyses of MWA data (e.g. Patwa2021). It
may be possible to model and remove these artefacts by subtracting out small
period components from the measured $C_{\ell}(\Delta\nu)$ using Gaussian
process regression or some similar technique Mertens2018; Elahi2023b. We plan
to address this in future work.

We identify a rectangular region $0.05\leq k_{\perp}\leq 0.16\,{\rm Mpc^{-1}}$
and $0.9\leq k_{\parallel{}}\leq 4.6\,{\rm Mpc^{-1}}$ to be relatively free of
foreground contamination, spikes, and other artefacts. We use the
$P(k_{\perp},k_{\parallel{}})$ estimates in this region to constrain the EoR
21-cm signal. The $P(k_{\perp},k_{\parallel{}})$ values, suitably scaled with
the uncertainties expected from the system noise, are found to be nearly
symmetrically distributed around 0. The distribution is well-fitted with a
Lorentzian profile. However, the standard deviation is roughly 100 times of
that expected from the contribution from system noise alone. The cause of this
excess variance is not known. We have used our measurements to place a
$2\sigma$ upper limit $\Delta^{2}(k)<(1.85\times 10^{4})^{2}\,{\rm mK^{2}}$ on
the mean squared 21-cm brightness temperature fluctuations at $k=1\,{\rm
Mpc}^{-1}$. The present results, which are restricted to a single pointing,
are limited by foregrounds and the systematics whose origins are not known to
us at present. It is necessary to address these issues before combining
multiple pointings. We plan to address these issues in future work.

##  Appendix A Preliminary results: MAPS

Figure 11: This figure shows $C_{\ell}(\Delta\nu)$ for the annotated $\ell$
values. The blue dashed curves show a polynomial fit on the range
$\Delta\nu<6$ MHz. The polynomial fit is subtracted from the measured
$C_{\ell}(\Delta\nu)$, and the residual $C_{\ell}(\Delta\nu)$ are shown in the
insets (red).

We have applied the TTGE on MWA drift scan observation to estimate the MAPS
$C_{\ell}(\Delta\nu)$ (Section 5). Figure 11 shows $C_{\ell}(\Delta\nu)$ for
different $\ell$ whose values are annotated in each panel. The estimated
$C_{\ell}(\Delta\nu)$ have values $\sim 2\times 10^{5}\,{\rm mK}^{2}$ at
$\Delta\nu=0$ and it gradually decreases with increasing $\Delta\nu$. This
gradual decorrelation is a typical behaviour of foregrounds, which has a
smooth frequency dependence. Further, the decorrelation becomes faster with
increasing $\ell$ due to baseline migration Pal2022. We note that point
sources, which are in the sidelobes of the telescopes, introduce oscillatory
features in $C_{\ell}(\Delta\nu)$, which are not prominent in these plots.
This is because the contribution from the sources within the field of view of
the primary beam is much larger than those appearing through the sidelobes. We
plan to identify and remove the point sources that are in the field of view of
the primary beam in future work.

We see no prominent features in the estimated $C_{\ell}(\Delta\nu)$ but find
features in the power spectrum (Figure 9), which leads us to investigate the
measured $C_{\ell}(\Delta\nu)$ closely. The idea is to fit the measured
$C_{\ell}(\Delta\nu)$ with a polynomial and subtract the polynomial fit to
look for any features in the residual $C_{\ell}(\Delta\nu)$. Here we have
chosen the range $\Delta\nu<6$ MHz and fitted the estimated
$C_{\ell}(\Delta\nu)$ with an even polynomial having 10 terms. The blue dashed
curves on top of the measured $C_{\ell}(\Delta\nu)$ show the polynomial fits,
and the residual $C_{\ell}(\Delta\nu)$ are shown in the insets (red).

## Acknowledgements

S. Chatterjee acknowledges support from the South African National Research
Foundation (Grant No. 84156) and the Department of Atomic Energy, Government
of India, under project no. 12-R&D-TFR-5.02-0700. S. Chatterjee would also
like to thank Dr. Devojyoti Kansabanik, for helpful discussions.

## Data Availability

The data sets were derived from sources in the public domain (the MWA Data
Archive: project ID G0031) at https://asvo.mwatelescope.org/.


Derangements and Continued Fractions for $e$

Peter<EMAIL_ADDRESS>
School of Mathematics & Statistics, University College Dublin

###### Abstract

Several continued fraction expansions for $e$ have been produced by an
automated conjecture generator (ACG) called _The Ramanujan Machine_. Some of
these were already known, some have recently been proved and some remain
unproven. While an ACG can produce interesting putative results, it gives very
limited insight into their significance. In this paper, we derive an elegant
continued fraction expansion, equivalent to a result from the Ramanujan
Machine, using the sequence of ratios of factorials to subfactorials or
derangement numbers.

## Arrangements and Derangements

Six students entering an examination hall place their cell-phones in a box.
After the exam, they each grab a phone at random as they rush out. What is the
likelihood that none of them gets their own phone? The surprising answer —
about 37% whatever the number of students — emerges from the theory of
derangements.

We may call any permutation of the elements of a set an arrangement. A
_derangement_ is an arrangement for which every element is moved from its
original position. Thus, a derangement is a permutation that has no _fixed
points._ The number of derangements of a set of $n$ elements is also called
the _subfactorial_ of $n$. Various notations are used for subfactorials: $!n$,
$d_{n}$ and $n$¡ are common; we will use $!n$ (read as ‘bang-en’).

Derangements were first considered by Pierre Reymond de Montmort. In 1713,
with help from Nicholas Bernoulli, he managed to find an expression for the
connection between $!n$ and $n!$. The answer, which he obtained using the
inclusion-exclusion principle (Zeilberger, 2008, pg. 560), is

$!n=n!\left(1-\frac{1}{1!}+\frac{1}{2!}-\frac{1}{3!}+\frac{1}{4!}-\cdots\pm\frac{1}{n!}\right)=n!\sum_{k=0}^{n}\frac{(-1)^{k}}{k!}\,.$
(1)

Of course, we see from this that $\lim_{n\to\infty}(!n)=n!/e$. In fact, we can
write a more precise connection between derangements and arrangements:

$!n=\left\lfloor\frac{n!+\textstyle{\frac{1}{2}}}{e}\right\rfloor\,.$

This means that $!n$ is the nearest whole number to $n!/e$.

The number $!n$ of derangements of an $n$-element set may be calculated using
a second-order recurrence relation:

$!n=(n-1)[!(n-1)+!(n-2)]$

with $!0=1$ and $!1=0$. The subfactorials also satisfy a first-order
recurrence relation,

$!n=n\times!(n-1)+(-1)^{n}\,,\qquad[\mbox{compare\ \ }n!=n\times(n-1)!]$

with initial condition $!0=1$. The first eight values of $!n$ are
$\\{1,0,1,2,9,44,265,1\,854\\}$.

There is an integral expression for the subfactorial,

$!n=\int_{0}^{\infty}(x-1)^{n}e^{-x}\mathrm{d}x\,,\qquad\left[\mbox{compare\ \
}n!=\int_{0}^{\infty}x^{n}e^{-x}\mathrm{d}x\right]\,.$ (2)

Expansion of (2) yields de Montmort’s result (1). It also allows extension of
the subfactorial function to non-integral arguments ($!x$) and analytic
continuation to the complex plane ($!z$).

## Continued Fractions and Convergents

The continued fraction expansion of an irrational number $x$ is written, in
expanded form (centre) and concise form (right), as

$x=a_{0}+\cfrac{1}{a_{1}+\cfrac{1}{a_{2}+\cfrac{1}{a_{3}+\ddots}}}=[a_{0};a_{1},a_{2},a_{3},\dots]$

where $a_{n}$ are integers. If $a_{n}$ is positive for $n\geq 1$ this is
called a _simple_ continued fraction.

The _generalized_ continued fraction expansion is written

$x=b_{0}+\cfrac{a_{1}}{b_{1}+\cfrac{a_{2}}{b_{2}+\cfrac{a_{3}}{b_{3}+\ddots}}}=b_{0}+\frac{a_{1}}{b_{1}+}\,\frac{a_{2}}{b_{2}+}\,\frac{a_{3}}{b_{3}+}\,\frac{a_{4}}{b_{4}+}\,\cdots\,.$

We assume that $a_{n}$ and $b_{n}$ are integers. Truncating the expansion at
various points, we obtain the _convergents_

$r_{n}=\frac{p_{n}}{q_{n}}=b_{0}+\frac{a_{1}}{b_{1}+}\,\frac{a_{2}}{b_{2}+}\,\frac{a_{3}}{b_{3}+}\,\frac{a_{4}}{b_{4}+}\,\cdots\,\frac{a_{n}}{b_{n}}$

where the numerators and denominators, $p_{n}$ and $q_{n}$, are integers. We
define the starting values

$p_{-1}=1\,,\qquad q_{-1}=0\,,\qquad p_{0}=b_{0}\,,\qquad q_{0}=1\,.$

Then, $p_{k}$ and $q_{k}$ for $k\geq 1$ are given by recurrence relations:

$p_{k}=b_{k}p_{k-1}+a_{k}p_{k-2}\,,\qquad q_{k}=b_{k}q_{k-1}+a_{k}q_{k-2}\,,$
(3)

which may be proved by induction (Jones & Thron, 1980, Pg. 20).

This process can be inverted: given a sequence of numerators $p_{n}$ and
denominators $q_{n}$ (or just their ratios, the convergents
$r_{n}=p_{n}/q_{n}$), we can solve (3) for $a_{n}$ and $b_{n}$:

$a_{n}=\frac{p_{n-1}q_{n}-p_{n}q_{n-1}}{p_{n-1}q_{n-2}-p_{n-2}q_{n-1}}\,,\qquad
b_{n}=\frac{p_{n}q_{n-2}-p_{n-2}q_{n}}{p_{n-1}q_{n-2}-p_{n-2}q_{n-1}}$ (4)

together with the starting values $b_{0}=p_{0}$, $a_{1}=(p_{1}-b_{0}q_{1})$
and $b_{1}=q_{1}$.

## Continued Fractions for $e$

Euler’s number is usually defined as the limit
$e=\lim_{n\to\infty}(1+1/n)^{n}$, which is the limit of the sequence

$\left\\{\frac{2^{1}}{1^{1}},\frac{3^{2}}{2^{2}},\frac{4^{3}}{3^{3}},\dots,\frac{(n+1)^{n}}{n^{n}},\dots\right\\}$

The terms may be regarded as the convergents of a continued fraction,

$r_{n}=\frac{p_{n}}{q_{n}}\,,\qquad\mbox{where}\qquad
p_{n}=(n+1)^{n}\quad\mbox{and}\quad q_{n}=n^{n}\,.$

We can generate a continued fraction by using (4). It begins as

$1+\frac{1}{1-}\,\frac{1}{5-}\,\frac{13}{10-}\,\frac{491}{196-}\,\frac{487903}{9952-}\,\frac{2384329879}{958144-}\,\cdots\,.$
(5)

The error of this expansion ($\log_{10}[r_{n}-e]$) as a function of truncation
is shown in Fig. 1 (dashed line). It is clear that the convergence is very
slow.

Euler made extensive studies of continued fractions. For example, his 50-page
paper, _Observations on continued fractions_ (Euler, 1750), contains numerous
original results. One of his best-known expansions is

$e=[2;1,2,1,1,4,1,1,6,1,1,8,\dots]$ (6)

The error of Euler’s expansion is shown in Fig. 1 (dotted line). It converges
much faster than (5). There is a clear signal of period 3, consistent with the
recurring pattern $(1,1,n)$ in (6).

## Continued fraction from derangement numbers

A beautiful continued fraction emerges from the relationship between
arrangements and derangements. We saw above that

$\left[\frac{\mbox{Arrangements of $n$ elements}}{\mbox{Derangements of $n$
elements}}\right]=\frac{n!}{!n}\to e$

If we define the numerators and denominators of convergents to be

$p_{n}=n!\quad\mbox{and}\quad q_{n}=\,!n\,,$

we can solve for the factors $a_{n}$ and $b_{n}$. The starting values
$p_{0}=1,p_{1}=1,q_{0}=1,q_{1}=0$ yield $a_{0}=0,b_{0}=1,a_{1}=1,b_{1}=0$.
Then (4) may be solved to yield $a_{n}=b_{n}=(n-1)$ for $n\geq 2$. Thus we get
the expansion

$e=1+\frac{1}{0+}\,\frac{1}{1+}\,\frac{2}{2+}\,\frac{3}{3+}\,\frac{4}{4+}\,\cdots\,.$

A small adjustment enables us to write this in the elegant form

$e=2+\frac{2}{2+}\,\frac{3}{3+}\,\frac{4}{4+}\,\frac{5}{5+}\,\frac{6}{6+}\,\cdots\,.$
(7)

The error of (7) is shown in Fig. 1 (solid line). Convergence is more rapid
than for the other two expansions.

Figure 1: Logarithm of the error $\log_{10}|r_{n}-e|$ in the continued
fraction expansions for $e$. Dashed line: $r_{n}=(1+1/n)^{n}$, Eq. (5) .
Dotted line: Convergents of Euler’s expansion (6). Solid line:
$r_{n}=(n+1)!/!(n+1)$, Eq. (7).

## The Ramanujan Machine

An Automated Conjecture Generator (ACG) called _The Ramanujan Machine_ 222G.
H. Hardy, in his Introduction to Ramanujan’s _Collected Papers_ (1927), wrote
that Ramanujan’s mastery of of continued fractions was “beyond that of any
mathematician in the world”. has been implemented by a team of mathematicians
at the Israel Institute of Technology. This ACG system is capable of producing
conjectures about mathematical (and physical) constants, expressed in the form
of continued fractions, using only numerical data as input. A paper describing
the system is available on the arXiv preprint server (Raayoni, et al., 2020).

The Ramanujan Machine comprises algorithms designed to discover new
conjectures, running on a network of computers. The goal of the project is to
formulate conjectures that may then be proved mathematically. The ACG has
already generated a number of very interesting new conjectures, as well as
reproducing several results that were already well known. The website
(http://www.ramanujanmachine.com/) enables researchers to submit proofs of
conjectures, code new algorithms and (if they wish) allow access to their
computers for distributed computation.

While the Ramanujan Machine generates conjectures but not proofs, it has
inspired a complementary project using _symbolic_ rather than numerical
computation. Dougherty-Bliss and Zeilberger (2020) describe a system that
generates automatic proofs of continued fraction expansions. Their system
produced some infinite families of expansions together with rigorous proofs of
their validity.

One of the continued fractions discovered by the Ramanujan Machine is

$\frac{1}{e-1}=\frac{1}{1+}\,\frac{2}{2+}\,\frac{3}{3+}\,\frac{4}{4+}\,\frac{5}{5+}\,\frac{6}{6+}\,\cdots\,,$
(8)

which is easily seen to be equivalent to (7) above. This is indicated in
Raayoni, et al. (2020) as a “known” result. A proof was presented by Kadyrov
and Mashurov (2019). Lu (2019) gave elementary proofs of other generalized
continued fraction formulae for $e$. However, the connection with derangement
numbers was not made by any of these authors.

Dougherty-Bliss and Zeilberger (2020) proved a generalized expansion of which
the Ramanujan Machine result is a special case. They noted the occurrence of
derangement numbers in their expansion, describing this as a “remarkable
coincidence”, and further commenting that “There does not seem to be any
immediate combinatorial reason for the derangement numbers to appear.” Our
above derivation of (7), starting from the ratio of factorials to
subfactorials, makes the connection clear.

## References

  * [1] Dougherty-Bliss, Robert and Doron Zeilberger, 2020: Automatic conjecturing and proving of exact values of some infinite families of infinite continued fractions. https://arxiv.org/pdf/2004.00090.pdf.
  * [2] Euler, L., 1750: De fractionibus continuis observationes. _Commentarii academiæ scientiarum Petropolitanæ_ , 11, 32–81. Reprinted in _Opera Omnia_ , Series 1, 14, 291–349. Translation by Alexander Aycock: _Observations on continued fractions_. https://arxiv.org/pdf/1808.07006.pdf
  * [3] Gorroochurn, Prakash, 2012: _Classic Problems of Probability_. Wiley, ISBN: 978-1-118-06325-5
  * [4] Jones, William B. and W. J. Thron, 1980: _Continued Fractions: Analytic Theory and Applications_. Encyclopedia of Mathematics and its Applications. No. 11. Addison-Wesley. ISBN 0-201-13510-8.
  * [5] Kadyrov, Shirali and Farukh Mashurov, 2019: Generalized continued fraction expansions for $\pi$ and $e$. https://arxiv.org/pdf/1912.03214.pdf
  * [6] Lu, Zhentao, 2019: Elementary proofs of generalized continued fraction formulae for $e$. https://arxiv.org/pdf/1907.05563.pdf
  * [7] Raayoni, Gal, Shahar Gottlieb, George Pisha, Yoav Harris, Yahel Manor, Uri Mendlovic, Doron Haviv, Yaron Hadad and Ido Kaminer, 2020: The Ramanujan Machine: Automatically generated conjectures on fundamental constants. https://arxiv.org/pdf/1907.00205.pdf
  * [8] Zeilberger, Doron, 2008: Enumerative and algebraic combinatorics. §IV.18 in _The Princeton Companion to Mathematics_ , Ed. Timothy Gowers. Princeton Unversity Press (pp. 550–561). ISBN: 978-0-6911-1880-2.