Uncovering the effects of array mutual coupling in 21-cm experiments with the SKA-Low radio telescope

ABSTRACT

We investigate the impact of mutual coupling (MC) between antennas on the time-delay power spectrum response of the core of the SKA-Low radio telescope. Using two in-house tools – Fast Array Simulation Tool (FAST; a fast full-wave electromagnetic solver) and Oxford Square Kilometre Array (OSKAR) (a GPU-accelerated radio telescope simulator) – we simulate station beams and compute visibilities for various station layouts (regular, sunflower, and random). Simulations are conducted in an epoch of reionization subband between 120–150 MHz, with a fine spectral resolution of 100 kHz, enabling the investigation of longer delays. Our results show that MC effects significantly increase foreground leakage into longer delays, especially for regular station layouts. For 21-cm science, foreground spill over into the 21-cm window extends beyond |$k_{\parallel } \sim 2\, \mathrm{ h}^{-1}$| Mpc for all station layouts and across all |$k_{\perp }$| modes, completely obscuring the detection window. We find that attempting to remove the foreground contribution from the visibilities using an approximated beam model – either based on the average embedded element pattern or by interpolating the embedded element patterns from a coarse channel rate of 781 kHz – results in residuals around |$\sim 10^{11}\,\mathrm{mK}^2\, \mathrm{ h}^{-3}\, \mathrm{Mpc}^3$|⁠. This is still seven orders of magnitude brighter than the expected level of the EoR signal (⁠|$\sim 10^{4}~\mathrm{mK}^2\, h^{-3}\, \mathrm{Mpc}^3$|⁠). We also find that station beam models with at least 4–5 significant digits in the far-field pattern and high spectral resolution are needed for effective foreground removal. Our research provides critical insights into the role of MC in SKA-Low experiments and highlights the computational challenges of fully integrating array patterns that account for MC effects into processing pipelines.

1 INTRODUCTION

Over the past decade, the detection of the faint cosmological 21-cm signal, arising from the hyperfine transition of neutral hydrogen atoms, has received increasing attention as one of the frontiers of observational radio astronomy. The signal promises a rich understanding of astrophysics, ranging from the earliest epochs of the Universe (the Dark Ages), through the first luminous sources (the Cosmic Dawn), and the eventual ionization of the intergalactic medium (epoch of reionization, or EoR). Measuring the spatial fluctuations of this signal would unlock a deeper understanding of the interplay between different heating and ionizing sources in the early Universe. The primary obstacle preventing such a detection is the isolation of the 21-cm signal from astrophysical foregrounds, such as Galactic synchrotron emission, radio galaxies, pulsars, and supernova remnants, which are up to five orders of magnitude brighter than the signal itself (Jelić et al. 2008). Foreground avoidance techniques, which exploit the rapid spectral variations of faint signals, enable their isolation at higher delays, helping to separate them from foreground emissions (Trott, Wayth & Tingay 2012; Parsons et al. 2012b; Liu, Parsons & Trott 2014a, b; Thyagarajan et al. 2015; Chapman et al. 2016). Alternatively, foreground removal techniques (Bowman, Morales & Hewitt 2006; Liu et al. 2009; Bonaldi & Brown 2015; Chapman et al. 2016; Pober et al. 2016) aim to subtract modelled foreground contributions from the observed data using knowledge of the telescope’s response and the celestial distribution of the foregrounds. In both foreground avoidance and removal approaches, precise models of the instrument’s beams and analogue chains are essential. Typically, two fundamental questions related to these models must be addressed:

• How accurately must we model the instrument’s response? This focuses on determining the necessary level of accuracy in matching the real instrument response with the model to ensure reliable data interpretation.

• In which regions of the observational domain are instrumental effects weakest? Identifying these regions will enable us to focus on areas of minimal instrumental uncertainty to preserve the data quality.

In addressing the accuracy criterion [question (i)], forecasts of the 21-cm power spectrum from CHIME (Bandura et al. 2014) indicate that achieving a 0.1 per cent accuracy in the primary beamwidth is necessary for reliable measurements (Shaw et al. 2015; Amiri et al. 2022). The REACH experiment (de Lera Acedo et al. 2022) has similarly determined that uncertainties in the antenna directivity must remain below |$-40$| dB to ensure robust global 21-cm detection (Cumner et al. 2024b). The MeerKAT Radio Telescope (Jonas & MeerKAT Team 2016) has evaluated the required pointing accuracy to use beam patterns with 1 per cent accuracy in power (Jonas & MeerKAT Team 2016; de Villiers & Cotton 2022). For SKA-Low (Dewdney et al. 2009), preliminary studies (Nasirudin et al. 2022) indicate that inaccuracies in beam knowledge – resulting from small positional deviations and up to 5 per cent of broken antennas – can produce residual errors as high as 1 per cent in the power spectrum EoR window and around 10 per cent in the foreground wedge. Similarly, for MWA (Tingay et al. 2013), 15–40 per cent of tiles typically have at least one broken dipole, leading to beam modelling errors that produce foreground residuals two orders of magnitude above the 21-cm signal after calibration (Joseph et al. 2019). With such error levels, systematic artefacts could well be misinterpreted as 21-cm structures, (i.e. Barry et al. 2016; Thyagarajan et al. 2016; Ewall-Wice et al. 2017).

For radio telescopes composed of densely packed antennas, addressing questions (i) and (ii) is complicated by a major issue: Mutual Coupling (MC) (Craeye & González-Ovejero 2011). In situ, strong electromagnetic interactions between antennas can induce sharp angular and spectral variations in the radiated far field, known as embedded element patterns (EEPs; Bird 2021; Bolli et al. 2022c; Anstey et al. 2024). Since each antenna operates in a distinct local environment within the array, the EEPs can vary significantly across elements in the station. In other words, MC effects are direction-dependent, frequency-dependent, and baseline-dependent. To date, MC is identified as a systematic source that potentially limits interferometric 21-cm power spectrum experiments such as the Low-Frequency Array (LOFAR) (van Haarlem et al. 2013; Mertens et al. 2020), the Hydrogen Epoch of Reionization Array (HERA) (DeBoer et al. 2017; Kern et al. 2019, 2020), the Precision Array for Probing the Epoch of Reionization (PAPER) (Parsons et al. 2010), and the Murchison Widefield Array (MWA) (Beardsley et al. 2016). LOFAR, NenuFAR (Zarka et al. 2012), and MWA have implemented full-Jones directional-dependent sky-based calibration algorithms SAGECAL (Kazemi et al. 2011; Yatawatta, Kazemi & Zaroubi 2012), DDECAL (van Diepen, Dijkema & Offringa 2018; Gan et al. 2023), and RTS (Mitchell et al. 2008). Incorporating realistic (simulated or measured) beam models into these calibration pipelines can help to reduce residual errors (Subrahmanyan 2021a; Chokshi et al. 2024). An alternative approach is mitigating MC effects using fringe rate filters. In this context, the HERA Collaboration modelled first-order dish–dish interactions with a semi-analytical formalism (Josaitis et al. 2022; Rath et al. 2024) to demonstrate the impact of MC at non-zero fringe rates. This informed the development of mitigation strategies and tailored fringe rate filters for foreground avoidance, building on Parsons et al. (2016).

