THE HYDROGEN EPOCH OF REIONIZATION ARRAY DISH. II. CHARACTERIZATION OF SPECTRAL STRUCTURE WITH ELECTROMAGNETIC SIMULATIONS AND ITS SCIENCE IMPLICATIONS

The time-dependent electromagnetic behavior of the antenna was modeled using the Microwave Studio, a commercial numerical simulation package produced by Computer Simulation Technology (CST). The model consists of an idealized 14 m diameter paraboloid reflector with a feed structure placed at the focus (4.5 m above the surface). The feed is a dual linear polarized sleeved dipole identical to that used as an element of PAPER but with a cylindrical skirt used in place of the angled planar reflectors. The dipole arms are modeled as copper tubing with aluminum disks and reflector surfaces. For simulation, the entire model was encapsulated in vacuum dielectric with radiative boundaries and discrete ports of 125 Ω impedance were defined at the terminals of the orthogonal dipoles. The simulation is conducted in a box encapsulating the dish geometry plus three wavelengths (at 100 MHz) on all sides with open boundary conditions.

The Microwave Studio divides the model into approximately 94 million discrete hexahedral cells, where the field is calculated over time using the Finite Integration Technique incorporating the Perfect Boundary Approximation¹⁴ (Workflow and Solver Overflow Document, CST Microwave Studio, 2015, Chapter 3) in response to a 100–200 MHz broadband pulse of about 40 ns in duration. The model was excited in two ways: (1) via an E–W polarized plane wave entering the model along the bore sight, and (2) from one of the discrete terminals (whose analysis we focus on in Section 4.5). The voltage responses at the terminals are monitored for these two cases over a duration of 500 ns as the excitation pulse travels throughout the model.

In Figure 3, we show the geometry of the electromagnetic simulation. A 150 MHz plane wave with a Gaussian envelope is initialized above the dish vertex traveling in the −z direction. It is reflected by the dish before entering the dipole feed, hidden below the cylindrical skirt in this figure. We record the electric field voltage of the plane wave at the feed output terminals as a function of time, plotted as a red line in Figure 3 along with the output voltage at the terminals for the feed polarized parallel to the plan wave (black line). The delay between the central envelope of the plane wave and the voltage output of the feed is ≈30 ns, which corresponds to the round trip travel time from the feed to the dish vertex and back. However, while the input plane wave, modulated by a Gaussian, falls off rapidly after the first ≈20 ns after its peak, we see that the voltage output decays far more slowly due to reflections and resonances within the antenna structure. We are able to get a qualitative feel for the amplitude of the reflections by inspecting the falloff of the time domain voltage response and see that after 60 ns it reaches ≈−25 dB.

Zoom In Zoom Out Reset image size

We can do much better than this. From Equation (1), we know that the voltage output results from the convolution of the plane wave input with the voltage gain of the antenna. Since we know the input wave, a straightforward application of the Fourier convolution theorem allows us to determine the voltage response.

Since our simulation is sampled in finite time steps, we will adopt a discretized notation for this section. In particular, our simulation consists of N samples, evenly spaced by dτ at times τ_n = n × dτ. We denote the output voltage at the feed terminals at time τ_n as ${\widetilde{v}}_{n}$. Rewriting the convolution in Equation (1) in discrete notation, we have

$$\begin{eqnarray}&&{\widetilde{v}}_{n}(\widehat{{\boldsymbol{s}}})=\displaystyle \sum _{m}\,{\widetilde{r}}_{m}(\widehat{{\boldsymbol{s}}}){\widetilde{s}}_{n-m}(\widehat{{\boldsymbol{s}}}).\end{eqnarray} \tag{ 10 }$$

We may undo this convolution by taking a discrete Fourier transform (DFT) of both $\widetilde{{\boldsymbol{v}}}$ and $\widetilde{{\boldsymbol{s}}}$ in time, dividing them in Fourier space, and taking an inverse DFT back. Symbolically,

$$\begin{eqnarray}&&\widetilde{r}(\widehat{{\boldsymbol{s}}})={{\boldsymbol{ \mathcal F }}}^{-1}\left[\displaystyle \frac{{\boldsymbol{ \mathcal F }}\widetilde{{\boldsymbol{v}}}(\widehat{{\boldsymbol{s}}})}{{\boldsymbol{ \mathcal F }}\widetilde{s}(\widehat{{\boldsymbol{s}}})}\right],\end{eqnarray} \tag{ 11 }$$

where ${\boldsymbol{ \mathcal F }}$ is the Fourier transform matrix for a 1D vector of length N.

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal F }}}_{{mn}}={e}^{-2\pi {imn}/N}.\end{eqnarray} \tag{ 12 }$$

In Figure 4, we show the amplitude of the Fourier transform of our Gaussian input, centered at 150 MHz along with the voltage response. Since our input only has support between ≈20 and 280 MHz, the direct ratio of our voltage response and input wave is dominated by numerical noise outside of this range. We eliminate these artifacts by multiplying our ratio by a Blackman–Harris window between 100 and 200 MHz and set our estimate to zero elsewhere. From a physical standpoint, this is sensible since 21 cm experiments only observe a limited bandwidth. PAPER's correlator, which will initially serve as the HERA backend samples over a 100 MHz instantaneous frequency interval and analog filtering is applied to prevent aliasing.

Zoom In Zoom Out Reset image size

