© The Authors 2025

1. Introduction

If we look at the sky at 150 MHz, we expect to find more than one source per deg² with a flux density lower than 1 Jy (Wilman et al. 2008; Intema et al. 2017). Radio sources become increasingly rare as their brightness increases (e.g. Franzen et al. 2019), with only a few having flux densities higher than 100 Jy. The brightest radio sources in constellations are named after their constellation with an ‘A’ suffix, forming the so-called ‘A-team’. The most powerful A-team source at 150 MHz is currently Cygnus A (Cyg A; 3C 405, RA ${19}^{\mathrm{h}}{59}^{\mathrm{m}}28\stackrel{\text{s}}{.}36$, Dec $+40°44\prime 02\stackrel{″}{.}10$ in J2000; Baumgartner et al. 2013; Charlot et al. 2020; Vinyaikin 2014), an FR-II radio galaxy with a flux density of about 11 kJy at 150 MHz, that, together with Cassiopeia A (Cas A), Taurus A (Tau A), and Virgo A (Vir A), forms the northern A-team (e.g. de Gasperin et al. 2020).

The high brightness of the A-team sources is a double-edged sword for radio interferometric observations. On the one hand, these sources are so easy to measure even for low sensitivity instruments that their fluxes are well known, making them good (amplitude) calibrators, especially when unresolved. On the other hand, they can affect observations of a given target field even when they are tens of degrees away. This poses a challenge for wide-field instruments, such as LOFAR¹ (van Haarlem et al. 2013) and MWA² (Lonsdale et al. 2009), which have strong primary beam sidelobes. Different techniques have been developed to reduce such a contribution in the targetted observations. Examples are fast flagging of the data where the predicted contamination is above a certain threshold (Shimwell et al. 2017), more computational demanding direction-dependent (DD) calibration (Gan et al. 2023), demixing (Van der Tol 2009) when the A-team source is more than 30° away, or peeling (Noordam 2004; Mitchell et al. 2008) when it is closer. In all of these cases, good spatial and spectral models of the A-team sources are required to minimise the residual contamination after the removal.

Minimising the A-team contamination is particularly important for those experiments that aim to probe the neutral hydrogen (HI) during the Epoch of Reionisation (EoR) and Cosmic Dawn (CD) via the redshifted 21-cm line (see e.g. Furlanetto et al. 2006; Pritchard & Loeb 2012; Liu & Shaw 2020, for some reviews on 21-cm cosmology). Because the current generation of instruments does not have the sensitivity to make direct images of the 21-cm signal, they mainly focus on estimating the spatial fluctuations of the 21-cm signal by measuring its power spectrum. Increasingly more stringent upper limits have been set on the 21-cm power spectrum from the EoR by radio interferometers such as LOFAR (Mertens et al. 2020) and MWA (Trott et al. 2020; Kolopanis et al. 2023) with the highest frequency band, GMRT (Paciga et al. 2013), PAPER³ (Parsons et al. 2010; Kolopanis et al. 2019), and HERA⁴ (DeBoer et al. 2017; Abdurashidova et al. 2022), and from the CD by instruments such as LOFAR (Gehlot et al. 2019) and MWA (Yoshiura et al. 2021) with the lowest frequency band, OVRO-LWA⁵ (Greenhill & Bernardi 2012; Eastwood et al. 2019), NenuFAR⁶ (Zarka et al. 2012; Munshi et al. 2024), and AARTFAAC⁷ (Gehlot et al. 2020). The next generation of instruments, such as HERA and SKA-Low⁸ (Dewdney et al. 2009; Koopmans et al. 2015) are expected to have the sensitivity to detect the cosmological signal via power spectra, and SKA-Low has the potential to directly make tomographic images (e.g. Giri et al. 2018a,b; Bianco et al. 2024).

In all of the aforementioned cases, the extra-galactic and Galactic foreground emission must be mitigated, being up to five orders of magnitude brighter than the cosmological 21-cm signal (e.g. Jelić et al. 2008; Bowman et al. 2009; Bernardi et al. 2009, 2010). In principle, a separation between the two signals can be made based on the spectral smoothness. The foreground emission, being dominated by synchrotron radiation, is expected to be smooth over tens of MHz, while the 21-cm signal, being a hyperfine spectral line, fluctuates rapidly over MHz-scales (Shaver et al. 1999; Morales & Hewitt 2004; Santos et al. 2005; Morales & Wyithe 2010). However, instrumental systematics and incorrect calibration can cause spatial and spectral fluctuations in the observed data and mix the signals, making it harder to separate them (e.g. Datta et al. 2010; Vedantham et al. 2012; Barry et al. 2016; Ewall-Wice et al. 2017; Mazumder et al. 2022; Gorce et al. 2023). For this reason, having a precise model of the sky emission, with correct spectral information, is vital to minimise such errors.

Foreground avoidance methods, such as the ones used by the HERA collaboration, are less dependent on sky modelling, because they mitigate the foreground by using the foreground avoidance technique, which discards the Fourier modes dominated by the foregrounds (e.g. Parsons et al. 2012; Thyagarajan et al. 2015). However, they might still need accurate sky models to improve the calibration (e.g. Li et al. 2018; Kern et al. 2020) or to access the most sensitive small Fourier modes (e.g. Kerrigan et al. 2018; Ghosh et al. 2020). On the other hand, sky-modelling subtraction methods, such as the ones used by the LOFAR team, require an accurate sky model (Mertens et al. 2020; Gan et al. 2022). The latest upper limits published by the LOFAR-EoR Key Science Project (KSP) showed an excess variance that, on the larger angular scales, reaches 22 times the thermal noise in the residual power of the North Celestial Pole (NCP) field after foreground removal (Mertens et al. 2020). There might be multiple causes for this excess power, such as polarisation leakage (Jelić et al. 2015; Asad et al. 2015, 2016), model incompleteness (Barry et al. 2016; Patil et al. 2016; Ewall-Wice et al. 2017), low-level radio frequency interference (RFI; Wilensky et al. 2019; Offringa et al. 2019a), and uncorrected ionospheric effects (Vedantham & Koopmans 2016). However, Gan et al. (2022) concluded that the excess seen in LOFAR data is likely dominated by residuals of bright sources outside the main lobe of the primary beam. The sources that mainly affect the NCP field are Cas A and Cyg A, which are approximately 30° and 50° away from the field centre, respectively. At that distance, the chromatic effects of the primary beam are stronger and can introduce rapid spectral fluctuations to the sources, mixing their emission with other foregrounds and the 21-cm signal (Pober et al. 2016). If the beam sidelobes have large spatial and spectral gradients and the source is extended enough, the beam could also imprint different spectral structures to different parts of such a source, making it even harder to calibrate for it (Cook et al. 2022).

Previous efforts to model bright and extended galaxies have been made, for example, by Procopio et al. (2017) for MWA observations. They demonstrated that when such sources are incorrectly modelled as point sources, their subtraction leaves strong residuals. As a result, an improved wide-field model was made by Lynch et al. (2021) for use in the MWA processing pipeline. Efforts to improve models of bright sources were also made by Line et al. (2020) on the Fornax A radio galaxy. In this case, shapelets (Refregier 2003) were employed to model the extended emission. Shapelets allow a source to be modelled with fewer components than when using a collection of point sources and 2D Gaussian functions, thereby reducing the computational time required for calibration. However, incorporating spatially varying spectral information into shapelets remains challenging and has not yet been tested in real scenarios. Line et al. (2020) showed that, while a theoretical improvement in the 21-cm power spectrum was achieved with the shapelet model, the improvement was negligible in real observations due to other systematics. Similar results were observed by Gehlot et al. (2022) when modelling diffuse Galactic emission in AARTFAAC observations. Shapelets have also been applied to LOFAR observations by Yatawatta et al. (2013), who made a model of 3C 61.1, the brightest source in the NCP field. This model was subsequently used by Mertens et al. (2020). However, the fixed spectrum of the shapelet components proved to be a limitation, and Mertens et al. (2025) showed that using a multi-scale model of 3C 61.1 with a spatially varying spectral index led to improved upper limits. McKinley et al. (2022) showed that combining point sources and Gaussian components can still produce an accurate, high-resolution model of a large source such as Centaurus A. All these efforts suggest that having a spatially and spectrally realistic model of A-team sources is a fundamental requirement to further investigate the cause of the excess power.

In this paper, we will focus on Cyg A, which is the brightest A-team source at 150 MHz. It is a complex source, with both extended and compact emission. At low-frequencies, its spectral energy distribution deviates from a pure power law. Using several measurements between 12 MHz and 15 GHz, McKean et al. (2016; henceforth cited as MK16) showed that the total flux density S at frequency ν is well described by a third-order (n = 4) logarithmic polynomial function of the form

$$\begin{array}{c}\hfill S(\nu )=S({\nu }_0){\displaystyle \prod _{i=1}^{n-1}}{10}^{c_i{log}_{10}^i(\nu /{\nu }_0)},\end{array}$$(1)

where c₁ corresponds to the spectral index α and c₂ to the spectral curvature β (see MK16 for the coefficients values). This total spectral shape does not describe the spatial variation of the Cyg A emission, whose spectral indices go from α ≈ −0.65 in the lobes to α ≈ −1.25 in the plume region, with a rapid turnover in the two brightest hotspots, the western of which is labelled A and the eastern D (Leahy et al. 1989; MK16). In particular, the spectral index of hotspot A, measured at 150 MHz, is α_A = +0.18 ± 0.01, with the turnover around 140 MHz, while for hotspot D it is α_D = +0.36 ± 0.02, with the turnover at 160 MHz (MK16).

The main goal of this work is to create an accurate low-frequency spatial and spectral model of Cyg A and assess its impact on the LOFAR 21-cm power spectrum. Both Mertens et al. (2020) and Gan et al. (2022) used models of Cyg A with a spatially fixed spectral index α = −0.8 at a resolution that is not enough to capture the finest structures of the radio galaxy. By using a better spatial and spectral model of Cyg A, we can investigate whether an improved model of a bright off-axis source reduces the observed excess power in the 21-cm signal power spectra. In order to embed physical spectral information into the Cyg A model, we take advantage of the forced-spectrum method during a multi-frequency deconvolution Ceccotti et al. (2023). This method uses an input spectral index map to constrain the spectral index of each component by fitting a logarithmic polynomial function during the deconvolution. Without this regularisation of the spectral behaviour, logarithmic (or ordinary) polynomial fitting can generate unrealistic spectral indices, especially when the frequency bandwidth is limited or systematics in data are high (Offringa et al. 2016). In this work, we used an updated version of the method that supports higher order terms, such as the curvature, as long as maps of the higher order terms are also provided. The spectrally varying model of Cyg A that we get from the forced-spectrum method is the first of its kind in the analysis of the LOFAR power spectrum⁹. Another important goal of this work is to provide a strategy to make a model with physical spectral information by using the forced-spectrum method, thereby making it easier to make accurate source models. The final aim is given by combining the two goals: assessing the importance of improved models of off-axis sources, including spectral information, and the role they play in the excess power that is observed in the 21-cm power spectra.

In Sect. 2, we describe how we obtained the new Cyg A model, and compare it with the model used in the current LOFAR-EoR processing pipeline. In Sect. 3, we introduce both the simulated and observed NCP datasets that we use to test the impact of the new model onto the 21-cm upper limits. In Sect. 4, we describe the NCP data processing, from calibration to power spectrum estimation. The results on the two datasets are presented in Sect. 5. Finally, in Sect. 6, we summarise and discuss the results, drawing the conclusions.