The Square Kilometre Array Low (SKA-Low), examined in this paper, is an interferometric array under construction in Inyarrimanha Ilgari Bundara, Western Australia. It will consist of 512 stations, each containing 256 dual-polarization SKA Log-periodic antennas (SKALA) antennas within a 19 m radius, spanning |$\sim 74$| km. Initial studies have investigated the spectral smoothness of the bandpass for individual SKALA (de Lera Acedo et al. 2017) and estimated preliminary calibration requirements (Trott et al. 2017). These studies concluded that the isolated antenna response is sufficiently smooth to recover the EoR signal through statistical methods (Trott & Wayth 2016). However, as noted in de Lera Acedo et al. (2017), these analyses did not include the effects of MC between antennas within a station. The impact of MC has been a long-standing concern (de Lera Acedo 2012), with variations in EEPs observed both experimentally (Raucy et al. 2013; de Lera Acedo, Pienaar & Fagnoni 2018) and through simulations (Gueuning et al. 2015; de Lera Acedo et al. 2016). Recent UAV measurements (Subrahmanyan 2021b; Paonessa et al. 2023) and electromagnetic simulations (Davidson et al. 2020; Bolli et al. 2022b; Anstey et al. 2024) confirm that MC introduces notches at frequencies (55, 78, and 125 MHz) within the EoR band (50–200 MHz). A significant milestone has been the incorporation of full-wave EEP simulations into calibration algorithms, demonstrating a substantial reduction in visibility residuals (Subrahmanyan 2021a). These results underscore the importance of numerical simulations for understanding and mitigating the effects of MC on the delay power spectrum of the SKA-Low telescope. In this context, our analysis extends the earlier studies (de Lera Acedo et al. 2017; Trott et al. 2017) by explicitly addressing the MC effects.

Accurately modelling SKA-Low beams requires full-wave simulations of the electromagnetic response of each of the 256 antennas embedded in large (19 metre radius), dense, and irregular arrays across a wide bandwidth (50–350 MHz) with high spectral resolution (below 1 MHz). The complex antenna geometry, small spacing (0.01 wavelengths at the lowest frequency), and large station size (up to 45 wavelengths at the highest frequency) result in a challenging multiscale Computational Electromagnetics (CEM) problem. General-purpose solvers can take days or weeks per frequency point, even when accelerated with the Multi-Level Fast Multipole Method (MLFMM; Engheta et al. 1992) on large workstations. To overcome this, state-of-the-art CEM methods tailored for large, irregular arrays have been developed (Maaskant, Mittra & Tijhuis 2008; Bui-Van et al. 2018; Chose et al. 2023; Conradie et al. 2023). The in-house solver, named Fast Array Simulation Tool (FAST) (Gueuning et al. 2025), simulates all the antenna responses in the station in about 10 min per frequency point on a current standard laptop, with re-simulation for a different station layout taking just 1 min, enabling fine spectral sampling and therefore evaluation of impulse responses up to late time delays.

Modelling SKA-Low visibilities also presents computational challenges. Each of the 512 stations generates beamformed signals, which are then cross-correlated to produce over |$130\, 000$| visibility pairs. To ensure accuracy on SKA-Low’s longest baselines (up to |$\sim 74$| km), high-resolution sky maps with tens of millions of pixels are essential for avoiding aliasing in the sky integrals. Furthermore, each station’s beam, derived from precomputed EEP data, must be re-evaluated across all these sky directions. To address this, the tool Oxford Square Kilometre Array (OSKAR) has been developed Dulwich (2020) to leverage GPU parallelization, and we have also implemented a GPU-parallelized approach for station beam evaluation, originally proposed in Bui-Van et al. (2018).

In this paper, we use the in-house tools, FAST and OSKAR, to investigate the accuracy required for the EEP models [question (i)] as well as the impact of MC effects on the time-delay power spectrum of foreground emissions [question (ii)] through numerical simulations. To address question (i), we estimate delay power spectrum residuals after foreground removal, using an approximate station beam model that is corrupted by random noise, or based on the average embedded element pattern (AEP) (Wijnholds et al. 2019) or spectrally interpolated EEPs. Regarding question (ii), we investigate MC effects by comparing beam and visibility time-delay responses across three station layouts – regular, sunflower, and random. We quantify foreground leakage into longer delays caused by MC, considering both rotated and non-rotated stations. This paper also highlights the computational costs associated with integrating array patterns, which fully account for MC effects, in the SKA-Low pipeline. This paper is organized as follows. Section 2 defines the conventions for the coordinates systems used, while Section 3 presents the mathematical formulation for the forward model of these visibilities, based on EEPs, array patterns, and sky brightness distribution. In Section 4, we outline the parameters used for the telescope simulations. Section 5 describes our simulation pipeline, emphasizing the recent advancements introduced in the tools, FAST and OSKAR. Section 6 presents the simulation results, including an analysis of the station beam impulse response across different layouts, the foreground spill over induced by MC, and the residual levels after foreground removal using several approximate beam models. Finally, we summarize our findings in Section 7.

2 CONVENTIONS

Before outlining our mathematical derivations, we would like to clarify the notation used by the coordinate systems in this paper.

Comoving distances are spatial coordinates that measure separation in cosmology. Per their definition, these distances remain constant over time despite the expansion of the universe, as they account for the changing scale at different epochs. Assuming a Cartesian coordinate system with unit vectors (⁠|$\hat{\boldsymbol {x}},\hat{\boldsymbol {y}},\hat{\boldsymbol {z}}$|⁠), the comoving coordinates will be denoted as follows:

The spatial frequency coordinates, which are associated with the cosmological Fourier domain and referred to as the k-space, are given by

where |$\boldsymbol {k}$| is the wavevector with cosmological wavenumber |$k=||\boldsymbol {k}||$|⁠.

Physical relative distances, known as baselines, are represented by the baseline vector |$\boldsymbol {b}$|⁠:

when expressed in metres, or by |$\boldsymbol {b}_\lambda = \boldsymbol {b} /\lambda _0$| when expressed in units of the free-space wavelength |$\lambda _0$|⁠. The vector |$\boldsymbol {u}_\lambda = u/\lambda _0 \, \hat{\boldsymbol {x}}+ v/\lambda _0\, \hat{\boldsymbol {y}}$| represents the in-plane projection of the baseline vector |$\boldsymbol {b}_\lambda$|⁠. In this paper, the baseline notation |$\boldsymbol {b}$| refers to the separation between station centres, while |$\boldsymbol {b}_{ant}$| (in metres) denotes the position of an antenna with respect to its station centre.

The directional unit vector |$\hat{\boldsymbol {s}}$|⁠, which represents a line of sight, is given by

where |$n = (1-l^2-m^2)^{1/2}$|⁠, ensuring that |$\hat{\boldsymbol {s}}$| is constrained to unit-length.

3 A FORWARD MODEL OF TELESCOPE VISIBILITIES

A forward model simulates how signals from a celestial distribution are transformed by the instrument, capturing the distortions introduced at each stage of the system. This model is essential for understanding how instrumental effects propagate through the digital processing pipeline and manifest in the final product, the delay power spectrum. In this section, we focus on modelling the measured visibilities, specifically including the impact of MC between antennas in the beam response of each station.

3.1 System overview