Applying Equation (11) to our simulation, we obtain estimates of the time-domain voltage response of the HERA dish toward zenith, which we plot in Figure 5. We also conduct a time-domain simulation of the voltage output in response to an identical input plane wave for the skirted dipole PAPER antenna (pictured above the HERA dish in Figure 1) in order to determine whether the presence of the parabolic dish introduces reflections and spectral structure in excess of previous successful antenna designs.¹⁵ We inspect the absolute value of $\widetilde{v}$ for both the PAPER and HERA antennas in Figure 5. Since non-zero values of $\widetilde{v}$ at negative delays violates causality, we assume such features are sourced by artifacts such as sidelobes of the zero-delay peak and/or numerical precision noise, which sets a limit of ≈−60 dB on the dynamic range of our method. We first note peaks in the HERA curve that are spaced by ≈35 ns and are absent in the PAPER simulation. Another significant difference between the two curves is a knee in the HERA gain at ≈120 ns that is not present in PAPER, leading to an increase in gain of approximately 20 dB at 200 ns. Because this knee does not exhibit nodes at 35 ns intervals, it is probably not sourced by reflections in the antenna-feed geometry but by reflections involving a geometric length that is not resolved by the 100 MHz bandwidth of the simulation.

Zoom In Zoom Out Reset image size

Our next step is to compute the power kernel, $\widetilde{R}$, given by Equation (7) by performing a convolution of the time-reversed voltage response with the complex conjugate voltage response. In Figure 6, we show the resulting power kernel for PAPER and HERA. Since both voltage gains drop rapidly with increasing delay, the approximation in Equation (9) holds quite well and we see that the kernels fall off at a rate similar to the voltage response.

Zoom In Zoom Out Reset image size

Because the 21 cm brightness temperature fluctuations evolve over redshift intervals of Δz ≳ 0.5 (Zaldarriaga et al. 2004), experiments will sub-divide their bands into ≈10 MHz intervals for power spectrum estimates at multiple redshifts. Thus the localization of reflections within the HERA band will determine whether all or only a subset of redshifts are affected. To determine whether the reflections are localized in frequency, we compute the voltage delay response and power kernel for three different subbands: 100–130 MHz, 130–160 MHz, and 160–190 MHz. In order to maintain a decent resolution of the kernel itself, we use frequency ranges that are several times larger than the actual subbands that will be used for EoR power spectrum estimation (≈10 MHz). The resulting power kernels are plotted in Figure 7. The central lobe of the subband kernels is significantly wider due to the the wider window functions incurred by the reduced bandwidth. We find that the shallow long-term falloff is only visible within the central 130–160 MHz band, indicating that long-term reflections are isolated around 150 MHz.

Zoom In Zoom Out Reset image size

To further illustrate the isolation of fine frequency structure in the center of the bandpass and to verify that our observations are not an artifact of our reduction of the simulation outputs, we fit 10 MHz intervals of the absolute value of the simulated gains to a sixth order polynomial and inspect the residuals in Figure 8. We find that the gain residuals on the sixth order fit are an order of magnitude greater on the 145–155 MHz subband than any other frequency interval.

Zoom In Zoom Out Reset image size

In Figure 5, we see that the long-term delay response of the HERA antenna differs from the PAPER design in two regards, the existence of node-like structures spaced every ∼35 ns, corresponding to the round trip delay between the dish and the feed, and a long-term knee that dominates the response function after 100 ns but does not exibit any periodicity that might be associated with the feed-dish geometry. Here we establish the origins of the long timescale structure using simulations. Reflectometry measurements conducted on the isolated feed in several different configurations (N. Patra et al. 2016, in preparation) independently discover and verify these origins. To determine whether this knee-like structure originates from feed-dish reflections, we perform a simulation of a plane wave incident along the bore sight the HERA feed without the dish and compare it to the time-domain response for the feed-and-dish system in Figure 10. The feed-only simulation lacks the reflections associated with the dish and feed configuration (present at small delays) but retains the large-delay knee. Thus, the dominant large timescale contamination of the HERA bandpass is not intrinsic to the arrangement of a feed suspended over the dish.

We isolate the source of the knee within the feed structure (shown in Figure 9) by running simulations of the dipole feed with only a backplane along with a sleeved dipole in isolation. We show both of these voltage responses in Figure 10. When the cylindrical skirt is removed from the feed, the knee vanishes. Because of the narrow-band nature of the knee (Section 4.2), its presence only when the cylinder is attached, and its exponential falloff (characteristic of a damped harmonic oscillator), we conclude that the cylindrical skirt is behaving as a resonant cavity and the resulting stored oscillating fields are responsible for the majority of the fine-scale spectral structure in the HERA antenna's gain. We estimate the quality factor of the resonance to be ≈6.5 from the time constant of the exponential falloff, which yields a resonance width of ≈20 MHz.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the termination of the dipole could be significantly different from the termination used in the simulation (100 Ω), we also investigate the impact of the termination impedance of the dipole on the delay response. we vary the termination impedance between 50 and 500 Ω for several plane-wave simulations and show the resulting voltage responses in Figure 11. We find that changing the termination impedance has a significant impact on the time-domain response below ≈150 ns in the region that is dominated by dish-feed reflections. Since the termination impedance does not affect resonance of the feed cavity, its effect is small beyond ≈200 ns. Only in the extreme, 500 Ω case do the reflections extend to appreciable delays. Since delays over ≈250 ns are responsible for contaminating the EoR window, termination impedance will only impact HERA's sensitivity for an extremely poor match.

