# Improved upper limits on the 21-cm signal power spectrum at $z = 17.0$ and $z = 20.3$ from an optimal field observed with NenuFAR

S. Munshi<sup>1,★</sup> F. G. Mertens,<sup>1,2</sup> J. K. Chege,<sup>1,3</sup> L. V. E. Koopmans<sup>1</sup>, A. R. Offringa,<sup>3</sup> B. Semelin,<sup>2</sup> R. Barkana,<sup>4</sup> J. Dhandha<sup>5</sup>,<sup>5</sup> A. Fialkov,<sup>5,6</sup> R. Mériot<sup>6</sup>,<sup>2,7</sup> S. Sikder,<sup>4</sup> A. Bracco,<sup>8,9</sup> S. A. Brackenhoff<sup>1</sup>,<sup>1</sup> E. Ceccotti,<sup>1,10</sup> R. Ghara<sup>1</sup>,<sup>11</sup> S. Ghosh,<sup>1</sup> I. Hothi,<sup>9</sup> M. Mevius<sup>1</sup>,<sup>3</sup> P. Ocvirk<sup>1</sup>,<sup>12</sup> A. K. Shaw<sup>1</sup>,<sup>13</sup> S. Yatawatta<sup>1,3</sup> and P. Zarka<sup>14,15</sup>

<sup>1</sup>Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, NL-9700 AV Groningen, the Netherlands

<sup>2</sup>LUX, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, F-75014 Paris, France

<sup>3</sup>ASTRON, PO Box 2, NL-7990 AA Dwingeloo, the Netherlands

<sup>4</sup>School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel

<sup>5</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

<sup>6</sup>Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK

<sup>7</sup>Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2AZ, UK

<sup>8</sup>INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

<sup>9</sup>Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France

<sup>10</sup>INAF – Istituto di Radioastronomia, Via P. Gobetti 101, I-40129 Bologna, Italy

<sup>11</sup>Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur 741 246, West Bengal, India

<sup>12</sup>Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg (ObAS), 11 rue de l'Université, Strasbourg, France

<sup>13</sup>Department of Computer Science, University of Nevada Las Vegas, 4505 S. Maryland Pkwy., Las Vegas, NV 89154, USA

<sup>14</sup>ORN, Observatoire de Paris, Université PSL, Univ Orléans, CNRS, F-18330 Nançay, France

<sup>15</sup>LIRA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

Accepted 2025 August 19. Received 2025 August 18; in original form 2025 July 14

## ABSTRACT

We report the deepest upper limits to date on the 21-cm signal power spectrum during the Cosmic Dawn (redshifts: $z > 15$), using four nights of observations with NenuFAR. The limits are derived from two redshift bins, centred at $z = 20.3$ and $z = 17.0$, with integration times of 26.1 and 23.6 h, from observations of an optimal target field chosen to minimize sidelobe leakage from bright sources. Our analysis incorporates improvements to the data processing pipeline, particularly in subtracting strong radio sources in the primary beam sidelobes and mitigating low-level radio frequency interference, yielding a 50-fold reduction in the excess variance compared to a previous analysis of the north celestial pole field. At $z = 20.3$, we achieve a best $2\sigma$ upper limit of $\Delta_{21}^2 < 4.6 \times 10^5 \text{ mK}^2$ at $k = 0.038 \text{ h cMpc}^{-1}$, while at $z = 17.0$, the best limit is $\Delta_{21}^2 < 5.0 \times 10^6 \text{ mK}^2$ at $k = 0.041 \text{ h cMpc}^{-1}$. These are the strongest constraints on the 21-cm power spectrum at the respective redshifts, with the limit at $z = 20.3$ being deeper by more than an order of magnitude over all previous Cosmic Dawn power spectrum limits. Comparison against simulated exotic 21-cm signals shows that while the $z = 20.3$ limits begin to exclude the most extreme models predicting signals stronger than the EDGES detection, an order-of-magnitude improvement would constrain signals compatible with EDGES. A coherence analysis reveals that the excess variance is largely incoherent across nights for the $z = 20.3$ redshift bin, suggesting that deeper integrations could yield significantly stronger constraints on the 21-cm signal from the Cosmic Dawn.

**Key words:** methods: data analysis – techniques: interferometric – dark ages, reionization, first stars – cosmology: observations.

## 1 INTRODUCTION

The 21-cm line of neutral hydrogen has emerged as one of the most powerful probes of the astrophysics and cosmology of the early Universe. The emerging field of 21-cm cosmology aims to detect the redshifted 21-cm signal from the early Universe, seen in absorption or emission against the background radiation, the cosmic microwave background (CMB). As the first luminous sources formed during the

Cosmic Dawn and subsequently reionized the intergalactic medium during the epoch of reionization (EoR), the astrophysical processes during these periods imprinted distinct signatures on the 21-cm brightness temperature distribution. Detection of these signatures thus offers a potential for direct experimental verification of our theories of early star formation and galaxy evolution (e.g. Furlanetto, Oh & Briggs 2006; Pritchard & Loeb 2012). Over the last decade, two complementary approaches to detecting the 21-cm signal from the early universe have gained momentum: global and interferometric 21-cm experiments.

<footnote>*E-mail: munshi@astro.rug.nl</footnote>
Global 21-cm experiments measure the sky-averaged 21-cm brightness temperature as a function of redshift. The claimed detection of an absorption trough in the global 21-cm signal during the Cosmic Dawn by EDGES<fnref n="1" /> (Bowman et al. 2018) has generated considerable interest in this epoch. The feature detected by EDGES has an unusual shape and is more than twice as deep as the absorption trough that is expected during the Cosmic Dawn based on standard models, which are incompatible with the detection. It thus becomes necessary to invoke exotic non-standard models of the Cosmic Dawn to explain the EDGES detection. These non-standard models fall under two categories. The first class of models proposes additional components to the background temperature, effectively boosting the contrast between the 21-cm signal and the background, which could result in the observed absorption trough (e.g. Ewall-Wice et al. 2018; Dowell & Taylor 2018a; Feng & Holder 2018a; Fialkov & Barkana 2019). Another class of models proposes that alternative mechanisms such as baryon–dark-matter interaction could lead to supercooling of the gas, lowering the spin temperature and enhancing the contrast against the background radiation, thus producing an unusually deep absorption trough (e.g. Barkana 2018; Berlin et al. 2018; Fialkov, Barkana & Cohen 2018; Muñoz & Loeb 2018; Liu et al. 2019). There have been no independent measurements with other instruments that have confirmed the EDGES detection. On the contrary, results from the SARAS<fnref n="2" /> experiment reject the EDGES detection with 95 per cent confidence (Singh et al. 2022). Several other global 21-cm experiments, such as BIGHORNS<fnref n="3" /> (Sokolowski et al. 2015), PRIZM<fnref n="4" /> (Philip et al. 2019), REACH<fnref n="5" /> (de Lera Acedo et al. 2022), MIST<fnref n="6" /> (Monsalve et al. 2024), and RHINO<fnref n="7" /> (Bull et al. 2024), aim to settle this debate, with many of them having already been deployed and recording data.

In contrast to the sky-averaged signal probed by global experiments, interferometers probe the three-dimensional spatial fluctuations in the 21-cm signal brightness temperature distribution and can potentially provide us with a wealth of information about the early universe. However, the signal-to-noise ratio (S/N) requirement for the detection of the 21-cm signal with interferometers is even more demanding than global 21-cm experiments, owing to the faint nature of the signal compared to instrumental thermal noise sensitivities reached in reasonable observation hours using the most sensitive currently operating radio interferometers. Thus, the current generation of experiments attempts to make a statistical detection of the 21-cm fluctuations by estimating its power spectrum. The upcoming SKA<fnref n="8" /> is expected to have the thermal noise sensitivity to not only detect the high-z 21-cm power spectrum, but also to generate tomographic maps of the EoR in reasonable observing hours (Koopmans et al. 2015; Bonaldi et al. 2025). Most interferometric experiments aiming to detect the high-redshift 21-cm signal have focused on the EoR, where they probe the distribution and evolution of ionized bubbles around the first stars and galaxies that eventually reionized the Universe. These efforts have been led by interferometers such as the GMRT<fnref n="9" /> (Paciga et al. 2013), LOFAR<fnref n="10" /> (Patil et al. 2017; Mertens et al.

2020, 2025; Ceccotti et al. 2025a), PAPER<fnref n="11" /> (Kolopanis et al. 2019), MWA<fnref n="12" /> (Barry et al. 2019; Li et al. 2019; Trott et al. 2020; Nunhokee et al. 2025), and HERA<fnref n="13" /> (Abdurashidova et al. 2022b; Adams et al. 2023). Though there has not been a detection of the signal yet, these experiments have progressively set stronger upper limits on the 21-cm power spectrum, enabling some scenarios of reionization to be ruled out (Ghara et al. 2020, 2025; Greig et al. 2021; Abdurashidova et al. 2022a; Adams et al. 2023).

Detection of the 21-cm signal power spectrum from the Cosmic Dawn ($z > 15$) is even more challenging than the EoR, primarily due to the higher system noise of instruments at these low radio frequencies, contributed mainly by the significantly higher sky temperature. The astrophysical foreground emission is much stronger at these low radio frequencies, owing to the power-law spectral shape of synchrotron emission, imposing extremely high requirements for the accurate subtraction of foregrounds to reach the instrument's thermal noise sensitivity. Additionally, the ionosphere and wide-field nature of the instruments present stronger data processing challenges compared to higher frequencies. Instruments that have attempted detection of the Cosmic Dawn 21-cm signal power spectrum include the MWA (Ewall-Wice et al. 2016; Yoshiura et al. 2021), LOFAR (Gehlot et al. 2019), AARTFAAC<fnref n="14" /> (Gehlot et al. 2020), OVRO-LWA<fnref n="15" /> (Eastwood et al. 2019; Garsden et al. 2021), and NenuFAR<fnref n="16" /> (Munshi et al. 2024). Currently, these power spectrum upper limits are a few orders of magnitude above the levels expected from the standard 21-cm models of the Cosmic Dawn. However, the exotic non-standard models of the Cosmic Dawn that have emerged to explain the EDGES detection also predict the fluctuations of the 21-cm brightness temperature to be much stronger than those expected from standard models (Feng & Holder 2018a; Fialkov et al. 2018), effectively boosting the corresponding power spectrum. Still, the strongest power spectrum upper limits on the 21-cm signal during the Cosmic Dawn are at least an order of magnitude beyond that expected from the most extreme exotic models at $z > 15$. Deep power spectrum limits on the 21-cm signal during the Cosmic Dawn can thus observationally confirm or rule out scenarios of the exotic physics suggested by the EDGES detection through an independent experimental approach.

Besides the challenging S/N requirements, the strong foregrounds impose steep technical requirements in instrument calibration necessary for a detection. These challenges include contamination from radio frequency interference (RFI), ionospheric scintillation, modelling of far-field sources, and accurate estimation of the instrumental primary beam. In HERA, for example, cross-talk between feeds in the closely packed redundant array layout has proved to be one of the current limiting factors (Abdurashidova et al. 2022b; Adams et al. 2023), which is mitigated to some extent using fringe-rate filtering (Kern et al. 2019, 2020; Charles et al. 2023, 2024; Garsden et al. 2024). LOFAR, MWA, and NenuFAR upper limits, on the other hand, suffer from excess noise above the theoretical noise sensitivity of the instrument (Mertens et al. 2020; Trott et al. 2020; Munshi et al. 2024). For LOFAR, inaccurate modelling of bright sources (Ceccotti et al. 2025b), the turbulent ionosphere (Brackenhoff et al. 2024), and incomplete sky models used in calibration (Höfer et al. 2025)

<footnote n="1">¹Experiment to Detect the Global EoR Signature</footnote>

<footnote n="2">²Shaped Antenna measurement of the background RAdio Spectrum</footnote>

<footnote n="3">³Broadband Instrument for Global HydrOgen ReioNisation Signal</footnote>

<footnote n="4">⁴Probing Radio Intensity at High-Z from Marion</footnote>

<footnote n="5">⁵Radio Experiment for the Analysis of Cosmic Hydrogen</footnote>

<footnote n="6">⁶Mapper of the IGM spin temperature</footnote>

<footnote n="7">⁷Remote HI eNvironment Observer</footnote>

<footnote n="8">⁸Square Kilometre Array</footnote>

<footnote n="9">⁹Giant Metrewave Radio Telescope</footnote>

<footnote n="10">¹⁰Low-Frequency Array</footnote>

<footnote n="11">¹¹Precision Array to Probe EoR</footnote>

<footnote n="12">¹²Murchison Widefield Array</footnote>

<footnote n="13">¹³Hydrogen Epoch of Reionization Array</footnote>

<footnote n="14">¹⁴Amsterdam ASTRON Radio Transients Facility And Analysis Center</footnote>

<footnote n="15">¹⁵Owens Valley Radio Observatory – Long Wavelength Array</footnote>

<footnote n="16">¹⁶New Extension in Nançay Upgrading LOFAR</footnote>
have been ruled out as the main contributors to the excess variance. Transient RFI from aeroplanes has also been shown to currently have a negligible contribution to the power spectra estimated using the widefield instrument AARTFAAC (Gehlot et al. 2024). As the calibration and analysis pipelines of the instruments continue to improve, several recent studies (Gan et al. 2022; Chokshi et al. 2024; Munshi et al. 2024; Brackenhoff et al. 2025) point to the inaccurate modelling of the instrumental primary beam near nulls and the leakage of power from bright sources via their point spread function (PSF) sidelobes as one of the major limiting factors in reaching the thermal noise limit and exploiting the full sensitivity of instruments.