The interferometer set-up comprises |$N_s$| phased array stations, each containing |$N_a$| antennas. After digitization and channelization of raw voltages from each antenna, the voltage from antenna n in station p is represented as a time-varying mono-frequency signal |$v_{\textrm {ant},np}(t, f)$|⁠. The |$N_a$| voltages in each station are then beamformed as |$v_{p}(f,t) = \sum _n w_{np}(f) v_{\textrm {ant},np}(f,t)$|⁠, with beamforming weights |$w_{np}(f)$| chosen to steer the beam or control its shape. Beamformed voltages are then cross-correlated to yield the complex visibility |$V_{pq}(f) = V(\boldsymbol {b}_{\lambda ,pq}, f)$|⁠, where |$\boldsymbol {b}_{\lambda ,pq}$| is the baseline vector linking the centres of stations p and q. A block diagram for a two-station, three-antenna configuration is shown in Fig. 1.

3.2 Mathematical model of the visibilities

The following derivations outline how the visibilities, measured by cross-correlating beamformed signals from different stations, are modelled using both celestial distribution and beam models.

Let us first assume a monochromatic incident plane wave propagating in the direction |$\hat{\boldsymbol {s}}$|⁠, with a time-varying electric field vector |$\boldsymbol {E}_{i}(f,t)$| representing its amplitude and polarization, where the subscript i denotes incident. Assuming and suppressing the time dependence |$e^{j2\pi ft}$|⁠, the electrical field of the incident plane wave is expressed as

where |$k_0$| is the free-space wavenumber. First, we consider single-polarized antennas (with a single feed port). According to Lorentz’s reciprocity theorem, the voltage |$v_{\textrm {ant},np}(f,t)$| appearing at the output port of a passive antenna n belonging to a station p illuminated by the plane wave can be obtained by treating antenna n as active while all the other antennas are passively terminated (Kildal 2015). Using the centre of antenna n as a phase reference in (5), this voltage is given by (Craeye & González-Ovejero 2011):

where |$\eta _0$| is the free-space characteristic impedance, |$Z_{L,np}$| is the input impedance of the amplifier connected to antenna n. The EEP, |$\mathbf {EEP}_{np}$| in equation (6), is expressed in volts and evaluated in the active case by feeding antenna n with a Thevenin equivalent of impedance |$Z_{L,np}$| and a source of unit voltage. In the passive case described in equation (6), the EEP is implicitly normalized by 1 volt, so that it becomes a dimensionless quantity. From there, the beamformed voltage associated with station p can be expressed as

where the array pattern (AP) is defined by

with the coefficients |$c_{pn} = w_{pn} Z_{L,np}/Z_0$|⁠. The phase reference is now taken at the centre of the station, introducing a phase factor in equation (8) based on the relative antenna distances |$\boldsymbol {b}_{\textrm {ant},np}$|⁠. The impedance |$Z_0$| in equation (7) has been introduced so that the AP shares the same units as the EEPs.

Considering now dual-polarized antennas with two feed ports, denoted as |$v_{p,X}$| and |$v_{p,Y}$|⁠, with indices X and Y referring to the orientation of the associated dipoles along |$\hat{\boldsymbol {x}}$|- or |$\hat{\boldsymbol {y}}$|-axis, we re-express the beamformed voltages using the |$2\times 2$| Jones matrix formalism (Jones 1941; Smirnov 2011). Omitting the explicit time (t) and frequency (f) dependence, this yields

where |$\boldsymbol {E}_i = E_{i,X}\, \hat{\boldsymbol {e}}_{X} + E_{i,Y}\, \hat{\boldsymbol {e}}_{Y}$|⁠, with |$\hat{\boldsymbol {e}}_X$| and |$\hat{\boldsymbol {e}}_Y$| representing the unit vectors corresponding to the horizontal and vertical components in the Ludwig-III polarization coordinate system (Ludwig 1973). In agreement with our previous derivation (equation 7), the elements of the Jones matrix are given by

where |$\mathbf {AP}_{p,X} = \textrm {AP}_{p,XX}\, \hat{\boldsymbol {e}}_{X} + \textrm {AP}_{p,XY}\, \hat{\boldsymbol {e}}_{Y}$| and |$\mathbf {AP}_{p,Y} = \textrm {AP}_{p,YX}\, \hat{\boldsymbol {e}}_{X} + \textrm {AP}_{p,YY}\, \hat{\boldsymbol {e}}_{Y}$| are the APs for feed ports X and Y, respectively. This formalism allows us to apply the radio interferometric measurement equation (RIME), as defined in Hamaker, Bregman & Sault (1996) and Smirnov (2011), to relate the Stokes I parameter of the measured visibilities |$V_{pq}$| to the spectral brightness distribution |$B(f,\hat{\boldsymbol {s}})$| and to the station beams:

with the instrument response:

where |$*$| denotes the complex conjugate. We refer to this function |$A_{pq}(f,\hat{\boldsymbol {s}})$| as the beam transfer function (BTF) while its time-delay Fourier counterpart |$\tilde{A}_{pq}(\tau ,\hat{\boldsymbol {s}})$| is called the beam impulse response (BIR). The factor |$\Omega _{pq}$| is called the integrated beam response and is given by

This factor ensures that the visibilities |$V_{pq}$| are expressed in Jy units when the brightness B is expressed in Jy sr⁻¹. To alleviate the notations, we will temporarily assume beams similar across stations, making the BTF effectively baseline independent, that is |$A_{pq} = A$|⁠.

3.3 Time-delay transform

As explained in Appendix  A, the power spectrum |$\tilde {\mathcal {P}}$| of an astronomical signal can be estimated from the amplitude of the time-delay Fourier transform |$\tilde{V}(\tau ,\boldsymbol {b})$|⁠, expressed in JyHz, of the visibilities |$V(f,\boldsymbol {b}_\lambda)$|⁠. In practice, the visibilities are only measured within a limited sub-band |$(f_{\min },f_{\max })$|⁠, a subset of the telescope’s full operational range. To minimize truncation effects, a window function |$W(f)$| is applied to the visibilities, resulting in

Substituting equation (11) into equation (14) and applying the convolution theorem, this expression becomes:

where |$\ast$| denotes a convolution product, |$c_0$| is the free-space speed of light, |$\tilde{W}(\tau)$| is the time-delay response of the window function W, and the spectrum |$\tilde{B}$| is the Fourier transform of the brigthness B. The spectral sampling step |$\Delta f$| is a critical factor in computing the delay Fourier transform (14). For time-limited functions over duration T, the Nyquist criterion requires |$\Delta f = \pi /(2T)$| to confine aliasing errors outside the observed delay range |$\tau \in [0, T]$|⁠.

3.4 Instrumental effects in the baseline-delay domain

Time-delay visibilities are analysed in baseline length-delay space, |$(b,\tau)$|⁠, by grouping visibilities into baseline length bins and averaging within each bin to obtain the average visibility, |$\tilde{V}_{av}(b,\tau)$|⁠. For simplicity, we consider the case of a single point source denoted by index i of flux density |$S(f) = B(f,\hat{\boldsymbol {s}}) \Delta l_i \Delta m_i/n_i$| in direction |$\hat{\boldsymbol {s}}_i$| and solid angle |$\Delta l_i \Delta m_i/n_i$|⁠. The relation (15) thus reduces to

The effect of the instrument appears as a cascade of convolution products in the time-delay domain, which we will analyse sequentially. All these effects are summarized and sketched in Fig. 2.

First, the baseline exponential in equation (11) causes a translation of sources flux density spectrum, denoted as |$\tilde{S}$|⁠, by a time delay that scales with the projection (Parsons et al. 2012b; Chapman et al. 2016),