Zoom In Zoom Out Reset image size

Simulations combining the dish with the analog signal-chain show that matching networks can reduce the levels feed-dish reflections below those observed in this paper, even under realistic termination conditions (Fagnoni & de Lera Acedo 2016).¹⁶

Since the feed properties dominate the antenna's long-term response, we investigate the voltage response (Figure 5) and the power kernel (Figure 6) of the HERA antenna without the cylindrical skirt and a bare dipole. Removing the cylinder leads to a ≈ 10 dB improvement in the antenna response after 100 ns and eliminating the backplane as well yields a 20 to 30 dB improvement. Hence, the delay response of the HERA antenna can potentially be improved by modifying the feed. There are several tradeoffs in polarization and sensitivity that are made when we remove feed components, which we now discuss.

We determined above that the long-delay knee in the gain of the HERA dish is primarily caused by resonance in the cylindrical skirt and backplane that surround the PAPER dipole. These structures are added to enhance the symmetry of its beam, to improve the effective area of the antenna, and reduce cross coupling (DeBoer et al. 2016). It is, therefore, worth examining the impact that removing these components has on the properties of the HERA beam.

We first examine the trade-off that is made in polarization when we remove elements of the original design. Faraday rotation is capable of imparting a fine-scale frequency structure into the Q/U polarized Stokes visibilities. Antennas in which the X and Y polarized beams are not identical can leak stokes Q/U visibilities and the spectral structure they contain into the Stokes-I visibilities from which delay power spectra are formed, potentially masking the 21 cm signal (Jelić et al. 2010; Moore et al. 2013, 2015). It is shown in Moore et al. (2013) and Moore et al. (2015) that this leakage is given approximately by

$$\begin{eqnarray}&&{P}_{{leak}}\approx \displaystyle \frac{{{ \mathcal A }}_{-}}{{{ \mathcal A }}_{+}}{P}_{P}\equiv \xi {P}_{P},\end{eqnarray} \tag{ 13 }$$

where P_P is the linearly polarized power spectrum, and

$$\begin{eqnarray}&&{{ \mathcal A }}_{\pm }=\int d{\rm{\Omega }}| {A}_{{xx}}(\widehat{{\boldsymbol{s}}})\pm {A}_{{yy}}(\widehat{{\boldsymbol{s}}}){| }^{2},\end{eqnarray} \tag{ 14 }$$

where A_xx is the x polarized beam and A_yy is the the y polarized beam. In Table 1, we give ξ at 150 MHz, computed for the three configurations of the HERA feed discussed in Section 4.3. We see that the feed with the cylinder has the smallest ellipticity with ξ ≈ 3 × 10⁻³ over much of the band. Polarization measurements from Asad et al. (2015) indicate a polarized power spectra of P_P ≈ 10³ mk² within the EoR window. Thus, an antenna with ξ ≈ 3 × 10⁻³ would produce leakage of several mk², which is roughly an order of magnitude below typical reionization peak amplitudes of several tens of mK². However, the feed variant in which the cylinder has been removed has an ellipticity that is roughly ≈4× larger and a bare dipole, ≈10× larger. These enhanced ellipticities would bring the level of the polarization leakage to or above the level predicted in Asad et al. (2015). On the other hand, it is possible that the LST binning and averaging over many nights with different rotation measures can suppress polarization leakage by a factor of ≈10³ in the power spectrum (Moore et al. 2015).

Table 1.  Ellipticity and Effective Area of the HERA Antenna for Different Feed Configurations

Feed Configuration	Ellipticity, ξ	Effective Area, ${A}_{\mathrm{eff}}$, (m2)
Cylindrical and Backplane	3.5 × 10−3	56.7
No Cylinder	1.5 × 10−2	97.2
No Cylinder or Backplane	2.8 × 10−2	41.9

Download table as:  ASCIITypeset image

The second impact that we consider is the effective area of the antenna, given by (Wilson et al. 2009),

$$\begin{eqnarray}&&{A}_{\mathrm{eff}}=\displaystyle \frac{{\lambda }^{2}}{\int d{\rm{\Omega }}A(\widehat{{\boldsymbol{s}}})}.\end{eqnarray} \tag{ 15 }$$

In Table 1, we see that the effective area at 150 MHz is actually enhanced by the removal of the cylinder but is negatively impacted by the removal of the backplane, which introduces a significant back-lobe, reducing the beam's directivity. As detailed in Neben et al. (2016), this reduction in the effective area of the dish lowers HERA's sensitivity by a factor of the order of unity.

The development of the HERA feed is still underway and it is possible that improvements in the feed will allow for the 150 MHz resonance to be eliminated while maintaining the illumination and beam-symmetry of the original feeds.

