Twentyonecentimetre signals from the Epoch of Reionization EoR are expected to be detected in the lowfrequency radio window by the nextgeneration interferometers particularly the Square Kilometre Array SKA. However precision data analysis pipelines are required to minimize the systematics within an infinitesimal error budget. Consequently there is a growing need to characterize the sources of errors in EoR analysis. In this study we identify one such error origin namely source blending which is introduced by the overlap of objects in the densely populated observing sky under SKA1Lows unprecedented sensitivity and resolution and evaluate its twofold impact in both the spatial and frequency domains using a novel hybrid evaluation HEVAL pipeline combining endtoend simulation with an analytic method to mimic EoR analysis pipelines. Sky models corrupted by source blending induce small but severe frequencydependent calibration errors when coupled with astronomical foregrounds impeding EoR parameter inference with strong additive residuals in the twodimensional power spectrum space. We report that additive residuals from poor calibration against sky models with blending ratios of 5 and 0.5 per cent significantly contaminate the EoR window. In contrast the sky model with a 0.05 per cent blending ratio leaves little residual imprint within the EoR window therefore identifying a blending tolerance at approximately 0.05 per cent. Given that the SKA observing sky is estimated to suffer from an extended level of blending strategies involving deblending frequencydependent error mitigation or a combination of both are required to effectively attenuate the calibration impact of sourceblending defects.
['Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies', '85-04', '85-08', '85-10', 'D.2.0', 'J.2', 'I.1.2', 'I.6.5', 'I.6.6']
An evaluation of sourceblending impact on the calibration of SKA EoR experiments
{'An evaluation of source-blending impact on the calibration of SKA EoR experiments': 'Chenxi Shan, 1 ★ Haiguang Xu, 1 ★ Yongkai Zhu, 1 Yuanyuan Zhao, 1 Sarah V. White, 2 Jack L. B. Line, 3 , 4 Dongchao Zheng, 1 Zhenghao Zhu, 5 Dan Hu, 6 Zhongli Zhang, 5 , 7 and Xiangping Wu 8  \n- 1 School of Physics & Astronomy, Shanghai Jiao Tong University, Shanghai, China\n- 2 Department of Physics and Electronics, Rhodes University, PO Box 94, Makhanda, 6140, South Africa\n- 3 International Centre for Radio Astronomy Research, Curtin University, Perth, WA 6102, Australia\n- 4 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO-3D)\n- 5 Shanghai Astronomical Observatory, CAS, 80 Nandan Road, Shanghai, China\n- 6 Department of Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University, Kotlářská 2, Brno, 611 37, Czech Republic\n- 7 Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, 20A Datun Road, Beĳing 100012, China\n- 8 National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Beĳing 100012, China  \nAccepted 2024 September 08. Received 2024 September 05; in original form 2023 December 30', 'ABSTRACT': "Twenty-one-centimetre signals from the Epoch of Reionization (EoR) are expected to be detected in the low-frequency radio window by the next-generation interferometers, particularly the Square Kilometre Array (SKA). However, precision data analysis pipelines are required to minimize the systematics within an infinitesimal error budget. Consequently, there is a growing need to characterize the sources of errors in EoR analysis. In this study, we identify one such error origin, namely source blending, which is introduced by the overlap of objects in the densely populated observing sky under SKA1-Low's unprecedented sensitivity and resolution, and evaluate its two-fold impact in both the spatial and frequency domains using a novel hybrid evaluation (HEVAL) pipeline combining end-to-end simulation with an analytic method to mimic EoR analysis pipelines. Sky models corrupted by source blending induce small but severe frequency-dependent calibration errors when coupled with astronomical foregrounds, impeding EoR parameter inference with strong additive residuals in the two-dimensional power spectrum space. We report that additive residuals from poor calibration against sky models with blending ratios of 5 and 0.5 per cent significantly contaminate the EoR window. In contrast, the sky model with a 0.05 per cent blending ratio leaves little residual imprint within the EoR window, therefore identifying a blending tolerance at approximately 0.05 per cent. Given that the SKA observing sky is estimated to suffer from an extended level of blending, strategies involving de-blending, frequency-dependent error mitigation, or a combination of both, are required to effectively attenuate the calibration impact of source-blending defects.  \nKey words: instrumentation: interferometers - dark ages, reionization, first stars - techniques: interferometric - software: simulations - radio continuum: general", '1 INTRODUCTION': "The low-frequency radio window opens up new opportunities for probing the high-redshift Universe through its unique capability to receive the redshifted 21-cm hyperfine line emission of neutral hydrogen from the early stages of the Universe, which can probe the untethered Cosmic Dawn ( 𝑧 ∼ 15 to 30) and the Epoch of Reionization (EoR; 𝑧 ∼ 6 to 15) with unprecedented precision (see Furlanetto et al. 2006; Pritchard & Loeb 2012; Zaroubi 2013; Furlanetto 2016, for reviews). The faint nature of the 21-cm signal demands instruments with extremely high sensitivity, wide field-of-view (FoV), and wideband coverage, making radio interferometers favoured for statistical detection of the 21-cm signal (Liu & Shaw 2020), presumably in the power spectrum (PS) space, given their flexibility to be designed  \nand constructed as an array that meets these technical demands while also offering precise baseline coverage of EoR-specific scales (0 . 1 ≲ 𝑘 ≲ 2 Mpc -1 ). Major interferometers in the low-frequency band include the Giant Metrewave Radio Telescope (GMRT 1 ; Swarup 1991), the Murchison Widefield Array (MWA 2 ; Bowman et al. 2013; Tingay et al. 2013), the LOw Frequency ARray (LOFAR 3 ; van Haarlem et al. 2013), the MIT Epoch of Reionization (MITEoR 4 ; Zheng et al. 2014), the Precision Array for Probing the Epoch of Reionization (PAPER 5 ; Parsons et al. 2010), and next-generation  \n- 1 http://www.ncra.tifr.res.in/ncra/gmrt\n- 2 https://www.mwatelescope.org\n- 3 https://www.astron.nl/telescopes/lofar\n- 4 https://space.mit.edu/home/tegmark/main\\_omniscope.html  \narrays, the Hydrogen Epoch of Reionization Array (HERA 6 ; DeBoer et al. 2017) and the Square Kilometre Array (SKA 7 ; Mellema et al. 2013; Koopmans et al. 2015). However, owing to the complex baseline designs and complicated instrumental effects of these advanced arrays, precision data analysis pipelines are required for EoR experiments to overcome calibration, imaging, and analysis challenges and to minimize exposure to systematics, such as calibration biases (Datta et al. 2009; Grobler et al. 2014; Wĳnholds et al. 2016; Patil et al. 2016; Gehlot et al. 2021), polarization leakages (Moore et al. 2017; Hales 2017; Dillon et al. 2018; Gehlot et al. 2018), wide-FoV and wide-band imaging errors (Bhatnagar et al. 2013; Jagannathan et al. 2015; Rau et al. 2016; Ye et al. 2022) and foreground residuals (Trott et al. 2012; Chapman et al. 2016; Nasirudin et al. 2020; Hothi et al. 2021). Therefore, it has become clear that EoR experiments must operate under an infinitesimal error budget with rigorous systematic characterizations to detect 21-cm signals successfully (see Liu & Shaw 2020; Shaw et al. 2023, for reviews).  \nAdvanced low-frequency interferometers, such as LOFAR and SKA, can reach new detection limits that transform the occupation of the observed low-frequency radio sky from sparsely distributed, relatively rare sources (e.g. bright radio quasars and radio galaxies) to densely distributed populations constituting the majority of the extragalactic sky. In particular, the upcoming SKA1-Low is estimated to probe extragalactic discrete radio source (EDRS) populations deep into the faint sub-mJy radio sky ( ∼ 4400 deg -2 above 0 . 1 mJy at 150 MHz) and detect the bulk of galaxies with radio power originating from star-forming processes, the central active galactic nucleus (AGN), and a composite of both (Padovani 2016). This increase in sensitivity comes with increasing chances of overlap between the surface brightness distributions of sources (known as blending), which introduces systematic effects due to contaminated measurements (see Melchior et al. 2021, for a review). Without proper identification and mitigation strategies, these blended measurements will bias the interferometric calibration and, therefore, contribute to the overall EoR error budget by introducing imperfections during the construction of sky models, which are required to overcome line-of-sight (LoS) effects and imperfect instrumental responses (Liu et al. 2010; Li et al. 2018; Liu & Shaw 2020; Shaw et al. 2023) by most, if not all, calibration strategies of EoR experiments, including sky-based (Pearson & Readhead 1984; Rau et al. 2009), array-based (Noordam &de Bruyn 1982; Wieringa 1992), and hybrid calibration methods (Sievers 2017; Kern et al. 2020; Zhang et al. 2020; Byrne et al. 2021).  \nOf all the contributions to the tight EoR error budget, calibration errors, especially chromatic residual gain errors originating from sky-model defects, are of great importance because biases from poor calibration 8 couple with contamination from strong astronomical foregrounds ( ∼ 10 4 -10 5 brighter than the 21-cm signal), propagate further down the analysis pipeline by leaving additive residuals in the measurement space, and ultimately impede the parameter inference of the EoR signal (Datta et al. 2010; Trott & Wayth 2016; Barry et al. 2016; Ewall-Wice et al. 2017; Byrne et al. 2019; Mazumder et al. 2022). Given the propagation nature of calibration errors, quantification of the impact of sky-model defects is achieved by mimicking full EoR analysis pipelines (see He et al. 2023, for an overview on EoR pipelines) and inferring the propagated residuals  \n- 6 http://reionization.org\n- 7 https://www.skao.int  \n8 Weuse the term 'poor calibration' to describe scenarios where an inaccurate or imperfect calibration is performed, in contrast to a perfect calibration. Poor calibrations result in an ill-calibrated instrument, thereby introducing calibration errors.  \nin the measurement space, particularly in the two-dimensional (2D) PS space. Although the exact path for establishing the calibration propagation effects in the measurement space varies among studies, we can primarily sort the existing efforts into two types based on the approach used to mimic an EoR analysis pipeline. The first type uses an end-to-end approach to simulate propagation effects (e.g. Barry et al. 2016, hereafter B16), whereas the second type directly drives the final PS biases via analytic analysis (e.g. Ewall-Wice et al. 2017, hereafter E17). Despite the different approaches adopted, both types of methods have demonstrated that chromatic residual gain errors originating from sky-model imperfections severely damage the detectability of the 21-cm signal by introducing strong contaminations owing to poor calibration (one to two orders of magnitude brighter than the EoR signal, E17) and require an extremely small tolerance for model defects ( ∼ 10 -5 of spectral features, B16). Therefore, even though sky models can never be perfect, defects should always be evaluated and factored into the overall EoR error budget.  \nWith its unprecedented sensitivity and baseline coverage in the low-frequency window, the SKA1-Low will contribute primary sky models for its observing sky (Trott & Wayth 2017). Thus, it is crucial to identify sky-model defects for the SKA and properly evaluate the impact of these defects on SKA EoR experiments under the telescope's observation specifications and conditions. Given the expected observing depth enabled by the SKA1-Low ( ∼ 0 . 1 mJy at 150 MHz, Prandoni & Seymour 2015), source blending will undoubtedly introduce one such defect to sky-model constructions of SKA EoR experiments and impact the EoR detections via blending-induced calibration errors. Since LOFAR, a pathfinder of SKA1-Low with comparable spatial resolution but lower sensitivity, has already suffered from source-blending effects to an extended level ( ∼ 3 per cent under a source density ∼ 3100 deg -2 , Kondapally et al. 2021), we expect the fraction of blended sources observed by the SKA1-Low to be at the same level or higher than that of the LOFAR observations, but lower than that of the optical surveys.  \nTo date, there is still a lack of investigations that evaluate the impact of source blending on interferometric calibration, estimate the blending-induced calibration residual error, and identify blending tolerance for SKA sky-model construction. In this study, we present the first systematic assessment of source blending in the low-frequency radio window for the upcoming SKA1-Low, identify a largely two-fold effect of blending-induced uncertainties in both the spatial and frequency domains, and investigate the impact of frequency-dependent errors originating from calibrations against a sky model corrupted by source blending (hereafter blending-corrupted sky model) for SKA EoR experiments. In contrast to the simulation and analytic approaches used to infer calibration effects by B16 and E17, respectively, we propose a novel hybrid evaluation (HEVAL) approach combining custom end-to-end simulations and analytic calibration analysis mimicking SKA EoR experiment pipelines: (i) the end-to-end modules of our hybrid evaluation pipeline simulate radio maps of the low-frequency radio sky down to the smallest scales resolvable by the SKA1-Low ( ∼ 6 arcsec at 196 MHz for the foregrounds and ∼ 0 . 245 cMpc for the EoR signal), introduce sky-model defects into skymodel construction, include SKA1-Low instrumental response with visibility synthesis, and reconstruct the radio sky with map-making; (ii) the analytic fraction of the HEVAL pipeline adopts the widely used logarithmic implementation of the calibration equations (Pearson & Readhead 1984; Wieringa 1992; Liu et al. 2010, hereafter the PWL method) to mimic the calibration processes of EoR experiments and estimates blending-induced calibration errors and propagation biases. By utilizing a pair of sky models, consisting of a perfect sky model and a corrupted one, HEVAL decouples the instrumental noise bias  \nfrom the calibration solution and derives the relative errors from poor calibration due to sky-model defects. The combination of end-to-end simulation and analytic analysis enables HEVAL to isolate the impact originating solely from source-blending defects and to quantify the blending-induced error budget under realistic SKA observations.  \nThis paper is organized as follows. Section 2 introduces the definition of source blending adopted for this study, identifies its underlying impact on sky-model construction, and evaluates the blending ratio expected for the SKA1-Low. Section 3 presents our HEVAL pipeline and the 2D PS measurement space for evaluating the impact of source-blending defects on the calibration of SKA EoR experiments. We introduce both the key components of our custom end-to-end simulation suite and the analytic formalism of blendinginduced residual gain error and its propagation under a sky-based calibration scheme. In Section 4, we present the propagated residual power in the 2D PS space and identify the blending tolerance of the SKA sky-model construction based on the residual power contamination. We then discuss the implications of these results for EoR experiments and provide potential strategies for mitigating blending-induced bias in Section 5. Throughout this work, we adopt a flat Λ CDM cosmology with 𝐻 0 = 100 ℎ = 67 . 66 km s -1 Mpc -1 , Ω 𝑚 = 0 . 3096, Ω Λ = 1 -Ω 𝑚 = 0 . 6904, Ω 𝑏 = 0 . 0489, 𝑛 𝑠 = 0 . 9665, and 𝜎 8 = 0 . 8102.", '2 SOURCE BLENDING': "Blending of astronomical sources often refers to the overlapping of structural components in the projected sky. To properly evaluate and simulate the source-blending effect, a clear definition of blending is required. In this study, source blending is defined as independent, unrelated objects (including compact sources, elements of multicomponent sources, and side lobes) distributed close enough in the projected sky plane so that they cannot be accurately measured independently under the telescope's resolving condition. Based on the spatial relationship between the two objects in the image space and the existing blending definitions (Arcelin et al. 2021; Sanchez et al. 2021), we identify three different blending cases:  \n- · Stacked : The projected radio emission contours of two objects coincide entirely, such that 100 per cent of the flux density distribution overlays each other.\n- · Overlapping : The projected radio emission contours of two objects partially intersect, with greater than 0 per cent but less than 100 per cent of the flux density distribution of one object overlapping the other.\n- · Proximity : The projected radio emission contours of two objects are close to each other but have a 0 per cent overlap in their flux density distribution. The angular separation of the peak flux density centroids can be up to several times the telescope's resolving beam.  \nIn the most common form, source blending in the radio domain occurs when two individual objects blend into an apparently elongated or reshaped source. For simplicity, we only consider the blending of two objects in our study. Although we note the existence of multisource blends, multi-source blending cases can be considered as the superposition of two-object blends. Using the observation sensitivity estimated by Prandoni & Seymour (2015) and the simulated source model from the SKA Designed Study Simulated Skies ( 𝑆 39 ; Wilman et al. 2008, hereafter W08), we can roughly estimate the fraction  \nof sources that are blended with other sources for upcoming SKA1Low observations. By calculating the per unit-area density of all the sources ( 𝐷 , arcsec -2 ) above the SKA1-Low's detection limit and projected extension of each source ( 𝑎 , arcsec 2 ), the likelihood of blending with another source can be inferred as 𝐿 = 𝑚𝑖𝑛 ( 𝑎 / 𝐷 / 𝑇, 1 ) , where 𝑇 ( 𝑇 ∈ (0, 1]) is the dimensionless threshold of blending. Subsequently, a quick and dirty estimation of the SKA blending ratio can be achieved by calculating the blending likelihood for each source above the SKA1-Low's detection limit. The total blending ratio of all the sources for the SKA1-Low is estimated to be approximately 5 -28 per cent considering different thresholds 𝑇 and projected extensions 𝑎 . As expected, this level of blending falls shy of the blending ratio of those optical surveys (see Melchior et al. 2021, and references therein), but still poses a significant challenge for SKA sky-model construction. In Fig. 1, we present the cumulative distribution function of the blending ratio under SKA1-Low's observation specifications and conditions at 150 MHz with a 𝑇 = 0 . 68 threshold. Given that the underlying W08 source-count models are characterized by Schechter parametrization, blending sources are unsurprisingly dominated by sources at the faint end.  \nFor the sky-model construction, defects due to source blending are introduced during the detection and measurement phases. In a typical catalogue-building pipeline, source-detection methods decompose radio emission into Gaussian components and group the associated components as detected sources (Mohan & Rafferty 2015; Hancock et al. 2018). The presence of source blending often leads to biased detection because the source finder may struggle to decompose the correct Gaussian components owing to blended structures. In addition, component associations can also lead to source-blending defects due to the complex morphology of EDRS populations (e.g. see the G4Jy Sample, White et al. 2020b,a). Without physical constraints from the host galaxy (links to optical counterparts through crossmatch techniques, such as Line et al. 2017), distinctively separated components, such as isolated lobes and nearby faint compact sources, are easily associated as one blended source. Therefore, source blending significantly affects the accuracy of sky-model construction and can cause uncertainties with a largely two-fold impact: (i) the spatial domain flux density and position error 10 , which originate from a mixture of originally independent flux density distributions, and (ii) the frequency-domain spectral index deviation , which is caused by the blending of intrinsically different spectral components. As most limiting factors of sky models (e.g. instrumental noise and defects from sky-model incompleteness) are often coupled with one another, methods with the ability to isolate the error contribution of one error of origin at a time are required to evaluate the impact of each sky-model defect properly. Section 3 introduces one such method to evaluate the sole impact of calibration against sky models with source-blending defects for SKA EoR experiments.", '3 METHODOLOGY': "This section presents a hybrid approach for evaluating the calibration impact of sky-model defects originating from source blending for future SKA EoR experiments. Fig. 2 outlines the data and software blocks of our HEVAL pipeline, in which the solid blue and dashed pink lines mark the end-to-end and analytic data flows, respectively.  \nFigure 1. The cumulative distribution function for blending ratio of extragalactic discrete objects under the SKA1-Low observation specifications and conditions. The x-axis represents the 150 MHz flux density of the discrete radio sources of the SKA1-Low observing sky, plotted on a logarithmic scale in the unit of Jy, whereas the y-axis represents the cumulative probability. The distribution is calculated by dividing the cumulative number of blended sources for each flux-density bin by the total number of blended sources. The shadow marks the 3-sigma uncertainty estimated using a Jackknife resampling technique with 60 subsamples. The lower and upper flux density are limited by the SKA1-Low's detection limit and the W08 simulation, respectively. For this estimation, the SKA observation condition is set using the synthesized beam size (8 arcsec) of the SKA1-Low at 150 MHz; the threshold of blending is set as 0.68; the total blending ratio of all the sources is estimated to be 7.25 per cent.  \n<!-- image -->  \nWe introduce the complete set of HEVAL modules 11 , including sky-model building, visibility synthesis, calibration error estimation, and map-making, and eventually arrive at the measurement space to quantify the blending-induced impact through 2D PS analysis. To better put our method and results into context, we attempt to summarize the assumptions adopted by the HEVAL pipeline in Table 1. We list the data-related mathematical notions in Table 2 and present the space transformation of the data cubes in Fig. 8.  \nFirstly, the HEVAL pipeline starts with the simulation in the image space by presenting the construction of sky models along with the introduction of sky-model defects, resulting in a pair of sky models, one with source-blending defects and one without, and the simulation of the low-frequency radio sky, which includes extragalactic discrete foregrounds, Galactic diffuse foregrounds, and 21-cm signals. All maps are simulated across frequency channels forming sky-map cubes [ 𝑆 ori ( 𝑙, 𝑚, Δ 𝑓 ) ] along the frequency axis. The modules used are from our custom F/g.pc21S/i.pc/m.pc+ 12 simulation suite, which is developed based on our previous simulation efforts (Wang et al. 2010, 2013; Li et al. 2019), notably the /f.pc/g.pc21/s.pc/i.pc/m.pc 13 simulator. As our next-generation simulation suite, F/g.pc21S/i.pc/m.pc+ is designed to simulate the low-frequency radio sky in the image space down to the arcsec level to cover the smallest angular scales resolved by high spatial resolution instruments, such as LOFAR and SKA1-Low. In addition, F/g.pc21S/i.pc/m.pc+ can introduce imperfections into the simulation, making it particularly useful for evaluating sky-model defects within data analysis pipelines.  \nSecondly, we detail image-to-visibility space transitions within the  \nHEVAL pipeline. Cross-correlation synthesis with the SKA1-Low array configuration is achieved by feeding the F/g.pc21S/i.pc/m.pc+-generated sky-map cubes of the sky-model pairs, foreground emissions, and the EoRsignals to the OSKAR array simulator, where the output visibility cubes mimic realistic instrumental responses under the SKA1-Low observation conditions. The array configuration and corresponding SKA1-Low specifications are summarized in Tables 3 and 4, respectively. Within the visibility space, we estimate the propagation effects during interferometric calibration under the sky-based calibration scheme by presenting an analytic analysis to infer the relative residual gain errors from poor calibration owing to sky-model defects and their subsequent propagation visibility biases. We use visibilities of the simulated sky model pairs to evaluate the blending-induced per-frequency per-antenna relative residual gain error and, subsequently, apply the relative error solutions to the visibilities of each sky component to infer the per-frequency per-baseline propagation visibility biases. Although no explicit foreground-removal module is implemented by HEVAL, the derived visibility biases are similar to the residual visibilities generated by pipelines explicitly performing foreground removal, owing to the nature of the relative calibration. Detailed analytic derivations are introduced in Section 3.4. In the visibility space, HEVAL produces two sets of visibility cubes: (i) the sky emission visibility cubes [ 𝑉 ori ( 𝑢, 𝑣, Δ 𝑓 ) ] representing the true response of the extragalactic discrete foregrounds, Galactic diffuse foregrounds, and EoR signals and (ii) the blending-induced visibility-bias cubes [ 𝑉 res ( 𝑢, 𝑣, Δ 𝑓 ) ] for each sky component.  \nThirdly, the HEVAL pipeline converts the visibilities back to the image space to reconstruct the radio sky through the map-making processes presented in Section 3.5, given our aim of inferring the residual PS impact under the 'reconstructed'-sky approach (Morales et al. 2019). The WSC/l.pc/e.pc/a.pc/n.pc imager is used for the Fourier transform, convolution, and de-convolution of the simulated visibilities. Both sets of visibility cubes [ 𝑉 ori ( 𝑢, 𝑣, Δ 𝑓 ) and 𝑉 res ( 𝑢, 𝑣, Δ 𝑓 ) ] are converted  \nFigure 2. The HEVAL pipeline for sky-model defect evaluation. The solid-blue lines mark the end-to-end simulation faction for sky-map generation, SKA1-Low visibility realization, and map making. The dashed pink lines mark the analytic analysis approach for the gain error estimation and visibility bias calculation. The solid blue lines mark the end-to-end simulation faction, including sky-map generation, SKA1-Low visibility realization, and map making, whereas the dashed pink lines mark the analytic analysis approach for the gain error estimation and visibility bias calculation. Although typical foreground removal is not explicitly included in the pipeline, this process is implicitly handled through our analytical fraction of the hybrid HEVAL pipeline directly in the visibility space, and the analytic derivation of visibility biases is analogous to the visibility residuals produced by pipelines subtracting the foregrounds. Coloured double diamonds and rounded grey squares represent data blocks and software blocks, respectively. From the data flow perspective, the upper and lower fractions denote the ill-calibrated residual and sky emission data flow, respectively. The simulation of sky model pairs, using ES/i.pc/m.pc (labelled 1) and B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc (labelled 2), is introduced in Section 3.1; the generation of evaluation data, using ES/i.pc/m.pc (labelled 1), GS/i.pc/m.pc (labelled 3), and 21/c.pc/m.pcS/i.pc/m.pc (labelled 4), are detailed in Section 3.2; the visibility realization, using OSKAR (labelled 5), is described in Section 3.3; the analytic analysis of the source-blending impact within the visibility space is introduced in Section 3.4; the step-by-step analytic derivation of gain error (labelled 6) and the analytic calculation of the propagation visibility bias (labelled 7) are detailed in Appendix B1 and Appendix B3, respectively; and the convolution and deconvolution of the radio maps, using WSC/l.pc/e.pc/a.pc/n.pc (labelled 8), are described in Section 3.5; assumptions adopted by the HEVAL pipeline are organized in Table 1; data-related mathematical notions are noted in Table 2; the array configuration and the simulation specifications used for this study are listed in Tables 3 and 4, respectively.  \n<!-- image -->  \ninto two sets of image cubes [ 𝐼 ori ( 𝑙, 𝑚, Δ 𝑓 ) and 𝐼 res ( 𝑙, 𝑚, Δ 𝑓 ) ] for each sky component, respectively.  \nFinally, we establish the PS measurement space for statistical EoR detection and transform the two sets of image cubes into two sets of three-dimensional (3D) PS [ 𝑃 ori ( 𝑘 𝑥 , 𝑘 𝑦 , 𝑘 𝑧 ) and 𝑃 res ( 𝑘 𝑥 , 𝑘 𝑦 , 𝑘 𝑧 ) ] and calculate the cylindrical-averaged 2D PS [ 𝑃 ori ( 𝑘 ⊥ , 𝑘 ∥ ) and 𝑃 res ( 𝑘 ⊥ , 𝑘 ∥ ) ] for each sky component, respectively. Being a key EoR figure of merit, 2D PS, along with EoR windows, is used to evaluate the eventual impact of blending-induced poor calibration.  \nThroughout this paper, we denote the F/g.pc21S/i.pc/m.pc+-generated sky emission, OSKAR-simulated cross-correlations, WSC/l.pc/e.pc/a.pc/n.pcconvoluted, and -deconvoluted images as sky maps, 'observed' visibilities, dirty maps, and clean maps, respectively. To clearly distinguish between the terms, we label the terms of sky emission with 'original' (e.g. 'original' image cubes and 'original' visibility cubes). In contrast, we label all the poor-calibration related terms with 'residual' (e.g. 'residual' visibility-bias cubes and 'residual' powers).", '3.1 Simulation of sky models and source-blending effect': "Wedetail the realization of sky-model construction in the image space. This study follows the convention and includes only EDRS populations in the sky model using the ES/i.pc/m.pc module, because most commonly used sky models for interferometric calibration are constructed using  \nEDRSpopulations, often from high-quality EDRS catalogues (Carroll et al. 2016; Procopio et al. 2017; White et al. 2020b,a; Kondapally et al. 2021) and dedicated surveys (Norris et al. 2011; Lane et al. 2014; Wayth et al. 2015; Heald et al. 2015; Jackson et al. 2016; Zheng et al. 2016; Lacy et al. 2016; Shimwell et al. 2017; Intema et al. 2017; Hurley-Walker et al. 2022). While diffuse emission contributes power to the radio sky at larger angular scales, the limited existing data and algorithmic challenges leave discrete sources as the primary ingredients for modelling the sky (Liu & Shaw 2020).  \nOur next-generation F/g.pc21S/i.pc/m.pc+ suite offers the flexibility to add defects to sky models in the image space through defect-simulation modules, such as B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc, which introduces additional sourceblending effects to the ES/i.pc/m.pc-simulated sky models. By combining ES/i.pc/m.pc and B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc, our simulation of sky-model construction and blending-induced defects arrives at a pair of sky models: (i) an ' ideal ' sky model representing the true flux density distribution of discrete sources without any blending of sources and (ii) a ' blended ' counterpart with some of the sources blended together resulting in flux density distribution deviations with respect to (w.r.t) their ' ideal ' counterparts. The data flow is reflected in the upper-left part of Fig. 2. With the ' ideal ' and ' blended ' sky models, our HEVAL pipeline can utilize a paired sky model estimation approach to derive the sole contribution of source blending during poor calibration, which is achieved by isolating the impact of blending defects from any  \nTable 1. Summary of assumptions adopted by the HEVAL pipeline.  \nTable 2. Data-related mathematical notions used in the paper.  \nother limiting factors of sky-model construction (e.g. baseline noises and other sky-model defects) and inferring the blending-induced relative residual gain errors. Given that the relative residual gain error and propagation bias are subject to the sky model differences in the visibility space, the simulation of sky models with end-to-end control and clearly defined procedures is crucial for estimating the source-blending effect properly.", '3.1.1 Simulation of EDRS populations': "EDRSpopulations (see Padovani 2016; Panessa et al. 2019; Blandford et al. 2019; O'Dea & Saikia 2021, for reviews) are key objects of interest in the radio window, not only due to their diversity in terms of morphological types but also due to the underlying physical processes, black-hole accretion and star formation, which govern how galaxies evolve. As the bulk of the extragalactic sky, EDRSs dominate the radio sky in terms of both number and brightness, making them a key foreground component of EoR experiments. Consequently, EDRS populations are critical for calibration, imaging, and subtraction in the field of 21-cm science. In light of this, we compose our sky model simulation by considering only extragalactic radio sources with radio emission originating from supermassive black hole (SMBH) activity, which comprises compact radio core emission and structural features (such as relativistic jets, plumes, lobes, and hotspots), and star  \nformation processes, which produce synchrotron radiation through relativistic plasma within the supernova remnants, associated with galaxies.  \nThe EDRS populations included in the ES/i.pc/m.pc module consist of (i) bright extended radio galaxies (RGs), which are Fanaroff-Riley (FR) type I and type II galaxies (Fanaroff & Riley 1974), (ii) compact radio-loud (RL) AGNs, which consist of both compact steep-spectrum (CSS) and peaked-spectrum (PS) AGNs, (iii) radio-quiet (RQ) AGNs, and (iv) star-forming (SF) and starburst (SB) galaxies. The underlying source count distribution follows the commonly used W08 source catalogue simulation. To generate sky maps of individual sources, we simulated the brightness distribution of compact sources (i.e. compact AGNs, SF, and SB galaxies) with a single 2D Gaussian model and extended sources (i.e. FR-I and FR-II RGs) with multipleGaussian models by adopting morphological properties from the W08 catalogue (e.g. angular size and position angle). In particular, the extended sources are created from the convolution of a 2D Gaussian kernel with a radial profile along the jet axis, which results in a superposition of multiple Gaussian distributions, representing the extended structures, such as the jet and lobes of radio galaxies (Fig. 3). For the spectral models, ES/i.pc/m.pc utilizes a power-law spectrum for most radio spectra, whereas a curvature model from Jarvis & Rawlings (2000) and spectral turnover discovered by O'Dea (1998) are implemented for compact core emission and PS sources, respectively.  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Illustration of ES/i.pc/m.pc-generated radio galaxies using a multiple-Gaussian model. (Upper) ES/i.pc/m.pc-simulated flux-density distributions of typical (a) FR-I and (b) FR-II sources. The two flux density distributions are presented using the same colour bar. The dashed white line marks the jet axis. (Lower) The corresponding radial profiles along the jet axis of the two sources above. The blue line marks the radial profile of the extended structures, whereas the red line marks the overall radial profile of the whole source. The two figures from the same column are plotted with the matching pixel coordinates.  \n<!-- image -->", '3.1.2 Simulation of the blending effect': "The B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc module is used to introduce source-blending effects in both the spatial (Fig. 4) and spectral domains (Fig. 5). The implementation of B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc contains all the EDRS populations in Section 3.1.1, resulting in 15 different blending types, such as RG-RQ AGN, AGN-SF, and RG-SB. To mimic the spatial displacement of brightness distribution, a beam with a limited resolution is used to convolute the ideal brightness distribution of the selected source pairs in the 'stacked', 'overlapping', or 'proximity' case by adjusting the distance between the centres of two objects. The width of the convolution beam should not be larger than that of the SKA synthesized beam, and we set the beam size to 6 arcsec for this study. To consider only spatial distribution deviations rather than flux density errors, B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc employs a blending convolution procedure that utilizes a normalized kernel to retain the same total flux density of the source pair before and after blending.  \nAfter the source pair is blended in the spatial domain, B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc adds the blending-induced error in the spectral domain. Because our sky-model building does not include spectral-index spatial variations for individual sources, the spectral domain error is introduced via the spectral index deviations of the integrated spectral indices of the blended pair. Given that most of the blended sources present either  \npower-law or power-law-like spectra within our simulated frequency band, a power-law 14 is used to model the integrated spectral index of the blended spectra (e.g. the blended spectra presented in Fig. 5). After the spectral model is assigned, B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc generates multi-frequency sky models using the ' blended ' monochromic sky model with the estimated blended spectral indices.  \nFor the simulation of blending defects, it is important to note that B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc does not introduce flux density errors, which are most likely unavoidable for blending sources in real data considering image noise and artefacts. Flux density errors will not only introduce their own impact on calibrations but also bias spectral modelling, leading to additional spectral index errors. The coupling of different error sources hinders the final interpretation of the impact in the measurement space. Given that our aim is to decouple the sole impact of blending, we explicitly choose not to simulate any flux density errors to ensure that the evaluation and interpretation target only the source-blending impact. While this omission of flux density errors in the simulation does mean that the errors of blending are most likely  \n14 Unlike spectral fitting with real data, our modelling processes are ideal cases with no noise or errors that may cause spectral fluctuations or distortions. A power-law model approach to mimic spectral index deviation is sufficient for this study.  \nFigure 4. Mostcommonlyoccurredblendingtypesbetweentwocomponentsin the image space (the cases of two extended structures fully stacked are omitted due to the rare occurrence). The rule-of-thumb illustrations demonstrate the possible blending-induced error in the image space during source detection with the presence of noises. The elongated and round shapes indicate extended and compact features, respectively. The colour-filled shapes represent the ground truth of the components, whereas the dashed black lines mark the detected components by source finders. For bright sources, the detection errors are dominated by blending. For faint sources, image noises contribute to flux density distribution bias in addition to blending-induced uncertainties.  \n<!-- image -->  \nunderestimated compared to real-world scenarios, it allows us to directly attribute the propagated biases to blending defects alone. This approach provides a clean and controlled evaluation, isolating the effects of a single source of origin, without the confounding influence of other error sources.", "3.1.3 The 'ideal' and 'blended' sky model": "To generate the sky map of the ' ideal ' sky model and its ' blended ' counterpart, a basis sample of 4,000 sources, ranging from 0.1 to 1000 mJy, are sampled from the W08 catalogue satisfying its underlying source count model with a minimum separation requirement of 100 arcsec within the estimated SKA1-Low FoV of 5 · . Therefore, the sparsely distributed sources are not biased by blending effects. With the sky model basis at hand, a sub-sample of the basis is selected as the host of the blended source pair and assigned to another source with a maximum separation of 6 arcsec, corresponding to the synthesized beamsize of the SKA1-Low. Using the ideal flux density distribution and the true spectral index of both the source pair and the basis sample, sky-map cubes of the ' ideal ' multi-frequency sky models are generated. As for the counterpart, the blended flux density distribution and deviated spectral index are simulated using B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc. By replacing the ideal model of the source pair with the blended model, we obtain sky-map cubes of the ' blended ' multifrequency sky models. In Fig. 6, we present a small patch of a pair of ' ideal ' and ' blended ' sky models. Adopting the 15 types of blending pairs described in Section 3.1.2, B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc assigns the types based on the likelihood using the W08 number count model.  \n(a) Compact-Compact Blend  \nFrequency [Hz]  \n<!-- image -->  \n(b) Compact-Extended Blend  \nFigure 5. Demonstration of blending defects in the spectral domain. (Top) Blended spectra of the two compact radio sources. (Bottom) Blended spectra of compact and extended sources. The two dashed lines mark the original spectra of the individual sources and the purple line marks the spectra of the blended pair. The green shaded area indicates the frequency range used in this study. In both cases, an apparent spectral index deviation exists after the two individual sources are blended.  \n<!-- image -->  \nFor a systematic estimation of the source-blending impact on the SKA EoR experiments, a set of 2, 20, and 200 blended source pairs are generated for the basis sample of 4,000 sources, representing source-blending levels of 0.05 per cent, 0.5 per cent, and 5 per cent, respectively. Because we aim to find an expectedly modest blending tolerance of sky-model construction for SKA EoR experiments, the 5 per cent minimum of the estimated SKA blending ratio is used as the maximum of the three levels. In the latter part of the manuscript, we refer to the three blending levels as 'mild', 'moderate', and 'high' blending ratios, respectively. For simplicity, sky models affected by the three blending ratio levels are marked as 'mildly'-, 'moderately'-, and 'highly'-corrupted sky models, respectively.", '3.2 Simulation of the evaluation data': "Evaluation data cubes are required to represent the low-frequency radio sky, including the EoR signal and its foregrounds, to evaluate the propagation effect of residual gain errors. The SKA1-Low has baselines across an extensive range of lengths that can probe the EoR signal across approximately three decades in scale. The ability to detect the EoR signal at a substantial number of scales, both large and small, gives the edge of the SKA1-Low over other EoR-aimed  \nFigure 6. Demonstration of a monochromic sky model pair simulated using our ES/i.pc/m.pc and B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc models. (Left) A small patch (1 . 5 arcmin × 1 . 5 arcmin) of an ' ideal ' sky model with the ideal source model pair. (Right) The same sky patch of the corresponding ' blended ' sky model with the blended source model pair. For a sky model pair, the differences between the ' blended ' and ' ideal ' sky models are the flux density distributions of these source model pairs marked by the dashed white lines, whereas the rest of the sources have exactly the same flux density distributions. Since B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc normalizes the blending convolution procedure, no flux density uncertainties are introduced for source model pairs. The two figures are plotted using the same colour bar. Multi-frequency sky model pairs are simulated by supplying the monochromatic version with spectral models. For this particular pair, the spectral models were similar to those presented in Fig. 5 (b).  \n<!-- image -->  \ninstruments. Even though the total power of the 21-cm signal is expected to be constituted mainly from large-scale contributions, the EoR signal at the medium and small scales pokes unique aspects of the Universe, particularly during the later stages of the EoR. For one thing, the 21-cm signal can trace the small-scale structure and its coupling to the large-scale structure (Lidz et al. 2007). For another, the EoR signal is expected to be highly non-Gaussian during the later stages of reionization (Mondal et al. 2016; Majumdar et al. 2018), and its non-Gaussianity is probed at the smaller spatial scales. Therefore, to thoroughly investigate the foreground coupling effect of the blending-induced calibration error for SKA EoR experiments, it is crucial to simulate the radio sky containing physical scales that are spatially resolvable by the SKA1-Low baselines. In the remainder of this section, we introduce the F/g.pc21S/i.pc/m.pc+ modules used to generate evaluation data containing simulated scales from the degree level down to the smallest scales resolvable by the SKA1-Low. The data flow of the low-frequency radio sky simulation is illustrated in the lower left part of Fig. 2. The corresponding F/g.pc21S/i.pc/m.pc+-generated sky maps are presented in Fig. 7.", '3.2.1 Extragalactic foregrounds': 'The extragalactic foregrounds, dominated by EDRS populations (see Section 3.1.1), predominate the contamination of the EoR signal mainly on smaller scales because they are discrete sources with flux density fluctuations at the arcsec to arcmin scales. Although other extragalactic sources, such as diffuse radio emission associated with galaxy clusters, present unique challenges to the EoR experiment and impose contamination at the arcmin scale (Li et al. 2019; Zhou et al. 2022), they are ultimately excluded from our model because of both the relatively faint flux, which is three to four orders of magnitude lower, and a significantly lower number density compared to the EDRS populations. We use the ES/i.pc/m.pc module to generate sky maps of the extragalactic foregrounds, adopting the simulation specifications from Section 3.3.2. Here, the full EDRS source-count model of the  \nW08 catalogue is implemented, in contrast to the sampling procedure when simulating sky models detailed in Section 3.1.', '3.2.2 Galactic foregrounds': 'The Galactic foregrounds mainly contain the diffuse Galactic synchrotron and free-free emission, with the synchrotron component being the dominant among the two. Our GS/i.pc/m.pc treats both the synchrotron and free-free components separately and simulates a patch of the sky using basis template maps: (i) For Galactic synchrotron emission, the basis synchrotron template map deployed by GS/i.pc/m.pc is the reprocessed Haslam 408 MHz all-sky map 15 , which is referred to as the HAS14 map, by Remazeilles et al. (2015) with significantly improved removal of both extragalactic sources and instrument artefacts. To calculate the sky maps in the SKA1-Low frequency via extrapolation with a power-law model from the 408 MHz, the all-sky spectral index map ( 𝛼 ∈ [2.50, 3.20], Giardino et al. 2002) of the Galactic synchrotron emission is used. (ii) As for the Galactic freefree emission, the basis template map is simulated through indirect estimations, as direct observations of the free-free component, which is overwhelmed by the strong synchrotron counterpart, is not possible. Owing to the common origin of free-free and H 𝛼 emission, a tight relation between the two has been demonstrated. It can be utilized to calculate the free-free emission (see Dickinson et al. 2003, and references therein). The GS/i.pc/m.pc module largely follows the prescription of Dickinson et al. (2003), employs the H 𝛼 survey from Finkbeiner (2003), corrects the dust absorption with a full-sky 100µ mdust map (Schlegel et al. 1998), and subsequently, simulates the free-free emission. Because the free-free and H 𝛼 relation is frequency dependent, the conversion from the H 𝛼 map to the free-free map is calculated for the desired SKA1-Low frequency window.', '(a) Galactic Emission (Synchrotron)': 'Figure 7. Simulated sky maps of evaluation data representing the low-frequency radio sky: the Galactic diffuse (a) synchrotron and (b) free-free emission, the (c) extragalactic discrete emission, and the (d) EoR signal at 158 MHz. All the images are zoomed in to the central 1 · × 1 · to show the small-scale structures. Each figure uses its own colour bar for the display. In particular, the EoR signal is plotted using unit [mK] different from the rest.  \n<!-- image -->  \n[K]  \n[mK]  \nLimited by existing observations, Galactic contributions to the radio sky are primarily at large scales. Under the scope of this work, flux density fluctuations at small scales are also required to properly estimate the impact at the longer baselines of the SKA1-Low. Our implementation in GS/i.pc/m.pc also includes the addition of Galactic small-scale fluctuations using the realization of a Gaussian random field (GRF). Our implementation of adding small-scale fluctuations is similar to the existing efforts based on the realization of GRFs (Miville-Deschênes et al. 2007; Delabrouille et al. 2013; Remazeilles et al. 2015; Hervías-Caimapo et al. 2016; Thorne et al. 2017) but with some enhancements. The fundamental idea is to extrapolate the angular PS ( C ℓ , parametrized by 𝛾 ) to the required small scales ( ℓ ) and generate the small-scale fluctuations ( 𝐺 ss , parametrized by 𝛼 and 𝛽 ) using a GRF corresponding to the extrapolated C ℓ . We detail the step-by-step realization of adding small-scale fluctuations in Appendix A.  \nTo generate sky maps of the galactic foregrounds, we set the simulation centre of sky maps at a high galactic latitude centre position (R.A., Dec. = 0 · , -27 · ) since future SKA observations are expected to point at high galactic latitudes (e.g. | 𝑏 | > 60 · ) to minimize contamination. As for the small-scale fluctuations, we use the best-fitting of the three required parameters: (i) the 𝛾 index of synchrotron and free-free component is fitted as -2.220 and -2.426, respectively, using the full sky template maps; (ii) the 𝛼 and 𝛽 index of the small-scale maps are fitted for the synchrotron ( 𝛼 = 0 . 0342,  \n𝛽 = 0 . 227) and free-free emission ( 𝛼 = 0 . 00785, 𝛽 = 0 . 526) using a patch of the template map centred at the selected simulation centre.', '3.2.3 EoR signal': "The final component of the radio sky is the 21-cm signal. To compensate for the need for high spatial resolution and a relatively large FoV, a fast semi-numerical approach is used in the 21/c.pc/m.pcS/i.pc/m.pc module by adopting the 21/c.pc/m.pcFAST 16 (Mesinger et al. 2011; Murray et al. 2020) package. The cosmic reionization process is simulated using the same physical parameters (A. Mesinger, priv. comm.) as the 'faint galaxies' model of the Evolution Of 21 cm Structure project 17 (Mesinger et al. 2016) with a box of 294 comoving Mpc and 1200 cells along the side. Our simulation run resolves the smallest physical scale at 0.245 cMpc, making the finest-resolution 21-cm simulation to date using 21/c.pc/m.pcFAST. From the simulated light-cone object, the designated sky maps of each frequency band are sliced w.r.t the corresponding redshift and tiled according to the simulation specifications listed in Table 4 to form sky-map cubes.  \n- 16 21/c.pc/m.pcFAST: https://github.com/21cmfast/21cmFAST\n- 17 Evolution Of 21 cm Structure project: http://homepage.sns.it/ mesinger/EOS.html  \nTable 3. Array specifications adopted for the SKA1-Low.", '3.3 The realization of simulated SKA observations': "This section presents the transformation from the image space to the visibility space. First, we detail the SKA array configuration used and the corresponding technical specifications of the given array. Then, we describe the simulation specifications under limited computational powers. Finally, we introduce the realization of visibility synthesis using an array simulator which outputs 'observed' visibilities with instrumental responses of the SKA1-Low. In the end, we transform all the sky-map cubes into visibility cubes for both the sky model pairs and the sky components, as presented in the middle-left part of Fig. 2.", '3.3.1 Array configuration': "As the most advanced and transformational project in radio astronomy to date, the upcoming SKA will be constructed in two phases, owing to challenges and difficulties in designing, constructing, and operating extremely large telescopes. Among the two telescopes of Phase 1, SKA1-Low in Australia and SKA1-Mid in South Africa, the low-frequency aperture array SKA1-Low operating in the 50 MHz to 350 MHz frequency range is the instrument intended for EoR detection and 21-cm cosmology.  \nBoth the SKA1-Low and SKA1-Mid instruments are yet to be fully built, although pathfinders, precursors, and verification systems have been built and operated as stepping stones for the operation of SKA1. For this study, we use the array configuration from the SKA1 System Baseline Design (Dewdney et al. 2013, 2016, 2019). These reports detail that SKA1-Low will contain 131072 log-period dipole antennas with band coverage from 50 MHz to 350 MHz, assembling 512 stations of diameter 35 -40 m each consisting of 256 antennas. Among all the stations, 224 stations are placed randomly within the 'core' region of 1000 m in diameter, and 288 stations are distributed into 'clusters' along the three array spiral arms, which extend the baselines up to 65 km. Following the SKA1-Low Design Baseline, the estimated SKA1-Low specifications are summarized in Table 3.", '3.3.2 Simulation specification estimation': 'The high spatial resolving power of the SKA1-Low advocates the capability of SKA EoR experiments to probe the structural scale ranging from several arcsecs to degrees via 21-cm signals. In accordance with the estimated spatial resolution of the SKA1-Low ( ∼ 6 arcsec at 196 MHz), the simulated sky maps are required to contain the smallest scale at the 6 arcsec level ( ∝ λ ) at the corresponding frequency, which expects the pixel size of the sky map to be approximately 2 arcsec.  \nWithin the designed SKA1-Low frequency window, three sparsely selected frequency bands, with a bandwidth of 8 MHz are chosen to limit the possible cosmological evolution of the EoR signal (e.g. Thyagarajan et al. 2013; Li et al. 2019), covering 120 - 128 MHz, 154 - 162 MHz, and 188 - 196 MHz.  \nSimulating the radio sky with both high spatial and frequency resolution across extensive sky coverage presents significant computational challenges. This is because each stage of the process requires substantial resources: (i) the cosmological simulation of the 21-cm signal is memory-intensive, (ii) the generation of multi-frequency sky maps demands a high-performance GPU, and (iii) the visibility synthesis and subsequent map-making processes require substantial CPU power. Consequently, the combined high demands on CPU, GPU, and memory inherently restrict the practical specifications of our simulations. Given these limitations, it is impractical to simulate the SKA observing sky in a way that achieves a large field of view (FoV) and high spatial resolution simultaneously. Since spatial resolution plays a more critical role under the scope of this study (e.g. EDRS-only sky models and determination of blending of sources), we prioritize spatial resolution over areal coverage. However, to still adequately cover the key EoR scales (0 . 1 ≲ 𝑘 ≲ 2 Mpc -1 ), we chose a reasonable balance to simulate sky patches of up to 2 · × 2 · . The simulation specifications for both sky models and evaluation data are presented in Table 4.', '3.3.3 Array simulator': "The realization of visibilities with practical instrumental effects of the SKA1-Low is archived with the OSKAR 18 array simulator (Mort et al. 2010) by generating SKA1-Low 'observed' visibilities with the SKA1-Low specification and configuration using GPU accelerations. The specific input telescope model is created using the SKA1-Low array layout with antenna coordinates 19 from Dewdney & Braun (2016). The three frequency bands (Section 3.3.2) are divided into 51 channels with a channel width of 160 kHz. To minimize the pointing effect, we choose to centre all the sky maps at the sky position of (R.A., Dec.) = (0 · , -27 · ), which is the expected SKA1-Low zenith.  \nSky-map cubes of both the ' ideal ' and ' blended ' sky models are directly passed to the OSKAR simulator without pre-processing processes under a sky coverage of 5 · × 5 · and a pixel size of 1 arcsec with a single 2-min snapshot exposure time. The sky-map cubes of the evaluation data, extragalactic discrete emission, Galactic diffuse emission, and EoR signal, are simulated with a sky coverage of 2 · × 2 · and a pixel size of 2 arcsec with deep observations for 6 h. For the EoR foregrounds, pre-processing progresses are made before the visibility synthesis. Being extended in nature, the Galactic components, synchrotron and free-free, are combined and tapered, which reduces the window effect of sky maps containing large-scale structures with a limited size before they are 'observed' by the simulator. Since the brightest sources of EDRS are already known to strongly bias the EoR observations (e.g. Datta et al. 2010), we assume these sources are adequately dealt with and removed during data analysis pipelines. Hence, our simulations mask the EDRS populations with a 158 MHz flux density above 50 mJy, as indicated  \nhttps://astronomers.skatelescope.org/wp-content/uploads/ 2016/09/SKA-TEL-SKO-0000422\\_02\\_SKA1\\_  \nLowConfigurationCoordinates-1.pdf (released on 2016 May 31)  \nTable 4. Simulation specifications of the sky models and evaluation datasets.  \nin Fig. 2. All the specifications and pre- and post-processing processes are listed in Table 4.", '3.4 Evaluation of blending-induced calibration errors': "This section introduces the analytic fraction of the HEAVL pipeline, which is shown in the central part of Fig. 2. A sky-based per-frequency per-antenna calibration scheme is first introduced as the basic framework. We then present the derivation of the frequency-dependent per-antenna gain error caused by a blending-corrupted sky model w.r.t the ' ideal ' sky model counterpart through analytic analysis based on the PWL logarithmic method. Unlike traditional methods aimed at a direct antenna-based calibration solver to infer the solution of the complex gain factors (Thompson & D'Addario 1982), we use an alternative approach to estimate the blending-induced relative gain error using a pair of ' ideal ' and ' blended ' sky models. By adopting the PWL logarithmic implementation, the analytic expressions arrive at two simple sets of linear systems of equations, resulting in an overdetermined linear regression problem mathematically. To solve for the estimated gain error, we use a singular value decomposition (SVD)-based method. With the solved relative residual gain error, we can derive an analytic expression of the propagation visibility bias of the sky signal in the visibility space induced by an ill-calibrated instrument. For conciseness, we state the critical steps in this section and leave the step-by-step mathematical derivation in Appendix B. A full description of the calibration-related mathematical notations is presented in Table B1.", '3.4.1 Calibration in a sky-based scheme': 'Considering an array with 𝑁 antennas, the 𝑖 th antenna measures the unpolarized voltage 𝑣 𝑖 at the antenna end at a single frequency channel 𝑓 within a limited time span Δ 𝑡 , and can be estimated as  \n𝑣 𝑖 = 𝑔 𝑖 𝑠 𝑖 + 𝑛 𝑖 , (1)  \nwhere 𝑔 𝑖 is the antenna complex gain factor, 𝑛 𝑖 denotes the intrinsic antenna noise, and 𝑠 𝑖 represents the true sky signal. As interferometers measure cross-correlations, each baseline correlates measurements from a pair of antennas, resulting in a time-averaged visibility of the baseline 𝑩 𝑖 𝑗 in the form of  \n𝑉 𝑖 𝑗 ≡ GLYPH<10> 𝑣 𝑖 𝑣 𝑗 GLYPH<11>  \n≊ 𝑔 𝑖 𝑔 𝑗 GLYPH<10> 𝑠 𝑖 𝑠 𝑗 GLYPH<11> + 𝑛 𝑖 𝑗 .  \n∗ (2a) = 𝑔 ∗ 𝑖 𝑔 𝑗 GLYPH<10> 𝑠 ∗ 𝑖 𝑠 𝑗 GLYPH<11> + 𝑔 ∗ 𝑖 GLYPH<10> 𝑠 ∗ 𝑖 𝑛 𝑗 GLYPH<11> + 𝑔 𝑗 GLYPH<10> 𝑛 ∗ 𝑖 𝑠 𝑗 GLYPH<11> + GLYPH<10> 𝑛 ∗ 𝑖 𝑛 𝑗 GLYPH<11> (2b) ∗ ∗ (2c)  \nBy assuming that only sky radio signals are correlated, the last three terms of equation (2b) can be reduced to noise 𝑛 𝑖 𝑗 specific to the  \nbaseline 𝑩 𝑖 𝑗 . Since each antenna gain is a complex factor, it can be parametrized by its amplitude 𝜂 and phase 𝜙 as  \n𝑔 𝑖 ≡ 𝑒 𝜂 𝑖 + i 𝜙 𝑖 . (3)  \nIf we use the parametric form of each antenna complex gain, the baseline measurement question, equation (2c), takes a new form:  \n𝑉 𝑖 𝑗 = exp GLYPH<2> GLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> GLYPH<3> GLYPH<10> 𝑠 ∗ 𝑖 𝑠 𝑗 GLYPH<11> + 𝑛 𝑖 𝑗 . (4)  \nIdeally, the calibration process solves both the true complex gain factor 𝑔 𝑖 of each antenna and the true sky signal cross-correlation 𝑆 𝑖 𝑗 ≡ GLYPH<10> 𝑠 ∗ 𝑖 𝑠 𝑗 GLYPH<11> of each baseline simultaneously. An array of 𝑁 antennas measures 𝐶 2 𝑁 visibilities, whereas the variables to solve for are 𝑁 true antenna gain amplitude 𝜂 𝑖 items, 𝑁 true antenna phase 𝜙 𝑖 items, and 𝐶 2 𝑁 ≡ 𝑁 ( 𝑁 -1 )/ 2 true sky cross-correlations 𝑆 𝑖 𝑗 . Therefore, the measurements are insufficient to directly solve for 𝐺 𝑖 𝑗 ≡ 𝑔 ∗ 𝑖 𝑔 𝑗 and recover 𝑆 𝑖 𝑗 . However, by adopting a sky-based calibration scheme, which supplies a sky model to replace the unknown sky cross-correlation, the calibration of 𝐺 𝑖 𝑗 becomes an overdetermined problem that can be solved using a least-squares solution. The sky cross-correlations, 𝑆 𝑖 𝑗 , can then be recovered by applying the solution of the estimated gain values. Calibrating all frequencies separately can recover a set of sky cross-correlations within the observed band.', '3.4.2 Analytic analysis of the calibration error': "As stated in Section 1, sky-based calibration techniques are ultimately restricted by the supplied sky model. Consequently, the precision and fidelity of these models becomes critical. Source blending, along with various limiting factors that couple with one another [such as model completeness (e.g. Barry et al. 2016; Gehlot et al. 2021), polarized emission (e.g. Moore et al. 2017; Asad et al. 2018), and diffuse components (e.g. Lanman et al. 2022; Sims et al. 2022)], collectively impacts the accuracy of calibration processes and subsequently dictates the detection ability of the EoR signal. To individually evaluate the source-blending impact on the overall calibration process, we propose a theoretical estimation that confines blending defects as the only sky model limiting factor using a pair of sky models consisting of an ideal sky model and its blended counterpart. The sky model pairs include only discrete sources and their simulation realization is presented in Section 3.1. Under this paired sky model estimation approach, we can use the measured visibility equation (equation 4) to infer the per-frequency per-antenna gain error induced by a ' blended ' sky model w.r.t its ' ideal ' counterpart within the sky-based calibration scheme. We refer to this type of calibration error as the relative residual gain error.  \nTo begin with, by using an ' ideal ' sky model 𝑉 ideal to calibrate the 𝑁 antenna composed array, we have the measured visibility and true  \ngain amplitudes and phases ( 𝜂 and 𝜙 , respectively) for the baseline 𝑩 𝑖 𝑗 in the form of  \n𝑉 𝑖 𝑗 = exp GLYPH<2> GLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> GLYPH<3> 𝑉 ideal 𝑖 𝑗 + 𝑛 𝑖 𝑗 . (5)  \nThen, by supplying a ' blended ' sky model ' 𝑉 blend , the same baseline measurement, albeit ill-calibrated, will take the form:  \n𝑉 𝑖 𝑗 = exp GLYPH<2> GLYPH<0> ' 𝜂 𝑖 + ' 𝜂 𝑗 GLYPH<1> + i GLYPH<0> ' 𝜙 𝑗 -' 𝜙 𝑖 GLYPH<1> GLYPH<3> ' 𝑉 blend 𝑖 𝑗 + 𝑛 𝑖 𝑗 , (6)  \nwhere ' 𝜂 and ' 𝜙 are the gain factors a ' blended ' sky model solves. Subsequently, by combining equations (5) with (6) under the PWL logarithmic method (see the step-by-step derivation in Appendix B1), we can obtain two sets of linear equations with the form  \nΔ 𝜂 𝑖 + Δ 𝜂 𝑗 = 𝑅 𝑖 𝑗 , (7a)  \nΔ 𝜙 𝑗 -Δ 𝜙 𝑖 = 𝐼 𝑖 𝑗 , (7b)  \nwhere both the relative gain errors ( Δ 𝜂 𝑖 ≡ 𝜂 𝑖 -' 𝜂 𝑖 and Δ 𝜙 𝑖 ≡ 𝜙 𝑖 -' 𝜙 𝑖 , respectively) and the visibility difference in terms of amplitude and phase ( 𝑅 𝑖 𝑗 ≡ ln GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑉 blend 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -ln GLYPH<12> GLYPH<12> GLYPH<12> 𝑉 ideal 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> and 𝐼 𝑖 𝑗 ≡ arg GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑉 blend 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -arg GLYPH<12> GLYPH<12> GLYPH<12> 𝑉 ideal 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> , respectively) are substituted for simplicity. For each baseline, both 𝑅 𝑖 𝑗 and 𝐼 𝑖 𝑗 are only subject to the visibility difference between the pair of sky models. Because we supplied simulated sky model pairs with source blending as the only limiting factor, the error contribution of a ' blended ' sky model w.r.t an ' ideal ' counterpart can be inferred by solving the amplitude and phase differences separately with those two simple sets of equations. In this way, we decouple the source-blending impact on the calibration of EoR experiments from both the baseline noises and other limiting factors of sky-model constructions.  \nAlthough recent efforts have been made to demonstrate solving the complex gain factor as a complex optimization for polarization and direction-dependent calibration (e.g. Smirnov & Tasse 2015; Grobler et al. 2018), we adopt the conventional approximation and treat the real and imaginary parts separately for our unpolarized case. Therefore, we can determine the blending-induced calibration error individually by solving the two overdetermined systems of equations and derive the per-frequency Δ 𝜂 𝑖 and Δ 𝜙 𝑖 for each antenna. A Python-based SVD solver is implemented to derive the solution following the exact derivation presented in Appendix B2. Because there is no contribution of baseline noise or possible biases from specific antenna-based calibration solvers, the solutions should be considered as a theoretical upper-limit estimation of the gain errors.  \nFinally, by supplying the per-frequency per-antenna relative residual gain errors back to equation (4) for each sky component (we note true visibility 𝑆 true 𝑖 𝑗 and ill-calibrated visibility ' 𝑆 ill 𝑖 𝑗 ), we can infer the per-frequency per-baseline propagation bias of the sky signal in the visibility space as:  \nΔ H 𝑖 𝑗 = GLYPH<2> exp GLYPH<0> Δ 𝜂 𝑖 + Δ 𝜂 𝑗 GLYPH<1> -1 GLYPH<3> GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> , (8a)  \nΔΦ 𝑖 𝑗 = Δ 𝜙 𝑗 -Δ 𝜙 𝑖 , (8b)  \nwhere Δ H 𝑖 𝑗 ≡ GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> and ΔΦ 𝑖 𝑗 ≡ arg GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -arg GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> denote the amplitude and phase bias, respectively. Based on equation (8), the derived amplitude and phase biases can be viewed as residuals after the removal of the corresponding foreground components without explicitly performing foreground subtraction. Because we use the true visibility of the foreground, the visibility bias is additive in nature, even after a perfect foreground removal. The detailed derivation is presented in Appendix B3. By utilizing the relative gain errors that are independently estimated for each frequency, the propagation bias can be inferred for all channels within each observation band. Thus, with the input visibility cubes [ 𝑉 ori , gal ( 𝑢, 𝑣, Δ 𝑓 ) , 𝑉 ori , ext ( 𝑢, 𝑣, Δ 𝑓 ) , and  \n𝑉 ori , eor ( 𝑢, 𝑣, Δ 𝑓 ) ], we can explicitly derive the visibility-bias cubes for each sky component, namely, 𝑉 res , gal ( 𝑢, 𝑣, Δ 𝑓 ) , 𝑉 res , ext ( 𝑢, 𝑣, Δ 𝑓 ) , and 𝑉 res , eor ( 𝑢, 𝑣, Δ 𝑓 ) , respectively.", '3.5 Deconvolution and imaging': "As the final fraction of the HEVAL pipeline, this section details mapmaking, which is a crucial final step to enter the measurement space under the 'reconstructed'-sky approach. To deal with the challenge of both wideband and wide-FOV imaging, the HEVAL pipeline uses the WSC/l.pc/e.pc/a.pc/n.pc 20 imager (Offringa et al. 2014) tailored for map-making in the low-frequency radio sky.  \nTo reconstruct 3D representations of the radio sky, image cubes along an observation frequency band are required. This is achieved by converting the 3D visibility cubes back to the image space. Both the 'original' visibility cubes (the lower part of Fig. 2) and the propagated 'residual' visibility-bias cubes (the upper part of Fig. 2) of each sky component are converted directly into 'original' [ 𝐼 ori ( 𝑙, 𝑚, Δ 𝑓 ) ] and 'residual' image cubes [ 𝐼 res ( 𝑙, 𝑚, Δ 𝑓 ) ], respectively. In particular, the direct imaging of foreground visibility-bias cubes, instead of ill-calibrated visibility cubes, allows an efficient map-making phase, given that imaging fidelity in the spatial domain is less a priority than spectral fidelity. Because we aim to evaluate the impact of frequency-dependent calibration errors, any spectral variations during the deconvolution phase need to be dealt with. The wideband deconvolution mode with multi-frequency weighting of the WSC/l.pc/e.pc/a.pc/n.pc (Offringa & Smirnov 2017) under a Briggs robustness weighting set to -0 . 5 (Briggs 1995) is used to account for spectral variations by griding the multi-frequency imaging weights together and cleaning all frequency channels jointly. In addition, a polynomial is fitted to the full bandwidth, reducing possible imaging noise and spectral artefacts that are typically caused by cleaning false peaks impersonated by side lobes or noise peaks. These measures allow HEVAL-produced foreground image cubes to mitigate possible spectral fluctuations, which might be coupled with the effect of frequency-dependent calibration errors in the measurement space. For the 21-cm signal, dirty images are used directly for analysis due to the insufficient deconvolution effect of the CLEAN algorithm to diffuse faint signals. All the final images are cropped to retain only the central 1 · × 1 · regions to mitigate imaging errors due to insufficient CLEANing of the marginal regions.", '3.6 The 2D power spectrum & EoR window': "As a redshifted line emission, observations of the EoR signal within a frequency span are three-dimensional, including two spatial dimensions across the sky and one vertical frequency dimension to the LoS. Statistical detection of the EoR signal under the 'reconstructed'-sky approach is expected to be achieved using the 3D nature of the EoR observation by converting reconstructed image cubes to 3D Fourier representations [ 𝑃 ( 𝑘 𝑥 , 𝑘 𝑦 , 𝑘 𝑧 ) ] in the measurement space via Fourier transform. The data flow from interferometric observations to the PS measurement space is presented in Fig. 8.  \nWithin the measurement space, PS analysis can be achieved by averaging 𝑃 ( 𝑘 𝑥 , 𝑘 𝑦 , 𝑘 𝑧 ) either in spherical shells of radii 𝑘 , due to its spherical symmetry, resulting in the 1D PS 𝑃 ( 𝑘 ) , or over the angular annuli of radii 𝑘 ⊥ ≡ √︃ 𝑘 2 𝑥 + 𝑘 2 𝑦 for each LoS plane 𝑘 ∥ ≡ 𝑘 𝑧 , owing to the independence between the spatial and frequency dimensions,  \nyielding the 2D PS 𝑃 ( 𝑘 ⊥ , 𝑘 ∥ ) . Theoretically, the cylindrical-averaged 𝑃 ( 𝑘 ⊥ , 𝑘 ∥ ) confines spectrally-smoothed foregrounds to be distributed within the low𝑘 ∥ region of the ( 𝑘 ⊥ , 𝑘 ∥ ) plane, leaving the rest of the space relatively free of foreground contamination. Thus, the 2D PS analysis is widely adopted as the most crucial figure of merit for experiments aimed to detect the EoR signal.  \nHowever, the reality of complicated instrumental and observational effects, such as chromatic primary beams, redistributes the foreground power to higher 𝑘 ∥ dimensions via mode mixing, resulting in an expanded wedge-shaped region known as the foreground wedge (Datta et al. 2010; Morales et al. 2012; Pober et al. 2014). Beyond the foreground wedge, a region known as the EoR window (for more details on the topic, please refer to Liu et al. 2014; Dillon et al. 2014), is relatively free of foreground power contamination and can be located at  \n𝑘 ∥ ≥ 𝐻 ( 𝑧 ) 𝐷 𝑀 ( 𝑧 ) ( 1 + 𝑧 ) 𝑐 GLYPH<20> 𝑘 ⊥ sin Θ + 2 𝜋𝑤 𝑓 21 -rf ( 1 + 𝑧 ) 𝐷 𝑀 ( 𝑧 ) 𝐵 GLYPH<21> (9)  \n(Thyagarajan et al. 2013; Li et al. 2019), where 𝐻 ( 𝑧 ) denotes the Hubble parameter at redshift 𝑧 , 𝐷 𝑀 ( 𝑧 ) measures the transverse comoving distance, 𝑐 is the speed of light, 𝐵 is the bandwidth of the frequency window, 𝑓 21 -rf denotes the frequency of the 21-cm line emission in the rest-frame, 𝑧 = 𝑓 21 -rf / 𝑓 𝑐 -1 is the signal redshift of the central frequency 𝑓 𝑐 within the measured frequency band, 𝑤 is the number of characteristic convolution widths of the spillover region due to frequency response variations, and Θ is the foreground source angular separation from the primary beam centre. Based on the angular separation Θ , commonly used EoR window boundaries include 'FoV'- and 'horizon'-bound windows, which are defined by supplying the FoV and horizon of the instrument.  \nFrequency-dependent calibration errors due to source blending are expected to aggravate the effect of mode mixing for 'reconstructed'sky PS, which will cause further foreground contamination by leaving additive 'residual' powers within EoR windows. Thus, to fully evaluate the impact of source blending on the EoR experiments, we discuss the 3D imprint (both the spatial and spectral deviations) of the blending-induced visibility bias with the 'original' [ 𝑃 ori ( 𝑘 ⊥ , 𝑘 ∥ ) ] and 'residual' 2D PS [ 𝑃 res ( 𝑘 ⊥ , 𝑘 ∥ ) ] for each sky component, which is calculated using the 'original' [ 𝐼 ori ( 𝑙, 𝑚, Δ 𝑓 ) ] and 'residual' image cubes [ 𝐼 res ( 𝑙, 𝑚, Δ 𝑓 ) ], respectively.", "4.1 'Residual' powers of EoR foregrounds": "We calculate the cylindrical-averaged 2D PS using image cubes produced by the HEVAL pipeline. In Fig. 9, we present the 'residual' 2D PS [ 𝑃 res , ext ( 𝑘 ⊥ , 𝑘 ∥ ) and 𝑃 res , gal ( 𝑘 ⊥ , 𝑘 ∥ ) ] transformed from foreground 'residual' image cubes [ 𝐼 res , ext ( 𝑙, 𝑚, Δ 𝑓 ) and 𝐼 res , gal ( 𝑙, 𝑚, Δ 𝑓 ) ] for each component owing to the propagation effects of blending-induced calibration errors. Given that these powers will still be present in the PS space even after the perfect removal of spectrally-smoothed foregrounds, we refer to them as 'residual' powers. 'Residual' powers originating from the 'highly'-corrupted sky model (5 per cent blending ratio) show a substantial occupation of the low 𝑘 ∥ modes for both the Galactic and extragalactic foregrounds. The strongest Galactic 'residual' powers reside at larger spatial scales ( 𝑘 ⊥ ≲ 0 . 2), while the strong extragalactic 'residual' powers spread across the large and intermediate spatial scales. For the high 𝑘 ∥ modes, we can see a clear excess of powers within the EoR window above both the dashed black 'FoV' line and solid black 'horizon' line. These 'residual' powers contaminate the expected foreground-free  \nregions. Hence, they pose a threat to the parameter inference of the EoR signal. For the 'moderate' case (0.5 per cent blending ratio), the 'residual' powers, albeit weaker, distribute similarly to that of the 'high' blending ratio case in the 2D PS space. Unlike the prior two cases, the 'mildly'-corrupted sky model (0.05 per cent blending ratio) introduces much weaker 'residual' powers with a slightly different vertical distribution, as there is no characteristic distribution of 'residual' powers occupying low 𝑘 ∥ modes for both the Galactic and extragalactic foregrounds. However, the contamination of the 'residual' power within the EoR window remains.  \nBytaking the 2D PS ratio between the ill-calibrated 'residual' power [ 𝑃 res , ext ( 𝑘 ⊥ , 𝑘 ∥ ) and 𝑃 res , gal ( 𝑘 ⊥ , 𝑘 ∥ ) ] and 'original' sky emission power [ 𝑃 ori , ext ( 𝑘 ⊥ , 𝑘 ∥ ) and 𝑃 ori , gal ( 𝑘 ⊥ , 𝑘 ∥ ) ] for each foreground component, we can identify the region with excess power caused by blending-induced calibration errors. Given the nature of frequencydependent calibration errors that cause relatively small fluctuations across the frequency channels, the region with excess power should be located at higher 𝑘 ∥ modes. In Fig. 10, we can see excess powers locate at the region covering important high 𝑘 ∥ modes ( 𝑘 ∥ ≳ 0 . 2), as expected. Among the three cases, the 'mild' blending case has the lowest excess power. All the 'residual' powers introduced by the 'mildly'-corrupted sky model are below 10 per cent of the 'original' sky emission powers. In contrast, the other two cases pose significant contamination by introducing excess foreground 'residual' powers, which are 10 to over 1000 times the 'original' sky emission powers, covering most of the EoR window. Outside the 'FoV' window, foreground 'original' powers still dominate the wedge region for all three cases. Given that we treat the two foreground components separately, we can also calculate the 2D PS ratio between the extragalactic and Galactic 'residual' powers to infer the dominant modes of each foreground component. Fig. 11 illustrates the spatial scales that divide the dominance. The Galactic 'residual' errors dominate the spatial scales larger than 𝑘 ⊥ ∼ 0 . 2, while the extragalactic 'residual' errors dominate the remaining 𝑘 ⊥ modes.", '4.2 Blending impact on the EoR PS and blending ratio tolerance': "The direct quantification of the blending-induced 'residual' power impact on EoR detection is achieved by calculating the 2D PS ratio between EoR 'residual' powers and 'original' powers. We can infer the specific impact by analysing the contamination within the estimated foreground-free EoR window. We plot the 𝑅 res , fg / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) ratio for the two EoR foregrounds for the three blending cases. From the plotted 𝑅 res , gal / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) in Fig. 12, the Galactic 'residual' powers dominate the large spatial scales, polluting the most important EoR 𝑘 ⊥ modes significantly. For both the 'high' (5 per cent) and 'moderate' (0.5 per cent) blending ratio cases, strong Galactic 'residual' powers, which are at least 1000 times the corresponding EoR powers, cover almost all the modes within the 'horizon'-bound region. Unanticipatedly, Galactic 'residual' powers also contaminate the EoR powers at the large 𝑘 ⊥ end with a wedge-shaped region at 𝑘 ∥ ≳ 7 and 𝑘 ∥ ≳ 9 for the 'high' and 'moderate' cases, respectively. Recall that the simulated Galactic diffuse foregrounds include small-scale fluctuations down to the SKA-Low's angular resolution, which enables PS analysis to infer the Galactic impact below the arcmin scale for the first time. Although the extragalactic 'residual' powers dominate at intermediate and small spatial scales, their impact within the EoR window remains significant viewing from the plotted 𝑅 res , ext / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) in Fig. 12. For the 'high' blending ratio case, all of the modes within the 'FoV'-bound region are not recoverable under the impact of either of the foreground components. Interestingly, there is a small window of opportunity for the 'moderate' case since there are modes,  \nFigure 8. The data flow of our 'reconstructed'-sky PS analysis pipeline. The multi-frequency visibilities of a given band are transformed into image cubes reconstructing the sky emission through the Fourier transform in the sky coordinates. Subsequently, the reconstructed image cubes are transformed into the measurement space through the Fourier transform in both the sky coordinates and the frequency dimension to analyse the spatial features of each sky component. Both 'original' visibility and 'residual' visibility-bias cubes for each sky component are transformed following the data flow. The 𝑉 ori 'original' visibility cubes are converted to 𝐼 ori and eventually transform into 𝑃 ori . The 𝑉 res 'residual' visibility-bias cubes follow the same trajectory to the image and PS space as 𝐼 res and 𝑃 res , respectively. This plot is inspired by fig. 2 of Morales & Hewitt (2004).  \n<!-- image -->  \nResidual Power of Extraglactic and Galactic ForegroundsFigure 9. (Left, a-i) Extragalactic 'residual' powers 𝑃 res , ext ( 𝑘 ⊥ , 𝑘 ∥ ) and (right, j-i) Galactic 'residual' powers 𝑃 res , gal ( 𝑘 ⊥ , 𝑘 ∥ ) originating from blending-induced calibration errors. Within each 3 × 3 grid, the three rows mark the 'high' (5 per cent), 'moderate' (0.5 per cent), and 'mild' (0.05 per cent) blending scenarios, respectively; the three columns mark the 188 - 196 MHz, 154 - 162 MHz, and 120 - 128 MHz frequency band, respectively. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. All the plots share the same logarithmic colour bar using unit [mK 2 Mpc 3 ].  \n<!-- image -->  \nwhich is located at 0 . 6 ≲ 𝑘 ⊥ ≲ 5 . 0 and 0 . 2 ≲ 𝑘 ∥ ≲ 1 . 0, relatively free of Galactic 'residual' power. However, this small area vanishes after combining the Galactic and extragalactic foregrounds. For the 'mild' (0.05 per cent) blending ratio case, significant less 'residual' powers are present in the measurement space. For one thing, only  \nthe largest spatial scales, which are located at 0 . 04 ≲ 𝑘 ⊥ ≲ 0 . 08 and 0 . 3 ≲ 𝑘 ∥ ≲ 1 . 0, are affected by the Galactic 'residual' power, leaving most of the 'FoV'-bound region relatively free of impact. For another, the extragalactic 'residual' power leaves all 'FoV'-bound regions free of contamination by polluting only the smaller spatial scales outside", "Extragalactic 'Residual' Power / Extragalactic 'Original' Power": "Figure 10. The 154 - 162 MHz 2D PS ratio between foreground 'residual' and 'original' powers 𝑅 res , fg / ori , fg ( 𝑘 ⊥ , 𝑘 ∥ ) . (Top) 𝑅 res , ext / ori , ext ( 𝑘 ⊥ , 𝑘 ∥ ) ratios under (a) 'high' (5 per cent), (b) 'moderate' (0.5 per cent), and (c) 'mild' (0.05 per cent) blending scenarios. (Bottom) Corresponding 𝑅 res , gal / ori , gal ( 𝑘 ⊥ , 𝑘 ∥ ) ratios under (d) 'high' (5 per cent), (e) 'moderate' (0.5 per cent), and (f) 'mild' (0.05 per cent) blending scenarios. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->", "Extragalactic 'Residual' Power / Galactic 'Residual' Power": "Figure 11. The 154 - 162 MHz 2D PS ratios between extragalactic and Galactic 'residual' powers 𝑅 res , ext / res , gal ( 𝑘 ⊥ , 𝑘 ∥ ) under (a) 'high' (5 per cent), (b) 'moderate' (0.5 per cent), and (c) 'mild' (0.05 per cent) blending scenarios. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->  \nthe EoR window. Although caution should be exercised when dealing with Galactic 'residual' pollution, the impact of 'residual' powers on EoR detection is insignificant.  \nThe HEVAL pipeline is designed to isolate source blending as the sole origin of errors to evaluate its impact on SKA EoR experiments under a sky-based calibration scheme. In the image space, the implementation of B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc simulates defects that can be attributed solely to source blending. In the visibility space, our analytical modules employ a relative calibration scheme to further segregate the calibration errors of blending-induced origin. During the map-making phase, wideband, multi-frequency, and spectral fitting procedures mitigate possible spectral artefacts that may cloud the impact evaluation in the measurement space. Given these measures of impact isolation, we can confidently conclude that the inference in the 2D PS space presented in this section originates from blending defects and propagates along the entire EoR analysis pipeline. Consider that the 'residual' powers are propagated from 'residual' visibility-biases that are present even after a perfect foreground subtraction, we can draw a blending ratio tolerance to guide the sky-model construction for the calibration of SKA EoR experiments. While caution should be taken when addressing Galactic 'residual' pollution, the blending ratio tolerance for the SKA EoR experiment can be drawn at or slightly below the 0.05 per cent level of the 'mildly'-blended case, which translates to containing only 5 pairs of blended sources per 10000 sources for sky models aimed at calibration.", '5.1 Limitations of the study': "Although this paper offers the first systematic insight into source blending in the upcoming SKA era, there are limitations to our study.  \nNoise consideration The simulation of the source-blending effect should be considered as the best-case scenario, as the noise impact of the sky-model construction is not considered. Given that noise plays a vital role during the detection phase of sky-model construction (i.e. the noise level determines the detection criteria), a mixture of components often results in a measurement error in the total flux density. Considering the accuracy of flux density measurements as its own source of sky-model defects, our implantation of source-blending effects only considers the flux density reallocation without affecting the total flux density of the blended source. Hence, the total flux density of the ' ideal ' and ' blended ' from the same pair remains the same.  \nEffects Decoupling One of the key considerations of this study is to infer the calibration errors solely from source-blending defects, given that sky-model defects often couple with each other, instrumental effects, or a combination of both in practical calibration pipelines. Hence, the realization of the HEVAL pipeline in this study decouples source-blending defects from other factors, such as polarimetric and side-lobe impact. In particular, source-blending defects from sources in the side lobes may introduce additional calibration-induced errors in the 2D PS space (see Li et al. 2019, for a discussion on side-lobe effects in the 2D PS space), which may contaminate different parts of the EoR window. We will consider these topics in future studies because of the high computational expenses required to evaluate these coupled effects.  \nTemporaleffects Themainaimofthisstudyistoinferthe frequencydependent gain error. The temporal effects of the complex gain are its  \nown topic (see Kumar et al. 2020, 2022, for detailed discussions on time-correlated gain effects). We omit consideration of the temporal effect of the complex gain while considering the calibration process. Therefore, our results should be viewed as a time-averaged scenario.  \nSmall-scale simulations The low-frequency radio data used in this study is simulated to include physical scales down to the estimated SKA1-Low spatial resolution. Although a commonly used practice, which effectively redistributes powers to the simulated small scales, is applied to add small-scale fluctuations to our Galactic components, Galactic observations to constrain the simulation at such scales are still lacking. For extragalactic foregrounds, our implemented Gaussian models of the EDRS populations indicate that fine structures below 6 arcsec are excluded from this study. We consider the fine structures of the extended sources as their own topic (i.e. the fidelity of the extended source model) and will discuss the impact in an upcoming project. Interested readers can refer to Procopio et al. (2017) for a discussion of foreground spatial fidelity impact on EoR foreground subtraction. As for the small-scale fluctuations of the EoR simulation, the accuracy of our data is subject to the limitations of the semi-analytic models used in 21/c.pc/m.pcFAST.", "5.2 Impact of EoR 'residual' powers": "Similar to the calculation of the 2D PS ratio between the foreground 'residual' and the EoR signal, we can plot the EoR 'residual' and 'original' power ratio 𝑅 res , eor / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) and infer the impact of EoR 'residual' powers on the EoR signal itself. As demonstrated in Fig. 13, only the 'high' (5 per cent) blending ratio case has introduced small excess 'residual' power within the EoR window. However, these powers are expected to have little impact on the EoR detection because they are smaller than one-tenth of the EoR signal and can be considered to be below the noise level. The only strong impact is outside the EoR window down to smaller scales, translating to little effect under the foreground avoidance strategy. Therefore, unlike the foreground 'residual' powers, the propagated EoR 'residual' powers originating from the 'highly'-corrupted sky model have little to no direct impact on the EoR detection. However, there might be impact on the accuracy of the EoR parameter inference (such as parameters related to the non-Gaussianity), depending on how precise the model can be constrained under the influence of the EoR 'residual' powers.", '5.3 Impact of imaging weighting on calibration errors': "Map-making is a critical step in the 'reconstructed'-sky PS analysis pipeline. Even with the advanced baseline coverage of the upcoming SKA1-Low, the interferometer can only discretely sample the sky without unlimited baselines. Thus, image weights are required to determine the filling of the telescope's sampling gap during the mapmaking phase. We discuss the impact of imaging-weight choice on the propagation effects of blending-induced calibration errors. We consider the 'moderate' (0.5 per cent) blending case, which is above the blending ratio tolerance, and add Briggs-weighted clean maps with different robustness, including -1, 0, 0.5, and 1, in addition to the -0.5 robustness used in Section 3.5. We retain the remaining imaging settings and change only the weight. Thus, we have 5 differently weighted maps for the 'residual' and 'original' image cubes of each foreground component. Fig. 14 plots the 2D PS ratio between foreground 'residual' powers and EoR 'original' powers across the 5 different weighting options. There are clear trends across the different imaging options: (i) there is a reduction in 'residual' power at the", "Extragalactic 'Residual' Power / EoR 'Original' Power": "Figure 12. The 154 - 162 MHz 2D PS ratios between foreground 'residual' powers and EoR 'original' powers 𝑅 res , fg / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) . (Top) 𝑅 res , ext / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) ratios under (a) 'high' (5 per cent), (b) 'moderate' (0.5 per cent), and (c) 'mild' (0.05 per cent) blending scenarios. (Bottom) Corresponding 𝑅 res , gal / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) ratios under (d) 'high' (5 per cent), (e) 'moderate' (0.5 per cent), and (f) 'mild' (0.05 per cent) blending scenarios. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->  \nEoR 'Residual' Power / EoR 'Original' PowerFigure 13. The 154 - 162 MHz 2D PS ratios 𝑅 res\\_eor / ori\\_eor ( 𝑘 ⊥ , 𝑘 ∥ ) (between EoR 'residual' and 'original' powers) under (a) 'high' (5 per cent), (b) 'moderate' (0.5 per cent), and (c) 'mild' (0.05 per cent) blending scenarios. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->  \nsmallest spatial scales with robustness varying from -1 to 1; (ii) the Galactic 'residual' power at the intermediate spatial scales increases with robustness varying from -1 to 1; and (iii) the 'residual' power at the largest spatial scales increases as the robustness varies from -1 to 1. Given these trends, we can identify a potential mitigation strategy for blending-induced 'residual' powers by choosing the suitable imaging weight for the particular observing sky during the map-making phase. As shown in Fig. 15, a viable approach is using the robustness -1 option for this particular simulated sky patch considering the reduced impact in the measurement space. There is an evident boost in the detection possibility at the largest spatial scales, which can be inferred from the weighting ratio figures in the second row of Fig. 15. From a general EoR detection point-of-view, a more holistic approach should be considered by balancing between different sources of errors across the analysis pipeline with each individual bias factor identified, evaluated, and mitigated.", '5.4 Mitigation of blending-induced effects': "To reduce the propagated 'residual' power discussed in this paper, upcoming SKA EoR experiments require dedicated de-blending efforts built directly within the sky-model construction pipeline. Current de-blending efforts based on LOFAR observations heavily involve human input (see Williams et al. 2019; Kondapally et al. 2021, for the LOFAR de-blending workflow). With the upcoming SKA1Low's improved sensitivity compared to the current observations, source blending could only be exacerbated. Therefore, automatic de-blending approaches, such as those applied in the optical domain (see Melchior et al. 2021, for a review), are required to address the blending issue both in decomposition and association, but the implementation is still an open question. To aid the development of de-blending pipelines and strategies, the developer community needs proper datasets through observations and simulations. From an observational perspective, obtaining higher-resolution radio observations of the SKA EoR fields through the international LOFAR Telescope (Vermeulen & van Haarlem 2011; Jackson et al. 2016) would both offer a direct blending effect deduction and help with the development of de-blending methods during decomposition. Multi-frequency synergies, such as joint observations with JWST 21 (McElwain et al. 2023), the Large Synoptic Survey Telescope (LSST 22 ; Ivezić et al. 2019), or the Chinese Survey Space Telescope (CSST 23 ; Zhan 2011), of the SKA observing sky, would also be crucial for the SKA to battle source blending during the source-association phase. From a simulation perspective, a high-fidelity multi-frequency (across the radio to infrared and optical bands) EDRS simulation is required to develop de-blending methods. An alternative approach to direct de-blending efforts is to treat blending-induced 'residual' errors in a similar manner to other frequency-dependent errors. General frequency-dependent error-mitigation strategies include the usage of spectrally-smoothed antennas (Barry et al. 2016), baseline weighting (Ewall-Wice et al. 2017), and simultaneous multi-channel calibration (Yatawatta 2015). In practice, a combination of both a direct de-blending pipeline and indirect error-mitigation strategies might be used to achieve an optimal solution that effectively attenuates the blending-induced impact.", '6 CONCLUSION': "Wemake the first attempt to systematically assess the impact of source blending in the low-frequency radio band for the upcoming SKA1Low. This study introduces a clear definition of source blending in the radio window, identifies the two-fold blending defect in both the spatial and frequency domains, and investigates the underlying impact of source blending on interferometric calibrations for SKA EoR experiments. We summarize our work as follows:  \n- (i) Using the quick-and-dirty method presented in Section 2, which utilizes the W08 source model and the latest SKA1-Low array configuration, we estimate an extended level of blending ( ∼ 5 -28 per cent sources are blended) for the upcoming SKA1-Low observing sky and identified source blending as one of the key imperfections impacting the sky-model construction for SKA EoR experiments.\n- (ii) ByusingourHEVALpipeline(whichsimulatesthelow-frequency radio sky containing physical scales spatially resolvable by the SKA1Low) to quantify the blending-induced calibration impact in a skybased scheme for EoR detections, we find that frequency-dependent calibration errors from poor calibration against blending-corrupted sky models coupled with strong foregrounds leave additive 'residual' powers within the EoR window in the 2D PS space and may significantly impede EoR detections. Our findings corroborate with both the B16 simulation approach and the E17 analytic analysis results, albeit focusing on three different sky-model defects. Furthermore, by considering the extragalactic and Galactic foregrounds separately, we conclude that the extragalactic and Galactic 'residual' powers contaminate EoR scales at larger 𝑘 ⊥ and smaller 𝑘 ⊥ , respectively, with a clear division around 𝑘 ⊥ ∼ 0 . 2.  \n(iii) To determine the blending tolerance for SKA EoR experiments, we perform three tests with relatively low levels of blending ratios (5, 0.5, and 0.05 per cent). The tests show that the blending defects at the former two levels introduce strong 'residual' powers in the measurement space, seriously contaminating the key EoR scales (0 . 1 ≲ 𝑘 ≲ 2 Mpc -1 ) within the EoR window, whereas sky models with a blending ratio at 0.05 per cent leave little impact in the measurement space. Hence, a blending ratio of approximately 0.05 per cent is identified as the blending tolerance for the upcoming SKA1-Low sky-model construction.  \n- (iv) Given that our estimated 5 -28 per cent blending ratio for the SKA1-Low has quite a large margin with a tolerance of 0.05 per cent, we discuss the possibility of mitigating and suppressing the blendinginduced impact through different imaging weights. For the designed baseline of the SKA1-Low, Briggs weighting with robustness -1 would be a better choice to mitigate the blending-induced impact at key EoR scales.  \n(v) To directly mitigate the impact of source-blending defects, deblending techniques must be involved in both the SKA source detection and sky-model construction to pass the required blending tolerance. From a software perspective, the source detection community should prioritize including de-blending or blending flagging during both the component-decomposition and -association phases. Additionally, priority should be set to port and adapt existing automatic de-blending methods, such as those used in the optical bands, to the radio bands. From an observation perspective, adding higher-resolution or lessblended observations to the SKA EoR fields via synergies with the international LOFAR Telescope and instruments from other wavelength bands will be crucial for the SKA to reduce the blending level through cross-match assessments directly. In practice, joint efforts from the software and observation perspectives are most likely required to contain blending-induced errors within the necessary overall EoR error budget.  \nFigure 14. The 154 - 162 MHz 2D PS ratios between foreground 'residual' powers and EoR 'original' powers 𝑅 res , fg / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) under the 'moderate' (0.5 per cent) blending scenario using different Briggs weighting options. (Top) 𝑅 res , ext / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) ratios using (a) -1, (b) -0.5, (c) 0, (d) 0.5, and (e) 1 robustness options. (Bottom) Corresponding 𝑅 res , gal / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) ratios using (f) -1, (g) -0.5, (h) 0, (i) 0.5, and (j) 1 robustness options. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->", 'ACKNOWLEDGEMENTS': "We thank Weitian Li for valuable discussions regarding /f.pc/g.pc21/s.pc/i.pc/m.pc, Andrei Mesinger, Yuxiang Qin, and Steven Murray for their help with 21/c.pc/m.pcFAST, Fred Dulwich for technical support in resolving issues with OSKAR, and André Offringa for guidance on map-making with WSC/l.pc/e.pc/a.pc/n.pc. We thank Hekun Lee, Benjamin McKinley, and Junhua Gu for scientific discussions, Sanjay Bhatnagar for reading the manuscript and providing insightful comments, and Bingxue Yang's help with the preparation of the final draft of this manuscript. This work was performed on the Gravity Cluster of Department of Astronomy (DOA) at Shanghai Jiao Tong University (SJTU). Some of the results in this paper have been derived using the healpy and HEALPix packages. This work is supported by the Ministry of Science and Technology of China (Grant No. 2020SKA0110201), and the National Natural Science Foundation of China (Grant Nos. 11973033 and 11835009). ZZH acknowledges the support from the National Science Foundation of China (Grant No. 12203085). SVW acknowledges the financial assistance of the South African Radio Astronomy Observatory (SARAO; https://www.sarao.ac.za). Parts of this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.", 'DATA AVAILABILITY': 'The end-to-end simulate suite F/g.pc21S/i.pc/m.pc+ is an open-source project developed at https://github.com/Fg21Sim/ . Owing to the large  \ndata volume, the data generated by this work will be shared under a reasonable request at https://github.com/Fg21Sim/Data . Datasets with the same specifications as in this study can be simulated using F/g.pc21S/i.pc/m.pc+ with the parameters provided by the authors.', 'AUTHOR CONTRIBUTION STATEMENT': "Here is the list of contributions:  \n- · C.Shan: project proposal, coding (blending ratio estimation and the majority of HEVAL pipeline), methodology, testing, validation, result analysis, and writing (original draft and modification).\n- · H.Xu: scientific feedback, results discussion, writing (feedback and modification), and funding support.\n- · Y.Zhu: coding (contribution to the ES/i.pc/m.pc and B/l.pc/e.pc/n.pc/d.pcS/i.pc/m.pc modules), validation, results analysis (feedback and discussion), and writing (feedback and modification).\n- · Y.Zhao: scientific feedback, results analysis (feedback and discussion), and writing (feedback and modification).\n- · S.White: scientific feedback, figure modification, and writing (feedback and modification).\n- · J.Line: scientific feedback, results analysis (feedback and discussion) and writing (feedback).\n- · D.Zheng: methodology (feedback and discussion) and coding (contribution to the GS/i.pc/m.pc).\n- · Z.Zhu: scientific feedback and writing (feedback).\n- · D.Hu: scientific feedback and writing (feedback).\n- · Z.Zhang: scientific feedback.\n- · X.Wu: scientific and funding support.  \nFigure 15. (Top) The 154 - 162 MHz 2D PS ratios between combined foreground 'residual' powers and EoR 'original' powers 𝑅 res , ext + gal / ori , eor ( 𝑘 ⊥ , 𝑘 ∥ ) under the 'moderate' (0.5 per cent) blending scenario using (a) -1, (b) -0.5, (c) 0, (d) 0.5, and (e) 1 Briggs robustness options, respectively. (Bottom) The 154 - 162 MHz 2D PS ratios between combined foreground 'residual' powers using various Briggs robustness options and the Briggs robustness 0 option 𝑅 res , Briggs ( i )/ res , Briggs ( 0 ) ( 𝑘 ⊥ , 𝑘 ∥ ) under the 'moderate' (0.5 per cent) blending scenario. The various robustness options include (f) -1, (g) -0.5, (h) 0, (i) 0.5, and (j) 1. Setting 𝑤 as 2, the black lines define the EoR windows with the dashed and solid lines marking the 'FoV' and the 'horizon' boundary, respectively. Figs. 10 to 15 share the same logarithmic colour bar, showing only six orders of magnitude for easy viewing.  \n<!-- image -->", 'REFERENCES': "Arcelin B., Doux C., Aubourg E., Roucelle C., LSST Dark Energy Science Collaboration 2021, MNRAS, 500, 531  \nAsad K. M. B., Koopmans L. V. E., Jelić V., de Bruyn A. G., Pandey V. N., Gehlot B. K., 2018, MNRAS, 476, 3051  \nBarry N., Hazelton B., Sullivan I., Morales M. F., Pober J. C., 2016, MNRAS, 461, 3135  \nBhatnagar S., Rau U., Golap K., 2013, ApJ, 770, 91  \nBlandford R., Meier D., Readhead A., 2019, ARA&A, 57, 467  \nBowman J. D., et al., 2013, Publ. Astron. Soc. Australia, 30, e031  \nBriggs D. S., 1995, PhD thesis, New Mexico Institute of Mining and Technology  \nBrunton S. L., Kutz J. N., 2019, Singular Value Decomposition (SVD). Cambridge University Press, p. 3-46  \nByrne R., et al., 2019, ApJ, 875, 70  \nByrne R., Morales M. F., Hazelton B. J., Wilensky M., 2021, MNRAS, 503, 2457  \nCarroll P. A., et al., 2016, MNRAS, 461, 4151  \nChapman E., Zaroubi S., Abdalla F. B., Dulwich F., Jelić V., Mort B., 2016, MNRAS, 458, 2928  \nDatta A., Bhatnagar S., Carilli C. L., 2009, ApJ, 703, 1851  \nDatta A., Bowman J. D., Carilli C. L., 2010, ApJ, 724, 526  \nDeBoer D. R., et al., 2017, PASP, 129, 045001  \nDelabrouille J., et al., 2013, A&A, 553, A96  \nDewdney P., Braun R., 2016, Technical report, SKA1-Low Configuration Coordinates - Complete Set. SKA Organisation, Design Report SKATEL-SKO-0000422  \nDewdney P., et al., 2013, Technical report, SKA1 System Baseline Design, Revision 1. SKA Organisation, Design Report SKA-TEL-SKO-DD-001 Dewdney P., et al., 2016, Technical report, SKA1 System Baseline Design V2.  \nSKA Organisation, Design Report SKA-TEL-SKO-0000002  \nDewdney P., et al., 2019, Technical report, SKA1: Design Baseline Description.  \nSKA Organisation, Design Report SKA-TEL-SKO-0001075  \nDickinson C., Davies R. D., Davis R. J., 2003, MNRAS, 341, 369  \nDillon J. S., et al., 2014, Phys. Rev. D, 89, 023002  \nDillon J. S., et al., 2018, MNRAS, 477, 5670  \nEwall-Wice A., Dillon J. S., Liu A., Hewitt J., 2017, MNRAS, 470, 1849  \nFanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P  \nFinkbeiner D. P., 2003, ApJS, 146, 407  \nForeman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306  \nFortet R., Mourier E., 1953, Annales scientifiques de l'École Normale Supérieure, 3e série, 70, 267  \nFurlanetto S. R., 2016, in Mesinger A., ed., Understanding the Epoch of Cosmic Reionization: Challenges and Progress. Springer International Publishing, Cham, pp 247-280  \nFurlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181  \nGehlot B. K., et al., 2018, MNRAS, 478, 1484  \nGehlot B. K., et al., 2021, MNRAS, 506, 4578  \nGiardino G., Banday A. J., Górski K. M., Bennett K., Jonas J. L., Tauber J., 2002, A&A, 387, 82  \nGórski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759  \nGrobler T. L., Nunhokee C. D., Smirnov O. M., van Zyl A. J., de Bruyn A. G., 2014, MNRAS, 439, 4030  \nGrobler T. L., Bernardi G., Kenyon J. S., Parsons A. R., Smirnov O. M., 2018, MNRAS, 476, 2410  \nHales C. A., 2017, AJ, 154, 54  \nHancock P. J., Trott C. M., Hurley-Walker N., 2018, Publ. Astron. Soc. Australia, 35, e011", 'APPENDIX A: SMALL-SCALE FLUCTUATION OF GALACTIC FOREGROUNDS': 'In this section, we introduce our method of adding small-scale fluctuations to the Galactic emission in detail. First, GS/i.pc/m.pc extrapolates the angular PS of the basis templates to the smaller scales and generates GRFs, 𝐺 grf , based on the extrapolated angular PS. We adopt the angular power spectrum model for the synchrotron and free-free emission as  \nC sync ℓ = ℓ 𝛾 h 1 -exp GLYPH<16> -ℓ 2 𝜎 2 temp GLYPH<17> i (A1)  \n(Remazeilles et al. 2015) and  \nC free ℓ = ℓ 𝛾 h exp GLYPH<16> -ℓ 2 𝜎 2 sim GLYPH<17> -exp GLYPH<16> -ℓ 2 𝜎 2 temp GLYPH<17> i (A2)  \n(Delabrouille et al. 2013), respectively, where 𝜎 temp and 𝜎 sim are the full width at half-maximum (FWHM) beam of the basis template and simulated sky maps, respectively, and the 𝛾 index is the power law index of the basis template. By adopting the two models, the GRF maps ( 𝐺 grf ) of the two Galactic diffuse components can be simulated using GS/i.pc/m.pc by converting C ℓ to maps with angular scales from degree level down to the required ℓ . To compensate for the high demand of both memory and CPU for generating high spatial resolution GRFs, we employ the GS/i.pc/m.pc module with a custom /c.pc/l.pc2/a.pc/l.pc/m.pc tool modified from the /h.pc/e.pc/a.pc/l.pc/p.pc/y.pc 24 package (Zonca et al. 2019) and /d.pc/u.pc/c.pc/c.pc0 25 with multi-threaded support. The former calculates the spherical harmonic coefficients ( 𝑎 𝑙𝑚 ) from the angular PS, and the latter offers map realization from the coefficients. Then, the small-scale fluctuation map ( 𝐺 ss ) is generated with a whitened 𝐺 grf , which has zero mean and unit variance, using  \n𝐺 ss = 𝛼𝐺 grf 𝐺 𝛽 temp , (A3)  \nwhere 𝛼 and 𝛽 , which dictate the mean and variance of the small-scale maps, respectively, are custom parameters to be determined. Finally, the simulated Galactic emission with the addition of small-scale fluctuations ( 𝐺 sim ) can be inferred with  \n𝐺 sim = 𝐺 temp + 𝐺 ss ,  \n= 𝐺 temp + 𝛼𝐺 grf 𝐺 temp .  \n(A4a) 𝛽 (A4b)  \nFor each simulation run, GS/i.pc/m.pc offers on-the-fly fitting of the three key parameters, 𝛼 , 𝛽 , and 𝛾 , to prescribe the generation of sky maps. The fitting process is based on the /e.pc/m.pc/c.pc/e.pc/e.pc 26 (Foreman-Mackey et al. 2013) Python library of the Markov Chain Monte Carlo (MCMC) method. To fit the power law index 𝛾 , GS/i.pc/m.pc employs the angular PS analysis of the Galactic all-sky basis template on a Hierarchical Equal Area and isoLatitude Pixelization (HEALPix 27 , Górski et al. 2005) sphere using the /a.pc/n.pc/a.pc/f.pc/a.pc/s.pc/t.pc 28 module. Subsequently, the 𝛾 index is inferred by fitting a power law to the large angular scale ( ℓ ≲ 500) of the computed C ℓ . As for the fitting of 𝛼 and 𝛽 , we use a set of statistical constraints, which require the preservation of the underlying statistical properties of the basis map with the addition of small-scale  \nfeatures (equation A4b). To achieve these constraints, GS/i.pc/m.pc uses a model that combines (i) an equivalence test, which is implemented through a two one-sided t-tests (TOST, Schuirmann 1987) procedure, (ii) a maximum mean discrepancy (MMD, Fortet & Mourier 1953) test, and (iii) difference minimization of the mean, variance, skewness, and kurtosis between the basis template ( 𝐺 temp ) and the feature-added map ( 𝐺 temp + 𝛼𝐺 grf 𝐺 𝛽 ).  \ntemp  \nFor the purpose of this work, the 𝛾 index of synchrotron and free-free component is fitted as -2.220 and -2.426, respectively, using C ℓ with ℓ ranging from 30 to 90. However, due to the high spatial resolution requirement of the small-scale features (a HEALPix map with Nside = 32,768 meets the SKA spatial resolving power), we fit 𝛼 and 𝛽 using partial sky coverage instead of the full sky map to avoid the extremely high CPU and memory costs. By fitting to the aimed sky coverage (R.A., Dec. = 0 · , -27 · ), we find the best-fitting result for the synchrotron ( 𝛼 = 0 . 0342, 𝛽 = 0 . 227) and free-free emission ( 𝛼 = 0 . 00785, 𝛽 = 0 . 526).', 'APPENDIX B: DERIVING THE GAIN ERROR SOLUTION & PER-BASELINE BIAS': 'This section details our formalism of the analytic analysis to evaluate the blending-induced calibration impact. Recall that we have established the sky-based per-frequency per-antenna calibration scheme in Section 3.4.1 and presented key equations in Section 3.4.2. Here, we start by deriving the linear systems of equations for the per-antenna gain amplitude and phase error in Appendix B1 using the PWL logarithmic implementation 29 . Then, Appendix B2 presents the SVDbased method used to infer the least-squares solution. Finally, the per-frequency per-baseline propagation visibility bias is estimated in Appendix B3.', 'B1 Deriving the per-antenna gain amplitude and phase error': "With our paired sky model approach, we can realize two equations using the measured visibility equation (equation 4) for applying an ' ideal ' sky model (equation 5) and a ' blended ' (equation 6) sky model, respectively. For each cross-correlation 𝑉 𝑖 𝑗 , we can determine the calibration error of the ill-calibrated scenario by combining equations (5) and (6) and solve  \nexp GLYPH<2> GLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> GLYPH<3> 𝑉 ideal 𝑖 𝑗 + 𝑛 𝑖 𝑗 = exp GLYPH<2> GLYPH<0> ' 𝜂 𝑖 + ' 𝜂 𝑗 GLYPH<1> + i GLYPH<0> ' 𝜙 𝑗 -' 𝜙 𝑖 GLYPH<1> GLYPH<3> ' 𝑉 blend 𝑖 𝑗 + 𝑛 𝑖 𝑗 . (B1)  \nHere, we effectively eliminated the impact of the baseline noise, since the noise item 𝑛 𝑖 𝑗 is specific to the baseline 𝐵 𝑖 𝑗 under the assumption made in equation (2c). By taking the logarithmic form of equation (B1), we have  \nGLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> + ln 𝑉 ideal 𝑖 𝑗 = GLYPH<0> ' 𝜂 𝑖 + ' 𝜂 𝑗 GLYPH<1> + i GLYPH<0> ' 𝜙 𝑗 -' 𝜙 𝑖 GLYPH<1> + ln ' 𝑉 blend 𝑖 𝑗 . (B2)  \nBy substituting the gain error ( Δ 𝜂 𝑖 ≡ 𝜂 𝑖 -' 𝜂 𝑖 and Δ 𝜙 𝑖 ≡ 𝜙 𝑖 -' 𝜙 𝑖 , respectively) and the visibility difference in terms of amplitude and  \nphase ( 𝑅 𝑖 𝑗 ≡ ln GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑉 blend 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -ln GLYPH<12> GLYPH<12> GLYPH<12> 𝑉 ideal 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> and 𝐼 𝑖 𝑗 ≡ arg GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑉 blend 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -arg GLYPH<12> GLYPH<12> GLYPH<12> 𝑉 ideal 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> , respectively), equation (B2) can be decoupled as two sets of linear equations presented in equation (7).", 'B2 Gain error solution': 'To solve the two sets of linear systems of equations (equation 7), we proceed by rewriting equation (7a) and equation (7b) in two matrix systems. Taking equation (7a) as an example, we have  \n<!-- image -->  \nwhere A 𝑅 , which will be referred to as the array configuration matrix, is a 𝑁 × 𝐶 2 𝑁 matrix, 𝒙 𝑅 is a 𝑁 × 1 matrix, and 𝒃 𝑅 is a 𝐶 2 𝑁 × 1 matrix. In general, there will not be enough degrees of freedom in 𝒙 𝑅 when solving 𝑁 unknown Δ 𝜂 with 𝐶 2 𝑁 measurements. Thus, the least-squares solution is solved and satisfies the optimization in equation (B3b). In a similar fashion, we can also write the set of the phase difference equations in a matrix form and solve the least-squares solution:  \nA 𝐼 𝒙 𝐼 = 𝒃 𝐼 , (B4a)  \nargmin ˜ 𝒙 𝐼 ∥ A 𝐼 ˜ 𝒙 𝐼 -𝒃 𝐼 ∥ 2 . (B4b)  \nTo solve those overdetermined systems of equations, we use the SVD (see Stewart 1993, for a review) to perform the pseudoinverse of the array configuration matrix and find the least-squares solution (see Brunton & Kutz 2019, for detailed demonstrations). According to the Moore-Penrose left pseudoinverse (Penrose 1955, 1956; Rohde 1965; Zlobec 1970), A matrix (referring to for both the amplitude and phase array configuration matrix) has a pseudoinverse A † , and the least-squares solution of the blending-induced amplitude and phase error, ˜ 𝒙 𝑅 and ˜ 𝒙 𝐼 , can be solved via  \n˜ 𝒙 = A † 𝒃 . (B5)', 'B3 Deriving the per-baseline propagation error': "The aftermath of the calibration against a sky model with sourceblending defects is further discussed here. As mentioned in Section 3.4.1, after solving each antenna gain factor, 𝑔 𝑖 , the calibration solution will be applied to the observed visibility to recover the true cross-correlation of the sky signal. In the ill-calibrated scenario, the blending-induced antenna gain errors, Δ 𝜂 s and Δ 𝜙 s, will propagate to the calibrated visibilities of the sky signal for each baseline and leave a per-frequency per-baseline propagation bias within the visibility space.  \nUsing the least-squares solution of both the amplitude and phase error matrix ( 𝒙 𝑅 and 𝒙 𝐼 , respectively) of the blending-induced antenna gain, we can further infer the per-frequency per-baseline propagation  \nbias upon the visibilities of the sky signal. To infer the poor-calibration impact on the SKA EoR experiments, a set of simulated data with instrumental responses and resolution fidelity, including extragalactic foregrounds, Galactic foregrounds, and the 21-cm signal, are utilized to estimate the propagation bias for each baseline in a realistic fashion. Here, taking Galactic foregrounds as an example, we can rewrite the measured visibility equation (equation 4) and its logarithmic form using the true visibility 𝑆 gal -true 𝑖 𝑗 and ill-calibrated visibility ' 𝑆 gal -ill 𝑖 𝑗 of Galactic emission as  \nexp GLYPH<2> GLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> GLYPH<3> 𝑆 gal -true 𝑖 𝑗 + 𝑛 𝑖 𝑗 = exp GLYPH<2> GLYPH<0> ' 𝜂 𝑖 + ' 𝜂 𝑗 GLYPH<1> + i GLYPH<0> ' 𝜙 𝑗 -' 𝜙 𝑖 GLYPH<1> GLYPH<3> ' 𝑆 gal -ill 𝑖 𝑗 + 𝑛 𝑖 𝑗 (B6)  \nand  \nGLYPH<0> 𝜂 𝑖 + 𝜂 𝑗 GLYPH<1> + i GLYPH<0> 𝜙 𝑗 -𝜙 𝑖 GLYPH<1> + ln 𝑆 gal -true 𝑖 𝑗 = GLYPH<0> ' 𝜂 𝑖 + ' 𝜂 𝑗 GLYPH<1> + i GLYPH<0> ' 𝜙 𝑗 -' 𝜙 𝑖 GLYPH<1> + ln ' 𝑆 gal -ill 𝑖 𝑗 , (B7)  \nrespectively. Similar to the decoupling of equation (B2), we further decompose equation (B7) into the real and imaginary parts separately and substitute the calibrated gain errors. Therefore, equation (B7) takes the form of  \nln GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 gal -ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -ln GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 gal -true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> = Δ 𝜂 𝑖 + Δ 𝜂 𝑗 , (B8a)  \narg GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 gal -ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -arg GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 gal -true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> = Δ 𝜙 𝑗 -Δ 𝜙 𝑖 . (B8b)  \nSince all the Δ 𝜂 s and Δ 𝜙 s were solved for each frequency individually using methods described in Appendix B2, the per-frequency perbaseline propagation amplitude and phase bias of Galactic emission can be inferred using equation (B8) by taking the exponential form of equation (B8a):  \nGLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 gal -ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 gal -true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> = GLYPH<2> exp GLYPH<0> Δ 𝜂 𝑖 + Δ 𝜂 𝑗 GLYPH<1> -1 GLYPH<3> GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 gal -true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> , (B9a)  \narg GLYPH<12> GLYPH<12> GLYPH<12> ' 𝑆 gal -ill 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> -arg GLYPH<12> GLYPH<12> GLYPH<12> 𝑆 gal -true 𝑖 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> = Δ 𝜙 𝑗 -Δ 𝜙 𝑖 . (B9b)  \nSimilarly, we also derive the amplitude and phase bias for the other components. Hence, the per-frequency per-baseline propagation bias of all the sky components in the visibility space can be estimated in the general form, as shown in equation (8).  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nTable B1. Calibration-related mathematical notions used in the paper."}