signal loss for which a bias correction was performed when calculating the power spectra (Sect. 8.4).

For 3C subtraction and NCP subtraction, we find that using a baseline cut of 40λ does result in higher residual power after subtraction compared to when we used a 20λ cut. However, we get negligible signal suppression for a 40λ cut and the increase in residual power is less than the inverse of the signal suppression factor for a 20λ cut. Thus, we used a 40λ cut for both the 3C subtraction and NCP subtraction steps. A frequency smoothing kernel of 2 MHz was used for the 3C subtraction step since it leads to a more accurate subtraction of the 3C sources, which are far from the phase center. For the NCP subtraction step, we used a 6 MHz smoothing kernel since it leads to better source subtraction than a 2 MHz kernel. This is because the phase center is at the NCP and a 6 MHz kernel is sufficient to account for the frequency behavior of the main lobe of the primary beam. Using a 6 MHz kernel also increases the S/N of the gain solutions and we find that the noise in the gain amplitude and phase solutions is smaller when a 6 MHz kernel is used for NCP subtraction. (Appendix B.2). The S/N becomes an important issue for the NCP subtraction since the total flux is divided into multiple clusters and it is necessary to make sure that the faintest cluster has enough S/N to yield reliable calibration solutions.

## 8.2. Robustness test on the full calibration pipeline

Once the results of the robustness tests for each DD subtraction step were assessed and the final settings to use in each processing and analysis step were selected, we performed a signal injection test on the entire calibration pipeline. This involved adding the scaled 21 cm simulated visibilities to the L2 data and passing it through all the steps of A-team subtraction and DI correction, 3C subtraction, and NCP subtraction one by one. The residual data cube without the injected signal was subtracted from the residual data cube with the injected signal to obtain the residual 21 cm signal data cube. The ratio of the power spectrum of this data cube with that of the injected signal then yields the suppression factor for the entire calibration pipeline. We note that the signal was injected before DI correction, so the calculated DI gains get applied to the injected signal in the course of passing the data through the pipeline. So for all steps after DI correction, to calculate the suppression factor, it is necessary to divide the power spectrum of the residual signal by the power spectrum of the injected signal with these DI gains applied.

The fourth column in Fig. 13 shows the results of injecting the signal to the L2 data and passing it through the pipeline with the final settings. We find that in addition to the A-team subtraction step, the DI correction step introduces a small additional suppression (2% of the signal). In the top panel, the additional suppression after the A-team subtraction step till the NCP subtraction step comes almost solely from the DI correction step. This additional factor is due to the fact that the calculated DI gains in the direction of the NCP are slightly higher when there is an injected signal. We repeated this signal injection test with twice the magnitude of the injected signal (Signal 2 in the top row, fourth and fifth columns of Fig. 9). We find that the signal suppression percentage in all DD calibration steps and in the DI correction step remains the same. This is in agreement with the results of Mouri Sardarabadi & Koopmans (2019) who conclude that in the linear regime where the signal is much weaker than the foregrounds, the suppression factor due to calibration is independent of the injected signal strength. In the bottom row of Fig. 13, we can see the factor by which the signal will be suppressed in the $k_⊥, k_‖$ space due to the calibration pipeline. It should be

noted that we only used the baseline range 20λ–40λ (the region between the vertical dotted lines) to construct the final spherical power spectra.

## 8.3. Robustness test on ML-GPR

Since we applied ML-GPR to the gridded visibility cubes, an additional signal injection test was done on the gridded data cube after NCP subtraction. This allows us to test the robustness of the ML-GPR step using a large variety of simulated 21 cm signals. This was not necessary for the calibration steps since the exact shape of the signal does not affect the signal injection results during calibration, as long as the injected signal is small compared to the foregrounds (Mouri Sardarabadi & Koopmans 2019). The decoder of the trained VAE kernel (described in Sect. 7) was used to generate power spectra corresponding to uniformly spaced points in the two-dimensional latent space. We chose five points between the values –2 and 2 for both $x_1$ and $x_2$, making a total of 25 different shapes of the 21 cm power spectrum. For each such signal, the power was scaled to be equal to 0.5, 1, 4, and 20 times the thermal noise power, thus covering a large range of intensities of the 21 cm signal for which we want to test the performance of ML-GPR. Each of these 100 mock 21 cm signals was separately injected into the data. To do this, one realization of the visibility cube was predicted per signal and added to the complex gridded data cube before ML-GPR. The data cube with the injected signal was next used to perform GPR with the same priors as was used for the actual data. The residual power spectrum of the data without an injected signal was then subtracted from the residual power spectrum obtained by applying ML-GPR to data with an injected signal. The workflow used in the signal injection tests for ML-GPR is shown on the right of Fig. 12. The recovered 21 cm power spectrum from the injection test ($\Delta_{\text{rec}}^2$) can be compared to the power spectrum of the injected signal ($\Delta_{\text{inj}}^2$) to see if the signal is absorbed in the procedure. We used two metrics to quantify the results of the signal injection simulations, the z-score and bias, defined as follows:

$$ 
\begin{aligned}
\text{z-score}(k) &= \frac{\Delta_{\text{rec}}^2(k) - \Delta_{\text{inj}}^2(k)}{\sigma_{\Delta_{\text{rec}}^2}(k)}, \\
\text{Bias}(k) &= \frac{\Delta_{\text{rec}}^2(k)}{\Delta_{\text{inj}}^2(k)}.
\end{aligned}
\qquad (7) 
$$

The z-score specifies the number of standard deviations by which the recovered signal power spectrum over or underestimates the input. A negative z-score indicates suppression of the signal and a value below –2 would thus indicate a suppression of the signal beyond the 2-sigma upper limits at that particular k bin. The z-scores for all signal injection tests are shown in the top panel of Fig. 14. The different k bins are indicated with different colored points. We find that none of the signals are suppressed beyond the 1-sigma error bars, and thus the 2-sigma upper limits are comfortably valid. The z-score for brighter signals has a large positive value since $\sigma_{\Delta_{\text{rec}}^2}(k)$ have very low values compared to the signal strength for stronger signals. This is because the contribution to the error bar from the uncertainty in the converged hyperparameter values is lower for stronger signals since the signals are estimated with higher precision. The bias tells us the factor by which the recovered signal power is higher than the injected signal. The bias at each k bin for all signal injection tests is illustrated in the bottom panel of Fig. 14. The signals for which the bias values are less than unity are all within the
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Results of the robustness tests performed on the ML-GPR step. Top: z-score for all injected signals at the different *k* bins. The vertical dashed lines demarcate the different strengths of the injected signals. The horizontal gray bands indicate the 1-sigma and 2-sigma levels. Bottom: Bias for all injected signals at the different *k* bins. The 25 cells in each panel correspond to the 25 different shapes of the injected signal. The different rows correspond to the different strengths of the injected signal and the different columns show the results in different *k* bins. The color scale indicates the value of the bias in each case.</figcaption></figure>

1 sigma limits and correspond to the points in the top panel that lie below zero. There is a positive bias, particularly for signals with high $x_2$ values. This suggests that there is some degeneracy between the foregrounds, excess, and the 21 cm components when modeled as Gaussian processes, but it only results in the absorption of the foregrounds into the excess or the 21 cm signal component, which are not subtracted from the data. This is acceptable as long as we set upper limits on the magnitude of the signal, and do not claim a detection.

## 8.4. Residual power spectra

The residual foreground subtracted data cube from ML-GPR was used to compute the spherical power spectra in logarithmically spaced bins in *k* space between $k = 0.035 \text{ h Mpc}^{-1}$ and $k = 0.5 \text{ h Mpc}^{-1}$. To compute the spherical power spectra, we used only the data in the uv range of $20\lambda-40\lambda$, in order to avoid the signal suppression in calibration due to the 3C subtraction and NCP subtraction steps. The spherical and cylindrical power spectra before and after the entire processing are shown in Fig. 15. We find that the residual power spectra are around two orders of magnitude higher than the thermal noise level at the lowest *k* bins, contributed mainly by the strong excess component that is left in the data after GPR. We note that the thermal noise power spectrum for “Residual” is slightly lower than that