We now assess the accuracy of our time-domain simulation framework by comparing simulations to measurements of the S₁₁ parameter of the HERA dish. Until now, our simulations have derived the voltage response of the dish using simulations of an incoming plane wave as is the case for radio signals arriving from objects at cosmological distances. It is possible to probe the gain of the dish using objects in the far field such as known radio sources (Thyagarajan et al. 2011; Pober et al. 2012; Colegate et al. 2015) or constellations of ORBCOMM satellites (Neben et al. 2015, 2016). However, natural radio sources are too weak to probe the dish response at the ≲10⁻⁴ level necessary to verify our simulations and the ORBCOMM technique can only be used to map the gain at the 137 MHz ORBCOMM transmission frequency. Work is currently underway to use broadband transmitters flown into the far field of the dish by drones (D. Jacobs et al. 2016, in preparation) but this system is still under development. Reflectometry of the dish with a VNA is a straightforward alternative used in N. Patra et al. (2016, in preparation) to estimate the gains directly. Rather than comparing their estimate of the gain with our predicted gains (which is done in their paper), we set up a time domain simulation to compute the S₁₁ parameter of HERA antenna and compare it to direct S₁₁ measurements.

S₁₁ refers to the complex ratio between a voltage signal transmitted into the feed terminals and the voltage reflected back as a function of frequency, ${S}_{11}(f)\equiv {v}_{\mathrm{trans}}(f)/{v}_{\mathrm{recv}}(f)$. Measurements, described in further detail in N. Patra et al. (2016, in preparation), of S₁₁ were taken on a prototype HERA dish at the National Radio Astronomy Observatory Green Bank facility using an Agilent 8753D VNA. The VNA was connected to the antenna's balanced feed terminal via a 30.5 m length of RG-8X-LL coaxial cable and M/A Com HH-128 180 degree hybrid junction. The VNA was calibrated at the balanced end of the hybrid junction using a set of termination standards having SMA connectors. A small adapter from SMA to the post terminals on the feed was not included in the calibration. The VNA excites the terminals with a band-limited (100–200 MHz) pulse, and the complex reflected voltage at the calibration plane is measured as a function of frequency. The reflection for the case of open terminals was measured to confirm the dynamic range and resolution of the measurement. While our simulations are terminated with an impedance of 125 Ω, the VNA measurements are terminated at 100 Ω, which only has a small effect on the large-delay response of the antenna (Figure 11).

Because S₁₁ is defined as a ratio in the frequency domain in a form identical to our voltage gains, we may run time-domain simulations similar to those described in Section 3.1 except, rather than simulating an incoming plane wave, we simulate the excitation of the feed terminals by a delta-gap impedance port and record output voltage as a function of time. We calculate S₁₁ from the simulation in the same manner that we obtained $\widetilde{r}$ using Equation (11) with ${\widetilde{v}}_{\mathrm{recv}}$ taking the place of $\widetilde{v}$ in the numerator inside the Fourier transform and ${\widetilde{v}}_{\mathrm{trans}}$ taking the place of $\widetilde{s}$ in the denominator.

In Figure 12, we show the simulated amplitude of S₁₁ as a function of frequency for the HERA dish, observing a distinctive two peaked structure before a steep die off in delay that transitions to a shallower falloff at ≈150 ns. The first peak is due to the reflection of the input wave off of the back of the feed while the second, roughly 35 ns later, arises from the transmitted component of the input wave reflecting off of the dish and arriving back at the feed. The ensuing long-term die off arises from reflections within the feed and dish structure and, for reasons that will be elaborated on in N. Patra et al. (2016, in preparation), corresponds very closely to our simulations of the dish gain itself (compare with Figure 5). We get a sense of the dynamic range of the measurement by unhooking the SMA adapter that attaches the sleeved dipole feed to the cable from our VNA, forming an open circuit that should ideally give a reflection coefficient of ≈1 at zero delay with no reflections at any other times. In this measurement, we find a noise-like structure at ≈−50 dB. Below ≈500 ns, we found that this structure does not integrate down with time, leading us to conclude that it is caused by systematics, likely uncalibrated low-level reflections in the VNA-feed cables. The level of this noise sets the systematic floor in our measurements, which we show as a gray shaded region in Figure 12. We find that in the region where the S₁₁ measurement is above the systematics floor, there is good agreement with our simulations (within several dB).

Zoom In Zoom Out Reset image size

Our simulations of the dish response only extend to ≈400 ns. However, interferometers are expected to be most sensitive to the 21 cm power spectrum at comoving scales between k ≈ 0.1–0.5 h Mpc⁻¹, corresponding to the delays between 180 and 900 ns at z = 8. To extrapolate out to the higher delays in this range, we assume that the response function is dominated by a sum of field oscillations within the feed's cylindrical skirt and reflections between the feed and the dish.

$$\begin{eqnarray}&&\widetilde{r}={\widetilde{r}}_{o}+{\widetilde{r}}_{r}\end{eqnarray} \tag{ 16 }$$

The long-term falloff after the knee, which is dominated by oscillations within the cylinder, appears as a line on a linear-log plot (Figure 5). This is the response we would expect for a damped harmonic oscillator. We thus model the voltage response from the cylindrical skirt as an exponential

$$\begin{eqnarray}&&{\widetilde{r}}_{o}={A}_{o}{X}_{o}^{-\tfrac{\tau 150\,\mathrm{MHz}}{2Q}}\end{eqnarray} \tag{ 17 }$$

where Q is the quality factor of the cavity, which we determine to be ≈6.5.

We see in Figure 5 that when the cylinder and backplane are removed, the response of the antenna falls off much faster, in line with the response below ≈130 ns. Like the response from the cylinder, this falloff also goes exponentially, which we now show is also indicative of reflections.