This relation indicates that, for planar baselines with |$\boldsymbol {b}=\boldsymbol {u}$|⁠, sources positioned at the zenith introduce no delay, while sources oriented in line with or opposite to the baseline direction yield either positive or negative delays. The average visibility |$\bar {V}_{av}(b,\tau)$| peaks around the delay |$\tau = b \sin \theta _i /c_0$| with the source zenith angle |$\theta _i$|⁠, thus following a line with slope |$\sin \theta _i/c_0$| in baseline length-delay space (b, |$\tau$|⁠). This delay reaches its maximum, the horizon limit, at |$\tau =b/c_0$|⁠, when the source is on the horizon (Morales et al. 2012). As the baseline length b increases, the separation between spectrum peaks at different elevation angles |$\theta _i$| increases, enhancing resolution in the time-delay domain. For sources with spatially smooth brightness distribution B (diffuse emission), the visibility |$\bar {V}_{av}(b,\tau)$| decreases rapidly with baseline length b and increases with delay |$\tau$| peaking at the horizon limit.

Secondly, the shifted source spectrum |$\tilde{S}$| is convolved with the BIR, |$\tilde{A}(\tau , \hat{\boldsymbol {s}}_i)$|⁠, in the direction |$\hat{\boldsymbol {s}}_i$|⁠. Sources within the main lobe thus appear brighter than those in the sidelobes, and higher sidelobe levels manifest in the |$(b, \tau)$| plots as increased peak values around |$\tau = \sin \theta _i b/c_0$|⁠. The beam limit, defined by |$\tau = b \sin \theta _i/c_0$|⁠, is based on the half-beamwidth |$\theta _i$| at the centre frequency. The BIR can decay slowly with delay |$\tau$|⁠, and when convolved with a more rapidly decaying signal spectrum |$\tilde{S}$|⁠, such as foreground emission, it spreads the signal power into longer delays, causing leakage that can extend far beyond the horizon limit, as shown in Section 6. This slow decay arises from unwanted electromagnetic interactions, such as multiple reflections or scattering, which prolong wave attenuation. Optimizing antenna geometries and array configurations during design can help minimize the BIR decay rate.

Finally, |$\tilde{W}(\tau)$| in equation (16) is the Fourier transform of the frequency window function |$W(f)$|⁠, used to reduce ringing artefacts from truncating the frequency band to |$(f_{\min }, f_{\max })$|⁠. A window with a narrow first lobe and high dynamic range is preferred (Lanman et al. 2020); here, we use the convolved Blackman–Harris window as in Parsons et al. (2014). The first lobe of the combined response |$\tilde{W} \ast \tilde{A}$| matches or exceeds that of |$\tilde{W}$|⁠, depending on beam chromaticity. However, the dynamic range is generally governed by the slower decay of the beam response |$\tilde{A}$|⁠.

3.5 Delay power spectrum

After applying the time-delay Fourier transform to the visibilities |$V_{pq}$|⁠, yielding the time-delay visibilities |$\tilde{V}_{pq}$|⁠, we can then evaluate the delay power spectrum from these time-delay visibilities using the following formula (Morales & Hewitt 2004):

where z is the redshift, |$X(z)$| and |$Y(z)$| are cosmological scaling factors, and |$\Omega _A$| is a normalization factor that depends on the instrument response A. The delay power spectrum is often represented as a function of the parallel and perpendicular components, |$\boldsymbol {k}_\perp$| and |$\boldsymbol {k}_\parallel$|⁠, of the cosmological wavevector |$\boldsymbol {k}$| defined in equation (2). The link to baseline time-delay coordinates is given by

For interested readers, a detailed reminder about the derivation of formulas (18) and (19), along with an expansion of the various constants, is provided in Appendix  A.

4 NUMERICAL EXPERIMENT: SKA-LOW

The spiral configuration of SKA-Low features expanding arms with baseline lengths up to |$\sim 74$| km for higher resolution, alongside a core of 224 densely packed stations that provide instantaneous uv-coverage with minimum baseline lengths starting from around |$b_{\min } = 40$| m to around |$b_{\max } = 1000$| m. One of SKA-Low’s key scientific goals is the detection of the 21-cm signal. In this section, we outline the parameters of a realistic scenario considered for the detection of this signal, with results to be presented in Section 6.

4.1 Observation parameters

The SKA Epoch of Reionization Science Working Group has declared their intent to target the 21-cm signal between the redshift of |$6\lt $| z |$\lt 30$| (Braun et al. 2015) and therefore an approximate frequency range of |$46~$| MHz < f <|$202~$| MHz. In this work, we focus on a sub-band from |$f_{\text{min}}=120$| MHz to |$f_{\text{max}}=150$| MHz at a 100 kHz spectral resolution (⁠|$\Delta f$|⁠) placing us firmly within the 21-cm band while allowing us to highlight the impact of MC-induced features (Davidson et al. 2020; Anstey et al. 2024) present in the beam. To minimize the ratio between the 21-cm signal and astrophysical foregrounds, observation times are typically chosen to ensure that the beam’s primary lobe and the field of view are relatively free from bright foreground emissions. This involves avoiding areas such as the Galactic centre (GC) and prominent point sources like Rigel. An example of a suitable observation field is the MWA EOR1, located at |$\mathrm{RA} = 60^{\circ }$| and |$\mathrm{Dec.} = -27^{\circ }$| (see Bowman et al. 2013; Beardsley et al. 2016). For our simulations, we select two observation times. The first is at 12:31:30 UTC on 2000 January 1, when the Galactic plane is far from the SKA-Low field of view and located along the horizon. The second is at 09:31:30 UTC on January 1, 2000, when the GC is below the horizon, minimizing any leakage through the far-out sidelobes. The sky map at 135 MHz for both observation times is shown in Fig. 3 and a description of the sky model content is outlined in Section 4.3.

4.2 Telescope model

4.2.1 Antenna geometry