2. Cygnus A modelling with LOFAR HBA

The current LOFAR-EoR processing pipeline uses the Cyg A model obtained from low-frequency VLA observations¹⁰ (see Appendix A for the catalogue of the model components), which we will refer to as the ‘old’ Cyg A model. This model has been used by standard LOFAR processing pipelines, such as LINC¹¹ (LOFAR Initial Calibration pipeline; de Gasperin et al. 2019), for calibration or demixing. Here we identify three major issues of this model when applied to HBA data:

1. Spectral extrapolation: the flux density of each component has been derived from a clean-based model of VLA observations at 73.8 MHz, and it is then extrapolated to higher frequencies using a given spectral index α;

2. Low resolution: the model consists of four point and eight Gaussian components, five of which have a zero-length minor axis, that are unable to model the complex structures of the source visible at higher frequencies (see Sect. 2.3.4);

3. Spatially constant spectral index: each component has an artificial, constant spectral index α = −0.8, even though the source has a strongly varying spectral index.

The total flux density given by the old Cyg A model at 73.8 MHz is 15.9 kJy, while the results of MK16 indicate it is 16.3 ± 2.2 kJy when extrapolated to that frequency using Eq. (1). Whereas this offset is acceptable within the 1σ error, issues (1) and (3) lead to a much larger difference at higher frequencies. Extrapolating the flux density of the old model to 150 MHz by using a first-order (n = 2) logarithmic polynomial with α = −0.8 results in a total flux of 9.0 kJy, which is 16% lower than the expected 10.7 kJy (MK16). The inconsistent flux density scale is another reason for improving the model of Cyg A. We will describe the new model in the next section.

2.1. The new Cygnus A model

LOFAR has the capabilities to operate between 110 and 250 MHz, by selecting two sub-ranges of the high band antenna (HBA) system, HBA-low in the 110–190 MHz range and HBA-high in the 210–250 MHz range. For each observation, a total of 488 sub-bands (SB) can be recorded, each with 195.3 kHz of bandwidth, which are further divided into 64 frequency channels. LOFAR currently consists of 24 core stations (CS) near the centre of the array, distributed within a diameter of ≈4 km and 14 remote stations (RS) across the Netherlands, with a maximum baseline of ≈120 km, forming what is called the Dutch array. Another 14 (international) stations are spread across Europe, providing baselines up to ≈2000 km. Each CS includes two HBA sub-stations (HBA0 and HBA1) that can either function independently (‘HBA Dual’ mode), improving the uv-coverage, or be individually selected (‘HBA Zero/One’ mode).

Given the large frequency coverage of the HBA system and the high spatial resolution provided by the Dutch array (approximately 3 arcsec), our goals with the improvement of the model for Cyg A are:

1. Obtain an accurate model at 120–160 MHz: this is the frequency range targeted by the LOFAR-EoR KSP and, therefore, we want to use HBA observations;

2. Increase resolution to 5 arcsec: hotspots are compact sources and we need a sufficient high resolution to not mix their emission with the surrounding lobes: 5 arcsec is the highest resolution provided by LOFAR HBA to get a model that is accurate for all Dutch baselines;

3. Capture the spatially varying spectral index: Cyg A is a morphologically complex radio galaxy, where different regions show very different spectral behaviours that have to be taken into account.

To satisfy goal (1), we use a combination of 12 h LOFAR HBA-low and HBA-high observations taken in 2014. The data, previously processed and calibrated (see Appendix B), consist of 412 SBs spanning 111–181 and 210–249 MHz, with a frequency resolution of 195.3 kHz (i.e. one channel per SB) and an integration time of 6 s. Other observational details are reported in Table 1. All modelling and data processing described in this paper has been performed on the ‘Dawn’ high-performance computing cluster (Pandey et al. 2020) at the University of Groningen.

In the next sections, we describe how we achieved goals (2) and (3) to get a new, more complete model of Cyg A. We used the WSCLEAN imager software (Offringa et al. 2014) to make images of the data. We start by describing the deconvolution methods that we employed.

2.2. MS-MF deconvolution and forced-spectrum fitting

Modelling the spectral structures of a radio source requires fitting some smooth function to the data from the frequency channels. This can be done during a multi-frequency (MF) clean approach, which models the spectral behaviour during deconvolution (Sault & Wieringa 1994). In the joined-channel deconvolution described by Offringa & Smirnov (2017) and implemented in WSCLEAN, a wide frequency band is divided into narrower output channels which are jointly imaged and deconvolved, and combined into a high dynamic range, frequency-integrated continuum image. The integrated image is used by the cleaning process to locate the peak of each model component. The frequency-varying brightness is measured across the images at each output frequency channel. This MF approach allows the components to capture the spectral variations of the source. Fitting the measured spectrum of each found component to an ordinary or logarithmic polynomial function ensures spectral smoothness of the model, which can be expected for sources dominated by synchrotron emission (e.g. Reich & Reich 1988; Tegmark et al. 2000; Jackson 2005; Jelić et al. 2008).

However, a clean-component spectral model might contain incorrect spectral information because, for instance, of the limited frequency bandwidth or when systematics in the data are large (e.g. Rau et al. 2016; Offringa et al. 2016). As a result, such spectral indices cannot be used for extrapolation. The forced-spectrum method (Ceccotti et al. 2023) solves this problem by using a spectral index map to assign to each clean-component a specific spectral index value based on the position of the component, fitting a logarithmic polynomial function of the form of Eq. (1). A modified weighted linear least-squares is then used to only determine the magnitude S(ν₀) of each clean-component. Ceccotti et al. (2023) used the forced-spectrum method to produce a model with only spectral indices. In this work, we additionally constrain higher order terms, such as the spectral curvature. This is necessary to model the spectral turnover of the Cyg A hotspots. When using higher order terms, it is important that the input spectral maps have been calculated at the same reference frequency at which WSCLEAN performs the MF spectral fitting. In Sect. 2.3 we will show a method to make spectral index and curvature maps from in-band data.

The forced-spectrum method, and in general the MF deconvolution, works well in combination with a multi-scale (MS) method (Cornwell 2008; Rich et al. 2008; Rau & Cornwell 2011). The MS deconvolution can incorporate diffuse emission by making a cleaned model as a summation of delta (i.e. point-like components) and basis functions (e.g. Gaussian or tapered quadratic). In this work, we used the MS deconvolution implemented in WSCLEAN (Offringa & Smirnov 2017) to model the extended emission with (circular) Gaussian functions that have an analytically defined Fourier transform, unlike tapered quadratic functions used by the algorithm described by Cornwell (2008). The combination of the MS and forced-spectrum methods result in a compactly described model, which is essential for calibration.

2.3. Imaging and spectral modelling of Cygnus A

To build the high-resolution spectral model of Cyg A, our imaging process was divided into four key steps: initial high-resolution imaging, imaging for accurate spectral modelling, spectral index extraction, and forced-spectrum fitting. This approach was designed to first capture the spatial details of the source, and then to accurately translate these observations into a comprehensive spectral model. The first step, the initial imaging, is optional and was done to check the data quality and the spectral features of hotspots A and D at the highest resolution allowed by the dataset.

2.3.1. High-resolution imaging and spectra of the hotspots

Our Cyg A dataset spans a wide bandwidth, from approximately 110 to 250 MHz, allowing us to create images with a synthesised beam (i.e. the point-spread function) of about 2 arcsec in size at the highest frequency. Although we cannot take advantage of this for reliable spectral index extraction, where maintaining a consistent beam size across the bandwidth is essential to ensure each frequency samples the same spatial scale, we can utilise it to generate high-resolution images for data quality assessment. However, if we clean deep enough so that the residual is mostly free of emission from the source, we can set the restoring beam to a fixed value and still get good spectral information from the restored images. This is the approach used in this section to extract the spectral behaviour of the hotspots.

For this initial step, we used the MF joined-channel deconvolution, splitting the full bandwidth into 100 independent (i.e. without spectral fitting) output channels. We did not use MF weighting, ensuring a better image quality for each output frequency by gridding the uv-plane of each channel on a separate grid. The output channels were weighted with a Briggs weighting scheme with a robust parameter of −2, resulting in an integrated synthesised beam with a full width at half-maximum (FWHM) of 3 × 2 arcsec. We used the MS deconvolution with point and Gaussian components of size 5, 10, 15, 20, 40, and 80 arcsec. We performed the cleaning to an initial threshold of 5σ using auto-masking and a final threshold of 1σ, where σ = 8.3 mJy/beam for the frequency integrated image, with a dynamic range of 1 : 19 540. We found that using a 1σ cleaning threshold, the residuals of the 100 output channels are noise-like at the hotpots location, but not completely on the rest of the source. This produced 10 times stronger residuals at the edges of the source at low frequencies and in the lobes at high frequencies. These are probably deconvolution artefacts caused by to the morphological complexity of the source. Because the condition of having the residual free of source emission is met for the brightest structures of Cyg A, we can still set a constant restoring beam to extract the spectra of hotspots A and D. We used a restoring circular beam with a FWHM of 3 arcsec, slightly larger than the frequency-integrated synthesised beam. Other imaging parameters are listed in the ‘High resolution’ column of Table 2.

The resulting 111–249-MHz frequency-integrated image is shown in the left panel of Fig. 1. At 3 arcsec resolution, the two brighter hotspots, A and D, are well resolved and we can extract their spectra from the 100 output images. We fitted their peak brightness using a non-linear least-squares method with a second-order (n = 3) logarithmic polynomial of the form

Fig. 1. High-resolution, frequency-integrated image of Cyg A between 111 and 249 MHz (left), and peak brightness of the hotspots A (top right) and D (bottom right), extracted from the high-resolution images of the 100 output channels (black solid line, with 1σ uncertainties shown as a grey shaded area). The fits using a second-order (blue dashed line) and a first-order logarithmic polynomial function (red dot-dashed line) are also shown, with the best fit parameters reported in the same colours. The latter has been fitted only between 111 and 181 MHz for a consistent comparison with MK16. The image on the left has a noise level of σ = 10 mJy/beam, with contours starting at 1 Jy/beam and increasing by a factor of 2.

$$\begin{array}{c}\hfill S(\nu )=S({\nu }_0){\left(\frac{\nu }{{\nu }_0}\right)}^{\alpha }{10}^{\beta {log}_{10}^2(\nu /{\nu }_0)},\end{array}$$(2)

which is equivalent to Eq. (1) when using three terms, where α and β are, respectively, the spectral index and the curvature terms evaluated at the reference frequency ν₀. We used 150 MHz as reference frequency to compare our spectral values with MK16. The panels on the right-hand side of Fig. 1 show these fits with a blue dashed line. Here, the black solid line represents the peak brightness as extracted from each output image, to which we have associated 1σ errors (grey shaded area). These are calculated by adding in quadrature the standard deviation of the noise of each image and a standard calibration error of 2% of the flux density (MK16), where the latter dominates the uncertainties. We find that the best fit values for the hotspot A are α_A = +0.1025 ± 0.0002 and β_A = −1.541 ± 0.009, and for the hotspots D are α_D = +0.3081 ± 0.0002 and β_D = −1.819 ± 0.009, confirming the presence of a turnover in both of them. Our values differ from those found by MK16 for a couple of reasons: they only had HBA-low data in the range 109–183 MHz, and they fitted a first-order logarithmic polynomial. Restricting our analysis in the same frequency range and fitting a simple power law without the curvature term (red dot-dashed line in Fig. 1), we find α_A = +0.1720 ± 0.0004 and α_D = +0.3917 ± 0.0004, which are consistent with the values reported by MK16, namely, α_A = +0.18 ± 0.01 and α_D = +0.36 ± 0.02. The 1σ uncertainties on α and β assume uncorrelated errors, whereas systematics such as calibration and deconvolution errors are typically correlated, and are therefore not captured by these uncertainties.

