# Ionospheric contributions to the excess power in high-redshift 21-cm power-spectrum observations with LOFAR

S.A Brackenhoff, M Mevius, L.V.E Koopmans, A Offringa, E Ceccotti, J.K Chege, B.K Gehlot, S Ghosh, C Höfer, F.G Mertens, et al.

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

## ► To cite this version:

S.A Brackenhoff, M Mevius, L.V.E Koopmans, A Offringa, E Ceccotti, et al. Ionospheric contributions to the excess power in high-redshift 21-cm power-spectrum observations with LOFAR. Monthly Notices of the Royal Astronomical Society, 2024, 533 (1), pp.632-656. 10.1093/mnras/stae1856. hal-04679348

HAL Id: hal-04679348

https://hal.science/hal-04679348v1

Submitted on 30 Aug 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Ionospheric contributions to the excess power in high-redshift 21-cm power-spectrum observations with LOFAR

S. A. Brackenhoff<sup>1,★</sup> M. Mevius<sup>2</sup> L. V. E. Koopmans<sup>1</sup> A. Offringa,<sup>1,2</sup> E. Cecotti<sup>1</sup> J. K. Chege,<sup>1</sup> B. K. Gehlot<sup>1</sup> S. Ghosh,<sup>1</sup> C. Höfer,<sup>1</sup> F. G. Mertens<sup>1,3</sup> S. Munshi<sup>1</sup> and S. Zaroubi<sup>1,4</sup>

Accepted 2024 July 29. Received 2024 July 29; in original form 2024 June 20

## ABSTRACT

The turbulent ionosphere causes phase shifts to incoming radio waves on a broad range of temporal and spatial scales. When an interferometer is not sufficiently calibrated for the direction-dependent ionospheric effects, the time-varying phase shifts can cause the signal to decorrelate. The ionosphere's influence over various spatiotemporal scales introduces a baseline-dependent effect on the interferometric array. We study the impact of baseline-dependent decorrelation on high-redshift observations with the Low Frequency Array (LOFAR). Data sets with a range of ionospheric corruptions are simulated using a thin-screen ionosphere model, and calibrated using the state-of-the-art LOFAR epoch of reionization pipeline. For the first time ever, we show the ionospheric impact on various stages of the calibration process including an analysis of the transfer of gain errors from longer to shorter baselines using realistic end-to-end simulations. We find that direction-dependent calibration for source subtraction leaves excess power of up to two orders of magnitude above the thermal noise at the largest spectral scales in the cylindrically averaged autopower spectrum under normal ionospheric conditions. However, we demonstrate that this excess power can be removed through Gaussian process regression, leaving no excess power above the 10 per cent level for a 5 km diffractive scale. We conclude that ionospheric errors, in the absence of interactions with other aggravating effects, do not constitute a dominant component in the excess power observed in LOFAR epoch of reionization observations of the North Celestial Pole. Future work should therefore focus on less spectrally smooth effects, such as beam modelling errors.

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

## 1 INTRODUCTION

