TY - JOUR AU - Fagnoni,, Nicolas AB - ABSTRACT The 21 cm hyperfine transition of neutral hydrogen offers a promising probe of the large-scale structure of the universe before and during the Epoch of Reionization (EoR), when the first ionizing sources formed. Bright radio emission from foreground sources remains the biggest obstacle to detecting the faint 21 cm signal. However, the expected smoothness of foreground power leaves a clean window in Fourier space where the EoR signal can potentially be seen over thermal noise. Though the boundary of this window is well defined in principle, spectral structure in foreground sources, instrumental chromaticity, and choice of spectral weighting in analysis all affect how much foreground power spills over into the EoR window. In this paper, we run a suite of numerical simulations of wide-field visibility measurements, with a variety of diffuse foreground models and instrument configurations, and measure the extent of contaminated Fourier modes in the EoR window using a delay-transform approach to estimate power spectra. We also test these effects with a model of the Hydrogen Epoch of Reionization Array (HERA) antenna beam generated from electromagnetic simulations, to take into account further chromatic effects in the real instrument. We find that foreground power spillover is dominated by the so-called pitchfork effect, in which diffuse foreground power is brightened near the horizon due to the shortening of baselines. As a result, the extent of contaminated modes in the EoR window is largely constant over time, except when the Galaxy is near the pointing centre. instrumentation: interferometers, methods: numerical, techniques: interferometric, dark ages, reionization, first stars 1 INTRODUCTION The Epoch of Reionization (EoR) comprises the period of the Universe’s history after recombination when the first luminous structures formed, emitting ultraviolet radiation that carved out ionized regions of the neutral intergalactic medium (IGM). The lack of bright sources and optically thick IGM during the EoR limit observational prospects through traditional high-redshift galaxy surveys and quasar absorption measurements. The hyperfine transition of neutral hydrogen, which emits a rest-frame photon with a wavelength of 21 cm, offers a promising probe of the IGM structure during the EoR and before, when neutral hydrogen was abundant in the Universe (see Pritchard & Loeb 2008; Furlanetto, Peng Oh & Briggs 2006, for a review]. As a forbidden transition, the IGM is largely transparent to 21 cm photons, so this signal is not obscured by intervening gas. As a line transition, the redshift of 21 cm emission/absorption directly corresponds to distance. In this way, mapping the distribution of redshifted 21 cm brightness can map the three-dimensional structure of neutral hydrogen during the EoR. Directly imaging structures of the EoR will likely require an instrument of the scale of the upcoming Square Kilometre Array (Furlanetto et al. 2006; Mellema et al. 2013), so current generation experiments are seeking statistical measures of the EoR signal, such as the power spectrum (these efforts are reviewed by Liu & Shaw 2020). Radio interferometers are well suited to power spectrum measurements, since they directly sample the angular Fourier transform of the sky, and can achieve high spectral resolution, which provides good sampling of line-of-sight Fourier modes. Many of these instruments feature highly redundant, compact layouts that maximize sensitivity to large angular scales, aiming to detect the EoR signal through their fine spectral resolution first (Parsons et al. 2012a). Though an actual detection of the power spectrum has yet to be made, experiments like the GMRT (Paciga et al. 2011, 2013), MWA (Bowman et al. 2013; Beardsley et al. 2016; Barry et al. 2019), PAPER (Parsons et al. 2010; Cheng et al. 2018; Kolopanis et al. 2019), and LOFAR (Patil et al. 2017) have placed upper limits on the power spectrum amplitude at several redshifts. The primary obstacle to these experiments is contamination of bright foreground sources, such as Galactic synchrotron emission, free–free emission, radio galaxies, pulsars, and supernova remnants, which are typically five orders of magnitude brighter than the expected EoR signal (Matteo et al. 2002; Oh & Mack 2003; Furlanetto et al. 2006). Fortunately, the smooth spectra exhibited by these sources distinguish them from the relatively complex structure of EoR signals, making it possible – in principle – to model and subtract them from data (Oh & Mack 2003; Liu & Tegmark 2012; Chapman et al. 2015; Carroll et al. 2016; Sims et al. 2016) or avoid them in Fourier space (Parsons et al. 2012a, b; Pober et al. 2013; Chapman et al. 2016). Since frequency maps directly to redshift and distance, this separation occurs in Fourier space parallel to the line of sight (k∥). Smooth-spectrum power will cluster around k∥ = 0, while EoR signal power will spread farther. Foreground mitigation strategies that rely on the spectral smoothness of foreground sources must also account for instrumental sources of spectral structure. One such effect is due to the natural chromaticity of an interferometer caused by the effective change of baseline length with frequency. A baseline samples angular scales on a scale inversely proportional to the baseline length in wavelengths, which means that a baseline will sample finer scales at higher frequencies. This is equivalent to the baseline drifting through Fourier space with frequency, which couples k⊥ and k∥, spreading power to higher k∥ modes (Datta, Bowman & Carilli 2010; Morales et al. 2012, 2019; Parsons et al. 2012b; Trott, Wayth & Tingay 2012; Vedantham, Udaya Shankar & Subrahmanyan 2012; Thyagarajan et al. 2013; Liu, Parsons & Trott 2014; Pober et al. 2014). Though baseline chromaticity itself is fundamental to the measurement, its effects are confined to a well-defined region of Fourier space called the foreground wedge, with a maximum k∥ set by the angular extent of the field of view, called the ‘horizon limit’. This leaves a region called the EoR window that should be uncontaminated by foregrounds. None the less, other sources of spectral structure, including the intrinsic spectra of foreground sources, the finite bandpass of the instrument, and the chromaticity of the antenna gains can spread foreground power beyond the wedge into the window. This suprahorizon power is usually accounted for in foreground avoidance approaches by avoiding k∥ modes within some buffer distance of the horizon line. Estimates of an appropriate buffer distance have been largely motivated by simple, conservative models and observations in data. Parsons et al. (2012b) argued that foreground power spillover in delay spectra is caused by the convolution of the spectral responses of the beam and sky brightness, and so is limited by the inverse bandwidth of the instrument. Through simplified foreground simulations, they demonstrate that this spillover does not exceed k∥ ∼ 0.13 Mpc−1. Thyagarajan et al. (2013) examined the effect of spectral windowing on this spillover using numerical formulations of point source and noise power spectra, finding that increased foreground suppression came at the expense of narrowing the ‘effective’ bandwidth, and hence increasing the extent to which foreground power spills into the EoR window. Pober et al. (2013) observed suprahorizon spillover in power spectra taken from PAPER data, finding that the spillover is not constant with k⊥ (equivalently, baseline length), and is larger on short baselines. Imaging the data after applying a high-pass filter to remove the wedge, they found that remaining contamination is dominated by unresolved diffuse emission. More recently, Thyagarajan et al. (2016) used wide-field visibility simulations with delay spectrum analysis to measure the levels of foreground power in the EoR window for baselines in the 19-element deployment of HERA. These simulations used a HERA beam model made through electromagnetic simulation of a HERA dish and receiver element, which incorporated realistic spectral structure due to reflections among the various components of the antenna. With these simulations, the authors determined the required performance of reflections in order for foreground power to be supressed below the EoR amplitude, arguing that HERA should be capable of detecting the EoR power spectrum for k∥ ≥ 0.2 h Mpc−1. Though these earlier works predicted that foreground power should extend as far as Δk∥ ∼ 0.2 h Mpc−1, published power spectrum limits usually choose narrower buffers when binning measurements in k. Dillon et al. (2015) used a moderate buffer of Δk∥ ∼ 0.02 h Mpc−1 when binning MWA power spectra in 1D, which noticeably left in some suprahorizon emission. Ali et al. (2015) applied a buffer of 15 ns when binning power spectra for PAPER-64 data, corresponding with Δk∥ ∼ 0.008 h Mpc−1 at their redshift of z ∼ 8.74. This buffer was set by the inverse of the bandwidth used for the analysis, and did not take into account any expectations of spectral structure in the instrument or sky. Both Beardsley et al. (2016) and the more recent Barry et al. (2019) account for suprahorizon emission by increasing the slope of the horizon line by 14 per cent, the horizon line being the line in (k⊥, k∥) space corresponding with the horizon limit. This choice makes for a wider buffer at large k⊥, which runs contrary to the expectation that foreground spillover is larger for small k⊥. Avoiding foreground-contaminated modes is essential to making a robust EoR power spectrum measurement, so it is advisable to set a wide buffer beyond the horizon limit. However, setting too wide a buffer means potentially throwing out the most sensitive measurements. It is therefore essential to have precise predictions of how foreground power spreads into the EoR window for different instruments in order to use foreground avoidance techniques. In this paper, we simulate visibility measurements from a zenith-pointing array with realistic, frequency-dependent antenna beams, observing a variety of diffuse foreground models. From these simulated observations, we make power spectra using delay-transform methods and measure how far beyond the horizon foreground power spills over as a result of the effects included in simulations. Such effects include the spectral structure of the primary beam or foreground model, baseline length, primary beam shape, and proximity of the brightest emission to the beam centre. We focus on diffuse foreground models because EoR power spectrum estimations typically only use the shortest baselines in the array, which are most sensitive to both the expected EoR signal and diffuse foregrounds. The following section gives some background on the 21 cm power spectrum, the delay-transform estimation approach, and the sources of foreground power spillover. We also discuss the theoretical explanation for the foreground pitchfork, and the choice of simulations to explore its impact on contamination. Section 3 describes the numerical simulator and the instrument and sky model configurations tested. Sections 4.1–4.4 present measurements of the effects of baseline foreshortening near the horizon, spectral windowing, source and beam spectral structure, and source position and baseline length. Lastly, we show the same results with a more realistic primary beam model in Section 4.5. 2 DIFFUSE FOREGROUND POWER The 21 cm power spectrum P21 is defined as $$\begin{eqnarray} \big\langle \widetilde{T}_{\mathrm{ b}} (\boldsymbol{k}) \widetilde{T}^*_{\mathrm{ b}} (\boldsymbol{k}^{\prime }) \big\rangle = (2 \pi)^3 \delta _\mathrm{ D}(\boldsymbol{k} - \boldsymbol{k}^{\prime }) P_{21}(\boldsymbol{k}) , \end{eqnarray}$$(1) where Tb is the brightness temperature contrast of the 21 cm line against the background radiation field, tildes denote the spatial Fourier transform, and δD is the three-dimensional Dirac delta function. The wave vector |$\boldsymbol{k}$| may be broken down into |$\boldsymbol{k}_\perp = (k_x, k_y)$|⁠, which is perpendicular to the line of sight, and k∥, which is parallel to the line of sight. We can form an estimator of the power spectrum from a measured volume |$\mathbb {V}$|⁠: $$\begin{eqnarray} \hat{P}(\boldsymbol{k}) \equiv \frac{\langle |\widetilde{T}(\boldsymbol{k})|^2\rangle }{\mathbb {V}} . \end{eqnarray}$$(2) 2.1 Delay transform There are a variety of methods available for estimating |$\widetilde{T}(\boldsymbol{k})$| from interferometric measurements, but they can be very generally categorized as either based on the delay transform on sky-reconstruction (see Morales et al. 2019 for a summary). In both categories, k∥ modes are accessed by Fourier-transforming over an axis probed by frequency. We will focus here on the delay spectrum approach (Parsons et al. 2012b, a, 2014). The complex visibility measured by a single baseline is given by the radio interferometer measurement equation (RIME) (Thompson, Moran & Swenson 2017) $$\begin{eqnarray} V(\boldsymbol{b}, \nu) = \int \nolimits _\text{sky} A(\hat{s},\nu) T_b(\hat{s},\nu) \mathrm{ e}^{-2 \pi i ((\hat{s} - \hat{s}_p) \cdot \boldsymbol{b}) \nu /c}\, \mathrm{ d}^2 s, \end{eqnarray}$$(3) where ν is the frequency, c is the speed of light, |$\boldsymbol{b}$| is the baseline vector, and A is the primary beam function giving directional gains of the measuring antennas. |$\hat{s}_p$| is a unit vector pointing to the phase centre, and the integral is carried out over all unit vectors |$\hat{s}$| over half of the unit sphere, denoted by ‘sky’. The phase centre term cancels out of the eventual estimator, so we omit it going forward. The delay transform is an inverse Fourier transform applied along the frequency axis: $$\begin{eqnarray} \widetilde{V}(\boldsymbol{b}, \tau) = \int \nolimits _{-\infty }^\infty \int \nolimits _\text{sky} \phi (\nu) A(\hat{s},\nu) T(\hat{s},\nu) \mathrm{ e}^{-2 \pi i (\tau _g - \tau) \nu } \,\mathrm{ d}^2s\:\, \mathrm{ d}\nu . \end{eqnarray}$$(4) The function ϕ(ν) is a spectral window, which both enforces the finiteness of the bandpass and allows us to applying weighting to each frequency. Here |$\tau _g = \hat{s} \cdot \boldsymbol{b} \nu /c$| is the geometric delay, corresponding to the time it takes a wavefront propagating from position |$\hat{s}$| to cross the baseline. Under a flat-sky approximation, the integral over the sky is equivalent to a 2D Fourier transform from |$\hat{s}$| to |$\boldsymbol{u}$|⁠, sampled at the baseline position |$\boldsymbol{u} = \boldsymbol{b} \nu /c$|⁠. The delay τ, dual to frequency, is then an approximate probe of line of sight Fourier modes.1 The delay-transformed visibility is thus related to the Fourier-transformed temperature |$\widetilde{T}_b$|⁠, convolved with the beam and taper function. With appropriate cosmological scaling, we can form a power spectrum estimator. $$\begin{eqnarray} \hat{P}(\boldsymbol{k}_\perp , k_\parallel) \equiv \frac{X^2 Y}{B_{pp} \Omega _{pp}} |\widetilde{V}(\tau)|^2 . \end{eqnarray}$$(5) The multiplicative factors in front approximate the inverse volume factor of equation (2). X and Y are redshift-dependent cosmological scaling factors, with units of length per angle and length per frequency, respectively. Bpp is an effective bandwidth, given by Bpp = ∫|ϕ(ν)|2 dν, which relates to the total line of sight extent observed. The primary beam squared integral |$\Omega _{pp} = \int |A(\hat{s})|^2 d^2 s$| gives the angular area. The exact derivation of this factor is given by appendix B of Parsons et al. (2014).2 The cosmological factors are $$\begin{eqnarray} X(z) = \chi (z) \qquad Y(z) = \frac{c (1+z)^2}{H(z) \nu _{21}}, \end{eqnarray}$$(6) where z is redshift, χ is the comoving distance, H(z) is the Hubble parameter, ν21 is the rest-frame 21 cm frequency, and c is the speed of light. The |$\boldsymbol{k}_\perp$| and k∥ modes corresponding to a given baseline |$\boldsymbol{u}$| and delay τ are given by $$\begin{eqnarray} \boldsymbol{k}_\perp = \frac{2\pi \boldsymbol{u}}{X} \qquad k_\parallel = \frac{2 \pi \tau }{Y} . \end{eqnarray}$$(7) For the above relations, X and Y are evaluated at the centre of the bandpass and assumed not to evolve much over the chosen bandpass. As mentioned before, the correspondence between k∥ and τ is only approximate, and gets worse for longer baselines. Throughout this paper, we will use a ΛCDM cosmological model with parameters measured by Planck Collaboration XIII (2016). 2.2 Foreground power spillover In the case of a single point source at position |$\hat{s}_p$|⁠, with corresponding geometric delay τp, and spectrum Ip(ν), equation (4) reduces to $$\begin{eqnarray} \widetilde{V}(\boldsymbol{b}, \tau) = (\widetilde{\phi } \ast \widetilde{A} \ast \widetilde{I}_p)(\tau - \tau _p) , \end{eqnarray}$$(8) where * denotes convolution in delay and tildes indicate a delay-transformed function. If the beam, source spectra, and window function are all flat in frequency, and the bandpass is infinite, then this becomes a delta function centred at the geometric delay of the point source. The maximum geometric delay is limited to the baseline length |$b = |\boldsymbol{b}|$| over the speed of light. $$\begin{eqnarray} |\tau _\mathrm{ g}| \le \frac{b}{c} . \end{eqnarray}$$(9) Under these conditions, the power from a point source on the horizon has the form of a delta function centred at this maximum delay. Hence, this maximum delay is known as the ‘horizon line’, and separates flat-spectrum power from spectrally structured emission. In cosmological terms, this horizon limit is given by $$\begin{eqnarray} k_\parallel \le \frac{H(z) \chi (z)}{c (1+z)} k_\perp , \end{eqnarray}$$(10) where |$k_\perp = |\boldsymbol{k}_\perp |$|⁠. This limit applies in the ideal point source case as discussed − infinite bandwidth and a flat spectrum. This power will spread out, possibly beyond the horizon limit, once these assumptions are relaxed. For example, assume the primary beam A and point source flux Ip have no spectral structure, but the window function has a limited domain. In this case, |$\widetilde{A}(\tau) = A_p \delta _\mathrm{ D}(\tau)$| and |$\widetilde{I}_p(\tau) = I_p \delta _\mathrm{ D}(\tau)$|⁠, where δD is the Dirac delta function, and equation (8) reduces to $$\begin{eqnarray} \widetilde{V}(\boldsymbol{b}, \tau) = I_p A_p \widetilde{\phi }(\tau - \tau _p) . \end{eqnarray}$$(11) If ϕ is a rectangular window function, then |$\widetilde{\phi }$| is a sinc function centred on τ = τp, which has sidelobes that are not confined to the horizon. The sinc function is known to lack the dynamic range needed to access the expected EoR amplitude, so most analyses use other common window functions such as Hamming, Blackman–Harris, or Blackman–Nuttall. These windows suppress their sidelobe amplitude at the expense of widening the main lobe. In any case, the finiteness of the bandpass necessarily causes even flat-spectrum power to spread. 2.3 Foreground pitchfork Resolved diffuse emission behaves in a fundamentally different way from unresolved point-like sources in interferometric measurements. The difference can be illustrated in angular Fourier space by taking the flat sky approximation of equation (3) such that the exponential fringe term acts as a 2D Fourier transform. Applying the convolution theorem, we can see that the visibility measured by a baseline is given by the convolution of the primary beam and sky brightness in Fourier space: $$\begin{eqnarray} V(u,v) \simeq \int A(l,m) T(l,m) e^{-2 \pi i (ul + vm)}\, \mathrm{ d}l \mathrm{ d}m \end{eqnarray}$$(12) $$\begin{eqnarray} \hphantom{V(u,v)}= (\widetilde{A} \ast \widetilde{T}) (u,v) \end{eqnarray}$$(13) Here, (l, m) are direction cosines, tildes denote the 2D Fourier transform and (u, v) are orthogonal components of the baseline vector measured in wavelengths. These are also the Fourier-duals to l and m. A point source has a surface brightness I in the form of a delta function, which is constant in Fourier space. A diffuse source, however, will be more compact in Fourier space, centred on (u, v ≈ 0. Short baselines will therefore measure more diffuse power than long baselines, while the measured power from point sources will be the same for both. The increased brightness of diffuse power on short baselines leads to an increase in measured brightness when diffuse sources are near the horizon, due to the shorter projected baseline length. This puts excess power near the horizon in delay space, leading to a characteristic pitchfork pattern, observed in simulations and in MWA data (Thyagarajan et al. 2015a, b). Since the projected baseline length can go to effectively zero at the horizon, even a monopole signal, normally undetectable to interferometers, produces a signal. The pitchfork effect can potentially change our intuition for how foreground power spreads beyond the horizon. For a zenith-tracking array, measured power is highest when bright sources are near zenith, because the primary antenna beam is strongest there. In delay space, however, this power is clustered around τ = 0, so is farther from the horizon limit . Power in the pitchfork is very close to the horizon, and so is more likely to spill over into the EoR window, but is attenuated by a weaker primary beam response. This raises a question – Does a bright patch of the sky contaminate the window more when it is at zenith or near the horizon? To answer this, we will look at how much power spreads beyond the horizon as a function of the hour angle of the Galactic centre, for a variety of primary beam models, in simulated data. Some of these beam models have effectively zero response near the horizon, which prevents any pitchfork from appearing. Comparing these, we can demonstrate that the pitchfork is the dominant source of power near the horizon for most observing time, while spectral structure in the beam or source affects how that power is spread. By taking both effects into consideration, it is possible to set bounds on the extent of foreground contamination. 3 SIMULATIONS Simulations were carried out by numerically integrating the RIME. The integration is done from horizon to horizon, taking into account the curvature of the sky. The simulator healvis is publicly available,3 previously used by Lanman & Pober (2019) and indexed by the ASCL (Lanman & Kern 2019). The brightness temperature |$T(\hat{s})$| describes the specific intensity as a function of angle. We can discretize the sky by dividing it into pixels p, each of which is defined by a function |$W_p(\hat{s})$| which is 1 when |$\hat{s}$| points into the pixel and 0 otherwise. The discretized form of the RIME (equation 3) has the form: $$\begin{eqnarray} V(\boldsymbol{b}, \nu) = \sum \limits _p \int \nolimits _\text{sky} W_p(\hat{s}) T(\hat{s}, \nu) A(\hat{s}) \mathrm{ e}^{-2\pi i \hat{s} \cdot \boldsymbol{b} \nu /c} \,\mathrm{ d}\Omega , \end{eqnarray}$$(14) $$\begin{eqnarray} \hphantom{V(\boldsymbol{b}, \nu)} \approx \sum \limits _p T(\hat{s}_p, \nu) A(\hat{s}_p)\, \mathrm{ e}^{-2\pi i \hat{s}_p \cdot \boldsymbol{b} \nu /c} \int \nolimits _\text{sky} W_p(\hat{s})\, \mathrm{ d}\Omega . \end{eqnarray}$$(15) The second step assumes that the fringe, surface brightness, and primary beam are all constant across the pixel, and so they can be moved outside the integral. The integral over each pixel window function then reduces to the solid angle of that pixel, ωp. $$\begin{eqnarray} V_b(\nu) = \sum \limits _{p} \omega _p T(\hat{s}_p, \nu) A(\hat{s}_p) \,\mathrm{ e}\mathrm{ }^{-2\pi i \hat{s}_p \cdot \boldsymbol{b} \nu /c} \end{eqnarray}$$(16) Equation (16) treats the diffuse sky model as a set of point sources at pixel centres with specific flux densities Tpωp. We thus refer to this as the point source approximation. As discussed in Section 2.3, the difference in behaviour between point sources and diffuse models is most significant on long baselines, so we need to ensure the map resolution is fine enough for this approximation to be valid. Appendix A discusses the limitations of the point source approximation and tests its validity for the healvis simulations discussed here. Precise and verified simulation is important for 21 cm cosmology, as the simulation technique itself can impart unexpected spectral structure on the simulated data, or fail to model an significant process. healvis has been verified against the reference simulator pyuvsim (Lanman et al. 2019), and is tested with analytic solutions for diffuse models (to be discussed in a forthcoming paper). Table 1 summarizes the various simulation configurations used, including antenna layouts, primary beam models, sky brightness distribution models, observing time ranges and time-step sizes. The same bandpass and channelization, covering 100–150 MHz at 97 kHz channel resolution, is used for all simulations. Generally, the ‘off-zenith’ time configuration is used whenever LST is not considered important to the result. This time set consists of only five 11 s integrations, from a time when the Galactic centre is in the sky but not near the beam centre (off by 4 h). This is considered as a typical sky brightness, an assumption supported by the tests with time dependence (using the ‘transit’ time set). More details on the sky and beam models are given in the subsequent subsections. Table 1. A summary of parameters used in simulations. Beam models Airy disc Radiation pattern of a circular aperture with a 14 m diameter (equation 17). Achromatic Gaussian |$\exp [-\theta ^2/(2 \sigma _0^2)]$|⁠, with σ0 = 7.37○. This is set to match the FWHM of the main lobe of the Airy beam (equation 18). Chromatic Gaussian The same as achromatic, but with σ = σ0(ν/ν0)α (equation 20). HERA CST Power beams made from far-field electric fields measured in CST simulations of a HERA dish with crossed-dipole feeds. Sky models GSM The global sky model in healpix, generated from three principal components per pixel (Zheng et al. 2017). Sym A model with the same angular power spectrum as the GSM, defined to be azimuthally symmetric about the Galactic centre. Each pixel is given a power-law frequency spectrum. Layout 37 hex A 37-element hexagon, with 14.6 m spacing. Imaging 80 antennas distributed randomly, with a Gaussian radial density profile with FWHM of 230 m. The longest baseline is 608 m. Line A line of 300 antennas, evenly spaced, east to west. Baselines are formed by pairing each of these antennas with a 301st antenna, such that the whole array provides east/west baselines with lengths ranging from 10 to 150 m. Time Off-zenith Five 11 s time-steps when the Galactic centre is 4 h off of zenith. Transit 24 h at 5 min time-step, such that the Galactic centre transits 6 h from the start. Bandpass 50 MHz 100–150 MHz with channels of width 97.65 kHz. The channel width is chosen to match that of the first generations of HERA. Beam models Airy disc Radiation pattern of a circular aperture with a 14 m diameter (equation 17). Achromatic Gaussian |$\exp [-\theta ^2/(2 \sigma _0^2)]$|⁠, with σ0 = 7.37○. This is set to match the FWHM of the main lobe of the Airy beam (equation 18). Chromatic Gaussian The same as achromatic, but with σ = σ0(ν/ν0)α (equation 20). HERA CST Power beams made from far-field electric fields measured in CST simulations of a HERA dish with crossed-dipole feeds. Sky models GSM The global sky model in healpix, generated from three principal components per pixel (Zheng et al. 2017). Sym A model with the same angular power spectrum as the GSM, defined to be azimuthally symmetric about the Galactic centre. Each pixel is given a power-law frequency spectrum. Layout 37 hex A 37-element hexagon, with 14.6 m spacing. Imaging 80 antennas distributed randomly, with a Gaussian radial density profile with FWHM of 230 m. The longest baseline is 608 m. Line A line of 300 antennas, evenly spaced, east to west. Baselines are formed by pairing each of these antennas with a 301st antenna, such that the whole array provides east/west baselines with lengths ranging from 10 to 150 m. Time Off-zenith Five 11 s time-steps when the Galactic centre is 4 h off of zenith. Transit 24 h at 5 min time-step, such that the Galactic centre transits 6 h from the start. Bandpass 50 MHz 100–150 MHz with channels of width 97.65 kHz. The channel width is chosen to match that of the first generations of HERA. Open in new tab Table 1. A summary of parameters used in simulations. Beam models Airy disc Radiation pattern of a circular aperture with a 14 m diameter (equation 17). Achromatic Gaussian |$\exp [-\theta ^2/(2 \sigma _0^2)]$|⁠, with σ0 = 7.37○. This is set to match the FWHM of the main lobe of the Airy beam (equation 18). Chromatic Gaussian The same as achromatic, but with σ = σ0(ν/ν0)α (equation 20). HERA CST Power beams made from far-field electric fields measured in CST simulations of a HERA dish with crossed-dipole feeds. Sky models GSM The global sky model in healpix, generated from three principal components per pixel (Zheng et al. 2017). Sym A model with the same angular power spectrum as the GSM, defined to be azimuthally symmetric about the Galactic centre. Each pixel is given a power-law frequency spectrum. Layout 37 hex A 37-element hexagon, with 14.6 m spacing. Imaging 80 antennas distributed randomly, with a Gaussian radial density profile with FWHM of 230 m. The longest baseline is 608 m. Line A line of 300 antennas, evenly spaced, east to west. Baselines are formed by pairing each of these antennas with a 301st antenna, such that the whole array provides east/west baselines with lengths ranging from 10 to 150 m. Time Off-zenith Five 11 s time-steps when the Galactic centre is 4 h off of zenith. Transit 24 h at 5 min time-step, such that the Galactic centre transits 6 h from the start. Bandpass 50 MHz 100–150 MHz with channels of width 97.65 kHz. The channel width is chosen to match that of the first generations of HERA. Beam models Airy disc Radiation pattern of a circular aperture with a 14 m diameter (equation 17). Achromatic Gaussian |$\exp [-\theta ^2/(2 \sigma _0^2)]$|⁠, with σ0 = 7.37○. This is set to match the FWHM of the main lobe of the Airy beam (equation 18). Chromatic Gaussian The same as achromatic, but with σ = σ0(ν/ν0)α (equation 20). HERA CST Power beams made from far-field electric fields measured in CST simulations of a HERA dish with crossed-dipole feeds. Sky models GSM The global sky model in healpix, generated from three principal components per pixel (Zheng et al. 2017). Sym A model with the same angular power spectrum as the GSM, defined to be azimuthally symmetric about the Galactic centre. Each pixel is given a power-law frequency spectrum. Layout 37 hex A 37-element hexagon, with 14.6 m spacing. Imaging 80 antennas distributed randomly, with a Gaussian radial density profile with FWHM of 230 m. The longest baseline is 608 m. Line A line of 300 antennas, evenly spaced, east to west. Baselines are formed by pairing each of these antennas with a 301st antenna, such that the whole array provides east/west baselines with lengths ranging from 10 to 150 m. Time Off-zenith Five 11 s time-steps when the Galactic centre is 4 h off of zenith. Transit 24 h at 5 min time-step, such that the Galactic centre transits 6 h from the start. Bandpass 50 MHz 100–150 MHz with channels of width 97.65 kHz. The channel width is chosen to match that of the first generations of HERA. Open in new tab 3.1 Primary beam models The primary beam |$A(\hat{s})$| gives the direction-dependent gains of the antenna receiver elements, normalized to the value at the antenna’s pointing centre. The primary beam of a given antenna can be found by considering the antenna as a transmitter, fed by a uniform voltage, and measuring the far-field electric field pattern, and normalizing to the value at zenith. The power beam is obtained by taking the squared amplitude of this electric field. We use three analytically defined beam models and one numerically calculated through electromagnetic simulations of a HERA dish. The analytic models are described in equations (17) to (20). The Airy beam is defined as the far-field radiation pattern of a uniformly illuminated circular aperture: $$\begin{eqnarray} A_\text{airy}(\theta , \nu) = \frac{\nu /c}{\pi D \sin {\theta }}\: J_1 \left(\frac{\pi D \sin {\theta }}{\nu /c} \right) , \end{eqnarray}$$(17) where θ is the angle from the pointing centre, J1 is the first-order Bessel function of the first kind, c is the speed of light, and D = 14 m is the antenna aperture. The Airy beam has a simple analytic form that behaves similarly to a physical antenna, in that it has a non-negligible response near the horizon and has nulls that move with frequency. The achromatic Gaussian beam is defined by a Gaussian function with a constant full width at half-maximum (FWHM) across all frequencies, defined by a width parameter σ0: $$\begin{eqnarray} A_\text{AG} (\theta) = \exp \left(-\frac{\theta ^2}{2\sigma _0^2}\right) \end{eqnarray}$$(18) This beam has no sidelobes and is effectively zero at the horizon. We choose σ0 = 7.37○, which sets the FWHM equal to the width of the main lobe of the Airy disc beam. The chromatic Gaussian beam is defined just as equation (18), but with a width that evolves according to a power law. $$\begin{eqnarray} A_\text{CG}(\theta , \nu) = \exp \left(-\frac{\theta ^2}{2\sigma (\nu)^2}\right) \end{eqnarray}$$(19) $$\begin{eqnarray} \sigma (\nu) = \sigma _0 \left(\frac{\nu }{\nu _0}\right)^\alpha \end{eqnarray}$$(20) The reference frequency ν0 is chosen to be the first frequency channel centre. Unless otherwise noted, the spectral index s is chosen to be −1, which resembles the frequency evolution of the Airy beam. The last beam model is defined numerically by a set of CST Suite4 simulations of a realistic HERA dish and receiver done by Fagnoni & Acedo (2016). The simulated antenna included a parabolic reflecting dish with a 14 m diameter, two 1.3 m copper dipole receivers with a pair of aluminium discs acting as sleeves, and a cylindrical mesh ‘skirt’ 36 cm high around the feed to reduce coupling with other antennas in the array (see Fagnoni & Acedo 2016 for more details). Reflections among the various components of the antenna introduce spectral structure into the beam, which is largely captured by the CST simulations. These simulated beam models are stored as a set of electric field components at altitude and azimuth positions on a spherical grid with steps of 1○ in both latitude and longitude, at frequencies between 50 and 250 MHz in steps of 1 MHz. For use in simulations of unpolarized (i.e. purely stokes-I) sky models, we convert this electric field beam into a power beam in the following way. The electric field beam of an antenna feed p may be expressed in a basis of unit vectors |$\hat{\theta }$| and |$\hat{\varphi }$| along the zenith angle and azimuthal directions, respectively: $$\begin{eqnarray} \boldsymbol{A}^p(\hat{s},\nu) = A^p_\theta (\hat{s}, \nu) \hat{\theta } + A^p_\varphi (\hat{s}, \nu) \hat{\varphi }. \end{eqnarray}$$(21) The components form the elements of a Jones matrix, describing the directional response of the antenna feed to the electric field components of the sky. This is described in detail in Kohn et al. (2019). For unpolarized sky models, the response of each feed is given by the power of the electric field: $$\begin{eqnarray} A^{pp}(\hat{s},\nu) = \big|A^p_\theta (\hat{s}, \nu)\big|^2 + \big|A^p_\varphi (\hat{s}, \nu)\big|^2. \end{eqnarray}$$(22) The CST simulated beams have two crossed-dipole feeds, one oriented east/west (designated X) and the other oriented north/south (Y). Under the same formalism, we may think of the analytic beams defined in this section as representing the response of an idealized single feed to an unpolarized sky model. 3.1.1 CST beam interpolation Our simulations are done at much higher spectral and angular resolution than the CST simulations, so we need to interpolate the CST simulated beam models values to the required angles and frequencies. The beam data is handled using the new pyuvdata class UVBeam (Hazelton et al. 2017).5 This class supports frequency interpolation at each pixel using the methods accessible to the scipy.interpolate.interp1d method, which includes linear, polynomial splines, and nearest-neighbour interpolation among other methods. UVBeam handles angular interpolation using a cubic bivariate spline. By interpolating the beam to higher frequency resolution, we are extrapolating the instrument’s response to higher delay modes. Since the main results of this paper concern the behaviour of foregrounds at high delays, we need to be sure that our interpolation does not introduce undue power. In this section, we interpolate the UVBeam to a finer frequency resolution using a variety of interpolation methods and compare these interpolated beams to the original beam data under a delay transform. The power beam is obtained from the E-field using equation (22), then the values at each frequency are normalized to their peaks. The peak-normalized power beam is then, at each pixel, interpolated to the 50 MHz bandpass (see Table 1) used for simulations. There are two general types of interpolation used: polynomial spline interpolation, which we do to linear, third, and fifth order, and Gaussian process regression (GPR) using the scikit-learn package. The GPR method uses a radial basis function (RBF) kernel, also known as a squared exponential, as the covariance function of the Gaussian process. For a fixed point on the sky, we can describe the covariance of the beam at two frequencies ν1 and ν2 as $$\begin{eqnarray} k(\nu _1,\nu _2|\ell) = \exp \left[-\frac{1}{2}(\nu _1 - \nu _2)\ell ^{-2}(\nu _1 - \nu _2)\right], \end{eqnarray}$$(23) where ℓ is a hyperparameter that sets the characteristic scale of the covariance function. Given a covariance function, the GPR maximum a posteriori estimate of the underlying signal constrained by the data can be calculated at any new frequency ν′ (equation 2.23 of Rasmussen & Williams 2006). The RBF kernel has the property of producing smoothly varying estimates with frequency structure set by the length-scale hyperparameter. By keeping this length-scale fixed (rather than solving for its optimal value given the data), the RBF kernel can act as a low-pass filter, filtering out structure at scales above the length-scale (e.g. Kern et al. 2019). Given the inherent Nyquist sampling rate of the CST output, we set ℓ = 1 MHz for our GPR interpolation method. Altogether, this leaves us with five methods of interpolation to test. To test the effect of interpolation on the beam, we construct a single metric to compare in delay space. At each frequency, the interpolated beam values are summed over pixels. This list of 512 sums is then is multiplied by a Blackman-Harris window function and then Fourier transformed from frequency to delay using an inverse DFT. The result of these steps is a quantity analogous to the power spectrum measured at k⊥ = 0. Fig. 1 shows this quantity for each of the interpolation types. The GPR results in the plot are shown with the chosen smoothing scales in parentheses. Overall, all of the interpolation schemes preserve the power on the known scales, marked with black dots, up to the Nyquist limit (shown by the vertical dashed line). Beyond that, the spline interpolations introduce noticeable structure. Further testing has shown that this is likely due to some sharp features in the spectra of pixels far off from the beam centre, which can cause aliasing under Fourier transform. The GPR-interpolation seems to smooth out these effects well, and so the curves for the GPR methods are relatively smooth beyond the Nyquist limit. Figure 1. Open in new tabDownload slide Comparison of the delay-transformed beam integral of the (black dots) non-interpolated beam to the interpolated beam with five methods of interpolation. The vertical dashed line marks 500 ns, the Nyquist sampling rate of the 1 MHz resolution beam. The inset zooms in on the region around the Nyquist limit. There is some disagreement among the curves at and just before this limit, so it seems the effects of interpolation are not strictly limited to interpolated delay modes. Figure 1. Open in new tabDownload slide Comparison of the delay-transformed beam integral of the (black dots) non-interpolated beam to the interpolated beam with five methods of interpolation. The vertical dashed line marks 500 ns, the Nyquist sampling rate of the 1 MHz resolution beam. The inset zooms in on the region around the Nyquist limit. There is some disagreement among the curves at and just before this limit, so it seems the effects of interpolation are not strictly limited to interpolated delay modes. To avoid any potential issues that may depend on the beam interpolation scheme, we avoid using delay modes beyond the Nyquist limit of 500 ns for any results. Simulations are still carried out to the full frequency resolution of HERA, with beams interpolated using the GPR method with a 1 MHz smoothing scale. This is a conservative choice, because the smoothing scale avoids accidentally smoothing over sub-Nyquist scales in the beam while suppressing the interpolation artefacts. It is of course possible that the HERA beam has structure on scales smaller than 1 MHz, but that cannot be determined from the existing CST simulation data. 3.2 Sky models Our simulations model the Stokes-I surface brightness temperature of the sky (in Kelvin) as a set of healpix6 maps (Gorski et al. 2005), one for each frequency channel. healpix divides the sky into equal-area pixels with a resolution set by a single parameter, Nside, such that the individual pixel area is ω = 4π/(12 × Nside2) sr. With the healvis simulator, there are trade-offs between resolution, run-time, and memory requirements, which are exacerbated by simulating wide fields of view and bandpasses. For this work we use maps at an Nside of 256, which corresponds with a resolution of 13.7 arcmin. The highest resolution probed by any of our simulation configurations is 14.13 arcmin, for the longest baseline and highest frequency, and shorter baselines are only sensitive to larger scales. Appendix A shows the effect of increasing resolution on the measured power spectrum of a GSM model is not noticeably improved by increasing the resolution above Nside 256. We use two diffuse foreground models. The first is the 2016 global sky model (GSM), which is compiled from 29 different sky maps covering a total bandpass of 10 MHz to 5TH (Zheng et al. 2017). The GSM is implemented in the Python package pygsm,7 which generates healpix GSM maps at an Nside of 1024 using a principal component decomposition to model the spectral structure at each pixel. We degrade this to an Nside of 256 by averaging together high-resolution pixels nested within low-resolution pixels, which can be done without interpolating or splitting pixel values due to the hierarchical structure of healpix maps. In our chosen bandpass, diffuse foreground power is dominated by Galactic synchrotron emission, which generally follows a power law with a spectral index of −2.5. Most pixels in the GSM model used here have power law spectra with index ∼−2.5. The second model is the ‘symmetric GSM’ or ‘sym’, which is constructed from the GSM but made to have azimuthal symmetry about the Galactic centre position. The sym model thus provides a well-defined locus of bright emission on the sky, which will make for a more meaningful test of how spillover relates to the ‘position of bright emission’. The sym model is constructed in the following way: Calculate the angular power spectrum $\mathcal {C}_l$ of the GSM at 100 MHz from the Nside 256 GSM. Generate a set of spherical harmonic expansion coefficients alm, with $a_{l0} = (2l + 1)\sqrt{\mathcal {C}_l}$ and alm = 0 for m ≠ 0. Apply a rotation to the alm using Wigner D-matrices so that the (θ, ϕ) = (0, 0) point is at the Galactic centre. Convert the alm to an Nside 256 map. Scale each pixel in this map by a power law in frequency: (ν/ν0)α for ν0 =100 MHz and spectral index α. For most tests, we use an achromatic symmetric model, with α = 0. The angular power spectrum calculation and spherical harmonic transforms are done using healpy’s anafast and alm2map functions, respectively. The rotation is done using rotate_map_alms. The GSM and symmetric GSM models at 100 MHz are shown in ICRS equatorial coordinates in Fig. 2, with a log scale on each plot. Figure 2. Open in new tabDownload slide The GSM (a) and symmetric GSM (b) models in a Mollweide projection of equatorial coordinates. Figure 2. Open in new tabDownload slide The GSM (a) and symmetric GSM (b) models in a Mollweide projection of equatorial coordinates. 3.3 Power spectrum estimation and measurements The form of the delay spectrum estimator was discussed in Section 2, in particular with equation (5). In this section, we will briefly discuss the steps of estimating power spectra from simulated visibilities, and define the power spillover and fiducial EoR amplitude that will be used to quantify the spread of power beyond the horizon. Rather than look at the power at specific k∥ modes, we are interested in the range of k∥ that are contaminated by foreground power. The simulator provides visibilities V[b, t, ν] for each baseline (b), time-step (t), and frequency (ν), in units of Jy, as well as the beam integral Ωpp at the zeroth frequency in units of steradians. The visibilities are converted from Jy to mK sr by weighting each channel with a frequency-dependent conversion factor, $$\begin{eqnarray} V_\text{mK sr}[b,t,\nu ] = 10^{-20} \left(\frac{c^2}{2 k_B \nu ^2} \right) V_\text{Jy}[b,t,\nu ] , \end{eqnarray}$$(24) where the factor of |$10^{-20} = 10^{-23}\, 10^{3}$| comes from the conversion from Jy to erg s−1 cm−2 Hz−1 combined with the conversion of K to mK. For each baseline and time, the delay transform is done by weighting each channel with the spectral window function and then doing an inverse discrete Fourier transform (IDFT), and taking the absolute square: $$\begin{eqnarray} \widetilde{V}[b, t, \tau _q] = \frac{B}{N_f}\sum \limits _p \phi (\nu _p) V_\text{mK} \mathrm{ e}^{2\pi \mathrm{ i} \nu _p \tau _q}, \end{eqnarray}$$(25) $$\begin{eqnarray} P[b, t, k_\parallel ] = \frac{X^2 Y}{B_{pp} \Omega _{pp}}|\widetilde{V}[b, t, \pm \tau _q] |^2 . \end{eqnarray}$$(26) In the above, νp is the pth frequency, and τq = q/(2B) for bandwidth B are the delay modes. Note that there are two delay modes for every k∥, since the IDFT goes to both positive and negative delays. Since the power spectrum is a real quantity, the positive and negative delay modes are equivalent, so we average them together. This leaves us with an estimator of the power spectrum for each baseline, time, and k∥. We average together power spectra from baselines with similar lengths (to within a half meter). Additional tests show no significant variation in the power spectrum with baseline orientation. For simulations using the off-zenith time configuration, averaging power spectra incoherently over the time-steps (55 s). The spectrum is not observed to change significantly over such a short period. For comparison with the foregrounds, we consider a fiducial EoR power spectrum amplitude of P21 = 292 mK2 Mpc3, which is approximately the amplitude at k ∼ 1 Mpc−1 at the band-centre redshift z ∼ 10.4. This amplitude is derived from a 21cmFAST EoR simulation with default settings (Mesinger, Furlanetto & Cen 2011). The true EoR power spectrum is stronger at smaller k at this redshift, so this is a conservatively low choice for comparison with the foreground power spectra. Using the fiducial EoR amplitude, we can define the power spectrum spillover Δk∥ as the distance in k∥ that the foreground power extends beyond the horizon. More precisely, the spillover is defined as the largest k∥ where the power spectrum crosses below a set threshold of P21, such that P(k∥) < P21 for all k∥ > Δk, less the horizon kh: $$\begin{eqnarray} \Delta k_\parallel (P_{21}) \equiv \max \left(k_\parallel | P(k_\parallel) = P_{21}\right) - k_h . \end{eqnarray}$$(27) Defining the spillover relative to the horizon makes it easier to compare the results among different length baselines, which have different horizon limits. Fig. 3 shows power spectra for a single baseline with several different beam models, with the threshold/horizon/spillover labelled. These spectra are from a simulation with the off-zenith time configuration, and with the 14.6 m baselines binned from the hex layout. Figure 3. Open in new tabDownload slide The delay spectra of the GSM for a single 14 m baseline with a variety of beam types. The fiducial EoR amplitude is shown by the horizontal line, the black vertical dashed line shows the horizon limit, and the foreground power spillover Δk∥ for the Gaussian beams is indicated by the horizontal green arrow. A vertical dashed red line marks the Nyquist limit for the 1 MHz resolution of the CST beam. For the results with the CST beam, we will use a higher threshold to avoid measuring spillovers past the Nyquist limit. Such measurements would be highly dependent on the frequency interpolation scheme used. Figure 3. Open in new tabDownload slide The delay spectra of the GSM for a single 14 m baseline with a variety of beam types. The fiducial EoR amplitude is shown by the horizontal line, the black vertical dashed line shows the horizon limit, and the foreground power spillover Δk∥ for the Gaussian beams is indicated by the horizontal green arrow. A vertical dashed red line marks the Nyquist limit for the 1 MHz resolution of the CST beam. For the results with the CST beam, we will use a higher threshold to avoid measuring spillovers past the Nyquist limit. Such measurements would be highly dependent on the frequency interpolation scheme used. 4 RESULTS 4.1 Pitchfork confirmation The pitchfork foreground signature in Fourier space was first found by simulations in PRISim (Thyagarajan et al. 2015a), and later confirmed in data form the Murchison Widefield Array (MWA) experiment (Thyagarajan et al. 2015b). PRISim was the first interferometer simulator that could handle diffuse sky models down to the horizon, and was thus sensitive to the effects of baseline foreshortening there. healvis is also sensitive to sources near the horizon, and so we expect to see pitchfork power in our simulations. The first set of simulations used a monopole sky signal to confirm the existence of the pitchfork, since a monopole should have minimal power away from the horizon for most baseline lengths. The simulated array consists of 300 east–west oriented baselines with lengths from 10 to 150 m, in even steps. This provides a range of k⊥ modes to study. We do see the pitchfork for randomly distributed antennas of all orientations as well (see Appendix A). We ran two simulations – first using the Airy beam, equation (17), and then with the achromatic Gaussian beam, equation (18). Fig. 4 shows the resulting delay spectra. The pitchfork appears for the Airy disc beam, since that has a non-zero response near the horizon. The Gaussian beam is effectively zero at the horizon, and so no pitchfork can be seen. We also observe a strong signal on the shortest baselines. This is evidence that short baselines of an interferometer may be made sensitive to a global signal, as discussed in Presley, Liu & Parsons (2015). Figure 4. Open in new tabDownload slide Delay spectra of a monopole signal for a variety of baseline lengths (top axis), and their corresponding k⊥ mode (bottom axis). (a) shows the results with an Airy beam, which has strong sidelobes and nonzero response down to the horizon. (b) is the same for a Gaussian beam, with width set to the width of the main lobe of the Airy beam. Black lines show the delay corresponding with sources on the horizon (b/c). Note the difference in colour scales. Figure 4. Open in new tabDownload slide Delay spectra of a monopole signal for a variety of baseline lengths (top axis), and their corresponding k⊥ mode (bottom axis). (a) shows the results with an Airy beam, which has strong sidelobes and nonzero response down to the horizon. (b) is the same for a Gaussian beam, with width set to the width of the main lobe of the Airy beam. Black lines show the delay corresponding with sources on the horizon (b/c). Note the difference in colour scales. The observed pitchfork appears to extend beyond the horizon by a consistent amount for all baseline lengths, suggesting that the pitchfork represents a constant minimum level of foreground spillover for antennas with a non-zero horizon response. Most realistic antennas used in EoR experiments have some response near the horizon, since the width of the beam varies inversely with the effective diameter of the antenna, and most of these experiments use relatively small antenna elements. As we will see, the pitchfork is an important factor of foreground spillover for most of the cases studied here. 4.2 Window function comparison As discussed in Section 2.2, the shape and width of the bandpass has an effect on the spread of power beyond the horizon, as the intrinsic power spectrum of the source is convolved in delay with the spectral window function and the beam. By choosing an appropriate window function to weight the frequency channels, one can reduce the amplitude of power beyond the horizon at the expense of widening the main lobe of the foreground signal. In this section, we explore the effects of several window functions on power spectra derived from a single symmetric-GSM simulation. The simulations used the Airy beam, off-zenith time configuration, and 37 hex layout. The window function is applied to the delay transform for each baseline separately, and the resulting power spectra were averaged over time and binned in k⊥ such that spectra from baselines of the same length to within 10 cm are binned together. The results presented here are for 14.6 m baselines. The most basic window function is the rectangular window, sometimes called a Dirichlet window, which applies even weighting to all samples in the bandpass. As mentioned before, in Fourier space this has the form of a sinc function, which has a limited dynamic range. We include it here for comparison, though it is not used in most 21 cm power spectrum analysis. Note that applying a window function, other than the rectangular, involves down-weighting channels towards the end of the bandpass. This reduces the effective bandwidth Bpp, and hence increases the width of the main lobe of the Fourier transform, which increases the basic spillover beyond the horizon. There is thus an implicit trade-off between the level of suppression in the foreground sidelobes and the loss of low-k modes to foreground spillover. We test two window functions, in addition to the unweighted rectangular spectral window. The first is the Blackman–Harris window, which is defined by a sum of three cosine functions optimized to reduce the level of sidelobes in Fourier space (see e.g. Harris (1978)). The Blackman–Harris window reduce has a dynamic range of about 120 dB in power, which is generally considered sufficient for 21 cm EoR experiments. Thyagarajan et al. (2016) applied a squared Blackman-Harris window to achieve additional 10 orders of magnitude in sidelobe suppression. The second window is the Dolph–Chebyshev, which is designed to minimize the sidelobe amplitude for a choice of main lobe width and thus has a tunable dynamic range (Rabiner & Gold 1975; Smith 2019). Functionally, it is constructed from Chebyshev polynomials in Fourier space and then inverse-transformed to direct space. Unlike the Blackman–Harris and rectangular windows, the sidelobes do not drop off with k, and so usually appear as flat lines in our power spectra. We use Dolph–Chebyshev windows set to 120, 150, and 180 dB ranges, which we refer to as DC120, DC150, and DC180, respectively. Fig. 5 shows the comparison of power spectra with the different window functions. The rectangular window function, which has the form of a sinc function in delay space, does not have sufficient dynamic range to keep spilled foreground power below the expected EoR amplitude. The three Dolph–Chebyshev windows demonstrate the expected trade-off between main lobe width and sidelobe suppression, as can be seen in the slight spread among the curves around P(k) ∼ 109 mK2 Mpc3). The Blackman–Harris window displays similar behaviour to the DC120 window, and continues to drop off at higher k∥ while the Dolph–Chebyshev window function stays flat. Figure 5. Open in new tabDownload slide Power spectra with rectangular (rect), Blackman–Harris (bh) window functions, as well as Dolph–Chebyshev window functions with dynamic range set to 120 and 150 dB (DC120 and DC150, respectively). The horizontal black line marks a fiducial EoR amplitude (see Section 4.2). The vertical dashed line marks the horizon for a 14.6 m baseline. The Blackman–Harris function alone does not provide sufficient dynamic range to push the leaked foreground power below the EoR level. By contrast, applying the right Dolph–Chebyshev window can suppress the leaked foreground power below the EoR threshold. Figure 5. Open in new tabDownload slide Power spectra with rectangular (rect), Blackman–Harris (bh) window functions, as well as Dolph–Chebyshev window functions with dynamic range set to 120 and 150 dB (DC120 and DC150, respectively). The horizontal black line marks a fiducial EoR amplitude (see Section 4.2). The vertical dashed line marks the horizon for a 14.6 m baseline. The Blackman–Harris function alone does not provide sufficient dynamic range to push the leaked foreground power below the EoR level. By contrast, applying the right Dolph–Chebyshev window can suppress the leaked foreground power below the EoR threshold. For the rest of this paper, we use the DC180 window for all delay spectrum estimation. The Dolph–Chebyshev is the only window tested that has sufficient dynamic range to reach the EoR amplitude without any assumed foreground subtraction,8 and does so with comparable foreground spillover to the Blackman-Harris used in other work (see e.g. Parsons et al. 2012a; Vedantham et al. 2012; Thyagarajan et al. 2013). 4.3 Spectral index comparison The next set of simulations examined the effect of power-law spectral structure, both in the sky model and in the beam. As mentioned in Section 3.2, diffuse power at low frequencies is dominated by Galactic synchrotron emission, which evolves as a power law. Testing the effects of a varied power law sky model gives some sense of how the intrinsic spectrum of the sky affects the power spectrum spillover. Alternatively, simulations with different beam spectral indices can show how much the evolution of the beam’s width couples angular and frequency modes together. Such coupling has been shown to be an important factor in the design of the EDGES experiment (Mozdzen et al. 2016). The angular width of a beam typically evolves linearly with frequency, as for the Airy disc, though antennas for 21 cm experiments are often chosen such that the beam evolves more slowly than that. The Phase I HERA dishes were built with sleeved dipole feeds after the design of PAPER (Ewall-Wice et al. 2016). These are currently being replaced with Vivaldi antenna feeds (Fagnoni et al., in preparation) for Phase II, which have nearly constant width over a wide bandpass (Gibson 1979). In this section, we parametrize the beam’s frequency evolution by letting the width evolve as a power law. Note that the structure is applied to the model and the beam in slightly different ways – the beam is normalized to its peak value at all frequencies, but its width changes with frequency. The sky model’s total brightness is evolving with frequency instead. We ran simulations in the following configurations: Chromatic Gaussian beamwidth power law index α = {−0.25, −0.5, −0.75, −1.0}, with a flat-spectrum sym model. Achromatic Gaussian beams with a sym model with spectral index {α = −3.0, −2.5, −1.5, −0.75}. These simulations were run with the 37 hex antenna layout with the off-zenith times (see Table 1). Power spectra were averaged in time and binned in baseline length. Fig. 6 shows the power spectra of 14.6 m baselines for these different sky and beam models, with Fig. 6(a) showing the power spectra for different beam spectral indices and Fig. 6(b) showing the effects of varying the spectral index of the sky. The black lines show the power spectrum for an achromatic Gaussian beam and spectrally flat sky. A lower panel on each plot shows the ratio of each power spectrum curve to this black line, to better emphasize the effects of the different spectral indices. Figure 6. Open in new tabDownload slide Power spectra versus k∥ for k⊥ = 0.046 Mpc−1 (14.6 m baselines), from simulated data using chromatic Gaussian beams and the sym sky model. In (a), the sky is given a flat spectrum and the beamwidth is allowed to evolve as a power law. In (b), the beam is constant with frequency and the sky brightness evolves with the given spectral index. The black lines in both show the spectrum for an achromatic beam and flat-spectrum sky. The bottom panels on each plot show the fraction of each curve over the black curve, to emphasize the deviations in power spectra for each case. The horizontal black line shows the EoR threshold and the vertical dashed line marks the horizon limit. Figure 6. Open in new tabDownload slide Power spectra versus k∥ for k⊥ = 0.046 Mpc−1 (14.6 m baselines), from simulated data using chromatic Gaussian beams and the sym sky model. In (a), the sky is given a flat spectrum and the beamwidth is allowed to evolve as a power law. In (b), the beam is constant with frequency and the sky brightness evolves with the given spectral index. The black lines in both show the spectrum for an achromatic beam and flat-spectrum sky. The bottom panels on each plot show the fraction of each curve over the black curve, to emphasize the deviations in power spectra for each case. The horizontal black line shows the EoR threshold and the vertical dashed line marks the horizon limit. For the range of spectral indices applied to the beam, the effects on the power spectrum are miniscule. At low k∥, the power is raised slightly, but no more than 1.5× the achromatic level. In the sidelobes of the Dolph–Chebyshev window, the power ripples around the black line but tapers off towards higher k∥. The effect of an increasingly steep spectral index applied to the sky seems to be that it tilts the spectrum more, suppressing the low-k modes and slightly slowing the drop-off in the DC-window sidelobes. In both plots, however, the change in the spillover at the EoR amplitude is less than the separation of the k∥ modes calculated. These results suggest that the monotonic frequency evolution of the antenna and the sky both have a negligible impact on the foreground spillover. Additional frequency structure appearing in the primary beam, due to reflections among components of the analogue system, are non-monotonic, and may have a significant impact on the spillover. These additional effects are included in the CST-simulated beam model. 4.4 Foreground spillover The foreground spillover is estimated by fitting a polynomial spline to the power spectrum measurements and finding the points where the spline crosses the threshold. The spline fit is done in log–log space where it is most stable. The largest k∥ of intersection, less the k∥ of the horizon, is taken as the measured foreground spillover Δk∥, as defined in equation (27). 4.4.1 Time We ran a set of simulations that covered 24 h of LST in steps of 5 min (the ‘transit’ time set), in order to measure how foreground power spillover on single baselines changes with the hour angle of the Galactic centre (recall that the ‘Galactic centre’ also refers to the centre of the sym model). As the Galactic centre moves away from the zenith, the brightest part of the sky is shifted toward a weaker part of the primary beam, but the delay spectrum peak moves closer to the horizon. Brightness near the horizon is further enhanced on all baselines by the pitchfork effect, so it is interesting to see how the power spillover versus time changes with different baseline lengths. Fig. 7 shows the power spectra versus Galactic centre hour angle for a 14.6 m baseline for the symmetric GSM model with Airy (top) and achromatic Gaussian (bottom) beams. The vertical dashed line marks the horizon and the horizontal line marks the chosen threshold Pthresh = 292 mK2 Mpc3. There is a clear increase in total power as the Galactic centre (hence, brightest part of the symmetric sky model) transits overhead. For the aGauss beam, this causes the sidelobe ripples of the Dolph-Chebyshev window to rise up through the threshold, quickly pushing the Δk∥ higher. For the Airy beam, the pitchfork effect is apparent as a persistent bump near the horizon line, and the spillover stays relatively constant except for when the Galactic centre is overhead. Figure 7. Open in new tabDownload slide 1D power spectra for a 38.6 m baseline versus the absolute value of the Galactic centre hour angle. Each curve corresponds with a different hour angle of the Galactic centre, for the achromatic sym model. The vertical dashed line marks the horizon limit, and the horizontal line marks the fiducial EoR amplitude. The top is for the Airy beam; the bottom is for the achromatic Gaussian beam. There Airy beam shows a persistent ‘shoulder’ near the horizon that is due to the foreground pitchfork, while power spectra with the Gaussian beam do not. The pitchfork keeps the spillover at the EoR threshold line very nearly constant for all hour angles, despite the changing amplitude within the horizon. Figure 7. Open in new tabDownload slide 1D power spectra for a 38.6 m baseline versus the absolute value of the Galactic centre hour angle. Each curve corresponds with a different hour angle of the Galactic centre, for the achromatic sym model. The vertical dashed line marks the horizon limit, and the horizontal line marks the fiducial EoR amplitude. The top is for the Airy beam; the bottom is for the achromatic Gaussian beam. There Airy beam shows a persistent ‘shoulder’ near the horizon that is due to the foreground pitchfork, while power spectra with the Gaussian beam do not. The pitchfork keeps the spillover at the EoR threshold line very nearly constant for all hour angles, despite the changing amplitude within the horizon. Fig. 8 shows the measured foreground spillover versus Galactic hour angle for the symmetric GSM (a) and the actual GSM (b) sky models for a 14.6 m baseline (left) and a 25.3 m baseline (right). As will be shown in the next section, the spillover is highest on short baselines because the foreground power is strongest there. Vertical dashed lines mark the rise and set times for the Galactic centre, which in both cases is the brightest part of the sky. The symmetric GSM model shows roughly the same behaviour as the GSM, with a few minor differences. The symmetric GSM is flatter and more symmetric with hour angle, as expected. The GSM is less peaked at hour angles near zero, probably because the full power is more spread out across the sky than in the sym model. For the Airy beam, which has a non-negligible response near the horizon and thus has a consistent pitchfork, the spillover is higher. It is clear that the pitchfork dominates suprahorizon emission during times when the Galaxy is out of the primary beam. The two Gaussian beams show very similar behaviour, with minimal spillover at times when the Galactic centre is away from the zenith. Figure 8. Open in new tabDownload slide The foreground spillover at the threshold shown in Fig. 7, for a 14.6 m baseline (left) and a 25.3 m baseline (right), plotted against the hour angle, for the three analytic beams. Part (a) shows the results with the sym model, while (b) shows the GSM results. Vertical dashed lines mark ±6 h, when the Galactic centre (or, symmetric GSM centre) is rising or setting. The spikes near ≃0 h arises from the Dolph–Chebyshev sidelobe ripples rising above the threshold. Overall, the spillover remains fairly constant with time, with the Airy disc beam showing consistently higher spillover due to the constant presence of the pitchfork. On a longer baseline, the overall amplitude is lower since longer baselines collect less power from diffuse emission. Figure 8. Open in new tabDownload slide The foreground spillover at the threshold shown in Fig. 7, for a 14.6 m baseline (left) and a 25.3 m baseline (right), plotted against the hour angle, for the three analytic beams. Part (a) shows the results with the sym model, while (b) shows the GSM results. Vertical dashed lines mark ±6 h, when the Galactic centre (or, symmetric GSM centre) is rising or setting. The spikes near ≃0 h arises from the Dolph–Chebyshev sidelobe ripples rising above the threshold. Overall, the spillover remains fairly constant with time, with the Airy disc beam showing consistently higher spillover due to the constant presence of the pitchfork. On a longer baseline, the overall amplitude is lower since longer baselines collect less power from diffuse emission. 4.4.2 Baseline length The last set of simulations looked at effect of baseline length on foreground power spillover. For this set we use the Line antenna layout, which provides 300 east–west baselines with lengths spanning 10–150 m. The 150 m baseline here is the longest baseline used in this paper for foreground spillover measurements, and is still within the range of baseline lengths where the point source approximation is good (see Appendix A). These simulations used the off-zenith time set, and power spectra were averaged over that time-span Fig. 9 shows the results with the symmetric GSM (a) and full GSM (b). In both cases, the foreground spillover is dominated by the pitchfork for the Airy beam above about 50 m, and is pushed up by the first sidelobe of the BH window on shorter baselines. For the Gaussian beams, which have no pitchfork, the spillover simply drops off with increasing baseline length. In all cases the spillover decreases with increasing baseline length, albeit more slowly for the Airy beam, due to the decrease in power measured with longer baselines, as is expected for diffuse models. Figure 9. Open in new tabDownload slide The foreground spillover at the threshold shown in Fig. 7, plotted against east–west baseline length, for the three analytic beams. The black line marks the horizon. These simulations were done with the off-zenith time configuration, when the spillover is maximized. Spillover increases steadily with decreasing baseline length. Note that a negative spillover means that foreground power is contained within the horizon limit at the EoR amplitude. Figure 9. Open in new tabDownload slide The foreground spillover at the threshold shown in Fig. 7, plotted against east–west baseline length, for the three analytic beams. The black line marks the horizon. These simulations were done with the off-zenith time configuration, when the spillover is maximized. Spillover increases steadily with decreasing baseline length. Note that a negative spillover means that foreground power is contained within the horizon limit at the EoR amplitude. 4.5 CST beam results The results of Sections 4.4.1 and 4.4.2 show that the foreground power spillover is dominated by the presence of a foreground pitchfork for the Airy disc beam, and is thus (for most Galactic centre hour angles) fairly constant. The analytically defined beams evolve monotonically with frequency, and even extreme levels of monotonic spectral structure were shown to have negligible effect on the spillover in Section 4.3. However, the primary beam of the HERA dishes does not evolve monotonically with frequency, due to reflections among the various physical components of the dish (Ewall-Wice et al. 2016; Patra et al. 2018; Fagnoni et al. 2019). These effects are largely captured by the CST-simulated beam model. Fig. 10 shows the power spectra versus Galactic hour angle for the achromatic sym model, simulated with the CST XX beam and hex layout. The vertical lines in each plot mark the horizon (left) and Nyquist limit of the CST beam (right), while the horizontal line shows the EoR amplitude threshold. The greater spectral structure of the CST beam throws excess power much farther than the analytic beams did (compare, e.g. the left plot of Fig. 10 to the top plot of Fig. 7). The persistent bump on the horizon line for the 38.6 m baseline power spectra shows that the pitchfork does appear for the CST beam model, and appears to set a minimum level of spillover. Figure 10. Open in new tabDownload slide Power spectra versus Galactic centre hour angle for the achromatic sym model simulated with the CST XX beam, averaged over 14.6 m baselines (left) and 38.6 m baselines (right) of the hex layout. The vertical dashed lines in each plot are the horizon and Nyquist limit (left and right, respectively, within each plot), while the horizontal line is the EoR amplitude threshold. The suprahorizon foreground power is not suppressed below the EoR threshold at any time within the Nyquist limit. From these plots, it is clear that the spillover is dominated for most baseline lengths by the existence of a pitchfork. Figure 10. Open in new tabDownload slide Power spectra versus Galactic centre hour angle for the achromatic sym model simulated with the CST XX beam, averaged over 14.6 m baselines (left) and 38.6 m baselines (right) of the hex layout. The vertical dashed lines in each plot are the horizon and Nyquist limit (left and right, respectively, within each plot), while the horizontal line is the EoR amplitude threshold. The suprahorizon foreground power is not suppressed below the EoR threshold at any time within the Nyquist limit. From these plots, it is clear that the spillover is dominated for most baseline lengths by the existence of a pitchfork. Importantly, the suprahorizon power in these plots does not drop below the EoR threshold for k∥ < kNyquist modes within the Nyquist limit of the beam. It is hard to draw any conclusions for the spectra in this region, since the values for k∥ > kNyquist are extrapolated and highly dependent on the choice of frequency interpolation, as shown in Section 3.1.1. The spectra appear to deviate from each other substantially from k ∼ kNyquist on, and it is unclear if this is a consequence of the interpolation or a real effect. To avoid relying on these extrapolated values, in this section we will look at the power spectral amplitude at chosen suprahorizon k∥ modes versus baseline length and Galactic hour angle, rather than try to measure the Δk∥ spillover at the EoR threshold. If the pitchfork is the dominant effect for the CST beam, we should see roughly constant amplitude versus time and baseline length for these modes, as we saw with the Airy beam pitchfork. Fig. 11 shows the amplitudes of different k∥ modes between the horizon and the Nyquist limit for simulations with the CST beam, the hex layout, and the sym or GSM models. Most of the suprahorizon modes stay largely constant in amplitude over time, brightening briefly when the Galactic centre reaches the horizon at ±6 h. This is clearest in the sym simulation results with the 38.6 m baselines. This is likely due to the power from the Galactic centre moving closer to the horizon, adding to the power in the pitchfork. In addition to these bumps near ±6 h, there are peaks around 0 h on the 14.6 m plots on k∥ modes near the horizon limit (blue), due to the increase in power as the Galactic centre crosses the beam centre. This does not appear as strongly with the GSM, and not at all for the sym model, the longer baseline, which is less sensitive to diffuse structure. Figure 11. Open in new tabDownload slide The power spectrum amplitude of various k∥ modes (colour) versus Galactic hour angle for simulated data with the CST XX polarization beam and the sym sky model (a) and GSM (b). The left plots are for the 14.6 m baselines in the hex layout, and the right are for 38.6 m baselines. Only k∥ modes outside the horizon and within the beam’s Nyquist sampling limit are plotted. For most k∥ modes outside the horizon, the power remains largely flat over time, as would be expected from the pitchfork effect. Figure 11. Open in new tabDownload slide The power spectrum amplitude of various k∥ modes (colour) versus Galactic hour angle for simulated data with the CST XX polarization beam and the sym sky model (a) and GSM (b). The left plots are for the 14.6 m baselines in the hex layout, and the right are for 38.6 m baselines. Only k∥ modes outside the horizon and within the beam’s Nyquist sampling limit are plotted. For most k∥ modes outside the horizon, the power remains largely flat over time, as would be expected from the pitchfork effect. Most of the suprahorizon k∥ modes are largely constant in time, with the exception of those particular times, as one would expect from a persistent pitchfork near the horizon. The modes that are close to the Nyquist limit, however, are much more structured with time (bright colours), especially with the sym model. This is reminiscent of the spread of the curves beyond k∥ ∼ 0.1 Mpc−1 in the top plot of Fig. 7 – though the pitchfork itself is largely independent of power within the wedge, the amplitude beyond the pitchfork (the side-ripples of the Dolph–Chebyshev window) varies with the power on sub-horizon modes. In short – the growth in power suggests that the pitchfork region ends before reaching the Nyquist limit. However, it is hard to draw a firm conclusion about behaviour on modes so close to the Nyquist limit in light of the tests in Section 3.1.1. The pitchfork dominance is clear in the 2D power spectra in Fig. 12, made from simulations with varied baseline lengths. The left plot is for the sym model, and the right is for the GSM. In both cases, the brightening along the horizon line (shown by the solid black line) is dominant across all baselines as the source of spillover. Beyond the Nyquist limit, there is clear aliasing that replicates the full pitchfork, showing distinct branches for aliases of the positive and negative delay modes. Discontinuities in the interpolated beam model on scales of the resolution Δν can create sharp peaks in delay space at delay modes 1/Δν. When the beam is convolved with the sky in delay space, this has the effect of replicating the foregrounds at higher delays. This plot suggests, therefore, that even the GPR-smoothed cubic interpolation is not sufficiently smooth for simulating foregrounds beyond the Nyquist limit. A beam model with higher frequency resolution is likely necessary. Figure 12. Open in new tabDownload slide 2D power spectra for simulations with the CST beam, with the sym sky model (left) and GSM (right). The solid black line marks the horizon limit, and the dashed black line marks the Nyquist limit set by the beam model’s frequency resolution. For both sky models, the pitchfork is wide, dominating power beyond the horizon. The entire foreground wedge is aliased beyond the Nyquist limit, causing a second pitchfork to appear. Figure 12. Open in new tabDownload slide 2D power spectra for simulations with the CST beam, with the sym sky model (left) and GSM (right). The solid black line marks the horizon limit, and the dashed black line marks the Nyquist limit set by the beam model’s frequency resolution. For both sky models, the pitchfork is wide, dominating power beyond the horizon. The entire foreground wedge is aliased beyond the Nyquist limit, causing a second pitchfork to appear. 5 SUMMARY Smooth-spectrum foreground sources are expected to be well-behaved in Fourier space, which leaves a clear and well-defined window in which to seek an EoR detection. Spectral structure of the foregrounds, chromatic effects of the baseline and primary beam, and chosen bandpass and window function all contribute to the breakdown of this simple picture and the spread of power into the EoR window. We carried out a series of simulations to characterize and measure this power spread for diffuse models in a variety of controlled cases. We used the GSM, and a symmetrized version of it with varied spectral steepness, in comparison to a fiducial EoR amplitude threshold. The analytic beams, antenna layouts, and frequency channelization and time-step sizes are all chosen to resemble the design of HERA. The choice of a spectral window function has the strongest effect on suppression of foreground power beyond the horizon; windowing along the frequency axis is required to reach the EoR level at all. Applying a window function effectively narrows the bandpass, which means that the main lobe of foreground power is made wider. The lowest k∥ modes are therefore lost to foreground spillover. This trade-off between foreground suppression and spillover is necessary for a detection without substantial foreground removal. For most of potential observing time, the most significant source of suprahorizon power appears to be the foreground pitchfork, which exists for primary beams with a response near the horizon. The pitchfork arises from the shortening of projected baselines near the horizon. As short baselines are more sensitive to diffuse structure, even faint monopole power near the horizon produces a strong response. This response can be observed for long baselines which are normally insensitive to diffuse power. It is worth noting that the simulations used here do not take into account mutual coupling effects among antennas, as discussed in Fagnoni et al. (2019). In a real array, the response of an antenna to sources near the horizon will be affected by the presence of the other antennas that interfere with the propagation of light. To take an extreme case, a bright source exactly in line with both antennas in a baseline will only be visible to the antenna nearer to it, since the farther antenna is in the near one’s shadow. In this case, the correlation between the two antennas is zero. However, it may be that the effect of baseline projection is still significant enough to produce a pitchfork when emission is near the horizon but not in the shadow of other antennas. The full effects of mutual coupling are too sophisticated an effect to model with healvis at this time, and so we leave this to future work. The Airy disc beam and CST simulated power beam both have a strong response near the horizon, and so have a constant foreground spillover due to this pitchfork for most times and baseline lengths. Avoiding LSTs when the Galactic centre is near the beam centre, and choosing an appropriate constant buffer beyond the foreground wedge, appears to be sufficient to avoid diffuse foreground power, absent other systematic effects. It is difficult to draw a conclusion about foreground suppression at the level of the EoR for the CST beam, since those power spectra do not reach the EoR amplitude within the Nyquist limit of the CST simulations. These results suggest that some degree of foreground subtraction is necessary for an avoidance strategy to work. ACKNOWLEDGEMENTS This work was supported by NSF grants 1613040 and 1636646, and NASA grant 80NSSC18K0389. Computational support was provided by the Brown Center for Computation and Visualization (CCV). We thank the developers of scipy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2011), matplotlib, astropy (Astropy Collaboration 2013; Price-Whelan et al. 2018), healpy (Zonca et al. 2019), and pyuvdata (Hazelton et al. 2017). Footnotes 1 " This correspondence is only approximate because the baseline length in wavelengths changes with frequency. Hence, the delay modes are not perfectly orthogonal to the k⊥ modes. 2 " See also Memos #27 and #43 at reionization.org (Parsons 2017; Liu 2018). 3 " https://github.com/rasg-affiliates/healvis 4 " https://www.cst.com 5 " https://github.com/RadioAstronomySoftwareGroup/pyuvdata 6 " Hierarchical equal-area isolatitude pixellization. 7 " https://github.com/telegraphic/PyGSM 8 " The Blackman–Harris window does drop below this flat EoR threshold at the high k∥ for some limited cases, but not enough to be useful for our results. REFERENCES Ali Z. S. et al. ., 2015 , ApJ , 809 , 61 10.1088/0004-637X/809/1/61 Crossref Search ADS Crossref Astropy Collaboration , 2013 , A&A , 558 , A33 10.1051/0004-6361/201322068 Crossref Search ADS Crossref Barry N. et al. ., 2019 , ApJ , 884 , 1 10.3847/1538-4357/ab40a8 Crossref Search ADS Crossref Beardsley A. P. et al. ., 2016 , ApJ , 833 , 102 10.3847/1538-4357/833/1/102 Crossref Search ADS Crossref Bowman J. D. et al. ., 2013 , Publ. Astron. Soc. Aust. , 30 , e031 10.1017/pas.2013.009 Crossref Search ADS Crossref Carroll P. A. et al. ., 2016 , MNRAS , 461 , 4151 10.1093/mnras/stw1599 Crossref Search ADS Crossref Chapman E. et al. ., 2015 , Proc. Sci., Cosmic Dawn and Epoch of Reionization Foreground Removal with the SKA . SISSA , Trieste , 5 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chapman E. , Zaroubi S., Abdalla F. B., Dulwich F., Jelic V., Mort B., 2016 , MNRAS , 458 , 2928 10.1093/mnras/stw161 Crossref Search ADS Crossref Cheng C. et al. ., 2018 , ApJ , 868 , 26 10.3847/1538-4357/aae833 Crossref Search ADS Crossref Datta A. , Bowman J. D., Carilli C. L., 2010 , ApJ , 724 , 526 10.1088/0004-637X/724/1/526 Crossref Search ADS Crossref Dillon J. S. et al. ., 2015 , Phys. Rev. D , 91 , 123011 10.1103/PhysRevD.91.123011 Crossref Search ADS Crossref Ewall-Wice A. et al. ., 2016 , ApJ , 831 , 196 10.3847/0004-637X/831/2/196 Crossref Search ADS Crossref Fagnoni N. , Acedo E. D. L., 2016 , Proceedings of the 2016 International Conference on Electromagnetics in Advanced Applications (ICEAA) . IEEE , Piscataway, New Jersey , p. 629 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Fagnoni N. et al. ., 2019 , preprint (arXiv:1908.02383) Furlanetto S. R. , Peng Oh S., Briggs F. H., 2006 , Phys. Rep. , 433 , 181 10.1016/j.physrep.2006.08.002 Crossref Search ADS Crossref Gibson P. , 1979 , 9th European Microwave Conference Proceedings . IEEE , Brighton, UK , p. 101 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gorski K. M. , Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005 , ApJ , 622 , 759 10.1086/427976 Crossref Search ADS Crossref Harris F. J. , 1978 , Proc. IEEE , 66 , 51 Crossref Search ADS Hazelton B. J. , Jacobs D. C., Pober J. C., Beardsley A. P., 2017 , J. Open Source Softw. , 2 , 140 10.21105/joss.00140 Crossref Search ADS Crossref Kern N. S. , Parsons A. R., Dillon J. S., Lanman A. E., 2019 , ApJ , 884 , 105 10.3847/1538-4357/ab3e73/meta Kohn S. A. et al. ., 2019 , ApJ , 882 , 58 10.3847/1538-4357/ab2f72 Crossref Search ADS Crossref Kolopanis M. et al. ., 2019 , ApJ , 883 , 133 10.3847/1538-4357/ab3e3a Crossref Search ADS Crossref Lanman A. E. , Kern N., 2019 , Astrophysics Source Code Library , record ascl:1907.002 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lanman A. E. , Pober J. C., 2019 , MNRAS , 487 , 5840 10.1093/mnras/stz1639 Crossref Search ADS Crossref Lanman A. , Hazelton B., Jacobs D., Kolopanis M., Pober J., Aguirre J., Thyagarajan N., 2019 , J. Open Source Softw. , 4 , 1234 10.21105/joss.01234 Crossref Search ADS Crossref Liu A. , 2018 , Updated Power Spectrum Normalization . Available at: http://reionization.org/wp-content/uploads/2013/03/Scalar_dev_memo43.pdf Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Liu A. , Shaw J. R., 2020 , PASP , 132 , 062001 10.1088/1538-3873/ab5bfd Liu A. , Tegmark M., 2012 , MNRAS , 419 , 3491 10.1111/j.1365-2966.2011.19989.x Crossref Search ADS Crossref Liu A. , Parsons A. R., Trott C. M., 2014 , Phys. Rev. D , 90 , 023018 10.1103/PhysRevD.90.023018 Crossref Search ADS Crossref Matteo T. D. , Perna R., Abel T., Rees M. J., 2002 , ApJ , 564 , 576 10.1086/324293 Crossref Search ADS Crossref Mellema G. et al. ., 2013 , Exp. Astron. , 36 , 235 10.1007/s10686-013-9334-5 Crossref Search ADS Crossref Mesinger A. , Furlanetto S., Cen R., 2011 , MNRAS , 411 , 955 10.1111/j.1365-2966.2010.17731.x Crossref Search ADS Crossref Morales M. F. , Hazelton B., Sullivan I., Beardsley A., 2012 , ApJ , 752 , 137 10.1088/0004-637X/752/2/137 Crossref Search ADS Crossref Morales M. F. , Beardsley A., Pober J., Barry N., Hazelton B., Jacobs D., Sullivan I., 2019 , MNRAS , 483 , 2207 10.1093/mnras/sty2844 Crossref Search ADS Crossref Mozdzen T. J. , Bowman J. D., Monsalve R. A., Rogers A. E. E., 2016 , MNRAS , 455 , 3890 10.1093/mnras/stv2601 Crossref Search ADS Crossref Oh S. P. , Mack K. J., 2003 , MNRAS , 346 , 871 10.1111/j.1365-2966.2003.07133.x Crossref Search ADS Crossref Paciga G. et al. ., 2011 , MNRAS , 413 , 1174 10.1111/j.1365-2966.2011.18208.x Crossref Search ADS Crossref Paciga G. et al. ., 2013 , MNRAS , 433 , 639 10.1093/mnras/stt753 Crossref Search ADS Crossref Parsons A. , 2017 , Power Spectrum Normalizations for HERA . Available at: http://reionization.org/wp-content/uploads/2013/03/Power_Spectrum_Normalizations_for_HERA.pdf Parsons A. R. et al. ., 2010 , AJ , 139 , 1468 10.1088/0004-6256/139/4/1468 Crossref Search ADS Crossref Parsons A. , Pober J., McQuinn M., Jacobs D., Aguirre J., 2012a , ApJ , 753 , 81 10.1088/0004-637X/753/1/81 Crossref Search ADS Crossref Parsons A. R. , Pober J. C., Aguirre J. E., Carilli C. L., Jacobs D. C., Moore D. F., 2012b , ApJ , 756 , 165 10.1088/0004-637X/756/2/165 Crossref Search ADS Crossref Parsons A. R. et al. ., 2014 , ApJ , 788 , 106 10.1088/0004-637X/788/2/106 Crossref Search ADS Crossref Patil A. H. et al. ., 2017 , ApJ , 838 , 65 10.3847/1538-4357/aa63e7 Crossref Search ADS Crossref Patra N. et al. ., 2018 , Exp. Astron. , 45 , 177 10.1007/s10686-017-9563-0 Crossref Search ADS Crossref Pedregosa F. et al. ., 2011 , J. Mach. Learn. Res. , 12 , 2825 Planck Collaboration XIII , 2016 , A&A , 594 , A13 10.1051/0004-6361/201525830 Crossref Search ADS Crossref Pober J. C. et al. ., 2013 , ApJ , 768 , L36 10.1088/2041-8205/768/2/L36 Crossref Search ADS Crossref Pober J. C. et al. ., 2014 , ApJ , 782 , 66 10.1088/0004-637X/782/2/66 Crossref Search ADS Crossref Presley M. , Liu A., Parsons A., 2015 , ApJ , 809 , 18 10.1088/0004-637X/809/1/18 Crossref Search ADS Crossref Price-Whelan A. M. et al. ., 2018 , AJ , 156 , 123 10.3847/1538-3881/aabc4f Crossref Search ADS Crossref Pritchard J. R. , Loeb A., 2008 , Phys. Rev. D , 78 , 103511 10.1103/PhysRevD.78.103511 Crossref Search ADS Crossref Rabiner L. R. , Gold B., 1975 , Theory and Application of Digital Signal Processing . Prentice-Hall , Englewood Cliffs, New Jersey Rasmussen C. E. , Williams C. K. I., 2006 , Gaussian Processes for Machine Learning . MIT Press , Cambridge, MA Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sims P. H. , Lentati L., Alexander P., Carilli C. L., 2016 , MNRAS , 462 , 3069 10.1093/mnras/stw1768 Crossref Search ADS Crossref Smith J. O. , 2019 , Spectral Audio Signal Processing . Available at: http://ccrma.stanford.edu/ jos/sasp/ Thompson A. R. , Moran J., Swenson G. W. Jr, 2017 , Interferometry and Synthesis in Radio Astronomy , 3rd edn. Springer International Publishing , Switzerland 10.1007/978-3-319-44431-4 Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Crossref Thyagarajan N. et al. ., 2013 , ApJ , 776 , 6 10.1088/0004-637X/776/1/6 Crossref Search ADS Crossref Thyagarajan N. et al. ., 2015a , ApJ , 804 , 14 10.1088/0004-637X/804/1/14 Crossref Search ADS Crossref Thyagarajan N. et al. ., 2015b , ApJ , 807 , L28 10.1088/2041-8205/807/2/L28 Crossref Search ADS Crossref Thyagarajan N. , Parsons A. R., DeBoer D. R., Bowman J. D., Ewall-Wice A. M., Neben A. R., Patra N., 2016 , ApJ , 825 , 9 10.3847/0004-637X/825/1/9 Crossref Search ADS Crossref Trott C. M. , Wayth R. B., Tingay S. J., 2012 , ApJ , 757 , 101 10.1088/0004-637X/757/1/101 Crossref Search ADS Crossref Vedantham H. , Udaya Shankar N., Subrahmanyan R., 2012 , ApJ , 745 , 176 10.1088/0004-637X/745/2/176 Crossref Search ADS Crossref Virtanen P. et al. ., 2020 , Nat. Methods , 17 , 261 10.1038/s41592-019-0686-2 Zheng H. et al. ., 2017 , MNRAS , 464 , 3486 10.1093/mnras/stw2525 Crossref Search ADS Crossref Zonca A. , Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019 , J. Open Source Softw. , 4 , 1298 10.21105/joss.01298 Crossref Search ADS Crossref APPENDIX: POINT SOURCE APPROXIMATION The simulator used in this paper assumes that variations in the primary beam, brightness, and interferometric fringe are negligible on the scale of a single pixel, such that the contribution of each pixel can be treated as that of a point source at the pixel centre. This appendix derives a rough bound on the error due to this point source approximation (PSA) for the sky and beam models used in this paper. A forthcoming paper will explore the validity of this simulator against analytic solutions for diffuse sky models, and will delve into more detail on the limitations of the PSA. We test the validity of our simulator here by performing a set of monopole simulations at different resolutions with an antenna layout featuring a broad range of baseline lengths. Our test array consists of 80 randomly distributed antennas, forming baselines with lengths spanning 2−608 m, with the Airy beam and a bandpass of 100−120 MHz at 97 kHz channel width. The monopole sky makes for a useful test because it should minimize the power within the foreground wedge, which makes any erroneous simulated power stand out. Fig. A1 shows the cylindrically binned power spectra for these three simulations, for maps with Nside 128, 256, and 512 from top to bottom. Note that the longest baseline used in this paper, for the ‘spillover versus baseline length’ results, was only 50λ. The pitchfork feature is clearly visible in each plot around the horizon (marked with the dashed line). At Nside 128, there is excess power on long baselines (equivalently, at large k⊥) which disappears at Nside 256. There is marginal improvement from raising the resolution again to Nside 512. Figure A1. Open in new tabDownload slide Power spectra for a monopole signal on healpix maps with Nside 128, 256, and 512 (from top to bottom). The top axis shows the baseline length corresponding with the k⊥ mode, and the dashed line marks the horizon. Excess power within the foreground wedge is erroneous for a monopole signal, and is clearly visibible at larger k⊥ for the Nside 128 simulations. Clearly errors are reduced by switching from Nside 128 to 256, but further resolution improvement does not yield much noticeable improvement. Note that the longest baseline used for any results in this paper is at 50λ, for the spillover versus baseline length tests. Further, the foreground pitchfork is not strongly influenced by the resolution of the map, as expected since it is a result of the projected baseline lengths being small near the horizon. Figure A1. Open in new tabDownload slide Power spectra for a monopole signal on healpix maps with Nside 128, 256, and 512 (from top to bottom). The top axis shows the baseline length corresponding with the k⊥ mode, and the dashed line marks the horizon. Excess power within the foreground wedge is erroneous for a monopole signal, and is clearly visibible at larger k⊥ for the Nside 128 simulations. Clearly errors are reduced by switching from Nside 128 to 256, but further resolution improvement does not yield much noticeable improvement. Note that the longest baseline used for any results in this paper is at 50λ, for the spillover versus baseline length tests. Further, the foreground pitchfork is not strongly influenced by the resolution of the map, as expected since it is a result of the projected baseline lengths being small near the horizon. The pitchfork is largely independent of map resolution, which is further evidence that it is caused by the projection of baselines near the horizon. The error in the point source approximation is worse for longer baselines, but the power in the pitchfork comes from the shortening of baselines near the horizon. We note also the presence of faint vertical striping, and non-zero power within the wedge for all resolutions. The exact visibility for a monopole sky signal with an Airy beam is given by the convolution of a circular disc (of diameter 14 m, in this case) with a sinc function of the baseline length. We suspect that the sidelobes of this sinc function are the cause of the stripes, as well as the broad spread of power at small k⊥. © 2020 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Quantifying EoR delay spectrum contamination from diffuse radio emission JF - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/staa987 DA - 2020-05-21 UR - https://www.deepdyve.com/lp/oxford-university-press/quantifying-eor-delay-spectrum-contamination-from-diffuse-radio-AWTZBFC0Z3 SP - 3712 VL - 494 IS - 3 DP - DeepDyve ER -