2.3.2. Spectral index extraction

The second step towards constructing an accurate spectral model of Cyg A was to generate images that could be used for a pixel-by-pixel spectral index extraction. As mentioned in the previous section, using the 100 high-resolution images was not ideal due to the ≈10σ residual emission, which could affect the spectral indices of fainter emission. While the clean model would be convolved with a consistent restoring beam across the bandwidth, the residual images retain the synthesised beam of their respective output channels. This might result in a mix of components at different spatial scales from low to high frequencies, preventing the extraction of a reliable source spectrum. To address this, we performed a separate imaging step with appropriate settings for an accurate full-source spectral index extraction.

Because we could not be certain that the residual images were noise-like across the entire bandwidth, we ensured a consistent synthesised beam (i.e. the same uv-coverage) by applying a strict uv-cut. We selected only the (projected) baselines in the range 160 − 45 500 λ, where the lower value corresponds to the shortest baseline at 249 MHz and the upper value to the longest baseline at 111 MHz. By removing baselines shorter than 160 λ, we lost sensitivity to structures larger than ≈20 arcmin, which was not a concern since our focus was only on Cyg A, which has a maximum extent of 2.3 arcmin. Conversely, the upper uv-limit reduced the resolution from 3 to 4.5 arcsec, leading us to fix the restoring beam to a circular Gaussian with a FWHM of 5 arcsec. Consequently, we increased the pixel scale from 0.25 arcsec (used for high-resolution imaging) to 1 arcsec, producing images of 256 × 256 pixels. We also reduced the number of output channels to 41 to get a better noise per output image. This number was still high enough to minimise any bias that would occur during the later fitting process. We measured a noise level of σ = 12.8 mJy/beam for the frequency-integrated image (dynamic range of 1 : 20 550). With these settings, we found that the MS deconvolution gave better results with the default spatial scales of 12, 25, 49, and 99 arcsec. The other WSCLEAN parameters were similar to those used in the high-resolution imaging. As we will apply similar settings in the next forced-spectrum fitting step, all parameters are listed in the ‘Forced spectrum’ column of Table 2, with specific parameters for this step provided in square brackets.

The resulting 41 output images were used to extract spectral index and curvature maps. The brightness of each pixel was fitted using Eq. (2) at the reference frequency of 190.18 MHz, which is the same frequency at which WSCLEAN extracted the model clean-components for our specific settings¹². In this operation, we considered only pixels brighter than 10 Jy/beam in the frequency-integrated image (this roughly corresponds to 3.6 Jy/beam in the frequency-integrated, high-resolution image of Fig. 1). The extracted spectral index and curvature maps are shown in Fig. 2. Although the image resolution is lower than in the preliminary imaging step, the large LOFAR HBA bandwidth and high image dynamic range still allow us to identify the main spectral features across the entirety of Cyg A. Hotspots A and D exhibit the flattest spectral indices and the lowest curvature, as expected due to their turnover. Additionally, in line with synchrotron spectral ageing models (e.g. Carilli et al. 1991; MK16), the spectral index map shows a steepening from the hotspots, where relativistic particles are continuously accelerated by the jets, to the outer regions of the lobe and counter-lobe, where the primary energy loss mechanism is synchrotron radiation. More energetic plasma, which emits at higher frequencies, loses energy faster than less energetic plasma, creating a break frequency that shifts to lower frequencies as the plasma ages (unless re-acceleration occurs). Because the plume of Cyg A has the steepest spectral index (α < −1), it consists of the oldest plasma in the source, indicating that this plasma underwent its last acceleration much earlier than in any other region.

Fig. 2. Forced-spectrum input spectral index (left) and curvature maps (right) at 190.18 MHz, extracted from the 5 arcsec resolution images. Contours in both figures are the same of the left panel of Fig. 1. The colour scale of the spectral index map spans the range −1.27 ≤ α ≤ −0.38, while the curvature is plotted in the range −1.04 ≤ β ≤ 1.04.

2.3.3. Forced-spectrum fitting and spectral rescaling

In both Sects. 2.3.1 and 2.3.2, the imaging was performed without any constraints on the spectral behaviour. In the next processing step, the extracted spectral index and curvature maps were used as inputs to the forced-spectrum method, acting as an imaging regularisation. We mostly applied the same WSCLEAN settings as in the previous step, as listed in the ‘Forced spectrum’ column of Table 2. The same uv-range was used to ensure that the sampled spatial scales were consistent with the spectral index map. While we could have used higher resolution to obtain a more detailed spatial model of Cyg A, the 5 arcsec resolution of the spectral index map meant that this would not have provided any additional benefit on the spectral modelling side. Furthermore, we preferred to maintain the same uv-cut and resolution to limit the number of clean-components in the final model, which consisted of 2082 components (1320 point sources and 762 Gaussians).

Due to the spectral constraints of the forced-spectrum method, the residual images may exhibit higher residuals, resulting from the clean-component overlap described in Ceccotti et al. (2023). This occurs because the spectral index is assigned based on the clean-component position, but for Gaussians, the wings overlap with other components in regions where the spectral index differs, causing a bias in the final model image and, consequently, in the residual image. For spectrally complex sources such as Cyg A, this bias can be significant, particularly in the extended regions where many Gaussians are used during cleaning. To mitigate this, we constrained the cleaning process using a spatial mask that included only pixels brighter than 1 Jy/beam, based on the previous frequency-integrated image (corresponding to approximately 0.4 Jy/beam in the high-resolution image of Fig. 1). This approach allowed us to push the initial cleaning threshold deeper, reducing residual emission and preventing artefacts from being cleaned beyond the source extension. As a result, clean components can still be selected between the spatial mask boundary and the spectral index map extension, which was defined down to 10 Jy/beam. To minimise the impact of the overlap bias in the spectrum of the final model within the 10 Jy/beam contour, we assigned values of α = 0 and β = 0 to the remaining areas of the spectral maps.

The output of our forced-spectrum imaging step is a new, high-resolution model of Cyg A, whose total flux density matches, between 110 and 250 MHz, the spectrum found by MK16 as shown in the top panel of Fig. 3. However, when extrapolated outside the LOFAR HBA frequency range, the spectral energy distribution diverges from MK16, likely because a sum of components with varying spectra would not necessarily follow the flux scale to which the data were calibrated. Moreover, the overlap bias might play a role. While extrapolation always requires caution, this new Cyg A model could be useful for applications at lower frequencies, such as calibrating or demixing LOFAR low band antenna (LBA) data, between 10 and 80 MHz. To address this, we evaluated and corrected the total flux density of our model at 10 frequencies between 10 and 300 MHz to match MK16. Each clean component was then adjusted by the found correction factors and fitted to a logarithmic polynomial to restore spectral index and curvature values at the model reference frequency. The rescaled model agrees with MK16 as shown by the black line in Fig. 3. The middle and bottom panels show the flux density of hotspots A and D, integrated within a circle of 5 arcsec diameter. Compared to Fig. 3, the hotspot spectra appear flatter, with the turnover shifted to lower frequencies. This is due to the lower resolution used for both the spectral maps and the final model, allowing more lobe emission to contaminate the hotspots. Such flattening is accentuated by the rescaling, although most effects impact frequencies outside the LOFAR HBA range.

Fig. 3. Flux density between 50 and 270 MHz measured from the Cyg A models for the total emission (top), hotspot A (middle), and hotspot D (bottom). The top panel shows the old model (solid blue line), the new model before (dashed green line) and after spectral rescaling (solid black line), and the third-order logarithmic polynomial used by MK16 (dashed orange line). In the middle and bottom panels, only the flux densities from the new Cyg A models before and after rescaling are shown. These values were extracted from a circular region with a 5 arcsec diameter, centred on the peak of each hotspot. The filled points indicate the spectral peaks, and the turnover frequency after rescaling is also reported. The green shaded area in all panels marks the LOFAR HBA range used to extract the new Cyg A model.

2.3.4. Old versus new Cygnus A model

The new, rescaled model of Cyg A and the old model are shown, respectively, in the top and bottom panels of Fig. 4, both rendered at 150 MHz. It is immediately clear that the new model is at a higher resolution than the old one. Hotspots A, B and D are modelled with a dominant point component, while the diffuse lobe and counter-lobe emission consists of a collection of Gaussian components of different sizes. Because WSCLEAN uses only circular Gaussians during the MS deconvolution, point components are required to shape these Gaussians into real source structures, such as the bow-shaped structure in the counter-lobe. This explains the presence of components with negative fluxes, that have no physical meaning but they are still needed to get models resembling real data. The final result is a model where also the fine structures of Cyg A are visible and easily distinguishable, fulfilling goal (2), mentioned in Sect. 2.1. On the other hand, we can hardly distinguish the different structures in the old model. Assuming that the two outermost Gaussian components model the hotspots A and D emission, a shift in the coordinates towards the north-west is present in the old model, the cause of which is unknown. This issue would be easily solved by shifting back each components to the correct location in the sky, but it is hard to find a reference point, such as hotspot A, in this model. When the old model is used on its own for calibration or demixing, such an offset, which is not too large, is absorbed into the gain phases.

Fig. 4. New (top) and old model (bottom) of Cyg A, rendered at 150 MHz. Contours are the same of Fig. 1. The zero-length axes of the old model Gaussians have been replaced with 1 arcsec axes for imaging purposes, resulting in the appearance of vertical spikes.

The new model also outperforms the old one in terms of spectral modelling. The forced-spectrum method transferred the pixel-by-pixel spectral information from the spectral index and curvature maps into the new model, whose components have different spectra depending on their location, as we saw in the middle and bottom panels of Fig. 3. Thus, the resulting model of Cyg A also satisfies goal (3). Plotting the old model in the top panel, extrapolated from 74.8 MHz to the HBA frequency range, shows the large difference in flux scale between the expected values, as already discussed at the beginning of Sect. 2.

3. NCP datasets

To test the performances of the new model of Cyg A against the old one, we used it within the LOFAR-EoR processing pipeline, comparing the results both from a simulated dataset and from a real observation of the NCP field. This allows us to assess the impact of this better model on the 21-cm power spectrum measurements, which is one of the main motivations of this work. We selected an average-quality night from Mertens et al. (2020) where Cyg A is above an elevation of ≈10° for most of the observing duration. The observational details of the selected dataset are reported in Table 3.

To produce the simulated dataset, we made a copy of the same measurement set and forward predicted (i.e. Fourier transformed from image space to visibility space) a sky model into it. When predicting visibilities, the calculation of the primary beam is computationally expensive. The computing time increases with the number of model components, because the primary beam has to be calculated at each component coordinate. Because the NCP sky model has more than 28 700 components (Yatawatta et al. 2013; Mertens et al. 2020), we used a simpler NCP model with only 684 point components with a flat spectrum, presented in Brackenhoff et al. (2024). This model extends up to 6° from the NCP and consists of 14 source clusters (i.e. directions), and its sources dominate the total power spectrum of the sky, making it sufficient for the purpose of this work. More details are reported in Appendix C. We added a low resolution model of Cas A, consisting of 24 Gaussian components with α = −0.8 constant across the source and ν₀ = 74 MHz, similar to the old Cyg A model. Finally, we add our new Cyg A model. Assuming that this new model is the true sky, we will be able to measure the theoretical excess power that might be caused by not having a Cyg A model with correct spectral indices.