NenuFAR, a new low-frequency radio interferometer in France, is one of the most sensitive instruments to the Cosmic Dawn 21-cm signal fluctuations at large scales, with the sensitivity to detect the predicted exotic model power spectra in about 100 observing hours. The first analysis of NenuFAR data obtained on the north celestial pole (NCP) field by Munshi et al. (2024), hereafter M24, produced upper limits on the 21-cm signal power spectrum at $z = 20.3$, the best limit being $\Delta_{21}^2 < 2.4 \times 10^7 \text{ mK}^2$ at $k = 0.041 \text{ h cMpc}^{-1}$. However, these limits were over two orders of magnitude above the thermal noise power spectrum due to excess variance contributed primarily by distant off-axis sky sources and local sources of RFI. Through a combination of analytical modelling and forward simulations, two separate studies by Munshi et al. (2025a) and Munshi et al. (2025b) demonstrated the impact of wide-field sources and local RFI sources, respectively, on the power spectrum estimated with phase tracking instruments and offered strategies to mitigate their impact.

In this paper, we analyse multiple nights of NenuFAR observations of an optimal target field at two redshift bins centred at $z = 17.0$ and $z = 20.3$, and implement improvements to the analysis pipeline that decrease the excess variance in the data caused by these two contaminating effects. We report $2\sigma$ upper limits on the 21-cm signal power spectrum at both redshift bins and compare the resulting limits against the power spectra predicted by exotic models of the Cosmic Dawn 21-cm signal. The paper is organized as follows: Section 2 introduces the instrument and the observation strategy. Section 3 describes the pre-processing and calibration-based sky-model subtraction. The power spectrum estimation is described in Section 4, and the residual foreground removal and robustness tests are described in Section 5. We present our results in Section 6 and discuss their implications and future prospects in Section 7. We summarize the main conclusions from this analysis in Section 8. In this paper, we have used a flat Lambda cold dark matter cosmology, with the cosmological parameters ($H_0 = 67.7, \Omega_m = 0.307$) obtained by Planck Collaboration XIII (2016).

# 2 NENUFAR COSMIC DAWN OBSERVATIONS

One of the main drivers of NenuFAR is its Cosmic Dawn Key Science Program (CD KSP), which aims to detect the redshifted 21-cm signal power spectrum from the Cosmic Dawn. NenuFAR operates at very low frequencies (10–85 MHz), and the CD KSP focuses on the 40–85 MHz frequency range corresponding to a redshift range of $z = 15–31$ (Mertens, Semelin & Koopmans 2021).

## 2.1 Instrument and observation status

NenuFAR is composed of inverted V-shaped dual-polarized antennas that are grouped into analogue steerable stations composed of 19 antennas in a hexagonal configuration, called Mini-Arrays (MAs

<figure><img src="image_1.png" /><figcaption>**Figure 1.** The baseline coverage of observations of the NT04 field. The different colours correspond to the four different nights of observation.</figcaption></figure>

hereafter). Most of the MAs are distributed within a 400 m diameter dense core, with a few remote MAs located at distances of a few kilometres from the core. We refer the reader to Zarka et al. (2012, 2015, 2020) for detailed descriptions of the instrument. The observations analysed in this paper used a maximum of 78 core MAs (out of the planned 96) and four remote MAs (out of the planned eight). The defining feature of NenuFAR for Cosmic Dawn studies is its extremely dense core, which provides a very high collecting area and filling factor. Fig. 1 shows the baseline coverage of the four nights of observation used in this analysis and illustrates the dense sampling at small uv by the baselines formed by pairs of core stations. Owing to this large collecting area, NenuFAR is one of the most sensitive instruments to the fluctuations in the 21-cm signal at these frequencies and the most sensitive instrument below 50 MHz, where both HERA and the upcoming SKA will not observe. The KSP has already garnered over 2000 h of observations across multiple target fields, focusing on two primary deep fields: the NCP field and the NT04 (RA = 7 h 20 min, Dec. = 35°) field.

## 2.2 Field and night selection

M24 showed that, after foreground subtraction, the residual data from NenuFAR observations of the NCP field are still dominated by a strong excess variance above the thermal noise power. This excess variance is contributed primarily by residuals from strong off-axis sources picked up by NenuFAR’s primary beam grating lobes and local sources of RFI near the NenuFAR core. Based on these results, we initiated a survey of six fields with NenuFAR (Mertens et al., in preparation), to identify optimal target fields where foreground contamination is minimized. The selection of the fields was performed based on two sets of criteria. The major criteria constitute the distance from bright A-team sources, the apparent flux from these sources as seen through the primary beam sidelobes, low emission from the Galactic plane, and phase centre transit close to the zenith. The minor selection criteria constitute the presence of a calibrator source within the field and year-round observability. The power spectra estimated
<figure><img src="image_2.png" /><figcaption>**Figure 2.** RFI flagging statistics at different stages of processing. Left: the flagged RFI percentage in the highest time and frequency resolution correlated data for the two spectral windows (rows) for the four nights (columns). Right: a scatter plot of NenuFAR MA locations showing the number of nights each MA has been identified and flagged (in Z20) based on bad solutions in the A-team subtraction step.</figcaption></figure>

<table><caption>**Table 1.** Observational parameters of the four nights of observation of the NT04 field used in this analysis.</caption><thead><tr><th rowspan="2">Id</th><th rowspan="2">Start date</th><th colspan="2">Calibrator</th><th colspan="2">Target</th><th colspan="2">RFI per cent (target)</th></tr><tr><th>Start (UTC)</th><th>Span (h)</th><th>Start (UTC)</th><th>Span (h)</th><th>Z20</th><th>Z17</th></tr></thead><tbody><tr><td>1</td><td>2023 Dec 8</td><td>21:00</td><td>0.5</td><td>21:30</td><td>9.5</td><td>1.6</td><td>10.7</td></tr><tr><td>2</td><td>2023 Dec 29</td><td>19:00</td><td>0.5</td><td>19:30</td><td>11.5</td><td>1.7</td><td>10.7</td></tr><tr><td>3</td><td>2024 Jan 29</td><td>17:00</td><td>0.5</td><td>17:30</td><td>10.5</td><td>1.6</td><td>9.9</td></tr><tr><td>4</td><td>2024 March 2</td><td>17:00</td><td>0.5</td><td>17:30</td><td>8.5</td><td>1.6</td><td>10.0</td></tr></tbody></table>

after subtraction of the brightest sources were used as a metric to identify the NT04 field as the optimal target field for 21-cm studies with NenuFAR. We discuss the field selection survey in detail in a forthcoming paper by Mertens et al. (in preparation). We have now initiated a deep survey of this field with NenuFAR, with over 600 h of observation approved over Cycles 3 to 6 of NenuFAR observations, out of which more than 500 h of observation have already been collected.

In this paper, we select four nights from observations of the NT04 field that were available from Cycle 3 of NenuFAR observations (2023 December to 2024 June) based on superior RFI statistics. The left panel of Fig. 2 shows the RFI statistics of the four selected nights. All four observations of the target field were preceded by a 30 min observation of the calibrator source Cassiopeia A (Cas A). In this analysis, we focus on the higher frequencies of 61–85 MHz, which we divide into two spectral windows: 61–72.5 MHz ($z = 20.3$, Z20 hereafter) and 73–84.5 MHz ($z = 17.0$, Z17 hereafter), which are analysed independently. This is done to minimize the light-cone effect caused by the evolution of the signal within the redshift range being analysed (Datta et al. 2012, 2014). The 0.5 MHz band between Z20 and Z17 is not used in order to avoid known RFI. The observational details of the four selected nights are summarized in Table 1. The baseline coverage of the four nights of observation shown in Fig. 1 illustrates how combining the four nights not only improves the sensitivity by sampling the same modes multiple times but also improves the sampling of the uv plane.

## 3 CALIBRATION AND SKY-MODEL SUBTRACTION

The data processing pipeline employed by the NenuFAR CD KSP can be divided into three main steps: pre-processing, calibration-based sky model subtraction, and power spectrum estimation (including residual foreground removal). The routine pre-processing step is performed in the nancep computing cluster at the Nançay radio observatory, while the remaining two steps of the processing are performed in the DAWN cluster in Groningen (Pandey et al. 2020). A workflow describing the sequence of steps in the full data processing pipeline is presented in Fig. 3. NenuFAR data are calibrated using DDECal (Gan et al. 2023), and we currently use the diagonal mode in which only the diagonal elements of the Jones matrices are estimated. The calibration parameters used at the different stages of the pipeline are summarized in Table 2.

### 3.1 Pre-processing

The correlated data of NenuFAR are available at a time resolution of 1 s, divided into sub-bands of width 195.3 kHz each, using a polyphase filter bank. Each sub-band is further channelized into 64 channels at the raw data stage, also called the L0 level. The data for both the target and calibrator fields for each night, at the highest time and frequency resolution, are flagged using AOFlagger (Offringa, Van De Gronde & Roerdink 2012) to remove intermittent strong RFI. The results of this RFI flagging step are shown in the left panel of Fig. 2. While Z20 is relatively clean, Z17 suffers from strong known
<figure><img src="image_3.png" /><figcaption>**Figure 3.** A flowchart describing the NenuFAR 21-cm cosmology processing pipeline. The green rectangles indicate the different data levels, and the blue regions correspond to the three main sections of the pipeline.</figcaption></figure>

<table><caption>**Table 2.** Calibration parameters used at the different stages of processing.</caption><thead><tr><th rowspan="2">Stage</th><th rowspan="2">Direction(s)</th><th colspan="2">Solution interval</th><th rowspan="2">Smoothness</th></tr><tr><th>Time</th><th>Frequency</th></tr></thead><tbody><tr><td>Bandpass calibration</td><td>Cas A</td><td>4 s</td><td>15.3 kHz</td><td>None</td></tr><tr><td>A-team subtraction</td><td>A-team, NT04</td><td>4 s–2 min (adaptive)</td><td>195.3 kHz</td><td>2 MHz</td></tr><tr><td>3C subtraction</td><td>3C, NT04</td><td>4 min</td><td>195.3 kHz</td><td>4 MHz</td></tr><tr><td>Target subtraction</td><td>NT04 clusters</td><td>6 min<sup>a</sup></td><td>195.3 kHz</td><td>6 MHz</td></tr></tbody></table>

Note. <sup>a</sup>For the first night, this value is 8 min, matching the duration of the time segments.

RFI in the middle of the band, which results in nearly 100 per cent flagging over a range of these frequencies. Next, two channels of width 15.3 kHz each are flagged at either edge of each 195.3 kHz sub-band to avoid edge effects of the polyphase filter. The data are then averaged into 4 s time resolution and 12 channels per sub-band frequency resolution, and compressed with *Dysco* (Offringa 2016; Chege et al. 2024) to produce the L1 visibilities.

The cables connecting the MAs in NenuFAR vary in length between 20 and 150 m, making it crucial to account for cable reflections using a bandpass calibration step to avoid artefacts at high delay, which can be limiting for 21-cm cosmology studies. The bandpass solutions for each MA are calculated from the calibrator observation with *DDECAL* (Table 2 lists the calibration parameters),

and the solutions are transferred to the target field observations.<fnref n="17" /> Since the difference in the beam response between the calibrator and the target direction is not corrected here, the data do not have the correct flux scale yet, and this is corrected only after the A-team subtraction and DI correction step (Section 3.2.1). The bandpass calibrated data, after an additional round of RFI flagging with *AOFlagger* and averaging into three channels per sub-band frequency resolution, form the L2 visibilities. The analysis in the rest of the paper is performed on the L2 data.

<footnote n="17">¹⁷See appendix B of M24 for an example of how this accounts for cable reflections and other direction-independent (DI) bandpass effects.</footnote>
## 3.2 Nenuflow

The faint nature of the 21-cm signal imposes stringent precision requirements on calibration algorithms that need to be tuned to the specific instrument so that they can account for the most prominent issues it suffers from. The first analysis by M24 presented such a calibration pipeline for NenuFAR, which consists of a series of steps that iteratively improve the sky model of the target field while removing the contribution of strong foregrounds from the visibilities through direction-dependent (DD) calibration. In this work, we have structured the calibration and sky model subtraction part of the pipeline into a streamlined workflow, called Nenuflow.<fnref n="18" /> Nenuflow is built on Nextflow (Di Tommaso et al. 2017), which provides key benefits such as reliability, reproducibility, robust error handling, and streamlined parallelization. Additionally, an automated pipeline is very beneficial for systematically processing data from a large number of nights, which is expected to be needed in 21-cm cosmology analyses to detect the signal power spectrum. While the pre-processing and power spectrum estimation sections of the pipeline could be incorporated in Nenuflow, in this work, it is used only for calibration-based sky model subtraction since this is where the robust structure and parallelization offer the most benefits. The pre-processed L2 data of the target field for each night and spectral window are taken as input by Nenuflow.

### 3.2.1 Bright A-team source subtraction

The regular hexagonal layout of the 19 antennas making up the MAs in NenuFAR gives rise to strong grating lobes in their primary beams.<fnref n="19" /> The residual power from bright off-axis radio sources passing through the grating lobes was determined by M24 to be one of the primary causes of the strong excess power in NenuFAR data. So, a major portion of our efforts to improve the calibration pipeline must go into the bright source subtraction step. The strongest sources in the radio sky (e.g. Cygnus A, Cas A, Taurus A) are collectively known as the A-team. In the case of NenuFAR, due to its wide-field nature and the very high flux densities of these sources, the contaminating effects of the A-team on the data are particularly strong. The PSF sidelobes of these sources can dominate the target field, making it necessary to remove their contribution from the visibilities before the gain of the instrument in the target field direction can be corrected. We achieve this through a DD calibration step using DDECal, which provides calibration solutions for each MA separately in the direction of each A-team source and the target field (as a single direction), as a function of time and frequency. Currently, the NenuFAR beam model is not integrated in EveryBeam,<fnref n="20" /> which is used by DDECal. Thus, the model visibilities, against which the calibration is performed, do not have an imprint of the NenuFAR primary beam. As a result, the calculated gains primarily capture the effect of the NenuFAR primary beam for each MA. For each A-team source, the predicted visibility contribution, corrupted by the corresponding calculated gains, is subtracted from the visibilities. This is followed by correction of the residual visibilities using the gains in the target direction. This essentially inverts and applies the station-based gain solutions in a single direction (towards the target field) to the visibilities, in a DI gain correction step.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Apparent fluxes of A-team and 3C sources as a function of LST for Z20, during Night 1, estimated using a simulated NenuFAR primary beam model. The top and bottom panels show the apparent fluxes for the A-team and 3C sources, respectively. The yellow shaded regions indicate the range of apparent fluxes above the flux threshold used to select sources to include in the DD calibration step for each segment. The vertical lines demarcate the different time segments.</figcaption></figure>