The 4th version of the SKA element, illustrated in Fig. 4, is a dual-polarized log-periodic antenna featuring 16 dipoles per arm (SKALA4, (de Lera Acedo et al. 2018). It operates over a |$7:1$| frequency range, spanning from 50 to 350 MHz. Each arm is differentially amplified using a low-noise amplifier (LNA) with an input impedance of 100 Ohms.

4.2.2 Station layout

Fig. 5 illustrates the three different layouts used in this numerical experiment. Each layout is arranged in a circular footprint with a radius of 19 m and a minimum distance of 1.7 m between the antennas.

• regular: A regularly spaced grid of elements within the circular station footprint, each element positioned 2.14 m apart.

• sunflower: A sunflower-head array, with the position of each element given in Vogel (1979) and evaluated for SKA-Low in Bolli et al. (2022a).

• random: A random layout which results from random perturbation of the sunflower layout (Anstey et al. 2024).

4.2.3 Array layout

The coordinates of the array centre and the projected spacing between these stations were obtained from Revision 3 of the SKA-Low configuration and illustrated in Fig. 6. We exclude all the long baselines and limit our simulations to the 224 core stations of SKA-Low, to ensure that the currently available sky model meets the Nyquist resolution. Note also that during our analysis, the baselines are irregularly binned based on a logarithmic length distribution specific to SKA-Low.

In the following, we consider the case where all stations in the core have an identical, un-rotated station layout. Additionally, we examine a scenario where each station (antenna orientation and position relative to the station centre) is rotated by a unique angle, as illustrated in Fig. 6. The rotation angles are assigned randomly to each station, following a uniform distribution between 0 and|$2\pi$|⁠.

4.3 Foreground model

To accurately model the localized and extended foreground emission accessible to the SKA-Low, a composite sky model is created by treating the extended emission as a set of unresolved point sources. The sky model components include:

• A diffuse Galactic foreground emission from the desourced and destriped Haslam 408 MHz survey, with |$\sim 5 \times 10 ^{7}$| pixels, equating to a 1.718 arcmin resolution (Remazeilles et al. 2015).

• A collection of point sources¹ from GaLactic and Extragalactic All-Sky MWA Survey (GLEAM), with |$\sim 3\times 10 ^{5}$| sources (Wayth et al. 2015).

At low frequencies, radio foreground observations may be approximated by a near-featureless power law across the frequency spectrum. The composite map can be frequency-scaled within OSKAR as needed using the equation:

In this equation, |$f_0$| denotes the reference frequency, and |$\beta$| signifies the spectral index. For the Haslam 408 MHz survey, |$\beta =-0.7$|⁠, whereas for the GLEAM survey, a distinct |$\beta$| is provided for each source.

5 ACCELERATED SIMULATION TOOLS

We present the simulation pipeline used to generate the delay power spectrum for large interferometric arrays. Accurately modelling visibilities for telescopes with many antennas and stations poses computational challenges, as it requires full-wave simulations of each antenna’s electromagnetic response, including MC effects, and the evaluation of APs over multiple directions and frequencies. This pipeline leverages two in-house tools: FAST for efficient electromagnetic simulations and the radio-telescope simulator OSKAR for parallelized interferometric calculations.

The pipeline, outlined in Fig. 7, begins with electromagnetic simulations for the |$2N_a$| EEPs in a frequency band |$(f_{\min },f_{\max })$| at resolution |$\Delta f$|⁠. Next, using sky catalogues and beam data from FAST, the OSKAR simulator computes the |$N_s(N_s-1)/2$| visibilities for |$N_s$| stations. During OSKAR visibility simulations, the beam computations are managed by the harp_beam library (using FAST-generated beam data). Once all visibilities |$V_{pq}(f)$| are generated, they are binned logarithmically by baseline length b and averaged. We then apply a window function and perform the time-delay Fourier transform to obtain the time-delay visibilities |$\tilde{V}_{av}(\tau ,b)$|⁠. Computation times for both simulation tools will be based on the simulation parameters provided in Section 4.

5.1 FAST

The recent advancements of FAST, (Gueuning et al. 2025), an accelerated direct solver based on the method of moments (MoM), have significantly reduced the simulation time required for the computation of all the EEPs of large, dense, and irregular arrays of identical but complex antennas. In the case of SKA-Low, it provides all the EEPs for a frequency point in |$\sim 10$| min on a standard current laptop (see Table 1 which details this benchmark). This represents a significant resource and time-saving improvement, as it operates approximately 25 times faster than FEKO’s Multilevel Fast Multipole Method on a 128-core workstation. Additionally, re-simulation for different station layouts with FAST takes under a minute, allowing for the optimization of station layout and analysis of MC effects at a fine resolution. The electromagnetic solver has been validated against several other numerical solvers, including CST,² WIPL³ (Bui-Van et al. 2018), and FEKO⁴ (Gueuning et al. 2022, 2025; Cavillot et al. 2024, showing agreement to within 2–3 digits of accuracy in the far-field pattern, provided the same mesh is used. Additionally, we have carefully validated our simulations against FEKO using smaller test-case stations to ensure that the chromaticity observed in our simulations results from MC and not from approximation errors related to the solver.

5.2 OSKAR

OSKAR is a GPU-accelerated telescope simulator specifically designed for the SKA (Dulwich 2020), utilizing the Radio Interferometric Measurement Equation (RIME). The RIME framework models the contributions of the complex interferometric visibility outlined in equation (11) through a discrete sum across all visible sources in the sky whilst accounting for any potential instrumental effects including MC. To obtain the visibilities, OSKAR requires a star catalogue to describe all the source properties and a description of the position of each station, known as the sky model and telescope model respectively. For SKA-Low, this tool allows us to obtain the visibilities in just 4.3 s per frequency for |$\sim 12.9\times 10^6$| sources (Table 2 which details this benchmark). The harp_beam library is integrated into OSKAR and efficiently calculates station beams directly from each antenna’s current distribution computed from FAST. On an NVIDIA A100 GPU, it computes SKA-Low station beams for a 16 000-pixel grid in 1 ms. Pre-computed APs save time for identical layouts, but differing layouts require re-computation. In such cases, OSKAR simulations show station beam evaluations dominate runtime, being nearly 1000 times slower than visibility integrations.

6 NUMERICAL RESULTS

We present the numerical results obtained using our pipeline, as outlined in Section 5, and the simulation parameters described in Section 4. First, we examine the variation in the BIR fall-off rate for three different station layouts, highlighting the impact of intra-station MC. Next, we explore a foreground avoidance scenario with the SKA telescope and discuss the size of the available detection window. Finally, we consider a foreground removal scenario in which the station beam is inaccurately modelled using either the AEP approximation, corrupted or interpolated EEPs. We then attempt to remove the foreground power using these approximate beam models.

6.1 Beam impulse responses

The station beams are plotted at 135 MHz the band centre, in Fig. 8 for the three different layouts illustrated in Fig. 5: regular, sunflower, and random. As can be expected (Clavier et al. 2014), at this frequency, all three layouts exhibit a similar main beam value, around 45 dBV, along with comparable levels for the first three sidelobes. The beam pattern for the regular layout exhibits large grating lobes (nearly 30 dBV). The beam for the random layout shows higher secondary sidelobes compared to those of both the sunflower and regular layouts.

Fig. 9 illustrates, in grey, all 256 individual EEPs for feed-H, evaluated at zenith and sampled at a 1 MHz rate, along with their corresponding BIRs. The solid black line represents the average embedded element pattern (AEP) across the elements, while the dashed line represents the isolated element pattern (IEP) of SKALA4. On the one hand, the IEP appears smooth across the frequency band, suggesting a nominal operation of the SKALA4 antenna. As a consequence, the BIR of the isolated antenna decays quickly reaching |$-50$| dB at 300 ns, with the main lobe width matching that of the window function |$\tilde{W}$|⁠. On the other hand, prominent notches in the EEPs can be observed at 124 and 132 MHz in the frequency response of both the regular and sunflower layouts. The drop at 124 MHz results from destructive interference between the field radiated by currents on the active antenna and the field radiated by currents induced on the passive antennas due to MC. This effect occurs when the inter-element distance approaches one wavelength (Cumner et al. 2024a). The depth and width of this notch have been observed to increase with the array regularity, i.e. the redundancy in baseline length distribution (Anstey et al. 2024). In the delay domain, these MC effects cause the first lobe of the BIRs of the EEPs to be wider for both the regular and sunflower layouts (extending to 200 ns) compared to the first lobe width of the random layout (which extends to 150 ns). The slower decay in the tail of the BIRs, at delays above 215 ns, suggests the presence of sharp, small-scale spectral variations in the BTFs that vary at a quick rate (faster than 1 MHz). These variations appear to be less pronounced at the AEP level (or equivalently here in the AP) in the sunflower and random layouts. Specifically, the BIR of the AEP is around |$-20$| dB for the regular layout, while it drops to |$-30/40$| dB in the sunflower and random layouts. This attenuation is likely due to the irregular arrangement of antennas, which may help minimise MC effects at these scales.

We analyse now the BTF and BIR of the AP for the random layout, sampled at 100 kHz, and examine the impact of cubic interpolation from a coarse spectral grid, with steps of |$\Delta f = {781}$| kHz, to match the expected SKA-Low station channel width (Trott & Wayth 2016). As shown in Fig. 10 (left), the zoomed region highlights rapid variations in the BTF that outpace the 781 kHz rate, resulting in interpolation errors of up to 0 dBV in the AP (roughly 2 significant digits in the E-field pattern, respectively). In Fig. 10 (right), this results in time-delay residual errors around 10 dB, leading to 3 per cent errors at short delays, and divergence from exact values at approximately 700 ns. These interpolated EEPs will be applied in our pipeline for the foreground removal scenario in Section 6.3 to evaluate the impact of interpolation errors on the delay power spectrum.

6.2 Foreground avoidance

Early studies suggest that the sensitivity of the SKA-Low telescopeshould be sufficient to detect the 21-cm power spectrum at scales defined by |$k \in [0.1,0.33]$|  cMpc|$^{-1}$| (Cohen, Fialkov & Barkana 2018; Barkana et al. 2023). However, this section illustrates, that MC effects may compromise these predictions.

To characterize foreground contamination, we define the delay spread, denoted as |$\tau _{\mathrm{ min}}(b)$|⁠, also known as the foreground spill over. This value represents the smallest delay at which the 21-cm signal power spectrum |$\tilde {\mathcal {P}}_{21}$| surpasses the foreground power spectrum |$\tilde {\mathcal {P}}_\mathrm{ F}$|⁠:

In our analysis, we calculate a fiducial 21-cm signal power spectrum amplitude (Lanman et al. 2020) of |$\tilde {\mathcal {P}}_{21} \approx 114^2~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| at the band’s central redshift |$z=9.52$|⁠.⁵ This was derived using 1.5 cGpc 21-cm brightness temperature cubes, generated with the semi-numerical code 21cmSPACE (Visbal et al. 2012; Fialkov et al. 2013).⁶

Using the telescope model, Fig. 11 shows the delay power spectrum for an observation at 09:31:30 UTC on 2000 January 1, with the GC along the horizon. The simulation spans 300 channels from 120 to 150 MHz, using an identical (un-rotated) random layout for all 224 stations of the SKA-Low core. The results are presented in two cases. In the first case, shown in the top plot (without MC), the AP is approximated using the IEP of SKALA4 multiplied by the array factor (AF). The AF is given by the equation

where unit coefficients |$c_{np}=1$| are assumed. In the second case, shown in the bottom plot (with MC), the exact AP is used, which incorporates all EEPs as defined in equation (8). The dashed white line marks the delay spread, |$\tau _{\min }$|⁠, outlining a 21-cm detection window shown as the upper left black triangular region in the top panel. In this region, the foreground signals are sufficiently attenuated, falling below the estimated level of the 21-cm signal. The solid black line indicates the horizon limit, defined by the maximum geometric delay (equation 17) for a source at the horizon. A bright band of power, present across all baseline lengths and confined to short delays, is identified as the intrinsic foreground. Additionally, along the horizon limit, a bright limb is observed, caused by an apparent increase in flux density (Thyagarajan et al. 2015). The dashed black line indicates the beam limit, determined by the primary beamwidth of |$4.2^{\circ }$| at the central frequency, as illustrated in Fig. 8. The slower decay rate of the BIR of EEPs, compared to that of the IEP, observed earlier in Fig. 9 translated into similar effects in the delay power spectrum plot shown in Fig. 11. The delay spread for the smallest baseline (⁠|$\sim 40$| m) generated with the AF approximation extends to |$700~$| ns, while it extends beyond the maximum delay for the exact AP. The chromaticity introduced by MC significantly increases the delay spread across all baselines, causing the foreground emission to leak and completely occlude the 21-cm detection window. As validation of our absolute power spectrum levels, we compare the peak delay power spectrum without MC effects (top plot in Fig. 11) against results from PAPER (Pober et al. 2013), MWA (Beardsley et al. 2016), HERA (Abdurashidova et al. 2023), and SKA-Low simulation (Trott & Wayth 2016). The agreement confirms consistent diffuse emission delay spreads and comparable peak values (around |$2 \ 10^{14}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$|⁠), despite differences in bandwidth and baseline binning. We will next investigate how observation times, station layouts, and random station rotations influence the spill over.

Using the same telescope model, Fig. 12 compares the delay power spectrum for two observations on 2000 January 1: at 09:31:30 UTC, when the Galactic plane is along the horizon, and at 12:31:30 UTC, when the GC is below the horizon. We decompose the delay power spectrum into separate contributions from point sources (using the GLEAM catalogue) and diffuse emission (using the HASLAM map). The results show that spill over is larger for diffuse emission (⁠|$\sim 10^{10}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| at 2000 ns) compared to point sources (⁠|$\sim 10^{6}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| at 2000 ns). This difference arises due to the increased chromaticity of the diffuse emission, resulting in an apparent broader delay spectrum |$S(\tau)$|⁠. The diffuse emission contribution also tapers off as baseline length increases. Fig. 12 clearly shows that, at both observation times, diffuse emission is concentrated around the horizon line, while point source emission predominantly originates from emissions in the main beam. Close inspection of the left-hand panels (GLEAM) exhibits faint lines (an example is denoted by a black arrow) of increased power running parallel to the horizon line. Each line originates from an increase in sensitivity due to the fourth sidelobe of the random layout appearing in Fig. 8. Furthermore, when the GC is below the horizon (bottom left plot), a notable increase in power is observed relative to when it is above the horizon (top left plot), extending to the beam limit. This increase is attributed to a higher source count and flux density within the primary beamwidth, as recorded by the GLEAM catalogue. In the central panels showing diffuse emission at both observation times, a decrease in diffuse power along the horizon is observed when the GC is below the horizon, compared to when it is at grazing angles, which explains the increased foreground spill over when the GC is up.

In Fig. 13, we analyse how the power spectrum varies with identical station layout (regular, sunflower, and random, as shown in Fig. 5) for the first observation time with the GC above the horizon. Given that the main beam and first sidelobe levels are similar across layouts (see Fig. 9), power levels remain consistent up to and slightly beyond the beam limit for both point source and diffuse emission. As the elevation angle increases to |$20^\circ$|⁠, the larger sidelobes in the random layout compared to those in the sunflower or regular layouts result in increased point source power levels within the foreground wedge, between the horizon and beam limit, and with a more prominent sidelobe line corresponding to the contribution of bright point sources. The grating lobe of the regular layout appears as a brighter region near the horizon line. Regarding MC-specific effects, in line with the discussions in Section 6.1, the regular and sunflower layouts exhibit a broader intrinsic foreground region, extending to approximately 200 ns. In contrast, the random layout, with its narrower first lobe width, shows a response extending only to around 150 ns. The delay spread is noticeably much larger for the regular layout for both point source and diffuse emission. The leakage of diffuse emission appears similar for the sunflower and random layouts. The spill over of point source emission is however slightly lower for the sunflower layout (around |$10^9~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| at 800 ns) compared to the random layout (around |$10^{10}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| at 800 ns). Overall, for all layouts, the entirety of the detection window (EoR window), up to |$k_{\parallel } \sim 2~\mathrm{h}^{-1} \, \mathrm{Mpc}$| (see bottom row of Fig. 13) is contaminated.

Each of the 224 SKA-Low core stations is now assigned a unique random rotation angle, uniformly distributed between |$0^\circ$| and |$360^\circ$| (see Fig. 6). This rotation is expected to attenuate contributions from sources in the sidelobes for station layouts with more irregular azimuthal distributions. Practically, this means that each station beam becomes unique, requiring the computation of 224 individual station beams instead of just one, which significantly increases computation time in OSKAR, as shown in Table 2. Fig. 14 compares the delay power spectrum slices for un-rotated and rotated random layouts at given baseline lengths |$b = [59,~198,~882]$| m. While the size of the intrinsic foreground region remains unchanged, the rotated case shows a roughly twofold decrease in power beyond the beam limit (sidelobe region), particularly along the brightened horizon limb. However, leakage at longer delays remains unaffected.

6.3 Foreground removal

The EoR program for the SKA-Low telescope (Koopmans et al. 2015) comprises three major experiments: power spectrum analysis, which uses visibilities averaged to around 100 kHz and 5 s; tomography, relying on imaging data formed from these visibilities; and the 21-cm Forest experiment, which employs high spectral resolution image cubes with finer spectral channels of |$\sim 4$| kHz. Several studies (Chapman et al. 2016; Sims et al. 2016; Li et al. 2019) suggest that a combination of removal and avoidance approaches will likely be employed for 21-cm science with SKA-Low. As previously discussed, the foreground removal approach aims to eliminate the foreground contribution from measured visibilities, leaving uncontaminated 21-cm signals. Achieving accurate foreground removal requires a precise model of the instrumental response on the sky. To estimate foreground residuals, we conduct simple removal tests using approximate station beam models. The residuals in the visibilities are computed as:

where |$\tilde {V}_{pq,~\mathrm{app}}(f)$| represents the visibility computed with the approximate beam models. The power spectra associated with these residuals are then computed using equations (14) and (A7). In this study, we consider three approximations to the beam model: station beams corrupted by random noise, the AEP approximation, and the spectral interpolation of EEPs from coarser channels. These tests allow us to assess the impact of beam model inaccuracies on foreground removal and the resulting residuals in the delay power spectrum. In the following examples, we consider un-rotated station layouts.

The first test consists of adding Gaussian white noise to the APs as follows:

where |$\mathcal {N}_q$| represents complex-valued Gaussian noise with zero mean and variance |$\sigma ^2$|⁠. This noise is a multivariable function of frequency f, the |$N_{\mathrm{ pix}}$| pixel directions |$\hat{\boldsymbol {s}}$|⁠, and the station indices p. We assume four different values for the variances, |$\sigma ^2 = {10^{0}, 10^{-2}, 10^{-4}, 10^{-6}}$|⁠. For a signal coming from zenith, where the beam value is around 45 dBV, this corresponds to signal-to-noise ratios |$\textrm {SNR} = {45, 65, 85, 105}$| dB. Since 20 dB corresponds to a single digit of accuracy in the E-field, the noise is approximately 2, 3, 4, and 5 digits below the received voltage for the respective SNR values. Fig. 15 illustrates a selection of delay power spectrum residuals across baseline bins with lengths of |$b = [59,~198,~882]$| m. Each additional significant figure of accuracy in the station model roughly corresponds to a reduction of two additional digits in the power spectrum residuals at short delays, with values around |$(10^9, 10^7, 10^5, 10^3)$| mK|$^2\, h^{-3}$| Mpc|$^3$| for the smallest baseline. Thus, to achieve effective foreground removal below our fiducial 21-cm signal, a 4–5-digit accuracy in the station beam model is required to keep the foreground residuals below the |$114^2~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| level of the 21-cm signal. Since the AP is the sum of the EEPs (see equation 8), the noise variance on the EEPs that results in similar noise on the AP is |$\sigma ^2/256$|⁠. We estimate the required accuracy on the EEPs to be around 3–4 significant digits.

Then, we analyse the residuals resulting from approximating the AP with the AEP multiplied by the AF (Wijnholds et al. 2019) as follows:

Fig. 16 shows slices of the delay power spectrum for |$b = [59,~198,~882]$| m, comparing the power spectrum using the exact AP, the AEP approximation, and the residual after attempting foreground removal in the visibilities. The AEP approximation is accurate within the intrinsic foreground region (up to 150 ns) with power spectrum residuals on the order of 1 per cent to 0.1 per cent at short delays. The AEP approximation is off by seven orders of magnitude, preventing the residuals from reaching the level of our fiducial 21-cm signal. As shown in Fig. 15, the residuals from the AEP approximation are one order of magnitude higher than those obtained by using station beams only accurate to 2 digits (assuming uncorrelated noise errors). Note that the tail of the estimated power spectrum decays faster than the true delay spectrum. This observation is consistent with the lower tail of the BIR associated with the AEP in Fig. 10.

In the final test, we examine the residuals resulting from interpolating the EEPs from a coarser spectral grid with a sampling rate of |$\Delta f = 781$| kHz to the fine sample rate of 100 kHz as discussed in Section 6.1. This coarse rate corresponds to the separation between the 384 coarse channels of SKA-Low, covering the 300 MHz bandwidth, which is formed at each station before applying beamforming and/or calibration weights (Comoretto et al. 2020). Each coarse channel contains around 200 fine 4 kHz channels. The calibration strategy outlined in Trott & Wayth (2016) proposes using a low-order polynomial fit for the calibration coefficients across three adjacent coarse channels. Inspired by this work, we interpolate the EEPs, simulated at the centre of each coarse channel within the |$120\!-\!150$| MHz band of interest, and then upsample these EEPs to a 100 kHz resolution. The AP is then formed using equation (8) from these interpolated EEPs, resulting in an approximate beam model. We have plotted the residuals for the three baseline lengths in Fig. 17. In agreement with the 30 dB residual observed in the BIR of the AP (Fig. 10), the power spectrum residuals are around 1 per cent. As shown in Fig. 17, the coarse interpolation of the EEPs leads to an estimated delay power spectrum that decays faster than the true spectrum at long delays due to spectral smoothing. Considering the results of previous tests, another conclusion from this experiment is that reducing residuals by employing different EEP models for each antenna is effective only when these models are computed at high spectral resolutions.

7 CONCLUSION

This paper analyses the impact of Mutual Coupling (MC) between antennas on the time-delay power spectrum of the SKA-Low telescope. Electromagnetic simulations reveal significant and rapid spectral variations in the station beam directivity, driven by energy travelling and resonating between antennas due to MC. For 21-cm science, these MC effects cause substantial foreground leakage into the detection window, compromising foreground avoidance techniques. Highly accurate station beam models (within 4–5 orders of magnitude) are required to effectively remove foregrounds. Brute-force approximations, such as these based on randomizing MC effects or on interpolating coarsely beam responses, will likely fall short. Although this study focuses on the direct evaluation of time-delay power spectrum, we expect image-based analyses, (see for instance Morales et al. 2019), to require similar beam model accuracy for image cleaning before power spectrum reconstruction. New calibration algorithms and measurement techniques will have to be developed to account for MC effects and enable precise beam modelling.

The analysis in this paper is made possible through the development of two key numerical tools:

• The FAST electromagnetic solver: All EEPs were simulated across 300 frequency points, completing in roughly two days on a standard laptop, with only a few additional hours needed for re-simulation in cases of new station layouts.

• The OSKAR visibility simulator: The simulator processes the 24 976 visibilities in the SKA core across all frequencies in just 10 min using 4 A100 GPUs and including diffuse emissions and point sources across approximately 10 million sky pixels. The beam evaluation method embedded into OSKAR directly calculates station beams from the aperture current distribution, achieving evaluation times around 2 s per station beam for the configuration above.

Using these simulation tools, we fully integrated exact array patterns (APs) directly into visibility calculations, enabling precise tracking of MC effects in the time-delay power spectrum analysis. This yielded the following key findings regarding intra-station MC impacts:

• Delay spread broadening: MC effects significantly extend the delay spread of the power spectrum. For EoR science, this can result in the complete obscuration of the detection window.

• Layout-dependent MC effects: while offering beneficial beam redundancy, regular station layouts are more prone to resonant wave effects that lead to sudden power drops across the field of view. In contrast, randomized layouts reduce these issues, providing more stable beam behaviour and a narrower delay spread in the impulse response.

• Residuals with AEP approximation: When approximating exact APs by using the AEP multiplied by the array factor, power spectrum residuals reach peak values around |$10^{11}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$| within the foreground wedge. This is about seven orders of magnitude higher than the expected level of the 21-cm power spectrum (⁠|$\sim 10^{4}~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$|⁠).

• Residuals with corrupted station beams: Testing various levels of accuracy (2, 3, 4, and 5 significant digits) resulted in residual power spectrum values of |$10^9$|⁠, |$10^7$|⁠, |$10^5$|⁠, and |$10^3~\mathrm{mK}^2\, h^{-3} \, \mathrm{Mpc}^3$|⁠, respectively. This demonstrates the need to model the station beams to |$4\!-\!5$| significant digits.

• Residuals with interpolated EEPs: Interpolating the EEPs from a coarse (781 kHz) to a finer (100 kHz) spectral grid before forming the AP results in power spectrum residuals around 1 per cent at short delays, highlighting the need for high spectral resolution in EEP models to reduce interpolation errors.

This work highlights the need to further develop station beam evaluation methods in processing pipelines, incorporate additional electromagnetic effects, and assess their impact across observational domains (baseline, delay, fringe rate, and image) to filter corrupted data and refine calibration algorithms. Future work should examine whether certain parts of the band are less affected by MC and extend these investigations to include image-based methodologies.

ACKNOWLEDGEMENTS

The first two authors, Oscar S.D. O’Hara and Quentin Gueuning, contributed equally to this work and should be considered joint first authors.

The authors thank Dr. Karel Adamek for his assistance in parallelizing the HARP beam library, Dr. Vladislav Stolyarov for benchmarking it. We also thank Pr. Christophe Craeye for supporting the development of the fast electromagnetic solver, Dr. Maciej Syralek and Dr. Robert Laing for their support. We also thank the anonymous reviewers for their comments and suggestions, which have improved the quality of this paper.

This research was supported by ESA, NPL, and the Science and Technology Facilities Council (STFC). OOH was supported by grant number G109464, ‘The Design of Highly Sensitive EM Sensors for Space Applications’, QG and FD by grant number ST/W00206X/1, EdLA by grant number ST/V004425/1 and DA and JC were supported by grant number ST/X00239X/1. AF is supported by a Royal Society University Research Fellowship #180523. JD acknowledges support from the Boustany Foundation and Cambridge Commonwealth Trust in the form of an Isaac Newton Studentship.

DATA AVAILABILITY

The data and software underlying this article will be shared on reasonable request to the corresponding authors.

Footnotes

REFERENCES

APPENDIX: THE DELAY POWER SPECTRUM

Here, we provide a brief overview of the link between the power spectrum associated with a temperature field, the statistics we aim to measure, and the visibilities, which are the observable quantities. This relationship is straightforward, enabling the instrument to function almost as a direct probe of the signal power spectrum. We begin with an estimation of the two-point correlation function, |$\xi (\boldsymbol {r})$| of a brightness temperature |$\mathcal {T}(\boldsymbol {r})$| as follows:

where |$\mathcal {V}$| is a comoving volume of interest (Hogg 2000). The power spectrum |$\mathcal {P}(\boldsymbol {k})$| associated with this temperature field, expressed in [mK|$^2\, h^{-3}$| Mpc|$^3$|], corresponds to the Fourier transform of the correlation function |$\xi (\boldsymbol {r})$|⁠, which transforms comoving distances |$\boldsymbol {r}$| to wavevectors |$\boldsymbol {k}$| (Parsons et al. 2012a, b):

According to Parseval’s theorem, we have |$\tilde {\mathcal {P}}(\boldsymbol {k}) =|\tilde{\mathcal {T}}(\boldsymbol {k})|^2$| where |$\tilde{\mathcal {T}}(\boldsymbol {k})$| denotes the Fourier transform of the temperature field |$\mathcal {T}(\boldsymbol {r})$| (Liu & Shaw 2020). The observed temperature |$\mathcal {T}$| is associated with a specific point in co-moving coordinates and corresponds to an epoch defined by a particular redshift z. We now express |$\mathcal {T}$| as a function of local observational coordinates. Assuming the observed volume |$\mathcal {V}$| is located at zenith and occupies a sufficiently small portion of the sky relative to the beamwidth, we can align the z-axis along the line of sight and apply the following transformation:

This change of variable maps co-moving coordinates |$(r_x,r_y,r_z)$| to observation coordinates |$(l,m,f)$|⁠. Here, |$X,Y$| depends on the redshift z and are defined as Lanman et al. (2020):

where |$\chi (z)$| is the co-moving distance, |$H(z)$| is the Hubble parameter and |$f_{21}$| is the rest-frame 21 cm frequency. The relationship in equation (A3) also establishes a mapping between the k-space and the observational baseline-delay domain: |$(u,v,\tau) = (k_x X(z)/\lambda _0,k_y X(z)/\lambda _0, k_zY(z))/(2\pi)$|⁠. We decompose the wavevector |$\boldsymbol {k}$| into components parallel and perpendicular to the line of sight |$\hat{\boldsymbol {z}}$|⁠:

Using the notation |$T(\hat{\boldsymbol {s}},f) = \mathcal {T}(\boldsymbol {r})$|⁠, we can express the Fourier transform with respect to the coordinates |$(l,m,f)$| as follows:

Here, we thus have |$\tilde{\mathcal {T}}(\boldsymbol {k}) = \tilde{T}(\tau ,\boldsymbol {u})/(X^2Y)$| for purely planar baseline (⁠|$\boldsymbol {b} = \boldsymbol {u}$|⁠). The relationship in equation (A6) is fundamental because |$\tilde{T}(\tau ,\boldsymbol {b})$| can be nearly directly estimated from the time-delay transform of the instrument observable: the visibilities |$V(f,\boldsymbol {b})$|⁠. The main distinction is that the visibilities account for the instrument’s response through the BTF, |$A(\hat{\boldsymbol {s}},f)$|⁠. When the width of the BIR, |$\tilde {A}(\hat{\boldsymbol {s}},\tau)$| (in both angular |$\hat{\boldsymbol {s}}$| and time dimensions |$\tau$|⁠), is significantly narrower than the temperature field |$\tilde {T}(\hat{\boldsymbol {s}},\tau)$|⁠, we can approximate the power spectrum |$\mathcal {P}$| as follows (Parsons et al. 2014):

where |$\Omega _{A}$| is a normalization factor (a full derivation can be found in appendix B of Parsons et al. 2014) that accounts for the instrument response and is evaluated here as:

The approximation (A7) thus serves as an estimator of the power spectrum of a given signal of brightness temperature T.