The combined sky model was predicted on to the visibilities of the previously selected night with SAGECAL¹³ (Yatawatta 2016), which allows using GPUs to speed up the process (see Sect. 4). The primary beam was applied during the prediction, and we finally added a per-visibility noise estimated from the radiometer equation (e.g. Thompson et al. 2017):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{vis}}=\frac{\mathrm{SEFD}}{\sqrt{2\mathrm{\Delta }t\mathrm{\Delta }\nu }},\end{array}$$(3)

where SEFD is the system equivalent flux density, Δt = 2 s is the visibility integration time and Δν = 61 kHz is the frequency channel width. The SEFD of CS and RS is expected to be different because RS have double collecting area (van Haarlem et al. 2013). However, our observations were carried with the ‘Inner’ setup, which sets the RS collecting area to the CS one by switching off the outer tiles of the station. For this reason, we assigned to every visibility the same SEFD of 4253 Jy as estimated by Mertens et al. (2020). To reduce the computing time, we kept only the first 6 h of the dataset, where the elevation of Cyg A peaks (77.9° after 35 min from the beginning) and stays above 35°.

The full-sky image of the NCP simulations is shown in Fig. 5, with a zoom-in of the NCP main field (top-right panel) and into the Cyg A field (bottom-right panel). Cyg A is 49.3° away from the NCP, and, because of that, the primary beam attenuates its flux by a factor of ≈10⁻⁴. However, its sidelobes are still strong and extend across the whole image, reaching the NCP main field and affecting its residuals if the source is not well modelled and subtracted. The same holds for Cas A, the effect of which could be even larger because it is closer to the NCP (31.2° away) and has a higher apparent brightness.

Fig. 5. All sky dirty image of the simulated NCP dataset (left), with zoom-in to the NCP main field (top right) and the Cyg A direction (bottom right). The all sky image has been obtained using a baseline range of 50 − 250λ, the NCP zoom-in with 50 − 500λ, while the Cyg A zoom-in with 50 − 4500λ. A natural weighting scheme has been used for the all sky and the NCP zoom-in images, whereas a uniform weighting has been chosen for the Cyg A zoom-in to make the source structures more visible. The contour in the bottom right panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

4. NCP data processing

Both datasets have been processed through an updated LOFAR-EoR production pipeline, which steps are extensively described by Mertens et al. (2020, 2025). The pipeline consists essentially of (1) pre-processing and RFI removal, (2) initial calibration, (3) DD solving and DD sky-model subtraction, (4) imaging, (5) removal of residual foregrounds, and (6) power spectrum estimation. In this work, we skipped step (1) because we did not include any RFI into the simulated measurement set and we start with the already pre-processed observed data from Mertens et al. (2020). Step (5), which is performed with the Gaussian process regression (GPR) method (Mertens et al. 2018, 2024; Hothi et al. 2021; Acharya et al. 2024) was also skipped for both simulated and observed datasets, because this step does not use a sky model and we therefore do not need it to make the final comparison. In the following, we briefly describe the steps used for this work.

4.1. Initial calibration

The main goal of the initial calibration is to set the correct flux scale and correcting the sky-averaged phase offset per station. For the NCP, this is done in an unconventional way, by obtaining calibration solutions for two directions. One direction covers 3C 61.1, which is the brightest source in the main field with a flux density of ≈38 Jy at 150 MHz (Baldwin et al. 1985) and close to the first null of the primary beam, making it necessary to solve the varying DD effects of the primary beam and imperfections in the LOFAR primary beam model. The other direction includes the rest of the field, where only sources brighter than 35 mJy are selected to reduce the computing time while keeping a high enough signal-to-noise ratio. This reduces the NCP sky model used for our real observations from 28 773 to 1416 sources, thereby reducing the computational cost. The remaining sources account for only 1% of the total flux model, and this has therefore only a small effect. For consistency, we performed a similar procedure for the sky model used for the simulations, but such a flux cut has a relatively small effect here, because the model has fewer sources (i.e. from 684 to 666 components). During the initial calibration, we only use baselines in the range 50 − 5000λ. The lower limit was applied in real observations to avoid poor uv-coverage and the diffuse Galactic emission which is not present in our model and would bias the calibration (Patil et al. 2017), while the upper limit minimises ionospheric effects and avoids resolving compact sources (Patil et al. 2016; Mevius et al. 2022).

The initial calibration was performed with SAGECAL, that uses the space alternating generalised expectation maximisation (SAGE) algorithm (Kazemi et al. 2011). It distributes the processing across multiple nodes using a consensus optimisation (Boyd et al. 2011). This method applies a frequency smoothness regularisation with a third-order Bernstein polynomial (Farouki & Rajan 1988; Yatawatta 2019), using an alternating direction method of multipliers (ADMM; Yatawatta 2015, 2016). The level of frequency smoothness for each iteration is determined by the regularisation parameter ρ, whose value depends on the apparent sky brightness for each source cluster and is therefore direction dependent (Mevius et al. 2022). More details about SAGECAL and its capabilities can be found in Yatawatta (2015, 2016, 2018), Yatawatta et al. (2017), and Gan et al. (2023). This smooth gain calibration was performed with a frequency and time solution interval of 195 kHz (i.e. one solution per SB) and 30 s, respectively, to address relatively fast direction-independent ionospheric effects on the longer baselines. We note, though, that we did not include any ionospheric effects in our simulations. Recent studies have shown that the effect of the ionosphere is currently still expected to be below the thermal noise level (Vedantham & Koopmans 2015, 2016; Brackenhoff et al. 2024). The data were corrected using the solutions from the main field.

After this high-temporal and smooth-spectral resolution calibration, we performed a low-temporal and high-spectral resolution bandpass calibration with the same sky model to solve for the fine frequency response of the instrument. This was done using a single solution per SB over a time interval of 3 h to mitigate the introduction of frequency structures due to an incomplete or inaccurate model as much as possible. We still used SAGECAL, but each SB was solved independently and no smoothness was applied to the solutions. As before, only the main field solutions were applied to the data.

After the bandpass calibration, all visibilities with amplitude higher than 70 Jy were flagged to prevent SAGECAL from diverging from the correct solutions in the subsequent DD calibration step (Mertens et al. 2020). We also averaged the data to a time resolution of 10 s to speed up the next step, which is the most demanding step of the pipeline.

4.2. Direction-dependent solving and source subtraction

The wide field of view of LOFAR introduces DD effects due to, for example, the time-varying primary beam and the ionospheric effects, that have to be taken into account to perform a good foreground subtraction. This is done by clustering the sky model components into a number of directions and ensuring a spectral smoothness using the consensus optimisation of SAGECAL. These clusters typically are 1 − 2° in size to reach a sufficient signal-to-noise with the sources inside the cluster and to minimise spatially varying DD effects. For the simulations, the NCP main field was divided into 14 clusters (see Table C.1), while for real observations it currently consists of 108 clusters, because the sky model extends up to 15° from the NCP. Two extra clusters were added for Cas A and Cyg A, because, similarly to the 3C 61.1 case, they are so bright and so far away that they need separate solutions. No flux limit was applied to the sky model during this DD step. The gains were solved for all clusters simultaneously using a frequency interval of 195 kHz. We set solution time intervals between 2.5 and 20 min for the main field directions and an interval of 2.5 min for Cas A and Cyg A to capture faster primary beam changes. All the sources in the sky models are included in this step.

To avoid over-fitting, a uv-cut at 250λ was introduced between baselines used for calibration and 21-cm signal extraction, where only baselines in the range 250 − 5000λ were used during the DD solving. This cut causes excess noise on the uncalibrated baselines, which is reduced by enforcing spectrally smooth solutions (Yatawatta 2016; Barry et al. 2016). Using this approach, it has been shown that there is no signal loss on the baselines of interest and very limited excess power (Mouri Sardarabadi & Koopmans 2019; Mevius et al. 2022).

In contrast to the initial calibration, after the DD gain solving we did not correct the data with the solutions. Instead, we used the gains when subtracting the sky model from the visibilities, to obtain the residuals that were finally used for creating power spectra.

In order to compare the performance of the new Cyg A model with the old one, we performed the DD calibration steps for both Cyg A models. In addition, we did a third run in which we modelled Cyg A as a point source with a spectral energy distribution following a third-order logarithmic polynomial with the coefficients given in MK16. This extreme case will demonstrate whether having a spatially resolved model is important, even when the source is not resolved on the baselines used for calibration. In fact, calibration can still be affected by the source structure if the signal-to-noise ratio is high enough, as in the case of Cyg A. The results from the simulations (see Sect. 5.2) that used the new model during the DD steps are only a benchmark for the calibration efficiency within the LOFAR-EoR pipeline, and do not say anything about how good our new model is. After the DD calibration, we considered only the first 6 h of the observed data, in order to better compare the results with the simulated dataset.

4.3. Imaging and conversion to Kelvin

After calibration and foreground subtraction, we gridded and imaged each SB of the residual visibilities using WSCLEAN, making an (l, m, ν) image cube. During the gridding, each visibility received equal weight and we applied a Kaiser-Bessel anti-aliasing filter with a kernel extent of 15 uv-cells, an oversampling factor of 4096, and 32 w-layers. These parameters ensure that any error arising from the gridding process is confined significantly below the expected 21-cm signal (Offringa et al. 2019b,a; Mertens et al. 2020).

We produced Stokes I and V images with a natural weighting scheme. Additionally, we generated even and odd time-step Stokes V images, the difference of which is used to estimate the thermal noise variance. We used a baseline range of 0 − 600λ during the imaging, setting a pixel scale of 0.5 arcmin and an image size of 512 × 512 pixels, which results in a field of view of 4.3 × 4.3°. After gridding, only the gridded uv-values between 50 and 250λ were used for the power spectrum.

The image cube, which has units of Jy/beam, has to be converted into units of brightness temperature, which is Kelvin. This is done by Fourier transforming the images into a gridded (u, v, ν) visibility cube and then using the conversion technique detailed by Offringa et al. (2019a). Making the visibility cube, we applied a spatial Tukey function with a 4° diameter to focus on the centre of the primary beam, which has a FWHM of ≈4.1° at 140 MHz, while also mitigating its uncertainties farther from the centre. The gridded visibilities (in units of Kelvin) and the number of visibilities in each (u, v, ν)-cell were used to estimate the power spectrum.

4.4. Power spectrum estimation

If T(x) is the brightness temperature of a signal at the physical coordinate x, its power spectrum as a function of the wavenumber k (units of h cMpc⁻¹) is defined as

$$\begin{array}{c}\hfill P(\mathit{k})=V|\stackrel{\sim }{T}{(\mathit{k})|}^2,\end{array}$$(4)

where V is the observed comoving volume and $\stackrel{\sim }{T}$ is the discrete Fourier transform of the brightness temperature. This 3D power spectrum is usually in units of K² h⁻³ cMpc³. The components of the wavenumber k perpendicular and along the line of sight are respectively defined as (Morales & Hewitt 2004; Vedantham et al. 2012):