The analysis by M24 divided the data into 52 min segments in order to parallelize calibration. Here, we refine this approach and divide the data into 12 min segments,<fnref n="21" /> such that the solution time intervals for all calibration steps are divisors of the time segment duration, avoiding residual narrower solution time intervals at the end of each segment. Additionally, the NenuFAR analogue beamformer repoints every 6 min, which is naturally accommodated by the 12 min calibration solution intervals. Apart from increasing the number of threads in parallelization, increasing the number of segments allows us to calibrate against an apparent sky model for each segment, where the simulated NenuFAR primary beam, calculated using nenupy<fnref n="22" /> (Loh et al. 2023) and averaged over time, frequency, and MAs, is applied to the intrinsic sky model before calibration. This provides a reasonable starting point for the calibration solutions, which need to capture the beam on time-scales shorter than a 12 min segment and across the 11.5 MHz frequency bandwidth. The top panel in Fig. 4 shows the A-team's apparent fluxes calculated using the simulated NenuFAR primary beam as a function of time for Z20 of Night 1. The ripples in the apparent fluxes with time are produced by the sources passing through the primary beam sidelobes. We impose a flux cut of 2 Jy such that a source is included in the calibration run for a segment

<footnote n="18">¹⁸https://github.com/kariukic/nenuflow</footnote>

<footnote n="19">¹⁹The different NenuFAR MAs are rotated by multiples of 10° with respect to each other, resulting in six different non-redundant sets of primary beam imprints due to the hexagonal symmetry.</footnote>

<footnote n="20">²⁰https://git.astron.nl/RD/EveryBeam</footnote>

<footnote n="21">²¹Only for Night 1 of Z20, which was processed first, we used 8 min segments instead.</footnote>