We let Γ_d represent the reflection coefficient of the dish vertex and Γ_f represent the reflection coefficient of the feed in the presence of each other. An electromagnetic wave incident on the feed, at t = 0, is accepted with amplitude $(1+{{\rm{\Gamma }}}_{f})$. The reflected component travels back to the dish and acquires an amplitude of (Γ_fΓ_d) before returning at time, τ_d, later where $(1+{{\rm{\Gamma }}}_{f})$ will be accepted and Γ_f will be reflected back toward the dish. Summing the infinite series of reflections, the time dependent voltage at the feed is

$$\begin{eqnarray}&&{\widetilde{v}}_{r}(t)=(1+{{\rm{\Gamma }}}_{f})\displaystyle \sum _{m}{({{\rm{\Gamma }}}_{f}{{\rm{\Gamma }}}_{d})}^{m}\widetilde{s}(t-m{\tau }_{d}),\end{eqnarray} \tag{ 18 }$$

which implies that

$$\begin{eqnarray}&&{\widetilde{r}}_{r}(\tau )=(1+{{\rm{\Gamma }}}_{f})\displaystyle \sum _{m}{({{\rm{\Gamma }}}_{f}{{\rm{\Gamma }}}_{d})}^{m}{\delta }_{D}(\tau -m{\tau }_{d}).\end{eqnarray} \tag{ 19 }$$

In practice, each delta function in Equation (21) is convolved with the window function arising from the finite bandwidth of the instrument. For our simulations, which have a bandwidth of 100 MHz, the width of the window function is ≈10 ns, partially filling in the gaps between each reflection. Since the number of reflections is m = t/τ_d, than we can write the long-term delay response in discrete form as

$$\begin{eqnarray}&&\widetilde{r}{({nd}\tau )\approx ({{\rm{\Gamma }}}_{f}{{\rm{\Gamma }}}_{d})}^{{nd}\tau /{\tau }_{d}}\end{eqnarray} \tag{ 20 }$$

which is an exponential in time. We thus model the response due to reflections between the feed and the dish as an exponential

$$\begin{eqnarray}&&{\widetilde{r}}_{r}={A}_{r}{X}_{r}^{(\tau /30\mathrm{ns})}.\end{eqnarray} \tag{ 21 }$$

To extrapolate beyond 400 ns, we fit the voltage response to the sum of two exponentials for ${\widetilde{r}}_{o}$ and ${\widetilde{r}}_{r}$ and set the feed response beyond the maximum time our simulation to this power-law sum. In our pedagogical treatment, we have assumed that the reflection coefficients are frequency independent, which is not the case in real life. The impact of frequency dependent reflection coefficients is to replace the summands in Equation (19) with the Fourier transform of ${({{\rm{\Gamma }}}_{f}{{\rm{\Gamma }}}_{d})}^{m}$ evaluated at $t=(\tau -m{\tau }_{d})$. As long as ${({{\rm{\Gamma }}}_{f}{{\rm{\Gamma }}}_{d})}^{m}$ is compact in delay space, which is the case if the reflection coefficients evolve smoothly with frequency, than our pedagogical approximation still yields the power law in Equation (21).

Given the chromaticity due to reflections and resonances in the HERA antenna, what Fourier modes will still be accessible with the delay filtering technique? To answer this question, we combine our extrapolated simulations of the HERA dish's spectral structure with simulations of foregrounds.

The foreground model is discussed in detail in Thyagarajan et al. (2015b), but we describe them briefly here for the reader's convenience. It consists of two major components: diffuse synchrotron emission from our Galaxy, the structure of which is described by the Global Sky Model (GSM) of de Oliveira-Costa et al. (2008), and a population of point sources including objects from the NRAO Sky Survey (NVSS; Condon et al. 1998) at 1.4 GHz and the Sydney University Mologolo Sky Survey (SUMMS; Bock et al. 1999) at 843 MHz. Their fluxes are extrapolated to the observed 100–200 MHz band using a spectral index of $\langle \alpha \rangle =-0.83$ determined in Mauch et al. (2003). Visibilities are computed from the diffuse and point-source models assuming an achromatic angular response for the HERA beam at 137 MHz described in Neben et al. (2016). We compute two sets of visibilities: one in which the spectral structure of the dish is assumed to be completely flat and another in which the beam is multiplied by the frequency dependent gain at zenith determined by our simulations.

The foreground filtering procedure employed by PAPER and HERA involves delay transforming the visibilities and performing a 1D CLEAN (Parsons & Backer 2009; Parsons et al. 2012b), which discovers and subtracts foregrounds within the horizon plus a small buffer, allowing for the suppression of foreground sidelobes in delay space. The level of foreground subtraction possible by this procedure is limited by the thermal noise level on the visibility, which in turn depends on the number of time steps and redundant baselines that are averaged before performing the CLEANing step. In this work, we assume that each visibility is CLEANed independently with a 20 minute cadence. Each short baseline is coherent for ≈10 s per night (Ali et al. 2015), so this is equivalent to each local sidereal time (LST) being integrated for ≈120 nights before CLEANing. If we were to average over redundant baselines this number would drop enormously. The standard deviation on the real and imaginary part of a single delay transformed visibility is (Morales & Hewitt 2004)