$$\begin{array}{c}\hfill {\mathit{k}}_{\perp }=\frac{2\pi }{D_M(z)}\mathit{u},k_{\Vert }=\frac{2\pi H_0{\nu }_{21}E(z)}{c{(1+z)}^2}\eta ,\end{array}$$(5)

where D_M(z) is the transverse comoving distance at redshift z, ν₂₁ = 1420 MHz, H₀ is the Hubble constant, E(z) is the dimensionless Hubble parameter, c is the speed of light, and η is the Fourier dual of the frequency ν (often referred as delay).

The power P(k) can be averaged in cylindrical shells of radius k_⊥ = |k_⊥| to get the cylindrically averaged (2D) power spectrum:

$$\begin{array}{c}\hfill P(k_{\perp },k_{\Vert })=\frac{{\displaystyle {\sum }_{\mathit{k}\in (k_{\perp },k_{\Vert })}}P(\mathit{k})}{N_{(k_{\perp },k_{\Vert })}},\end{array}$$(6)

where N_{(k⊥, k∥)} is the number of k voxels falling into the (k_⊥, k_∥)-bin. If we choose spherical shells of radius $k=|\mathit{k}|=\sqrt{k_{\perp }^2+k_{\Vert }^2}$, a spherically averaged (1D) power spectrum can also be obtained:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^2(k)=\frac{k^3}{2{\pi }^2}\frac{{\displaystyle {\sum }_{\mathit{k}\in k}}P(\mathit{k})}{N_k},\end{array}$$(7)

where N_k is the number of k voxels falling into the k-bin. This is the power spectrum estimator used to set 21-cm upper limits because it is in units of mK². Both cylindrical and spherical power spectra are weighted by the gridded visibility thermal noise noise (see Mertens et al. 2020, for more details).

Because k_∥ is essentially the Fourier dual of the frequency, a spectrally smooth signal is confined to low k_∥ values. This is the case for Galactic and extra-galactic foreground emission, which is expected to be dominated by synchrotron radiation at low frequencies and spectrally smooth over tens of MHz. On the other hand, the 21-cm signal is fluctuating rapidly in frequency, and its power extends over all k_∥ (e.g. Santos et al. 2005). However, every interferometer is intrinsically chromatic because a single uv-cell is crossed by different baselines at different frequencies, and this effect becomes stronger at longer baselines, corresponding to larger k_⊥ (Morales et al. 2012, 2019). Because of this, the foreground power is spread at higher k_∥, generating a characteristic shape known as the ‘foreground wedge’ in the cylindrically averaged power spectrum (Datta et al. 2010; Liu et al. 2014a,b). The maximum extent of the wedge in k_∥ is given by the horizon line, that can be accurately calculated in a non-flat sky formalism knowing the horizon delay (Munshi et al. 2025). In the same way, we could find the delay lines for every direction in the sky (Munshi et al. 2025). We then expect that all the sky emission lies in the wedge below the horizon line, leaving a foreground-free region called ‘EoR window’. However, the finite bandwidth of the data causes the Fourier transformation to leak some power above the wedge, requiring some frequency window function to reduce such a leakage (e.g. Vedantham et al. 2012). Also calibration errors might introduce rapid frequency fluctuations that contaminate the EoR window. This effect can be limited by using frequency smoothed gains, as we did in the first calibration step and during the DD solving.

In this work, we estimated the power spectra using the power spectrum pipeline PSPIPE¹⁴. We used only gridded visibilities in the 50 − 250λ baseline range and applied a Hann window function to the frequency bandwidth. The reported power spectrum uncertainties were estimated from the sample variance, as described by Mertens et al. (2020).

5. Results

In this section, we show the results from both the simulated and observed datasets, comparing the residuals after the DD subtraction and their power spectra. We also compare results from the three Cyg A models described previously: the old, new, and point source models. All the cylindrical power spectra are shown as a ratio between the Stokes I and the thermal noise power spectrum, where the thermal noise was evaluated from time differences of Stokes V images. This ratio removes most of the effects that are caused by the sensitivity response of the instrument, allowing us to focus on the differences between the different Cyg A models.

In the next section, we sketch how errors in the sky model could propagate into the foreground-subtracted visibilities after the DD subtraction. This is an important aspect to better understand the following results.

5.1. DD gain errors due to model incompleteness

In the following sections, the effect of model incompleteness for Cyg A, after its subtraction, is illustrated both in the image domain via residual images and in the power spectrum domain. The latter is particularly relevant for 21-cm cosmology experiments on baselines where the 21-cm signal is expected and illustrates that the effect in the image domain can be extremely subtle because continuum images contain all baselines and are averaged over frequency, thereby averaging out frequency-dependent fluctuations. In the cylindrical power spectra (Sects. 5.2.1 and 5.3.1), however, the effects as a function of baseline (proportional to k_⊥) and frequency (or its Fourier dual, ‘delay’, proportional to k_∥) become much more evident.

In general, when determining LOFAR station-based gains, especially their DD gains, the incompleteness of the sky model can have a significant impact on the quality of the gains as a function of baselines and frequency (Patil et al. 2016; Ewall-Wice et al. 2017; Mouri Sardarabadi & Koopmans 2019; Mevius et al. 2022). When solving for multiple directions, the gains can modify the point spread function (PSF) used for subtraction in such a way that its sidelobes absorb unmodelled structures in the sky, thereby removing diffuse foregrounds and 21-cm emission on larger scales (see Patil et al. 2016 for a demonstration). This is why spectrally smooth gain regularisation was introduced (e.g. Yatawatta 2015; Gan et al. 2023), as it suppresses the ability of gains to remove unmodelled source structures.

Similarly, very bright sources far from the phase centre (e.g. Cas A and Cyg A) can have a dominant impact via their PSF sidelobes, leaking power into the main target field (e.g. NCP or 3C 196), as demonstrated by Gan et al. (2022) and Munshi et al. (2024). Hence, an incomplete sky model for these bright sources might lead to incorrect DD solutions in their direction. When subtracting them with their DD gains applied, the PSF-convolved source model will be incorrectly subtracted, leaving residuals in the main field of interest (NCP in our case). For fainter sources, this effect is much smaller and averages out when there are many sources in various directions and at different distances. However, Cyg A and Cas A are two extremely dominant sources, and their effect does not average away (Gan et al. 2022; Munshi et al. 2024).

The errors on these source models can then enter the power spectrum in several different and complex ways: (i) even if the DD solutions in the direction of a bright source were known perfectly, applying them to an incorrect source model would leave residuals in the main field via the sidelobes, as discussed above; (ii) in practice, however, the incomplete sky model (not just the bright source model) also leads to complex gain errors, even in the case of strong gain regularisation as a function of frequency (see e.g. Yatawatta 2015). The latter effect can redistribute power across baselines depending on which baselines exhibit the largest differences between data and model, as well as on which baselines have the most visibilities (both contributing to the minimisation of the penalty function for the gain solutions).

In the following sections, we illustrate these effects for the case of an incomplete model for Cyg A.

5.2. Results from simulations

After the DD subtraction step, we imaged the time and frequency averaged residuals at the Cyg A position in the simulated dataset for all the three tested models. In the top row of Fig. 6, from second to fourth panel, we show these residuals divided by the apparent peak brightness of Cyg A in the initially calibrated data, shown in the first panel. This gives an estimate of the relative errors using different models during the DD solving. In the second image, the model used for prediction and calibration was the same. The residuals therefore do not reflect the quality of the model, but the effectivity of calibration. The third and fourth panels show the extra residual power caused by (i) not modelling spectral indices, and (ii) assuming Cyg A is a point source. In the bottom row we also report the differences between the three residual images, to better highlight the structures affected by the different Cyg A models.

Fig. 6. Dirty images of the simulated NCP dataset in the direction of Cyg A. The top row shows, from left to right: (i) the image after initial calibration, and residual images after DD subtraction using the (ii) new, (iii) old and (iv) point source model, divided by the peak brightness of the initially calibrated image. The bottom row shows the differences between residual images obtained after DD subtraction for the three Cyg A models: Old − New (left), Point − New (middle), and Point − Old (right). All the images have been obtained with a natural weighting scheme and using a baseline range of 50 − 4500λ, resulting in the synthesised beam shown in the bottom-left corner of the top-left panel. The contour in every panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

Although the old model is incomplete and spectrally incorrect, it leaves residual similar to the new model, both with a maximum ratio of approximately 0.13 close to the source centre. The point source model shows higher residuals, with a maximum of 0.17. The removal of all the three models generated ripples far from the sources and directed toward the NCP, eventually affecting the power spectrum. The intensity of these ripples are slightly larger for the extreme point source model, which also shows stronger small-scale fluctuations in the perpendicular direction. It is interesting to note that residuals for the new model are not fully consistent with the noise, even though the same model was used for both prediction and DD calibration. This highlights that even in the perfect scenario, calibration errors still affect the resulting residuals. This can be partly attributed to the added noise and the non-linear behaviour of calibration. However, in our experience, sources in the centre of the primary beam, are easier to subtract and leave lower residuals. A possible explanation for the relatively large residuals with a perfect model might be that the source passes through nulls in the primary beam and that not all stations have sufficient signal-to-noise at all times to accurately solve the source.

Fig. 8. Power spectra along the Cyg A direction after DD subtraction (top panel) and their ratios (middle and bottom panels) for the simulated NCP dataset. The top panel shows Stokes IP(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) models, along with the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the ratio New/Old (black dots), while the bottom panel shows the ratios New/Point (pink squares), and Old/Point (purple diamonds). Filled markers represent the ratio of the cylindrical power spectra P(k⊥, k∥) for each (k⊥, k∥)-cell within the Cyg A delay lines, while the white-faced markers indicate the ratios of the power spectra P(k). For the latter, the horizontal error bars indicate the k-bin extension, while the vertical error bars indicate the 1σ standard deviation over all (k⊥, k∥)-cells of the ratios shown in Fig. 7. A small offset in k has been added to the white-faced pink squares and purple diamonds in the bottom panel to avoid overlapping of error bars. Fig. 7. Stokes I cylindrical power spectra divided by the thermal noise power spectrum at different stages of the processing of the simulated NCP dataset. Top row, from left to right: power spectrum ratio after initial calibration, after DD subtraction using the new, old, and point source model of Cyg A. The bottom row shows ratios of results from processing with different models. From left to right: the initial calibration results over new model results, new model results over old model results, new model results over point source model results, and old model results over point source model results. ‘DD New/Point’ and ‘DD Old/Point’ share the same colour bar. In all the panels, the thick line indicates the horizon delay line (foreground wedge), whereas the thin solid and dashed line pairs indicate the delay ranges where we expect most of the power of Cyg A and Cas A, respectively.

5.2.1. Cylindrical power spectra and ratios

The cylindrical power spectra P(k_⊥, k_∥), centred on the NCP, were estimated both after initial calibration and DD subtraction, and show that the DD steps removed most of the foreground emission within the wedge. The top row of Fig. 7 shows the Stokes I power spectra divided by the thermal noise, estimated from time-differenced Stokes V. The thick line represents the horizon limit, above which the EoR window should be foreground-free. The other two pairs of lines, the solid and the dashed ones, show the range where we can maximally find Cyg A and Cas A, respectively. Even with error-free simulations, the residual power within the wedge is quite high after the DD subtraction, up to 10 times the thermal noise. This might be due to calibration errors occurring during the DD solving or inaccuracies in the primary beam prediction, which also leaks power to higher k_∥, into the EoR window. Most of the residual power is along the Cyg A direction, while Cas A seems better subtracted. The maximum residual brightness of Cas A is 1% of the peak of Cas A after the initial calibration and it is the same irrespective of which Cyg A model was used in the calibration. This shows that using spatially and spectrally different models of Cyg A during the DD calibration step does not considerably affect the gain calibration of other bright sources such as Cas A. The bottom-left panel of Fig. 7 shows that the DD subtraction removes most of the power within the horizon. The residual power is usually taken into account with GPR, a step skipped in this work.