Detecting the spatial variations in the redshifted 21-cm line from neutral hydrogen during the epoch of reionization (EoR) and cosmic dawn is an important objective in observational cosmology, because it offers a unique glimpse into the evolution of the early Universe, its radiating sources, and how they interact with their environment (e.g. Furlanetto, Peng Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012; Zaroubi 2013). Inhomogeneity in the density of neutral hydrogen in the early Universe causes the 21-cm signal to vary spatially. The aim of interferometric experiments is to capture these spatial variations in a power spectrum. However, achieving a detection of this high-redshift 21-cm field remains technically challenging for many reasons (see Liu & Shaw 2020, for example). A number of low-frequency arrays have been built to detect this 21-cm signal over the past two decades, or are currently under construction. Arrays that have set upper limits on the 21-cm power spectrum are the GMRT (Giant Metrewave Radio Telescope, Swarup

1990; Paciga et al. 2013),<fnref n="1" /> the MWA (Murchison Widefield Array, Bowman et al. 2013; Barry et al. 2019; Li et al. 2019; Trott et al. 2020),<fnref n="2" /> PAPER (Donald C. Backer Precision Array to Probe EoR, Parsons et al. 2010; Kolopanis et al. 2019),<fnref n="3" /> HERA (Hydrogen Epoch of Reionization Array, DeBoer et al. 2017; Abdurashidova et al. 2022),<fnref n="4" /> LOFAR (Low-Frequency Array, Haarlem et al. 2013; Patil et al. 2017; Mertens et al. 2020),<fnref n="5" /> OVRO-LWA (Owen's Valley Radio Observatory – Long Wavelength Array, Hallinan 2015; Eastwood et al. 2019; Garsden et al. 2021), and NenuFAR (New Extension in Nançay Upgrading LOFAR, Zarka et al. 2012; Munshi et al. 2024b).<fnref n="6" /> Furthermore, the construction of SKA (Square Kilometre Array, Dewdney et al. 2009; Koopmans et al. 2015)<fnref n="7" />, an instrument

<footnote n="1"><sup>1</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Netherlands Institute for Radio Astronomy (ASTRON), PO Box 2, 7990 AA Dwingeloo, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, PO Box 808, Ra'anana 4353701, Israel</footnote>

<footnote>* E-mail: brackenhoff@astro.rug.nl</footnote>

<footnote n="1">¹http://gmrt.ncra.tifr.res.in</footnote>

<footnote n="2">²http://www.mwatelescope.org</footnote>

<footnote n="3">³http://eor.berkeley.edu</footnote>

<footnote n="4">⁴http://reionization.org</footnote>

<footnote n="5">⁵http://www.lofar.org</footnote>

<footnote n="6">⁶https://nenufahr.obs-nancay.fr</footnote>

<footnote n="7">⁷http://www.skatelescope.org</footnote>
that will have a much higher sensitivity than its predecessors, has begun.

Galactic and extragalactic foreground sources that are several orders of magnitude brighter obscure the 21-cm signal (Shaver et al. 1999). However, the foregrounds are more spectrally smooth than the EoR signal, such that the two can be separated during data analysis. This requires high-precision calibration of the instrument. Issues can arise, for example, due to sky model errors and incompleteness (Ewall-Wice et al. 2017), or because the gain errors vary too fast as a function of time, direction, and baseline to be solved for. In addition, effects such as correlated radio-frequency interference (RFI, Offringa, Mertens & Koopmans 2019a) and instrumental leakage between Stokes components or cable reflections (Jelić et al. 2010; Gasperin et al. 2018; Offringa et al. 2019b) can leave residuals in the data after calibration. Finally, there are ionospheric errors that introduce chromatic errors (Koopmans 2010; Vedantham & Koopmans 2016). Combined, these effects lead to an 'excess variance' that biases the recovered 21-cm power spectrum (Mertens et al. 2020). The origins of this excess variance have remained hard to identify (see e.g. Gan et al. 2023). The goal of this work is to assess the part that the ionosphere plays in this excess variance.

Recent work using LOFAR and NenuFAR suggests that the excess variance may partially be attributed to residual foregrounds from bright off-axis sources (Ceccotti et al. in preparation) and to unflagged low-level RFI emitted near the interferometer (Munshi et al. 2024b). The removal of off-axis sources poses a challenge because the primary beam of the instrument is highly spatially, spectrally and to a lesser degree temporally varying far from the target direction (see e.g. Yatawatta et al. 2013; Cook, Seymour & Sokolowski 2021). Bright sources moving through side lobes in this variable beam cause an excess variance in the data, because of the inability to model the beam with sufficient accuracy (Asad et al. 2018). The difficulty of modelling this excess flux is further exacerbated by stochastic corruptions of the flux of these bright sources in the ionosphere (e.g. Koopmans 2010; Vedantham & Koopmans 2015). It is of great importance to identify the impact of each possible excess variance source on the power spectrum, such that they may be mitigated.

The ionosphere is a known source of difficulty in creating radioastronomical images, and many methods for correcting for the ionosphere have been investigated (see, for example, Intema et al. 2009; Weeren et al. 2016; Rioja, Dodson & Franzen 2018; Albert et al. 2020; Gasperin et al. 2020; Tasse et al. 2021; Chege et al. 2022). The ionosphere can affect the data in various ways, such as through its own thermal emission, refraction, and diffraction, and in its lower layers also through absorption. While thermal emission and absorption is of great concern in global 21-cm experiments (Sokolowski et al. 2015; Datta et al. 2016; Shen et al. 2022), an interferometer is more sensitive to refraction and diffraction (depending on the lengths of its baselines). In the context of the foreground-avoidance approach followed by HERA, work has been done on estimating the ionospheric attenuation on polarized foregrounds (Martinot et al. 2018). The MWA is an instrument that is relatively similar to LOFAR. However, during calibration, the MWA uses only its shorter baselines, and assumes that the effect of the ionosphere can be corrected with source position-dependent refractive shifts (Lonsdale 2005; Chege et al. 2021) if the most active ionospheric time-bins are discarded (Trott et al. 2018; Waszewski, Morgan & Jordan 2022). Work on MWA EoR measurements by Jordan et al. (2017) and Chege et al. (2022) has shown that the refractive effects of the ionosphere are not a dominant contributor to their excess variance and that updating the model based on ionospheric refraction does not improve results.

Because the LOFAR calibration strategy uses longer baselines for calibration, the variation of ionosphere across the array cannot be ignored. This can cause sources to scintillate as well as refract, which has not been simulated in the context of MWA.

Analytic work on the expected variance of measured visibilities has been done by Vedantham & Koopmans (2015). In a subsequent paper, Vedantham & Koopmans (2016) have computed the impact of scintillation noise on power spectra with an idealized source flux density distribution. However, the LOFAR EoR calibration strategy employs a method in which a longer subset of baselines is used for calibration and a shorter subset for power-spectrum analysis. This separation between baseline lengths is not considered in the work by Vedantham & Koopmans (2016).

In this work, we present simulations of ionospheric effects on LOFAR EoR measurements taken with the High-Band Array (HBA), in which we can control the behaviour of the ionosphere and separate it from other *in situ* effects. This provides a novel result, because the longer baselines of LOFAR cause us to be in a different ionospheric regime than earlier MWA-based work. Furthermore, we analyse whether ionospheric calibration errors from longer baselines are transferred to shorter baselines. We present, for the first time, a full simulation of ionospheric impact on the LOFAR EoR pipeline and assess the impact on the final power spectrum.

This work is structured as follows. First, we discuss the effect of the ionosphere on a minimally redundant array such as LOFAR, in Section 2. In Section 3, we describe our end-to-end simulation and calibration pipelines. In these, we first generate data with ionospheric corruptions and then process these data in a pipeline that replicates the LOFAR EoR data processing pipeline. In Section 4, we show the resulting power spectra, study the impact of the ionosphere on the gain solutions, and discuss how we can extrapolate the results to longer observations. In Section 5, we discuss the main conclusions of this work.

## 2 THE IONOSPHERE

The ionosphere is a partially ionized layer of turbulent plasma in the upper atmosphere. This layer causes electromagnetic propagation delays, such that incoming wavefronts are perturbed, which predominantly introduces phase errors (and in some cases amplitude errors). These perturbations vary in time, direction, and frequency and also depend on the point where the electromagnetic field is measured on the ground. The errors are difficult to calibrate for due to their rapid time decorrelation. This section describes the effects of the ionosphere on radio-interferometric data. We focus on phase errors only and leave higher order effects, such as differential Faraday rotation and the three-dimensional nature of the ionosphere, out of the scope of this work, because they are expected to affect the data to a much lesser extent (Gasperin et al. 2018).

### 2.1 Frequency behaviour

The propagation delays and resulting phase errors are caused by the free electron density $N_e$ in the column of the ionosphere that the incoming wavefront traverses. This is quantified in the Total Electron Content (TEC), measured in TEC Units ($1 \text{ TECU} = 10^{16} \text{ m}^{-2}$),

$$ \text{TEC} \approx \int N_e(s) \, ds \, [\text{m}^{-2}]. \qquad (1) $$

In this equation, $s$ denotes the line of sight along which the TEC is measured. For LOFAR HBA frequencies, the phase delay $\Delta\phi$ in
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Projected baseline density of a LOFAR observation (left axis) and phase structure function (right axis) for both a diffractive scale of 5 and 10 km. The baseline lengths are binned in 40 logarithmically spaced bins between 20 m and 20 km. The baseline lengths of 50, 250, and 5000 λ are also indicated for a frequency of 140 MHz.</figcaption></figure>

radians at frequency $v$ can be computed using

$$ \Delta\phi = \frac{e^2}{2\pi vcm_e\epsilon_0} \left[ \frac{\text{TEC}}{\text{TECU}} \right] \approx 8.4 \times 10^9 \cdot \left[ \frac{\text{TEC}}{\text{TECU}} \right] \frac{1}{v} \quad (2) $$

Here, $e$ is the electron charge, $v$ the frequency of the incident wave in Hz, $c$ the speed of light in vacuum, $m_e$ the rest mass of an electron, and $\epsilon_0$ the vacuum permittivity. The phase delay can be predicted across the entire waveband if the TEC is known.

## 2.2 Spatiotemporal behaviour

In addition to its spectral variation, the TEC varies spatially over a range of length-scales and in time. In this paper, we use the frozen screen approximation, where the differential TEC that causes ionospheric turbulence is modelled as an unchanging screen that moves over the telescope (Taylor 1938). This movement causes the differential TEC observed by the stationary array to evolve with time. The statistical properties of the phase errors are modelled by the phase structure function $D_\phi$, which describes the phase variance between two piercing points through the ionospheric screen as a function of the distance between them (Narayan 1992; van der Tol 2009; Vedantham & Koopmans 2016). A commonly used phase structure function is that of a Kolmogorov screen of the form

$$ D_{\phi}(r) = \left(\frac{r}{r_{\text{diff}}}\right)^{\beta} \quad (3) $$

Here, $D_\phi$ is given in rad$^2$, $r$ is the scale being probed, and $r_{\text{diff}}$ the diffractive scale, defined as the scale at which $D_\phi(r) = 1$ rad$^2$. We use 150 MHz as the frequency for which we define the diffractive scale.<fnref n="8" /> A smaller diffractive scale means more phase variance on all scales (see the slanted lines in Fig. 1). A turbulent ionosphere is described by this equation using the Kolmogorov index $\beta$ of 5/3 (Rufenach 1972).

Observations show that this structure function closely matches the ionospheric perturbations experienced by LOFAR (Mevius et al.

2016; Gasperin et al. 2019; Gan et al. 2022; Waszewski et al. 2022). A notable difference is that the true index $\beta$ has been shown to vary somewhat, reaching levels of nearly 2 (Mevius et al. 2016). This index indicates a more structured ionosphere, for example, due to travelling ionospheric disturbances (TIDs) (Wandzura 1980; Loi et al. 2016). Since the effects of large-scale TIDs are easier to calibrate for (Gasperin et al. 2018), we choose to adopt the more commonly used value of the index of $\beta = 5/3$. Some ionospheric models also include an inner and an outer scale in their phase structure functions, where the screen behaves differently from a pure Kolmogorov screen. The inner scale is of the order of metres (Booker 1979; Vedantham & Koopmans 2015). The outer scale at 150 MHz was observationally found to be at least ~80 km (Mevius et al. 2016). Like in real LOFAR EoR observations, we use baselines between 50λ and 5000λ for calibration and power spectrum analysis, corresponding to physical lengths between ~100 m and ~11.2 km in the simulated frequency band. These baseline lengths are also indicated in Fig. 1. Since all of our baselines are sufficiently longer than the inner scale and shorter than the outer scale, these scales are omitted in the description used here. Finally, we do not include simulations of ionospheric ducts or other anisotropies, such as those discussed by Trott et al. (2018).

## 2.3 The role of baseline length in sensitivity to ionospheric errors

We use the thin-screen approximation and the pierce-point approximation to model ionospheric effects, because the thickness of the ionospheric layer is a second-order effect. The thin-screen approximation means that we do not model the full three-dimensional structure of the ionosphere, but instead simulate a two-dimensional screen. The ionosphere is simulated using the TEC at the pierce-points and scaled based on the angle under which incident radiation would travel through the ionospheric slab (see Edler, de Gasperin & Rafferty 2021). A pierce-point is defined as the point at which the line of sight of a station observing in a direction crosses the thin ionospheric screen. The ionospheric phase difference that a baseline experiences in a direction is given by the difference in TEC value between its pierce points. Under these assumptions, the pierce-point separation is the projected baseline length on to the phase screen.

According to the phase structure function in equation (3), a longer baseline experiences a larger phase variance. The impact of the largest variations on those baselines is mitigated by the fact that long baselines can be better calibrated for ionospheric disturbances due to their large time coherence.<fnref n="9" />

Lonsdale (2005) defined four ionospheric regimes, based on the size of the array and its field of view (FoV) compared to ionospheric scales. LOFAR is in the regime where both its FoV and its longest baselines are large compared to the scales of the ionospheric irregularities (van der Tol, Jeffs & van der Veen 2007). Its FoV means that the phase delay varies over the FoV, and therefore direction-dependent (DD) calibration solutions are needed to correct for the ionosphere. Furthermore, because the stations are far apart geographically, such a set of DD solutions is needed for each station separately.

Vedantham & Koopmans (2015) show that for baselines shorter than the Fresnel scale $r_F$, and for a source in the target direction, the decorrelation time $t_d$ becomes $t_d = 2r_F/v$, where $v$ is the bulk wind velocity. For baselines longer than the Fresnel scale, the decorrelation time of a baseline of physical length $b$ varies between $t_d = 2b/v$ and

<footnote n="8">⁸This frequency was chosen since it is the customary reference frequency for LOFAR EoR observations. Therefore, using 150 MHz makes it easier to compare our results to previous work, and we choose it despite it being outside of our simulation band (see Table 1).</footnote>

<footnote n="9">⁹These long baselines are also perturbed by the smaller scale variations between their pierce points, however.</footnote>
$t_d = 4b/v$, depending on the angle between the baseline and the wind direction.<fnref n="10" /> We use a frozen screen moving at 100 m s⁻¹. For baselines shorter than the Fresnel scale of ~300 m (which occurs for a thin-screen height of 300 km), this implies a decorrelation time of ~6 s. The longest baselines used in this work are up to ~11 km (5000 λ). The decorrelation time of the largest ionospheric scales they are sensitive to can go up to a few minutes, but they are also influenced by the more rapidly decorrelating small-scale variations. Solution intervals during gain calibration, on the other hand, range from 30 s for initial calibration up to 20 min for DD calibration in the LOFAR EoR project (see Mertens et al. 2020). Hence, rapid small-scale variations within a solution interval are not corrected for on these baselines, leading to phase errors on the data (Koopmans 2010). This is especially true for the small-scale ionospheric variations, where the solution intervals far exceed the coherence time and hence, all ionospheric errors remain in the data after calibration. This also holds for any other instrument, such as MWA, HERA, and SKA. Some 'speckle noise', i.e. a halo-structure caused by the ionosphere around point-like sources, may also occur (Koopmans 2010). This is often not seen in low-frequency images, because this speckle noise is strongly baseline-dependent and often below the thermal or confusion noise. Furthermore, it is often dominated by the decorrelation of the point spread function (PSF) of the array (Vedantham & Koopmans 2015). However, in 21-cm experiments, it may show up as additional excess variance in the power spectrum due to its high sensitivity. These errors on the station-based gain solutions are not straightforward to predict and, furthermore, are strongly direction-dependent. Each station is part of baselines of various lengths and therefore experiences a range of wave-modes and coherence time-scales simultaneously.

## 2.4 Impact on the power spectrum

Gan et al. (2022) have correlated observed excess variance with ionospheric activity parameters as measured over 3-h time bins. No clear correlation between these parameters and the observed excess variance was shown. A stronger correlation between the excess variance and local sidereal time (LST) was found instead. However, this does not exclude the possibility that the ionosphere still plays a major role in the excess variance for two reasons. First of all, ionospheric errors make the removal of bright foregrounds more challenging. This implies that errors in foreground subtraction would also correlate with the on-sky position of these sources and therefore LST. This would mean an ionospheric effect may still show a predominant correlation with local time. Furthermore, the activity level of the ionosphere is difficult to measure with LOFAR, because short bursts of high ionospheric activity can dominate ionospheric variance estimates. A time bin may therefore have data that are only occasionally strongly affected by the ionosphere, whereas most data are only weakly affected by ionospheric errors. In this work, we assume that the ionospheric behaviour remains ergodic over time, such that there are no temporal changes in the phase-structure functions and therefore no ionospheric activity bursts.

The expected impact of ionospheric scintillation on the LOFAR 21-cm power spectrum has been estimated analytically by Vedantham & Koopmans (2016). They have shown that for an instrument similar to LOFAR, in the case of a simplified sky power spectrum and calibration strategy, the excess variance above thermal noise follows

a wedge-like structure that is dependent on the ionospheric activity level. This structure leaves a part of the power spectrum (often called the 'EoR window') mostly unaffected. However, a full analysis of a realistic sky and calibration pipeline was outside of the scope of their work. Ewall-Wice et al. (2017) have shown that modelling errors can cause foreground leakage into the EoR window. Since ionospheric errors cause distortions of foreground sources and therefore have a similar effect to modelling errors, this supports the idea that ionospheric errors may leak into other parts of the power spectrum.

DD gain calibration in the LOFAR EoR pipeline (see Patil et al. 2017; Mertens et al. 2020) employs a calibration step where the full set of baselines is divided into a group of longer and a group of shorter baselines. The long baselines are used to calibrate gains. These gains are then applied to the sky model and this model is subtracted on shorter baselines. This avoids signal suppression on shorter baselines (Patil et al. 2016; Mevius et al. 2022), but may transfer gain-calibration errors due to the ionosphere from the longer to the shorter baselines, if the calibration intervals on the long baselines are longer than their coherence time. Although this transfer of gain errors has been investigated for incomplete sky models and beam errors (Barry et al. 2016; Ewall-Wice et al. 2017; Mouri Sardarabadi & Koopmans 2019), to our knowledge the impact of the transfer of ionospheric errors from longer to shorter baselines during calibration, and the impact of speckle noise, has not been investigated in detail via simulation but only estimated theoretically.

# 3 SIMULATION AND CALIBRATION PIPELINES

To gauge the effects of the ionosphere on LOFAR HBA data and data processing, we create realistic mock data sets via forward simulations, both including and excluding ionospheric errors, and calibrate those using a pipeline nearly identical to the state-of-the-art LOFAR EoR standard pipeline (Mertens et al. in preparation) that is described in Appendix B (from here on referred to as the 'standard pipeline'). The pipeline for simulation is described in Section 3.1 and the pipeline for calibration and power spectrum estimation in Section 3.2.

## 3.1 Simulation pipeline

To isolate the impact of the ionosphere from other sources of errors and systematics, our forward simulations contain a limited number of effects. The included components are a sky model, ionospheric phase-screen, the station beam models, and thermal noise contributions. The simulations are limited to 12 h, a full rotational synthesis. Because the 21-cm signal is much weaker than the thermal noise level in observations of this duration, there is no need to include a 21-cm signal in the simulations. The simulation process is described below and shown schematically in Fig. 2.

### 3.1.1 Sky model

The sky model is a representation of the North Celestial Pole (NCP) field at 145 MHz. The NCP field has been used to produce upper limits on the EoR signal with LOFAR (Patil et al. 2017; Mertens et al. 2020) and NenuFAR (Munshi et al. 2024b). Vedantham & Koopmans (2015) show that the square of the scintillation noise of a set of point sources is roughly equal to the sum of squares of the scintillation noises of the individual sources. As such, the scintillation noise due to all foregrounds can be well approximated by only simulating the

<footnote n="10">¹⁰This holds for a frozen screen model. The decorrelation time is longer if the baseline is oriented in the direction of the motion of the ionosphere.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Pipeline used to create the simulated visibilities. The input sky model is denoted with a pentagon pointing leftwards. Processing steps are shown with a rounded rectangle. The title of the step is shown in bold font, the used software package is in mono-spaced font, and any further details are noted in regular font. The circles with a '+' sign denote an addition in visibility space. The pentagons pointing to the right denote output simulated visibilities saved in the MS.</figcaption></figure>

brightest sources. We limit the sky model to a restricted model of 684 unpolarized flat-spectrum sources, rather than using the 28 778 source components model used in the calibration of real LOFAR EoR data by Mertens et al. (2020). The restricted model was produced with WSCLEAN (Offringa et al. 2014; Offringa & Smirnov 2017)<fnref n="11" /> on a small LOFAR data set with a 10 × 10° FoV. Clean-components within the scale of the PSF were merged into single point sources. We set the intrinsic flux scale of the model using the predicted apparent flux of ~100 sources. The restricted model and full model agree to a 1 per cent level in apparent flux. The cylindrically averaged power spectra of the models agree to within 10 per cent, sufficient to assess the impact of the ionosphere in a realistic scenario.

Because the spatial coherence scale of the ionosphere is roughly 3.5 arcmin at the frequencies of interest on the baselines used to create the power spectra (Vedantham & Koopmans 2016), one would ideally compute a gain solution for each patch of this spatial scale on the sky. This is not possible, however, because this would create too many degrees of freedom during calibration and therefore increase the risk of 21-cm signal suppression. Furthermore, the signal-to-noise ratios of the patches would become too low to reach accurate gain solutions (Kazemi, Yatawatta & Zaroubi 2013). Instead, the 684 sources in our model are divided into only 24 clusters for gain calibration. This is a factor of 4–5 fewer than used in the standard pipeline, because we do not solve for the instrument beam during calibration and because the coherence time of the ionosphere on the power-spectrum baselines is much smaller than the calibration time interval. The bright source 3C 61.1 is located near the first beam null in the NCP field. The apparent brightness of this source varies strongly, because the station beam rotates around the NCP. Due to its distinctive role in calibration strategy (see Section 3.2.2), we assign 3C 61.1 to its own cluster<fnref n="12" /> and divide the remaining sources into 23 clusters based on their positions and apparent flux using Voronoi tesselation with LSMTOOL.<fnref n="13" /> Our clustered model is shown in Fig. 3 and has been included as supplementary material to this paper. An excerpt of the model is shown in Appendix A.

Due to its extreme apparent brightness, the bright far-field source Cassiopeia A (Cas A) has been included in the simulations.<fnref n="14" /> While Cygnus A (Cyg A) is also taken into consideration when the standard pipeline is applied to real data, it has a peak apparent brightness that

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic representation of the sky model used in this work. Each combination of shape and colour represents a calibration cluster for which a single gain solution is inferred and applied during sky-model subtraction. Due to its large angular separation from the NCP, Cas A is not shown in this figure, although it is included in the simulations.</figcaption></figure>

is an order of magnitude lower than Cas A in the simulated LST range. Its impact on the power spectrum, which scales with the flux squared, is therefore much lower. See, for example, Cook, Trott & Line (2022), who show the impact of unsubtracted far-field sources on the power spectrum. In their figs 6 and 7, they show the removal of a single bright source, the effect of which dominates over other, weaker off-axis sources through most of the power spectrum. In a similar way, Cas A dominates over the far-field foregrounds, and we use it as a representative source to describe the general behaviour of bright far-field sources on the power spectrum and omit Cyg A.

### 3.1.2 Measurement set

Measurement sets (MS) are created using the LOFAR Simulation Tool (LOSITO).<fnref n="15" /> An MS contains properties of the telescope array, such as station positions, uv-coverage, and frequency channels. The

<footnote n="11">¹¹W-Stacking CLEAN, https://wsclean.readthedocs.io/en/latest/index.html</footnote>

<footnote n="12">¹²This was done based on source position alone, such that the actual cluster consists of a few sources as can be seen in Fig. 3.</footnote>

<footnote n="13">¹³LOFAR Sky Model Tool, https://lsmtool.readthedocs.io/en/latest</footnote>

<footnote n="14">¹⁴Using the LOFAR Low-Band Array model taken from https://github.com/lofar-astron/prefactor/tree/master/skymodels</footnote>

<footnote n="15">¹⁵https://losito.readthedocs.io</footnote>
MS is adjusted to a specified phase centre (target direction) and observation time. We simulate the Dutch HBA stations used in the EoR project (Patil et al. 2017; Mertens et al. 2020) with a realistic pattern of inactive tiles (groups of 4×4 antennas) within the station. Two types of stations are included: core stations (CS), which are concentrated within a few kilometres diameter area near the LOFAR core, and remote stations (RS), which predominantly provide longer baselines. RS have twice the number of tiles that CS have and are beamformed to closely match the station layout of CS by assigning zero weights to a sub-set of the tiles. In real data, baselines sharing a common electronics cabinet lead to crosstalk and are therefore removed, so they are removed in the simulation as well.

A time resolution of 2 s is used. This is the highest time resolution data stored for LOFAR EoR observations and captures temporal variations at the same scale as available in measured data. Furthermore, this time interval is shorter than the temporal coherence at the Fresnel scale of ~6 s (Vedantham & Koopmans 2015), such that it is smaller than the shortest ionospheric time-scales. The frequency resolution is set to 195 kHz, the same spacing as the smallest spectral gain calibration intervals used in real data. As described in equation (2), the spectral behaviour of the ionosphere is smooth, such that a higher spectral resolution is not needed for this simulation. Because we aim to assess the impact of the ionosphere, there is no need to simulate spectral effects introduced by the receiver chain of the telescope. The total simulated bandwidth is 12.6 MHz, between 134.1 and 146.7 MHz (the same as used for $z = 9.1$ observations). We apply a baseline filter removing baselines of > 20 km as they are much longer than the longest baselines used during data processing (~ 11 km). Due to the distribution of stations, this leaves our longest baseline at 18 km, reducing the size of the TEC-screen needed.

### 3.1.3 TEC screen

We use the TEC method in LOSITO as described in Edler et al. (2021) to model the ionosphere. A frozen two-dimensional TEC screen is created for a given diffractive scale at 150 MHz at a given height (300 km in this work). This is done using a combination of a large-scale and small-scale screen as proposed by Buscher (2016), such that both scales can be computed in an efficient manner. Subsequently, pierce-points are computed for each source and station separately. The TEC values at these pierce-points are found through interpolation. A geometric factor is applied to the TEC value at each pierce-point to correct for the angle at which the ionosphere is pierced. The geometric factor takes into account both the angle of incidence and the curvature of the Earth. Hence, lines of sight near the horizon experience a higher TEC and therefore smaller diffractive scale than points near the zenith. This can lead to bright off-axis sources dominating the ionospheric errors when they are low in the sky as shown by Gehlot et al. (2018) and is only mitigated by beam attenuation.

### 3.1.4 Thermal noise

Thermal noise is characterized by the system equivalent flux density (SEFD). The used SEFD is 4.2 kJy, based on the estimated SEFD for LOFAR EoR observations listed by Mertens et al. (2020). Thermal noise is generated by drawing the real and imaginary component of the noise component of the visibility separately from a Gaussian distribution with a standard deviation of

<table><caption>**Table 1.** Properties of the simulated MS.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Telescope</td><td>LOFAR HBA</td></tr><tr><td>Pointing</td><td>RA = 0h00m00s, Dec. = +90°00′00″</td></tr><tr><td>Sky model</td><td>684 sources in the central 10 × 10° Cassiopeia A</td></tr><tr><td>Bandwidth</td><td>134.1–146.7 MHz</td></tr><tr><td>Frequency resolution</td><td>195 kHz</td></tr><tr><td>Time resolution</td><td>2 s</td></tr><tr><td>Total duration</td><td>12 h</td></tr><tr><td>System equivalent flux density</td><td>4.157 kJy</td></tr><tr><td>Baseline range</td><td>42 m–18 km</td></tr><tr><td>Height of ionospheric layer</td><td>300 km</td></tr></tbody></table>

$$ \sigma_{\text{thermal}} = \frac{\text{SEFD}}{\sqrt{2\Delta t \Delta v}} \text{ [Jy].} \quad (4) $$

Here, $\Delta t$ and $\Delta v$ represent the temporal and spectral resolution of the visibilities in the MS as listed in Table 1.

The thermal noise is stored separately, such that the same noise realization can be added to several simulated data sets. This allows for a better comparison between calibration results for different simulations while still including the effect of noise.

### 3.1.5 Visibility prediction

Forward prediction of the visibilities is done using DP3 (van Diepen, Dijkema & Offringa 2018).<fnref n="16" /> LOSITO has been designed to interface with DP3, such that the TEC screen can easily be applied in a direction-dependent way during the prediction step. DP3 allows for control over the simulation process, such that we can choose not to include temporal and spectral smearing in our simulations and therefore use a lower spectral resolution for the simulations (see Table 1). DP3 includes a realistic model of the LOFAR beam including sidelobes. This is important for simulating the behaviour of Cas A, which can pass through such sidelobes during the observations. We produce three data products:

(i) TRUE VISIBILITIES – This set of visibilities is a forward prediction of the sky model with LOFAR beam effects and contains no noise or ionospheric effects. It is considered the ‘ideal case’ and can be compared to other data products to quantify their ionospheric errors.

(ii) REFERENCE VISIBILITIES – This set is created by adding a thermal noise realization to the TRUE VISIBILITIES. We use this set to assess the impact of an ionosphere-free to a full model, which includes ionospheric effects.

(iii) FULL VISIBILITIES – This set is generated by adding phase perturbations to the predicted visibilities of each source separately using the TEC-screen described in Section 3.1.3 to create DD ionospheric effects. Furthermore, the same thermal noise realization is added as to the REFERENCE VISIBILITIES.

## 3.2 Calibration and power spectrum estimation pipeline

The calibration and power spectrum estimation pipeline largely follows the LOFAR EoR standard pipeline (Mertens et al. in preparation, see Appendix B). Because our input data are simplified compared to

<footnote n="16">¹⁶Default Pre-Processing Pipeline</footnote>
observational LOFAR EoR data (such as the omission of a 21-cm signal, bandpass errors, RFI, etc.), some parts of the standard pipeline can be omitted, resulting in the pipeline shown in Fig. 4. However, there are a few differences between this pipeline and the standard pipeline:

(i) The calibration pipeline presented in this paper utilizes the DDECal (Gan et al. 2023) submodule of DP3 instead of SAGECAL-CO (Yatawatta 2015). As mentioned in Section 3.1.5, there are advantages to using DP3 for simulating the data. However, because SAGECAL-CO and DP3 use slightly different beam models, we prefer to use DP3 in both simulation and calibration, as it allows us to use the exact same beam model in both. Therefore, any ambiguity between beam modelling mismatches and ionospheric effects is avoided. As shown by Gan et al. (2023), calibration pipelines using DDECal and SAGECAL-CO produce comparable power spectra, such that the impact of this difference in calibration algorithm will minimally affect the resulting power spectra. Additionally, differential experiments have shown that the difference in calibration algorithm is a second-order effect in the calibration result. The differences in the calibration results between SAGECAL-CO and DDECal are a factor $10^{-7}$ smaller than the source powers (Gehlot, private communications).

(ii) Assigning solution intervals that differ between source clusters is difficult, because of technical differences between DDECal and SAGECAL-CO.<fnref n="17" /> Therefore, we use the two-step method presented by Gan et al. (2023), where bright off-axis sources are subtracted at a higher time resolution in a separate step. This step is followed by a lower time resolution subtraction of the weaker sources near the target direction.

(iii) Radio frequency interference flagging and bandpass calibration steps are omitted, because these effects are not present in the simulated data.

(iv) We only consider a single night of observation at a time for power spectrum estimation.

### 3.2.1 Initial calibration

The settings in all gain calibration steps are listed in Table 2. In the first gain-calibration step, a common gain across the field is computed for each station. In the case of the NCP field, this round of calibration is complicated by the bright source 3C 61.1 near the first beam-null, which can cause time- and direction-dependent effects during calibration. Initial calibration is therefore done in two directions, such that 3C 61.1 can be isolated. The gains computed in the direction of 3C 61.1 are discarded, and the gain solutions computed for the rest of the field are applied to the full data set. Cas A is not included in any of the initial calibration, because its apparent brightness is highly temporally variable due to beam sidelobes far from the target direction.

Initial calibration is performed on a solution interval of 30 s, for baselines ranging from 50 to 5000 λ. We use the spectral smoothness constraint, with a width of the Gaussian smoothing kernel of 4 MHz recommended by Gan et al. (2023). Because the sky model is unpolarized, instrumental polarization errors are not introduced, and dispersive effects of the ionosphere produce a scalar

gain perturbation.<fnref n="18" /> Only gain solutions on the diagonal of the Jones matrices are computed.

### 3.2.2 DD calibration

The DD calibration step removes the sources in the sky model from the data. Gains are computed in the directions of the sources, and then the contributions to the visibilities from each source are predicted and subtracted using those gains. Before DD gain calibration, the data are first averaged to a time resolution of 10 s. Only the longer baselines of length > 250 λ are used during this calibration step to avoid subtraction of the 21-cm signal of interest during the calibration process (Patil et al. 2016).

We first subtract Cas A in a separate calibration step with a time resolution of 2.5 min. To do this, the sky model is divided into two clusters: a cluster for Cas A and a cluster containing all other sources. After calibration, the gains computed for Cas A are used to subtract it from the data, whereas the other gains are discarded. Subsequently, a second DD-calibration step is performed to remove the remaining 684 sources near the target direction. We perform DD-calibration towards 24 clusters of sources on a 10-min time resolution. This process is summarized in the right grey box in Fig. 4.

### 3.2.3 A-team flagging

After DD gain-calibration and sky model subtraction, the dynamic range in the data is drastically reduced and the data can be gridded in uv-space. Because the baselines migrate through uv-space as a function of time and frequency and multiple baselines are gridded into a single uv-cell, this step mixes residual post-calibration artefacts (e.g. from the ionosphere) between different baselines. Therefore, these artefacts can no longer be removed after gridding. In real data, a few bright A-team sources still create large residuals in the power spectra at this point in the pipeline. They must be flagged before further analysis in a step called 'uv-space masking' or 'A-team flagging'. This is done by creating a mask in uv-space where the A-team sources are expected to be brightest (see Munshi et al. 2024a). A line in the direction of each A-team source on the uv-plane is calculated, and all gridded visibilities within 2 pixels of this line are masked. This is illustrated in Fig. 5. Although the A-team flagging step removes the strongest effects of these sources, not all of their impact is removed.

### 3.2.4 Residual foreground removal

After DD-source subtraction, grinding, and A-team flagging, residual foregrounds can remain in the data. The spectral variability of the foregrounds is increased, due to a process called 'mode-mixing' (Morales et al. 2012), where the chromatic nature of the interferometer mixes spatial fluctuations into the spectral domain, such that they obscure the much fainter 21-cm signal. As a result, the foregrounds dominate over the 21-cm signal in a larger region of the power spectrum than just the most spectrally smooth modes (Datta, Bowman & Carilli 2010; Vedantham, Udaya Shankar & Subrahmanyan 2012). In addition, chromatic effects due to the ionosphere (Koopmans 2010; Vedantham & Koopmans 2016; Trott et al. 2018) and instrument beam can cause power from the foreground sources

<footnote n="17">¹⁷This is because DDECAL performs gain calibration by iterating over stations rather than calibration directions. While DDECal does offer the functionality of DD calibration intervals, it comes at the cost of sacrificing precision and robustness by treating all directions independently within the update step.</footnote>

<footnote n="18">¹⁸Higher order effects such as differential Faraday rotation are not simulated and therefore do not need to be corrected for.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Pipeline used to calibrate the simulated visibilities. The input visibilities simulated using the pipeline shown in Fig. 2 are denoted with a pentagon pointing to the left. Processing steps are shown with a rounded rectangle. The title of the step is shown in bold font, the used software package in mono-spaced font, and any further details are noted in regular font. The circles with a '-' sign denote a subtraction in visibility space. The pentagons pointing to the right denote output power spectra shown in Section 4.</figcaption></figure>

<table><caption>**Table 2.** Settings for initial and DD calibration.</caption><thead><tr><th>Property</th><th>Value in initial calibration</th><th>Value in DD-calibration for Cas A</th><th>Value in DD-calibration for near phase-centre sources</th></tr></thead><tbody><tr><td>Operation</td><td>Apply gains</td><td>Subtract with gains</td><td>Subtract with gains</td></tr><tr><td>Minimum baseline (λ)</td><td>50</td><td>250</td><td>250</td></tr><tr><td>Maximum baseline (λ)</td><td>5000</td><td>5000</td><td>5000</td></tr><tr><td>Number of directions</td><td>2<sup>a</sup></td><td>2<sup>b</sup></td><td>24</td></tr><tr><td>Solution time interval (s)</td><td>30</td><td>150</td><td>600</td></tr><tr><td>Solution frequency interval (kHz)</td><td>195</td><td>195</td><td>195</td></tr><tr><td>Number of iterations</td><td>100</td><td>200</td><td>200</td></tr><tr><td>Jones matrix elements calibrated</td><td>Diagonal</td><td>Diagonal</td><td>Diagonal</td></tr><tr><td>Smoothness Constraint</td><td>4 MHz</td><td>4 MHz</td><td>4 MHz</td></tr></tbody></table>

Notes. <sup>a</sup>3C 61.1 is placed in one cluster, the remainder of the central model is placed in the other.

<sup>b</sup>Cas A is placed in one cluster, the full central model in the other.

to leak to other modes. Because of this, these residual foregrounds need to be removed in order to recover the 21-cm signal.

To remove the residual foregrounds from the data, Mertens, Ghosh & Koopmans (2018) introduced a method based on Gaussian process regression (GPR). Gaussian processes are used to model the components in the data through fitting frequency covariances (kernels) as a function of uv-coordinates. These components can describe a specific process (such as the 21-cm signal, mode-mixing, thermal noise, etc.), but also more general residual foregrounds or systematics (such as excess variance). The covariance kernels are described by a general functional form that depends on the frequency separation and several hyperparameters. These hyperparameters can describe e.g. the power normalization or spectral coherence scale. Markov Chain Monte Carlo or nested sampling is used to obtain samples from the posterior probability distribution of the hyperparameters given the observed (or simulated) data. To compute a power spectrum of a component, a number of realizations are drawn from the posterior distribution for its kernel. Those realizations are used to compute a probability distribution of the power spectrum of the component.

Its power spectrum is calculated from realizations drawn from this distribution.

In real data, the residual foregrounds are comprised of foreground sources not included in the model (such as diffuse emission or extragalactic emission near the confusion noise limit), and foreground sources that are not well subtracted (for example due to calibration errors). We expect to need fewer and slightly different components to model the data in our simulations than needed in real data. This is because our data contain no 21-cm signal, beam errors, or sky-model incompleteness. Fewer foreground components may therefore be needed to describe the residual foregrounds after DD-subtraction. However, foreground contamination is still expected, due to ionospheric variations on time-scales shorter than the calibration intervals which result in source subtraction errors.

In this paper, we have therefore chosen a simplified set of kernels that subtract the foregrounds to a satisfactory level. This set might not be optimal or unique but demonstrates the behaviour of ionospheric errors in GPR. The used GPR components are the intrinsic foregrounds, mode-mixing contaminants, and uncorrelated
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Illustration of uv-space masking. The coloured pixels show the SEFD estimated from the spectrally differenced simulated visibilities of the central model and Cas A. A line in the direction of Cas A is shown in black and the pixels highlighted in red represent the pixels discarded during (u, v)-space flagging.</figcaption></figure>

thermal noise. The intrinsic foreground kernel is intended to capture smooth emission that has not been accurately removed during sky-model subtraction inside the image field used for gridding. Similar to Mertens et al. (2020), we use a Median Radial Basis Function (MRBF) with a large length-scale (150 MHz, much larger than the total bandwidth of the simulations) to model this emission. The mode-mixing kernel describes spectrally rapidly fluctuating effects introduced by the instrument, such as the imprint of the PSF and beam chromaticity on residual foregrounds both inside and outside the gridding image. The latter is of a lesser concern in this work, due to the limited angular extent of our sky model. Mode-mixing is modelled with a 5/2 Matérn covariance kernel, which is narrower than the MRBF kernel. Finally, thermal noise is modelled to be uncorrelated in frequency.

A noise estimate is made through the Stokes V temporal differences in the visibilities. The noise variance is scaled with respect to this estimate. All other variances are normalized to the data variance. The hyperparameters are listed in Table 3. A nested sampler with 500 live points is used to optimize the evidence for each simulation.

### 3.2.5 Power spectra

Due to the extreme sensitivity needed to detect the 21-cm signal, several modes are averaged in Fourier space in order to create a suitably sensitive power spectrum. A version of this power spectrum often used to gauge the effectiveness of data analysis pipelines is the cylindrically averaged power spectrum. In this spectrum, the Fourier dual of the frequency channels of an observation are used on one axis ($k_\parallel$). The other axis ($k_\perp$) averages over the Fourier dual of spatial variations on the sky. As such, the pixels of the cylindrically averaged power spectrum represent annuli in Fourier space. Ideally, the spectrally smooth foregrounds would be confined to the lowest $k_\parallel$-modes, leaving the remainder of the power spectrum free for 21-cm signal analysis. However, the aforementioned effect of 'mode-mixing' causes the foregrounds to dominate over the 21-cm signal in a larger wedge-shaped region of the power spectrum, called the 'foreground-wedge'.

The power spectra in this work are computed using the code PSPIPE<fnref n="19" /> with the conventional definitions by Morales & Hewitt (2004). Only the shortest baselines (50–250λ) which are the most sensitive to the 21-cm signal are used for this. A power spectrum is produced at each calibration step, to investigate how the ionosphere impacts that step.

Assessing errors made during the DD-calibration step can be difficult, because the DD-calibration step also includes sky-model subtraction. This means that, even if calibration errors are introduced, the foreground power is typically reduced by orders of magnitude. Therefore, we have differenced the FULL VISIBILITIES and TRUE VISIBILITIES, both before and after initial calibration of the FULL VISIBILITIES, and created power spectra of the residual as well. Although this cannot be done with real data, it provides a fairer comparison of the ionospheric effects on the data.

## 4 RESULTS AND DISCUSSION

In this section, we analyse the effects of the ionosphere on the simulations described in Section 4.1. We first inspect how ionospheric errors appear in the cylindrically averaged power spectrum in different steps of the calibration and data processing pipeline, in Section 4.2. The impact of the level of ionospheric turbulence (i.e. different diffractive scales) is compared in Section 4.3. A deeper inspection of the effects A-team flagging is performed in Section 4.4. The effects of phase and amplitude calibration, respectively, are analysed in Section 4.5. Finally, the correlation between nights is inspected in Section 4.6, to investigate whether ionospheric errors average down when multiple nights of data are integrated.

### 4.1 Simulations

Five simulations created using the simulation pipeline described Section 3.1 are summarized in Table 4. The first two simulations ($r_{10a}$ and $r_{10b}$) assume a diffractive scale of 10 km, and are the most representative of LOFAR EoR data (see Mevius et al. 2016; Gan et al. 2022). These simulations are used to investigate the coherence of ionospheric residuals in observations at different LST. A different LST means that the baseline orientations and station beams compared to the sky will be different. However, there may still be some coherence between the Fourier modes probed within a gridded uv-bin across multiple nights, such that the residuals of strong sources remain correlated. As such, these simulations make it possible to estimate whether the residuals are incoherent across multiple nights and can therefore be suppressed through processing more data.

Two simulations represent an active ionosphere and have a diffractive scale of 5 km. This scale is the limit at which LOFAR EoR data are considered to be of sufficient quality. One ($r_{5b}$) of these two simulations has a unique TEC screen, and one ($r_{5a}$) has the TEC screen used in one of the 10 km diffractive scale simulations ($r_{10a}$), with all TEC values scaled up to a diffractive scale of 5 km. This allows us to differentiate between effects due to stronger ionospheric activity for the exact same screen and between different random realizations.

Finally, we have performed a simulation that does include thermal noise, but no ionospheric errors (Null). This enables differentiation between effects due to the ionosphere and effects due to data processing, calibration (e.g. gain-solution time-scale), sky-model

<footnote n="19">¹⁹Power-Spectrum PIPEline, https://gitlab.com/flomertens/pspipe</footnote>
<table><caption>**Table 3.** List of used hyperparameters in GPR. In this table *U* and *LU* indicate a uniform and log-uniform prior, respectively. This is a simplified set of kernels in comparison to the kernels used in the standard pipeline, because many effects present in real data have been omitted from our simulations to isolate ionospheric effects. More details on these kernels can be found in Mertens et al. (2018) and Mertens, Bobin & Carucci (2024).</caption><thead><tr><th>Component</th><th>Kernel</th><th>Hyperparameter</th><th>Prior</th></tr></thead><tbody><tr><td>Intrinsic foregrounds</td><td>MRBF</td><td>Data normalized variance</td><td>U(10<sup>-4</sup>, 10<sup>2</sup>)</td></tr><tr><td></td><td></td><td>Length-scale</td><td>Fixed at 150</td></tr><tr><td>Mode-mixing foreground</td><td>5/2 Matérn kernel</td><td>Data normalized variance</td><td>U(10<sup>-4</sup>, 10<sup>2</sup>)</td></tr><tr><td></td><td></td><td>Length-scale</td><td>U(0.1, 20)</td></tr><tr><td>Noise</td><td>Uncorrelated</td><td>Scaling with respect to Stokes V time difference</td><td>U(0.75, 1.3)</td></tr></tbody></table>

<table><caption>**Table 4.** Overview of the simulated data used in this work. A diffractive scale of 5 km is the cut-off for real data, below which data are not used for processing due to its high distortion level.</caption><thead><tr><th>Label</th><th>r<sub>diff</sub></th><th>Start time (UTC)</th><th>TEC realization</th><th>Noise realization</th></tr></thead><tbody><tr><td>r10a</td><td>10 km</td><td>2013-11-27 18:00:00</td><td>Nominal model</td><td>Nominal model</td></tr><tr><td>r10b</td><td>10 km</td><td>2014-03-06 18:00:00</td><td>Unique</td><td>Unique</td></tr><tr><td>r5a</td><td>5 km</td><td>2013-11-27 18:00:00</td><td>Scaled TEC-values of r10a</td><td>Same as r10a</td></tr><tr><td>r5b</td><td>5 km</td><td>2013-11-27 18:00:00</td><td>Unique</td><td>Same as r10a</td></tr><tr><td>Null</td><td>&infin;</td><td>2013-11-27 18:00:00</td><td>No ionospheric errors</td><td>Same as r10a</td></tr></tbody></table>

subtraction, and GPR (e.g. chosen kernels). All simulations that share a start time also share the same thermal noise realization. This minimizes ambiguity in whether different results in the various simulations can be attributed to the ionospheric or thermal noise realizations. To create a realistic second 'night' of simulation, the simulation for another start time (r10b) has been given a unique thermal noise realization.

## 4.2 The propagation of ionospheric errors

To assess how ionospheric errors in the calibration process evolve, first, we analyse the errors in the nominal simulation r10a separately and compare it to the ionosphere-free case in the Null simulation.

### 4.2.1 Simulated data and initial calibration

Fig. 6 shows (residual) images made using simulation r10a. The top row shows images after initial calibration. The bright source near the first null, 3C 61.1, causes the ring structures coming from the bottom-left corner in the initial-calibrated target field image. The overall behaviour of the simulated data in these images resembles that of real data.<fnref n="20" />

Fig. 7 shows the residual power spectra created at different stages of the calibration process for simulations r10a and Null (we present the same figure for the other simulations in Appendix C). Note that the power spectra shown here are all residual power spectra, i.e. the TRUE VISIBILITIES have been subtracted from the raw and initial-calibrated data, in order to show the ionospheric errors rather than the source power. The ratio between the power spectrum of the simulations and their injected noise is presented in the third and fourth rows. They show the 'excess' variance as any residual power above the thermal noise level. The bottom row, which shows the ratio between the power spectrum of simulations r10a and Null, discriminates between systematic errors due to the ionosphere, and any other possible systematic errors.

Up to the delay line of Cas A (the source farthest away from the phase centre), the residuals follow the wedge and are typically a few times the thermal noise (third row, panels one and two). We expect the residuals to be coherent as a function of frequency, such that they can be removed in the GPR step. There is a slight excess in the ratio with thermal noise after initial-calibration above the wedge (the red hue starting just above Cas A and extending into the EoR window in Fig. 7 columns 2 panel 3). This is due to the application of the gains from initial calibration, which change the visibilities and therefore also this region of the power spectrum. The ratio then changes, because the thermal noise we are comparing to is not calibrated.

In Fig. 8, changes in the cylindrically averaged power spectrum throughout the calibration pipeline are presented. We show the ratio of power spectra extracted at each subsequent stage of the pipeline. This figure highlights where power is added or removed and therefore we use the FULL VISIBILITIES directly for these plots (rather than creating a residual like in Fig. 7). In the leftmost panel, the ratio between un-calibrated and initial-calibrated data is shown. This ratio does not appreciably deviate from unity anywhere in the power spectrum. Large deviations would mainly occur for large phase variations that are coherent across the full FoV (tilt corrections). However, the ionospheric phase variations in the simulation are small, because only baselines up to 250λ are used in the power spectra. Therefore, the changes between the raw and initial-calibrated data being minimal indicates that no major systematic errors or corrections are introduced at this point in the pipeline.

### 4.2.2 DD calibration

The Cas A model is subtracted in the next step of the pipeline, followed by the NCP source model. The residuals in the image plane can be seen in the bottom two panels of Fig. 6. While DD-calibration can remove most of the source flux, some power remains visible near Cas A in the wide field image. However, the impact of these residuals on the power spectrum is best estimated in the central 4 × 4° FoV, because this is where the power spectrum is estimated from. When comparing the variance of this part of the image to an image made with the same simulation, but with Cas A omitted from both the simulation and calibration, we find a ~0.5 per cent change in variance in the full-band image and a ~0.1 per cent change in

<footnote n="20">²⁰See for example the right side of fig. 2 in Mertens et al. (2020), which is plotted at the same angular scale as our figure.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Images created using the full nominal simulation r10a with WSCLEAN. The top panels are created with initial-calibrated data and the bottom panels with DD-source subtracted data. The left panels represent a 68° × 68° FoV with a 7 arcmin resolution made using baselines between 50 and 1000λ. The right panels are a zoomed-in version showing only the central 4° × 4°, with a 42 arcsec resolution and baselines between 50 and 5000 λ. A few important sources are labelled. Power near the NCP and around Cas A is clearly visible in the initial-calibrated wide field image (top left). The flux of both Cas A and the 684 source model is reduced after DD-calibration, although residual power is still visible.</figcaption></figure>

variance in a single-channel image. We conclude that the impact of Cas A is limited on images in this FoV. In the target field image (bottom-right panel), the brightest sources leave a ring-like structure behind. This could either be an artefact of the PSF convolved with the residual of the source, or due to source scintillation creating a halo-like structure around the original source position directly (Koopmans 2010; Vedantham & Koopmans 2016).

The effect of removing Cas A and the NCP sources on the power spectra is shown separately in Fig. 8. During the Cas A-subtraction step, power is removed mostly around the delay region of Cas A, with a slight extension above it at low $k_⊥$. Removal of the NCP sources results in power reduction at low $k_‖$, extending above the expected delay of the central sky model. This cannot be due to the spectral dependence of ionospheric phase errors on the short baselines directly, as such power should be visible in Fig. 8. Therefore, we attribute this to the transfer of errors from longer baselines to shorter ones, which introduces stronger spectral variations in the source power, extending the foreground power to higher $k_‖$-modes. Any transfer of errors happens along constant delay

lines up to the maximum delay causing a so-called ‘brick’ in the power spectrum (Ewall-Wice et al. 2017). Additionally, there is a vertical feature at $k_⊥$ = 0.06 in the third panel of Fig. 8 (a similar feature is present in Figs C2 and C3), attributed to low uv-coverage in this part of the spectrum. Baselines at this separation have been removed because $k_⊥$ = 0.06 $hcMpc^{-1}$ corresponds to the separation between a pair of stations that share an electronics cabinet (see Ewall-Wice et al. 2017, Section 3.2). There is a gap in the uv-coverage around this $k_⊥$, because we remove these baselines in the simulation, similar to how we treat real data.

Additional power introduced during the DD-calibration steps are more clearly seen comparing the initial-calibrated and DD-source subtracted residuals in Fig. 7. Because the TRUE VISIBILITIES have been subtracted from the visibility data shown in the second column of this figure, the change between these columns represents errors introduced by imperfect gain calibration during the DD-step. Note that this is not the same as a comparison between gains that perfectly compensate for the ionosphere and gains that do not, because the TRUE VISIBILITIES ignore the ionosphere and
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Overview of the residual power in cylindrical power spectra of nominal simulation r10a and the ionosphere-free simulation Null. The columns represent the different calibration steps. From left to right: before any calibration, after initial calibration, after DD calibration with source subtraction, and after GPR. The TRUE VISIBILITIES have been subtracted from the data in the first two columns, such that all columns represent residuals, even before the source subtraction step in the pipeline. From top to bottom, the rows represent the power spectra for r10a and Null, respectively, the ratio between these power spectra and a thermal noise power spectrum, and the ratio between the power spectrum of simulation r10a and that of simulation Null. In the first three columns of the middle row, the thermal noise is the injected thermal noise. In the last column, this is replaced by a thermal noise cube drawn from GPR (as described in Section 3.2.4). In this cylindrically averaged power spectrum and in all such power spectra throughout this work we plot four delay lines, from top to bottom: the horizon line, two lines defining the delay area where the peak power from Cas A should lie, and the delay below which the peak power of a source at 5° from the phase centre should lie. These lines are based on improved delay calculations that take the sky curvature and phase-referencing away from zenith into account by Munshi et al. (2024a).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The ratio between cylindrically averaged power spectra between calibration steps for simulation r10a to illustrate the impact of each calibration step. From left to right, the columns represent the ratio between raw data and initial-calibrated data, the ratio between initial-calibrated data and the same data after subtraction of Cas A, the ratio between the data after Cas A has been subtracted and after the central sky model has been subtracted, the ratio between this residual and the same data with Cas A masked out in A-team flagging, and the ratio between this data product and the same after residual foreground removal.</figcaption></figure>

<table><caption>**Table 5.** Mean and maximum power spectrum ratios between simulations and the thermal noise or Null simulation, respectively. The shown simulations all share the same LST as the Null simulation. The lowest level subcolumn represents the point in the pipeline; after initial calibration, after DD-source subtraction, and after GPR, respectively. For the ratio between GPR and thermal noise, two ratios are shown: the first is the ratio with the noise cube drawn from GPR (the same as used in the noise ratios in the fourth column of Fig. 7), the second is the ratio with the input thermal noise realization (as in the other thermal noise ratios). This second value suffers from sampling variance after GPR, making the maximum more prone to outliers. We therefore refer to the first ratio when discussing excess variance, and provide the second one for completeness.</caption><thead><tr><th rowspan="3">Label</th><th colspan="8">Ratio with thermal noise</th><th colspan="8">Ratio with Null simulation</th></tr><tr><th colspan="4">Maximum</th><th colspan="4">Mean</th><th colspan="4">Maximum</th><th colspan="4">Mean</th></tr><tr><th>Initial</th><th>DD</th><th>GPR<br>versus<br>estimate</th><th>GPR<br>versus<br>input</th><th>Initial</th><th>DD</th><th>GPR<br>versus<br>estimate</th><th>GPR<br>versus<br>input</th><th>Initial</th><th>DD</th><th>GPR</th><th>Initial</th><th>DD</th><th>GPR</th></tr></thead><tbody><tr><td>r10a</td><td>12.6</td><td>22.3</td><td>1.0</td><td>1.9</td><td>1.2</td><td>1.6</td><td>1.0</td><td>1.0</td><td>13.0</td><td>13.2</td><td>1.4</td><td>1.2</td><td>1.4</td><td>1.0</td></tr><tr><td>r5a</td><td>39.8</td><td>454.5</td><td>1.1</td><td>2.7</td><td>1.7</td><td>9.7</td><td>1.0</td><td>1.1</td><td>41.2</td><td>237.3</td><td>2.6</td><td>1.7</td><td>6.4</td><td>1.1</td></tr><tr><td>r5b</td><td>30.7</td><td>801.2</td><td>1.0</td><td>2.1</td><td>1.7</td><td>16.9</td><td>1.0</td><td>1.1</td><td>32.5</td><td>558.1</td><td>1.8</td><td>1.7</td><td>11.3</td><td>1.1</td></tr></tbody></table>

assume gains of unity. Power in the wedge is more pronounced in the power spectra of the residuals after DD-calibration (compared to the idealized foreground subtraction in the previous column), both at the delays of the NCP source and at the delay lines of Cas A. In the power spectrum itself (the first panel in the third column of Fig. 7), this is mainly visible at the lowest k-modes and for Cas A around $k_{\parallel} = 0.2 \text{ hCpc}^{-1}$. When inspecting the ratio with thermal noise in the middle panel in the third column of Fig. 7, the excess of Cas A is strongly pronounced. Interestingly, however, the excess power near the delay lines of Cas A is weak in the ratio between simulations Null and r10a (bottom row). The same effect is clear in Table 5, where the mean and maximum values of the ratio between r10a and its thermal noise is higher than the ratio between r10a and Null. Because the residuals of Cas A are stronger if there are ionospheric errors, but do not disappear if no ionospheric errors are present, we conclude that the increase in excess power at this stage in the data processing is linked to ionospheric errors, but cannot be completely attributed to them.

### 4.2.3 A-team flagging and residual foreground removal

During the final data processing steps, the uv-cells where bright off-axis sources are expected to be strongest are masked and GPR model-fitting is applied to remove residual foregrounds present in the data after calibration and model subtraction. We show the posterior distributions of the GPR hyperparameters in Figs 9 and 10 (the posteriors of the remaining distributions are shown in

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Corner plot of the posterior probability distributions of GPR on the Null simulation. This figure illustrates how the components are fitted in GPR. The vertical dashed lines in the marginalized distributions indicate the 68 per cent confidence level.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The same as Fig. 9, but with simulation r10a shown in blue and simulation r10b shown in red. Despite having the same ionospheric activity level, the fitted kernels are slightly different between simulations r10a and r10b due to different ionospheric and thermal noise realizations and due to the LST difference between the simulations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Cylindrically averaged power spectrum representation of the fitted covariance kernels produced by drawing 1000 power spectra from each component to illustrate where they capture power. These realizations were generated using the hyperparameters found for simulation r10a.</figcaption></figure>

Fig. D1 in Appendix D). The cylindrically averaged power spectra of the different components fitted for simulation r10a are shown in Fig. 11.

The masking of Cas A only has a very minor effect in the fourth panel of Fig. 8, indicating that the source has been removed well during DD calibration. GPR has removed power near the target direction, as is evident in the lowest $k_\parallel$-modes. This is in line with expectations from the foreground kernels visualized in Fig. 11. Most notably, we do not see an appreciable excess variance due to ionospheric errors in the final power spectrum of simulation r10a. Residual foreground removal effectively removes ionospheric errors from the cylindrically averaged power spectra to the thermal noise level. This is also evident from Fig. 7 and Table 5, where no clear excess over thermal noise is observed anymore after GPR. There are some minor effects visible in the wedge region in the comparison between simulations r10a and Null. However, these are small and

randomly distributed around unity, possibly speckle noise, indicating an increased variance in the power spectrum estimate, rather than a systematic bias.

Recall that the noise scaling hyperparameter is defined compared to the noise estimate from Stokes V time difference, and we therefore expect the noise power to be unity. However, there is a positive bias in noise scaling of the posteriors on the level of a few per cent. This issue persists when the priors are chosen differently, and when all variances are estimated in linear space rather than log-space. Therefore, the bias indicates that there is another power component being absorbed into the noise kernel. Because there is no clear systematic difference between the noise scaling in the Null simulation and the other four simulations (that do contain ionospheric effects), we do not expect this bias to be an ionospheric effect. We do not investigate it further here, because the resulting error is well below the level where it would significantly influence results in the analysis of real 21-cm data, and because the bias does not appear to result from the ionosphere.

When comparing the GPR components, thermal noise is dominant outside of the lowest $k_\parallel$-modes. The foreground residuals are mostly captured by the mode-mixing kernel, but also to a lower level by the intrinsic foreground kernel. This is the decomposition of the ionospherically corrupted residuals in a spectrally smooth and more rapidly fluctuating component. That these are ionospheric effects is supported by the results for the Null simulation, where the variance of the intrinsic foreground kernel hits the lower bound of its prior and the mode-mixing component is much weaker than for the simulations affected by the ionosphere. The absence of ionospheric errors makes source subtraction more effective, hence fewer residual foregrounds are present before GPR. The spectral coherence scales found for the mode-mixing kernel are low: around or below 3 MHz. Mertens et al. (2018) use 3 MHz as the limit below which some signal loss is expected, if the 21-cm signal itself is either not a part of the GPR model or incorporated in an approximate manner. The newly developed ML-GPR method proposed by Mertens et al. (2024) and tested by Acharya et al. (2024) has an improved kernel for modelling the 21-cm component, such that it can be distinguished from spectrally varying foreground power more easily and the 3 MHz cut-off has become less stringent. As such, the low spectral coherence scale may not lead to 21-cm signal subtraction in real data.

Although the current set of kernels is able to remove the ionospheric residuals to the noise level, two kernels are required to fully fit the ionospheric signature. For real data, this may pose a challenge in the presence of other errors, which also need to be removed by these kernels. Furthermore, the mode-mixing spectral coherence scale drops to a low spectral width, which could potentially lead to some signal loss. These issues may be remedied by creating a more physically motivated kernel for either instrument-induced mode-mixing or ionospheric effects.

The residual power spectra clearly differ from those obtained using LOFAR observations, as described by Mertens et al. (2020), both in shape and amplitude. For a similar length of observation and after similar data processing, the foreground wedge is still clearly visible in real data, whereas it can be removed in our simulations. This is despite real data having a similar or larger diffractive scale than the simulations, such that ionospheric effects should be smaller in real data. Furthermore, observed data also show an increase in power for the lowest $k_\perp$-modes, which we do not reproduce in our simulations. We therefore exclude ionospheric errors as the dominant cause of the excess variance found in LOFAR EoR observations under normal ionospheric conditions of $r_{\text{diff}} = 10$ km.
## 4.3 Ionospheric conditions

Because simulations r10a and r5a are created using the same thermal noise realization and TEC screen, be it scaled in amplitude in the latter case, the impact of the level of ionospheric distortions can be assessed without being affected by differences in the noise and ionospheric sample variance between the two simulations. A diffractive scale of 5 km leads to an increase in residual power during calibration for two reasons. First of all, a smaller diffractive scale means more ionospheric phase errors on the short baselines on time-scales shorter than the solution time interval for the gains, which is directly reflected in the residuals in the power spectrum. In conjunction with this, larger ionospheric variations on large spatial scales reduce the decorability of long baselines and are therefore detrimental to the gain solutions which are applied to the shorter baselines as well. However, even for a diffractive scale of 5 km, complete decorrelation is not reached on the majority of baselines, because they are below the diffractive scale.<fnref n="21" />

In Fig. 12, we present the ratios of the power spectra of r10a and r5a at different steps in the calibration pipeline. In the raw and initial-calibrated data, the effect of the larger ionospheric errors is visible in the foreground wedge, most prominently in the lower $k_\parallel$-modes. The difference between the two data sets becomes more pronounced after the DD-calibration step, as is expected due to the calibration limitations outlined above. The excess power at this point in the pipeline is about an order of magnitude higher for simulation r5a compared to simulation r10a, as is shown in Table 5. This table also shows the same ratios for simulation r5b, which has a unique ionospheric TEC-screen. Although the evolution of excess variance of r5a and r5b through the pipeline follows a similar trend and the excess is clearly higher than in r10a, their excess variances are significantly different. We attribute the difference to sampling variance, but do conclude that a lower diffractive scale leads to larger foreground residuals. This demonstrates that the impact of the ionospheric activity level is a non-linear effect that can drastically affect the gains, motivating the need for an ionospheric activity cut-off.

In the rightmost column of Fig. 12, slightly more residual power is present near the delay line of Cas A for the data set with a more active ionosphere. This is expected to be a result of the non-optimal (i.e. hard flagging) uv-plane masking of Cas A. The uv-plane has been masked in the same way in both data sets, but because we only mask in a 2-pixel width around the direction of Cas A, whereas its true footprint exists throughout the uv-plane, some residuals of Cas A can remain in the data after (u, v)-space flagging. Furthermore, because GPR only removes power at the delay lines near the target direction rather than in the area of Cas A, it is not able to remove the residual of Cas A and a larger residual before (u, v)-space flagging also results in a larger residual after. This is true for both $r_{\text{diff}} = 5$ km simulations (as can be seen in the third panels of the fourth columns in Figs C2 and C3). However, this residual has very low power and is only visible here because the simulations have an identical thermal noise realization. In real data, we do not probe down to the level of this residual within a single night, which can be seen from Table 5. Therefore, based on these results ionospheric errors can even be

removed in the worst ionospheric conditions used for data processing ($r_{\text{diff}} = 5$ km), with excess of at most 10 per cent.

## 4.4 Necessity of A-team flagging

To investigate how well GPR can capture the residuals of Cas A, we test the effects of uv-space masking. We use the same gridded visibilities but do not perform (u, v)-space flagging before applying GPR for one case while using the regular pipeline including uv-space masking for the other. The ratios between the results for simulation r5a are shown in Fig. 13. We use this simulation rather than r10a, because the residuals of Cas A are stronger such that the effect of flagging is more visible.

On the left panel of Fig. 13, there is no clear residual of Cas A visible, suggesting that the error is not dominating the error budget before flagging. On the right panel, flagging removes power in two regions: one near the expected delay of Cas A, and the other in a small region near $k_\perp = 0.06$ hcMpc$^{-1}$ at high $k_\parallel$. We attribute the power reduction in the latter region to having very few baselines at these $k_\perp$-values. Flagging the uv-cells at which Cas A is expected to dominate may lead to a significant loss in the number of visibilities going into these uv-cells, leading to the loss of power here. We do not expect this to be an ionospheric effect, because the region is affected on a similar level for simulation r10a, which has weaker ionospheric errors than the shown simulation r5a.

In the region where Cas A is expected, the change in residual power upon flagging does depend on the ionospheric activity. We conclude that at least for some ionospheric activity levels, far-field source removal is not completely effective, even if both the sky and the beam are perfectly modelled. We expect such far-field sources to be more difficult to remove during GPR model-fitting as the frequency-dependent beam sidelobes and nulls in the far-field imprint a spectral structure on the source during the observation (see e.g. Cook et al. 2021). GPR mainly removes spectrally smooth power, to avoid subtracting the 21-cm signal, such that these residuals are not removed. We expect the same to be true for other A-team sources, because they are also bright and far from the target direction. However, as shown by Gan et al. (2022), Cas A and Cyg A are by far the most dominant far-field sources in LOFAR HBA observations of the NCP, and therefore dominate the off-axis foreground power (Cook et al. 2022). We analyse the spectral impact of the beam on bright far-field sources in more detail in future work.

## 4.5 Amplitude and phase error transfer

One of the main motivations for this work is to assess whether ionospheric errors may lead to gain errors on long baselines during calibration, resulting in an error transfer to the short baselines used for power spectrum estimation. Because the solution intervals (being 30 s for the initial calibration, 10 min for the Cas A DD-calibration, and 2.5 min for the phase centre DD-calibration step) are longer than the simulation time resolution, a decorrelation in the visibilities occurs. This results in a change in both the amplitude and phase of the inferred gains. To investigate this further, we analyse the estimated gains from initial and DD-calibration on a baseline level. For stations i and j with station gains $g_i$ and $g_j$, we define the baseline gain $g_i g_j^*$, where $g^*$ indicates the complex conjugate of g. Each pair of stations is used only once, such that if baseline i, j is included, baseline j, i is omitted. The differences between the ionospheric effects and gains for all baselines at 140.3 MHz (the central frequency in our simulation band) are visualized in Fig. 14 for a single polarization. Two cases are shown here: off-axis source

<footnote n="21">²¹The longest baselines do experience a higher level of decorrelation, because both the NCP sources and Cas A are positioned away from the zenith for at least part of the observation. Due to the increased airmass away from the zenith, the experienced diffractive scale is lower than that set for the TEC screen.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The ratio between cylindrically averaged residuals (computed in the same way as in Fig. 7) of simulations r10a and r5a. These simulations share the same ionospheric phase screen, barring a re-scaling of all TEC values, such that the diffractive scale of simulation r5a simulation is half of that of simulation r10a. Their ratio illustrates the impact of the activity level of the ionosphere on each step of data processing. A green hue indicates more residual power in the simulation with a larger ionospheric variance. Note that outside the wedge, little to no difference is observed between the two power spectra as they share a thermal noise realization.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** The effect of flagging Cas A in simulation r5a. Left: the SEFD computed using the spectrally differenced residuals after subtraction of all sources. Masked pixels are highlighted in red. Right: the ratio of a cylindrically averaged power spectrum with uv-space masking of Cas A and a similar power spectrum without masking Cas A. Both power spectra used in this ratio are a residual after GPR has been performed on the gridded data cube.</figcaption></figure>

Cas A, and NVSS J011732+892848, a bright source at 30 arcmin from the target direction (see Fig. 6) that is also used as an amplitude calibrator for observed LOFAR EoR data. Both sources dominate the flux density in their respective DD-calibration source clusters.

### 4.5.1 Behaviour of gain phases

The left panels of Fig. 14 show the difference between the phase shift as a result of ionospheric dispersion and the phase correction performed during calibration. The values shown are statistics over time.<fnref n="22" /> In the case of a perfect correction for the ionospheric phase errors, this difference would vanish. The split in the mean values occurs due to baseline orientation, a sign change occurs between a baseline oriented parallel to the movement of the frozen TEC screen compared to a similar baseline oriented antiparallel to this movement. The red envelope describes the expected standard deviation of phase

errors for a source in the zenith at $r_{\text{diff}} = 10$ km, namely $\sqrt{D_{\phi}(\mathbf{b})}$ rad.

The top row demonstrates that overfitting is not a significant issue for target field sources, because the distribution of the phase errors does not evolve significantly during subsequent processing steps. Although the mean values of the phase differences shift, the standard deviations remain at approximately the same level as given by the structure function, despite it being impossible to solve for the rapid temporal variations of the ionosphere on the shortest temporal scales. This implies that neither a large bias nor an increase in gain errors are imprinted during data processing. The observed errors are the direct result of the ionosphere itself on time-scales that are too short to solve for on short baselines. For Cas A, on the other hand, the variance of the gain solutions drastically increases in the DD-calibration step on all baselines.

We conclude that this occurs when Cas A passes through nulls in the station beams. When this happens, there is very little signal power to constrain the gain. Due to this low signal-to-noise ratio, some stations obtain gains that deviate strongly from the true values. In these short intervals the overall standard deviation increases. The

<footnote n="22">²²Within the simulation, this means that the averages are taken with to a changing diffractive scale for Cas A, since its elevation changes over time.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Difference between the ionospheric distortions and baseline gain at 140.3 MHz at different stages of data processing as a function of baseline length. Top row: bright calibrator source NVSS J011732+892848 (also visible in the images of Fig. 6), bottom row: Cas A. From left to right the first three columns show the phase shifts of the TEC in the direction of the source (without including gains), the difference between the ionospheric phase shifts and the initial gain solutions, and the difference between the ionospheric phase shifts and the combination of both initial and DD-phase solutions. The right two columns show the amplitude levels of the gain solutions for initial calibration baseline gains, and the combination of initial and DD-gains, respectively. The black dots represent the median value across time, and their blue error bars show the standard deviation. These points split apart for longer baselines because of sign changes due to baseline orientation. The bottom-right panel is plotted in a different colour and on a different scale. The vertical black dashed lines show the 250 λ baseline length. Finally, the red lines show the expected behaviour for a source in the zenith at a 10 km diffractive scale.</figcaption></figure>

visibilities at these nulls do not dominate the residual power of Cas A, because in or near the null, the erroneous gains are strongly attenuated by the synthesized beam as the source traverses a null during visibility prediction. The overall impact of the off-axis sources on the power spectrum is therefore limited when they traverse a beam null (Cook et al. 2022).

### 4.5.2 Behaviour of gain amplitudes

The right-hand panels in Fig. 14 show the gain amplitude per baseline averaged over time. Because the ionosphere does not directly influence the amplitude of incident radiation in our simulations, these amplitudes approach unity in the ideal case. However, the ionosphere does indirectly influence the observed source amplitudes. This is due to uncorrected phase errors due to the ionosphere, which leads to ‘seeing’ determined by the phase structure function around these sources and a lowered peak flux compared to the original point source. The flux is spread out more and the speckle noise increases random errors in the data. The seeing-halo is given by $\langle V(\mathbf{b}) \rangle = V(\mathbf{b})\exp[-\frac{1}{2}D_{\phi}(\mathbf{b})]$ (Koopmans 2010; Vedantham & Koopmans 2015), where $\langle \cdot \rangle$ indicates the expectation value when averaged over sufficient time or well-separated baselines, and **b** a baseline vector. The shape of this amplitude drop-off is shown in red on the right side of Fig. 14.

To compensate for the lower peak flux compared to the sky model as a result of this ‘seeing’, the amplitude of the average gain is lowered during calibration. This effect is stronger for longer baselines, because they experience larger ionospheric errors. In the direction of Cas A, the variance of the amplitudes is very large, due to the aforementioned issues near nulls.

### 4.5.3 Gain phases and amplitudes in the power spectrum

To test the effects of both types of corrections separately, we compute cylindrically averaged power spectra using only phase and only amplitude corrections, respectively. We do not recompute the gains,

but instead apply either the phase or the amplitude part of the solutions. We do this both for the gain correction during initial calibration and for source subtraction during DD-calibration. Fig. 15 shows the results for simulation x10a. The top panel shows the ratio between a power spectrum created using only amplitude corrections and a power spectrum created using gains of unity. The middle panel shows the same, but for phase instead of amplitude and the bottom panel shows the full gains versus gains of unity.

In the top panel, the ratio only deviates from unity inside the wedge, where foregrounds and their mode-mixing effects reside. All deviations are larger than unity, indicating that, in the absence of instrumental errors, amplitude corrections are detrimental to foreground removal. In the modes strongly affected by A-team flagging (the wedge area of Cas A around $k_{\perp} = 0.15-0.19 \, h \, \text{Mpc}^{-1}$), we see increased errors due to amplitude solutions, but not due to phase solutions. In the plot which uses only phase solutions (middle panel), stronger residuals are introduced at low $k$-modes. However, a slight reduction in residual power appears in other places in the wedge. This coincides with the larger $k_{\perp}$ modes on the delay lines of the NCP source and Cas A, indicating that in some calibration directions, the phase corrections reduce the residual power. Overall, however, we still see an excess as a result of phase calibration. We conclude that we are only able to solve for the ionosphere in the phase part of the gain of some of the longest baselines in power estimation, but are unable to do so for most of the power spectrum. Furthermore, when gain amplitudes are also computed, we cannot solve for the errors on any power-spectrum baselines anymore.

## 4.6 Time coherence

The time coherence of residuals after calibration is one of the most important metrics to determine whether the LOFAR EoR data processing is successful. A low or absent time correlation indicates that integrating more data will lead to a reduction in excess variance, and therefore make a detection of the 21-cm-signal possible. A caveat to this is that an uncorrelated error on top of the thermal noise can
<figure><img src="image_17.png" /><figcaption>Figure 15. The ratio between the correction and sky-model subtraction using only part of the gain solutions and an uncalibrated residual computed using gains of unity for sky-model subtraction. We perform both initial calibration and DD-calibrated source subtraction. Top: only amplitude solutions are used, middle: only phase solutions are used, bottom: both phase and amplitude solutions are used (as in the regular pipeline). Amplitude solutions without phase solutions show an increased residual throughout the wedge. Phase solutions without amplitude show a slight reduction in power near the delay of Cas A and at the largest $k_{\perp}$ at the lowest $k_{\parallel}$, but an increase at the lower $k$-modes. A combination of both parts of the gains behaves like a sum of both excesses.</figcaption></figure>

still make the required time integration so large that it becomes prohibitive. If strong time correlations are present, however, any residual power from the foregrounds could pose a systematic limit on the detectable signal power. On the other hand, strong correlation may make it easier to model the residual foregrounds such that they can be removed more effectively.

To gauge the correlation between foreground residuals in different MSs, two metrics are used: the cross-power and the cross-coherence. The cross-power between simulations i and j is computed as

$$P_{ij}(k_{\perp}, k_{\parallel}) = \mathcal{V}_c \langle \tilde{I}_i^*(u, v, k_{\parallel}) \tilde{I}_j(u, v, k_{\parallel}) \rangle [\text{K}^2\text{h}^{-3}\text{cMpc}^3]. \quad (5)$$

Here, $\mathcal{V}_c$ is the observed comoving cosmological volume, $\langle \cdot \rangle$ indicates a cylindrical average in $(u, v)$-space to transfer to cylindrically averaged coordinates, and $\tilde{I}(u, v, k_{\parallel})$ refers to gridded brightness temperature cube at voxel $(u, v, k_{\parallel})$. We also use the cross-coherence between simulations i and j,

$$C_{ij}(k_{\perp}, k_{\parallel}) = \frac{\langle \tilde{I}_i^*(u, v, k_{\parallel}) \tilde{I}_j(u, v, k_{\parallel}) \rangle}{\sqrt{\langle |\tilde{I}_i(u, v, k_{\parallel})|^2 \rangle \langle |\tilde{I}_j(u, v, k_{\parallel})|^2 \rangle}}. \quad (6)$$

Since the cross-coherence is normalized by the power spectra, its value is between -1 and 1 and its magnitude is large if there is a strong (anti)correlation, regardless of the power level. We determine both the cross-power and cross-coherence to illustrate not only the level of correlation ($C_{ij}$) but also the impact that would have on a power-spectrum measurement ($P_{ij}$).

We have not included the residual foreground removal step in this analysis, because the residuals after GPR are already below the thermal noise floor. Additionally, we do not perform uv-space masking, to be able to interpret how the residual of Cas A behaves. Finally, thermal noise has not been included in this analysis. The reason for this is, first of all, that realistic thermal noise is uncorrelated over time, and will therefore artificially decrease the coherence and increase the sample variance. Furthermore, we have chosen to perform several simulations with the same thermal noise realization to avoid random realization differences between simulations. This means that the thermal noise would dominate the cross-power between these simulations if it is not removed.

The noiseless cross-spectra are shown in Fig. 16. In the first two columns, we show the correlation between two different LST ranges. The diffractive scales in the leftmost column are 10 km for both simulations, whereas the observation in the middle column has a scaled TEC screen with $r_{\text{diff}} = 5$ km (but otherwise the exact same ionospheric realization). In both simulations, the cross-coherence has a peak at the delay of Cas A around $k_{\perp} = 0.15-0.19$ hMpc$^{-1}$. The location of this peak is comparable between the two different ionospheric activity levels and coincides with the pixels where A-team flagging has a large effect (see the right panel of Fig. 13). The peak in the cross-power is slightly higher for the higher ionospheric activity level. For real data, this means that after combining nights, some correlation between different observations is expected to be present at the top of the wedge and the coherent power level is expected to depend on the residual power level in the individual observations. However, as A-team flagging and GPR are still performed after the steps we show here, this correlated power is likely removed later.

The third and fourth columns, which show the cross-coherence between two simulations at the same LST and date, show a large increase in cross-coherence at delay modes near and above Cas A. The third column has the same diffractive scales as used in the second column, whereas the fourth column uses the simulation with the same ionospheric realization at a higher activity level again. Interestingly, the correlation region is narrower in the fourth column, even though both simulations in this spectrum use the same diffractive scale. The opposite would be expected from ionospheric speckle noise, because the 'seeing-halo' around the sources would take the same form at the same LST and diffractive scale, whereas a higher diffractive scale would lead to a narrower halo and therefore less correlation. Therefore, we attribute the increased correlation for the same LST to a gain bias introduced during calibration. This is most likely a gain amplitude-based effect, because the gain amplitudes cause an increase in the foreground residuals in the relevant parts of the power spectra (see Fig. 15).
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Cylindrically averaged cross-power spectra (top) and cross-coherences (bottom). Left: correlation between simulations r10a and r10b (both with r<sub>diff</sub> = 10 km, but at a different LST). Second column: correlation between simulations r5a and r10b. The TEC screens are the same as on the left, but the TEC-values of r5a have been scaled to have r<sub>diff</sub> = 5 km. Third column: correlation between simulations r10a simulation and r5b. These simulations also have diffractive scales of r<sub>diff</sub> = 10 km and r<sub>diff</sub> = 5 km, like in the previous column, but these are at the same LST, such that the baselines observe the same sky. Final column: correlation between simulations r5a and r5b (both with r<sub>diff</sub> = 5 km and at the same LST). Cross-power in all plots is concentrated in the foreground wedge, approximately up to the upper delay line of Cas A. In the cross-coherence, the lower part of the wedge shows a weaker signal, indicating that the cross-coherence levels in this area are low, despite the cross-power being high. There is no clear difference between the two leftmost columns, whereas the rightmost columns show a larger coherence. This supports the idea that the coherence between different simulations is not dependent on the ionosphere itself, but more so on LST.</figcaption></figure>

We conclude that the coherence between different simulations is more dependent on LST separation than on ionospheric activity level. Simulations at different LST mostly have a high coherence at pixels that are strongly affected by A-team sources, which will be masked out later in the pipeline. As such, we expect ionospheric errors to be removed when combining a sufficient number of nights at different LSTs.

## 5 SUMMARY AND CONCLUSIONS

In this work, we analyse the impact of ionospheric effects on the 21-cm signal power spectrum extracted from simulated LOFAR HBA data sets that are processed in a very similar way to real observations. To do so, we create simulations that isolate ionospheric errors from other possible sources of error and process these simulated data sets using the LOFAR EoR data analysis pipeline. This data processing pipeline includes initial gain-calibration, multidirection gain calibration for source subtraction, uv-space masking of bright A-team sources, and residual foreground removal through GPR (Mertens et al. 2020). Crucially, these steps are performed using the same calibration settings as used on real data, including a 50–250λ and 250–5000λ baseline split in which longer baselines are used for calibration and shorter baselines for 21-cm signal estimation. Our simulations of a typical LOFAR HBA data set of the NCP at z = 9.1 include ionospheric phase errors representing one nominal, and one more extreme case, at the limit where real data are typically discarded. We conclude the following:

(i) Ionospheric errors that occur under normal observing conditions with LOFAR HBA ($r_{\text{diff}} > 5$ km) are not a dominant component of the excess power as seen in the 21-cm power spectrum. Unlike the results for a single night of real observations shown by Mertens et al. (2020), we are able to remove the foregrounds to below the thermal noise, despite having a lower diffractive scale (Gan et al. 2023).<sup>23</sup> After GPR, no excess variance remains at a level exceeding 10 per cent in any of our simulations, considerably lower than the excess variance of factors a few in current observations. We conclude that the propagation of ionospheric errors from longer baselines, which are the most affected by induced phase errors, to shorter baselines, is not a dominant contributor to the excess variance.

(ii) The ionosphere makes signal extraction more difficult. We have shown that configurations that remove the ionospheric errors exist. However, our configuration shows a positive noise bias of a few per cent, and the fitted component for mode-mixing varies on spectral scales that could approach those of the 21-cm signal itself. This can potentially cause signal suppression, so more extensive tests of the GPR settings, or more physically motivated foreground kernels, are needed in the case of real data.

(iii) The gain solutions that try to correct for ionospheric errors imprint low-level coherent errors on to the data after sky-model subtraction. This even occurs when data sets with uncorrelated ionospheric phase screens are used. This is mostly likely due to the bias in the gain solution as a result of the decorrelation of the visibilities within a solution interval.

<footnote n="23"><sup>23</sup>The pipeline used in this work is slightly different, however.</footnote>
(iv) Removal of bright off-axis sources remains more difficult than sources in the target field (the NCP). Errors remain present near the delay of Cas A in our simulations. For a smaller (5 km) diffractive scale, (u, v)-space masking is needed to remove the source down to the thermal noise level. Furthermore, the gains computed in the direction of this source show a high variance, possibly due to the crossing of nulls and side lobes of the primary beam. We show that these residuals also occur when the modelling and calibration beams are identical and therefore devoid of beam modelling errors. Although these subtraction errors are exacerbated by the ionosphere, we cannot completely attribute them to ionospheric errors, as they also occur in our simulation without an ionosphere.

Based on our findings, there is no need to exclude ionospheric activity levels currently used in processing. However, a physically motivated ionospheric covariance kernel may need to be added to our GPR step as higher sensitivity levels are approached, both for LOFAR and the future SKA. However, at the moment, the imprint of the ionosphere is well below that of the thermal noise after full data processing including GPR. Therefore, the ionosphere is not the cause of the observed excess variance in LOFAR HBA 21-cm signal data sets. It might be possible to improve the multidirection calibration step by taking speckle noise into account during forward prediction, such that gain amplitude suppression is avoided, especially in the case of off-axis sources such as Cas A. Because the sky contains several bright off-axis sources, the higher sensitivity of SKA may lead to many off-axis sources being detected above the noise level. It then becomes infeasible to mask all such sources, because this can lead to a significant loss of uv-coverage. These improvements may be especially useful for future processing with SKA, which has much higher sensitivity levels and may therefore suffer from ionospheric errors to a larger extent, although this may be mitigated somewhat by the better uv-coverage of SKA.

## ACKNOWLEDGEMENTS

SAB, LVEK, JKC, BKG, SG, CH, and SM 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'). EC acknowledges support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen. FGM acknowledges the support of the PSL Fellowship.

## DATA AVAILABILITY

The simulation and data processing pipelines underlying this article are available at https://doi.org/10.5281/zenodo.12188371. The simulations underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

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

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

* Albert J. G., Weeren R. J. v., Intema H. T., Röttgering H. J. A., 2020, A&A, 635, A147

* Asad K. M. B., Koopmans L. V. E., Jelić V., de Bruyn A. G., Pandey V. N., Gehlot B. K., 2018, MNRAS, 476, 3051

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

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

* Booker H. G., 1979, J. Atmos. Terr. Phys., 41, 501

* Bowman J. D. et al., 2013, Publ. Astron. Soc. Aust., 30, e031

* Buscher D. F., 2016, Opt. Express, 24, 23566

* Chege J. K., Jordan C. H., Lynch C., Line J. L. B., Trott C. M., 2021, Publ. Astron. Soc. Aust., 38, e028

* Chege J. K., Jordan C. H., Lynch C., Trott C. M., Line J. L. B., Pindor B., Yoshiura S., 2022, Publ. Astron. Soc. Aust., 39, e047

* Cook J. H., Seymour N., Sokolowski M., 2021, Publ. Astron. Soc. Aust., 38, e063

* Cook J. H., Trott C. M., Line J. L. B., 2022, MNRAS, 514, 790

* Datta A., Bowman J. D., Carilli C. L., 2010, ApJ, 724, 526

* Datta A., Bradley R., Burns J. O., Harker G., Komjathy A., Lazio T. J. W., 2016, ApJ, 831, 6

* DeBoer D. R. et al., 2017, PASP, 129, 045001

* Dewdney P. E., Hall P. J., Schilizzi R. T., Lazio T. J. L. W., 2009, Proc. IEEE, 97, 1482

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

* Edler H. W., de Gasperin F., Rafferty D., 2021, A&A, 652, A37

* Ewall-Wice A., Dillon J. S., Liu A., Hewitt J., 2017, MNRAS, 470, 1849

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

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

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

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

* Gasperin F. d., Mevius M., Rafferty D. A., Intema H. T., Fallows R. A., 2018, A&A, 615, A179

* Gasperin F. d. et al., 2019, A&A, 622, A5

* Gasperin F. d. et al., 2020, A&A, 642, A85

* Gehlot B. K. et al., 2018, MNRAS, 478, 1484

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

* Hallinan G., 2015, American Astronomical Society Meeting. American Astronomical Society Publishing, Washington, DC 20006 USA

* Intema H. T., Tol S. v. d., Cotton W. D., Cohen A. S., Gemmel I. M. v., Röttgering H. J. A., 2009, A&A, 501, 1185

* Jelić V., Zaroubi S., Labropoulos P., Bernardi G., de Bruyn A. G., Koopmans L. V. E., 2010, MNRAS, 409, 1647

* Jordan C. H. et al., 2017, MNRAS, 471, 3974

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

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

* Koopmans L. V. E., 2010, ApJ, 718, 963

* Koopmans L. V. E. et al., 2015, Proc. Sci., The Cosmic Dawn and Epoch of Reionization with the Square Kilometre Array. SISSA, Trieste, PoS#001

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

* Liu A., Shaw J. R., 2020, PASP, 132, 062001

* Loi S. T. et al., 2016, J. Geophys. Res., 121, 1569

* Lonsdale C. J., 2005, in Kassim N. E., Perez M. R., Junor W., Henning P. A., eds, ASP Conf. Ser. Vol. 345, From Clark Lake to the Long Wavelength Array: Bill Erickson's Radio Science. Astron. Soc. Pac., San Francisco, p. 399

* Martinot Z. E., Aguirre J. E., Kohn S. A., Washington I. Q., 2018, ApJ, 869, 79

* Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640

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

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

* Mevius M. et al., 2016, Radio Sci., 51, 927

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

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

* Morales M. F., Wyithe J. S. B., 2010, ARA&A, 48, 127

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

* Mouri Sardarabadi A., Koopmans L. V. E., 2019, MNRAS, 483, 5480

* Munshi S. et al., 2024a, preprint arXiv:2407.10686

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

* Narayan R., 1992, Phil. Trans. R. Soc. A, 341, 151

* Offringa A. R., Smirnov O., 2017, MNRAS, 471, 301

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

* Offringa A. R., Mertens F., Koopmans L. V. E., 2019a, MNRAS, 484, 2866

* Offringa A. R., Mertens F., Tol S. v. d., Veenboer B., Gehlot B. K., Koopmans L. V. E., Mevius M., 2019b, A&A, 631, A12

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

* Parsons A. R. et al., 2010, AJ, 139, 1468
Patil A. H. et al., 2016, MNRAS, 463, 4317

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

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

Rioja M. J., Dodson R., Franzen T. M. O., 2018, MNRAS, 478, 2337

Rufenach C. L., 1972, J. Geophys. Res., 77, 4761

Shaver P. A., Windhorst R. A., Madau P., de Bruyn A. G., 1999, A&A, 345, 380

Shen E., Anstey D., de Lera Acedo E., Fialkov A., 2022, MNRAS, 515, 4565

Sokolowski M. et al., 2015, ApJ, 813, 18

Swarup G., 1990, Indian J. Radio Space Phys., 19, 493

Tasse C. et al., 2021, A&A, 648, A1

Taylor G. I., 1938, Phil. Trans. R. Soc. A, 164, 476

van der Tol S., 2009, PhD thesis, Delft Univ. Technology, https://repository.tudelft.nl/islandora/object/uuid

van der Tol S., Jeffs B. D., van der Veen A.-J., 2007, IEEE Trans. Signal Process., 55, 4497

van Diepen G., Dijkema T. J., Offringa A., 2018, Astrophysics Source Code Library, record ascl:1804.003

Trott C. M. et al., 2018, ApJ, 867, 15

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

Vedantham H. K., Koopmans L. V. E., 2015, MNRAS, 453, 925

Vedantham H. K., Koopmans L. V. E., 2016, MNRAS, 458, 3099

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

Wandzura S. M., 1980, J. Opt. Soc. Am., 70, 745

Waszewski A., Morgan J., Jordan C. H., 2022, Publ. Astron. Soc. Aust., 39, e036

Weeren R. J. v. et al., 2016, ApJS, 223, 2

Yatawatta S., 2015, MNRAS, 449, 4506

Yatawatta S. et al., 2013, A&A, 550, A136

Zarka P., Girard J. N., Tagger M., Denis L., 2012, in Boissier S., de Laverny P., Nardetto N., Samadi R., D. Valls-Gabaud, Wozniak H., eds, SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. p. 687

Zaroubi S., 2013, in Wiklind T., Mobasher B., Bromm V., eds, The First Galaxies: Theoretical Predictions and Observational Clues. Springer, Berlin, Heidelberg, p. 45, https://doi.org/10.1007/978-3-642-32362-1_2

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: SKY MODEL

An excerpt from the supplementary sky-model file is tabulated in Table A1. Note that the right ascensions of sources within a calibration cluster can have a large spread, as this does not necessarily lead to a large angular separation near the NCP.
<table><caption>**Table A1.** Excerpt from the file that is provided as supplementary material to this work. All listed numbers are given with the same sign and number of significant digits as used in the file itself. Because all sources in the sky model are unpolarized, flat-spectrum point sources, we do not show the superfluous fields for polarized components, spectral shape, and characteristics of Gaussian sources.</caption><thead><tr><th>Type</th><th>Calibration cluster</th><th>Source label</th><th>RA</th><th>Dec</th><th>Flux density (Jy)</th></tr></thead><tbody><tr><td>cluster definition</td><td>3C61.1</td><td></td><td>+02<sup>h</sup>33<sup>m</sup>06<sup>s</sup>8705</td><td>+86<sup>°</sup>30<sup>′</sup>45<sup>″</sup>2960</td><td></td></tr><tr><td>cluster definition</td><td>Patch_0</td><td></td><td>+01<sup>h</sup>26<sup>m</sup>43<sup>s</sup>8591</td><td>+89<sup>°</sup>27<sup>′</sup>48<sup>″</sup>4997</td><td></td></tr><tr><td>cluster definition</td><td>Patch_1</td><td></td><td>-02<sup>h</sup>02<sup>m</sup>25<sup>s</sup>5563</td><td>+89<sup>°</sup>14<sup>′</sup>43<sup>″</sup>1056</td><td></td></tr><tr><td>...</td><td>...</td><td></td><td>...</td><td>...</td><td></td></tr><tr><td>point source</td><td>3C61.1</td><td>s0c125</td><td>+02<sup>h</sup>22<sup>m</sup>34<sup>s</sup>186</td><td>+86<sup>°</sup>17<sup>′</sup>19<sup>″</sup>596</td><td>37.493542737027255</td></tr><tr><td>point source</td><td>3C61.1</td><td>s0c199</td><td>+02<sup>h</sup>46<sup>m</sup>56<sup>s</sup>127</td><td>+86<sup>°</sup>18<sup>′</sup>30<sup>″</sup>985</td><td>0.2866027935546541</td></tr><tr><td>point source</td><td>3C61.1</td><td>s0c204</td><td>+02<sup>h</sup>17<sup>m</sup>34<sup>s</sup>851</td><td>+86<sup>°</sup>44<sup>′</sup>34<sup>″</sup>757</td><td>0.19606612284926137</td></tr><tr><td>point source</td><td>Patch_0</td><td>s0c1063</td><td>-02<sup>h</sup>35<sup>m</sup>35<sup>s</sup>333</td><td>+89<sup>°</sup>41<sup>′</sup>09<sup>″</sup>287</td><td>0.03310936307446051</td></tr><tr><td>point source</td><td>Patch_0</td><td>s0c1064</td><td>-02<sup>h</sup>36<sup>m</sup>52<sup>s</sup>226</td><td>+89<sup>°</sup>41<sup>′</sup>17<sup>″</sup>052</td><td>0.27724967346183066</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

## APPENDIX B: REFERENCE DATA PROCESSING PIPELINE

The data processing (or calibration) pipeline used in this work stems from the pipeline used to produce LOFAR EoR upper limits on the 21-cm signal, as discussed by Mertens et al. (2020) (see their fig. 1 in particular). We have chosen to base our work on a recently updated version of this pipeline, which is used to produce an improved version of these upper limits by Mertens et al. (in preparation). The advantage of using this improved pipeline is that our results reflect the ionospheric impact on the most competitive LOFAR EoR results. However, as the improved pipeline is not yet described in a publication, we list the updates with respect to the pipeline presented by Mertens et al. (2020) that are relevant to our work below:

(i) The time interval in initial calibration has been increased to 30 s rather than 10 s.

(ii) A-team flagging has been introduced after gridding. We describe how this is done in more detail in Section 3.2.3.

(iii) The Gaussian process models used for residual foregrounds removal have been updated following the introduction of ML-

GPR (Acharya et al. 2024; Mertens et al. 2024). As discussed in Section 3.2.5, the kernels used in the analysis of our simulations cannot be the same as those used on real data. We therefore use a set of kernels that resembles those investigated in the improvement of the data processing pipeline.

## APPENDIX C: ERROR EVOLUTION PER SIMULATION

Here, we provide propagation of ionospheric errors through the pipeline of simulations r10b (Fig. C1), r5a (Fig. C2), and r5b (Fig. C3). We cannot plot the ratio with the Null simulation and thermal noise realization without performing additional simulations for r10b, because this simulation uses a different LST and a thermal noise realization. As such, we only show the autopower spectrum for that simulation. Figs C2 and C3 show a larger residual of Cas A in the ratio between the simulation with an ionosphere and the Null simulation. Similar to what is discussed in Section 4.3, this effect is expected to be below the thermal noise level, and only visible here because the power spectra share the same LST and thermal noise realization.
<figure><img src="image_19.png" /><figcaption>**Figure C1.** Overview of the residual power in cylindrical power spectra of simulation r10b. The columns represent the different calibration steps in the same way as Fig. 7. The TRUE VISIBILITIES have been subtracted from the data in the first two columns, such that all columns represent residuals, even before the sky-model subtraction step in the pipeline.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure C2.** The same as Fig. 7, but for simulation r5a. To make the comparison between these figures easier, we use the same colour scale. From top to bottom, the rows show the power spectrum of the residuals, the ratio between this power spectrum and the input thermal noise power spectrum, and the ratio between the power spectrum and the same spectrum produced using the Null simulation. In the first three columns the input noise realization is used as thermal noise, whereas in the last column, this is replaced by a thermal noise cube drawn from GPR.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure C3.** Overview of the calibration of simulation r5b simulation. The figure is plotted in the same way as Fig. C2.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure D1.** Corner plot of the posterior probability distributions of GPR on the *r*<sub>diff</sub> = 5 km simulations. The vertical dashed lines in the marginalized distributions indicate the 68 per cent confidence level.</figcaption></figure>

APPENDIX D: POSTERIOR DISTRIBUTIONS
OF RESIDUAL FOREGROUND REMOVAL

In this appendix, we show the posterior distributions from GPR run
on the simulations with a 5 km diffractive scale. These simulations
have a more active ionosphere than those shown in Fig. 10, such
that the foreground variances are higher (less negative on the log-
scale).

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