<footnote n="22">²²Nenupy (https://github.com/AlanLoh/nenupy) uses a dipole beam model obtained from electromagnetic simulations to compute the full NenuFAR primary beam model using the array factor of each MA.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Comparison between calibration with a fixed and adaptive solution time interval. The top row shows the XX gain amplitudes for an example station, in the direction of Cyg A, during a 1 h portion of Night 1. The middle row shows a zoomed-in view of a single 8 min time segment, corresponding to the dashed rectangular region of the top row. The bottom row shows the standard deviation of the real part of the visibilities for the full night, after A-team subtraction and DI correction, for both XX and YY correlations.</figcaption></figure>

only if it has an apparent flux above this threshold (yellow shaded region). This is slightly different to the approach adopted by M24, where all A-team sources located above the horizon were included. The flux cut is an input parameter to Nenuflow, allowing multiple calibration runs to ensure calibration solutions have sufficient S/N and describe the primary beam well.

Another improvement we make to the A-team source subtraction compared to M24 is the implementation of adaptive solution time intervals. In the DD calibration run for any given segment, the apparent flux for the different directions can be very different, as seen in Fig. 4. A direction with a lower apparent flux needs a longer solution interval to garner enough S/N. However, not accounting for the primary beam on such time-scales could lead to a large residual flux in a direction with high apparent flux. We address this by adapting the solution time interval as a function of direction. The number of sub-solutions needed per direction is calculated internally by Nenuflow, based on the estimated apparent fluxes for each segment. The directions with higher flux are allocated a higher number of sub-solutions per interval, and the higher total flux ensures that there is enough S/N within a sub-solution interval to obtain reliable solutions. The different A-team directions are allocated solution intervals between 4 s (the time resolution of L2 data) and 2 min, depending on their apparent flux in a 12 min time segment, as opposed to a fixed solution interval of 2 min used by M24. The target direction is solved with the full 2 min time interval, since the beam variations in the target direction are expected to be small on such time-scales. The top and middle rows of Fig. 5 illustrate

the effect of using adaptive solution intervals for an example MA by comparing the gain amplitude solutions in the direction of Cygnus A (Cyg A) for the XX visibility. The top panels show the solutions for a 1-h duration during Night 1, while the middle panels give a zoomed-in view of a single 8 min time segment. We find that Cyg A passes through multiple primary beam sidelobes and nulls, which are captured considerably better through the adaptive solution interval approach. Within the highlighted time segment, Cyg A direction is allocated the smallest possible (4 s) solution interval, which can better capture the nulls in both time (around 3.86 h) and frequency (around 65 MHz). The bottom row of Fig. 5 compares the standard deviation of the residual visibilities of Night 1 after A-team source subtraction and DI correction using the two approaches. Using an adaptive interval reduces the residuals significantly, especially during the middle of the observation when the phase centre is closer to the zenith and the sensitivity is maximum.

In total, we include five A-team sources, Cas A, Taurus A (Tau A), Virgo A (Vir A), Hydra A (Hyd A), and Cyg A, in the A-team subtraction step. The sky model for the target field is initially an intrinsic catalogue model obtained from the global sky model<sup>23</sup> in a 15-deg radius around the NT04 field phase centre. This model is composed of compact sources modelled using delta functions and Gaussians, and was created using data from The GMRT Sky Survey (TGSS; Intema et al. 2017), NRAO VLA Sky Survey (NVSS; Condon et al. 1998), WEsterbork Northern Sky Survey (WENSS; Rengelink et al. 1997), and VLA Low-frequency Sky Survey (VLSS; Cohen et al. 2007). M24 used a 20λ baseline cut in calibration and showed, through signal injection tests, that this induces signal suppression of ≈3 per cent in the A-team subtraction step due to overfitting in the baselines used in calibration. Here, we avoid this suppression by performing the calibration and the final power spectrum estimation with two mutually exclusive sets of baselines (Patil et al. 2016, 2017; Mouri Sardarabadi & Koopmans 2019; Mevius et al. 2022). This is possible due to the presence of an extra remote MA compared to M24, which forms a large number of additional baselines with the core stations. Only baselines longer than 40λ are used in calibration, and the resulting station-based gains are applied to the predicted models and subtracted from all visibilities. We note that this baseline cut, selected by M24, is smaller than the 250λ cut used by LOFAR, since the baselines in NenuFAR are shorter than LOFAR, and a higher baseline cut would drastically reduce the number of baselines available during calibration. After A-team subtraction and DI correction of the target field, an apparent sky model (effectively beam-averaged over the full night) of the target field is constructed by performing multifrequency synthesis imaging and deconvolution of the visibilities from all time segments within a night using WSClean (Offringa et al. 2014). The spectral information in the model is captured by fitting a second-order polynomial along frequency for each clean component during deconvolution. The sources outside a fixed radius (11° for Z20 and 9° for Z17) from the phase centre are filtered out to only keep genuine sources within the main lobe of the primary beam. Fig. 6 shows the extent of the sky model after this filtering, which remains confined within the first null of the primary beam, as indicated by the contours corresponding to 4 per cent of the peak primary beam amplitude. The A-team subtraction, DI correction, and imaging are then repeated using this improved sky model on the original (L2) visibilities, following the approach taken by M24, to produce an updated model of the target field. Fig. 6 shows the cleaned images constructed from

<footnote>23https://lcs165.lofar.eu/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Cleaned images of the target field after DI correction, with the sky models from Night 1 overplotted. The two panels correspond to the two spectral windows. The overplotted white curves indicate the clusters of sources, which are used as separate directions when the target field subtraction is performed in a DD calibration and subtraction step.</figcaption></figure>

all four nights, with the black circles indicating the location of the clean components of the updated sky model for Night 1.

The A-team source subtraction and DI correction are performed in a strongly parallelized fashion, over 15 nodes of the DAWN cluster, with up to five calibration runs performed in parallel on each node. This enables multiple runs that are necessary to identify and flag stations with bad calibration solutions over a series of iterations.<fnref n="24" /> The stations that are completely flagged in this step for each night are illustrated in the right panel of Fig. 2. This calibration solution-based station flagging results in the removal of 15 stations from all nights, and these were identified to have electronic malfunctions.<fnref n="25" /> A few more stations are flagged on individual nights based on artefacts in the calibration solutions. Fig. 7 shows wide-field dirty images constructed before and after sky-model subtraction for both spectral windows. Cyg A is more than 90 deg from the phase centre, so its contribution cannot be distinguished in these images. Cas A is the strongest source in both spectral windows, and for Z17, Hyd A is picked up more strongly by the grating lobes, which are closer to the phase centre compared to Z20. We find that the A-team sources are subtracted well, and their PSF sidelobe ripples that run across the entire sky are significantly reduced after the A-team source subtraction step.

### 3.2.2 3C source subtraction

NenuFAR's grating lobes pick up a large number of bright radio sources (top row of Fig. 7). M24 adopted a manual approach of identifying the sources that have the most significant contribution in the uv plane and performing a DD calibration and subtraction. Here, we use a more systematic approach, similar to what is used for the

A-team subtraction. A large number of 3C sources (Bennett & Simth 1962) visible near the primary beam grating lobes in the full-sky images are included in an intrinsic 3C source model. This is given as an input to Nenuflow, which applies the simulated primary beam for each segment and includes those sources above a minimum flux limit. The bottom panel of Fig. 4 shows the simulated apparent fluxes of the selected 3C sources for Z20 during Night 1. The sources that do not make the flux threshold for any segment are shown with faded colours. The effective S/N within a calibration solution interval for a certain direction depends primarily on the apparent flux received from this direction. However, if the DD calibration is performed for a large number of directions, the number of free parameters being solved for increases, effectively increasing the solver noise and raising the S/N requirement. We find that including all sources above a 2 Jy flux limit causes such an overfitting where, even though the contribution of the 3C sources is subtracted on the baselines used in calibration (> 40λ), the shorter baselines used for power spectrum estimation exhibit an excess variance. We thus adopt a more conservative approach of imposing a much higher flux limit of 10 Jy, such that fewer than six directions are included for any time segment. Thus, the source contributions are still subtracted when they have a very high apparent flux, and the subtraction is accurate and avoids overfitting. The choice of the flux cut, and hence the maximum number of directions to include in the 3C subtraction step, depends on the trade-off between the fraction of the flux that can be reliably subtracted, and the excess variance caused on the shorter baselines due to overfitting. The adaptive solution interval is not used in the 3C subtraction since the apparent flux of the strongest 3C sources is much lower than that of the strongest A-team sources, so the errors made within the solution time interval of 4 min are less significant. Comparing the top row of Fig. 7 with the bottom row, we see that while the 3C sources are not completely subtracted in this more conservative approach, the peak fluxes of the brightest sources decrease in the process.

<footnote n="24">²⁴Typically, fewer than three iterations are sufficient to identify and flag such unstable stations.</footnote>

<footnote n="25">²⁵This has since been fixed.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Wide-field dirty images before and after calibration-based sky model subtraction. The two columns correspond to the two spectral windows. The top and bottom rows show the images before and after sky model subtraction, respectively. The NenuFAR primary beam is simulated at a time resolution of 1 h for all four nights and averaged over time, frequency, and MAs for both spectral windows. The red contours indicate 4 per cent and 40 per cent levels of the primary beam.</figcaption></figure>

### 3.2.3 Target field subtraction

The apparent sky model of the target field, after filtering out sources outside the main lobe from the phase centre, is divided into multiple clusters (10 for Z20, eight for Z17) using a spatial *k*-means clustering algorithm (MacQueen 1967; Lloyd 1982). The choice of the number of clusters ensures sufficient S/N for the faintest cluster within the solution interval. The clusters for Night 1 for both spectral windows are shown in white convex hulls in Fig. 6. Each cluster forms a separate direction in DD calibration and subtraction, which is performed on the L3 visibilities in the same manner as M24, using a 6 MHz frequency smoothness constraint (see Table 2). This results in

the subtraction of sources till the confusion noise limit. The resulting full-sky dirty images are shown in the bottom row of Fig. 7. The maximum flux densities within the main lobe of the primary beam calculated from the residual dirty images are 0.74 Jy (for Z20) and 0.51 Jy (for Z17), which are close to the expected confusion noise levels of 0.70 Jy/PSF (for Z20) and 0.47 Jy/PSF (for Z17) estimated from equation (6) of van Haarlem et al. (2013). Comparing the same measured band in M24 to Z20, the residual peak flux for Z20 is lower than the value of 1.05 Jy reached by M24 due to a higher maximum baseline length of 3.3 km as opposed to that of 1.4 km for M24, resulting in a narrower PSF and thus higher angular resolution and lower confusion noise. The residual visibilities after target field subtraction
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Power spectra estimated at different data levels. The two rows correspond to the two spectral windows. The first column corresponds to the L2 data. The second and third columns show the level of power subtracted in the off-axis source subtraction and target subtraction, respectively. The fourth column shows the power spectrum after sky-model subtraction, and the fifth column shows its ratio with the thermal noise power spectrum. The black dashed lines indicate the full-sky horizon limit, and the white dashed lines indicate the source wedge of Cas A. The vertical dotted lines indicate the baseline cut. We note that the power spectra presented in this figure are estimated before the steps described in Sections 4.2 and 4.3 are applied to the data.</figcaption></figure>

are used to estimate the power spectrum after a final foreground removal step using Gaussian process regression (GPR; Section 5).

# 4 POWER SPECTRUM ESTIMATION

We estimate power spectra using *pspipe*<sup>26</sup> (Mertens et al. 2020, 2025). The visibility data are composed of multiple measurement sets for a large number of time segments for each night and spectral window. Gridded data cubes for each night and spectral window are produced by constructing dirty image cubes using *wsclean* (Offringa et al. 2014, 2019), converting into Kelvin units, applying a spatial Hann window, and spatial Fourier transforming into the *uvv* domain. The combined data cube for all four nights is obtained by computing a weighted average of the visibility cubes of the individual nights using the corresponding weight cubes obtained from the naturally weighted PSFs.

Though the final power spectrum estimation for setting the upper limits on the 21-cm signal is performed after residual foreground removal using GPR, we also estimate the power spectrum from the visibility data at each stage of calibration and foreground subtraction for diagnostic purposes. The visibility cubes are Fourier transformed along the frequency direction with a Blackman–Harris window, modulus squared, and averaged in spherical and cylindrical shells to estimate the spherical and cylindrical power spectra, respectively. We refer the reader to M24 for more details on the power spectrum estimation. For Z17, a large number of frequency channels in the middle of the frequency band are flagged (left panel of Fig. 2). Thus, as long as there is significant power in the lowest $k_\parallel$ modes, the Fourier transformation along the frequency axis causes strong

artefacts in the cylindrical power spectrum. Thus, the power spectra before GPR-based foreground removal for Z17 are computed only in the frequency range 79 to 85 MHz, for illustration purposes.

## 4.1 Residual power spectra

Before constructing image cubes and power spectra, we perform an additional flagging step based on the standard deviation of the residual visibilities per time segment, which is discussed in Appendix A. We then compute the cylindrical power spectra from the Stokes *I* visibilities at the different data levels L2, L3, and L4. The thermal noise power spectrum is estimated from the time-differenced Stokes *V* visibilities. Fig. 8 shows the power spectra at different data levels for both spectral windows, and illustrates the $k_\perp$, $k_\parallel$ modes at which the power is subtracted by computing the ratio of the power spectra from the combined data cube at successive data levels. The black dashed lines indicate the full-sky horizon limit of the foreground wedge (Munshi et al. 2025a) for the four nights. The white lines indicate the source wedge of Cas A, the region where the peak power of Cas A is expected to dominate (Munshi et al. 2025a). The vertical dotted line indicates the 40$\lambda$ baseline cut. Only $k_\perp$ modes below this threshold are used to set power spectrum upper limits, as they are unaffected by signal suppression during calibration. For both spectral windows, we find that off-axis A-team and 3C source subtraction (L2/L3) removes over two orders of magnitude of power, particularly in the $k_\perp$, $k_\parallel$ modes occupied by Cas A. We note that the slight increase in power in the EoR window between L2 and L3 data for Z17 is caused by time steps near the start and end of observations for each night, with higher RFI contamination (Fig. A1). These time steps are eventually flagged in the next two sections (Sections 4.2 and 4.3). The target field subtraction (L3/L4), as expected, removes power at the lowest $k_\parallel$ modes corresponding

<footnote n="26"><sup>26</sup>https://gitlab.com/flomertens/pspipe</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** The signature of RFI in the data. The two columns correspond to the two spectral windows. The top row shows the frequency channel differenced noise in the *uv* plane. The grey dashed circles indicate the baseline cut. The bottom row shows the ratio of the power spectrum estimated before and after five *uv* cells in the vicinity of the *u* = 0 line (indicated by vertical dashed lines in the top row) are flagged. In the bottom row, the black dashed lines indicate the full-sky horizon limit, and the vertical dotted lines indicate the baseline cut.</figcaption></figure>

to the region near the phase centre. After sky model subtraction, the L4 data are still more than four orders of magnitude beyond the thermal noise power at the lowest $k_\parallel$ modes due to confusion-limited foreground sources. The high power at the lowest $k_\perp$, $k_\parallel$ modes is possibly caused by diffuse Galactic foreground emission. At high $k_\parallel$, the residual power approaches the thermal noise power for both spectral windows. However, the excess power above the thermal noise is significantly higher for Z17, possibly because of higher RFI contamination. For Z20, we find a region of high residual power at low $k_\perp$ and high $k_\parallel$ modes. Z17 also exhibits high residual power at the lowest $k_\perp$ modes across the entire $k_\parallel$ range.

## 4.2 RFI and bright source flagging in the uv plane

We estimate the frequency difference noise in the gridded data cube after sky model subtraction for each uv cell, to identify residual contributions to the visibilities coming from different directions. The noise in the uv plane, in units of system equivalent flux density, is shown in the top row of Fig. 9. For both spectral windows, we find a region of higher residual rms along the u = 0 line. Such a signature can be caused by stationary RFI with respect to the array and has been demonstrated through forward simulations of local RFI by Munshi et al. (2025b). Since the celestial poles are the only points in the sky whose delays are fixed with respect to the array during a synthesis observation, stationary RFI adds up coherently at the NCP. For phase

centres other than the NCP, this resembles an RFI-induced source towards the north of the image, which results in a region of higher power along the u = 0 line due to sampling in uv tracks and gridding in the uv plane (Munshi et al. 2025a).<fnref n="27" /> We see that not only is the peak RFI variance higher for Z17 than for Z20 due to higher RFI contamination, but the variance across the rest of the uv plane is also increased. We flag five uv cells centred on the u = 0 line for all v values before constructing the power spectrum (indicated by the black dashed lines in the top panels of Fig. 9). The bottom row of Fig. 9 shows the ratio of the power spectra estimated before and after such a flagging step. We find that the higher power at the small $k_\perp$ modes is reduced by more than an order of magnitude through this flagging in both the spectral windows. We note that such a reduction of RFI power by flagging a region of the uv plane is only possible for phase centres away from the NCP field, where the RFI power is spread throughout the uv plane (Munshi et al. 2025b). As a result, to decrease the contamination from RFI in the NCP field, M24 identified and flagged baselines strongly affected by RFI by inspecting their delay power spectra. A second (f fainter) line of higher variance in the top panels of Fig. 9, along the south-west-north-east direction, is likely caused by the residuals of Cas A. Five uv cells in the vicinity of this direction are also flagged, to suppress the contribution from Cas A residuals.

## 4.3 Optimal elevation range

Over the course of each night, the elevation of the NT04 field changes with time. The foreground power in the cylindrical power spectrum is confined to a region known as the foreground wedge. The boundary of the foreground wedge, called the horizon line, has conventionally been assumed to be independent of the observation parameters. However, accounting for curved sky effects, Munshi et al. (2025a) showed that the slope of the horizon line depends on the elevation of the phase centre. These full-sky horizon line equations suggest that the horizon line is steepest when the phase centre is at the lowest elevation. Thus, selecting only a subset of the LST range for power spectrum estimation could significantly reduce the foreground contamination, albeit at the cost of a loss of sensitivity due to a smaller integration time. To reduce the foreground contamination at high $k_\parallel$, we estimate the power spectrum including only the LST range where the elevation of the NT04 field is higher than 50° for all four nights. This is compared against the power spectrum estimated from the data for the full LST range in Fig. 10. The horizon delay line for both cases is plotted in black dashed lines. We find that confining the elevation range results in a reduction of power in the $k_\perp$, $k_\parallel$ modes between the two horizon lines for both spectral windows, with the reduction for Z20 being over an order of magnitude in some modes.<fnref n="28" /> We find that for Z17, this also reduces power above the full LST range horizon line. This power could possibly be caused by RFI sources with a non-smooth spectrum, during the LST range when the phase centre is below 50° elevation. The exact origin of the higher contamination on baselines with lengths close to the baseline cut is currently unknown. As expected, we find an increase in power in other modes caused by increased thermal noise power due to the smaller integration

<footnote n="27">²⁷A slight deviation from the *u* = 0 direction seen in Fig. 9 could be attributed to the *uv* tracks not being perfect circles (see appendix B of Munshi et al. 2025a).</footnote>

<footnote n="28">²⁸Since the horizon delay lines indicate the maximum extent of the wedge, this extra power need not be confined completely between the horizon lines. Specifically, this power corresponds to the modes occupied by the foreground sources when the target field is below 50° elevation.</footnote>
Full LST range / Elevation > 50°

<figure><img src="image_10.png" /><figcaption>**Figure 10.** The ratio of power spectra estimated from the full data to that estimated from the LST range where the phase centre elevation is above 50°. The two columns correspond to the two spectral windows. The black dashed lines indicate the full-sky horizon limit, corresponding to the two cases. The vertical dotted lines indicate the baseline cut.</figcaption></figure>

time. Still, the significant reduction in foreground power at high $k_\parallel$ modes far outweighs this minor reduction of sensitivity, especially since high $k_\parallel$ modes are the hardest to clean in later steps. We note that the NCP is at a fixed elevation, and the horizon line is fixed throughout the observation duration, making it impossible to perform such an optimization of the LST range to confine the foreground contamination to low $k_\parallel$.

# 5 RESIDUAL FOREGROUND REMOVAL

After the calibration-based sky-model subtraction, RFI flagging in the uv plane, and elevation-based data selection, the data are still dominated by unmodelled foreground emission that is several orders of magnitude stronger than the thermal noise level. In NenuFAR data, this is constituted primarily by foreground emission within the main lobe of the primary beam below the confusion noise limit and residuals of bright off-axis sources, which still contaminate the target field through their PSF sidelobes. These foregrounds can be disentangled from the 21-cm signal based on their spectral smoothness, compared to the intrinsic spectral fluctuations of the 21-cm signal along the line of sight. We perform GPR (Mertens, Ghosh & Koopmans 2018) on the combined gridded visibility cube in the uvv space to model and subtract such residual foreground contributions from the data.

## 5.1 Gaussian process regression

Following M24, we use machine-learning-enhanced GPR (ML-GPR; Acharya et al. 2024; Mertens, Bobin & Carucci 2024) for residual foreground removal. In the context of foreground removal in 21-cm cosmology, GPR models the gridded data in the uvv space (d) as a sum of Gaussian processes along frequency that describe foregrounds (f<sub>fg</sub>), the 21-cm signal (f<sub>21</sub>), and thermal noise (n). Often, an additional excess variance component (f<sub>ex</sub>) is needed to account for power at small frequency coherence scales that cannot be described by the foreground or the 21-cm signal components. The GPR model can then be written as

$$d = f_{\text{fg}}(\nu) + f_{\text{21}}(\nu) + f_{\text{ex}}(\nu) + n(\nu). \quad (1)$$

<table><caption>**Table 3.** The covariance model used in ML-GPR for both spectral windows. The prior columns indicate the prior bounds in log space as a fraction of data variance for all $\sigma^2$ parameters and in linear space in units of MHz for all l parameters. In the posterior column, parameters that are unconstrained within the prior range are indicated with a ‘-’, and for those that reach the lower prior boundary, only the upper limit is listed.</caption><thead><tr><th rowspan="2">Kernel</th><th rowspan="2">Parm.</th><th colspan="2">Z20</th><th colspan="2">Z17</th></tr><tr><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>K<sub>int</sub> (RBF)</td><td>σ<sup>2</sup></td><td>-0.5, 0.5</td><td>0.104<sup>+0.014</sup><sub>-0.016</sub></td><td>-0.5, 0.5</td><td>0.046<sup>+0.011</sup><sub>-0.013</sub></td></tr><tr><td></td><td>l</td><td>20,60</td><td>27.171<sup>+0.563</sup><sub>-0.681</sub></td><td>20,60</td><td>30.705<sup>+0.857</sup><sub>-0.721</sub></td></tr><tr><td>K<sub>mix</sub> (RBF)</td><td>σ<sup>2</sup></td><td>-2, -1</td><td>-1.657<sup>+0.003</sup><sub>-0.004</sub></td><td>-2, -1</td><td>-1.157<sup>+0.008</sup><sub>-0.009</sub></td></tr><tr><td></td><td>l</td><td>0.1, 1</td><td>0.503<sup>+0.001</sup><sub>-0.001</sub></td><td>0.1, 1</td><td>0.582<sup>+0.003</sup><sub>-0.003</sub></td></tr><tr><td>K<sub>21</sub> (VAE)</td><td>x<sub>1</sub></td><td>-4, 4</td><td>-</td><td>-4, 4</td><td>-</td></tr><tr><td></td><td>x<sub>2</sub></td><td>-4, 4</td><td>-</td><td>-4, 4</td><td>-</td></tr><tr><td></td><td>σ<sup>2</sup></td><td>-6, -1</td><td>&lt;-3.208</td><td>-6, -1</td><td>&lt;-2.749</td></tr><tr><td>K<sub>ex1</sub> (Matern 3/2)</td><td>σ<sup>2</sup></td><td>-5, -2</td><td>-4.586<sup>+0.048</sup><sub>-0.062</sub></td><td>-5, -2</td><td>-3.163<sup>+0.012</sup><sub>-0.014</sub></td></tr><tr><td></td><td>l</td><td>0.02, 0.1</td><td>0.045<sup>+0.004</sup><sub>-0.003</sub></td><td>0.02, 0.1</td><td>0.057<sup>+0.002</sup><sub>-0.002</sub></td></tr><tr><td>K<sub>ex2</sub> (Matern 3/2)</td><td>σ<sup>2</sup></td><td>-5, -2</td><td>-3.812<sup>+0.067</sup><sub>-0.050</sub></td><td>-4, -1</td><td>-2.159<sup>+0.088</sup><sub>-0.092</sub></td></tr><tr><td></td><td>l</td><td>0.1, 1</td><td>0.251<sup>+0.025</sup><sub>-0.020</sub></td><td>0.1, 1</td><td>0.561<sup>+0.045</sup><sub>-0.048</sub></td></tr></tbody></table>

Each component is characterized by a frequency-frequency covariance function that describes its spectral behaviour. The total covariance kernel of the data (K) is then a sum of the covariance kernels of the foreground (K<sub>fg</sub>), 21-cm signal (K<sub>21</sub>), noise (K<sub>n</sub>), and excess (K<sub>ex</sub>) components:

$$K = K_{\text{fg}} + K_{\text{21}} + K_{\text{ex}} + K_{\text{n}}. \quad (2)$$

The covariance kernels for the foreground and excess components are generated from analytical covariance functions that are parametrized by a set of hyperparameters, such as the frequency coherence scale and variance of the components. We use covariance functions from the Matern class, defined by

$$K_{\text{Matern}}(\Delta\nu) = \sigma^2 \frac{2^{1-\eta}}{\Gamma(\eta)} \left( \frac{\sqrt{2\eta}\Delta\nu}{l} \right)^{\eta} K_{\eta} \left( \frac{\sqrt{2\eta}\Delta\nu}{l} \right). \quad (3)$$

Here, $\Gamma$ is the Gamma function, $K_\eta$ is the modified Bessel function of the second kind, $\eta$ defines the nature of the Matern class function, $\Delta\nu$ is the absolute frequency separation, and $\sigma^2$ and $l$ are hyperparameters that describe the variance and frequency coherence scale of the component, respectively. For the 21-cm signal covariance kernel, however, instead of an analytical covariance function, a variational auto-encoder (VAE) kernel trained on simulations of the 21-cm signal is used in order to describe the signal of interest using physically motivated models. Table 3 summarizes the GPR covariance model. All $\sigma^2$ values are specified as fractions of the data variance in log space, and the $l$ values are in units of MHz. Next, we describe the covariance kernels that are used for the different components in ML-GPR.

**Foregrounds:** We divide the foreground component further into an intrinsic component f<sub>int</sub> with covariance K<sub>int</sub> and a mode-mixing component f<sub>mix</sub> with covariance K<sub>mix</sub>. The intrinsic component describes the power at the lowest $k_\parallel$, which is constituted by residual foregrounds within the field of view and has extremely large-frequency coherence scales. The mode-mixing component captures the off-axis foreground sources that occupy the rest of the foreground wedge created by the chromatic and non-ideal uv coverage. For both the intrinsic and mode-mixing foreground components, we use radial basis function (RBF) covariance kernels, which are obtained by setting $\eta = \infty$ in equation (3). RBF covariances have a rapid decay in power along $k_\parallel$ and thus give rise to frequency-smooth functions that
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Power spectrum shapes captured by the 21-cm covariance kernel in ML-GPR. The normalized power spectra of the 1000 21cmFAST simulations for both redshifts are shown in grey curves. The normalized power spectra generated by the trained VAEs at a set of uniformly spaced points in the latent space are overplotted in brown.</figcaption></figure>

are ideal for modelling foregrounds, enabling easier separation from the 21-cm signal and other components with spectral fluctuations. We distinguish between the intrinsic and mode-mixing foregrounds based on the priors on $l$, with the intrinsic component having priors at much higher $l$ values. In future, this will be replaced by a physically motivated foreground covariance kernel that naturally produces the mode-mixing covariance using the primary beam models and foreground source distribution functions.

**21-cm signal:** The 21-cm covariance kernels are VAEs trained on simulations of the 21-cm signal. The code 21cmFAST (Mesinger, Furlanetto & Cen 2011; Murray Steven et al. 2020) is used to perform simulations of the 21-cm signal at $z = 17$ and 20 by varying a set of astrophysical parameters. We reuse the trained kernel that was used by M24 for $z = 20$. For $z = 17$, we perform an additional 1000 simulations of the 21-cm signal, estimate the power spectra, and train a VAE kernel (see M24 for details about the simulations and training). We note that the VAE is trained to capture only the shape of the 21-cm power spectrum, described by the two latent space dimensions $x_1$ and $x_2$, while the variance is left as a free parameter optimized during GPR. This enables the trained VAE to describe not only the standard 21-cm signal models but also extreme and exotic models of the 21-cm signal that can have a significantly boosted power spectrum<fnref n="29" /> Fig. 11 shows the normalized power spectra obtained from the 21cmFAST simulations at both redshifts (in grey) and the power spectra generated from the trained VAEs at 25 uniformly spaced points in the two-dimensional latent space (in brown). We see that the VAE is able to capture the variety of power spectrum shapes at both redshifts.

**Excess:** To account for residual power in the data at smaller frequency coherence scales than the foregrounds, we need to include an excess component. The excess component is described by a Matern covariance kernel with $\eta = 3/2$, which gives it a covariance structure that decays much slower than foregrounds along the $k_\parallel$ direction, enabling it to capture rapidly fluctuating spectral structures that are still correlated in frequency. Here, we divide the excess further into two components: $f_{ex1}$ (with covariance $K_{ex1}$), which

captures very small coherence scale structures in the data and accounts for almost noise-like excess variance above the estimated thermal noise, and $f_{ex2}$ (with covariance $K_{ex2}$), which has a larger coherence scale, and accounts for residual power in the data possibly caused by low-level systematics. The priors on $l$ are chosen in such a way as to reflect the corresponding frequency coherence scales.

**Thermal noise:** The thermal noise covariance is a diagonal matrix, with its variance estimated from the time-differenced Stokes V visibility cube.

We apply GPR separately to the visibility cubes of the two spectral windows. The posterior distribution of the hyperparameters of the Gaussian process covariance model is estimated with Nested sampling using dynasty (Speagle 2020), and is shown in Fig. B1 of Appendix B. We find that all parameters are constrained within compact regions in the parameter space within the prior bounds, except the 21-cm component. Here, the parameters controlling the 21-cm power spectrum shape, $x_1$ and $x_2$, remain unconstrained, and the variance $\sigma^2$ reaches the lower prior boundary, indicating that the 21-cm signal is not detected. The length-scale of $K_{int}$ converges to values of 27.2 and 30.7 MHz for Z20 and Z17, respectively, and thus describes extremely frequency smooth components, with length-scales well above the bandwidth of 11.5 MHz. Small differences in the converged intrinsic foreground length-scales between the two spectral windows are thus expected, first because they describe different portions of the foreground source spectrum, and secondly because the primary beams are different. $K_{mix}$ and $K_{ex1}$ also converge to similar length-scales for both Z17 and Z20. We note, however, that $K_{ex2}$ converges to significantly different length-scales for Z20 and Z17, with its variance being a higher fraction of the data variance for Z17 compared to Z20 (Table 3). This is likely because $K_{ex2}$ describes the excess variance due to the RFI contamination, which can have different spectral signatures in Z17 and Z20. The higher excess variance in Z17 is also expected due to the significantly more severe RFI contamination in Z17 compared to Z20 (Fig. 2) and was evident in the cylindrical power spectra before residual foreground removal (Fig. 8).

For each sample in the posterior distribution, the predictive mean and covariance of the foreground component can now be estimated

$$ 
\begin{aligned}
\mathcal{E}(\mathbf{f}_{\text{fg}}) &= \mathbf{K}_{\text{fg}}\mathbf{K}^{-1}\mathbf{d}, \\
\text{cov}(\mathbf{f}_{\text{fg}}) &= \mathbf{K}_{\text{fg}} - \mathbf{K}_{\text{fg}}\mathbf{K}^{-1}\mathbf{K}_{\text{fg}}
\end{aligned}
\qquad (4) 
$$

A realization of the foreground component is then subtracted from the visibility cube to yield the residual data cube. Following the approach described in section 3.8.1 of Mertens et al. (2025), an ensemble of power spectra estimated from the residual data cubes corresponding to 1000 hyperparameters sampled from their posterior distribution is then used to estimate the median residual power spectrum with uncertainties. The cylindrical power spectra estimated before and after GPR for the two spectral windows are shown in Fig. 12, as ratios over the corresponding thermal noise power spectra. To minimize the artefacts in the cylindrical power spectra before GPR caused by the flagged channels in Z17, we interpolate the gridded data values at the flagged channels using the converged data covariance ($K$) following the description in section 3.8.2 of Mertens et al. (2025). This interpolation is performed only for visualization and is not used anywhere else, as post-GPR power spectra inherently suppress such artefacts due to the reduced dynamic range along $k_\parallel$. Additionally, here the Fourier transformation is performed using least squares spectral analysis (Vaníček 1969) to further suppress these artefacts. We find that the power at low $k_\parallel$ is reduced significantly through GPR, leaving excess power above the thermal noise level,

<footnote n="29">²⁹While this naturally accounts for excess radio backgrounds, if the power spectrum boost arises from effects that alter the 21-cm power spectrum shapes, the framework assumes that such shapes can still be described by the variety of shapes accommodated by the trained VAE.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Cylindrical power spectra, expressed as a ratio over the thermal noise power spectrum, before and after GPR-based foreground subtraction. The two rows correspond to the two spectral windows.</figcaption></figure>

which is much stronger in Z17 than Z20, likely caused by stronger RFI contamination, as discussed earlier. The residual power near $k_\parallel = 1.8 \text{ h cMpc}^{-1}$ (for Z20) and $k_\parallel = 2.0 \text{ h cMpc}^{-1}$ (for Z17) are artefacts caused by the polyphase filterbank used in channelization into 195.3 kHz sub-bands.

## 5.2 Signal injection tests

To verify the robustness of the foreground removal with GPR, we perform a suite of signal injection tests. The test consists of injecting the visibilities corresponding to the 21-cm signal into the gridded data cube before GPR, running GPR in the same manner as done on the data, and comparing the recovered 21-cm signal against the injected signal to determine the level of bias. For a detailed description of how the signal injection tests are performed, we refer the reader to M24. For each spectral window, the injection test is repeated for 25 different shapes of the signal sampled uniformly within the range $x_1, x_2 \in [-2, 2]$ from the 2D latent space of the VAE. For each shape, we repeat the test for four different strengths of the input 21-cm signal corresponding to 0.25, 0.5, 1, and 2 times the thermal noise variance. For each of these 100 injected signals, we use the z-score as a metric for the level of bias that is incurred from GPR. The recovered power spectrum, corresponding to 1000 different points from the posterior distribution of the GPR hyperparameters, is estimated for each injected signal, and the z-score at each $k$ bin is obtained by calculating the inverse CDF of the $p$-value under the assumption of a Gaussian error distribution. The z-score represents the number of

standard deviations by which the recovered signal overestimates or underestimates the injected signal, with values below -2 indicating a 21-cm signal absorption beyond 2σ error bars.

The results from all 100 injection tests for each spectral window are summarized in Fig. 13. The left panel shows the z-scores at each $k$ bin, separately for the different strengths of the injected signal. The spread in the z-score values is caused by the results from the different shapes of the injected signals. For $z = 20.3$, we find that there is a positive bias at high $k$ bins, and here, several recovered signals are more than 3σ above the injected signal. At the lower $k$ bins, the z-score distributions are relatively well centred around zero, as we would expect if the recovery is unbiased. There are a few $k$ bins for which the z-score goes below -2, but these outliers constitute 1.33 per cent out of all computed z-score values. While there is a positive bias also at large $k$ bins for $z = 17.0$, the z-scores are closer to zero compared to $z = 20.3$, indicating that the injected signals are better constrained within the recovered signal uncertainties. This is due to the larger error bars and the stronger excess component in $z = 17.0$, which can more easily absorb a portion of the 21-cm signal when the VAE kernel fails to do so. The right panel of Fig. 13 illustrates the dependence of the average z-score on the shape of the injected signal for each of the four different injected signal strengths. Flatter signal power spectrum shapes have a lower z-score on average compared to steeper power spectra, which are better recovered with GPR, while often incurring a positive bias. We also find that such a positive bias is stronger in injected signals with more power (also seen in the left panel for higher $k$ values). This occurs because stronger signals are recovered with greater precision, resulting in a smaller standard deviation and, consequently, a higher estimated z-score. None the less, since our analysis focuses on setting upper limits on the 21-cm power spectrum, this positive bias is not a limiting factor.

# 6 RESULTS

In this section, we discuss the main results from this analysis and set upper limits on the 21-cm signal power spectrum. We then investigate the possible origin of the excess variance seen in the residual data.

## 6.1 Comparison with the NCP field

A key distinction in our analysis is the shift from the NCP field used by M24 to the NT04 field, selected through a targeted survey of multiple candidate fields (Mertens et al. in preparation), and the improvements in reducing the excess variance that this enabled (Section 4). In Fig. 14, we compare the cylindrical power spectrum from Z20 after calibration and sky model subtraction to that of the NCP field at the same stage. Since the NT04 field integration is ~2.3 times longer than the NCP observation, to leave the difference in noise power out of the comparison, we subtract the corresponding thermal noise power spectrum estimated from time-differenced Stokes V data from both power spectra before plotting. The modes with negative power correspond to the noise-dominated regions within the EoR window and are indicated in white. The foreground power seems to be much better confined to low $k_\parallel$ modes in the NT04 field analysis, and the residual wedge-like power in the NT04 data is possibly caused by the large number of radio sources in the grating lobes of the NenuFAR primary beam. In the right-most panel of Fig. 14, the ratio of the noise-bias-subtracted power spectra for the NCP and the NT04 analyses is shown. We find a reduction in power by at least an order of magnitude in the entire EoR window, and by over two orders of magnitude for some $k_\perp, k_\parallel$.
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Results of the signal injection tests performed on GPR for each spectral window. The two rows correspond to the two redshift bins. The left panels show the distribution of the z-scores as a function of *k* values and injected signal variances. The black crosses with the numbers alongside them indicate the number of signals for which the recovery is beyond the 3σ limit. The right panel shows the mean z-score as a function of the shape and strength of the injected signal power spectrum.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 14.** Power spectrum of the Z20 data before GPR, compared against the results obtained from the same spectral window from observations of the NCP field by M24 at the same stage of analysis. The three panels show the noise-bias-subtracted cylindrical power spectra, corresponding to this (NT04) analysis, the NCP analysis, and the ratio of the NCP and NT04 power spectra.</figcaption></figure>

modes. This suggests that the improved choice of the field, combined with improvements in data processing methods and the ability to flag RFI contributions and confine the foreground wedge, significantly reduces the excess variance that was observed in the analysis by M24. The power in the foreground wedge remains at a roughly similar level for the two fields. The reduction in power at the lowest $k_\parallel$ modes in the NT04 power spectrum could be attributed to the better target subtraction using more clusters in this analysis, enabled by the higher angular resolution due to an additional remote MA being used during the NT04 observation. The factor of ≈ 2 reduction in the power in these modes matches the value expected from the confusion noise values reported in Section 3.

## 6.2 Power spectrum upper limits

After GPR-based foreground removal, the region outside the uv range of 10λ–40λ is filtered out from the visibility cube. The spherical

<figure><img src="image_15.png" /><figcaption>**Figure 15.** Upper limits on the 21-cm signal power spectrum. The results from this analysis on the NT04 field are shown in brown, while those from the NCP analysis by M24 are shown in blue. The two panels show the two redshift bins. The black dashed and dotted lines indicate the thermal noise power spectrum and 2σ sensitivity, respectively.</figcaption></figure>

power spectrum and uncertainties are computed in logarithmically spaced bins in *k* space, ranging from the smallest accessible *k* value to a maximum value of $k = 0.5 \, \text{h cMpc}^{-1}$. Next, the thermal noise power spectrum is subtracted to yield the noise-bias-subtracted residual power spectrum ($\Delta_{21}^2$). We set upper limits on the 21-cm signal at a 2σ level above the residual power spectrum at both redshifts ($\Delta_{\text{ul}}^2$). Unlike M24, a calibration bias correction step is not required here since the baselines used in calibration ($> 40\lambda$) are not used in estimating the upper limits. Fig. 15 shows the resulting noise-bias-subtracted power spectra. For $z = 20.3$, the corresponding limits obtained by M24 are also shown. The vertical error bars indicate the 2σ errors on $\Delta_{21}^2$, and the horizontal error bars correspond to the *k* bin range. The dashed and dotted lines are the estimated thermal noise power spectrum and the 2σ uncertainties on the thermal noise, respectively, the latter describing the sensitivity of the instrument and indicating the minimum value that the upper limits can have. We find that for $z = 20.3$, the data power spectrum is within an order of
<table><caption>**Table 4.** 2σ upper limits on the 21-cm signal power spectrum derived from this analysis.</caption><thead><tr><th rowspan="2">k<br>h cMpc<sup>-1</sup></th><th colspan="3">z = 20.3</th><th colspan="3">z = 17.0</th></tr><tr><th>Δ<sup>2</sup><sub>21</sub>(k)<br>[mK<sup>2</sup>]</th><th>Δ<sup>2</sup><sub>ul</sub>(k)<br>[mK<sup>2</sup>]</th><th>k<br>h cMpc<sup>-1</sup></th><th>Δ<sup>2</sup><sub>21</sub>(k)<br>[mK<sup>2</sup>]</th><th>Δ<sup>2</sup><sub>ul</sub>(k)<br>[mK<sup>2</sup>]</th></tr></thead><tbody><tr><td>0.038</td><td>2.1 × 10<sup>5</sup></td><td>4.6 × 10<sup>5</sup></td><td>0.041</td><td>2.5 × 10<sup>6</sup></td><td>5.0 × 10<sup>6</sup></td></tr><tr><td>0.066</td><td>9.0 × 10<sup>5</sup></td><td>1.8 × 10<sup>6</sup></td><td>0.071</td><td>1.2 × 10<sup>7</sup></td><td>2.3 × 10<sup>7</sup></td></tr><tr><td>0.107</td><td>3.2 × 10<sup>6</sup></td><td>5.6 × 10<sup>6</sup></td><td>0.104</td><td>3.0 × 10<sup>7</sup></td><td>5.8 × 10<sup>7</sup></td></tr><tr><td>0.166</td><td>9.3 × 10<sup>6</sup></td><td>1.5 × 10<sup>7</sup></td><td>0.170</td><td>8.3 × 10<sup>7</sup></td><td>1.4 × 10<sup>8</sup></td></tr><tr><td>0.267</td><td>2.6 × 10<sup>7</sup></td><td>3.7 × 10<sup>7</sup></td><td>0.271</td><td>1.4 × 10<sup>8</sup></td><td>2.1 × 10<sup>8</sup></td></tr><tr><td>0.408</td><td>4.0 × 10<sup>7</sup></td><td>4.9 × 10<sup>7</sup></td><td>0.405</td><td>1.8 × 10<sup>8</sup></td><td>2.5 × 10<sup>8</sup></td></tr></tbody></table>

magnitude of the thermal noise power at all k bins, while the upper limits are still up to two orders of magnitude above the sensitivity. The upper limits derived from this analysis are an improvement to those obtained by M24 at the same redshift bin by over a factor of 50. For z = 17.0, the data power spectrum is still more than an order of magnitude above the thermal noise power at all but the largest k modes.

Table 4 summarizes the upper limits on the 21-cm signal obtained from this analysis. At z = 17.0, the best 2σ limit on the 21-cm signal power spectrum is $\Delta_{21}^2 < 5.0 \times 10^6 \text{ mK}^2$ at $k = 0.041 \text{ h cMpc}^{-1}$. This limit suffers from excess variance, primarily caused by RFI contamination of NenuFAR observations in this spectral window. Still, this limit is considerably deeper than all Cosmic Dawn power spectrum limits set by other instruments. At z = 20.3, the best 2σ upper limit on the 21-cm signal power spectrum is $\Delta_{21}^2 < 4.6 \times 10^5 \text{ mK}^2$ at $k = 0.038 \text{ h cMpc}^{-1}$. This limit is deeper than all previous Cosmic Dawn power spectrum upper limits by more than an order of magnitude (Ewall-Wice et al. 2016; Eastwood et al. 2019; Gehlot et al. 2019, 2020; Garsden et al. 2021; Yoshiura et al. 2021; Munshi et al. 2024).

Fig. 16 presents a compilation of all 21-cm power spectrum upper limits from the Cosmic Dawn ($z > 15$). The 2σ sensitivities of the full NenuFAR array configuration at $k = 0.04 \text{ h cMpc}^{-1}$ for 10 and 1000 h of integration are simulated and indicated in the figure. The predicted range of power spectra of standard models is obtained from the Evolution of 21 cm Structure (EOS) simulation suite (Mesinger, Greig & Sobacchi 2016) and shown in blue. An approximate range for exotic models is also shown in grey. A more comprehensive exploration of the exotic model power spectra and their comparison against the upper limits derived in this analysis is presented in Section 7.1.

## 6.3 Coherence of the excess variance

In both spectral windows, we find excess variance that limits our ability to reach the theoretical sensitivity, with Z17 having a much stronger excess than Z20 due to RFI contamination. It is important to understand the nature of this excess variance in order to remove it or mitigate its impact and assess the future prospects of longer integrations of hundreds of hours of NenuFAR observations at both spectral windows. For this purpose, we use *pspipe* to compute gridded visibility cubes in brightness temperature units from the sky-model subtracted residual data of each 12 min time segment for all four nights. Next, these visibility cubes are combined by computing a weighted average using the corresponding weight cubes, taking an increasing number of cubes starting from the first cube of night 1, and the corresponding cylindrical power spectra are estimated. To understand the nature of the excess variance as a function of regions in the $k_\perp, k_\parallel$ space, we define three regions:

<figure><img src="image_16.png" /><figcaption>**Figure 16.** Cosmic dawn power spectrum upper limits from different interferometers. The range of power spectra expected from standard models (blue shaded region) and an approximate range of power spectra expected from exotic models (grey shaded region) are indicated. The k values mentioned in this figure are in units of $h \text{ cMpc}^{-1}$.</figcaption></figure>

(i) Main lobe: The region within an angular distance of 15° from the phase centre, which contains the main lobe of the NenuFAR primary beam. This is given by the region in the $k_\perp, k_\parallel$ space below the 15° delay line derived under a flat sky approximation (Morales et al. 2012).$^{30}$

(ii) Sidelobes: The region between the primary beam main lobe and the horizon. In the $k_\perp, k_\parallel$ space, this corresponds to the region between the primary beam delay line and the steepest full sky horizon delay line across the four nights.

(iii) EoR window: The region in the $k_\perp, k_\parallel$ space beyond the full-sky horizon delay line.

We note that the NenuFAR primary beam main lobe extends up to ≈ 11° from the phase centre for Z20 and ≈ 9° for Z17. Thus, the region within 15° from the phase centre that is chosen here also includes the first nulls and sidelobe. However, the finite bandwidth imposes a minimum $k_\parallel$ value that we can measure, preventing the use of smaller angular distances where enough modes can be averaged below the delay line. For the sidelobe and wedge regions, we impose a minimum $k_\parallel$ value of $0.2 \text{ h cMpc}^{-1}$ as a buffer, to better isolate the behaviour within the wedge and the EoR window from spectral leakage from the primary beam main lobe. For all regions, a maximum $k_\parallel$ cut-off at $1 \text{ h cMpc}^{-1}$ is imposed to avoid the artefact caused by the polyphase filter.

The estimated powers in all $k_\perp, k_\parallel$ bins within each region are averaged separately. The entire analysis is repeated for the time-differenced Stokes V data cubes, which give us an estimate of the thermal noise power in each region. The results are shown in the top row of Fig. 17. The noise power is consistent across the three regions and exhibits the expected inverse dependence on the integration time.

<footnote>$^{30}$Though this delay line equation is not strictly correct for phase centres far from zenith (Munshi et al. 2025a), they are still valid for such small angular distances where the flat-sky approximation holds ($\sin(\pi/12) \approx \pi/12$).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 17.** Coherence of the excess variance across LSTs and nights in different regions of the $k_{\perp}$, $k_{\parallel}$ space. The different columns correspond to three regions: the main lobe, the sidelobes, and the EoR window. The top row shows the average data and noise power in the cylindrical power spectrum as a function of integration time, with the vertical solid lines demarcating the four nights. The middle and bottom rows show the cross-coherence values calculated for the two spectral windows.</figcaption></figure>

In the main lobe, the reduction in data power with integration time in the beginning could be caused by the non-axisymmetric primary beam main lobe picking up different sources as it rotates against the sky. The confusion noise-dominated power then adds up coherently with increasing integration time. In the sidelobes, we find that the power adds up relatively incoherently within a night, which was also observed by M24. This is a combination of two factors: noise dominating in the first few segments and the grating lobes picking up different groups of sources as the sky rotates. For the subsequent nights, the overall sidelobe data power stops decreasing further since the power from off-axis sources in the grating lobes at the same LST for the different nights still adds up coherently. The slight increase in the data power at the beginning of each night is possibly caused by the stronger foreground contamination in this specified region of the

$k_{\perp}$, $k_{\parallel}$ space due to steeper horizon lines corresponding to the low elevation of the phase centre.

In the EoR window, we find that for Z20, the power decreases monotonically with integration time, though there is an excess above the thermal noise power. However, for Z17, after the first night of integration, the data power starts to average down significantly slower than incoherent noise, indicating that a portion of the excess variance is coherent and starts to dominate the noise power. Residual power in the EoR window is a limiting factor for 21-cm cosmology analyses since the foreground subtraction algorithms cannot easily remove this power, without also suppressing the 21-cm signal. Thus, unless improvements are made to the processing strategies that can mitigate this coherent excess in Z17, possibly caused by RFI contamination, deeper integrations are unlikely to improve significantly upon the
current upper limits. However, Z20 looks promising since the excess power seems to integrate down incoherently, roughly like thermal noise, at least within these four nights of observation.

We can investigate the nature of the excess variance in more detail through a cross-coherence analysis. The cylindrically averaged cross coherence between two data cubes $T_i$ and $T_j$ corresponding to the ith and jth segments is defined as

$$C_{i,j}(k_{\perp}, k_{\parallel}) = \frac{\langle |T_i^*(k)T_j(k)| \rangle^2}{\langle |T_i(k)|^2 \rangle \langle |T_j(k)|^2 \rangle}, \quad (5)$$

where $\langle \cdots \rangle$ represents a weighted average over all k values lying within the $k_{\perp}$, $k_{\parallel}$ bin. The cross-coherence is computed for all pairs of 12 min gridded data cubes across the four nights.

The computed cross coherence in all $k_{\perp}$, $k_{\parallel}$ bins within each region is averaged separately. The results for both spectral windows are shown in the middle and bottom rows of Fig. 17, where each pixel corresponds to one such average cross-coherence value. We note that the values of the diagonals in the cross-coherence between the same nights are unity by definition. For the cross-coherence between different nights, the higher values near the diagonal indicate that similar LSTs are more coherent across nights due to the apparent sky being identical and the PSF being similar across nights for the same LST values. The values of cross-coherence at the same LST for different nights are highest within the main lobe and decrease as we go to the sidelobes, ultimately reducing to less than 10 per cent in the EoR window. Away from the diagonal, the coherence decreases while still remaining above 30 per cent within the main lobe for Z20. The coherence in the main lobe is lower for Z17, possibly because of the narrower primary beam, resulting in a higher decoherence due to thermal noise within the chosen 15° angular radius. Additionally, only the frequency range above 79 MHz is used for Z17 to avoid artefacts due to frequency gaps. This reduces the resolution along $k_{\parallel}$ and, consequently, the number of $k_{\parallel}$ bins that lie below the delay line for which the computed coherence values are averaged, increasing the sampling variance in the computed coherence. In the sidelobes, away from the diagonals, there is a higher coherence at the beginning and the end of the observations for both spectral windows. This is possibly caused by the fact that the horizon line is steeper during these periods because the phase centre is at a lower elevation. Thus, the power from off-axis sources can add up more coherently in this region of the $k_{\perp}$, $k_{\parallel}$ space. The residual coherence in the EoR window for both spectral windows can be attributed to leakage from foregrounds or low-level RFI that is partially coherent at the same LST across nights. This coherence is lower for nights that are more separated in time, as we move towards the bottom right square in each panel. We note that the residual coherence away from the diagonal in the EoR window is higher for Z17 for several pairs of nights. A possible origin could be the contaminating effects of RFI, which also result in the much stronger excess in Z17. The coherence analysis supports the insights from the analysis of how the power integrates down with time, in that Z20 is significantly cleaner, with the excess variance being more incoherent with time compared to Z17.

## 7 DISCUSSION

In this section, we discuss the astrophysical relevance of the upper limits derived in this analysis by comparing them against simulated models. We also highlight the main limitations of the current analysis and discuss prospects for future improvement.

## 7.1 Astrophysical interpretation

We focus only on the upper limits at $z = 20.3$, which are the more stringent limits from this analysis. Since these upper limits remain above the predictions of most models, including exotic ones, we do not perform a full inference analysis. None the less, it is instructive to illustrate how much these limits must be reduced in order to place meaningful constraints on theoretical models. This is what is shown in Fig. 18.

The left panel shows models computed with the simulation code 21-cm Semi-numerical Predictions Across Cosmological Epochs (21cmSPACE; e.g. Fialkov et al. 2012, 2013; Visbal et al. 2012; Fialkov, Barkana & Visbal 2014),<fnref n="31" /> which include the contribution from radio-emitting high-redshift galaxies. Such models (Feng & Holder 2018b; Mirocha & Furlanetto 2019; Reis, Fialkov & Barkana 2020; Sikder et al. 2024a) were originally inspired by the tentative detection of the sky-averaged (global) 21-cm signal from cosmic dawn (Bowman et al. 2018) by the EDGES experiment. Here, we use the latest version of the code as described by Liu et al. (2025).<fnref n="32" /> Since the astrophysical parameters of high-redshift galaxies are highly uncertain, we fix the star formation efficiency of Pop II stars to $f_{\star,II} = 0.5$ (not unlikely for bursty or dense compact star-forming galaxies at high-$z$, e.g. Dhandha et al. 2025; Yung, Somerville & Iyer 2025), and Pop III stars to $f_{\star,III} = 0.05$ with a fast recovery period after Pop III supernovae (Magg et al. 2022). In the plot, only the parameter to which the 21-cm signal at cosmic dawn is most sensitive, namely the X-ray efficiency $f_X$, is varied within the range $[10^{-3}, 10^3]$. We allow for low values still not completely ruled out (Dhandha et al. 2025) since they produce the highest 21-cm fluctuations. This parameter is normalized so that a value of unity corresponds to the ratio of X-ray luminosity to star formation rate ($L_X/SFR$) of local starburst galaxies with low-metallicity high-mass X-ray binaries (Fragos et al. 2013a, b). For each model, we fix the radio efficiency $f_R$ (i.e. the ratio of low-frequency radio synchrotron luminosity to star formation rate (SFR) normalized to low-redshift galaxies; Hardcastle et al. 2016; Gürkan et al. 2018) to be the maximum allowed according to the observed intensity of the extragalactic background, based on the measurements by ARCADE-2 (Fixsen et al. 2011) and LWA-1 (Dowell & Taylor 2018b). At $z = 20.3$, these constraints yield $f_R \sim 140 000$ following the methodology of Sikder, Barkana & Fialkov (2024b).

The right panel of Fig. 18 shows a set of models derived from the Loreli II simulations (Meriot, Semelin & Cornu 2025) obtained using the LICORICE radiative transfer code (Semelin, Combes & Baek 2007; Baek et al. 2009, 2010; Semelin 2016; Semelin et al. 2017). While the original Loreli II simulations are standard models of the early universe, a $\times 2000$ boost to the SFR was applied compared to the original simulations. This leads to a star formation rate density (SFRD) at $z \sim 20$ in the $2 \times 10^{-6} - 8 \times 10^{-3} M_{\odot} yr^{-1} cMpc^{-3}$ range, which matches those presented, for instance, by Cang et al. (2025). Then, similarly to 21cmSPACE models, high-redshift galaxies were considered the source of intense radio emissions, creating a fluctuating radio

<footnote n="31">³¹https://www.cosmicdawnlab.com/21cmSPACE</footnote>

<footnote n="32">³²We note that the cosmology assumed in 21cmSPACE is Planck 2013 (Planck Collaboration XVI 2014), while the rest of the work uses Planck 2015 parameters. Although important for precision inference, at the current level of the limits and given that the errors on seminumerical simulation approaches of $\sim 20$ per cent (Zahn et al. 2011; Hutter 2018), this should have a negligible impact on the conclusions.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 18.** Comparison of the upper limits derived in this analysis at $z = 20.3$ against simulated exotic 21-cm signals. The left panel shows the power spectra of signals simulated using 21cmSPACE and the right panel shows the power spectra obtained from the LICORICE simulations. The colourbar indicates the corresponding global 21-cm signal depth.</figcaption></figure>

background that dominates over the CMB. The radio background was computed using a modified version of the SPINTER code (Semelin et al. 2023), following Sikder et al. (2024a), and using a radio efficiency parameter of $f_r = 10^3$ or $f_r = 10^4$. We assumed negligible X-ray heating (i.e. $f_X = 0$) and ionization. For a given value of $f_r$, the different power spectra shown in Fig. 18 correspond to different pairs of minimum halo mass for star formation and gas conversion time-scale, spanning the full range described in (Meriot et al. 2025), and resulting in the range of SFRD mentioned above.

While our goal here is not to compare the two sets of models but rather to show a diversity of possible signals, we briefly comment on the observed differences between the 21cmSPACE and LICORICE models. While their power spectra exhibit similar slopes and span a similar range of amplitudes, the LICORICE models reach somewhat higher amplitudes, even surpassing the NT04 upper limits at large scales. A plausible explanation is that, while both sets of models reach a similar SFRD, the emissivity is not distributed in the same way. LICORICE models can only form stars in haloes more massive than $10^8 M_\odot$ (thus considering molecular cooling to be ineffective), while 21cmSPACE models allow for star formation in haloes down to $10^6 M_\odot$, although Lyman-Werner feedback is allowed to limit the efficiency in the least massive haloes. This results in a more diffuse emissivity in 21cmSPACE, for both the Wouthuysen-Field coupling and the radio background, and thus in a decreased power spectrum amplitude.

As we can see in Fig. 18, some models exceed the NT04 upper limits at large scales only when stringent conditions are met: a very strong radio emission by high-redshift galaxies, a high SFRD to generate sufficient Wouthuysen-Field coupling, and an efficient suppression of star formation in haloes below $10^8 M_\odot$ to generate a strongly fluctuating emissivity field. These models have a corresponding global signal in the 3000 mK range, much stronger than the EDGES signal, and at a higher redshift. An improvement of a factor of $\gtrsim 10$ in the upper limits will allow us to constrain models with a power spectrum amplitude similar to the EDGES signal.

## 7.2 Limitations and future work

As hundreds of hours of observation with NenuFAR are integrated in the future, there are a few key areas where the current processing pipeline could be improved.

**Local RFI sources:** The current approach to mitigating contamination from local RFI sources in NenuFAR utilizes the non-stationarity of the target field with respect to terrestrial RFI sources to flag their contribution in the *uv* plane. However, the RFI power continues to be present in the rest of the *uv* plane, but at a lower level due to suppression by the gridding kernel. Detailed modelling of the spectral and temporal behaviour of the local RFI sources in NenuFAR has already been investigated (Munshi et al. 2025b), and approaches to subtract these near-field models can be further developed. The NenuFAR team has also improved the shielding of one of the identified RFI sources, which has been seen to reduce the contamination in more recent NenuFAR observations.

**Off-axis source contamination:** While the current approach of a carefully chosen target field significantly improves the A-team contamination compared to previous analyses, the signature of the residuals from the strongest A-team source (Cas A) can still be seen in the gridded *uv* plane (Fig. 9). Additionally, the 3C subtraction in this analysis followed a conservative approach to avoid overfitting, which results in 3C sources in the primary beam grating lobes dominating the power spectra after sky model subtraction. The inclusion of a primary beam model in calibration could improve the off-axis source subtraction, where the calibration solutions do not need to describe the complete beam variations but only need to capture deviations of the simulated beam from the true beam. The solutions can thus be constrained better due to the lower number of free parameters. The NenuFAR MA primary beam is analogue steered, and inactive antennas within an MA could drastically change the primary beam pattern and impact gain calibration (Brackenhoff et al. 2025). Thus, it could be important to have prior knowledge of the broken antennas for each MA for effective inclusion of the primary beam in calibration.

**Faster processing mode:** Currently the computational time required by the entire calibration and sky-model subtraction pipeline
implemented in Nenuflow (Fig. 3) is about 12 h per night of
observation (with duration ≈ 10 h) when run in a parallelized manner
over 15 computational nodes of the DAWN cluster, running up to five
processes parallelly in each node (one per CPU). A large fraction of
this runtime (≈ 61 per cent) is spent on the imaging and sky model
iteration step. In the future, we are planning to use a deep intrinsic sky
model from these calibrated four nights of data, effectively reducing
the total runtime to ≈ 4.7 h per night. In practice, the processing
of a single night of NenuFAR observations takes longer than the
runtime of 12 h, since the first A-team subtraction step needs to be
repeated multiple times to identify faulty stations based on calibration
solutions (right panel of Fig. 2). This step could be automated in
the future by using outlier detection algorithms on the calibration
solutions.

Excess variance bias correction: The excess variance in the data above the thermal noise level has not been subtracted in this analysis, to set conservative upper limits on the 21-cm power spectrum. In forthcoming work, we will investigate the nature of the (in)coherence of the excess noise in greater detail, since making a bias correction for the incoherent component (as for the thermal noise) could lead to a considerably stronger upper limit.

8 CONCLUSIONS

In this paper, we have analysed four nights of observations of the
NT04 field (RA = 7 h 20 min, Dec. = 35°) with NenuFAR to set
upper limits on the 21-cm power spectrum from the Cosmic Dawn
at two redshift bins centred at *z* = 20.3 and *z* = 17.0. Our main
conclusions and insights are summarized below.

Choice of a new field: We have switched from the NCP field analysed by M24 to the NT04 target field, selected carefully through a survey of multiple fields. First, the new field selection is optimized to reduce PSF sidelobe leakage from A-team sources, one of the main contributors to the excess variance observed by M24. Secondly, switching to a field other than the NCP improves our ability to mitigate power from terrestrial RFI and strong off-axis sources. RFI manifests as a vertical line in the uv plane (Munshi et al. 2025b), which we mitigate to an extent by flagging the specific region in the uv plane. Off-axis sources are most problematic when the phase centre is at a low elevation due to the foreground wedge being steeper (Munshi et al. 2025a), and this foreground contamination at high *k*<sub>‖</sub> modes can be reduced by implementing an elevation threshold in the power spectrum estimation. Each of these steps decreases the contamination at the relevant *k*<sub>⊥</sub>, *k*<sub>‖</sub> modes by more than an order of magnitude. Compared to the NCP field used in M24, the foreground wedge in the NT04 field is significantly better confined, reducing leakage into the EoR window by up to two orders of magnitude (Section 6.1).

Improved A-team subtraction: Since the PSF sidelobe contam-
ination from A-team sources is one of the primary causes of the
excess variance in NenuFAR data, we have implemented two key
improvements in the A-team subtraction step of the calibration
pipeline. First, the data are divided into small time segments and
calibrated against apparent sky models of sources above a flux cut
for each segment. Secondly, the implementation of adaptive solution
intervals in calibration improves the residuals in the direction of
strong A-team sources. We organize all calibration and sky-model
construction steps into a robust pipeline Nenuflow, which performs
these tasks in a heavily parallelized fashion. This is particularly
important as we look into the future, when we plan to integrate
hundreds of hours of this new field, where computational time will
become an important constraint.

Power spectrum upper limits: After sky model subtraction and residual foreground removal, we set upper limits on the 21-cm signal power spectrum at two redshift bins. At $z = 20.3$, the best $2\sigma$ upper limit, derived from 26.1 h of observations, is $\Delta_{21}^2 < 4.6 \times 10^5 \text{ mK}^2$ at $k = 0.038 \text{ h cMpc}^{-1}$, which is the strongest Cosmic Dawn power spectrum limit to date and improve upon previous results by over an order of magnitude. At $z = 17.0$, the best $2\sigma$ limit, from 23.6 h of observations, is $\Delta_{21}^2 < 5.0 \times 10^6 \text{ mK}^2$ at $k = 0.041 \text{ h cMpc}^{-1}$. While this limit is less stringent than the one at $z = 20.3$, it is still deeper than Cosmic Dawn power spectrum limits obtained from all other instruments. A comparison of the upper limits at $z = 20.3$ against simulated exotic 21-cm signals indicates that only the most extreme models, which predict a global signal significantly stronger than the absorption trough observed by EDGES, are currently excluded by these upper limits. An order-of-magnitude improvement in these limits would allow us to probe models with signal strengths comparable to the EDGES detection.

Coherence of excess variance: The new redshift bin we analysed at *z* = 17.0 is severely affected by RFI and exhibits significantly stronger excess variance than *z* = 20.3. Investigation of the nature of the excess variance at both redshift bins suggests that while the excess in the *z* = 20.3 bin is largely incoherent, the excess in the *z* = 17.0 bin has a coherent component that begins to dominate the power in the EoR window after several hours of integration. The redshift bin centred at *z* = 20.3 thus offers a better window into 21-cm cosmology studies for the NenuFAR Cosmic Dawn KSP, and will be the primary focus of future efforts of the KSP.

This paper demonstrates how improved observation and analysis
strategies can help address the challenges of 21-cm power spectrum
detection from the Cosmic Dawn. The derived power spectrum limits,
for the first time, allow us to reach astrophysically relevant 21-cm
signal levels during the Cosmic Dawn, demonstrating the capability
of NenuFAR to provide strong constraints on the 21-cm signal from
the Cosmic Dawn as more observations are integrated in the future.

ACKNOWLEDGEMENTS

SM, JKC, LVEK, SAB, and SG acknowledge the financial support
from the European Research Council (ERC) under the European
Union's Horizon 2020 research and innovation programme (Grant
agreement No. 884760, 'CoDEX'). FGM acknowledges support
from the I-DAWN project, funded by the DIM-ORIGINS pro-
gramme. RB and SS acknowledge the support of the Israel Science
Foundation (grant no. 1078/24). AB acknowledges financial support
from the INAF initiative 'IAF Astronomy Fellowships in Italy' (grant
name MEGASKAT). RG acknowledges support from SERB, DST
Ramanujan Fellowship no. RJF/2022/000141. AKS acknowledges
support from National Science Foundation (grant no. 2206602). This
paper is based on data obtained using the NenuFAR radiotelescope.
NenuFAR has benefitted from the following funding sources: CNRS-
INSU, Observatoire de Paris, Station de Radioastronomie de Nançay,
Observatoire des Sciences de l'Univers de la Région Centre, Région
Centre-Val de Loire, Université d'Orléans, DIM-ACAV and DIM-
ACAV + de la Région Ile de France, Agence Nationale de la
Recherche. We acknowledge the Nançay Data Center resources used
for data reduction and storage.