Ratios of DD subtracted power spectra for the three Cyg A models are shown in the bottom row of Fig. 7. Overall, the new model leaves the least residual power compared to the other two models. This is expected because the new model is the same model used during the prediction step. Most differences are visible around the Cyg A lines (thin solid lines), where the new model outperforms the old one, particularly at k_⊥ ≲ 0.12 h cMpc⁻¹. The point source model, however, leaves the highest residuals across all k_⊥ − k_∥ values within the Cyg A lines.

The top panel of Fig. 8 shows the 1D power spectrum P(k) for the three models, extracted from the cylindrical power spectra (i.e. top row of Fig. 7, but not divided by the thermal noise) along the Cyg A direction and averaged within k-bins of size 0.07 h cMpc⁻¹.¹⁵ The middle and bottom panels show the ratios of these P(k) pairs (white-faced markers), alongside the ratio values from bottom row of Fig. 7 for each (k_⊥, k_∥)-cell within the Cyg A delay lines (filled markers). We plotted the ratio P_new(k)/P_old(k) separately in the middle panel to make it easier to see where the new model has improved residuals compared to the old model. The 1σ uncertainties associated with P(k) have been calculated as

$$\begin{array}{c}\hfill P_{\mathrm{err}}(k)=\sqrt{\frac{{\displaystyle \sum _{(k_{\perp },k_{\Vert })\in k}}{P}_{\mathrm{err}}^2(k_{\perp },k_{\Vert })}{N_k^2}},\end{array}$$(8)

where N_k is the number of (k_⊥, k_∥)-cells falling in k-bin and P_err(k_⊥, k_∥) is the error on the cylindrical power spectrum, estimated from the sample variance. For the ratio P_i(k)/P_j(k), the sample variance is not a reliable uncertainty estimator, as all three datasets were drawn from the same sample and have the same noise. Therefore, the uncertainties of P_i(k)/P_j(k) were set to the standard deviation over all (k_⊥, k_∥)-cells of the ratios shown in the bottom row of Fig. 7. The New/Old ratio reaches a minimum value of P_new(k)/P_old(k) = 0.80 ± 0.03 at k = 0.30 h cMpc⁻¹ and is significantly less than one for all k < 0.4 h cMpc⁻¹. At higher k values, the ratio is close to one, suggesting that the new model does not offer improvement over the old model at smaller frequency and spatial scales. On the other hand, both P_new(k)/P_point(k) and P_old(k)/P_point(k) are consistently lower than one. The inverse-variance weighted averages across all k-bins are ⟨P_new(k)/P_old(k)⟩ = 0.95 ± 0.05, ⟨P_new(k)/P_point(k)⟩ ≈ ⟨P_old(k)/P_point(k)⟩ = 0.7 ± 0.1, indicating that the point source model performs worse than both the new and old models. This highlights that, although the point source model has correct spectral information, even a rough spatial characterisation of the source is necessary to improve DD calibration results.

5.2.2. Spherical power spectra and differences

Although the power spectra P(k) along the Cyg A direction give a robust estimation of the impact that the three different models have on the residual power, we saw in the third and fourth bottom panels of Fig. 7 that differences are present also in different regions of the cylindrical power spectrum. We can then use Eq. (7) to estimate the spherical power spectrum Δ²(k) in k-bins spanning the whole (k_⊥, k_∥)-space and assess the overall impact of the three models on the 21-cm upper limits.

The total Stokes I spherical power spectra are shown in the top panel of Fig. 9, evaluated within ten logarithmically spaced k-bins in the range 0.08 ≤ k ≤ 1.15 h cMpc⁻¹, with a bin size of dk/k ≈ 0.3. The differences between model pairs i − j and the associated 1σ uncertainties have been calculated as

Fig. 9. Spherical power spectra after DD subtraction (top panel) and their difference (middle and bottom panels) for the simulated NCP dataset. The top panel shows Stokes IΔ2(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) Cyg A model, and the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the difference Δold2 − Δnew2 (black dots), while the bottom panel shows the differences Δpoint2 − Δnew2 (pink squares), and Δpoint2 − Δold2 (purple diamonds), with the associated 1σ uncertainties (values reported in Table D.1). The horizontal error bars indicate the k-bin extension. We added a small offset in k to the pink and purple markers in the bottom panel to avoid overlapping of error bars. The grey shaded area delimits the range in k where Δ2 upper limits are usually extracted from LOFAR data, being the most sensitive range for LOFAR.

$$\begin{array}{c}\hfill D_{\mathit{ij}}={\mathrm{\Delta }}_i^2-{\mathrm{\Delta }}_j^2,D_{ij,\mathrm{err}}=\sqrt{{({\mathrm{\Delta }}_{i,\mathrm{err}}^2)}^2+{({\mathrm{\Delta }}_{j,\mathrm{err}}^2)}^2},\end{array}$$(9)

where Δ_err²(k) is estimated from the sample variance of the power spectrum (see P_err in Eq. 8). In this case, the sample variance is used also for the uncertainties of the differences, because D_ij is calculated directly from Δ²(k). The values are reported in Table D.1, with also the inverse-variance weighted mean of the differences calculated over two representative k-bin ranges. The difference Δ_old² − Δ_new² (middle panel of Fig. 9) is larger than zero across most of the k-bins, with just a few exceptions. However, all these values are consistent with zero when the sample variance error is taken into account, except at k = 0.35 h cMpc⁻¹ where Δ_old² > Δ_new² significantly. We observe the strongest differences at this k-bin for all the model pairs, consistent with Fig. 8. Looking at the differences with the point source model (bottom panel of Fig. 9), we find that Δ_point² is higher than Δ_new² and Δ_old² across all k, with significant values at 0.19 ≤ k ≤ 0.85 h cMpc⁻¹.

5.3. Results from real observations

For the actual observations, we follow a similar route compared to the simulated data. In the top row of Fig. 10, from second to fourth panel, we show the residuals at the Cyg A position for all the three tested models, divided by the apparent peak brightness of the initial calibrated Cyg A, shown in the first panel. Differences between the residual images are shown in the bottom row. The residuals left from the new and old model are similar, as we saw in the simulations, with a maximum ratio within the Cyg A region of 0.09 and 0.11, respectively. The point source model shows the highest residuals, with a maximum ratio of 0.23. Although the measured residuals for the new and old models are smaller than those in the simulations, the fluctuations originated from the subtraction of these models are stronger even far from the source and more complex. The subtraction of a single point component leaves the strongest sidelobes and ripples, both close and far from the source location.

Fig. 10. Dirty images of the real observed NCP dataset in the direction of Cyg A. The top row shows, from left to right: (i) the image after initial calibration, and residual images after DD subtraction using the (ii) new, (iii) old and (iv) point source model, divided by the peak brightness of the initially calibrated image. The bottom row shows the differences between residual images obtained after DD subtraction for the three Cyg A models: Old − New (left), Point − New (middle), and Point − Old (right). All the images have been obtained with a natural weighting scheme and using a baseline range of 50 − 4500λ, resulting in the synthesised beam shown in the bottom-left corner of the top-left panel. The contour in every panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

5.3.1. Cylindrical power spectra

The top row of Fig. 11 shows the Stokes I cylindrical power spectra after the initial calibration and the DD subtraction, all divided by the thermal noise power spectrum. In this real case scenario, the DD subtraction is not able to remove some foreground emission at low k_∥. The power spectrum reaches 10⁴ times the thermal noise level, because of systematics such as calibration errors and sky model incompleteness, including diffuse Galactic emission not present in the model. However, the decrease in power after the DD subtraction step is a factor of 10³ at low k_∥. Some residual power is also observed along the Cas A direction, which might be the cause of some of the ripples visible in Fig. 10.

Fig. 11. Stokes I cylindrical power spectra divided by the thermal noise power spectrum at different stages of the processing of the real observed NCP dataset. Top row, from left to right: power spectrum ratio after initial calibration, after DD subtraction using the new, old, and point source model of Cyg A. The bottom row shows ratios of results from processing with different models. From left to right: the initial calibration results over new model results, new model results over old model results, new model results over point source model results, and old model results over point source model results. ‘DD New/Point’ and ‘DD Old/Point’ share the same colour bar. In all the panels, the thick line indicates the horizon delay line (foreground wedge), whereas the thin solid and dashed line pairs indicate the delay ranges where we expect most of the power of Cyg A and Cas A, respectively.

In the bottom row of Fig. 11 (second to fourth panel), we show the ratios of the DD subtracted power spectra for the three Cyg A models. In all three cases, most differences are confined along the Cyg A direction. Although residuals from Cas A are observed in the top panels, the three different Cyg A models do not affect them differently, and the power along the Cas A direction cancels out in the ratio.

The 1D power spectra P(k) along the Cyg A direction for the three models are shown in the top panel of Fig. 12. The ratio P_new(k)/P_old(k) is shown in the middle panel, while P_new(k)/P_point(k) and P_old(k)/P_point(k) are shown in the bottom panel. In this real case, the 1σ uncertainties associated with ratio P_i(k)/P_j(k) have been estimated from the sample variance as

Fig. 12. Power spectra along the Cyg A direction after DD subtraction (top panel) and their ratios (middle and bottom panels) for the real observed NCP dataset. The top panel shows Stokes IP(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) models, along with the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the ratio New/Old (black points), while the bottom panel shows the ratios New/Point (pink squares), and Old/Point (purple diamonds). Filled markers represent the ratio of the cylindrical power spectra P(k⊥, k∥) for each (k⊥, k∥)-cell within the Cyg A delay lines, while the white-faced markers indicate the ratios of the power spectra P(k). For the latter, the horizontal error bars indicate the k-bin extension, while the vertical error bars indicate the 1σ uncertainties of the ratio. A small offset in k has been added to the white-faced pink squares and purple diamonds in the bottom panel to avoid overlapping of error bars.

$$\begin{array}{c}\hfill {\left[\frac{P_i}{P_j}\right]}_{\mathrm{err}}=\frac{P_i}{P_j}\sqrt{{\left(\frac{P_{i,\mathrm{err}}}{P_i}\right)}^2+{\left(\frac{P_{j,\mathrm{err}}}{P_j}\right)}^2},\end{array}$$(10)

where P_err is from Eq. (8). In real observations, the new model results in a residual power spectrum similar to the old model at low k. The ratio P_new(k)/P_old(k) reaches a minimum of 0.8 ± 0.1 at k = 0.89 h cMpc⁻¹, with an inverse-variance weighted average of 0.91 ± 0.05 over the entire k range. In contrast, the point source model shows higher power at all scales, with ⟨P_new(k)/P_point(k)⟩ = 0.67 ± 0.03 and ⟨P_old(k)/P_point(k)⟩ = 0.74 ± 0.04. These values are consistent with those measured from the simulations. Interestingly, the new model leaves lower residuals than the old one at high k, whereas in the simulations we observed the opposite trend, with larger differences at small k. This discrepancy results from the DD calibration, which successfully removed most of the Cyg A power at k_⊥ < 0.12 h cMpc⁻¹ in the observed dataset (see Fig. 11), but not in the simulated one, where the cylindrical power spectra where noise-like at higher k_⊥ (see Fig. 7).

