MNRAS 538, 31–48 (2025) https://doi.org/10.1093/mnras/staf264 Advance Access publication 2025 February 12 Unco v ering the effects of array mutual coupling in 21-cm experiments with the SKA-Low radio telescope Oscar S. D. O’Hara , 1 , 2 ‹ Quentin Gueuning, 1 , 2 ‹ Eloy de Lera Acedo, 1 , 2 Fred Dulwich, 1 John Cumner , 1 , 2 Dominic Anste y , 1 , 2 Anthon y Brown, 3 Anastasia Fialkov, 2 , 4 Jiten Dhandha , 2 , 4 Andrew Faulkner 1 and Yuchen Liu 1 , 2 1 Cavendish Astrophysics, University of Cambridge, Cambridge CB3 0HE, UK 2 Kavli Institute for Cosmology in Cambridge, University of Cambridge, Cambridge CB3 0HA, UK 3 Queen Mary, University of London, London E1 4NS, UK 4 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Accepted 2025 February 7. Received 2025 January 28; in original form 2024 December 2 A B S T R A C T 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 ‖ ∼ 2 h −1 Mpc for all station layouts and across all k ⊥ modes, completely obscuring the detection window. We find that attempting to remo v e 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 ∼ 10 11 mK 2 h −3 Mpc 3 . This is still seven orders of magnitude brighter than the expected level of the EoR signal ( ∼ 10 4 mK 2 h −3 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 ef fecti v e fore ground remo val. Our research pro vides 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. Key words: scattering – instrumentation: interferometers – dark ages, reionization, first stars. 1 O c a o s U C m fl t e i s s � b t e f e 2 t B a d d r a r © P C p D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 I N T RO D U C T I O N ver the past decade, the detection of the faint cosmological 21- m signal, arising from the hyperfine transition of neutral hydrogen toms, has received increasing attention as one of the frontiers of bservational radio astronomy. The signal promises a rich under- tanding of astrophysics, ranging from the earliest epochs of the niverse (the Dark Ages), through the first luminous sources (the osmic Dawn), and the eventual ionization of the intergalactic edium (epoch of reionization, or EoR). Measuring the spatial uctuations of this signal would unlock a deeper understanding of he interplay between different heating and ionizing sources in the arly Universe. The primary obstacle preventing such a detection s the isolation of the 21-cm signal from astrophysical foregrounds, uch as Galactic synchrotron emission, radio galaxies, pulsars, and upernova remnants, which are up to five orders of magnitude E-mails: osdo2@cam.ac.uk (OSDO); qdg20@cam.ac.uk (QG) f t i 2025 The Author(s). ublished by Oxford University Press on behalf of Royal Astronomical Society. Th ommons Attribution License ( https:// creativecommons.org/ licenses/ by/ 4.0/ ), whic rovided the original work is properly cited. righter than the signal itself (Jeli ́c et al. 2008 ). F ore ground a v oidance echniques, which exploit the rapid spectral variations of faint signals, nable their isolation at higher delays, helping to separate them rom foreground emissions (Trott, Wayth & Tingay 2012 ; Parsons t al. 2012b ; Liu, Parsons & Trott 2014a , b ; Th yag arajan et al. 015 ; Chapman et al. 2016 ). Alternativ ely, fore ground remo val echniques (Bowman, Morales & Hewitt 2006 ; Liu et al. 2009 ; onaldi & Brown 2015 ; Chapman et al. 2016 ; Pober et al. 2016 ) im to subtract modelled foreground contributions from the observed ata using knowledge of the telescope’s response and the celestial istribution of the foregrounds. In both foreground a v oidance and emoval approaches, precise models of the instrument’s beams and nalogue chains are essential. Typically, two fundamental questions elated to these models must be addressed: (i) How accurately must we model the instrument’s response? This ocuses on determining the necessary level of accuracy in matching he real instrument response with the model to ensure reliable data nterpretation. is is an Open Access article distributed under the terms of the Creative h permits unrestricted reuse, distribution, and reproduction in any medium, http://orcid.org/0009-0006-3633-5816 http://orcid.org/0000-0001-8479-3746 http://orcid.org/0000-0003-1742-7417 http://orcid.org/0000-0002-1481-0907 mailto:osdo2@cam.ac.uk mailto:qdg20@cam.ac.uk https://creativecommons.org/licenses/by/4.0/ 32 O. S. D. O’Hara et al. M e o q 2 t i e h m ( M t M ( i p p E S t e a s a E a M s s a A l a d T l L e ( P a L f S 2 e ( c e f fi ( o m a p I e ( s f p c t 2 a s A ( a e P e M t i r ( n o t e l e a ( w t i G e ( t h e s 2 1 s fi u l a e h f b a K t G p t ( s s D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 (ii) In which regions of the observational domain are instrumental ffects weakest? Identifying these regions will enable us to focus n areas of minimal instrumental uncertainty to preserve the data uality. In addressing the accuracy criterion [question (i)], forecasts of the 1-cm power spectrum from CHIME (Bandura et al. 2014 ) indicate hat achieving a 0.1 per cent accuracy in the primary beamwidth s necessary for reliable measurements (Shaw et al. 2015 ; Amiri t al. 2022 ). The REACH experiment (de Lera Acedo et al. 2022 ) as similarly determined that uncertainties in the antenna directivity ust remain below −40 dB to ensure robust global 21-cm detection Cumner et al. 2024b ). The MeerKAT Radio Telescope (Jonas & eerKAT Team 2016 ) has e v aluated the required pointing accuracy o use beam patterns with 1 per cent accuracy in power (Jonas & eerKAT Team 2016 ; de Villiers & Cotton 2022 ). For SKA-Low Dewdney et al. 2009 ), preliminary studies (Nasirudin et al. 2022 ) ndicate that inaccuracies in beam knowledge – resulting from small ositional deviations and up to 5 per cent of broken antennas – can roduce residual errors as high as 1 per cent in the power spectrum oR window and around 10 per cent in the foreground wedge. imilarly, for MWA (Tingay et al. 2013 ), 15–40 per cent of tiles ypically have at least one broken dipole, leading to beam modelling rrors that produce foreground residuals two orders of magnitude bo v e the 21-cm signal after calibration (Joseph et al. 2019 ). With uch error levels, systematic artefacts could well be misinterpreted s 21-cm structures, (i.e. Barry et al. 2016 ; Th yag arajan et al. 2016 ; wall-Wice et al. 2017 ). For radio telescopes composed of densely packed antennas, ddressing questions (i) and (ii) is complicated by a major issue: utual Coupling (MC) (Craeye & Gonz ́alez-Ovejero 2011 ). In situ , trong electromagnetic interactions between antennas can induce harp angular and spectral variations in the radiated far field, known s embedded element patterns (EEPs; Bird 2021 ; Bolli et al. 2022c ; nstey et al. 2024 ). Since each antenna operates in a distinct ocal environment within the array, the EEPs can vary significantly cross elements in the station. In other words, MC effects are irection-dependent, frequency-dependent, and baseline-dependent. o date, MC is identified as a systematic source that potentially imits interferometric 21-cm power spectrum experiments such as the o w-Frequency Array (LOFAR) (v an Haarlem et al. 2013 ; Mertens t al. 2020 ), the Hydrogen Epoch of Reionization Array (HERA) DeBoer et al. 2017 ; Kern et al. 2019 , 2020 ), the Precision Array for robing the Epoch of Reionization (PAPER) (Parsons et al. 2010 ), nd the Murchison Widefield Array (MWA) (Beardsley et al. 2016 ). OF AR, NenuF AR (Zarka et al. 2012 ), and MWA have implemented ull-Jones directional-dependent sky-based calibration algorithms AGECAL (Kazemi et al. 2011 ; Yatawatta, Kazemi & Zaroubi 012 ), DDECAL (van Diepen, Dijkema & Offringa 2018 ; Gan t al. 2023 ), and RTS (Mitchell et al. 2008 ). Incorporating realistic simulated or measured) beam models into these calibration pipelines an help to reduce residual errors (Subrahmanyan 2021a ; Chokshi t al. 2024 ). An alternative approach is mitigating MC effects using ringe rate filters. In this context, the HERA Collaboration modelled rst-order dish–dish interactions with a semi-analytical formalism Josaitis et al. 2022 ; Rath et al. 2024 ) to demonstrate the impact f MC at non-zero fringe rates. This informed the development of itigation strategies and tailored fringe rate filters for foreground v oidance, b uilding on Parsons et al. ( 2016 ). The Square Kilometre Array Lo w (SKA-Lo w), examined in this aper, is an interferometric array under construction in Inyarrimanha lgari Bundara, Western Australia. It will consist of 512 stations, NRAS 538, 31–48 (2025) ach containing 256 dual-polarization SKA Log-periodic antennas SKALA) antennas within a 19 m radius, spanning ∼ 74 km. Initial tudies hav e inv estigated the spectral smoothness of the bandpass or individual SKALA (de Lera Acedo et al. 2017 ) and estimated reliminary calibration requirements (Trott et al. 2017 ). These studies oncluded that the isolated antenna response is sufficiently smooth o reco v er the EoR signal through statistical methods (Trott & Wayth 016 ). Ho we ver, as noted in de Lera Acedo et al. ( 2017 ), these nalyses did not include the effects of MC between antennas within a tation. The impact of MC has been a long-standing concern (de Lera cedo 2012 ), with variations in EEPs observed both experimentally Raucy et al. 2013 ; de Lera Acedo, Pienaar & Fagnoni 2018 ) nd through simulations (Gueuning et al. 2015 ; de Lera Acedo t al. 2016 ). Recent UAV measurements (Subrahmanyan 2021b ; aonessa et al. 2023 ) and electromagnetic simulations (Davidson t al. 2020 ; Bolli et al. 2022b ; Anstey et al. 2024 ) confirm that C introduces notches at frequencies (55, 78, and 125 MHz) within he EoR band (50–200 MHz). A significant milestone has been the ncorporation of full-wave EEP simulations into calibration algo- ithms, demonstrating a substantial reduction in visibility residuals Subrahmanyan 2021a ). These results underscore the importance of umerical simulations for understanding and mitigating the effects f MC on the delay power spectrum of the SKA-Low telescope. In his context, our analysis extends the earlier studies (de Lera Acedo t al. 2017 ; Trott et al. 2017 ) by explicitly addressing the MC effects. Accurately modelling SKA-Low beams requires full-wave simu- ations of the electromagnetic response of each of the 256 antennas mbedded in large (19 metre radius), dense, and irregular arrays cross a wide bandwidth (50–350 MHz) with high spectral resolution below 1 MHz). The complex antenna geometry, small spacing (0.01 avelengths at the lowest frequency), and large station size (up o 45 wavelengths at the highest frequency) result in a challeng- ng multiscale Computational Electromagnetics (CEM) problem. eneral-purpose solvers can take days or weeks per frequency point, ven when accelerated with the Multi-Level Fast Multipole Method MLFMM; Engheta et al. 1992 ) on large workstations. To o v ercome his, state-of-the-art CEM methods tailored for large, irregular arrays av e been dev eloped (Maaskant, Mittra & Tijhuis 2008 ; Bui-Van t al. 2018 ; Chose et al. 2023 ; Conradie et al. 2023 ). The in-house olv er, named F ast Array Simulation Tool (FAST) (Gueuning et al. 025 ), simulates all the antenna responses in the station in about 0 min per frequency point on a current standard laptop, with re- imulation for a different station layout taking just 1 min, enabling ne spectral sampling and therefore e v aluation of impulse responses p to late time delays. Modelling SKA-Low visibilities also presents computational chal- enges. Each of the 512 stations generates beamformed signals, which re then cross-correlated to produce o v er 130 000 visibility pairs. To nsure accuracy on SKA-Low’s longest baselines (up to ∼ 74 km), igh-resolution sky maps with tens of millions of pixels are essential or a v oiding aliasing in the sk y inte grals. Furthermore, each station’s eam, derived from precomputed EEP data, must be re-e v aluated cross all these sky directions. To address this, the tool Oxford Square ilometre Array (OSKAR) has been developed Dulwich ( 2020 ) o leverage GPU parallelization, and we have also implemented a PU-parallelized approach for station beam e v aluation, originally roposed in Bui-Van et al. ( 2018 ). In this paper, we use the in-house tools, FAST and OSKAR, o investigate the accuracy required for the EEP models [question i)] as well as the impact of MC effects on the time-delay power pectrum of foreground emissions [question (ii)] through numerical imulations. To address question (i), we estimate delay power Mutual coupling effects in 21-cm experiments 33 s s o 2 i d a c T i S d S m b f p t s a a a S 2 B c r c a C c r T c g k w t b w o r p c w g ŝ w l Figure 1. Block diagram of a radio telescope composed of phased array sta- tions. Channelized voltages v ant ,np ( f , t) from each antenna are beamformed within each station p, and the resulting beamformed voltages are then cross- correlated. 3 V A t i f d d t b 3 T c r s v a w v V c a 3 T c m i r d D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 pectrum residuals after fore ground remo val, using an approximate tation beam model that is corrupted by random noise, or based n the average embedded element pattern (AEP) (Wijnholds et al. 019 ) or spectrally interpolated EEPs. Regarding question (ii), we nvestigate MC effects by comparing beam and visibility time- elay responses across three station layouts – regular , sunflower , nd random. We quantify foreground leakage into longer delays aused by MC, considering both rotated and non-rotated stations. his paper also highlights the computational costs associated with ntegrating array patterns, which fully account for MC effects, in the KA-Low pipeline. This paper is organized as follows. Section 2 efines the conventions for the coordinates systems used, while ection 3 presents the mathematical formulation for the forward odel of these visibilities, based on EEPs, array patterns, and sky rightness distribution. In Section 4 , we outline the parameters used or the telescope simulations. Section 5 describes our simulation ipeline, emphasizing the recent advancements introduced in the ools, FAST and OSKAR. Section 6 presents the simulation re- ults, including an analysis of the station beam impulse response cross different layouts, the foreground spill over induced by MC, nd the residual levels after foreground removal using several pproximate beam models. Finally, we summarize our findings in ection 7 . C O N V E N T I O N S efore outlining our mathematical deri v ations, we would like to larify the notation used by the coordinate systems in this paper. Comoving distances are spatial coordinates that measure sepa- ation in cosmology. Per their definition, these distances remain onstant o v er time despite the e xpansion of the univ erse, as the y ccount for the changing scale at different epochs. Assuming a artesian coordinate system with unit vectors ( ̂ x , ̂ y , ̂ z ), the comoving oordinates will be denoted as follows: = r x ̂ x + r y ̂ y + r z ̂ z . (1) he spatial frequency coordinates, which are associated with the osmological Fourier domain and referred to as the k-space, are iven by = k x ̂ x + k y ̂ y + k z ̂ z , (2) here k is the wav ev ector with cosmological wavenumber k = || k || . Physical relati ve distances, kno wn as baselines, are represented by he baseline vector b : = u ̂ x + v ̂ y + w ̂ z , (3) hen expressed in metres, or by b λ = b /λ0 when expressed in units f the free-space wavelength λ0 . The vector u λ = u/λ0 ̂ x + v/λ0 ̂ y epresents the in-plane projection of the baseline vector b λ. In this aper, the baseline notation b refers to the separation between station entres, while b ant (in metres) denotes the position of an antenna ith respect to its station centre. The directional unit vector ˆ s , which represents a line of sight, is iven by = l ̂ x + m ̂ y + n ̂ z , (4) here n = (1 − l 2 − m 2 ) 1 / 2 , ensuring that ˆ s is constrained to unit- ength. A F O RWA R D M O D E L O F TELESCOPE ISIBILITIES forward model simulates how signals from a celestial distribu- ion are transformed by the instrument, capturing the distortions ntroduced at each stage of the system. This model is essential or understanding how instrumental effects propagate through the igital processing pipeline and manifest in the final product, the elay power spectrum. In this section, we focus on modelling he measured visibilities, specifically including the impact of MC etween antennas in the beam response of each station. .1 System o v er view he interferometer set-up comprises N s phased array stations, each ontaining N a antennas. After digitization and channelization of aw voltages from each antenna, the voltage from antenna n in tation p is represented as a time-varying mono-frequency signal ant ,np ( t, f ). The N a voltages in each station are then beamformed s v p ( f , t) = ∑ n w np ( f ) v ant ,np ( f , t), with beamforming weights np ( f ) chosen to steer the beam or control its shape. Beamformed oltages are then cross-correlated to yield the complex visibility pq ( f ) = V ( b λ,pq , f ), where b λ,pq is the baseline vector linking the entres of stations p and q. A block diagram for a two-station, three- ntenna configuration is shown in Fig. 1 . .2 Mathematical model of the visibilities he following derivations outline how the visibilities, measured by ross-correlating beamformed signals from different stations, are odelled using both celestial distribution and beam models. Let us first assume a monochromatic incident plane wave propagat- ng in the direction ̂ s , with a time-varying electric field vector E i ( f , t) epresenting its amplitude and polarization, where the subscript i enotes incident. Assuming and suppressing the time dependence MNRAS 538, 31–48 (2025) 34 O. S. D. O’Hara et al. M e E w p r p t a t g v w i i c Z i i v v w A w n i i s d o t ( f [ w v t a t [ w A f e a v t V w A w A F ( g � T u n m 3 A i F V a o W V S c V w o W B t d a 3 T ( w s b a V T u A o s 2 τ D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 j2 πf t , the electrical field of the incident plane wave is expressed as ( f , t, b ) = E i ( f , t) e −jk 0 ̂ s ·b , (5) here k 0 is the free-space wavenumber. First, we consider single- olarized antennas (with a single feed port). According to Lorentz’s eciprocity theorem, the voltage v ant ,np ( f , t) appearing at the output ort of a passive antenna n belonging to a station p illuminated by he plane wave can be obtained by treating antenna n as active while ll the other antennas are passively terminated (Kildal 2015 ). Using he centre of antenna n as a phase reference in ( 5 ), this voltage is iven by (Craeye & Gonz ́alez-Ovejero 2011 ): ant ,np ( f , t) = 2 λ0 jη0 Z L,np ( f ) EEP np ( f , ̂ s ) · E i ( f , t) , (6) here η0 is the free-space characteristic impedance, Z L,np is the input mpedance of the amplifier connected to antenna n . The EEP, EEP np n equation ( 6 ), is expressed in volts and e v aluated in the active ase by feeding antenna n with a The venin equi v alent of impedance L,np and a source of unit voltage. In the passive case described n equation ( 6 ), the EEP is implicitly normalized by 1 volt, so that t becomes a dimensionless quantity. From there, the beamformed oltage associated with station p can be expressed as p ( f , t) = 2 λ0 jη0 Z 0 AP p ( f , ̂ s ) · E i ( f , t) , (7) here the array pattern (AP) is defined by P p ( f , ̂ s ) = N a ∑ n = 1 c pn ( f ) EEP np ( f , ̂ s ) e jk 0 ̂ s ·b ant ,np , (8) ith the coefficients c pn = w pn Z L,np /Z 0 . The phase reference is ow taken at the centre of the station, introducing a phase factor n equation ( 8 ) based on the relative antenna distances b ant ,np . The mpedance Z 0 in equation ( 7 ) has been introduced so that the AP hares the same units as the EEPs. Considering now dual-polarized antennas with two feed ports, enoted as v p,X and v p,Y , with indices X and Y referring to the rientation of the associated dipoles along ̂ x - or ̂ y -axis, we re-express he beamformed voltages using the 2 × 2 Jones matrix formalism Jones 1941 ; Smirnov 2011 ). Omitting the explicit time ( t) and requency ( f ) dependence, this yields v p,X v p,Y ] = [ J p ,X X J p ,X Y J p,YX J p,YY ][ E i,X E i,Y ] (9) here E i = E i,X ̂ e X + E i,Y ˆ e Y , with ̂ e X and ̂ e Y representing the unit ectors corresponding to the horizontal and vertical components in he Ludwig-III polarization coordinate system (Ludwig 1973 ). In greement with our pre vious deri v ation (equation 7 ), the elements of he Jones matrix are given by J p ,X X J p ,X Y J p,YX J p,YY ] = 2 λ0 jη0 Z 0 [ AP p ,X X AP p ,X Y AP p,YX AP p,YY ] (10) here AP p,X = AP p ,X X ̂ e X + AP p ,X Y ˆ e Y and AP p,Y = AP p,YX ̂ e X + P p,YY ˆ e Y are the APs for feed ports X and Y , respectively. This ormalism allows us to apply the radio interferometric measurement quation (RIME), as defined in Hamaker, Bregman & Sault ( 1996 ) nd Smirnov ( 2011 ), to relate the Stokes I parameter of the measured isibilities V pq to the spectral brightness distribution B( f , ̂ s ) and to he station beams: pq ( f , b λ) = “ A pq ( f , ̂ s ) B( f , ̂ s ) e −j2 π ˆ s ·b λ,pq d l n d m, (11) NRAS 538, 31–48 (2025) ith the instrument response: pq ( f , ̂ s ) = 1 �pq ( f ) ( J p ,X X J ∗ q,XX + J p ,X Y J ∗ q,XY + J p,YX J ∗ q,YX + J p,YY J ∗ q,YY ) , (12) here ∗ denotes the complex conjugate. We refer to this function pq ( f , ̂ s ) as the beam transfer function (BTF) while its time-delay ourier counterpart ˜ A pq ( τ, ̂ s ) is called the beam impulse response BIR). The factor �pq is called the integrated beam response and is iven by pq ( f ) = “ | A pq ( f , ̂ s ) | d l n d m (13) his factor ensures that the visibilities V pq are expressed in Jy nits when the brightness B is expressed in Jy sr −1 . To alleviate the otations, we will temporarily assume beams similar across stations, aking the BTF ef fecti vely baseline independent, that is A pq = A . .3 Time-delay transform s explained in Appendix A , the power spectrum ˜ P of an astronom- cal signal can be estimated from the amplitude of the time-delay ourier transform ˜ V ( τ, b ), expressed in JyHz, of the visibilities ( f , b λ). In practice, the visibilities are only measured within limited sub-band ( f min , f max ), a subset of the telescope’s full perational range. To minimize truncation ef fects, a windo w function ( f ) is applied to the visibilities, resulting in ˜ ( τ, b ) = ∫ f max f min W ( f ) V ( f , b λ) e −2 πjf τ d f . (14) ubstituting equation ( 11 ) into equation ( 14 ) and applying the onvolution theorem, this expression becomes: ˜ ( τ, b ) = ˜ W ( τ ) ∗ “ ˜ A ( τ, ̂ s ) ∗ ˜ B ( τ − ˆ s · b /c 0 , ̂ s ) d l n d m, (15) here ∗ denotes a convolution product, c 0 is the free-space speed f light, ˜ W ( τ ) is the time-delay response of the window function , and the spectrum ˜ B is the Fourier transform of the brigthness . The spectral sampling step �f is a critical factor in computing he delay Fourier transform ( 14 ). For time-limited functions over uration T , the Nyquist criterion requires �f = π/ (2 T ) to confine liasing errors outside the observed delay range τ ∈ [0 , T ]. .4 Instrumental effects in the baseline-delay domain ime-delay visibilities are analysed in baseline length-delay space, b, τ ), by grouping visibilities into baseline length bins and averaging ithin each bin to obtain the average visibility, ˜ V av ( b, τ ). For implicity, we consider the case of a single point source denoted y index i of flux density S( f ) = B( f , ̂ s ) �l i �m i /n i in direction ̂ s i nd solid angle �l i �m i /n i . The relation ( 15 ) thus reduces to ˜ ( τ, b ) = ˜ W ( τ ) ∗ ˜ A ( τ, ̂ s i ) ∗ ˜ S ( τ − ˆ s i · b /c 0 ) , (16) he effect of the instrument appears as a cascade of convolution prod- cts in the time-delay domain, which we will analyse sequentially. ll these effects are summarized and sketched in Fig. 2 . First, the baseline exponential in equation ( 11 ) causes a translation f sources flux density spectrum, denoted as ˜ S , by a time delay that cales with the projection (Parsons et al. 2012b ; Chapman et al. 016 ), = ̂ s · b /c 0 . (17) Mutual coupling effects in 21-cm experiments 35 Figure 2. A log space illustration of instrumental effects in the delay-baseline domain. T p i o d a d t l e d B b l B a m s h w s p t u o g t f f n 2 P m H d 3 A V t f P w t r f t b k F f c 4 T b c c 4 g t t 4 T d o f f a 2 i i a e f s R l B t G a w t t i D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 his relation indicates that, for planar baselines with b = u , sources ositioned at the zenith introduce no delay, while sources oriented n line with or opposite to the baseline direction yield either positive r ne gativ e delays. The av erage visibility V̄ av ( b, τ ) peaks around the elay τ = b sin θi /c 0 with the source zenith angle θi , thus following line with slope sin θi /c 0 in baseline length-delay space ( b, τ ). This elay reaches its maximum, the horizon limit , at τ = b/c 0 , when he source is on the horizon (Morales et al. 2012 ). As the baseline ength b increases, the separation between spectrum peaks at different le v ation angles θi increases, enhancing resolution in the time-delay omain. For sources with spatially smooth brightness distribution (diffuse emission), the visibility V̄ av ( b, τ ) decreases rapidly with aseline length b and increases with delay τ peaking at the horizon imit. Secondly, the shifted source spectrum ˜ S is convolved with the IR, ˜ A ( τ, ̂ s i ), in the direction ̂ s i . Sources within the main lobe thus ppear brighter than those in the sidelobes, and higher sidelobe levels anifest in the ( b, τ ) plots as increased peak values around τ = in θi b/c 0 . The beam limit, defined by τ = b sin θi /c 0 , is based on the alf-beamwidth θi at the centre frequency. The BIR can decay slowly ith delay τ , and when convolved with a more rapidly decaying ignal spectrum ˜ S , such as foreground emission, it spreads the signal ower into longer delays, causing leaka g e that can extend far beyond he horizon limit, as shown in Section 6 . This slow decay arises from nwanted electromagnetic interactions, such as multiple reflections r scattering, which prolong wave attenuation. Optimizing antenna eometries and array configurations during design can help minimize he BIR decay rate. Finally, ˜ W ( τ ) in equation ( 16 ) is the Fourier transform of the requency window function W ( f ), used to reduce ringing artefacts rom truncating the frequency band to ( f min , f max ). A window with a arrow first lobe and high dynamic range is preferred (Lanman et al. 020 ); here, we use the convolved Blackman–Harris window as in arsons et al. ( 2014 ). The first lobe of the combined response ˜ W ∗ ˜ A atches or exceeds that of ˜ W , depending on beam chromaticity. o we ver, the dynamic range is generally governed by the slower ecay of the beam response ˜ A . .5 Delay power spectrum fter applying the time-delay Fourier transform to the visibilities pq , yielding the time-delay visibilities ˜ V pq , we can then e v aluate he delay power spectrum from these time-delay visibilities using the ollowing formula (Morales & Hewitt 2004 ): ˜ ( k ⊥ , k ‖ ) � X 2 ( z ) Y ( z ) �A ‖ ̃ V ( τ, b ) | 2 , (18) here z is the redshift, X( z) and Y ( z) are cosmological scaling fac- ors, and �A is a normalization factor that depends on the instrument esponse A . The delay power spectrum is often represented as a unction of the parallel and perpendicular components, k ⊥ and k ‖ , of he cosmological wav ev ector k defined in equation ( 2 ). The link to aseline time-delay coordinates is given by ⊥ = 2 πu λ X( z) , k ‖ = 2 πτ Y ( z) . (19) or interested readers, a detailed reminder about the deri v ation of ormulas ( 18 ) and ( 19 ), along with an expansion of the various onstants, is provided in Appendix A . N U M E R I C A L EXPERI MENT: SKA-LOW he spiral configuration of SKA-Low features expanding arms with aseline lengths up to ∼ 74 km for higher resolution, alongside a ore of 224 densely packed stations that provide instantaneous uv- o v erage with minimum baseline lengths starting from around b min = 0 m to around b max = 1000 m. One of SKA-Low’s key scientific oals is the detection of the 21-cm signal. In this section, we outline he parameters of a realistic scenario considered for the detection of his signal, with results to be presented in Section 6 . .1 Obser v ation parameters he SKA Epoch of Reionization Science Working Group has eclared their intent to target the 21-cm signal between the redshift f 6 < z < 30 (Braun et al. 2015 ) and therefore an approximate requency range of 46 MHz < f < 202 MHz. In this work, we ocus on a sub-band from f min = 120 MHz to f max = 150 MHz at 100 kHz spectral resolution ( �f ) placing us firmly within the 1-cm band while allowing us to highlight the impact of MC- nduced features (Davidson et al. 2020 ; Anstey et al. 2024 ) present n the beam. To minimize the ratio between the 21-cm signal and strophysical foregrounds, observation times are typically chosen to nsure that the beam’s primary lobe and the field of view are relatively ree from bright foreground emissions. This inv olves a v oiding areas uch as the Galactic centre (GC) and prominent point sources like igel. An example of a suitable observation field is the MWA EOR1, ocated at RA = 60 ◦ and Dec . = −27 ◦ (see Bowman et al. 2013 ; eardsley et al. 2016 ). For our simulations, we select two observation imes. The first is at 12:31:30 UTC on 2000 January 1, when the alactic plane is far from the SKA-Low field of view and located long the horizon. The second is at 09:31:30 UTC on January 1, 2000, hen the GC is below the horizon, minimizing any leakage through he far-out sidelobes. The sky map at 135 MHz for both observation imes is shown in Fig. 3 and a description of the sky model content s outlined in Section 4.3 . MNRAS 538, 31–48 (2025) 36 O. S. D. O’Hara et al. M Figure 3. The composite sky model illustrating the source flux density as Dirac delta functions at 135 MHz, observed at two times: (top) 2000-01-01 09:31:30 [UTC] with the GC along the horizon and (bottom) 2000-01-01 12:31:30 [UTC] with the GC below the horizon. 4 4 T d ( f d i 4 F e Figure 4. Geometry of the SKALA4 (de Lera Acedo et al. 2018 ). a a s e e t 4 T t c b L N a s c e p a t 2 4 T a t T d t D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 .2 Telescope model .2.1 Antenna geometry he 4th version of the SKA element, illustrated in Fig. 4 , is a ual-polarized log-periodic antenna featuring 16 dipoles per arm SKALA4, (de Lera Acedo et al. 2018 ). It operates o v er a 7 : 1 requency range, spanning from 50 to 350 MHz. Each arm is ifferentially amplified using a low-noise amplifier (LNA) with an nput impedance of 100 Ohms. .2.2 Station layout ig. 5 illustrates the three different layouts used in this numerical xperiment. Each layout is arranged in a circular footprint with NRAS 538, 31–48 (2025) radius of 19 m and a minimum distance of 1.7 m between the ntennas. (i) regular: A regularly spaced grid of elements within the circular tation footprint, each element positioned 2.14 m apart. (ii) sunflower: A sunflower-head array, with the position of each lement given in Vogel ( 1979 ) and e v aluated for SKA-Low in Bolli t al. ( 2022a ). (iii) random: A random layout which results from random per- urbation of the sunflower layout (Anstey et al. 2024 ). .2.3 Array layout he coordinates of the array centre and the projected spacing between hese stations were obtained from Revision 3 of the SKA-Low onfiguration and illustrated in Fig. 6 . We exclude all the long aselines and limit our simulations to the 224 core stations of SKA- ow, to ensure that the currently available sky model meets the yquist resolution. Note also that during our analysis, the baselines re irregularly binned based on a logarithmic length distribution pecific to SKA-Low. In the following, we consider the case where all stations in the ore have an identical, un-rotated station layout. Additionally, we xamine a scenario where each station (antenna orientation and osition relative to the station centre) is rotated by a unique angle, s illustrated in Fig. 6 . The rotation angles are assigned randomly o each station, following a uniform distribution between 0 and π . .3 For egr ound model o accurately model the localized and extended foreground emission ccessible to the SKA-Low, a composite sky model is created by reating the extended emission as a set of unresolved point sources. he sky model components include: (i) A diffuse Galactic foreground emission from the desourced and estriped Haslam 408 MHz surv e y, with ∼ 5 × 10 7 pixels, equating o a 1.718 arcmin resolution (Remazeilles et al. 2015 ). Mutual coupling effects in 21-cm experiments 37 Figure 5. The antenna positions within a SKA-Low station are shown for the ‘regular,’ ‘sunflower,’ and ‘random’ layouts. Figure 6. The centres of the 224 stations in the SKA-Low core, each with a unique rotation applied. Baseline lengths range from approximately 40 m to 1 km. A e m T n B 1 A p Figure 7. Flow-chart of the simulation pipeline. I t w s 5 W s v c e t D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 (ii) A collection of point sources 1 from GaLactic and Extragalactic ll-Sk y MWA Surv e y (GLEAM), with ∼ 3 × 10 5 sources (Wayth t al. 2015 ). At low frequencies, radio foreground observations may be approxi- ated by a near-featureless power law across the frequency spectrum. he composite map can be frequency-scaled within OSKAR as eeded using the equation: ( f , ̂ s ) = B( f 0 , ̂ s ) [ f f 0 ]β . (20) Note that the GLEAM catalogue contains A-team sources, such as Fornax , which exhibit extended spatial scales that are challenging to model using oint sources. For more accurate modelling, refer to Line et al. ( 2020 ). p m p s n this equation, f 0 denotes the reference frequency, and β signifies he spectral inde x. F or the Haslam 408 MHz surv e y, β = −0 . 7, hereas for the GLEAM surv e y, a distinct β is provided for each ource. AC C E L E R A TED SIMULA T I O N TO O L S e present the simulation pipeline used to generate the delay power pectrum for large interferometric arrays. Accurately modelling isibilities for telescopes with many antennas and stations poses omputational challenges, as it requires full-wave simulations of ach antenna’s electromagnetic response, including MC effects, and he e v aluation of APs o v er multiple directions and frequencies. This ipeline leverages two in-house tools: FAST for efficient electro- agnetic simulations and the radio-telescope simulator OSKAR for arallelized interferometric calculations. The pipeline, outlined in Fig. 7 , begins with electromagnetic imulations for the 2 N a EEPs in a frequency band ( f min , f max ) at MNRAS 538, 31–48 (2025) 38 O. S. D. O’Hara et al. M Table 1. Computation times to e v aluate 512 EEPs with FAST per frequency point. Solver: Fast Station layout Random Machine Intel Core i7-13800H, 2.50 GHz, 14 cores, 32.0 GB RAM Pre-computations 10.1 min Matrix filling 0.23 min Solve time 0.51 min Total time 10.8 min Update time per new layout 0.7 min Peak memory 1.3 GB r F f c g t W t t p 5 T a h c i a l a a M d f a a V e t w s o r 5 O s r t i 2 3 4 Table 2. Computation times to simulate visibilities with OSKAR. OSKAR: common properties Machine 2 × AMD EPYC 7763 64-Core Processor 1.8 GHz, 1 TB RAM, 4 × NVIDIA A100-SXM-80 GB GPUs Source catalogue Composite Sky model: # ∼ 12 . 9 × 10 6 Number of frequency channels 300 Number of time-steps 1 Station layout Randomised-sunflower Antenna element Embedded SKALA4 Number of stations 224 Simulation Station beamDuplication Unique stationRotation Number of unique stations 1 224 Beam e v aluation time [min] 12.5 2268.8 Visibility calculation [min] 8.5 9.1 Time per frequency [sec] 4.3 455.6 Total simulation time [min] 21.5 2277.9 i i c t m t ( i d F b f I d i 6 W o S f s t w w a r 6 T t r D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 esolution �f . Next, using sky catalogues and beam data from AST , the OSKAR simulator computes the N s ( N s − 1) / 2 visibilities or N s stations. During OSKAR visibility simulations, the beam omputations are managed by the harp beam library (using FAST- enerated beam data). Once all visibilities V pq ( f ) are generated, hey are binned logarithmically by baseline length b and averaged. e then apply a window function and perform the time-delay Fourier ransform to obtain the time-delay visibilities ˜ V av ( τ, b). Computation imes for both simulation tools will be based on the simulation arameters provided in Section 4 . .1 FAST he recent advancements of F AST , (Gueuning et al. 2025 ), an ccelerated direct solver based on the method of moments (MoM), ave significantly reduced the simulation time required for the omputation of all the EEPs of large, dense, and irregular arrays of dentical but complex antennas. In the case of SKA-Low, it provides ll the EEPs for a frequency point in ∼ 10 min on a standard current aptop (see Table 1 which details this benchmark). This represents significant resource and time-saving impro v ement, as it operates pproximately 25 times faster than FEKO’s Multilevel Fast Multipole ethod on a 128-core workstation. Additionally, re-simulation for ifferent station layouts with FAST takes under a minute, allowing or the optimization of station layout and analysis of MC effects t a fine resolution. The electromagnetic solver has been validated gainst several other numerical solvers, including CST, 2 WIPL 3 (Bui- an et al. 2018 ), and FEKO 4 (Gueuning et al. 2022 , 2025 ; Cavillot t al. 2024 , showing agreement to within 2–3 digits of accuracy in he far-field pattern, provided the same mesh is used. Additionally, e have carefully validated our simulations against FEKO using maller test-case stations to ensure that the chromaticity observed in ur simulations results from MC and not from approximation errors elated to the solver. .2 OSKAR SKAR is a GPU-accelerated telescope simulator specifically de- igned for the SKA (Dulwich 2020 ), utilizing the Radio Interferomet- ic Measurement Equation (RIME). The RIME framework models he contributions of the complex interferometric visibility outlined n equation ( 11 ) through a discrete sum across all visible sources NRAS 538, 31–48 (2025) See https:// www.3ds.com/ products/ simulia/ cst- studio- suite . See https:// wipl-d.com/ . See https:// altair.com/ feko . t w f T c n the sky whilst accounting for any potential instrumental effects ncluding MC. To obtain the visibilities, OSKAR requires a star atalogue to describe all the source properties and a description of he position of each station, known as the sky model and telescope odel respectiv ely. F or SKA-Low, this tool allows us to obtain he visibilities in just 4.3 s per frequency for ∼ 12 . 9 × 10 6 sources Table 2 which details this benchmark). The harp beam library is ntegrated into OSKAR and efficiently calculates station beams irectly from each antenna’s current distribution computed from AST . On an NVIDIA A100 GPU, it computes SKA-Low station eams for a 16 000-pixel grid in 1 ms. Pre-computed APs save time or identical layouts, but differing layouts require re-computation. n such cases, OSKAR simulations show station beam e v aluations ominate runtime, being nearly 1000 times slower than visibility ntegrations. N U M E R I C A L RESULTS e present the numerical results obtained using our pipeline, as utlined in Section 5 , and the simulation parameters described in ection 4 . First, we examine the variation in the BIR fall-off rate or three different station layouts, highlighting the impact of intra- tation MC. Next, we explore a foreground a v oidance scenario with he SKA telescope and discuss the size of the available detection indow . Finally , we consider a foreground removal scenario in hich the station beam is inaccurately modelled using either the AEP pproximation, corrupted or interpolated EEPs. We then attempt to emo v e the foreground power using these approximate beam models. .1 Beam impulse responses he station beams are plotted at 135 MHz the band centre, in Fig. 8 for he three different layouts illustrated in Fig. 5 : regular , sunflower , and andom. As can be expected (Clavier et al. 2014 ), at this frequency, all hree layouts exhibit a similar main beam value, around 45 dBV, along ith comparable levels for the first three sidelobes. The beam pattern or the regular layout exhibits large grating lobes (nearly 30 dBV). he beam for the random layout shows higher secondary sidelobes ompared to those of both the sunflower and regular layouts. https://www.3ds.com/products/simulia/cst-studio-suite https://wipl-d.com/ https://altair.com/feko Mutual coupling effects in 21-cm experiments 39 Figure 8. H-plane cut of the SKA-Low station beams for the ‘regular’, ‘sunflower’, and ‘random’ layouts at 135 MHz. e c e d O b A r o i r 1 r b o ( o i d t t ( a s 1 l l r r r s s f e s i e t t a 7 t o 6 E s a B t s T p τ I a a 1 s 2 s w f f i t a A w s i l s t b i ( p i l d t t r e 5 Each cosmological cube, generated at an integer redshift z, contains 21- cm brightness temperature samples T 21 , defined on 512 3 pixel grid, with each pixel having a side length of 3 cMpc. To obtain the fiducial 21-cm signal power spectrum amplitude, we begin by calculating the wavevector k , given by equation ( 2 ). A three-dimensional discrete Fourier Transform is then computed on each simulation cube to obtain, T 21 and the power spectrum calculated using, ˜ P 21 ( k, z) = | T 21 | 2 × V 2 pix /V cube av eraged o v er 100 logarithmically spaced k-bins in the range [7 × 10 −3 , 2]. 6 Information regarding the astrophysical parameters used to model the 21-cm cubes with 21cm SPACE can be found in section 2.2 of O’Hara et al. ( 2024 ). D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 Fig. 9 illustrates, in grey, all 256 individual EEPs for feed-H, v aluated at zenith and sampled at a 1 MHz rate, along with their orresponding BIRs. The solid black line represents the average mbedded element pattern (AEP) across the elements, while the ashed line represents the isolated element pattern (IEP) of SKALA4. n the one hand, the IEP appears smooth across the frequency and, suggesting a nominal operation of the SKALA4 antenna. s a consequence, the BIR of the isolated antenna decays quickly eaching −50 dB at 300 ns, with the main lobe width matching that f the window function ˜ W . On the other hand, prominent notches n the EEPs can be observed at 124 and 132 MHz in the frequency esponse of both the regular and sunflower layouts. The drop at 24 MHz results from destructive interference between the field adiated by currents on the active antenna and the field radiated y currents induced on the passive antennas due to MC. This effect ccurs when the inter-element distance approaches one wavelength Cumner et al. 2024a ). The depth and width of this notch have been bserved to increase with the array regularity, i.e. the redundancy n baseline length distribution (Anstey et al. 2024 ). In the delay omain, these MC effects cause the first lobe of the BIRs of the EEPs o be wider for both the regular and sunflower layouts (extending o 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, t delays abo v e 215 ns, suggests the presence of sharp, small-scale pectral variations in the BTFs that vary at a quick rate (faster than MHz). These variations appear to be less pronounced at the AEP evel (or equivalently here in the AP) in the sunflower and random ayouts. Specifically, the BIR of the AEP is around −20 dB for the egular layout, while it drops to −30 / 40 dB in the sunflower and andom layouts. This attenuation is likely due to the irregular ar- angement of antennas, which may help minimise MC effects at these cales. We analyse now the BTF and BIR of the AP for the random layout, ampled at 100 kHz, and examine the impact of cubic interpolation rom a coarse spectral grid, with steps of �f = 781 kHz, to match the xpected SKA-Low station channel width (Trott & Wayth 2016 ). As hown in Fig. 10 (left), the zoomed region highlights rapid variations n the BTF that outpace the 781 kHz rate, resulting in interpolation rrors of up to 0 dBV in the AP (roughly 2 significant digits in he E-field pattern, respectively). In Fig. 10 (right), this results in ime-delay residual errors around 10 dB, leading to 3 per cent errors t short delays, and divergence from exact values at approximately 00 ns. These interpolated EEPs will be applied in our pipeline for he foreground removal scenario in Section 6.3 to e v aluate the impact f interpolation errors on the delay power spectrum. .2 For egr ound avoidance arly studies suggest that the sensitivity of the SKA-Low tele- copeshould be sufficient to detect the 21-cm power spectrum t scales defined by k ∈ [0 . 1 , 0 . 33] cMpc −1 (Cohen, Fialkov & arkana 2018 ; Barkana et al. 2023 ). Ho we ver, this section illustrates, hat MC effects may compromise these predictions. To characterize foreground contamination, we define the delay pread, denoted as τmin ( b), also known as the foreground spill over. his value represents the smallest delay at which the 21-cm signal ower spectrum ˜ P 21 surpasses the foreground power spectrum ˜ P F : min ( b) = min [ τ | ˜ P F ( b, τ ) < ˜ P 21 ( b, τ ) ] . (21) n our analysis, we calculate a fiducial 21-cm signal power spectrum mplitude (Lanman et al. 2020 ) of ˜ P 21 ≈ 114 2 mK 2 h −3 Mpc 3 t the band’s central redshift z = 9 . 52. 5 This was derived using .5 cGpc 21-cm brightness temperature cubes, generated with the emi-numerical code 21cmSPACE (Visbal et al. 2012 ; Fialkov et al. 013 ). 6 Using the telescope model, Fig. 11 shows the delay power pectrum for an observation at 09:31:30 UTC on 2000 January 1, ith the GC along the horizon. The simulation spans 300 channels rom 120 to 150 MHz, using an identical (un-rotated) random layout or all 224 stations of the SKA-Low core. The results are presented n two cases. In the first case, shown in the top plot (without MC), he AP is approximated using the IEP of SKALA4 multiplied by the rray factor (AF). The AF is given by the equation F ( f , ̂ s ) = N a ∑ n = 1 e jk 0 ̂ s ·b ant ,np , (22) here unit coefficients c np = 1 are assumed. In the second case, hown in the bottom plot (with MC), the exact AP is used, which ncorporates all EEPs as defined in equation ( 8 ). The dashed white ine marks the delay spread, τmin , outlining a 21-cm detection window hown as the upper left black triangular region in the top panel. In his region, the foreground signals are sufficiently attenuated, falling elow the estimated level of the 21-cm signal. The solid black line ndicates the horizon limit, defined by the maximum geometric delay equation 17 ) for a source at the horizon. A bright band of power, resent across all baseline lengths and confined to short delays, is dentified as the intrinsic foreground. Additionally, along the horizon imit, a bright limb is observed, caused by an apparent increase in flux ensity (Th yag arajan et al. 2015 ). The dashed black line indicates he beam limit, determined by the primary beamwidth of 4 . 2 ◦ at he central frequency, as illustrated in Fig. 8 . The slower decay ate of the BIR of EEPs, compared to that of the IEP, observed arlier in Fig. 9 translated into similar effects in the delay power MNRAS 538, 31–48 (2025) 40 O. S. D. O’Hara et al. M Figure 9. The beam transfer function (left) and beam impulse response (right) for a SKA-Low station with regular (top), sunflower (middle), and random (bottom) layouts, e v aluated at zenith. The results, shown for all EEPs (solid grey lines), the IEP (dashed black lines), and the AEP (solid black lines), are e v aluated o v er 30 frequenc y channels ranging from 120 to 150 MHz. s b 7 A d t v p a e s t ( a t o p U U d ( H e s d i c D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 pectrum plot shown in Fig. 11 . The delay spread for the smallest aseline ( ∼ 40 m) generated with the AF approximation extends to 00 ns, while it extends beyond the maximum delay for the exact P. The chromaticity introduced by MC significantly increases the elay spread across all baselines, causing the foreground emission o leak and completely occlude the 21-cm detection window. As alidation of our absolute power spectrum levels, we compare the eak delay power spectrum without MC effects (top plot in Fig. 11 ) gainst results from PAPER (Pober et al. 2013 ), MWA (Beardsley t al. 2016 ), HERA (Abdurashidova et al. 2023 ), and SKA-Low imulation (Trott & Wayth 2016 ). The agreement confirms consis- ent diffuse emission delay spreads and comparable peak values around 2 10 14 mK 2 h −3 Mpc 3 ), despite differences in bandwidth nd baseline binning. We will next investigate how observation NRAS 538, 31–48 (2025) imes, station layouts, and random station rotations influence the spill v er. Using the same telescope model, Fig. 12 compares the delay ower spectrum for two observations on 2000 January 1: at 09:31:30 TC, when the Galactic plane is along the horizon, and at 12:31:30 TC, when the GC is below the horizon. We decompose the elay power spectrum into separate contributions from point sources using the GLEAM catalogue) and diffuse emission (using the ASLAM map). The results show that spill o v er is larger for diffuse mission ( ∼ 10 10 mK 2 h −3 Mpc 3 at 2000 ns) compared to point ources ( ∼ 10 6 mK 2 h −3 Mpc 3 at 2000 ns). This difference arises ue to the increased chromaticity of the diffuse emission, resulting n an apparent broader delay spectrum S( τ ). The diffuse emission ontribution also tapers off as baseline length increases. Fig. 12 Mutual coupling effects in 21-cm experiments 41 Figure 10. The beam transfer function (BTF, left) and beam impulse response (BIR, right) simulated at 100 kHz (solid black line) and a cubic spline interpolated to 100 kHz from data sampled at 781 kHz (dashed black line) for a random layout SKA-Low station, evaluated at zenith. The BIR is truncated to 1500 ns. The BTF residual was initially calculated in volts and subsequently converted to dBV, whereas the BIR residual represents the Fourier transform of the power difference. Figure 11. The delay power spectrum at 09:31:30 UTC on 2000 January 1, without MC (top) and with MC (bottom), for the un-rotated random layout of the SKA-Low core across 300 channels from 120 to 150 MHz. The solid black line represents the horizon limit, the dashed black line indicates the beam limit, and the dashed white line shows the foreground spill over. F ore ground contamination induced by MC completely masks the accessible Fourier modes of the detection window (EoR window). c c p i e p s i l i i t c a G w u i i h s u a l s l l c l M r r l o f l r s 8 a D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 learly shows that, at both observation times, diffuse emission is oncentrated around the horizon line, while point source emission redominantly originates from emissions in the main beam. Close nspection of the left-hand panels (GLEAM) exhibits faint lines (an xample is denoted by a black arrow) of increased power running arallel to the horizon line. Each line originates from an increase in ensitivity due to the fourth sidelobe of the random layout appearing n Fig. 8 . Furthermore, when the GC is below the horizon (bottom eft plot), a notable increase in power is observ ed relativ e to when it s abo v e the horizon (top left plot), e xtending to the beam limit. This ncrease is attributed to a higher source count and flux density within he primary beamwidth, as recorded by the GLEAM catalogue. In the entral panels showing diffuse emission at both observation times, decrease in diffuse power along the horizon is observed when the C is below the horizon, compared to when it is at grazing angles, hich explains the increased foreground spill over when the GC is p. In Fig. 13 , we analyse how the power spectrum varies with dentical station layout (regular , sunflower , and random, as shown n Fig. 5 ) for the first observation time with the GC abo v e the orizon. Given that the main beam and first sidelobe levels are imilar across layouts (see Fig. 9 ), po wer le vels remain consistent p to and slightly beyond the beam limit for both point source nd diffuse emission. As the ele v ation angle increases to 20 ◦, the arger sidelobes in the random layout compared to those in the unflower or regular layouts result in increased point source power evels within the foreground wedge, between the horizon and beam imit, and with a more prominent sidelobe line corresponding to the ontribution of bright point sources. The grating lobe of the regular ayout appears as a brighter region near the horizon line. Regarding C-specific effects, in line with the discussions in Section 6.1 , the egular and sunflower layouts exhibit a broader intrinsic foreground e gion, e xtending to approximately 200 ns. In contrast, the random ayout, with its narrower first lobe width, shows a response extending nly to around 150 ns. The delay spread is noticeably much larger or the regular layout for both point source and diffuse emission. The eakage of diffuse emission appears similar for the sunflower and andom layouts. The spill o v er of point source emission is ho we ver lightly lower for the sunflower layout (around 10 9 mK 2 h −3 Mpc 3 at 00 ns) compared to the random layout (around 10 10 mK 2 h −3 Mpc 3 t 800 ns). Overall, for all layouts, the entirety of the detection MNRAS 538, 31–48 (2025) 42 O. S. D. O’Hara et al. M Figure 12. The delay power spectrum for different sky map compositions (GLEAM, Haslam 408 MHz, and the composite sky map) is shown for two snapshots: at 09:31:30 UTC on 2000 January 1, when the Galactic plane is along the horizon, and at 12:31:30 UTC on 2000 January 1, when the Galactic centre is below the horizon. Each observation was simulated with MC, using an un-rotated random station layout and the SKA-Low core across 300 channels from 120 to 150 MHz. w F r ( f a b b t p a t s ( H 6 T 2 w p 2 i ( t b d f t r e a c R w p r s b t t f s l a A D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 indo w (EoR windo w), up to k ‖ ∼ 2 h −1 Mpc (see bottom row of ig. 13 ) is contaminated. Each of the 224 SKA-Low core stations is now assigned a unique andom rotation angle, uniformly distributed between 0 ◦ and 360 ◦ see Fig. 6 ). This rotation is expected to attenuate contributions rom sources in the sidelobes for station layouts with more irregular zimuthal distributions. Practically, this means that each station beam ecomes unique, requiring the computation of 224 individual station eams instead of just one, which significantly increases computation ime in OSKAR, as shown in Table 2 . Fig. 14 compares the delay ower spectrum slices for un-rotated and rotated random layouts t given baseline lengths b = [59 , 198 , 882] m. While the size of he intrinsic foreground region remains unchanged, the rotated case hows a roughly twofold decrease in power beyond the beam limit sidelobe region), particularly along the brightened horizon limb. o we ver, leakage at longer delays remains unaffected. .3 For egr ound r emo v al he EoR program for the SKA-Low telescope (Koopmans et al. 015 ) comprises three major experiments: power spectrum analysis, hich uses visibilities averaged to around 100 kHz and 5 s; tomogra- hy, relying on imaging data formed from these visibilities; and the 1-cm Forest experiment, which employs high spectral resolution mage cubes with finer spectral channels of ∼ 4 kHz. Several studies Chapman et al. 2016 ; Sims et al. 2016 ; Li et al. 2019 ) suggest NRAS 538, 31–48 (2025) hat a combination of removal and a v oidance approaches will likely e employed for 21-cm science with SKA-Low. As previously iscussed, the fore ground remo val approach aims to eliminate the oreground contribution from measured visibilities, leaving uncon- aminated 21-cm signals. Achieving accurate foreground removal equires a precise model of the instrumental response on the sky. To stimate foreground residuals, we conduct simple removal tests using pproximate station beam models. The residuals in the visibilities are omputed as: es V pq ( f ) = V pq ( f ) − V pq, app ( f ) . (23) here ˜ V pq, app ( f ) represents the visibility computed with the ap- roximate beam models. The power spectra associated with these esiduals are then computed using equations ( 14 ) and ( A7 ). In this tudy, we consider three approximations to the beam model: station eams corrupted by random noise, the AEP approximation, and he spectral interpolation of EEPs from coarser channels. These ests allow us to assess the impact of beam model inaccuracies on ore ground remo val and the resulting residuals in the delay power pectrum. In the following examples, we consider un-rotated station ayouts. The first test consists of adding Gaussian white noise to the APs s follows: P p, noise ( f , ̂ s ) = AP p ( f , ̂ s ) + N q ( f , ̂ s ) , (24) Mutual coupling effects in 21-cm experiments 43 Figure 13. Delay power spectrum across different sky map compositions (GLEAM, Haslam 408 MHz, and composite) and un-rotated station layouts (regular, sunflower, and random) including MC, for a snapshot at 09:31:30 UTC on 2000 January 1, with the SKA-Low core across 300 channels from 120 to 150 MHz. w m o d 1 t n a 2 S r E r t ( t s k t 8 A a A f A F [ t r t r A t A D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 here N q represents complex-valued Gaussian noise with zero ean and variance σ 2 . This noise is a multi v ariable function f frequency f , the N pix pixel directions ˆ s , and the station in- ices p. We assume four different values for the variances, σ 2 = 0 0 , 10 −2 , 10 −4 , 10 −6 . For a signal coming from zenith, where he beam value is around 45 dBV, this corresponds to signal-to- oise ratios SNR = 45 , 65 , 85 , 105 dB. Since 20 dB corresponds to single digit of accuracy in the E-field, the noise is approximately , 3, 4, and 5 digits below the received voltage for the respective NR values. Fig. 15 illustrates a selection of delay power spectrum esiduals across baseline bins with lengths of b = [59 , 198 , 882] m. ach additional significant figure of accuracy in the station model oughly corresponds to a reduction of two additional digits in he 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, o achie ve ef fecti ve foreground remov al belo w our fiducial 21-cm ignal, a 4–5-digit accuracy in the station beam model is required to eep the foreground residuals below the 114 2 mK 2 h −3 Mpc 3 level of he 21-cm signal. Since the AP is the sum of the EEPs (see equation ), the noise variance on the EEPs that results in similar noise on the P is σ 2 / 256. We estimate the required accuracy on the EEPs to be round 3–4 significant digits. Then, we analyse the residuals resulting from approximating the P with the AEP multiplied by the AF (Wijnholds et al. 2019 ) as ollows: P p ( f , ̂ s ) ≈ EEP av,p ( f , ̂ s ) AF ( f , ̂ s ) . (25) ig. 16 shows slices of the delay power spectrum for b = 59 , 198 , 882] m, comparing the power spectrum using the exact AP, he AEP approximation, and the residual after attempting foreground emoval in the visibilities. The AEP approximation is accurate within he intrinsic foreground region (up to 150 ns) with power spectrum esiduals on the order of 1 per cent to 0.1 per cent at short delays. The EP approximation is off by seven orders of magnitude, preventing he residuals from reaching the level of our fiducial 21-cm signal. s shown in Fig. 15 , the residuals from the AEP approximation are MNRAS 538, 31–48 (2025) 44 O. S. D. O’Hara et al. M Figure 14. A delay power spectrum slice for b = [59 , 198 , 882] m was generated using 224 identical (un-rotated) stations and 224 uniquely rotated stations, both arranged in a random layout. o b N t t l � S t w a c s o a E 1 1 t h I A A a s r t a s 7 T a t s t c i t o f e f t f F r D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niv ne order of magnitude higher than those obtained by using station eams only accurate to 2 digits (assuming uncorrelated noise errors). ote that the tail of the estimated power spectrum decays faster than he true delay spectrum. This observation is consistent with the lower ail of the BIR associated with the AEP in Fig. 10 . In the final test, we examine the residuals resulting from interpo- ating the EEPs from a coarser spectral grid with a sampling rate of f = 781 kHz to the fine sample rate of 100 kHz as discussed in ection 6.1 . This coarse rate corresponds to the separation between he 384 coarse channels of SKA-Low, co v ering the 300 MHz band- NRAS 538, 31–48 (2025) igure 15. Delay power spectrum slices of constant b = [59 , 198 , 882] m, were s andom station layout. The residuals were computed with respect to the visibilities idth, which is formed at each station before applying beamforming nd/or calibration weights (Comoretto et al. 2020 ). Each coarse hannel contains around 200 fine 4 kHz channels. The calibration trategy outlined in Trott & Wayth ( 2016 ) proposes using a low- rder polynomial fit for the calibration coefficients across three djacent coarse channels. Inspired by this work, we interpolate the EPs, simulated at the centre of each coarse channel within the 20 –150 MHz band of interest, and then upsample these EEPs to a 00 kHz resolution. The AP is then formed using equation ( 8 ) from hese interpolated EEPs, resulting in an approximate beam model. We ave plotted the residuals for the three baseline lengths in Fig. 17 . n agreement with the 30 dB residual observed in the BIR of the P (Fig. 10 ), the power spectrum residuals are around 1 per cent. s shown in Fig. 17 , the coarse interpolation of the EEPs leads to n estimated delay power spectrum that decays faster than the true pectrum at long delays due to spectral smoothing. Considering the esults of previous tests, another conclusion from this experiment is hat reducing residuals by employing different EEP models for each ntenna is ef fecti ve only when these models are computed at high pectral resolutions. C O N C L U S I O N his paper analyses the impact of Mutual Coupling (MC) between ntennas on the time-delay power spectrum of the SKA-Low elescope. Electromagnetic simulations reveal significant and rapid pectral variations in the station beam directi vity, dri ven by energy ravelling and resonating between antennas due to MC. For 21- m science, these MC effects cause substantial foreground leakage nto the detection window, compromising foreground a v oidance echniques. Highly accurate station beam models (within 4–5 orders f magnitude) are required to ef fecti v ely remo v e fore grounds. Brute- orce approximations, such as these based on randomizing MC ffects or on interpolating coarsely beam responses, will likely all short. Although this study focuses on the direct e v aluation of ime-delay power spectrum, we expect image-based analyses, (see or instance Morales et al. 2019 ), to require similar beam model imulated with an AP accuracy limited to 2, 3, 4, and 5 digits for an un-rotated generated with the exact AP. ersity of C am bridge user on 28 January 2026 Mutual coupling effects in 21-cm experiments 45 Figure 16. Delay power spectrum slices of constant b = [59 , 198 , 882] m for the exact AP, the AEP approximation and their absolute residual (Shaded grey region) for an un-rotated random station layout. Figure 17. Delay power spectrum slices of constant b = [59 , 198 , 882] m using the exact AP, the AP with spline interpolated EEPs (upsampled from 781 to 100 kHz) and their absolute residual (Shaded grey region) for an un-rotated random station layout. a N t m m a a s 2 1 p e b t p t T i d i D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 ccuracy for image cleaning before power spectrum reconstruction. ew calibration algorithms and measurement techniques will have o be developed to account for MC effects and enable precise beam odelling. The analysis in this paper is made possible through the develop- ent of two key numerical tools: (i) The FAST electromagnetic solver: All EEPs were simulated cross 300 frequency points, completing in roughly two days on standard laptop, with only a few additional hours needed for re- imulation in cases of new station layouts. (ii) The OSKAR visibility simulator: The simulator processes the 4 976 visibilities in the SKA core across all frequencies in just 0 min using 4 A100 GPUs and including diffuse emissions and oint sources across approximately 10 million sky pixels. The beam v aluation method embedded into OSKAR directly calculates station eams from the aperture current distribution, achieving evaluation imes around 2 s per station beam for the configuration abo v e. Using these simulation tools, we fully inte grated e xact array atterns (APs) directly into visibility calculations, enabling precise racking of MC effects in the time-delay power spectrum analysis. his yielded the following key findings regarding intra-station MC mpacts: (i) Delay spread broadening: MC effects significantly extend the elay spread of the power spectrum. For EoR science, this can result n the complete obscuration of the detection window. MNRAS 538, 31–48 (2025) 46 O. S. D. O’Hara et al. M r e c s r e s w t ( o p r t a t d m e e d d w t A T c a p b t S a h a b S S J b a C D T r R A A A B B B B B B B B B B B B B C C C C C C C C C C C D d d d d d d D D D E E D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 (ii) Layout-dependent MC effects: while offering beneficial beam edundanc y, re gular station layouts are more prone to resonant wave ffects that lead to sudden power drops across the field of view. In ontrast, randomized layouts reduce these issues, providing more table beam behaviour and a narrower delay spread in the impulse esponse. (iii) Residuals with AEP approximation: When approximating xact APs by using the AEP multiplied by the array factor, power pectrum residuals reach peak values around 10 11 mK 2 h −3 Mpc 3 ithin the foreground wedge. This is about seven orders of magni- ude higher than the expected level of the 21-cm power spectrum ∼ 10 4 mK 2 h −3 Mpc 3 ). (iv) Residuals with corrupted station beams: Testing various levels f accuracy (2, 3, 4, and 5 significant digits) resulted in residual ower spectrum values of 10 9 , 10 7 , 10 5 , and 10 3 mK 2 h −3 Mpc 3 , espectively. This demonstrates the need to model the station beams o 4 –5 significant digits. (v) Residuals with interpolated EEPs: Interpolating the EEPs from coarse (781 kHz) to a finer (100 kHz) spectral grid before forming he AP results in power spectrum residuals around 1 per cent at short elays, highlighting the need for high spectral resolution in EEP odels to reduce interpolation errors. This work highlights the need to further develop station beam v aluation methods in processing pipelines, incorporate additional lectromagnetic effects, and assess their impact across observational omains (baseline, delay, fringe rate, and image) to filter corrupted ata and refine calibration algorithms. Future work should examine hether certain parts of the band are less affected by MC and extend hese investigations to include image-based methodologies. C K N OW L E D G E M E N T S he first two authors, Oscar S.D. O’Hara and Quentin Gueuning, ontributed equally to this work and should be considered joint first uthors. The authors thank Dr. Karel Adamek for his assistance in arallelizing the HARP beam library, Dr. Vladislav Stolyarov for enchmarking it. We also thank Pr. Christophe Craeye for supporting he development of the fast electromagnetic solver, Dr. Maciej yralek and Dr. Robert Laing for their support. We also thank the nonymous re vie wers for their comments and suggestions, which av e impro v ed the quality of this paper. This research was supported by ESA, NPL, and the Science nd Technology Facilities Council (STFC). OOH was supported y grant number G109464, ‘The Design of Highly Sensitive EM ensors for Space Applications’, QG and FD by grant number T/W00206X/1, EdLA by grant number ST/V004425/1 and DA and C were supported by grant number ST/X00239X/1. AF is supported y a Royal Society University Research Fellowship #180523. JD cknowledges support from the Boustany Foundation and Cambridge ommonwealth Trust in the form of an Isaac Newton Studentship. ATA AVA ILA BILITY he data and software underlying this article will be shared on easonable request to the corresponding authors. E FEREN C ES bdurashidova T. H. C. Z. et al., 2023, ApJ , 945, 124 miri M. et al., 2022, ApJS , 261, 29 NRAS 538, 31–48 (2025) nstey D. et al., 2024, in 2024 18th European Conference on Antennas and Propagation (EuCAP). IEEE, New Yerk, p. 1 andura K. et al., 2014, in Proc. SPIE Conf. Ser. Vol. 9145, Ground-based and Airborne Telescopes V. SPIE, Bellingham, p. 20 arkana R. , Fialkov A., Liu H., Outmezguine N. J., 2023, Phys. Rev. D , 108, 063503 arry N. , Hazelton B., Sulli v an I., Morales M. F., Pober J. C., 2016, MNRAS , 461, 3135 eardsley A. et al., 2016, ApJ , 833, 102 ird T. S. , 2021, Mutual Coupling Between Antennas , 177 olli P. , Davidson D. B., Braun R., Di Ninni P., Ung D., 2022a, 2022 International Conference on Electromagnetics in Advanced Applications (ICEAA). IEEE, New York, p. 200 olli P. , Davidson D. B., Bercigli M., Ninni P. D., Labate M. G., Ung D., Virone G., 2022b, J. Astron. Telesc. Instrum. Syst. , 8, 011017 olli P. , Bercigli M., Ninni P. D., Mezzadrelli L., Virone G., 2022c, J. Astron. Telesc. Instrum. Syst. , 8, 011023 onaldi A. , Brown M. L., 2015, MNRAS , 447, 1973 raun R. , Bourke T. L., Green J. A., Keane E. F ., W agg J., 2015, Proc. Sci., Advancing Astrophysics with the Square Kilometre Array. SISSA, Trieste, PoS(AASKA14)174 owman J. D. , Morales M. F., Hewitt J. N., 2006, ApJ , 638, 20 owman J. D. et al., 2013, Publ. Astron. Soc. Aust. , 30, e031 ui-Van H. , Abraham J., Arts M., Gueuning Q., Raucy C., Gonz ́alez-Ovejero D., de Lera Acedo E., Craeye C., 2018, IEEE T. Antenn. Propag. , 66, 1805 avillot J. , Tihon D., Gueuning Q., de Lera Acedo E., Craeye C., 2024, IEEE T. Antenn. Propag. , 72, 7560 hapman E. , Zaroubi S., Abdalla F. B., Dulwich F., Jeli ́c V., Mort B., 2016, MNRAS , 458, 2928 hokshi A. , Barry N., Line J. L. B., Jordan C. H., Pindor B., Webster R. L., 2024, MNRAS , 534, 2475 hose M. , Conradie A. S., Cilliers P. I., Botha M. M., 2023, Microw. Opt. Techn. Lett. , 65, 2359 la vier T. , Raza vi-Ghods N., Glineur F., Gonz ́alez-Ovejero D., de Lera Acedo E., Craeye C., Alexander P., 2014, IEEE T. Antenn. Propag. , 62, 1596 ohen A. , Fialkov A., Barkana R., 2018, MNRAS , 478, 2193 omoretto G. et al., 2020, in Proc. SPIE. Conf. Ser. Vol. 11445, Ground-based and Airborne Telescopes VIII. SPIE, Bellingham, p. 14 onradie A. S. , Chose M., Cilliers P. I., Botha M. M., 2023, IEEE T. Antenn. Propag. , 71, 5199 raeye C. , Gonz ́alez-Ovejero D., 2011, Radio Sci. , 46, 1 umner J. et al., 2024a, 2024 International Conference on Electromagnetics in Advanced Applications (ICEAA). IEEE, Lisbon, Portugal, p. 171 umner J. , Pieterse C., De Villiers D., de Lera Acedo E., 2024b, MNRAS , 531, 4734 avidson D. B. et al., 2020, 2020 XXXIIIrd General Assembly and Scientific Symposium of the International Union of Radio Science. IEEE, New York, p. 1 e Lera Acedo E. , 2012, 2012 International Conference on Electromagnetics in Advanced Applications. IEEE, New York, p. 353 e Lera Acedo E. , Arts M., Craeye C., Van H. B., 2016, 2016 International Conference on Electromagnetics in Advanced Applications (ICEAA). IEEE, New York, p. 683 e Lera Acedo E. et al., 2017, MNRAS , 469, 2662 e Lera Acedo E. , Pienaar H., Fagnoni N., 2018, in 2018 International Conference on Electromagnetics in Advanced Applications (ICEAA). p. 636 e Lera Acedo E. et al., 2022, Nat. Astron. , 6, 984 e Villiers M. S. , Cotton W. D., 2022, AJ , 163, 135 eBoer D. R. et al., 2017, PASP , 129, 045001 ewdney P. E. , Hall P. J., Schilizzi R. T., Lazio T. J. L. W., 2009, Proc. IEEE , 97, 1482 ulwich F. , 2020, OSKAR 2.7.6 . ngheta N. , Murphy W. D., Rokhlin V., Vassiliou M. S., 1992, IEEE T. Antenn. Propag. , 40, 634 wall-Wice A. , Dillon J. S., Liu A., Hewitt J., 2017, MNRAS , 470, 1849 http://dx.doi.org/10.3847/1538-4357/acaf50 http://dx.doi.org/10.3847/1538-4365/ac6fd9 http://dx.doi.org/10.1103/PhysRevD.108.063503 http://dx.doi.org/10.1093/mnras/stw1380 http://dx.doi.org/10.3847/1538-4357/833/1/102 http://dx.doi.org/10.1002/9781119565048.ch7 http://dx.doi.org/10.1117/1.JATIS.8.1.011017 http://dx.doi.org/10.1117/1.JATIS.8.1.011023 http://dx.doi.org/10.1093/mnras/stu2601 http://dx.doi.org/10.1086/498703 http://dx.doi.org/10.1017/pas.2013.009 http://dx.doi.org/10.1109/TAP.2018.2806222 http://dx.doi.org/10.1109/TAP.2024.3445132 http://dx.doi.org/doi.org/10.1093/mnras/stw161 http://dx.doi.org/10.1093/mnras/stae2264 http://dx.doi.org/10.1002/mop.33696 http://dx.doi.org/10.1109/TAP.2013.2284816 http://dx.doi.org/10.1093/mnras/sty1094 http://dx.doi.org/10.1109/TAP.2023.3268726 http://dx.doi.org/10.1029/2010RS004518 http://dx.doi.org/10.1093/mnras/stae1475 http://dx.doi.org/10.1093/mnras/stx904 http://dx.doi.org/ 10.1038/s41550-022-01709-9 http://dx.doi.org/10.3847/1538-3881/ac460a http://dx.doi.org/10.1088/1538-3873/129/974/045001 http://dx.doi.org/10.1109/JPROC.2009.2021005 http://dx.doi.org/10.5281/zenodo.3758491 http://dx.doi.org/10.1109/8.144597 http://dx.doi.org/10.1093/mnras/stx1221 Mutual coupling effects in 21-cm experiments 47 F G G G G H H J J J J J K K K K K L L L L L L L L M M M M M M N O P P P P P P P P R R R S S S S S T T T T T T v v V V W W Y Z A H s m T f b o ξ w s i c t P A T ( s d l l r s ( D ow nloaded from https://academ ic.oup.com /m nras/article/538/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 ialkov A. , Barkana R., Pinhas A., Visbal E., 2013, MNRAS , 437, L36 an H. et al., 2023, A&A , 669, A20 ueuning Q. , Craeye C., Colin-Beltran E., de Lera Acedo E., 2015, 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA). IEEE, New York, p. 1214 ueuning Q. , de Lera Acedo E., Brown A. K., Craeye C., 2022, IEEE T. Antenn. Propag. , 70, 9511 ueuning Q. , Acedo E. d. L., Brown A. K., Craeye C., O’Hara O., 2025, IEEE T. Antenn. Propag. , 1 amaker J. , Bregman J., Sault R., 1996, A&AS , 117, 137 ogg D. W. , 2000, preprint ( arXiv:astro-ph/9905116 ) eli ́c V. et al., 2008, MNRAS , 389, 1319 onas J. , MeerKAT Team , 2016, Proc. Sci., The MeerKAT Radio Telescope . SISSA, Trieste, PoS(MeerKAT2016)001 ones R. C. , 1941, J. Opt. Soc. Am. , 31, 488 osaitis A. T. , Ewall-Wice A., Fagnoni N., de Lera Acedo E., 2022, MNRAS , 514, 1804 oseph R. C. , Trott C. M., Wayth R. B., Nasirudin A., 2019, MNRAS , 492, 2017 azemi S. , Yatawatta S., Zaroubi S., Lampropoulos P., De Bruyn A., Koopmans L., Noordam J., 2011, MNRAS , 414, 1656 ern N. S. , Parsons A. R., Dillon J. S., Lanman A. E., Fagnoni N., de Lera Acedo E., 2019, ApJ , 884, 105 ern N. S. et al., 2020, ApJ , 888, 70 ildal P.-S. , 2015, Foundations of antenna engineering: a unified approach for line-of-sight and multipath. Artech House oopmans L. V. E. et al., 2015, Proc. Sci., The Cosmic Dawn and Epoch of Reionization withthe Square Kilometre Array. SISSA, Trieste, p. 001 anman A. E. , Pober J. C., Kern N. S., de Lera Acedo E., DeBoer D. R., Fagnoni N., 2020, MNRAS , 494, 3712 i W. et al., 2019, ApJ , 879, 104 ine J. L. B. et al., 2020, Publ. Astron. Soc. Aust. , 37, e027 iu A. , Shaw J. R., 2020, PASP , 132, 062001 iu A. , Tegmark M., Bowman J., Hewitt J., Zaldarriaga M., 2009, MNRAS , 398, 401 iu A. , Parsons A. R., Trott C. M., 2014a, Phys. Rev. D , 90, 023018 iu A. , Parsons A. R., Trott C. M., 2014b, Phys. Rev. D , 90, 023019 udwig A. , 1973, IEEE T. Antenn. Propag. , 21, 116 aaskant R. , Mittra R., Tijhuis A., 2008, IEEE T. Antenn. Propag. , 56, 3440 ertens F. G. et al., 2020, MNRAS , 493, 1662 itchell D. A. , Greenhill L. J., Wayth R. B., Sault R. J., Lonsdale C. J., Cappallo R. J., Morales M. F., Ord S. M., 2008, IEEE J. Sel. Top. Signal Process. , 2, 707 orales M. F. , Hewitt J., 2004, ApJ , 615, 7 orales M. F. , Hazelton B., Sulli v an I., Beardsley A., 2012, ApJ , 752, 137 orales M. F. , Beardsley A., Pober J., Barry N., Hazelton B., Jacobs D., Sulli v an I., 2019, MNRAS , 483, 2207 asirudin A. , Prelogovic D., Murray S. G., Mesinger A., Bernardi G., 2022, MNRAS , 514, 4655 ’Hara O. S. , Dulwich F., de Lera Acedo E., Dhandha J., Gessey-Jones T., Anstey D., Fialkov A., 2024, MNRAS , 533, 2876 aonessa F. , Ciorba L., Kyriakou G., Bolli P., Virone G., 2023, IEEE Antennas Wirel. Propag. Lett. , 22, 2735 arsons A. R. et al., 2010, AJ , 139, 1468 arsons A. , Pober J., McQuinn M., Jacobs D., Aguirre J., 2012a, ApJ , 753, 81 arsons A. R. , Pober J. C., Aguirre J. E., Carilli C. L., Jacobs D. C., Moore D. F., 2012b, ApJ , 756, 165 arsons A. R. et al., 2014, ApJ , 788, 106 arsons A. R. , Liu A., Ali Z. S., Cheng C., 2016, ApJ , 820, 51 ober J. C. et al., 2013, ApJ , 768, L36 ober J. et al., 2016, ApJ , 819, 8 ath E. et al., 2024, preprint ( arXiv:2406.08549 ) aucy C. , de Lera Acedo E., Razavi-Ghods N., Gonz ́alez-Ov ejero D., Crae ye C., 2013, 2013 7th European Conference on Antennas and Propagation (EuCAP). IEEE, New York, p. 661 emazeilles M. , Dickinson C., Banday A., Bigot-Sazy M.-A., Ghosh T., 2015, MNRAS , 451, 4311 haw J. R. , Sigurdson K., Sitwell M., Stebbins A., Pen U.-L., 2015, Phys. Rev. D , 91, 083514 ims P. H. , Lentati L., Alexander P., Carilli C. L., 2016, MNRAS , 462, 3069 mirnov O. M. , 2011, A&A , 527, A106 ubrahmanyan R. , 2021a, AAVS2 Antenna Total Powers Versus Time at 110 MHz. CSIRO-CASS, Curtin-CIRA ubrahmanyan R. , 2021b, Analysis of Total Power Spectra from AAVS2 Antennas. CSIRO-CASS, Curtin-CIRA h yag arajan N. et al., 2015, ApJ , 804, 14 h yag arajan N. , Parsons A. R., DeBoer D. R., Bowman J. D., Ewall-Wice A. M., Neben A. R., Patra N., 2016, ApJ , 825, 9 ingay S. J. et al., 2013, Publ. Astron. Soc. Aust. , 30, e007 rott C. M. , Wayth R. B., 2016, Publ. Astron. Soc. Aust. , 33, e019 rott C. M. , Wayth R. B., Tingay S. J., 2012, ApJ , 757, 101 rott C. M. , de Lera Acedo E., Wayth R. B., Fagnoni N., Sutinjo A. T., Wakley B., Punzalan C. I. B., 2017, MNRAS , 470, 455 an Diepen G. , Dijkema T. J., Offringa A., 2018, Astrophysics Source Code Library, record ascl:1804.003 an Haarlem M. P. et al., 2013, A&A , 556, A2 isbal E. , Barkana R., Fialkov A., Tseliakhovich D., Hirata C. M., 2012, Nature , 487, 70 ogel H. , 1979, Math. Biosci. , 44, 179 ayth R. B. et al., 2015, Publ. Astron. Soc. Aust. , 32, 12 ijnholds S. J. , Arts M., Bolli P., Di Ninni P., Virone G., 2019, 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA). IEEE, New York, p. 0437 atawatta S. , Kazemi S., Zaroubi S., 2012, 2012 Innov ati ve Parallel Comput- ing (InPar). IEEE, New York, p. 1 arka P. , Girard J. N., Tagger M., Denis L., 2012, in Boissier S., de Lav ern y P., Nardetto N., Samadi R., Valls-Gabaud D., Wozniak H., eds, SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. p. 687 PPENDI X: T H E DELAY POWER SPECT RUM ere, we provide a brief overview of the link between the power pectrum associated with a temperature field, the statistics we aim to easure, and the visibilities, which are the observable quantities. his relationship is straightforward, enabling the instrument to unction almost as a direct probe of the signal power spectrum. We egin with an estimation of the two-point correlation function, ξ ( r ) f a brightness temperature T ( r ) as follows: ( r ′ ) ≈ 1 V • V T ( r ) T ( r + r ′ ) d r x d r y d r z , (A1) here V is a comoving volume of interest (Hogg 2000 ). The power pectrum P( k ) associated with this temperature field, expressed n [mK 2 h −3 Mpc 3 ], corresponds to the Fourier transform of the orrelation function ξ ( r ), which transforms comoving distances r o wav ev ectors k (P arsons et al. 2012a , b ): ( k ) = •∞ −∞ ξ ( r ) e −j k ·r d r x d r y d r z . (A2) ccording to Parse v al’s theorem, we have ˜ P ( k ) = | ̃ T ( k ) | 2 where ˜ ( k ) denotes the Fourier transform of the temperature field T ( r ) Liu & Shaw 2020 ). The observed temperature T is associated with a pecific point in co-moving coordinates and corresponds to an epoch efined by a particular redshift z. We now express T as a function of ocal observational coordinates. Assuming the observed volume V is ocated at zenith and occupies a sufficiently small portion of the sky elative to the beamwidth, we can align the z -axis along the line of ight and apply the following transformation: l, m, f ) = ( r x X( z) , r y X( z ) , r z Y ( z )) . (A3) MNRAS 538, 31–48 (2025) http://dx.doi.org/10.1093/mnrasl/slt135 http://dx.doi.org/ 10.1051/0004-6361/202244316 http://dx.doi.org/10.1109/TAP.2022.3177465 http://dx.doi.org/10.1109/TAP.2025.3528766 http://dx.doi.org/https://doi.org/10.1051/aas:1996146 http://arxiv.org/abs/astro-ph/9905116 http://dx.doi.org/10.1111/j.1365-2966.2008.13634.x http://dx.doi.org/10.22323/1.277.0001 http://dx.doi.org/10.1364/JOSA.31.000488 http://dx.doi.org/10.1093/mnras/stac916 http://dx.doi.org/10.1093/mnras/stz3375 http://dx.doi.org/10.1111/j.1365-2966.2011.18506.x http://dx.doi.org/10.3847/1538-4357/ab3e73 http://dx.doi.org/10.3847/1538-4357/ab5e8a http://dx.doi.org/10.1093/mnras/staa987 http://dx.doi.org/10.3847/1538-4357/ab21bc http://dx.doi.org/10.1017/pasa.2020.18 http://dx.doi.org/10.1088/1538-3873/ab5bfd http://dx.doi.org/10.1111/j.1365-2966.2009.15156.x http://dx.doi.org/10.1103/PhysRevD.90.023018 http://dx.doi.org/10.1103/PhysRevD.90.023019 http://dx.doi.org/10.1109/TAP.1973.1140406 http://dx.doi.org/10.1109/TAP.2008.2005471 http://dx.doi.org/10.1093/mnras/staa327 http://dx.doi.org/10.1109/JSTSP.2008.2005327 http://dx.doi.org/10.1086/424437 http://dx.doi.org/10.1088/0004-637X/752/2/137 http://dx.doi.org/10.1093/mnras/sty2844 http://dx.doi.org/10.1093/mnras/stac1588 http://dx.doi.org/10.1093/mnras/stae1952 http://dx.doi.org/10.1109/LAWP.2023.3290967 http://dx.doi.org/10.1088/0004-6256/139/4/1468 http://dx.doi.org/10.1088/0004-637X/753/1/81 http://dx.doi.org/10.1088/0004-637X/756/2/165 http://dx.doi.org/10.1088/0004-637X/788/2/106 http://dx.doi.org/10.3847/0004-637X/820/1/51 http://dx.doi.org/10.1088/2041-8205/768/2/L36 http://dx.doi.org/10.3847/0004-637X/819/1/8 http://arxiv.org/abs/2406.08549 http://dx.doi.org/10.1093/mnras/stv1274 http://dx.doi.org/10.1103/PhysRevD.91.083514 http://dx.doi.org/10.1093/mnras/stw1768 http://dx.doi.org/10.1051/0004-6361/201016082 http://dx.doi.org/10.1088/0004-637X/804/1/14 http://dx.doi.org/10.3847/0004-637X/825/1/9 http://dx.doi.org/10.1017/pasa.2012.007 http://dx.doi.org/10.1017/pasa.2016.18 http://dx.doi.org/10.1088/0004-637X/757/1/101 http://dx.doi.org/10.1093/mnras/stx1224 http://dx.doi.org/ 10.1051/0004-6361/201220873 http://dx.doi.org/ 10.1038/nature11177 http://dx.doi.org/https://doi.org/10.1016/0025-5564(79)90080-4 http://dx.doi.org/10.1017/pasa.2015.26 48 O. S. D. O’Hara et al. M T o z X w r s k ( t z k U F f T H ( T f d t a t P P w a r � T s T his change of variable maps co-moving coordinates ( r x , r y , r z ) to bservation coordinates ( l, m, f ). Here, X, Y depends on the redshift and are defined as Lanman et al. ( 2020 ): ( z) ≈ χ ( z) , Y ( z) ≈ c 0 (1 + z) 2 H ( z) f 21 (A4) here χ ( z) is the co-moving distance, H ( z) is the Hubble pa- ameter and f 21 is the rest-frame 21 cm frequency. The relation- hip in equation ( A3 ) also establishes a mapping between the -space and the observational baseline-delay domain: ( u, v, τ ) = k x X( z) /λ0 , k y X( z ) /λ0 , k z Y ( z )) / (2 π ). We decompose the wav ev ec- or k into components parallel and perpendicular to the line of sight ˆ : ⊥ = 2 πu λ X( z) , k ‖ = 2 πτ Y ( z) . (A5) sing the notation T ( ̂ s , f ) = T ( r ), we can express the ourier transform with respect to the coordinates ( l, m, f ) as ollows: ˜ ( τ, b ) = • T ( ̂ s , f ) e −j2 π( ̂ s ·b λ+ τf ) d ld m d f . (A6) ere, we thus have ˜ T ( k ) = ˜ T ( τ, u ) / ( X 2 Y ) for purely planar baseline b = u ). The relationship in equation ( A6 ) is fundamental because ˜ ( τ, b ) can be nearly directly estimated from the time-delay trans- NRAS 538, 31–48 (2025) Published by Oxford University Press on behalf of Royal Astronomical Society. This is an ( https://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reus orm of the instrument observable: the visibilities V ( f , b ). The main istinction is that the visibilities account for the instrument’s response hrough the BTF, A ( ̂ s , f ). When the width of the BIR, ˜ A ( ̂ s , τ ) (in both ngular ̂ s and time dimensions τ ), is significantly narrower than the emperature field ˜ T ( ̂ s , τ ), we can approximate the power spectrum as follows (Parsons et al. 2014 ): ˜ ( k ⊥ , k ‖ ) � X( z ) 2 Y ( z ) �A ‖ ̃ V ( τ, b ) | 2 , (A7) here �A is a normalization factor (a full deri v ation can be found in ppendix B of Parsons et al. 2014 ) that accounts for the instrument esponse and is e v aluated here as: A = ∫ f max f min “ | W ( f ) A ( f , ̂ s ) | 2 d l n d m d f (A8) he approximation ( A7 ) thus serves as an estimator of the power pectrum of a given signal of brightness temperature T . his paper has been typeset from a T E X/L A T E X file prepared by the author. © 2025 The Author(s). Open Access article distributed under the terms of the Creative Commons Attribution License e, distribution, and reproduction in any medium, provided the original work is properly cited. D ow nloaded from https://academ ic.oup.com /m nras/article/53 8/1/31/8010856 by U niversity of C am bridge user on 28 January 2026 https://creativecommons.org/licenses/by/4.0/ 1 INTRODUCTION 2 CONVENTIONS 3 A FORWARD MODEL OF TELESCOPE VISIBILITIES 4 NUMERICAL EXPERIMENT: SKA-LOW 5 ACCELERATED SIMULATION TOOLS 6 NUMERICAL RESULTS 7 CONCLUSION ACKNOWLEDGEMENTS DATA AVAILABILITY REFERENCES APPENDIX: THE DELAY POWER SPECTRUM