DATA AVAILABILITY

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

Abdurashidova Z. et al., 2022a, ApJ, 924, 51

Abdurashidova Z. et al., 2022b, ApJ, 925, 221

Acharya A. et al., 2024, MNRAS, 527, 7835

Adams T. et al., 2023, ApJ, 945, 124

Baek S., Di Matteo P., Semelin B., Combes F., Revaz Y., 2009, A&A, 495, 389

Baek S., Semelin B., Di Matteo P., Revaz Y., Combes F., 2010, A&A, 523, A4

Barkana R., 2018, Nature, 555, 71

Barry N. et al., 2019, ApJ, 884, 1

Bennett A., Simth F., 1962, MNRAS, 125, 75

Berlin A., Hooper D., Krnjaic G., McDermott S. D., 2018, Phys. Rev. Lett., 121, 011102

Bonaldi A. et al., 2025, preprint (arXiv:2503.11740)

Bowman J. D., Rogers A. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67

Brackenhoff S. et al., 2024, MNRAS, 533, 632

Brackenhoff S. et al., 2025, MNRAS, 541, 3993

Bull P. et al., 2024, preprint (arXiv:2410.00076)

Cang J., Mesinger A., Murray S. G., Breitman D., Qin Y., Trotta R., 2025, A&A, 698, A152