$$\begin{eqnarray}&&{\sigma }_{V}=\displaystyle \frac{\sqrt{2B}{k}_{{\rm{B}}}{T}_{\mathrm{sys}}}{{A}_{e}\sqrt{\tau }},\end{eqnarray} \tag{ 22 }$$

where A_e is the effective area of the dish, B is the bandwidth, T_sys is the system temperature, τ is the integration time, and k_B is the Boltzmann constant. For T_sys, we use the equation ${T}_{\mathrm{sys}}=100{\rm{K}}+{T}_{\mathrm{sky}}$, where 100 K is the temperature of the PAPER receiver and ${T}_{\mathrm{sky}}=60{(\lambda /1{\rm{m}})}^{2.55}$ is the sky temperature (Rogers & Bowman 2008; Fixsen et al. 2011). For A_e, we use the value of 98 m² determined in Neben et al. (2016). We CLEAN down to five times the thermal noise level. Subtracting a model of the sky is an alternative to CLEANing that could reach below the noise level on a single visibility but also assumes a priori knowledge of the foregrounds and the instrumental response. In Figure 13, we compare the delay transform of visibilities before and after CLEANing at the LST of 4 hr. While CLEANing is able to remove structure within the horizon, it does not reduce any of the power leaked outside of the horizon by the chromaticity of the dish. If we knew the spectral response of the dish perfectly, then we might be able to CLEAN with this kernel and remove the supra-horizon structure though the true direction-dependent nature of the chromaticity would complicate this task.

Zoom In Zoom Out Reset image size

To form estimates of the 21 cm power spectrum, we split each visibility into Blackman–Harris windowed subbands centered at redshift intervals of Δz = 0.5 and each with a noise equivalent bandwidth of 10 MHz.¹⁷ The flat sky approximation allows us to Fourier transform each interval in frequency, square, and multiply by a set of prefactors to obtain a power spectrum estimate (Parsons et al. 2012a, 2014),

$$\begin{eqnarray}&&\widehat{P}({\boldsymbol{k}})={\left(\displaystyle \frac{2{k}_{{\rm{B}}}}{{\lambda }^{2}}\right)}^{2}\displaystyle \frac{{X}^{2}Y}{{B}_{{pp}}{{\rm{\Omega }}}_{{pp}}}| \widetilde{V}({\boldsymbol{u}}){| }^{2}.\end{eqnarray} \tag{ 23 }$$

Here, λ is the central wavelength of the observation, k_B is the Boltzmann constant, B_pp is the integral of the square of the bandpass, and Ω_pp is the integral of the primary beam squared over solid angle. X and Y are linear factors converting between native interferometry and cosmological coordinates, defined through the relation $2\pi {\boldsymbol{u}}=2\pi (u,v,\eta )=({{Xk}}_{x},{{Xk}}_{y},{{Yk}}_{z})$.

As a drift scan instrument, HERA will observe the sky at many LSTs within the declination stripe that passes through its primary beam, averaging over the power spectrum estimate at each LST. It is well documented that foreground power varies significantly over LST (Thyagarajan et al. 2015b), hence such an estimate will either filter or weight LSTs in a way that minimizes the impact of the most contaminated observations. For our analysis, we focus on a single LST of 4 hr, which is a patch of sky with low Galactic contamination. Such a patch is representative of the kind that HERA will focus most of its observing time on.

Computing the power spectra, we inspect the amplitude of foregrounds given the delay-response and angular pattern of the HERA dish for baselines of two different lengths in Figure 14. In both baselines, we find that the residuals after CLEANing tend to be at similar levels except at the subband centered at z = 8.5 (150 MHz), where the cylinder resonance is present. In bands outside of the resonance, foreground power outside of the horizon is dominated by sidelobes of the finite window function used in the Fourier transform. Hence, if we remove power inside of the horizon through cleaning, the level of this leakage is reduced. Power that is not reduced by the 1D clean is caused by intrinsic super-horizon structure in the antenna gain. Over the resonance, foreground residuals remain above the signal level out to k_∥ = 0.23 h Mpc⁻¹. Several other baselines, especially those oriented entirely in the E–W direction, which we do not show, contain up to two orders of magnitude greater galactic contamination as is noted in Thyagarajan et al. (2015a). In the data analysis, these baselines would be down-weighted or discarded so that they do not bias our final estimate. Outside of the resonance, beam chromaticity is dominated by beam-antenna reflections and since the level of foreground contamination is nearly identical to the achromatic beam, we conclude that feed-dish reflections do not pose a significant limitation on HERA's ability to measure the EoR power spectrum. Dishes are therefore a viable strategy for scaling the collecting area of 21 cm interferometry experiments at low cost.

Zoom In Zoom Out Reset image size

The level of the foregrounds in Figure 14 is conservative since no attempt has been made to apply inverse covariance weighting techniques (Tegmark 1997; Liu & Tegmark 2011; Dillon et al. 2013, 2015a, 2015b; Liu et al. 2014a, 2014b; Parsons et al. 2014; Trott et al. 2016) or fringe rate filtering (Parsons et al. 2015). Applications of inverse covariance filters in recent PAPER observations have yielded reductions in foreground power by greater than an order of magnitude (Ali et al. 2015); hence, our simulations show that even with the presence of the enhanced spectral structure from reflections, HERA will be able to isolate foregrounds well below the level of thermal noise in most of the EoR window.