5.3.2. Spherical power spectra and differences

The total Stokes I spherical power spectra Δ²(k) are shown in the top panel of Fig. 13. As already seen in the top panel of Fig. 11, the power is a few orders of magnitude higher than the thermal noise at k < 0.3 h cMpc⁻¹, which corresponds to approximately k_∥ < 0.3 h cMpc⁻¹ in the cylindrical power spectrum. Here, the foreground emission within and close to the main lobe of the primary beam has not been completely removed, because of either calibration errors or an incomplete sky model (e.g. extended Galactic emission is not removed and may affect low k_⊥). At higher k, Δ² is closer to the thermal noise, even if some excess power is present because of poorly subtracted off-axis sources, mainly Cyg A and Cas A. In fact, the difference between model pairs, shown in the middle and bottom panels of Fig. 13 and evaluated with Eq. (9), increases in these k-bins, reaching peaks of ≈(1600 mK)² for Δ_old² − Δ_new² and ≈(4500 mK)² for both Δ_point² − Δ_new² and Δ_point² − Δ_old² at k = 0.85 h cMpc⁻¹. The real observed data do not prefer one model at k ≤ 0.35 h cMpc⁻¹, where all the differences have large uncertainties and are consistent with zero. This is in contrast with the simulations, where differences and uncertainties were smaller. This happens because in the observed data there are other systematics, such as model incompleteness, calibration errors, and inaccuracies in the primary beam modelling, that play a bigger role during the DD subtraction. Given the current situation of the systematics, we do not expect to see larger differences by using different Cyg A models. However, an overall improvement is achievable by using the new model instead of the old one in the range 0.63 − 1.15 h cMpc⁻¹, where ⟨Δ_old² − Δ_new²⟩=(1259 ± 742 mK)². At lower k, the average difference is consistent with zero, but the expected value is ≈(257 mK)² in the range 0.08 − 0.63 h cMpc⁻¹. This shows that, in principle, the new Cyg A model is capable of reducing the excess power by hundreds of mK, even if Δ_new² is still 2–3 times higher than the thermal noise level.

Fig. 13. Spherical power spectra after DD subtraction (top panel) and their difference (middle and bottom panels) for the real observed NCP dataset. The top panel shows Stokes IΔ2(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) Cyg A model, and the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the difference Δold2 − Δnew2 (black dots), while the bottom panel shows the differences Δpoint2 − Δnew2 (pink squares), and Δpoint2 − Δold2 (purple diamonds), with the associated 1σ uncertainties (values reported in Table D.1). The horizontal error bars indicate the k-bin extension. We added a small offset in k to the pink and purple markers in the bottom panel to avoid overlapping of error bars. The grey shaded area delimits the range in k where Δ2 upper limits are usually extracted from LOFAR data, being the most sensitive range for LOFAR.

6. Discussion and conclusions

The latest upper limits on the 21-cm power spectrum from the LOFAR-EoR KSP show an excess power in Δ²(k) with respect to thermal noise (Mertens et al. 2020). Gan et al. (2022) analysed various causes of such an excess, concluding that it might be dominated by residuals of bright sources, especially those in the far field. Their power spectra, even after GPR foreground removal, show non-negligible power in the Cas A and Cyg A direction, hinting that the modelling and subtraction of these A-team sources has been performed poorly. The models of these sources, also used by the current LOFAR-EoR processing pipeline, are low resolution, with no spectral variation and incorrect spectral indices. These models have been obtained from VLA data at ≈74 MHz and assigned a fixed spectral index of −0.8 to each component, which is used to extrapolate values to the LOFAR HBA band where the 21-cm power spectrum is calculated.

In this work, we focused on Cyg A, the brightest of the A-team at 150 MHz, and used it as a test case to assess the impact of a better spatial and physical model on the 21-cm power spectrum. In Sect. 2, we introduced what we called the ‘old’ Cyg A model and showed that extrapolating it to the HBA frequency range results in 16% missing flux. Not including any spectral variation across the source is incorrect, because we know from MK16 that the brightest hotspots (A and D) show a frequency turnover at around 140–150 MHz. To solve these issues, we made a new model directly from HBA observations (110–250 MHz), taking advantage of the high spatial resolution allowed by the Dutch LOFAR configuration (i.e. ≈5 arcsec) and using a second-order logarithmic polynomial to model the spectral behaviour of each component of the source (Sect. 2.1). We used the forced-spectrum method of WSCLEAN to transfer this spectral information into the model during a multi-scale and multi-frequency deconvolution (Ceccotti et al. 2023). This makes our new Cyg A model the first test case of such a method for improving 21-cm power spectrum results.

In Sect. 5, we compared the performances of the old model with our new one within the LOFAR-EoR processing pipeline (Sect. 4), using both simulated and observed data of the NCP field. Additionally, we considered an extreme Cyg A model, constituted by a single point source at the Cyg A location but with correct spectral information. On the simulated data, the most significant differences appear in P(k) along the Cyg A direction at k < 0.4 h cMpc⁻¹, where the new model performs better than both the old and point source models. These differences are washed out in Δ²(k), where all (k_⊥, k_∥)-cells are considered, resulting in Δ_new ≈ Δ_old across most k bins. For the real observed data, differences between the three models in both P(k) along the Cyg A direction and Δ²(k) are less pronounced at k < 0.4 h cMpc⁻¹, where uncertainties make their power spectra consistent with one another. Larger differences are observed at higher k, where the point source model results in significantly higher power in both the power spectrum estimators. The lowest residuals are achieved with the new model at k > 0.6 h cMpc⁻¹, where the old and point source models show an average excess of more than (1250 mK)² and (2300 mK)², respectively.

As we did not inject any systematics in the simulations, apart from adding the thermal noise, we expect all the differences to come from the different Cyg A models. In the real observations, where other systematics are present, these systematics have a larger impact on the residuals, hiding the effect of the different models, especially at the lower k where upper limits on the 21-cm power spectrum are usually extracted. We know that Cyg A is the strongest source at 150 MHz, but not in apparent flux for NCP observations. There are many other bright sources outside the main lobe of the primary beam that could leave subtraction residuals strong enough to outweigh those left from Cyg A, most notably Cas A, which has an apparent flux approximately three times higher than Cyg A. The DD subtracted plots in Fig. 11 show that some residuals associated with Cas A are present, highlighting the need of a better model for this source as well. However, as long as we focus our analysis at k < 0.6 h cMpc⁻¹, as in the analysis of real data, we do not expect significant improvements until the effects of other systematics are lowered.

In conclusion, we have demonstrated that the forced-spectrum method can generate spectrally accurate models of complex sources such as Cyg A, provided that physical spectral index maps are available or can be extracted from the data. This approach ensures that robust spectral information is embedded in a high-resolution model, which is fundamental for using such a model at different frequencies from those at which the model was originally extracted. For example, we showed that while the old model has a correct flux density at 74 MHz, its incorrect spectral index results in a much lower extrapolated flux density at 150 MHz than the expected value. Moreover, the new model is described as a standard sky model, with a low number of point and Gaussian components with spectral information, and can therefore directly be used in calibration. The new Cyg A model presented in this work, and more generally, the forced-spectrum fitting method, could greatly benefit high-resolution and wide-field interferometers operating at low frequencies, such as the upcoming SKA-Low. This advancement will be beneficial for enhancing the accuracy and reliability of astronomical observations at low frequencies.

Because no GPR analysis was applied before power spectrum estimation, all the improvements that we measured in the 21-cm power spectrum come from the DD calibration, the only step where we used Cyg A. Alongside with assessing the impact of a better model on the power spectrum, we demonstrated that the calibration itself can benefit from a forced-spectrum model. We cannot make a one-to-one comparison with Mertens et al. (2020) about the excess power in the LOFAR power spectrum, because we did not apply a GPR technique. However, we can conclude that improving the modelling of Cyg A can reduce the excess power up to (250 mK)² over the k = 0.08 − 0.63 h cMpc⁻¹ range, as measured in the real observations (see Table D.1). This reduction can be smaller after GPR, because some of the extra power that we observed would be absorbed in the foreground mode-mixing kernel, making the difference between the three models less strong. Our simulation scenario can be expected to be comparable to the results of a GPR analysis, because we have not included the systematics that GPR removes. Here, indeed, the measured differences are smaller, with ⟨Δ_old² − Δ_new²⟩≈(17 mK)² over the k = 0.08 − 0.63 h cMpc⁻¹ range. What we showed is that using both spatially and spectrally accurate models of bright off-axis sources is important to lower the upper limits on the 21-cm power spectrum.

In order to see a significant impact of a better spectral modelling in real data, other systematics must be reduced first. Because we flagged data after the DD subtraction, we expect even low-level RFI to be removed and they would therefore not be a dominant source of excess power. Another possible cause of the excess power might be ionospheric effects that have not been completely solved for in the initial calibration step. Although, Gan et al. (2022) and Brackenhoff et al. (2024) concluded that the ionosphere is currently not the dominant source of the excess power. In the analysis of Brackenhoff et al. (2024), it is shown that rapid spatial and temporal variations in the primary beam model could, however, result in incorrect DD calibration solutions. This is especially true for directions far from the field centre, where primary beam sidelobes move faster in time and have steep edges with deep nulls. Whereas we do not see incorrect DD gains, we see bright residuals associated with Cyg A even in the simulations where the calibration model is the same as the one used in the prediction step (first panel of Fig. 6). This could mean that the DD calibration gains are not able to capture the finest spatial, temporal, and spectral variations of the source due to the primary beam sidelobes. Understanding and characterising the time and frequency behaviour of the LOFAR HBA primary beam appears to be the next important step in trying to reduce the excess power in the NCP field. Processing and analysing a different sky region, such as the 3C 196 field, could be helpful as well in understanding the source of the systematics (Ceccotti et al., in prep.).

Data availability

The component list of the new model of Cygnus A is available at https://doi.org/10.5281/zenodo.14863290. The Cygnus A dataset (Project code: LC1_013) and both the observed and simulated NCP datasets (Observation ID: L246309) underlying this article will be shared upon reasonable request to the corresponding author.

Acknowledgments

EC, ARO and LVEK would like to acknowledge support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen. EC acknowledges support from the Ministry of Universities and Research (MUR) through the PRIN project ‘Optimal inference from radio images of the epoch of reionization’. SAB, SG, LVEK 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’). FGM acknowledges support from a PSL Fellowship. RG acknowledges support from SERB, DST Ramanujan Fellowship no. RJF/2022/000141. The post-doctoral contract of IH was funded by Sorbonne Université in the framework of the Initiative Physique des Infinis (IDEX SUPER). LOFAR, the Low Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. In this work, we made use of the KVIS (Gooch 1996) and DS9 (Joye & Mandel 2003) FITS file image viewers, and the ASTROPY (Astropy Collaboration 2022), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), PANDAS (McKinney 2010), SCIPY (Virtanen et al. 2020) PYTHON packages.

References

Appendix A: Old Cygnus A model components

In Table A.1 we report the components of the old Cyg A model used by the current LOFAR-EoR processing pipeline. The model catalogue is available at https://raw.githubusercontent.com/lofar-astron/prefactor/master/skymodels/Ateam_LBA_CC.skymodel.

Appendix B: Processing of the Cygnus A dataset