Ceccotti E. et al., 2025a, preprint (arXiv:2504.18534)

Ceccotti E. et al., 2025b, A&A, 696, A56

Charles N., Kern N., Bernardi G., Bester L., Smirnov O., Fagnoni N., Acedo E. d. L., 2023, MNRAS, 522, 1009

Charles N. et al., 2024, MNRAS, 534, 3349

Chege J. et al., 2024, A&A, 692, A211

Chokshi A., Barry N., Line J., Jordan C., Pindor B., Webster R., 2024, MNRAS, 534, 2475

Cohen A., Lane W., Cotton W., Kassim N., Lazio T., Perley R., Condon J., Erickson W., 2007, AJ, 134, 1245

Condon J. J., Cotton W., Greisen E., Yin Q., Perley R. A., Taylor G., Broderick J., 1998, AJ, 115, 1693

Datta K. K., Mellema G., Mao Y., Iliev I. T., Shapiro P. R., Ahn K., 2012, MNRAS, 424, 1877

Datta K. K., Jensen H., Majumdar S., Mellema G., Iliev I. T., Mao Y., Shapiro P. R., Ahn K., 2014, MNRAS, 442, 1491

Dhandha J. et al., 2025, preprint (arXiv:2503.21687)

Di Tommaso P., Chatzou M., Floden E. W., Barja P. P., Palumbo E., Notredame C., 2017, Nat. Biotechnol., 35, 316