A primary near-term goal of 21 cm EoR observations is to obtain information about the nature of the sources that drove reionization. Since the amplitude of the 21 cm signal is largest at smaller k values, a loss of large-scale signal due to foreground contamination eliminates the modes on which HERA would otherwise have the greatest signal-to-noise detections, reducing its overall sensitivity. In this section, we use the Fisher Matrix formalism to estimate the impact of HERA's intra-dish reflections on its sensitivity and ability to constrain the astrophysics of reionization. The Fisher Matrix allows us to forecast the covariances and errors on reionization parameters given errors on power spectrum observations due to the uncertainties caused by cosmic variance and thermal noise, which is in turn determined by the uv coverage and observing time of the interferometer. The covariance between the parameters, ${\boldsymbol{\theta }}$, of an astrophysical model is given by the inverse of the Fisher matrix, ${\boldsymbol{F}}$, which for Gaussian and independently determined power spectrum bins may be written approximately as (Pober et al. 2014),

$$\begin{eqnarray}&&{F}_{{ij}}\approx \displaystyle \sum _{k,z}\displaystyle \frac{1}{{\sigma }^{2}(k,z)}\displaystyle \frac{\partial {{\rm{\Delta }}}^{2}(k,z)}{\partial {\theta }_{i}}\displaystyle \frac{\partial {{\rm{\Delta }}}^{2}(k,z)}{\partial {\theta }_{j}},\end{eqnarray} \tag{ 24 }$$

where ${{\rm{\Delta }}}^{2}(k,z)$ is the power spectrum amplitude for some k − z bin and σ²(k, z) is the variance of the power spectrum estimate. We write this equation as approximate since it ignores additive terms arising from the dependence of σ²(k, z) on model parameters due to cosmic variance, which only contribute at the ≈1% level (Ewall-Wice et al. 2016b).

To simulate Δ², we use the publicly available 21cmFAST ¹⁸ code (Mesinger et al. 2011), which generates realizations of the 21 cm brightness temperature field using the excursion set formalism of Furlanetto et al. (2004). We employ a popular three parameter model of reionization (Mesinger et al. 2012) with the following variables.

• 1.  
  ${\boldsymbol{\zeta }}$: the "ionization efficiency" is defined in Furlanetto et al. (2004) to be the inverse of the mass collapse fraction necessary to ionize a region and is computed from a number of other physical parameters including the fraction of collapsed baryons that form stars and the UV photon escape fraction. Because ζ acts as an efficiency parameter, its primary effect is to change the timing of reionization. We choose a fiducial value of ζ = 20, though expected values range anywhere between 5 and 50 (Songaila & Cowie 2010).
• 2.  
  ${{\boldsymbol{R}}}_{\mathrm{mfp}}$: the presence of Lyman limit systems and other potential absorbers within H ii regions causes UV photons to have a finite mean-free path denoted by R_mfp. In the 21cmFAST framework, H ii regions cease to grow after reaching the radius of R_mfp, primarily impacting the morphology of the signal. We choose a fiducial value of R_mfp = 15 Mpc, which is in line with recent simulations accounting for the subgrid physics of absorption (Sobacchi & Mesinger 2014).
• 3.  
  ${{\boldsymbol{T}}}_{\mathrm{vir}}^{\min }$: the minimal mass of dark matter halos that hosted ionizing sources. While in principle, halos with virial temperatures as small as 10² K are thought to be able to form stars (Haiman et al. 1996; Tegmark et al. 1997), thermal and mechanical feedback have been seen to raise this limit to as high as 10⁵ K (Springel & Hernquist 2003; Mesinger & Dijkstra 2008; Okamoto et al. 2008). We choose a fiducial value of ${T}_{\mathrm{vir}}^{\min }=1.5\times {10}^{4}$ K, which is set by the atomic line cooling threshold.

In order to account for the degeneracies in the power spectrum between heating from X-rays and reionization from UV photons, we also marginalize over three additional parameters that describe the impact of heating from early X-ray sources as explored in Ewall-Wice et al. (2016b). These are the X-ray heating efficiency, f_X; the maximal energy of X-ray photons that are self absorbed by the ISM of early galaxies, ν_min; and the spectral slope, α, which are taken to have fiducial values of 1, 0.3 keV, and −1.2 respectively. We choose to parameterize our model in terms of the fractional difference of each variable from its fiducial value so that, for example, ${\theta }_{\zeta }=(\zeta -{\zeta }_{\mathrm{fid}})/{\zeta }_{\mathrm{fid}}$ and compute the derivatives in Equation (24) by performing a linear fit to realizations of the 21 cm power spectrum calculated by 21cmFAST at ${\theta }_{i}=\pm {10}^{-2},\pm 5\times {10}^{-2},\pm {10}^{-1}$, and ±2 × 10⁻².

For each measurement in the uv plane, the standard deviation of a power spectrum measurement (σ²(k, z)) is given by the direct sum of sample variance and thermal noise (McQuinn et al. 2006) which in turn depends on the primary beam of the instrument and the time spent sampling each uv cell. For our analysis, we assume that the uv plane is sampled by circular apertures with effective areas of 98 m² and that $\tau ({\boldsymbol{k}})$ is determined by a drift scan in which non-instantaneously redundant baselines are combined within each uv cell. We compute the standard deviation of each (k, z) bin for the proposed 350-element deployment of HERA (HERA-350) using the public 21cmSense ¹⁹ code (Pober et al. 2013b, 2014).