The initial dataset consisted of 399 and 203 SBs for HBA-low and HBA-high, respectively, with frequency and time resolutions of 48.8 kHz (i.e. 4 channels per SB) and 2 s. In the following, we will briefly describe the flagging, averaging, and calibration that have been previously performed using CASA (Common Astronomy Software Applications; McMullin et al. 2007).

Firstly, the dataset has been averaged to a single channel per SB (i.e. 195.3 kHz of frequency resolution) and to an integration time of 2 s. The SB at the edges of both bands are flagged because of strong RFI. Then, the calibration of the HBA-low data has been performed using the best model from McKean et al. (2016), solving both phase and amplitude for each channel (i.e. one solution every 195.3 kHz). The calibrated data has been imaged with the CASA multiscale multifrequency synthesis deconvolution, decomposing the spectral variation with two terms of the logarithmic polynomial function (see Eq. 1). This model is used to calibrate the HBA-high data, which are then concatenated to the HBA-low data.

The full dataset has been self-calibrated with a two Taylor term model obtained from the combined bands and finally flagged, resulting in a total of 412 SBs spanning 111–181 and 210–249 MHz.

Appendix C: Simple NCP sky model

In Table C.1, we report the central position and total flux density of the each of the 14 clusters that constitute the simple NCP sky model. We used this model for the simulated dataset.

Appendix D: Final power spectrum differences

In Table D.1, we report the differences D_ij = Δ_i² − Δ_j² and the associated 1σ uncertainties evaluated using Eq. (9). Both the results from simulations (Sect. 5.2.2) and real observations (Sect. 5.3.2) are shown. In the last two rows, we also report the inverse-variance weighted mean of the differences D_ij over the k-bins ranges and the associated 1σ uncertainties, evaluated as

$$\begin{array}{c}\hfill ⟨D_{\mathit{ij}}⟩=\frac{{\displaystyle {\sum }_k}(D_{\mathit{ij}}(k)·D_{ij,\mathrm{err}}^{-2}(k))}{{\displaystyle {\sum }_k}{D}_{ij,\mathrm{err}}^{-2}(k)},{⟨D_{\mathit{ij}}⟩}_{\mathrm{err}}=\sqrt{\frac{1}{{\displaystyle \sum _k}{D}_{ij,\mathrm{err}}^{-2}}}.\end{array}$$(D.1)

All Tables

All Figures

Fig. 1. High-resolution, frequency-integrated image of Cyg A between 111 and 249 MHz (left), and peak brightness of the hotspots A (top right) and D (bottom right), extracted from the high-resolution images of the 100 output channels (black solid line, with 1σ uncertainties shown as a grey shaded area). The fits using a second-order (blue dashed line) and a first-order logarithmic polynomial function (red dot-dashed line) are also shown, with the best fit parameters reported in the same colours. The latter has been fitted only between 111 and 181 MHz for a consistent comparison with MK16. The image on the left has a noise level of σ = 10 mJy/beam, with contours starting at 1 Jy/beam and increasing by a factor of 2.

In the text

	Fig. 2. Forced-spectrum input spectral index (left) and curvature maps (right) at 190.18 MHz, extracted from the 5 arcsec resolution images. Contours in both figures are the same of the left panel of Fig. 1. The colour scale of the spectral index map spans the range −1.27 ≤ α ≤ −0.38, while the curvature is plotted in the range −1.04 ≤ β ≤ 1.04.
In the text

Fig. 3. Flux density between 50 and 270 MHz measured from the Cyg A models for the total emission (top), hotspot A (middle), and hotspot D (bottom). The top panel shows the old model (solid blue line), the new model before (dashed green line) and after spectral rescaling (solid black line), and the third-order logarithmic polynomial used by MK16 (dashed orange line). In the middle and bottom panels, only the flux densities from the new Cyg A models before and after rescaling are shown. These values were extracted from a circular region with a 5 arcsec diameter, centred on the peak of each hotspot. The filled points indicate the spectral peaks, and the turnover frequency after rescaling is also reported. The green shaded area in all panels marks the LOFAR HBA range used to extract the new Cyg A model.

In the text

	Fig. 4. New (top) and old model (bottom) of Cyg A, rendered at 150 MHz. Contours are the same of Fig. 1. The zero-length axes of the old model Gaussians have been replaced with 1 arcsec axes for imaging purposes, resulting in the appearance of vertical spikes.
In the text

Fig. 5. All sky dirty image of the simulated NCP dataset (left), with zoom-in to the NCP main field (top right) and the Cyg A direction (bottom right). The all sky image has been obtained using a baseline range of 50 − 250λ, the NCP zoom-in with 50 − 500λ, while the Cyg A zoom-in with 50 − 4500λ. A natural weighting scheme has been used for the all sky and the NCP zoom-in images, whereas a uniform weighting has been chosen for the Cyg A zoom-in to make the source structures more visible. The contour in the bottom right panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

In the text

Fig. 6. Dirty images of the simulated NCP dataset in the direction of Cyg A. The top row shows, from left to right: (i) the image after initial calibration, and residual images after DD subtraction using the (ii) new, (iii) old and (iv) point source model, divided by the peak brightness of the initially calibrated image. The bottom row shows the differences between residual images obtained after DD subtraction for the three Cyg A models: Old − New (left), Point − New (middle), and Point − Old (right). All the images have been obtained with a natural weighting scheme and using a baseline range of 50 − 4500λ, resulting in the synthesised beam shown in the bottom-left corner of the top-left panel. The contour in every panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

In the text

Fig. 8. Power spectra along the Cyg A direction after DD subtraction (top panel) and their ratios (middle and bottom panels) for the simulated NCP dataset. The top panel shows Stokes IP(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) models, along with the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the ratio New/Old (black dots), while the bottom panel shows the ratios New/Point (pink squares), and Old/Point (purple diamonds). Filled markers represent the ratio of the cylindrical power spectra P(k⊥, k∥) for each (k⊥, k∥)-cell within the Cyg A delay lines, while the white-faced markers indicate the ratios of the power spectra P(k). For the latter, the horizontal error bars indicate the k-bin extension, while the vertical error bars indicate the 1σ standard deviation over all (k⊥, k∥)-cells of the ratios shown in Fig. 7. A small offset in k has been added to the white-faced pink squares and purple diamonds in the bottom panel to avoid overlapping of error bars. Fig. 7. Stokes I cylindrical power spectra divided by the thermal noise power spectrum at different stages of the processing of the simulated NCP dataset. Top row, from left to right: power spectrum ratio after initial calibration, after DD subtraction using the new, old, and point source model of Cyg A. The bottom row shows ratios of results from processing with different models. From left to right: the initial calibration results over new model results, new model results over old model results, new model results over point source model results, and old model results over point source model results. ‘DD New/Point’ and ‘DD Old/Point’ share the same colour bar. In all the panels, the thick line indicates the horizon delay line (foreground wedge), whereas the thin solid and dashed line pairs indicate the delay ranges where we expect most of the power of Cyg A and Cas A, respectively.

In the text

Fig. 7. Stokes I cylindrical power spectra divided by the thermal noise power spectrum at different stages of the processing of the simulated NCP dataset. Top row, from left to right: power spectrum ratio after initial calibration, after DD subtraction using the new, old, and point source model of Cyg A. The bottom row shows ratios of results from processing with different models. From left to right: the initial calibration results over new model results, new model results over old model results, new model results over point source model results, and old model results over point source model results. ‘DD New/Point’ and ‘DD Old/Point’ share the same colour bar. In all the panels, the thick line indicates the horizon delay line (foreground wedge), whereas the thin solid and dashed line pairs indicate the delay ranges where we expect most of the power of Cyg A and Cas A, respectively.

In the text

Fig. 9. Spherical power spectra after DD subtraction (top panel) and their difference (middle and bottom panels) for the simulated NCP dataset. The top panel shows Stokes IΔ2(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) Cyg A model, and the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the difference Δold2 − Δnew2 (black dots), while the bottom panel shows the differences Δpoint2 − Δnew2 (pink squares), and Δpoint2 − Δold2 (purple diamonds), with the associated 1σ uncertainties (values reported in Table D.1). The horizontal error bars indicate the k-bin extension. We added a small offset in k to the pink and purple markers in the bottom panel to avoid overlapping of error bars. The grey shaded area delimits the range in k where Δ2 upper limits are usually extracted from LOFAR data, being the most sensitive range for LOFAR.

In the text

Fig. 10. Dirty images of the real observed NCP dataset in the direction of Cyg A. The top row shows, from left to right: (i) the image after initial calibration, and residual images after DD subtraction using the (ii) new, (iii) old and (iv) point source model, divided by the peak brightness of the initially calibrated image. The bottom row shows the differences between residual images obtained after DD subtraction for the three Cyg A models: Old − New (left), Point − New (middle), and Point − Old (right). All the images have been obtained with a natural weighting scheme and using a baseline range of 50 − 4500λ, resulting in the synthesised beam shown in the bottom-left corner of the top-left panel. The contour in every panel is the 1 Jy/beam level from Fig. 1 to show the position and extension of Cyg A.

In the text

Fig. 11. Stokes I cylindrical power spectra divided by the thermal noise power spectrum at different stages of the processing of the real observed NCP dataset. Top row, from left to right: power spectrum ratio after initial calibration, after DD subtraction using the new, old, and point source model of Cyg A. The bottom row shows ratios of results from processing with different models. From left to right: the initial calibration results over new model results, new model results over old model results, new model results over point source model results, and old model results over point source model results. ‘DD New/Point’ and ‘DD Old/Point’ share the same colour bar. In all the panels, the thick line indicates the horizon delay line (foreground wedge), whereas the thin solid and dashed line pairs indicate the delay ranges where we expect most of the power of Cyg A and Cas A, respectively.

In the text

Fig. 12. Power spectra along the Cyg A direction after DD subtraction (top panel) and their ratios (middle and bottom panels) for the real observed NCP dataset. The top panel shows Stokes IP(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) models, along with the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the ratio New/Old (black points), while the bottom panel shows the ratios New/Point (pink squares), and Old/Point (purple diamonds). Filled markers represent the ratio of the cylindrical power spectra P(k⊥, k∥) for each (k⊥, k∥)-cell within the Cyg A delay lines, while the white-faced markers indicate the ratios of the power spectra P(k). For the latter, the horizontal error bars indicate the k-bin extension, while the vertical error bars indicate the 1σ uncertainties of the ratio. A small offset in k has been added to the white-faced pink squares and purple diamonds in the bottom panel to avoid overlapping of error bars.

In the text

Fig. 13. Spherical power spectra after DD subtraction (top panel) and their difference (middle and bottom panels) for the real observed NCP dataset. The top panel shows Stokes IΔ2(k) after the DD subtraction step using the new (blue), old (orange), and point source (green) Cyg A model, and the thermal noise level (dashed grey line). The shaded areas represent the 1σ uncertainties. The middle panel shows the difference Δold2 − Δnew2 (black dots), while the bottom panel shows the differences Δpoint2 − Δnew2 (pink squares), and Δpoint2 − Δold2 (purple diamonds), with the associated 1σ uncertainties (values reported in Table D.1). The horizontal error bars indicate the k-bin extension. We added a small offset in k to the pink and purple markers in the bottom panel to avoid overlapping of error bars. The grey shaded area delimits the range in k where Δ2 upper limits are usually extracted from LOFAR data, being the most sensitive range for LOFAR.

In the text