Dowell J., Taylor G. B., 2018a, ApJ, 858, L9

Dowell J., Taylor G. B., 2018b, ApJ, 858, L9

Eastwood M. W. et al., 2019, AJ, 158, 84

Ewall-Wice A. et al., 2016, MNRAS, 460, 4320

Ewall-Wice A., Chang T.-C., Lazio J., Doré O., Seiffert M., Monsalve R., 2018, ApJ, 868, 63

Feng C., Holder G., 2018a, ApJ, 858, L17

Feng C., Holder G., 2018b, ApJ, 858, L17

Fialkov A., Barkana R., 2019, MNRAS, 486, 1763

Fialkov A., Barkana R., Tseliakhovich D., Hirata C. M., 2012, MNRAS, 424, 1335

Fialkov A., Barkana R., Visbal E., Tseliakhovich D., Hirata C. M., 2013, MNRAS, 432, 2909

Fialkov A., Barkana R., Visbal E., 2014, Nature, 506, 197

Fialkov A., Barkana R., Cohen A., 2018, Phys. Rev. Lett., 121, 011101

Fixsen D. J. et al., 2011, ApJ, 734, 5

Fragos T. et al., 2013a, ApJ, 764, 41

Fragos T., Lehmer B. D., Naoz S., Zezas A., Basu-Zych A., 2013b, ApJ, 776, L31

Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181