We have seen in Figure 14 that the simulated chromaticity of the dish leaks foregrounds beyond the horizon to varying degrees depending on the subband with the resonance occuring at z = 8.5(150 MHz). In this subband, foregrounds exceed the signal out to ∼380 ns, which for a 14.6 m baseline is ≈330 ns beyond the horizon. At all other redshifts, the leakage due to reflections only extends to ≈250 ns beyond the horizon. To determine the impact of the observed reflections on HERA's ability to constrain the astrophysics of reionization, we vary the cutoff in delay below which cosmological Fourier modes will be unobservable. To determine this cutoff, we consider three different scenarios for beam chromaticity that capture a range of possibilities informed by our simulations.

• 1.  
  Optimistic: the cylinder resonance is mitigated. In our most optimistic scenario, we assume that a more advanced feed design is able to mitigate the 150 MHz resonance in the cylinder while preserving polarization symmetry and effective area. In this case, the foregrounds pass below the level of the signal at 250 ns beyond the horizon at all redshifts between z = 7 and 12, consistent with what is observed in the bands where the reflections are less severe or when they are not present at all. Since the achromatic beam stays above the foregrounds below ≈250 ns as well, this number is set by the width of the Blackman–Harris window function and more optimal estimators of the power spectrum may reduce it substantially (see Ali et al. 2015). Hence, even our optimistic scenario is on the conservative side, making it more pessimistic than previous Fisher analyses (Pober et al. 2014; Ewall-Wice et al. 2016b; Liu & Parsons 2016; Liu et al. 2016). In this scenario, we assume that only modes within 250 ns of the horizon are unobservable.
• 2.  
  Moderate: the simulations accurately capture the chromaticity of the dish. In this case, we assume that the impedance match is sufficient for the contamination from reflections to pass below the level of the signal at 250 ns beyond the horizon except at redshift 8.5, where the feed resonance causes foregrounds to extend above the signal to 350 ns beyond the horizon. In this scenario, we assume that modes within 250 ns of the horizon are unobservable except for redshift 8.5, where modes below 350 ns beyond the horizon are unobservable. This is the scenario that would follow directly from our simulation results.
• 3.  
  Pessimistic: structure is present in all subbands. In this scenario, we assume that the resonance is still present in the neighborhood of 150 MHz and that a very poor impedance match (similar to the 500 Ω case in Figure 11) is obtained, allowing for significant contamination from reflections over the entire band, beyond the resonance. For this case, structure would be present out to ≈350 ns over the entire 100–200 MHz frequency range covered by HERA. In this scenario, we assume that modes within 350 ns of the horizon are unobservable across all redshift intervals.

In Figure 15, we compare the level of 1σ thermal noise for our optimistic and pessimistic scenarios to the amplitude of the 21 cm signal at several different redshifts. The smallest k modes lost to foreground contamination are the highest signal-to-noise measurements that HERA is expected to obtain. Their absence impacts sensitivity in two ways. By leading to a reduction in the maximal signal to noise ratio of ≈1.5 and a reduction in the total number of modes that the instrument is able to measure.

Zoom In Zoom Out Reset image size

Folding our calculations of thermal noise and the derivatives of Δ² into Equation (24) and inverting, we obtain the covariance matrix for model parameters. We show the 95% confidence ellipses for the reionization parameters in Figure 16. The presence of the resonance within a limited subband about z = 8.5 leads to an almost neglible increase in the extent of our confidence intervals while its presence across the entire band causes the lengths and widths of our confidence ellipses to increase by a factor of ≈2. The diagonal elements of our covariance give us error bars on each parameter, which we plot in Figure 17. We see that, similar to our confidence regions, the error bars on reionization parameters for our optimistic and moderate scenarios are nearly indestinguishable. In our worst case, where the resonance is present across the entire band, we see an increase in our error bars by a factor of ≈2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The error bars on reionization parameters, even for our most optimistic model, are a factor of a few larger than the "moderate" errors predicted in previous works using likelihood analyses such as Pober et al. (2014), Greig & Mesinger (2015), Greig et al. (2015a, 2015b) Ewall-Wice et al. (2016b), Liu et al. (2016), and Liu & Parsons (2016). There are several reasons for this. First, in our most optimistic scenario, we assumed that foregrounds cause the signal to be inaccessible below 250 ns beyond the edge of the wedge, which corresponds to k_min ≈ 0.15 h Mpc⁻¹ at z = 8.5, while in previous works a minimal comoving k_min ≈ 0.1 h Mpc⁻¹ was used. Second, previous studies assumed a fully illuminated HERA aperture, which for a 14 m diameter dish predicts an effective area of ≈155 m². Electromagnetic simulations and ORBCOMM mapping of the angular beam pattern of the HERA dish show that the effective area of the antenna element is actually ≈98 m² at 137 MHz, which leads to an increase in the overall thermal noise levels by a factor of ≈1.5 (Neben et al. 2016).

Although conservative, our analysis shows that the level of the HERA antenna's delay response is not an insurmountable obstacle for the 21 cm power spectrum measurement. Even if the additional systematics introduced by the feed resonance are present over the entire band, HERA will obtain a ≳10σ detection of the power spectrum and be capable of establishing precision constraints on the properties of the sources that reionized the IGM.