<table><caption>**Table 6.** $2\sigma$ upper limits on the 21 cm power spectrum $\Delta_{\text{ul}}^2(k)$ at different *k* bins. $\Delta_{21}^2(k)$ is the residual power, $\Delta_{21,\text{err}}^2(k)$ is the error on the residual power, and $\Delta_N^2(k)$ is the thermal noise power.</caption><thead><tr><th>k bin<br>(h cMpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21, err</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>ul</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>N</sub>(k)<br>(K<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.036–0.055</td><td>17.81</td><td>2.90</td><td>23.60</td><td>0.12</td></tr><tr><td>0.055–0.086</td><td>71.30</td><td>11.11</td><td>93.52</td><td>0.55</td></tr><tr><td>0.086–0.134</td><td>261.63</td><td>30.96</td><td>323.54</td><td>2.01</td></tr><tr><td>0.134–0.207</td><td>456.33</td><td>49.54</td><td>555.41</td><td>8.10</td></tr><tr><td>0.207–0.322</td><td>731.75</td><td>57.18</td><td>846.12</td><td>31.29</td></tr><tr><td>0.322–0.5</td><td>1157.37</td><td>85.06</td><td>1327.49</td><td>104.66</td></tr></tbody></table>

for “Data” since the post-calibration RFI flagging removes 7.5% of the data.

The noise bias at each *k* bin was calculated by passing the time-differenced Stokes-V image cube (see Sect. 6.1) through the power spectrum pipeline, and accounting for the factor of 2. This was subtracted from the residual power spectrum to obtain the noise bias subtracted dimensionless power spectrum. In Sect. 8.2, we found that there is a maximum of 11% suppression in the signal as a cumulative effect of the entire calibration pipeline, coming mainly from the A-team subtraction and DI correction steps in which the $20\lambda$ to $40\lambda$ baselines were used (Fig. 13). To account for this, we performed a bias correction, which amounts to applying the inverse of the calculated average cumulative signal suppression factors at each *k* bin to the residual power spectra. We place the upper limits at a level of $2\sigma$ above these values at the different *k* modes. The power spectrum values with the $2\sigma$ error bars after noise bias subtraction and after bias correction are shown in Fig. 15 in the rightmost panel. The $2\sigma$ upper limits are listed in Table 6. The best $2\sigma$ upper limit is $2.4 \times 10^7 \text{ mK}^2$ at $k = 0.041 \text{ h Mpc}^{-1}$.

# 9. Discussion and next steps

In this section, we discuss the results and the limitations of the current analysis and present the plans for future improvements to the analysis pipeline.

## 9.1. Nature of the “excess”

The residual data after foreground subtraction is seen to be dominated by excess power up to two orders of magnitude beyond the thermal noise limit. This excess power is not possible to subtract from the data using GPR since its spectral behavior is similar to what we expect from the 21 cm signal.

In order to understand the behavior of the excess as a function of time, we utilized the fact that the visibility data were divided into 13 time segments of 52 min duration each (with a 56 min last segment). The gridded data cube for each of these segments was obtained following Sects. 6.1 and 6.2. A separate GPR was performed on each cube, with the length scales fixed to the values obtained for the full data, and the variances left as free parameters, which were optimized in GPR. Fixing the length scales prevents overfitting of the data, and optimizing the variances allows the time dependence of the strength of the GPR components to be modeled over timescales longer than a time segment. The maximum a posteriori solutions of the foreground cubes were subtracted from the data cube of each time segment to give the residual data cubes, which were used to construct the
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Results of the robustness tests performed on the ML-GPR step. Top: z-score for all injected signals at the different *k* bins. The vertical dashed lines demarcate the different strengths of the injected signals. The horizontal gray bands indicate the 1-sigma and 2-sigma levels. Bottom: Bias for all injected signals at the different *k* bins. The 25 cells in each panel correspond to the 25 different shapes of the injected signal. The different rows correspond to the different strengths of the injected signal and the different columns show the results in different *k* bins. The color scale indicates the value of the bias in each case.</figcaption></figure>

1 sigma limits and correspond to the points in the top panel that lie below zero. There is a positive bias, particularly for signals with high $x_2$ values. This suggests that there is some degeneracy between the foregrounds, excess, and the 21 cm components when modeled as Gaussian processes, but it only results in the absorption of the foregrounds into the excess or the 21 cm signal component, which are not subtracted from the data. This is acceptable as long as we set upper limits on the magnitude of the signal, and do not claim a detection.

## 8.4. Residual power spectra

The residual foreground subtracted data cube from ML-GPR was used to compute the spherical power spectra in logarithmically spaced bins in *k* space between $k = 0.035 \text{ hMpc}^{-1}$ and $k = 0.5 \text{ hMpc}^{-1}$. To compute the spherical power spectra, we used only the data in the uv range of $20\lambda-40\lambda$, in order to avoid the signal suppression in calibration due to the 3C subtraction and NCP subtraction steps. The spherical and cylindrical power spectra before and after the entire processing are shown in Fig. 15. We find that the residual power spectra are around two orders of magnitude higher than the thermal noise level at the lowest *k* bins, contributed mainly by the strong excess component that is left in the data after GPR. We note that the thermal noise power spectrum for “Residual” is slightly lower than that

<table><caption>**Table 6.** $2\sigma$ upper limits on the 21 cm power spectrum $\Delta_{ul}^2(k)$ at different *k* bins. $\Delta_{21}^2(k)$ is the residual power, $\Delta_{21,\text{err}}^2(k)$ is the error on the residual power, and $\Delta_N^2(k)$ is the thermal noise power.</caption><thead><tr><th>k bin<br>(h cMpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21, err</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>ul</sub>(k)<br>(K<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>N</sub>(k)<br>(K<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.036–0.055</td><td>17.81</td><td>2.90</td><td>23.60</td><td>0.12</td></tr><tr><td>0.055–0.086</td><td>71.30</td><td>11.11</td><td>93.52</td><td>0.55</td></tr><tr><td>0.086–0.134</td><td>261.63</td><td>30.96</td><td>323.54</td><td>2.01</td></tr><tr><td>0.134–0.207</td><td>456.33</td><td>49.54</td><td>555.41</td><td>8.10</td></tr><tr><td>0.207–0.322</td><td>731.75</td><td>57.18</td><td>846.12</td><td>31.29</td></tr><tr><td>0.322–0.5</td><td>1157.37</td><td>85.06</td><td>1327.49</td><td>104.66</td></tr></tbody></table>

for “Data” since the post-calibration RFI flagging removes 7.5% of the data.

The noise bias at each *k* bin was calculated by passing the time-differenced Stokes-V image cube (see Sect. 6.1) through the power spectrum pipeline, and accounting for the factor of 2. This was subtracted from the residual power spectrum to obtain the noise bias subtracted dimensionless power spectrum. In Sect. 8.2, we found that there is a maximum of 11% suppression in the signal as a cumulative effect of the entire calibration pipeline, coming mainly from the A-team subtraction and DI correction steps in which the $20\lambda$ to $40\lambda$ baselines were used (Fig. 13). To account for this, we performed a bias correction, which amounts to applying the inverse of the calculated average cumulative signal suppression factors at each *k* bin to the residual power spectra. We place the upper limits at a level of $2\sigma$ above these values at the different *k* modes. The power spectrum values with the $2\sigma$ error bars after noise bias subtraction and after bias correction are shown in Fig. 15 in the rightmost panel. The $2\sigma$ upper limits are listed in Table 6. The best $2\sigma$ upper limit is $2.4 \times 10^7 \text{ mK}^2$ at $k = 0.041 \text{ hMpc}^{-1}$.

# 9. Discussion and next steps

In this section, we discuss the results and the limitations of the current analysis and present the plans for future improvements to the analysis pipeline.

## 9.1. Nature of the “excess”

The residual data after foreground subtraction is seen to be dominated by excess power up to two orders of magnitude beyond the thermal noise limit. This excess power is not possible to subtract from the data using GPR since its spectral behavior is similar to what we expect from the 21 cm signal.

In order to understand the behavior of the excess as a function of time, we utilized the fact that the visibility data were divided into 13 time segments of 52 min duration each (with a 56 min last segment). The gridded data cube for each of these segments was obtained following Sects. 6.1 and 6.2. A separate GPR was performed on each cube, with the length scales fixed to the values obtained for the full data, and the variances left as free parameters, which were optimized in GPR. Fixing the length scales prevents overfitting of the data, and optimizing the variances allows the time dependence of the strength of the GPR components to be modeled over timescales longer than a time segment. The maximum a posteriori solutions of the foreground cubes were subtracted from the data cube of each time segment to give the residual data cubes, which were used to construct the
<figure><img src="image_15.png" /></figure>

**Fig. 15.** Spherical and cylindrical power spectra at some key stages of the analysis pipeline. The left-most panel shows the spherical power spectra after preprocessing (“Data”), after sky model subtraction (“Skymodel Sub”), and after GPR (“GPR Residual”), along with the thermal noise power spectrum (“Thermal Noise”). For the cylindrical power spectra (second and third panels), the ratio with respect to the thermal noise power spectrum is shown. The white dashed lines indicate the horizon limit. The spherical power spectra after noise bias subtraction and suppression factor correction are shown in the rightmost panel.

respective power spectra, shown in the left panel of Fig. 16. The excess is stronger in the first six segments, and the last seven segments have a similar power level. We note that the last segment has the lowest power because it is 56 min long and has a lower thermal noise power level than the other 52 min segments.

Possible contributors to this excess are the residuals from Cyg A, which dominate the *uv* plane after sky model subtraction as seen in Fig. 6. Also, Cyg A goes to lower elevations as time progresses and the beam gain in its direction has very low values after the first 6 segments (Fig. 4). The residuals from Cyg A subtraction are expected to be directly related to its apparent brightness and this is confirmed by wide-field images of the data (for example, the images on the left panels of Fig. 5) where the residuals at its location are seen to go down with time. Thus, Cyg A residuals are expected to be a major contributor to the excess. It is important to note that this excess power is a result of the variance on the NCP field caused by the residual PSF sidelobes of Cyg A. Therefore, even if the peak flux at the position of the Cyg A is reduced to be below the thermal noise limit, it does not necessarily remove the effect of its PSF sidelobes on the NCP field where the power spectrum is computed. This can only be improved in the future with careful time and frequency modeling of the beam, ionosphere, and any other DD effects, during calibration.

To determine if the excess power is coherent with time within a single night, the power spectrum was constructed with data from an increasing number of time segments included. This was done by computing the weighted average (using the *uv* weights) of the residual data cubes for those segments and constructing the power spectra. It should be noted that the weighted average of all the time segments retrieves the exact full data cube since the gridding is linear. The residual power is seen to go down as we integrate more data (right panel of Fig. 16). We find that the power spectrum indicated by “*N* = 13” has slightly lower power than that of “full data”. This is because, in the case of *N* = 13, GPR has been applied to the individual data cubes with the variances of the individual GPR components allowed to vary in time instead of having a single variance representing the strength of the GPR components for the duration of the entire

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

**Fig. 16.** Spherical power spectra of the data after applying GPR on individual time segments. The left panel shows the power spectra of the GPR residuals for individual time segments of the data and the right panel shows the power spectra with increasing volume of data, in the order of observation.

observation as is the case for full data, thus allowing slightly better foreground modeling and subtraction. The ratio between the power of the first segment and that of the full data is seen to be more than a factor of 20, as opposed to a factor of 13 that we expect from using 13 times as much data if the excess is incoherent between time segments. This is likely due to the first few segments having more power from the Cyg A residuals, as discussed above. The fact that the power goes down by a factor of 13 or more as we integrate 13 times as much data suggests that the excess power, which dominates the residuals,
<figure><img src="image_16.png" /></figure>

**Fig. 15.** Spherical and cylindrical power spectra at some key stages of the analysis pipeline. The left-most panel shows the spherical power spectra after preprocessing (“Data”), after sky model subtraction (“Skymodel Sub”), and after GPR (“GPR Residual”), along with the thermal noise power spectrum (“Thermal Noise”). For the cylindrical power spectra (second and third panels), the ratio with respect to the thermal noise power spectrum is shown. The white dashed lines indicate the horizon limit. The spherical power spectra after noise bias subtraction and suppression factor correction are shown in the rightmost panel.

respective power spectra, shown in the left panel of Fig. 16. The excess is stronger in the first six segments, and the last seven segments have a similar power level. We note that the last segment has the lowest power because it is 56 min long and has a lower thermal noise power level than the other 52 min segments.

Possible contributors to this excess are the residuals from Cyg A, which dominate the *uv* plane after sky model subtraction as seen in Fig. 6. Also, Cyg A goes to lower elevations as time progresses and the beam gain in its direction has very low values after the first 6 segments (Fig. 4). The residuals from Cyg A subtraction are expected to be directly related to its apparent brightness and this is confirmed by wide-field images of the data (for example, the images on the left panels of Fig. 5) where the residuals at its location are seen to go down with time. Thus, Cyg A residuals are expected to be a major contributor to the excess. It is important to note that this excess power is a result of the variance on the NCP field caused by the residual PSF sidelobes of Cyg A. Therefore, even if the peak flux at the position of the Cyg A is reduced to be below the thermal noise limit, it does not necessarily remove the effect of its PSF sidelobes on the NCP field where the power spectrum is computed. This can only be improved in the future with careful time and frequency modeling of the beam, ionosphere, and any other DD effects, during calibration.

To determine if the excess power is coherent with time within a single night, the power spectrum was constructed with data from an increasing number of time segments included. This was done by computing the weighted average (using the *uv* weights) of the residual data cubes for those segments and constructing the power spectra. It should be noted that the weighted average of all the time segments retrieves the exact full data cube since the gridding is linear. The residual power is seen to go down as we integrate more data (right panel of Fig. 16). We find that the power spectrum indicated by “*N* = 13” has slightly lower power than that of “full data”. This is because, in the case of *N* = 13, GPR has been applied to the individual data cubes with the variances of the individual GPR components allowed to vary in time instead of having a single variance representing the strength of the GPR components for the duration of the entire

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

**Fig. 16.** Spherical power spectra of the data after applying GPR on individual time segments. The left panel shows the power spectra of the GPR residuals for individual time segments of the data and the right panel shows the power spectra with increasing volume of data, in the order of observation.

observation as is the case for full data, thus allowing slightly better foreground modeling and subtraction. The ratio between the power of the first segment and that of the full data is seen to be more than a factor of 20, as opposed to a factor of 13 that we expect from using 13 times as much data if the excess is incoherent between time segments. This is likely due to the first few segments having more power from the Cyg A residuals, as discussed above. The fact that the power goes down by a factor of 13 or more as we integrate 13 times as much data suggests that the excess power, which dominates the residuals,
is incoherent. Thus, as we integrate more nights of observations, it is possible that this excess could go down significantly with time if a part of the excess is incoherent across different nights. This needs further investigation by processing multiple nights of data. However, even if this is confirmed, the final power spec- tra would remain well above what can be expected from thermal noise alone, which also averages down in a similar manner.

## 9.2. Impact of wide-field sources

The NenuFAR primary beam has very strong grating lobes roughly 60 degrees away from the phase center. Wide-field images of the data reveal a large number of strong radio sources falling in these grating lobes, which remain in the data after sky model subtraction. In the 3C subtraction step, seven of the brightest among these sources were subtracted through DD calibration. However, the remaining sources have too little S/N to obtain reliable calibration solutions. Though these sources are far away from the phase center, the PSF sidelobes of these sources will still pass through the NCP field and add a significant contribution to the excess power. To understand the impact of these sources far away from the phase center, we computed wide-field MFS dirty images from the data for each of the time segments into which the data had been divided. Each image was divided into four annular regions around the NCP, labeled “Main Lobe”, “Nulls”, “Grating Lobes”, and “Edge”, based on the amplitude of the radial profile of the NenuFAR primary beam (bottom panel in Fig. 17). For the pixels within the nth annulus, the cross-coherence between the ith and jth time segments is given by

$$C_{ij}(n) = \frac{\langle I_i(l, m)I_j(l, m) \rangle}{\sqrt{\langle I_i^2(l, m) \rangle\langle I_j^2(l, m) \rangle}}, \quad (8)$$

where $\langle \cdots \rangle$ indicates an average over all $l, m$ lying within the nth annulus. Figure 17 shows the cross-coherence at four different calibration and foreground subtraction stages. It should be noted that the cross-coherence of a segment with itself is unity by definition.

Before A-team subtraction, at the “data” level, the nulls annulus is dominated by flux from Cas A, and the grating lobes annulus is dominated by Cyg A. Cas A is high in the sky throughout the observation, but Cyg A goes below the horizon halfway through the observation, resulting in lower cross-coherence values for the later time segments. After A-team subtraction, the grating lobes annulus is dominated by the strong 3C sources. The 3C subtraction step removes seven of the brightest sources, but the remaining sources have a significant cross-coherence in the grating lobes annulus. There is also higher cross-coherence for the separation of five time segments, corresponding to a time difference of 4.3 h. This is possibly because of the six-fold symmetry of the hexagonal MAs, which results in the grating lobes in the sky that overlap exactly once every 4 h. In the main lobe, the cross-coherence values are close to unity until before NCP subtraction. NCP subtraction decreases the cross-coherence in the main lobe and the modulation seen in the cross-coherence with increasing time difference is because there are large residuals from sources at the edge of the main lobe. As the elongated main lobe rotates against the sky, it completes close to a 180 degree rotation in the duration of the observation as a result of which the time segments separated by 13 segments again have high cross-coherence values. We note that this modulation due to the rotation of the non-axisymmetric main lobe of the beam is also seen

<figure><img src="image_17.png" /><figcaption>Fig. 17. Cross-coherence of the data between different time segments. The top panel shows the cross-coherence computed in four different regions of the primary beam (different columns), for the data at four different stages of sky model subtraction (different rows). The bottom panel shows what the four regions correspond to in the radial profile of the time-averaged primary beam.</figcaption></figure>

before NCP subtraction, but the effect is lower since the field has a more uniform distribution of flux before NCP subtraction. The edge region is dominated by noise, so there is negligible cross-coherence between time segments.

These cross-coherence results show that the timescale of coherence is heavily dependent on the distance from the phase center because of the primary beam. For the main lobe, this timescale is nearly the duration of the full night, while for the nulls, the coherence drops to 20% or less for consecutive segments. For the grating lobes, the coherence is about 10% for the separation of between one and six segments and falls to less than 2% for longer timescales. If the edge region is a reasonable estimate of the noise, we can state that the coherence of the nulls and grating lobes (and, of course, the main lobe) are still dominated by real sources and not too much diluted by noise. These sources behave incoherently in their flux or position from segment to segment due to the primary beam fixed to the Earth. Thus, the power from these sources, coupled with the non-smooth spectral nature of the primary beam at the nulls and grating lobes, could have a large contribution to the excess power at high $k_\parallel$ seen in the residual data, which averages down with time like incoherent noise (as seen in Sect. 9.1), on timescales that exceed their own coherence time (as seen in Fig. 17).
is incoherent. Thus, as we integrate more nights of observations, it is possible that this excess could go down significantly with time if a part of the excess is incoherent across different nights. This needs further investigation by processing multiple nights of data. However, even if this is confirmed, the final power spectra would remain well above what can be expected from thermal noise alone, which also averages down in a similar manner.

## 9.2. Impact of wide-field sources

The NenuFAR primary beam has very strong grating lobes roughly 60 degrees away from the phase center. Wide-field images of the data reveal a large number of strong radio sources falling in these grating lobes, which remain in the data after sky model subtraction. In the 3C subtraction step, seven of the brightest among these sources were subtracted through DD calibration. However, the remaining sources have too little S/N to obtain reliable calibration solutions. Though these sources are far away from the phase center, the PSF sidelobes of these sources will still pass through the NCP field and add a significant contribution to the excess power. To understand the impact of these sources far away from the phase center, we computed wide-field MFS dirty images from the data for each of the time segments into which the data had been divided. Each image was divided into four annular regions around the NCP, labeled “Main Lobe”, “Nulls”, “Grating Lobes”, and “Edge”, based on the amplitude of the radial profile of the NenuFAR primary beam (bottom panel in Fig. 17). For the pixels within the nth annulus, the cross-coherence between the ith and jth time segments is given by

$$C_{ij}(n) = \frac{\langle I_i(l, m)I_j(l, m) \rangle}{\sqrt{\langle I_i^2(l, m) \rangle\langle I_j^2(l, m) \rangle}}, \quad (8)$$

where $\langle \cdots \rangle$ indicates an average over all $l, m$ lying within the nth annulus. Figure 17 shows the cross-coherence at four different calibration and foreground subtraction stages. It should be noted that the cross-coherence of a segment with itself is unity by definition.

Before A-team subtraction, at the “data” level, the nulls annulus is dominated by flux from Cas A, and the grating lobes annulus is dominated by Cyg A. Cas A is high in the sky throughout the observation, but Cyg A goes below the horizon halfway through the observation, resulting in lower cross-coherence values for the later time segments. After A-team subtraction, the grating lobes annulus is dominated by the strong 3C sources. The 3C subtraction step removes seven of the brightest sources, but the remaining sources have a significant cross-coherence in the grating lobes annulus. There is also higher cross-coherence for the separation of five time segments, corresponding to a time difference of 4.3 h. This is possibly because of the six-fold symmetry of the hexagonal MAs, which results in the grating lobes in the sky that overlap exactly once every 4 h. In the main lobe, the cross-coherence values are close to unity until before NCP subtraction. NCP subtraction decreases the cross-coherence in the main lobe and the modulation seen in the cross-coherence with increasing time difference is because there are large residuals from sources at the edge of the main lobe. As the elongated main lobe rotates against the sky, it completes close to a 180 degree rotation in the duration of the observation as a result of which the time segments separated by 13 segments again have high cross-coherence values. We note that this modulation due to the rotation of the non-axisymmetric main lobe of the beam is also seen

<figure><img src="image_18.png" /><figcaption>Fig. 17. Cross-coherence of the data between different time segments. The top panel shows the cross-coherence computed in four different regions of the primary beam (different columns), for the data at four different stages of sky model subtraction (different rows). The bottom panel shows what the four regions correspond to in the radial profile of the time-averaged primary beam.</figcaption></figure>

before NCP subtraction, but the effect is lower since the field has a more uniform distribution of flux before NCP subtraction. The edge region is dominated by noise, so there is negligible cross-coherence between time segments.

These cross-coherence results show that the timescale of coherence is heavily dependent on the distance from the phase center because of the primary beam. For the main lobe, this timescale is nearly the duration of the full night, while for the nulls, the coherence drops to 20% or less for consecutive segments. For the grating lobes, the coherence is about 10% for the separation of between one and six segments and falls to less than 2% for longer timescales. If the edge region is a reasonable estimate of the noise, we can state that the coherence of the nulls and grating lobes (and, of course, the main lobe) are still dominated by real sources and not too much diluted by noise. These sources behave incoherently in their flux or position from segment to segment due to the primary beam fixed to the Earth. Thus, the power from these sources, coupled with the non-smooth spectral nature of the primary beam at the nulls and grating lobes, could have a large contribution to the excess power at high $k_\parallel$ seen in the residual data, which averages down with time like incoherent noise (as seen in Sect. 9.1), on timescales that exceed their own coherence time (as seen in Fig. 17).
## 9.3. Limitations and future improvement

This first analysis on a single night of NenuFAR observations reveals the limitations of the current analysis pipeline. Here, we describe these limitations and discuss future improvements in the processing pipeline necessary to overcome them.

**Near-field RFI sources.** Near-field imaging (i.e., imaging of sources on the ground plane in the vicinity of the array) of the data after calibration and sky model subtraction revealed that two prominent sources of local RFI contaminate the data. The approach adopted in this analysis to mitigate its impact was by identifying the baselines where these RFI sources have a strong signature in the delay power spectra and flagging those baselines completely. However, this is an overly simplistic approach, and all un-flagged baselines still have some power due to these local RFI sources. It was seen that flagging the affected baselines significantly reduced excess power at high delay modes. So, local RFI in the unflagged baselines could make up much of the excess power we see in the data. A more targeted approach in mitigating the RFI sources could be by filtering them out from delay or fringe rate power spectra of all baselines, utilizing the fact that the RFI sources have a localized signature in the delay-fringe rate space. We defer that to a future analysis. We note that the impact of the near-field RFI sources is expected to be maximum for NCP observations since the NCP is the only point in the northern sky that is fixed with respect to the Earth, and the local stationary RFI sources add up coherently at the NCP even as the sky rotates with time. This is less the case for other fields, not including the NCP, and the contaminating effects of these local RFI sources on observations of these fields could be significantly less.

**Impact of diffuse (polarized) foregrounds.** NenuFAR, with its dense core consisting of many small baselines, is very sensitive to diffuse Galactic foregrounds. These (polarized) foregrounds are not present in the sky model used in calibration, and the power from the diffuse foregrounds can be absorbed in the station gains, leading to calibration errors. We used a 20λ baseline cut in calibration for the A-team subtraction step and a 40λ cut for the other calibration steps. This avoids the power due to diffuse foregrounds at spatial modes measured on baselines shorter than 20λ and 40λ, respectively. However, using a baseline cut results in excess power on the excluded baselines. Thus, both 3C and NCP subtraction result in an excess noise on the baselines smaller than 40λ, which are used in power spectrum generation. The errors in calibration due to diffuse foregrounds can be improved by including diffuse foreground models. However, the ideal way to model the diffuse foregrounds, perhaps using shapelets decomposition (Yatawatta 2011; Gehlot et al. 2022), must still be tested. A future step could also be to build a broadband spectral model of diffuse foregrounds by combining AARTFAAC High Band Antenna and NenuFAR observations.

**Beam model in calibration.** The beam model of NenuFAR is currently not integrated into the calibration algorithm. As a result, we performed calibration against an apparent sky model obtained by multiplying the intrinsic sky model with the average beam for each time segment of data. The effect of the beam on the apparent flux of sources on timescales more than the calibration solution interval is accounted for by the calibration gain solutions. However, the variation of the beam on time and frequency scales shorter than a solution interval is currently not accounted for. The high residuals due to Cyg A (seen in the middle-left panel of Fig. 5 and the middle and right panels of

Fig. 6), which are expected to be a major contributor to the excess variance (Sect. 9.1), is partly a consequence of this. The solution interval used in the A-team subtraction step is not small enough to model the primary beam accurately at the grating lobes in which Cyg A lies. This could be addressed to some extent when the primary beam model is integrated into the calibration software, and the level of improvement will depend on the accuracy of the beam model.

For the NCP subtraction step, including a beam model in the calibration step will allow us to use longer solution intervals since this will partly mitigate the requirement of the gain solutions to represent the rotating non-axisymmetric beam. This will lower the expected thermal noise for each solution interval and thus make it possible to use smaller clusters with enough flux to be significantly beyond the thermal noise level. Thus, it will be possible to account for DD effects on smaller spatial scales than currently possible. Including a beam model in the calibration software will also mean that the gain solutions do not need to account for the frequency behavior of the primary beam. Therefore, even with a large value for the smoothing scale, it will be possible to fit and subtract the sources well, thus avoiding increasing the risk of overfitting by decreasing the width of the frequency smoothing kernel.

Another advantage of having the beam in the calibration step is using a wider sky model of the NCP field. This is because the spatial gradient of the source fluxes due to attenuation by the beam will now be partly accounted for by the beam model, thus allowing us to use wider clusters with higher S/N. We see in the bottom-right panel of Fig. 5 and the bottom-left panel of Fig. 17 that after NCP subtraction, there are residual sources near the edge of the main lobe of the primary beam. Using an extended sky model, along with the beam model, will allow us to subtract these sources better.

Finally, including a beam model will allow us to subtract many sources in the grating lobes of the NenuFAR primary beam. Currently, these sources cannot be subtracted because to model the beam, we needed to calibrate and absorb the beam in the gain solutions. Performing such a calibration run is impossible since the sources do not have enough S/N to yield reliable calibration solutions within the time and frequency solution interval needed to represent the beam faithfully. Once a beam model is available, the visibilities for these sources with the primary beam attenuation can be predicted and subtracted directly from the data. The accuracy of this subtraction will depend on the accuracy of the simulated beam model. We note that all sources that are not subtracted and the residuals from the subtracted sources will continue to have the effect of the rotating beam and will still contribute to the excess power.

**Confusion noise limit.** The NCP subtraction step is limited by confusion noise due to the maximum angular resolution of NenuFAR, which is given by the baseline formed by two remote MAs. Once more remote MAs are available, a 5 km maximum baseline in the future will improve the angular resolution by a factor of ≈3.4 and allow us to construct a deeper sky model due to a lower confusion limit (by a factor of ≈6.7). This will allow us to use smaller clusters in the NCP subtraction step since smaller clusters will have a high enough S/N to give reliable calibration solutions.

**Polarization leakage.** In this analysis, we have used diagonal mode in calibration with DDECal and Stokes-I and Stokes-V data for all steps. However, we expect some polarization leakage in NenuFAR since it is a wide-field phased array
### 9.3. Limitations and future improvement

This first analysis on a single night of NenuFAR observations reveals the limitations of the current analysis pipeline. Here, we describe these limitations and discuss future improvements in the processing pipeline necessary to overcome them.

**Near-field RFI sources.** Near-field imaging (i.e., imaging of sources on the ground plane in the vicinity of the array) of the data after calibration and sky model subtraction revealed that two prominent sources of local RFI contaminate the data. The approach adopted in this analysis to mitigate its impact was by identifying the baselines where these RFI sources have a strong signature in the delay power spectra and flagging those baselines completely. However, this is an overly simplistic approach, and all un-flagged baselines still have some power due to these local RFI sources. It was seen that flagging the affected baselines significantly reduced excess power at high delay modes. So, local RFI in the unflagged baselines could make up much of the excess power we see in the data. A more targeted approach in mitigating the RFI sources could be by filtering them out from delay or fringe rate power spectra of all baselines, utilizing the fact that the RFI sources have a localized signature in the delay-fringe rate space. We defer that to a future analysis. We note that the impact of the near-field RFI sources is expected to be maximum for NCP observations since the NCP is the only point in the northern sky that is fixed with respect to the Earth, and the local stationary RFI sources add up coherently at the NCP even as the sky rotates with time. This is less the case for other fields, not including the NCP, and the contaminating effects of these local RFI sources on observations of these fields could be significantly less.

**Impact of diffuse (polarized) foregrounds.** NenuFAR, with its dense core consisting of many small baselines, is very sensitive to diffuse Galactic foregrounds. These (polarized) foregrounds are not present in the sky model used in calibration, and the power from the diffuse foregrounds can be absorbed in the station gains, leading to calibration errors. We used a 20λ baseline cut in calibration for the A-team subtraction step and a 40λ cut for the other calibration steps. This avoids the power due to diffuse foregrounds at spatial modes measured on baselines shorter than 20λ and 40λ, respectively. However, using a baseline cut results in excess power on the excluded baselines. Thus, both 3C and NCP subtraction result in an excess noise on the baselines smaller than 40λ, which are used in power spectrum generation. The errors in calibration due to diffuse foregrounds can be improved by including diffuse foreground models. However, the ideal way to model the diffuse foregrounds, perhaps using shapelets decomposition (Yatawatta 2011; Gehlot et al. 2022), must still be tested. A future step could also be to build a broadband spectral model of diffuse foregrounds by combining AARTFAAC High Band Antenna and NenuFAR observations.

**Beam model in calibration.** The beam model of NenuFAR is currently not integrated into the calibration algorithm. As a result, we performed calibration against an apparent sky model obtained by multiplying the intrinsic sky model with the average beam for each time segment of data. The effect of the beam on the apparent flux of sources on timescales more than the calibration solution interval is accounted for by the calibration gain solutions. However, the variation of the beam on time and frequency scales shorter than a solution interval is currently not accounted for. The high residuals due to Cyg A (seen in the middle-left panel of Fig. 5 and the middle and right panels of

Fig. 6), which are expected to be a major contributor to the excess variance (Sect. 9.1), is partly a consequence of this. The solution interval used in the A-team subtraction step is not small enough to model the primary beam accurately at the grating lobes in which Cyg A lies. This could be addressed to some extent when the primary beam model is integrated into the calibration software, and the level of improvement will depend on the accuracy of the beam model.

For the NCP subtraction step, including a beam model in the calibration step will allow us to use longer solution intervals since this will partly mitigate the requirement of the gain solutions to represent the rotating non-axisymmetric beam. This will lower the expected thermal noise for each solution interval and thus make it possible to use smaller clusters with enough flux to be significantly beyond the thermal noise level. Thus, it will be possible to account for DD effects on smaller spatial scales than currently possible. Including a beam model in the calibration software will also mean that the gain solutions do not need to account for the frequency behavior of the primary beam. Therefore, even with a large value for the smoothing scale, it will be possible to fit and subtract the sources well, thus avoiding increasing the risk of overfitting by decreasing the width of the frequency smoothing kernel.

Another advantage of having the beam in the calibration step is using a wider sky model of the NCP field. This is because the spatial gradient of the source fluxes due to attenuation by the beam will now be partly accounted for by the beam model, thus allowing us to use wider clusters with higher S/N. We see in the bottom-right panel of Fig. 5 and the bottom-left panel of Fig. 17 that after NCP subtraction, there are residual sources near the edge of the main lobe of the primary beam. Using an extended sky model, along with the beam model, will allow us to subtract these sources better.

Finally, including a beam model will allow us to subtract many sources in the grating lobes of the NenuFAR primary beam. Currently, these sources cannot be subtracted because to model the beam, we needed to calibrate and absorb the beam in the gain solutions. Performing such a calibration run is impossible since the sources do not have enough S/N to yield reliable calibration solutions within the time and frequency solution interval needed to represent the beam faithfully. Once a beam model is available, the visibilities for these sources with the primary beam attenuation can be predicted and subtracted directly from the data. The accuracy of this subtraction will depend on the accuracy of the simulated beam model. We note that all sources that are not subtracted and the residuals from the subtracted sources will continue to have the effect of the rotating beam and will still contribute to the excess power.

**Confusion noise limit.** The NCP subtraction step is limited by confusion noise due to the maximum angular resolution of NenuFAR, which is given by the baseline formed by two remote MAs. Once more remote MAs are available, a 5 km maximum baseline in the future will improve the angular resolution by a factor of ≈3.4 and allow us to construct a deeper sky model due to a lower confusion limit (by a factor of ≈6.7). This will allow us to use smaller clusters in the NCP subtraction step since smaller clusters will have a high enough S/N to give reliable calibration solutions.

**Polarization leakage.** In this analysis, we have used diagonal mode in calibration with DDECal and Stokes-I and Stokes-V data for all steps. However, we expect some polarization leakage in NenuFAR since it is a wide-field phased array
instrument. Additionally, due to its dense core, NenuFAR is par-
ticularly sensitive to diffuse foregrounds, which are expected to
have a polarized component. The leakage of Faraday rotated
polarized foregrounds into Stokes-I can result in the foregrounds
assuming a non-smooth spectral behavior. This can be tested
with simulations once the full Stokes beam model is available,
and the inclusion of such a beam model in calibration can poten-
tially decrease instrumental polarization.

**Ionospheric effects.** An important effect that can impact the power spectrum is ionospheric scintillation noise (Vedantham & Koopmans 2015, 2016; Gehlot et al. 2018). A turbulent ionosphere introduces phase shifts on the incoming electromagnetic waves, which can cause the sources to move around in the sky with time, thus affecting the accuracy of their subtraction through calibration. In this analysis, we have selected a relatively quiet ionospheric night, and the calibration solution time intervals used at each stage are longer than the expected timescales of ionospheric effects, making it impossible to solve for these effects. A detailed simulation of the ionosphere is necessary to quantify its impact on NenuFAR data and will be done in the future.

In summary, a few specific effects could be the dominant
contributors to the excess power we see in the data. Targeted
improvements in the processing pipeline can account for, miti-
gate, and possibly eliminate these effects.

10. Summary

This paper describes the first end-to-end analysis of the spectral
window of 61-72 MHz from an 11.4 h observation of the NCP
field with NenuFAR. We draw the following main conclusions
from this analysis:

- We set the best 2σ upper limits on the 21 cm power spectrum at $2.4 \times 10^7$ mK$^2$ at $k = 0.041$ $h$ cMpc$^{-1}$ and $z = 20.3$. These upper limits are two orders of magnitude higher than the thermal noise power spectrum.

- Localized RFI sources strongly impact the data and are a major contributor to the excess power. The RFI is mitigated to an extent by baseline flagging, and a more targeted approach to filtering these local RFI sources could significantly mitigate the excess power.

- The excess power seems largely incoherent with time within a single night and averages down like noise as we include more data to construct the power spectra. Processing of multiple nights is necessary to understand whether coherent averaging of multiple nights of data will decrease this excess power.

- Many sources lying in the nulls and grating lobes of the NenuFAR primary beam are left in the data after sky model subtraction. Coupled with the frequency fluctuations of the primary beam in these regions, these sources could strongly contribute to the excess power seen in the data. Including a simulated primary beam model in the calibration step would allow us to subtract these sources to some extent and partly mitigate this effect.

We will continue to investigate the origin of this excess power
and implement improvements in the analysis pipeline to mitigate
it. This will allow us to exploit the full sensitivity of NenuFAR.
In parallel, we plan to process multiple nights of observations
and obtain upper limits closer to the thermal noise and expected
21 cm signal levels, which will allow us to constrain the under-
lying astrophysical models.

**Acknowledgements.** SM, LVEK, and BKG 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 the support of the PSL Fellowship. AB acknowledges support from the European Research Council through the Advanced Grant MIST (FP7/2017-2022, No. 742719). EC would like to acknowledge support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen. AF is supported by the Royal Society University Research Fellowship. The postdoctoral contract of IH was funded by Sorbonne Université in the framework of the Initiative Physique des Infinis (IDEX SUPER). We are grateful to Bradley Greig for his valuable input and guidance in setting up 21cmFAST simulations for cosmic dawn. 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.

References

Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221

Adams, T., Aguirre, J. E., Alexander, P., et al. 2023, ApJ, 945, 124

Asad, K., Koopmans, L., Jelić, V., et al. 2015, MNRAS, 451, 3709

Asad, K., Koopmans, L., Jelić, V., et al. 2016, MNRAS, 462, 4482

Asad, K., Koopmans, L., Jelić, V., et al. 2018, MNRAS, 476, 3051

Barkana, R. 2018, Nature, 555, 71

Barkana, R., & Loeb, A. 2005, ApJ, 624, L65

Barry, N., Hazelton, B., Sullivan, I., Morales, M., & Pober, J. 2016, MNRAS, 461, 3135

Barry, N., Wilensky, M., Trott, C., et al. 2019, ApJ, 884, 1

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

Bevins, H. T., Heimersheim, S., Abril-Cabezas, I., et al. 2024, MNRAS, 527, 813

Bowman, J. D., Rogers, A. E., & Hewitt, J. N. 2008, ApJ, 676, 1

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

Bowman, J. D., Rogers, A. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018b, Nature, 564, E35

Broekema, P. C., Mol, J. J. D., Nijboer, R., et al. 2018, Astron. Comput., 23, 180

Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, MNRAS, 423, 2518

Chapman, E., Abdalla, F. B., Bobin, J., et al. 2013, MNRAS, 429, 165

Charrier, D., & the CODALEMA Collaboration 2007, in 2007 IEEE Antennas and Propagation Society International Symposium (IEEE), 4485

Chen, X., & Miralda-Escudé, J. 2004, ApJ, 602, 1

Cohen, A., Lane, W., Cotton, W., et al. 2007, AJ, 134, 1245

Condon, J. J., Cotton, W., Greisen, E., et al. 1998, AJ, 115, 1693

DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001

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

Eastwood, M. W., Anderson, M. M., Monroe, R. M., et al. 2019, AJ, 158, 84

Ewall-Wice, A., Dillon, J. S., Hewitt, J. N., et al. 2016, MNRAS, 460, 4320

Ewall-Wice, A., Chang, T.-C., Lazio, J., et al. 2018, ApJ, 868, 63

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

Fialkov, A., & Barkana, R. 2014, MNRAS, 445, 213

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

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

Field, G. B. 1958, Proc. IRE, 46, 240

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

Gan, H., Mertens, F., Koopmans, L., et al. 2023, A&A, 669, A20

Garsden, H., Greenhill, L., Bernardi, G., et al. 2021, MNRAS, 506, 5802

Gehlot, B. K., Koopmans, L. V., de Bruyn, A., et al. 2018, MNRAS, 478, 1484

Gehlot, B., Mertens, F., Koopmans, L., et al. 2019, MNRAS, 488, 4271

Gehlot, B., Mertens, F., Koopmans, L., et al. 2020, MNRAS, 499, 4158

Gehlot, B., Koopmans, L., Offringa, A., et al. 2022, A&A, 662, A97

Gessey-Jones, T., Sartorio, N. S., Fialkov, A., et al. 2022, MNRAS, 516, 841

Hibbard, J. J., Mirocha, J., Rapetti, D., et al. 2022, ApJ, 929, 151

Hicks, B. C., Paravastu-Dalal, N., Stewart, K. P., et al. 2012, PASP, 124, 1090

Hills, R., Kulkarni, G., Meerburg, P. D., & Puchwein, E. 2018, Nature, 564, E32

Hothi, I., Chapman, E., Pritchard, J. R., et al. 2021, MNRAS, 500, 2264

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

Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2008, MNRAS, 389, 1319

Joseph, R. 2012, Bachelor's Thesis, Faculty of Science and Engineering

Kazemi, S., Yatawatta, S., & Zaroubi, S. 2013, MNRAS, 430, 1457

Kern, N. S., & Liu, A. 2021, MNRAS, 501, 1463
instrument. Additionally, due to its dense core, NenuFAR is par-
ticularly sensitive to diffuse foregrounds, which are expected to
have a polarized component. The leakage of Faraday rotated
polarized foregrounds into Stokes-I can result in the foregrounds
assuming a non-smooth spectral behavior. This can be tested
with simulations once the full Stokes beam model is available,
and the inclusion of such a beam model in calibration can poten-
tially decrease instrumental polarization.

**Ionospheric effects.** An important effect that can impact the power spectrum is ionospheric scintillation noise (Vedantham & Koopmans 2015, 2016; Gehlot et al. 2018). A turbulent ionosphere introduces phase shifts on the incoming electromagnetic waves, which can cause the sources to move around in the sky with time, thus affecting the accuracy of their subtraction through calibration. In this analysis, we have selected a relatively quiet ionospheric night, and the calibration solution time intervals used at each stage are longer than the expected timescales of ionospheric effects, making it impossible to solve for these effects. A detailed simulation of the ionosphere is necessary to quantify its impact on NenuFAR data and will be done in the future.

In summary, a few specific effects could be the dominant
contributors to the excess power we see in the data. Targeted
improvements in the processing pipeline can account for, miti-
gate, and possibly eliminate these effects.

10. Summary

This paper describes the first end-to-end analysis of the spectral
window of 61-72 MHz from an 11.4 h observation of the NCP
field with NenuFAR. We draw the following main conclusions
from this analysis:

- We set the best 2σ upper limits on the 21 cm power spectrum at $2.4 \times 10^7$ mK$^2$ at $k = 0.041$ $h$ cMpc$^{-1}$ and $z = 20.3$. These upper limits are two orders of magnitude higher than the thermal noise power spectrum.

- Localized RFI sources strongly impact the data and are a major contributor to the excess power. The RFI is mitigated to an extent by baseline flagging, and a more targeted approach to filtering these local RFI sources could significantly mitigate the excess power.

- The excess power seems largely incoherent with time within a single night and averages down like noise as we include more data to construct the power spectra. Processing of multiple nights is necessary to understand whether coherent averaging of multiple nights of data will decrease this excess power.

- Many sources lying in the nulls and grating lobes of the NenuFAR primary beam are left in the data after sky model subtraction. Coupled with the frequency fluctuations of the primary beam in these regions, these sources could strongly contribute to the excess power seen in the data. Including a simulated primary beam model in the calibration step would allow us to subtract these sources to some extent and partly mitigate this effect.

We will continue to investigate the origin of this excess power
and implement improvements in the analysis pipeline to mitigate
it. This will allow us to exploit the full sensitivity of NenuFAR.
In parallel, we plan to process multiple nights of observations
and obtain upper limits closer to the thermal noise and expected
21 cm signal levels, which will allow us to constrain the under-
lying astrophysical models.

Acknowledgements. SM, LVEK, and BKG 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 the support of the PSL Fellowship. AB acknowledges support from the European Research Council through the Advanced Grant MIST (FP7/2017-2022, No. 742719). EC would like to acknowledge support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen. AF is supported by the Royal Society University Research Fellowship. The postdoctoral contract of IH was funded by Sorbonne Université in the framework of the Initiative Physique des Infinis (IDEX SUPER). We are grateful to Bradley Greig for his valuable input and guidance in setting up 21cmFAST simulations for cosmic dawn. 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.

References

Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221

Adams, T., Aguirre, J. E., Alexander, P., et al. 2023, ApJ, 945, 124

Asad, K., Koopmans, L., Jelić, V., et al. 2015, MNRAS, 451, 3709

Asad, K., Koopmans, L., Jelić, V., et al. 2016, MNRAS, 462, 4482

Asad, K., Koopmans, L., Jelić, V., et al. 2018, MNRAS, 476, 3051

Barkana, R. 2018, Nature, 555, 71

Barkana, R., & Loeb, A. 2005, ApJ, 624, L65

Barry, N., Hazelton, B., Sullivan, I., Morales, M., & Pober, J. 2016, MNRAS, 461, 3135

Barry, N., Wilensky, M., Trott, C., et al. 2019, ApJ, 884, 1

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

Bevins, H. T., Heimersheim, S., Abril-Cabezas, I., et al. 2024, MNRAS, 527, 813

Bowman, J. D., Rogers, A. E., & Hewitt, J. N. 2008, ApJ, 676, 1

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

Bowman, J. D., Rogers, A. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018b, Nature, 564, E35

Broekema, P. C., Mol, J. J. D., Nijboer, R., et al. 2018, Astron. Comput., 23, 180

Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, MNRAS, 423, 2518

Chapman, E., Abdalla, F. B., Bobin, J., et al. 2013, MNRAS, 429, 165

Charrier, D., & the CODALEMA Collaboration 2007, in 2007 IEEE Antennas and Propagation Society International Symposium (IEEE), 4485

Chen, X., & Miralda-Escudé, J. 2004, ApJ, 602, 1

Cohen, A., Lane, W., Cotton, W., et al. 2007, AJ, 134, 1245

Condon, J. J., Cotton, W., Greisen, E., et al. 1998, AJ, 115, 1693

DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001

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

Eastwood, M. W., Anderson, M. M., Monroe, R. M., et al. 2019, AJ, 158, 84

Ewall-Wice, A., Dillon, J. S., Hewitt, J. N., et al. 2016, MNRAS, 460, 4320

Ewall-Wice, A., Chang, T.-C., Lazio, J., et al. 2018, ApJ, 868, 63

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

Fialkov, A., & Barkana, R. 2014, MNRAS, 445, 213

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

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

Field, G. B. 1958, Proc. IRE, 46, 240

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

Gan, H., Mertens, F., Koopmans, L., et al. 2023, A&A, 669, A20

Garsden, H., Greenhill, L., Bernardi, G., et al. 2021, MNRAS, 506, 5802

Gehlot, B. K., Koopmans, L. V., de Bruyn, A., et al. 2018, MNRAS, 478, 1484

Gehlot, B., Mertens, F., Koopmans, L., et al. 2019, MNRAS, 488, 4271

Gehlot, B., Mertens, F., Koopmans, L., et al. 2020, MNRAS, 499, 4158

Gehlot, B., Koopmans, L., Offringa, A., et al. 2022, A&A, 662, A97

Gessey-Jones, T., Sartorio, N. S., Fialkov, A., et al. 2022, MNRAS, 516, 841

Hibbard, J. J., Mirocha, J., Rapetti, D., et al. 2022, ApJ, 929, 151

Hicks, B. C., Paravastu-Dalal, N., Stewart, K. P., et al. 2012, PASP, 124, 1090

Hills, R., Kulkarni, G., Meerburg, P. D., & Puchwein, E. 2018, Nature, 564, E32

Hothi, I., Chapman, E., Pritchard, J. R., et al. 2021, MNRAS, 500, 2264

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

Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2008, MNRAS, 389, 1319

Joseph, R. 2012, Bachelor's Thesis, Faculty of Science and Engineering

Kazemi, S., Yatawatta, S., & Zaroubi, S. 2013, MNRAS, 430, 1457

Kern, N. S., & Liu, A. 2021, MNRAS, 501, 1463
Kingma, D. P., & Welling, M. 2013, arXiv e-prints [arXiv:1312.6114]

Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, arXiv e-prints
[arXiv:1505.07568]

Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133

Li, W., Pober, J., Barry, N., et al. 2019, ApJ, 887, 141

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

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429

McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815

Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2020, MNRAS, 493, 1217

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

Mertens, F. G., Mevius, M., Koopmans, L. V., et al. 2020, MNRAS, 493, 1662

Mertens, F. G., Semelin, B., Koopmans, L. V. E., et al. 2021, in SF2A-2021: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. A. Siebert, K. Baillié, E. Lagadec, et al., 211

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

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

Mevius, M. 2018, RMextract: Ionospheric Faraday Rotation calculator, Astrophysics Source Code Library [record ascl:1806.024]

Mevius, M., Mertens, F., Koopmans, L., et al. 2022, MNRAS, 509, 3693

Monsalve, R. A., Fialkov, A., Bowman, J. D., et al. 2019, ApJ, 875, 67

Morales, M. F., & Hewitt, J. 2004, ApJ, 615, 7

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

Munoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657

Murray, S. G., Greig, B., Mesinger, A., et al. 2020, J. Open Source Software, 5, 2582

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

Offringa, A., & Smirnov, O. 2017, MNRAS, 471, 301

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

Offringa, A., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606

Offringa, A., Mertens, F., Van der Tol, S., et al. 2019a, A&A, 631, A12

Offringa, A., Mertens, F., & Koopmans, L. 2019b, MNRAS, 484, 2866

Paciga, G., Chang, T.-C., Gupta, Y., et al. 2011, MNRAS, 413, 1174

Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639

Pandey, V., Koopmans, L., Tiesinga, E., Albers, W., & Koers, H. 2020, Astronomical Data Analysis Software and Systems XXIX, 527, 473

Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933

Patil, A. H., Zaroubi, S., Chapman, E., et al. 2014, MNRAS, 443, 1113

Patil, A. H., Yatawatta, S., Zaroubi, S., et al. 2016, MNRAS, 463, 4317

Patil, A., Yatawatta, S., Koopmans, L., et al. 2017, ApJ, 838, 65

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

Pritchard, J. R., & Furlanetto, S. R. 2007, MNRAS, 376, 1680

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

Reis, I., Fialkov, A., & Barkana, R. 2020, MNRAS, 499, 5993

Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, A&AS, 124, 259

Semelin, B., Mériot, R., Mertens, F., et al. 2023, A&A, 672, A162

Shaver, P., Windhorst, R., Madau, P., & De Bruyn, A. 1999, arXiv e-prints
[arXiv:astro-ph/9901320]

Singh, S., Subrahmanyan, R., Shankar, N. U., et al. 2017, ApJ, 845, L12

Singh, S., Nambissan, T., Subrahmanyan, R., et al. 2022, Nat. Astron., 1

Smirnov, O., & Tasse, C. 2015, MNRAS, 449, 2668

Subrahmanyan, R., Somashekar, R., Shankar, N. U., et al. 2021, Exp. Astron., 51, 193

Trott, C. M., Jordan, C., Midgley, S., et al. 2020, MNRAS, 493, 4711

van Diepen, G., Dijkema, T. J., & Offringa, A. 2018, Astrophysics Source Code Library [record ascl:1804.0003]

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

Vedantham, H., & Koopmans, L. 2015, MNRAS, 453, 925

Vedantham, H., & Koopmans, L. 2016, MNRAS, 458, 3099

Vedantham, H., Shankar, N. U., & Subrahmanyan, R. 2012, ApJ, 745, 176

Venkatesan, A., Giroux, M. L., & Shull, J. M. 2001, ApJ, 563, 1

Wouthuysen, S. 1952, AJ, 57, 31

Yatawatta, S. 2011, in 2011 XXXth URSI General Assembly and Scientific Symposium (IEEE), 1

Yoshiura, S., Pindor, B., Line, J., et al. 2021, MNRAS, 505, 4775

Zarka, P., Girard, J., Tagger, M., & Denis, L. 2012, in SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 687

Zarka, P., Tagger, M., Denis, L., et al. 2015, 2015 International Conference on Antenna Theory and Techniques (ICATT) (IEEE), 1

Zarka, P., Denis, L., Tagger, M., et al. 2020, in URSI GASS 2020
Kingma, D. P., & Welling, M. 2013, arXiv e-prints [arXiv:1312.6114]

Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, arXiv e-prints
[arXiv:1505.07568]

Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133

Li, W., Pober, J., Barry, N., et al. 2019, ApJ, 887, 141

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

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429

McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815

Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2020, MNRAS, 493, 1217

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

Mertens, F. G., Mevius, M., Koopmans, L. V., et al. 2020, MNRAS, 493, 1662

Mertens, F. G., Semelin, B., Koopmans, L. V. E., et al. 2021, in SF2A-2021: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. A. Siebert, K. Baillié, E. Lagadec, et al., 211

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

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

Mevius, M. 2018, RMextract: Ionospheric Faraday Rotation calculator, Astrophysics Source Code Library [record ascl:1806.024]

Mevius, M., Mertens, F., Koopmans, L., et al. 2022, MNRAS, 509, 3693

Monsalve, R. A., Fialkov, A., Bowman, J. D., et al. 2019, ApJ, 875, 67

Morales, M. F., & Hewitt, J. 2004, ApJ, 615, 7

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

Munoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657

Murray, S. G., Greig, B., Mesinger, A., et al. 2020, J. Open Source Software, 5, 2582

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

Offringa, A., & Smirnov, O. 2017, MNRAS, 471, 301

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

Offringa, A., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606

Offringa, A., Mertens, F., Van der Tol, S., et al. 2019a, A&A, 631, A12

Offringa, A., Mertens, F., & Koopmans, L. 2019b, MNRAS, 484, 2866

Paciga, G., Chang, T.-C., Gupta, Y., et al. 2011, MNRAS, 413, 1174

Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639

Pandey, V., Koopmans, L., Tiesinga, E., Albers, W., & Koers, H. 2020, Astronomical Data Analysis Software and Systems XXIX, 527, 473

Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933

Patil, A. H., Zaroubi, S., Chapman, E., et al. 2014, MNRAS, 443, 1113

Patil, A. H., Yatawatta, S., Zaroubi, S., et al. 2016, MNRAS, 463, 4317

Patil, A., Yatawatta, S., Koopmans, L., et al. 2017, ApJ, 838, 65

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

Pritchard, J. R., & Furlanetto, S. R. 2007, MNRAS, 376, 1680

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

Reis, I., Fialkov, A., & Barkana, R. 2020, MNRAS, 499, 5993

Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, A&AS, 124, 259

Semelin, B., Mériot, R., Mertens, F., et al. 2023, A&A, 672, A162

Shaver, P., Windhorst, R., Madau, P., & De Bruyn, A. 1999, arXiv e-prints
[arXiv:astro-ph/9901320]

Singh, S., Subrahmanyan, R., Shankar, N. U., et al. 2017, ApJ, 845, L12

Singh, S., Nambissan, T., Subrahmanyan, R., et al. 2022, Nat. Astron., 1

Smirnov, O., & Tasse, C. 2015, MNRAS, 449, 2668

Subrahmanyan, R., Somashekar, R., Shankar, N. U., et al. 2021, Exp. Astron., 51, 193

Trott, C. M., Jordan, C., Midgley, S., et al. 2020, MNRAS, 493, 4711

van Diepen, G., Dijkema, T. J., & Offringa, A. 2018, Astrophysics Source Code Library [record ascl:1804.0003]

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

Vedantham, H., & Koopmans, L. 2015, MNRAS, 453, 925

Vedantham, H., & Koopmans, L. 2016, MNRAS, 458, 3099

Vedantham, H., Shankar, N. U., & Subrahmanyan, R. 2012, ApJ, 745, 176

Venkatesan, A., Giroux, M. L., & Shull, J. M. 2001, ApJ, 563, 1

Wouthuysen, S. 1952, AJ, 57, 31

Yatawatta, S. 2011, in 2011 XXXth URSI General Assembly and Scientific Symposium (IEEE), 1

Yoshiura, S., Pindor, B., Line, J., et al. 2021, MNRAS, 505, 4775

Zarka, P., Girard, J., Tagger, M., & Denis, L. 2012, in SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 687

Zarka, P., Tagger, M., Denis, L., et al. 2015, 2015 International Conference on Antenna Theory and Techniques (ICATT) (IEEE), 1

Zarka, P., Denis, L., Tagger, M., et al. 2020, in URSI GASS 2020
# Appendix A: RFI statistics

<figure><img src="image_18.png" /><figcaption>**Fig. A.1.** RFI percentage at the L1 data level. The two panels show the flagging percentage as a function of time and frequency (left) and as a function of baseline (right).</figcaption></figure>

RFI flagging was performed at different stages of the data analysis pipeline. In Fig. A.1, we present the percentage of flagged data by AOFlagger until the L1 data level as a function of time, frequency, and baseline.

# Appendix B: Calibration solutions

## B.1. Bandpass solutions

Bandpass calibration is a crucial step to account for the frequency ripples produced in the data due to cable reflections. In NenuFAR the different MAs have different cable lengths, and the cable reflections are clearly visible in the autocorrelations before bandpass calibration. This is illustrated in Fig. B.1 where the delay power spectra of the autocorrelations of the NCP data before (middle panel) and after (right panel) bandpass calibration are plotted along with the delay power spectra of the bandpass solutions obtained separately from the Cas A observation (left panel), for a few example MAs. The vertical dashed lines indicate the expected delays due to the cable reflections for the MAs calculated from their respective cable lengths. We see a peak near the expected delay in both the bandpass solutions from the Cas A observation and the autocorrelations of the NCP data before bandpass calibration. The peaks are not present in the autocorrelations after the bandpass solutions have been applied to the data.

In the peak-normalized delay power spectra of the autocorrelations before and after bandpass calibration (middle and right panels in Fig. B.1), we find that even though the peaks due to the cable reflections are accounted for by bandpass calibration, the relative power at high delay values increases after bandpass calibration. To understand the impact of this effect on the power spectra, we performed a test in which the bandpass solutions were smoothed by a Gaussian filter of 2 MHz width, separately on the real and imaginary parts, before applying them to the data. The delay spectra of the autocorrelations with the smoothed bandpass show the peaks at the expected delays corresponding to the cable lengths but have a lower noise floor. However, in the final power spectra, we found negligible differences in the power levels of the autocorrelations with the actual bandpass against the autocorrelations with the smoothed bandpass. This is because the thermal noise level is well above the noise floor of the bandpass gain solutions. For the NCP field, the total flux

<figure><img src="image_19.png" /><figcaption>**Fig. B.1.** Bandpass calibration results for some example MAs. Left column: Delay power spectra of bandpass solutions obtained from the Cas A observation. Middle column: Delay power spectra of autocorrelations for the NCP observation before bandpass calibration. Right column: Delay power spectra of autocorrelations for the NCP observation after bandpass calibration. The cable lengths for the respective MAs are indicated by vertical dashed lines.</figcaption></figure>

is not high enough to give a dynamic range, which is needed for these effects to show up in the cylindrical and spherical power spectra.

## B.2. DD gain errors

We find that for the NCP subtraction step, a 6 MHz smoothing kernel results in better subtraction of sources compared to a 2 MHz smoothing kernel. Here, we compare the noise in the gain solutions in these two cases. The difference in the gain solutions at successive time steps is used as a proxy for the error in the gain solutions under the assumption that the gains do not change substantially over two solution intervals. Fig. B.2 shows the errors in the amplitude and phase of the DD gain solutions
# Appendix A: RFI statistics

<figure><img src="image_19.png" /><figcaption>**Fig. A.1.** RFI percentage at the L1 data level. The two panels show the flagging percentage as a function of time and frequency (left) and as a function of baseline (right).</figcaption></figure>

RFI flagging was performed at different stages of the data analysis pipeline. In Fig. A.1, we present the percentage of flagged data by AOFlagger until the L1 data level as a function of time, frequency, and baseline.

# Appendix B: Calibration solutions

## B.1. Bandpass solutions

Bandpass calibration is a crucial step to account for the frequency ripples produced in the data due to cable reflections. In NenuFAR the different MAs have different cable lengths, and the cable reflections are clearly visible in the autocorrelations before bandpass calibration. This is illustrated in Fig. B.1 where the delay power spectra of the autocorrelations of the NCP data before (middle panel) and after (right panel) bandpass calibration are plotted along with the delay power spectra of the bandpass solutions obtained separately from the Cas A observation (left panel), for a few example MAs. The vertical dashed lines indicate the expected delays due to the cable reflections for the MAs calculated from their respective cable lengths. We see a peak near the expected delay in both the bandpass solutions from the Cas A observation and the autocorrelations of the NCP data before bandpass calibration. The peaks are not present in the autocorrelations after the bandpass solutions have been applied to the data.

In the peak-normalized delay power spectra of the autocorrelations before and after bandpass calibration (middle and right panels in Fig. B.1), we find that even though the peaks due to the cable reflections are accounted for by bandpass calibration, the relative power at high delay values increases after bandpass calibration. To understand the impact of this effect on the power spectra, we performed a test in which the bandpass solutions were smoothed by a Gaussian filter of 2 MHz width, separately on the real and imaginary parts, before applying them to the data. The delay spectra of the autocorrelations with the smoothed bandpass show the peaks at the expected delays corresponding to the cable lengths but have a lower noise floor. However, in the final power spectra, we found negligible differences in the power levels of the autocorrelations with the actual bandpass against the autocorrelations with the smoothed bandpass. This is because the thermal noise level is well above the noise floor of the bandpass gain solutions. For the NCP field, the total flux

<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Bandpass calibration results for some example MAs. Left column: Delay power spectra of bandpass solutions obtained from the Cas A observation. Middle column: Delay power spectra of autocorrelations for the NCP observation before bandpass calibration. Right column: Delay power spectra of autocorrelations for the NCP observation after bandpass calibration. The cable lengths for the respective MAs are indicated by vertical dashed lines.</figcaption></figure>

is not high enough to give a dynamic range, which is needed for these effects to show up in the cylindrical and spherical power spectra.

## B.2. DD gain errors

We find that for the NCP subtraction step, a 6 MHz smoothing kernel results in better subtraction of sources compared to a 2 MHz smoothing kernel. Here, we compare the noise in the gain solutions in these two cases. The difference in the gain solutions at successive time steps is used as a proxy for the error in the gain solutions under the assumption that the gains do not change substantially over two solution intervals. Fig. B.2 shows the errors in the amplitude and phase of the DD gain solutions
## Appendix C: A-team subtraction trials

<figure><img src="image_20.png" /><figcaption>**Fig. B.2.** Gain calibration errors at the NCP subtraction step. The top and bottom panels show the errors in the amplitude and phase of the gain solutions, respectively. The dashed lines correspond to the solutions obtained using 2 MHz smoothing while the solid lines correspond to those obtained using 6 MHz smoothing. The thick solid lines show the beam factor and are to be read using the vertical axis on the right-hand side. The different colors refer to the gains in the direction of the different clusters.</figcaption></figure>

obtained in this manner, by taking the standard deviation of the time-differenced gains across frequencies and stations. Clusters 1 to 7 are arranged in the decreasing order of flux and we see that fainter clusters have higher calibration errors. For each cluster, the calibration errors for a 6 MHz smoothing kernel are smaller compared to a 2 MHz smoothing kernel. The modulation of the errors in the gain solutions toward Clusters 6 and 7 with time is due to the fact that the central lobe of the NenuFAR primary beam is not axisymmetric and is elongated in one direction. As the beam rotates with respect to the sky with time, Cluster 7 overlaps with the beam for the first few hours and then moves out of the beam. So the errors are lower in the first few hours. The opposite is true for Cluster 6 and hence we see the reverse trend in the gain errors. The average primary beam for each data segment at the center of Cluster 6 and Cluster 7 are plotted in Fig. B.2 and the anticorrelation of the gain calibration errors with the primary beam can be seen clearly as a function of time.

<figure><img src="image_21.png" /><figcaption>**Fig. C.1.** Standard deviation of the real part of the visibilities along frequency and baselines, as a function of time for different calibration trials. At each successive trial, only one setting was changed in order to understand its impact on the calibration. Trial 1 is the result of a fulljones calibration with Cyg A and Cas A subtracted when they are at elevations above 10 degrees. The subsequent changes in the other trials are described as follows. Trial 2: MAs 62, 64, 72 flagged, Trial 3: A-team included till 0 degree elevation, Trial 4: Tau A included, Trial 5: MA 6 flagged, Trial 6: MA 18 flagged, Trial 7: Vir A included, Trial 8: Calibration done in the diagonal mode.</figcaption></figure>

We varied the number of stations to be flagged, the number of A-team sources to be included, the elevation range for which an A-team is included, and the calibration mode (fulljones or diagonal) in these trial runs and converged to the optimized calibration parameters by monitoring the calibration solutions, data statistics, and images at the end of each trial. As an example, Fig. C.1 shows the standard deviation of the real parts of the XX and YY calibrated visibilities, integrated along all frequencies and baselines, for these different trial runs. We see that it improves progressively as we converge to the final calibration settings. Many more trials were performed with more MAs flagged (MAs 1, 9, 22, 30, 56, 59, and 63), but they did not lead to any further improvement in the solutions, statistics, or images; hence, these stations were not flagged in the final data used for the analysis.
Appendix C: A-team subtraction trials

<figure><img src="image_21.png" /><figcaption>**Fig. B.2.** Gain calibration errors at the NCP subtraction step. The top and bottom panels show the errors in the amplitude and phase of the gain solutions, respectively. The dashed lines correspond to the solutions obtained using 2 MHz smoothing while the solid lines correspond to those obtained using 6 MHz smoothing. The thick solid lines show the beam factor and are to be read using the vertical axis on the right-hand side. The different colors refer to the gains in the direction of the different clusters.</figcaption></figure>

obtained in this manner, by taking the standard deviation of the
time-differenced gains across frequencies and stations. Clusters
1 to 7 are arranged in the decreasing order of flux and we see that
fainter clusters have higher calibration errors. For each cluster,
the calibration errors for a 6 MHz smoothing kernel are smaller
compared to a 2 MHz smoothing kernel. The modulation of the
errors in the gain solutions toward Clusters 6 and 7 with time
is due to the fact that the central lobe of the NenuFAR primary
beam is not axisymmetric and is elongated in one direction. As
the beam rotates with respect to the sky with time, Cluster 7
overlaps with the beam for the first few hours and then moves
out of the beam. So the errors are lower in the first few hours.
The opposite is true for Cluster 6 and hence we see the reverse
trend in the gain errors. The average primary beam for each data
segment at the center of Cluster 6 and Cluster 7 are plotted in Fig.
B.2 and the anticorrelation of the gain calibration errors with the
primary beam can be seen clearly as a function of time.

<figure><img src="image_22.png" /><figcaption>**Fig. C.1.** Standard deviation of the real part of the visibilities along frequency and baselines, as a function of time for different calibration trials. At each successive trial, only one setting was changed in order to understand its impact on the calibration. Trial 1 is the result of a *fulljones* calibration with Cyg A and Cas A subtracted when they are at elevations above 10 degrees. The subsequent changes in the other trials are described as follows. Trial 2: MAs 62, 64, 72 flagged, Trial 3: A-team included till 0 degree elevation, Trial 4: Tau A included, Trial 5: MA 6 flagged, Trial 6: MA 18 flagged, Trial 7: Vir A included, Trial 8: Calibration done in the *diagonal* mode.</figcaption></figure>

We varied the number of stations to be flagged, the number of
A-team sources to be included, the elevation range for which
an A-team is included, and the calibration mode (**fulljones** or
**diagonal**) in these trial runs and converged to the optimized
calibration parameters by monitoring the calibration solutions,
data statistics, and images at the end of each trial. As an exam-
ple, Fig. C.1 shows the standard deviation of the real parts of
the XX and YY calibrated visibilities, integrated along all fre-
quencies and baselines, for these different trial runs. We see that
it improves progressively as we converge to the final calibra-
tion settings. Many more trials were performed with more MAs
flagged (MAs 1, 9, 22, 30, 56, 59, and 63), but they did not lead
to any further improvement in the solutions, statistics, or images;
hence, these stations were not flagged in the final data used for
the analysis.