Gan H. et al., 2022, A&A, 663, A9

Gan H. et al., 2023, A&A, 669, A20

Garsden H. et al., 2021, MNRAS, 506, 5802

Garsden H. et al., 2024, MNRAS, 535, 3218

Gehlot B. et al., 2019, MNRAS, 488, 4271

Gehlot B. et al., 2020, MNRAS, 499, 4158

Gehlot B. et al., 2024, A&A, 681, A71

Ghara R. et al., 2020, MNRAS, 493, 4728

Ghara R. et al., 2025, A&A, 699, A109

Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021, MNRAS, 500, 5322

Gürkan G. et al., 2018, MNRAS, 475, 3010

Hardcastle M. J. et al., 2016, MNRAS, 462, 1910

Höfer C. et al., 2025, preprint (arXiv:2504.03554)

Hutter A., 2018, MNRAS, 477, 1549

Intema H. T., Jagannathan P., Mooley K. P., Frail D. A., 2017, A&A, 598, A78

Kern N. S., Parsons A. R., Dillon J. S., Lanman A. E., Fagnoni N., de Lera Acedo E., 2019, ApJ, 884, 105

Kern N. S. et al., 2020, ApJ, 888, 70

Kolopanis M. et al., 2019, ApJ, 883, 133

Koopmans L. et al., 2015, preprint (arXiv:1505.07568)

de Lera Acedo E. et al., 2022, Nat. Astron., 6, 984

Li W. et al., 2019, ApJ, 887, 141

Liu H., Outmezguine N. J., Redigolo D., Volansky T., 2019, Phys. Rev. D, 100, 123011

Liu B. et al., 2025, MNRAS, 541, 3113

Lloyd S., 1982, IEEE Trans. Inf. Theory, 28, 129

Loh A., coutouly, Girard J. N., EmilieMauduit, Viou C., 2023, Alan-Loh/nenupy: 2023 release, Zenodo, https://doi.org/10.5281/zenodo.7994526

MacQueen J., 1967, in Le Cam L. M., Neyman J., eds, Proc. Fifth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1: Statistics, Univ. California Press, p. 281

Magg M. et al., 2022, MNRAS, 514, 4433

Meriot R., Semelin B., Cornu D., 2025, A&A, 698, A80

Mertens F., Ghosh A., Koopmans L., 2018, MNRAS, 478, 3640

Mertens F. G. et al., 2020, MNRAS, 493, 1662

Mertens F. G., Semelin B., Koopmans L. V. E., 2021, in Siebert A. et al., eds, SF2A-2021: Proc. Annual meeting of the French Society of Astronomy and Astrophysics, Société Française d'Astronomie et d'Astrophysique, p. 211

Mertens F. G., Bobin J., Carucci I. P., 2024, MNRAS, 527, 3517

Mertens F. et al., 2025, A&A, 698, A186

Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955

Mesinger A., Greig B., Sobacchi E., 2016, MNRAS, 459, 2342

Mevius M. et al., 2022, MNRAS, 509, 3693

Mirocha J., Furlanetto S. R., 2019, MNRAS, 483, 1980

Monsalve R. et al., 2024, MNRAS, 530, 4125

Morales M. F., Hazelton B., Sullivan I., Beardsley A., 2012, ApJ, 752, 137

Mouri Sardarabadi A., Koopmans L., 2019, MNRAS, 483, 5480

Muñoz J. B., Loeb A., 2018, Nature, 557, 684

Munshi S. et al., 2024, A&A, 681, A62

Munshi S. et al., 2025a, A&A, 693, A276

Munshi S. et al., 2025b, A&A, 697, A203

Murray Steven G., Bradley G., Andrei M., Muñoz Julian B., Yuxiang Q., Jaehong P., Watkinson Catherine A., 2020, Open Source Softw, 5, 2582

Nunhokee C. et al., 2025, preprint (arXiv:2505.09097)

Offringa A., 2016, A&A, 595, A99

Offringa A., Van De Gronde J., Roerdink J., 2012, A&A, 539, A95

Offringa A. et al., 2014, MNRAS, 444, 606

Offringa A., Mertens F., Van der Tol S., Veenboer B., Gehlot B., Koopmans L., Mevius M., 2019, A&A, 631, A12

Paciga G. et al., 2013, MNRAS, 433, 639

Pandey V., Koopmans L., Tiesinga E., Albers W., Koers H., 2020, ASP Conf. Ser. Vol. 527, Astronomical Data Analysis Software and Systems XXIX. Astron. Soc. Pac., San Francisco, p. 473

Patil A. H. et al., 2016, MNRAS, 463, 4317

Patil A. et al., 2017, ApJ, 838, 65

Philip L. et al., 2019, J. Astron. Instrum., 08, 1950004

Planck Collaboration XVI, 2014, A&A, 571, A16

Planck Collaboration XIII, 2016, A&A, 594, A13

Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys., 75, 086901

Reis I., Fialkov A., Barkana R., 2020, MNRAS, 499, 5993
Rengelink R. B., Tang Y., de Bruyn A. G., Miley G. K., Bremer M. N.,
Roettgering H. J. A., Bremer M. A. R., 1997, A&AS, 124, 259

Semelin B., 2016, MNRAS, 455, 962

Semelin B., Combes F., Baek S., 2007, A&A, 474, 365

Semelin B., Eames E., Bolgar F., Caillat M., 2017, MNRAS, 472, 4508

Semelin B. et al., 2023, A&A, 672, A162

Sikder S., Barkana R., Fialkov A., Reis I., 2024a, MNRAS, 527, 10975

Sikder S., Barkana R., Fialkov A., 2024b, ApJ, 970, L25

Singh S. et al., 2022, Nat. Astron., 6, 607

Sokolowski M. et al., 2015, Publ. Astron. Soc. Aust., 32, e004

Speagle J. S., 2020, MNRAS, 493, 3132

Trott C. M. et al., 2020, MNRAS, 493, 4711

van Haarlem M. P. et al., 2013, A&A, 556, A2

Vanfček P., 1969, Ap&SS, 4, 387

Visbal E., Barkana R., Fialkov A., Tseliakhovich D., Hirata C. M., 2012,
Nature, 487, 70

Yoshiura S. et al., 2021, MNRAS, 505, 4775

Yung L. Y. A., Somerville R. S., Iyer K. G., 2025, preprint (arXiv:2504.18618)

Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011,
MNRAS, 414, 727

Zarka P., Girard J., Tagger M., Denis L., 2012, in Boissier S. et al., eds,
SF2A-2012: Proc. Annual meeting of the French Society of Astronomy
and Astrophysics, Société Française d'Astronomie et d'Astrophysique,
p. 687

Zarka P. et al., 2015, in International Conference on Antenna Theory and
Techniques (ICATT), IEEE, p. 1

Zarka P. et al., 2020, in URSI GASS 2020

APPENDIX A: STATISTICS-BASED SEGMENT
FLAGGING

Before power spectrum estimation, we perform a round of flagging
of a selection of 12 min segments based on visibility statistics. We
compute the standard deviation of the residual visibilities after A-
team subtraction and DI correction as a function of time. Fig. A1
shows this for all four nights of the two spectral windows. The
decrease in variance during the middle of the observation is caused
by the increased sensitivity as the target field moves closer to the
zenith. The vertical dashed lines indicate the LST range where the
target field is above 50° elevation, which is ultimately used in upper
limits estimation. We find that Z20 is significantly cleaner than Z17,
which exhibits strong peaks in the residual standard deviation at
different time segments for each night. We first detrend the standard

<figure><img src="image_19.png" /><figcaption>**Figure A1.** Standard deviation of residual visibilities after A-team subtraction and DI correction. The vertical dashed lines indicate the LSTs where the phase centre crosses 50° elevation. The grey bands indicate the time segments that are flagged before power spectrum estimation.</figcaption></figure>

deviation curve using a running mean and identify and flag the
time segments containing 3σ outliers. The time segments that were
flagged in this step are indicated with grey vertical bands. Flagging
these time segments completely before power spectrum estimation
significantly reduces the excess variance for Z17.

APPENDIX B: POSTERIOR DISTRIBUTION OF
GPR HYPERPARAMETERS

The posterior distribution of the hyperparameters in the Gaussian
process covariance model used in both redshift bins are presented in
Fig. B1.
<figure><img src="image_20.png" /></figure>

$z = 20.3$

<figure><img src="image_21.png" /></figure>

$z = 17.0$

**Figure B1.** Corner plots showing the posterior distribution of the Gaussian process model hyperparameters for the two redshift bins. The contours in the 2D histograms show the 68 per cent and 95 per cent confidence intervals, and the vertical dashed lines in the 1D histograms indicate the 68 per cent confidence intervals.

This paper has been typeset from a TeX/LATEX file prepared by the author.