The following article is Free article

Ionospheric Attenuation of Polarized Foregrounds in 21 cm Epoch of Reionization Measurements: A Demonstration for the HERA Experiment

Zachary E. Martinot, James E. Aguirre, Saul A. Kohn, and Immanuel Q. Washington

Published 2018 December 13 • © 2018. The American Astronomical Society. All rights reserved.
The Astrophysical Journal, Volume 869, Number 1Citation Zachary E. Martinot et al 2018 ApJ 869 79DOI 10.3847/1538-4357/aaeac6

PDF Opens in a new tab.

ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

Article metrics

1043 Total downloads
0 Video abstract views

Permissions

Get permission to re-use this article

Dates

Received 2018 August 11
Revised 2018 October 18
Accepted 2018 October 21
Published 2018 December 13

Keywords

atmospheric effects; cosmology: observations; dark ages, reionization, first stars; polarization; techniques: interferometric

Journal RSS

Create or edit your corridor alerts

What are corridors?

Abstract

Foregrounds with polarization states that are not smooth functions of frequency present a challenge to H i Epoch of Reionization (EOR) power spectrum measurements if they are not cleanly separated from the desired Stokes I signal. The intrinsic polarization impurity of an antenna’s electromagnetic response limits the degree to which components of the polarization state on the sky can be separated from one another, leading to the possibility that this frequency structure could be confused for H i emission. We investigate the potential of Faraday rotation by Earth’s ionosphere to provide a mechanism for both mitigation of and systematic tests for this contamination. Specifically, we consider the delay power spectrum estimator, which relies on the expectation that foregrounds will be separated from the cosmological signal by a clearly demarcated boundary in Fourier space and is being used by the Hydrogen Epoch of Reionization Array (HERA) experiment. Through simulations of visibility measurements that include the ionospheric Faraday rotation calculated from real historical ionospheric plasma density data, we find that the incoherent averaging of the polarization state over repeated observations of the sky may attenuate polarization leakage in the power spectrum by a factor of 10 or more. Additionally, this effect provides a way to test for the presence of polarized foreground contamination in the EOR power spectrum estimate.

Export citation and abstractBibTeX RIS

Previous article in issue

Next article in issue

1. Introduction

Experiments seeking to observe the redshifted H i signal from the Epoch of Reionization (EOR) must contend with foregrounds that are ∼10⁴ times brighter than the cosmological signal by employing foreground removal or avoidance strategies (e.g., Santos et al. 2005; Bernardi et al. 2009, 2010; Pober et al. 2013; Dillon et al. 2014). These techniques rely on the smooth frequency structure of the foreground emission, in contrast to the spectrally structured cosmological signal (e.g., Datta et al. 2010; Morales et al. 2012; Trott et al. 2012; Pober et al. 2014; Liu et al. 2014a, 2014b; Thyagarajan et al. 2015a, 2015b). While the total intensity (Stokes I) of foreground radiation is spectrally smooth, Faraday rotation during propagation through our galaxy produces frequency structure in the linear polarization state (Stokes Q and U) at low frequencies (Jelić et al. 2010). Although extragalactic point sources appear largely depolarized, the large-scale synchrotron emission within the Milky Way appears to retain a significant level of polarization by the time it reaches an observer on Earth (Bernardi et al. 2013; Lenc et al. 2016).

The cosmological signal is expected to be effectively unpolarized given current experimental sensitivities (Babich & Loeb 2005; Hirata et al. 2018) and thus will be detected by measurements of Stokes I on the sky. On its own, frequency structure in the polarization state would not seem to be a concern when the objective is a measurement of the total intensity. However, the dipole antenna elements used in low radio frequency interferometers generally have significant sensitivity over the full sky when compared to the faintness of the EOR emission—bright foreground emission away from the antenna's boresight may still be relatively bright compared to the cosmological emission along the antenna’s boresight. Additionally, these dipole antennae are necessarily imperfect polarimeters over the full sky and do not naturally produce measurements of the incident radiation field in an orthogonal basis—a necessary condition to properly measure Stokes I. This imperfection in the measurements, commonly referred to as “polarization leakage,” means that even though we would like to obtain a pure measurement of the Stokes I intensity field on the sky, the measured visibilities will always involve a coupling to the polarization state of incident radiation. Although the sky is thought to be largely depolarized in the low-frequency radio spectrum (Farnes et al. 2014), even a polarization fraction of p ≈ 10⁻², which implies a polarized brightness that is “small” compared to total intensity, is not necessarily negligible compared to the cosmological signal, and therefore has the potential to produce contamination that is comparable to the EOR signal. This coupling must be understood and appropriately addressed to ensure that frequency structure in the polarization state of astrophysical foregrounds will not be mistaken for the cosmological power spectrum.

As a successor to the PAPER experiment (Parsons et al. 2010), the Hydrogen Epoch of Reionization Array (HERA) experiment (DeBoer et al. 2017) plans to use a delay-spectrum-based estimator (Parsons et al. 2012) to make measurements of the EOR power spectrum. In contrast to other efforts to observe the EOR that pursue imaging-based methods, the delay spectrum analysis approach does not involve precision imaging and thus has not included detailed modeling and subtraction of polarized foregrounds. This makes potential contamination due to polarization leakage particularly concerning for the HERA experiment. However, in Moore et al. (2017) it was proposed that the natural density fluctuations of the plasma in Earth’s ionosphere will produce a kind of polarization filter that can attenuate the coupling of visibility measurements to the polarization state of the sky.

In this paper we seek to understand the magnitude of this ionospheric attenuation effect in visibility measurements and the derived power spectra. We simulate interferometric visibilities based on models that include the wide-field effect of the ionosphere on the polarization state of diffuse foregrounds and the full-polarization instrumental response of an early HERA antenna design. This paper is organized as follows: In Section 2, we review the relevant mathematical description of polarization in interferometric measurements including ionospheric Faraday rotation and present a pedagogical picture of its attenuating effect on the measured polarized power. In Section 3, we discuss our implementation of this formalism, which involves modeling of the instrumental response, the diffuse polarized foreground emission on the sky, and calculations using archival data of Faraday rotations based on real ionospheric behavior. Section 4 presents the results of our simulations and analysis of the effect of ionospheric behavior on HERA observations. We conclude in Section 5.

2. Preliminary Formalism

2.1. Ionospheric Variation and Polarization Attenuation

The ionosphere is a turbulent upper region of Earth’s atmosphere that is ionized by solar radiation (e.g., Kintner & Seyler 1985; Loi et al. 2015). The permeation of this ionized medium by the persistent magnetic field of Earth then produces a magnetized plasma that will induce a rotation in the linear polarization state of electromagnetic plane waves propagating through it—the effect known as Faraday rotation. The rotation angle of the electric vector is φλ², where $\lambda =c/\nu$ is the wavelength and φ is the rotation measure (RM), which is given—in SI units—by Thompson et al. (2017):

where we have written the position vector ${\boldsymbol{s}}=s\hat{{\boldsymbol{s}}}$ and φ has units of rad m⁻². Here the integral is taken along the line of sight $\hat{{\boldsymbol{s}}}$ , the function ${\rho }_{e}(\hat{{\boldsymbol{s}}},s)$ is the free electron density at a radial distance s through the ionosphere, and ${\boldsymbol{B}}(\hat{{\boldsymbol{s}}},s)$ is the geomagnetic field.

The primary effect of the ionosphere that has concerned EOR power spectrum measurements so far is the refractive effect of the ionosphere (Vedantham & Koopmans 2015a, 2015b), and the fact that the ionosphere is a source of thermal radio emission has also been considered in the context of 21 cm global signal measurements (Sokolowski et al. 2015). Here we are instead concerned with the effect of the ionosphere on the polarization state of the sky and with the short baselines (approximately tens of wavelengths) of the compact HERA array that are most sensitive to the large-scale cosmological signal. The polarized emission at low radio frequencies appears to be dominated by large-scale diffuse Galactic emission rather than many unresolved point sources, and the effect of small refractive shifts on such spatially smooth emission is thus expected to be negligible. We focus here on the changing Faraday rotation over months-long timescales due to variations in the ionospheric free electron density.

Driven by the heating from the Sun, the free electron density ρ_e varies quasi-cyclically with the rotation of Earth at any fixed geographic location—the plasma density increases when the Sun is up and decreases at night, but the ionosphere will not return to exactly the same state. This means that observations of a polarized source on the sky made on different days will always be made through an ionospheric screen that is at least slightly different than on the previous day.

Moore et al. (2017) proposed that the effect of the ionospheric Faraday rotation on polarization leakage in a visibility could be estimated by approximating the RM over the sky as a constant $\varphi (\hat{{\boldsymbol{s}}})\approx \bar{\varphi }$ and additionally that the level of polarization leakage attenuation could be estimated without regard for the details of the instrumental response. While this simple approximation turns out to be an inadequate description of real visibilities, it is equivalent to considering the effect of the ionosphere for a single source on the sky and is a good way to build some intuition. This will be useful for interpreting the results of the detailed simulations in Section 4.

Suppose you used a good polarimeter to repeatedly observe a polarized source with a linear polarization state (Q, U) on each of N different days. Propagating through the ionosphere on the nth day will rotate the polarization state by an angle $2{\varphi }_{n}{\lambda }^{2}$ so that the observed polarization state is

where φ_n ∈ {φ₁, …, φ_N} is a sequence of different ionospheric RMs toward the source on the nth day. If one then averages over all these observations, the resulting quantity would be

While this would not be a sensible thing to do if one were actually interested in measurements of the polarization state, this averaging process is realized in the standard processing of HERA visibility data for power spectrum estimation. It is straightforward to see that this decreases the magnitude of the polarization $\overline{L}=| \overline{Q}+i\overline{U}|$ since the magnitude of the sum in the second line is always ≤N. The ratio of the intrinsic polarized power L² = Q² + U² to the power of the incoherently averaged polarization state ${\overline{L}}^{2}={\overline{Q}}^{2}+{\overline{U}}^{2}$ is then

We can think of the varying RM as defining a set of steps in a 2D plane with unit-length steps where the position after the Nth step is

Then, the attenuation is the (squared) ratio of the actual distance traveled from the origin $| Z(N)| $ to the maximum distance N that could have been traveled,

Examples of walks for several distributions of the φ_n values are shown in Figure 1 along with the resulting attenuation curves as a function of N. From the top, the first panel shows the walk when φ_n is simply a linear function

for some slope α. In this case the expression for the attenuation can be simplified by summing the geometric series

The attenuation factor is then

and the top right panel plots A² with α = 0.02 rad m⁻².

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Examples of four different types of normalized walks Z(N)/N in the unit disk described in the text. The left panels show the path of the walk over 100 steps, while the right panels show the (squared) distance from the origin at each step—i.e., the attenuation of the polarization amplitude for a source observed through an ionospheric variation characterized by the walk in the left panel. From top to bottom the walks are generated by a linear phase angle, an uncorrelated Gaussian distributed sequence of angles, a correlated Gaussian distributed sequence, and finally the same correlated Gaussian sequence with an additional linear trend added. In the middle column the solid circle is the circle of radius 1, while the dashed circle denotes the final radius of the walk at the position Z(100)/100, which is marked by a red dot.
Download figure:
Standard image High-resolution image

The second panel from the top shows the random walk generated by a realization of a sequence of N uncorrelated Gaussian-random variables ${\varphi }_{n}\sim { \mathcal N }(0,{\sigma }^{2})$ with σ = 0.2 rad m⁻². In this case it is straightforward to compute the expectation value of the attenuation factor in Equation 4(b), which is

While the real ionospheric RM sequences of interest to us are not necessarily well described as a Gaussian-random variable or a purely linear trend, the features of these simple models are worth noting. In the case of the linear trend, the attenuation oscillates as the walk passes near the origin, but it is bounded by a $\sim \tfrac{1}{{N}^{2}}$ envelope. For the Gaussian distribution with finite variance σ², the attenuation eventually approaches an asymptote ${A}^{2}\to {e}^{-4{\sigma }^{2}{\lambda }^{4}}$ as $N\to \infty$ . On the other hand, taking ${\sigma }^{2}\to \infty$ produces the limit $\langle A\rangle \to \tfrac{1}{N}$ , corresponding to a uniform distribution over the angle $2{\varphi }_{n}{\lambda }^{2}$ .

The third panel from the top then shows the walk generated by a sequence of correlated Gaussian-random variables $({\varphi }_{1},\,\ldots ,\,{\varphi }_{N})\in { \mathcal N }(0,{\boldsymbol{\Sigma }})$ , where the covariance matrix ${\boldsymbol{\Sigma }}$ is given by

Then, the panel at the bottom shows the result of adding a linear trend to the exact same φ_n sequence as in the third panel, ${\varphi }_{n}\to {\varphi }_{n}+\alpha n$ with α = 0.007 rad m⁻².

These walks and the resulting attenuation curves provide some intuition for the features of the attenuation curves obtained from the more complex visibility simulations—varying degrees of smoothness, discontinuities, oscillation, and lack of consistent monotonicity in frequency. The attenuation of polarization in a visibility can then be seen as a function of each of the random walks taken by the polarization in each direction on the sky—and this function is explicitly the visibility measurement equation.

2.2. Polarization in Interferometric Visibilities and the Delay Spectrum

The fundamental measurement made by interferometric arrays is the correlation function of the electric field. The van Cittert–Zernike theorem (Born & Wolf 1999), suitably generalized to include the four possible two-point correlation functions for pairs of isolated and identical dual-feed dipole antennae (Carozzi & Woan 2009; Smirnov 2011a, 2011b), relates the polarized intensity distribution on the sky to a measured visibility matrix by

where the integral is taken over the unit sphere ${{\mathbb{S}}}^{2}=\{\hat{{\boldsymbol{s}}}:\hat{{\boldsymbol{s}}}\in {{\mathbb{R}}}^{3},\parallel \hat{{\boldsymbol{s}}}\parallel =1\}$ (aka “the sky”). Here the vector ${\boldsymbol{b}}$ denotes a baseline between the two antennas, ν is the sampled frequency, ${\boldsymbol{ \mathcal J }}$ is a direction-dependent Jones matrix, a and b label two different antenna feed orientations, and ${\boldsymbol{ \mathcal C }}$ is the polarized brightness field on the sky expressed as the rank-2 coherency tensor field

The symbol ⨂ denotes a tensor product of the unit vectors, and the $\langle .\rangle$ denotes an ensemble average of the incoherent celestial radiation fields. The functions ${{ \mathcal E }}_{\delta },{{ \mathcal E }}_{\alpha }$ are the projections of the complex-valued electric vector amplitude

of a plane wave with wavevector $\propto \hat{{\boldsymbol{s}}}$ and components specified in the normalized tangent basis $\{{\hat{{\boldsymbol{e}}}}_{\delta },{\hat{{\boldsymbol{e}}}}_{\alpha }\}$ induced by the equatorial coordinates R.A. α and decl. δ. As a Hermitian matrix the coherency matrix is by definition specified by the frequency- and direction-dependent Stokes parameters I, Q, U, V so that

where the ${{\boldsymbol{\sigma }}}_{{ \mathcal S }}$ matrices are the Pauli matrices

Equation or symbol description not available

For our purposes here it is useful to adopt the point of view of a fixed observer under a rotating sky. Therefore, we will think of ${\boldsymbol{ \mathcal C }}$ as a time t dependent function, while the instrumental response of a drift-scanning antenna is fixed with respect to t. The coherency tensor is then a periodic function of the time t of the observation with a period T, which is the rotational period of Earth,

On the other hand, the ionospheric RM is only quasi-cyclic and thus not periodic in t. Since visibility data are averaged over multiple days of observation at the same local sidereal time (LST), it is useful to break the time variable into the LST t ∈ [0, T) and an integer n that indexes sidereal days. So from here on we will use the t variable to refer only to the LST of an observation. Then we write ${\boldsymbol{ \mathcal V }}(n,t,\nu ,{\boldsymbol{b}})$ as the visibility matrix observed on the nth sidereal day at the LST t.

The effect of the ionospheric Faraday rotation is described by a Jones matrix

which describes the rotation ${\boldsymbol{ \mathcal E }}\to {{\boldsymbol{R}}}_{n}{\boldsymbol{ \mathcal E }}$ of the polarization vector ${\boldsymbol{ \mathcal E }}$ of a plane wave upon propagation through the ionosphere. In this work we compute the visibilities resulting from

where ${\boldsymbol{J}}$ is the instrumental Jones matrix that describes the response of the instrument to polarized plane wave excitations. We do not include further direction-independent Jones matrices, so that our analysis concerns an idealized limit of data that has been calibrated for the direction-independent receiver-chain effects, nor do we include a thermal noise term, in order to isolate the effect of the ionospheric Faraday rotation.

From the visibilities we may form linear combinations analogous to the Stokes parameters, which we will refer to as “Vokes” parameters (e.g., “Vokes I parameter”). The Vokes parameters are

As noted in Section 1, the cosmological signal we are interested in detecting is thought to be effectively unpolarized, so Vokes I provides the highest sensitivity to the cosmological Stokes I signal, even though it is not generally a pure measurement of Stokes I. Thus, the quantity used for estimation of the power spectrum in the delay spectrum estimator is

where

are the instrumental Mueller matrix elements as shown in Figures 2 and 3, and Q_n, U_n are the linear polarization components after undergoing Faraday rotation in the ionosphere (Equation (2)).

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Mueller matrix elements for the $\{{\hat{{\boldsymbol{e}}}}_{\delta },{\hat{{\boldsymbol{e}}}}_{\alpha }\}$ basis (see Appendix B) in Equation 19(c) derived from a simulation of a transmitting HERA antenna’s far-field electric vector fields. Each image is a Lambert equal-area projection of the function on the hemisphere centered on the antenna’s boresight—the inscribed circle is the antenna’s local horizon. The rows are different frequencies 110, 130, 150, 170, and 190 MHz, from top to bottom. Since the M_IQ, M_IU elements take values in a range that is symmetric about 0, the color scale is a symmetric log₁₀ that spans six orders of magnitude on each side, with linearized values in the range of 10⁻⁶ to 10⁻⁸. All values are <1 in absolute value, so the sign is unambiguous. Overlaid on each image in the last column on the right is a unit tensor field that shows the geometric orientation of the polarization leakage. The plotted vector is headless and thus symmetric under a rotation by an angle π. This reflects the symmetry of the polarization state Q + iU of the *source* of incident radiation under a rotation by π.
Download figure:
Standard image High-resolution image

**Figure 2.** Mueller matrix elements for the $\{{\hat{{\boldsymbol{e}}}}_{\delta },{\hat{{\boldsymbol{e}}}}_{\alpha }\}$ basis (see Appendix B) in Equation 19(c) derived from a simulation of a transmitting HERA antenna’s far-field electric vector fields. Each image is a Lambert equal-area projection of the function on the hemisphere centered on the antenna’s boresight—the inscribed circle is the antenna’s local horizon. The rows are different frequencies 110, 130, 150, 170, and 190 MHz, from top to bottom. Since the M_IQ, M_IU elements take values in a range that is symmetric about 0, the color scale is a symmetric log₁₀ that spans six orders of magnitude on each side, with linearized values in the range of 10⁻⁶ to 10⁻⁸. All values are <1 in absolute value, so the sign is unambiguous. Overlaid on each image in the last column on the right is a unit tensor field that shows the geometric orientation of the polarization leakage. The plotted vector is headless and thus symmetric under a rotation by an angle π. This reflects the symmetry of the polarization state Q + iU of the *source* of incident radiation under a rotation by π.
Download figure:
Standard image High-resolution image

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Mueller matrices at 150 MHz as defined in Equation (20) for the HERA model (top) and the analytically defined Airy beam dipole (bottom). As in Figure 2, the matrix elements are shown in the $\{{\hat{{\boldsymbol{e}}}}_{\delta },{\hat{{\boldsymbol{e}}}}_{\alpha }\}$ basis (see Appendix B). The rows are the kernels for each of the Vokes parameters, i.e., the first row of each matrix is M_II, M_IQ, M_IU, M_IV from left to right, etc. The color scales are as described in Figure 2, but note that the diagonal elements have a different range than the off-diagonal components.
Download figure:
Standard image High-resolution image

As an aside, it may also be useful to note that in the same way that the linear polarization on the sky may be described as a spin-2 field (in either Cartesian or polar form),

the associated Mueller matrix elements describing the polarization impurity are also the components of a spin-2 field,

The Vokes I polarization leakage terms can then be thought of as an inner product between the polarization impurity of the instrument and the polarization state Q, U of the incident radiation, which is

The right column of Figure 2 shows the scalar function M_IL overlaid with a unit tensor field that shows the orientation on the sky of the instrumental impurity—geometrically, when the unit tensors of the instrumental impurity and the polarization state of the sky are at 45° (or 135° measured the other way), the inner product is 0. When the angle is 90°, the inner product is negative and minimized, i.e., M_IQQ + M_IUU = −M_ILL. The visualization in Figure 2 permitted by this representation may be useful for understanding the effect of the ionospheric Faraday rotation, which will be examined in more detail in Section 4.

The power spectrum can then be estimated from ${{ \mathcal V }}_{I}$ for each baseline through the delay transform (Parsons et al. 2012)

where ${ \mathcal B }$ is the frequency band selected to estimate the power spectrum at a given redshift and $W(\nu ,{ \mathcal B })$ is a windowing function that accounts for the finite bandwidth and any further choice of tapering function. Then, an estimator $\widehat{P}(k)$ for the spherically averaged power spectrum P(k) is obtained by averaging over time samples of the delay spectra (see, e.g., Ali et al. 2015 for the additional analysis complexities required with real measurements, and Liu et al. 2016 for more in-depth theoretical considerations),

This method of estimating the power spectrum motivates the analysis in this paper, but since we will only be considering ratios of different power spectra, the precise proportionality is unimportant here.

3. Visibility Simulation Components

We compute the visibility matrix in Equation 12(a) by a quadrature on a HEALPix¹ pixelization of the sky (Górski et al. 2005). The functions ${\boldsymbol{J}}(\nu ,\hat{{\boldsymbol{s}}})$ , ${\boldsymbol{ \mathcal C }}(t,\nu ,\hat{{\boldsymbol{s}}})$ , and ${{\boldsymbol{R}}}_{n}(t,\nu ,\hat{{\boldsymbol{s}}})$ are evaluated for each $\hat{{\boldsymbol{s}}}={\hat{{\boldsymbol{s}}}}_{p}$ in the set of HEALPix pixels ${\{{\hat{{\boldsymbol{s}}}}_{p}\}}_{p=1}^{{N}_{p}}$ . Explicitly, Equation 12(a) is estimated as

where ${\rm{\Delta }}{\rm{\Omega }}=\tfrac{4\pi }{N}$ . In this section we discuss our definition and evaluation of these three functions, as well as particular parameters of our simulations. The visibility simulation code is publicly available.²

3.1. Calculation of the Ionospheric RM from Archival Total Electron Content

Over the past three decades, methods to measure the total electron content (TEC) using global positioning system (GPS) dual-frequency receivers have been developed and improved (e.g., Royden et al. 1984; Lanyi & Roth 1988; Mannucci et al. 1998; Schaer 1999. Commission géodésique 1999; Iijima et al. 1999; Komjathy et al. 2005; Erdogan et al. 2016). These methods utilize the TEC-induced time delay between the arrival of radio waves of two closely spaced frequencies to estimate the TEC value of the ionosphere above a GPS station. Repeating this method for stations around the world and interpolating spatially provides an estimate of the TEC above any location on Earth.

Meanwhile, many generations of the International Geomagnetic Reference Field (IGRF; e.g., Finlay et al. 2010) have continually improved the model of Earth’s magnetic field. This model is composed by spatial interpolation of magnetic field measurements (in up to 13th-order spherical harmonic coefficients) reported by institutions around the world.

Based on the IonFR³ package of Sotomayor-Beltran et al. (2013), we have developed radionopy,⁴ a python package to calculate ionospheric RM values (the function φ in Section 2). Like IonFR, radionopy uses GPS-derived TEC maps (in the IONosphere Map EXchange format; ionex) and the IGRF to estimate the value of φ at a given latitude, longitude, and date. Unlike its predecessor, radionopy is written to calculate φ( $\hat{s}$ ) over an arbitrary point set of directions on the sky, allowing images of the full sky. Additionally, radionopy implements the temporal interpolation scheme recommended in Schaer et al. (1998) to obtain full-sky maps for arbitrary times between the 2 hr time resolution of the provided IONEX data, such that the resulting φ(n, t, $\hat{{\boldsymbol{s}}}$ ) is a fairly smooth function of t. The interpolation scheme is as follows. Let η denote universal time (UT) and i index the times at which the TEC maps

are available as a function of the geocentric latitude θ and longitude ϕ. Then the interpolated TEC map at an arbitrary time η such that η_i ≤ η ≤ η_i+1 is a linear interpolation of the forward- and backward-rotated preceding and succeeding maps given by

where ϕ_k = ϕ + ω_⊕(η − η_k), with ω_⊕ the angular speed of Earth. The RM function φ(n, t(η), $\hat{{\boldsymbol{s}}}$ (θ, ϕ)) is then computed by the approximation of Equation (1) described in Sotomayor-Beltran et al. (2013). Figure 4 shows an example of radionopy output for the RM function φ(n, t, $\hat{{\boldsymbol{s}}}$ ) evaluated in altitude/azimuth coordinates at the location of the HERA array, and Figure 18 shows example output TEC maps over Earth.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Images of $\varphi (\hat{{\boldsymbol{s}}})$ at a fixed LST of 2.5 hr. The first two maps are consecutive sidereal days; the third is 30 days later. At the top of each panel is the civil time and date of the RM snapshot. The images are horizon to horizon in a Lambert equal-area projection centered on longitude = +214283, latitude = −307215. The square, diamond, and circle indicate the points for which RM sequences over the day index n are shown in Figure 5. The images are oriented so that north is up and east is to the right. The places where $\varphi \to 0$ correspond to points where ${\boldsymbol{B}}\cdot \hat{{\boldsymbol{s}}}=0$ ; the null to the north is near the equator.
Download figure:
Standard image High-resolution image

**Figure 4.** Images of $\varphi (\hat{{\boldsymbol{s}}})$ at a fixed LST of 2.5 hr. The first two maps are consecutive sidereal days; the third is 30 days later. At the top of each panel is the civil time and date of the RM snapshot. The images are horizon to horizon in a Lambert equal-area projection centered on longitude = +214283, latitude = −307215. The square, diamond, and circle indicate the points for which RM sequences over the day index n are shown in Figure 5. The images are oriented so that north is up and east is to the right. The places where $\varphi \to 0$ correspond to points where ${\boldsymbol{B}}\cdot \hat{{\boldsymbol{s}}}=0$ ; the null to the north is near the equator.
Download figure:
Standard image High-resolution image

While the intrinsic time and spatial resolution of the resulting RM maps is relatively low, the RM computed in this way has recently been validated as being reasonably close to more precise measurements using pulsar timing dispersion (Malins et al. 2018).

Figure 5 shows a selection of RM sequences over 100 sidereal days at a fixed LST hour of 2.5 for several different years. The RM due to the ionosphere is generally a random function over time, with the underlying random variable being the TEC whose variation is driven by solar radiation. The RM varies randomly from day to day but follows a clear trend over the course of 100 days.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Rotation measure sequences over 100 days at three different points on the sky at fixed LST for each year from 2011 to 2014. The error bars are the 1σ error bars propagated from TEC uncertainties provided with the IONEX data. There is a clear trend along with the random variation. A significant solar event is observable as the large jump in 2014 (note the scale on the vertical axis of the bottom row of panels differs from the top three rows). This appears to correspond to a relatively large solar flare that was observed on 2014 October 19 by NASA’s *Solar Dynamics Observatory*, which was followed by several weeks of abnormally high solar activity.
Download figure:
Standard image High-resolution image

The cause of this trend can be understood broadly by noting that the magnitude of the RM goes inversely as the time since the Sun went down. As the season progresses, a given LST transit occurs progressively closer to the previous sunset. Since the Sun is the driver of ionization in the atmosphere, as this proximity increases, the ionosphere has had less time for recombination to occur since it was last heated, resulting in a higher free electron density and thus a higher magnitude of RM.

Careful inspection of Figure 5 would reveal that each of the three different points in the RM maps is not exactly a rescaling of a common function of n, i.e., φ is not a separable function of n and $\hat{{\boldsymbol{s}}}$ . However, it is clear in Figure 4 that there is a distinct average shape to the function $\varphi (\hat{{\boldsymbol{s}}})$ , which is largely due to the increasing path length through the ionosphere with increasing zenith angle, as well as the projection of the geomagnetic field, along different lines of sight. These observations are quantified somewhat by considering the matrix

where

which measures the similarity in spatial structure between different days. The integral is taken over the observed hemisphere ${{\mathbb{S}}}_{+}^{2}$ , and the absolute value of φ is taken because the sign does not vary between days. We can also compute this correlation for the TEC by replacing φ(n, t, $\hat{{\boldsymbol{s}}}$ ) with ${\overline{\rho }}_{e}(n,t,\hat{{\boldsymbol{s}}})$ in Equations (30) and (31). A representative example of the matrices c_kl for both the RM and TEC are shown in Figure 6, where we see that the shape of the RM field is not as variable between days as the underlying TEC field.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Left: matrices as defined in Equation (29) at t = LST 2.5 hr for 100 days starting on 2011 September 10. Right: same function as in the left panel, but applied to the TEC column density, i.e., with $\varphi (n,t,\hat{{\boldsymbol{s}}})$ replaced by ${\overline{\rho }}_{e}(n,t,\hat{{\boldsymbol{s}}})$ .
Download figure:
Standard image High-resolution image

**Figure 6.** Left: matrices as defined in Equation (29) at t = LST 2.5 hr for 100 days starting on 2011 September 10. Right: same function as in the left panel, but applied to the TEC column density, i.e., with $\varphi (n,t,\hat{{\boldsymbol{s}}})$ replaced by ${\overline{\rho }}_{e}(n,t,\hat{{\boldsymbol{s}}})$ .
Download figure:
Standard image High-resolution image

3.2. Antenna Response Model

The instrumental Jones matrix ${\boldsymbol{J}}(\nu ,\hat{{\boldsymbol{s}}})$ is derived from electromagnetic simulations of a HERA antenna element using the commercial software CST, which solves for the far-field electric field radiated by an array antenna element when operated in transmission (Fagnoni & Acedo 2016). By Lorentz reciprocity these electric field functions define the polarized response of the instrument to the incident plane waves produced by celestial sources (De Hoop 1968; Born & Wolf 1999; Potton 2004; Balanis 2005) and thus define the instrumental Jones matrix as described in Appendix B.

We take the data for the fields output by CST and interpolate to a HEALPix map. Only the simulation for a single feed is available, so we use the assumption that the antenna structure is symmetric under 90° rotations about the antenna boresight to derive the response of the second feed. Since the simulation was computed at 1 MHz resolution in frequency (which was deemed sufficient to capture the frequency structure), an interpolation in frequency is performed by cubic spline fit to the components of the spherical harmonic transforms of the electric field components and then synthesizing the fields at the desired frequency and spatial resolution. The fields simulated by CST are defined over the full sphere. We apply a hard cut to the fields at the local horizon, which is defined as the zenith angle of π/2, and set the instrumental response ${\boldsymbol{J}}$ to zero below this horizon.

Additionally, for comparison with a simplified model we use a Hertzian dipole with an Airy disk directivity taper. In terms of a set of Cartesian coordinates (x, y, z) such that ${\hat{{\boldsymbol{e}}}}_{z}$ is along the antenna’s boresight and $\theta (\hat{{\boldsymbol{s}}})={\cos }^{-1}({\hat{{\boldsymbol{e}}}}_{z}\cdot \hat{{\boldsymbol{s}}})$ , this Jones matrix is

where J₁(x) is the Bessel function of the first kind of order one and 2a = 14.6 m is the diameter of a HERA dish. The purpose of this simpler and less realistic model is to illustrate how the details of the instrumental response affect the simulations. The HERA antenna simulation used here is that of an early model that is still in the process of development. We expect that the final model of the antenna far-field response will differ slightly from the one available to us now, but not significantly so. The difference will be much smaller than the difference that can be seen in Figure 3 between the current HERA model and this simple Airy dipole construction. Further, although the polarization properties of the antenna beam have not been measured, other, more accessible properties have been measured, and their agreement with the CST simulated model suggests that it is quite realistic. Therefore, despite the similarity with the HERA model, the Airy dipole should probably be thought of as a large change to the instrument model, rather than a small one.

3.3. Sky Model

While a good model of the actual diffuse Stokes I emission is available in the form of the Global Sky Model, there is currently no equivalent full-sky model of the polarization state of the diffuse galactic synchrotron emission. Therefore, the best method available is to use a random realization generated from a statistical model with constrained parameters. We use the CORA⁵ code to generate a set of Stokes I, Q, and U diffuse maps. The CORA package was developed for use in Shaw et al. (2015), and the details of its physical motivation and implementation are discussed there. Briefly, CORA generates random realizations of the polarization state of the diffuse emission by RM synthesis (Brentjens & De Bruyn 2005; Jelić et al. 2010)

where F describes the distribution and polarization angle (F is a complex-valued function) of polarized emission as a function of the Faraday depth ϕ. The function F is in turn decomposed into spherical harmonic components as

The components f_lm(ϕ) are Gaussian-random complex-valued functions, while $w(\phi ,\hat{{\boldsymbol{s}}})$ is a fixed function of ϕ and $\hat{{\boldsymbol{s}}}$ . Realizations of diffuse linear polarization components Q, U are thus generated by drawing realizations of the components f_lm(ϕ). Figure 7 shows an example of the diffuse polarized power $L(\hat{{\boldsymbol{s}}})=\sqrt{{Q}^{2}(\hat{{\boldsymbol{s}}})+{U}^{2}(\hat{{\boldsymbol{s}}})}$ and polarization orientation tensor field generated by this model.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Image of the polarized power $L=\sqrt{{Q}^{2}+{U}^{2}}$ in the sky model over half of the sky as would be observed by an antenna instantaneously, meaning that the edge of the image is the local horizon 90° from zenith. The linear color scale is normalized to the peak of the image, and both panels show the same map. Overlaid is a unit tensor field that shows the orientation of the linear polarization state (Q, U). The tensor field in the left panel shows the initial polarization orientation field, while the right panel shows the polarization orientation after ionospheric Faraday rotation at 150 MHz (i.e., the polarization orientation field of the polarization state (Q_n, U_n) in Equation 19(c)). Since the RM field $\varphi (\hat{{\boldsymbol{s}}})$ is spatially smooth, the orientation field after Faraday rotation maintains its initial spatial correlation.
Download figure:
Standard image High-resolution image

**Figure 7.** Image of the polarized power $L=\sqrt{{Q}^{2}+{U}^{2}}$ in the sky model over half of the sky as would be observed by an antenna instantaneously, meaning that the edge of the image is the local horizon 90° from zenith. The linear color scale is normalized to the peak of the image, and both panels show the same map. Overlaid is a unit tensor field that shows the orientation of the linear polarization state (Q, U). The tensor field in the left panel shows the initial polarization orientation field, while the right panel shows the polarization orientation after ionospheric Faraday rotation at 150 MHz (i.e., the polarization orientation field of the polarization state (Q_n, U_n) in Equation 19(c)). Since the RM field $\varphi (\hat{{\boldsymbol{s}}})$ is spatially smooth, the orientation field after Faraday rotation maintains its initial spatial correlation.
Download figure:
Standard image High-resolution image

Stokes I is generated by CORA as an extrapolation of the Haslam map at 408 MHz, but we subtract the Stokes I term from the Vokes parameters when analyzing the simulated visibilities in Section 4. The one exception is in Figure 8, where, for context, we show the simulated ${{ \mathcal V }}_{I}(t,\nu )$ function including the Stokes I term.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Amplitude and phase of ${{ \mathcal V }}_{I}(t,\nu )$ (left) and ${{ \mathcal V }}_{{IL}}(t,\nu )$ (right) for the intrinsic Vokes I visibility (φ = 0) computed in our fiducial visibility simulation. The amplitude plotted in both figures is relative to the maximum of $| {{ \mathcal V }}_{I}(t,\nu )|$ .
Download figure:
Standard image High-resolution image

**Figure 8.** Amplitude and phase of ${{ \mathcal V }}_{I}(t,\nu )$ (left) and ${{ \mathcal V }}_{{IL}}(t,\nu )$ (right) for the intrinsic Vokes I visibility (φ = 0) computed in our fiducial visibility simulation. The amplitude plotted in both figures is relative to the maximum of $| {{ \mathcal V }}_{I}(t,\nu )|$ .
Download figure:
Standard image High-resolution image

While the constraints on CORA's model parameters are, in the author’s words, “crude,” the model is sufficiently realistic to capture the important characteristic features of diffuse polarized emission, namely, the unsmooth frequency structure and spatial correlation. In particular, our results are by construction independent of the absolute level of polarized power present in the sky model and somewhat insensitive to the particular frequency spectrum of the polarization. We purposefully avoid speculation about the absolute level of polarization leakage that may be observed with HERA.

3.4. Simulation Parameters

1.
We compute the visibility matrix for a 30 m east–west-oriented baseline.
2.
The two dipole-feed orientations a and b are east–west and north–south, as is the case for the HERA antenna elements.
3.
The visibilities are computed for each of 201 equally spaced frequency points ν = ν_j (i.e., a 0.5 MHz channel width) in the band 100–200 MHz; the smallest frequency is 100 MHz, and the largest is 200 MHz. From this band five 20 MHz sub-bands are used:
4.
In the delay transform we use a Blackman–Harris window function.
5.
We compute visibilities using the historical ionospheric data for the 100-day sequence starting on September 10 in each year of interest; this is the sequence for which RM data are shown in Section 3.1. The LST hour range 1–4 was chosen so that all local times in this range are between sunset and sunrise for each of the 100 days at the geographic location of the HERA array.
6.
Fiducial simulations: We picked a single realization of the sky model to use for a set of fiducial simulations. In these simulations N_t = 96 equally spaced time samples t = t_l in the LST hour range 1–4 were computed. Since each of the functions modeled in our simulation is, by construction, smooth on the scale of our sky pixelization, this is sufficient to completely sample the time dependence of the visibilities. Visibilities were computed using the historical ionosphere data from the years 2009, 2011, 2012, and 2014 and for both instrumental response models.
7.
Sky model variance simulations: We also performed simulations using many realizations of the statistical sky model. In order to save computational time in these simulations, six equally spaced time samples were computed in the same LST range. For each of 12 yr from 2003 to 2014, visibilities for 100 different realizations of the sky model were computed, i.e., the 100 realizations are different for each year. The reduced cadence of the time sampling has an effect on the results, but we found from resamplings of the fiducial simulations that it was negligible compared to the change due to the sky model.

4. Results from Simulations

4.1. Attenuation of Linear Polarization in a Vokes Parameter Delay Spectrum

Power spectrum estimators based on the delay spectrum will generally average visibility measurements taken over multiple days at fixed (t, ν) in order to attenuate thermal noise. We follow this procedure by averaging the simulated visibilities over a set S_k of N sidereal days. The simulated visibilities ${\boldsymbol{ \mathcal V }}(n,t,\nu )$ are computed on an LST grid for each of N_d consecutive days indexed by $n\in S=\{1,2,\,\ldots ,\,{N}_{d}\}$ . We can then choose a subset ${S}_{k}\subset S$ and compute the average at fixed t as

and then the corresponding averaged Vokes parameters are

Here we find the averaging process for visibilities that was alluded to in Section 2. Now, instead of the polarization state of a single source, the averaged quantity is ${\overline{{ \mathcal V }}}_{{ \mathcal S }}$ , which can be considered a functional over the sequence of functions $\{\varphi (n)\}{}_{n\in {S}_{k}}$ . For the simple model in Section 2 of a single polarized point source the notion of attenuation of the polarized power was clear, but extending the idea to the delay spectrum of visibilities warrants some additional consideration.

For each individual day n the effect of the ionospheric Faraday rotation of the polarization state on the sky is a small change to the frequency spectrum ${\boldsymbol{ \mathcal V }}(\nu )$ . We can make the consequences more apparent by considering Equations (37) and (38) in greater detail. The instrumental response ${\boldsymbol{J}}$ is independent of n, so we take the sum in Equation (37) inside the integral defining ${\boldsymbol{ \mathcal V }}(\nu )$ in Equation 12(a):

The sum in parentheses may always be expressed as

since the cumulative effect of summing the N different rotations may be described by a single rotation by an angle 2μ where

and an amplitude factor

which define the resultant matrix

Note that ${{\boldsymbol{ \mathcal T }}}^{2}$ is the matrix representation of the complex number Z(N)/N (Equation (5)). Then, for each ${ \mathcal S }\in \{I,Q,U,V\}$ the averaged Vokes ${ \mathcal S }$ is

This form exposes the fact that the ionospheric Faraday rotation need not reduce the Vokes I polarization leakage; in fact, it can increase it. Suppose that M_IQQ + M_IUU = 0 for some ν and some $\hat{{\boldsymbol{s}}}$ so that the polarization leakage term is—by cosmic accident of alignment—intrinsically zero. Then, any rotation by a small angle $2\mu \ne 0$ will make the polarization leakage term nonzero. On the other hand, the rotation by itself can reduce the polarization leakage. Suppose now that ${M}_{{IQ}}Q+{M}_{{IU}}U\ne 0$ . Then, there is always a choice of rotation angle 2μ that will null the polarization leakage term given by

Of course, generally the change in the magnitude of the leakage terms at each point $\hat{{\boldsymbol{s}}}$ will fall between these two extremes, and the change in the visibility will be the result of integrating over all such changes. This may be visualized by comparing the polarization orientation of the model sky in Figure 7 and the instrumental response in Figure 2. The spatial coherence of the fields means that merely changing the polarization angle over the whole sky can have dramatic effects on the polarization leakage terms.

For small N where the ionosphere does not change very much between different days, the amplitude factor A is generally fairly close to unity, but the resultant effective rotation of the polarization state by the angle 2μ can produce a stronger (or weaker) instrumental coupling to the polarization state if the original state (Q, U) was not maximally (or minimally) aligned with the instrument. While A and $\cos (2\mu ),\sin (2\mu )$ are fairly smooth functions of frequency, the change in the frequency spectrum due to the realignment term (the term proportional to $\sin (2\mu )$ ) including nonsmooth Q, U will tend to make the power in any given τ mode of the leakage delay spectrum fluctuate slightly. As N increases and the variation between successive ionospheric Faraday screens becomes significant, the amplitude factor A decreases enough to attenuate the polarization leakage regardless of the relative orientation of the sky’s polarization state to the instrumental response. The result is that the precise attenuation may be somewhat variable as a function of N for different τ modes, but as N increases, it should tend to converge to an overall trend. For this reason, taking the ratio of a mean over modes of the delay spectrum of the polarization leakage provides a good summary measure of the attenuation.

With these considerations in mind, we define a metric to assess the overall level of attenuation of the polarization leakage. Since we are interested in how polarization will affect power spectrum measurements, we define the attenuation in terms of the delay spectra that would be used in such measurements.

The attenuation of polarization leakage due to the ionosphere is defined as the ratio of the total power in the leakage function after averaging visibilities over different ionospheric Faraday rotations to the total power in the intrinsic leakage that would occur if the observation was made in the absence of ionosphere rotation

where ${ \mathcal B }$ is the band over which the delay transform is computed. Then, the leakage delay power spectrum ${{ \mathcal L }}_{I}$ is computed as

and ${\widetilde{{ \mathcal V }}}_{{IL}}(n,t,\tau ,{ \mathcal B })$ is the delay transform (Equation (24)) of the leakage terms ${{ \mathcal V }}_{{IL}}(n,t,\nu )$ in the Vokes I visibility for the nth day:

The intrinsic leakage function ${{ \mathcal L }}_{I,\mathrm{int}}$ is formally obtained from Equation (47) by setting the ionospheric RM as φ = 0, so

is the Vokes I polarization leakage that would be observed without ionospheric interference and is constant as a function of n, which we think of as the intrinsic leakage.

Figures 8 shows the simulated ${{ \mathcal V }}_{I}(t,\nu )$ including the intrinsic leakage term, as well as ${{ \mathcal V }}_{{IL},\mathrm{int}}(t,\nu )$ on its own, and Figure 9 shows some examples of the effect of the ionospheric Faraday rotation on the frequency and delay spectra of the visibilities.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Left: LST-averaged square magnitude of the simulated Vokes I polarization leakage as a function of frequency ν, e.g., the result of averaging over t in the top right panel of Figure 8. Summing the black “Intrinsic” and pink “Averaged” curves over any of the sub-bands ${ \mathcal B }$ would produce the numerator and denominator for S_k = S₁₀₀ in Equation (46). Right: delay power spectra over the full 100–200 MHz band of the intrinsic Vokes I polarization leakage and a representative sample of the leakage when Faraday rotation for a single day is included. In both ν and τ representations there is a characteristic amplitude and shape, but the ionospheric Faraday rotation for a single day perturbs the spectra. The spectra resulting from averaging the visibilities over 100 days of different Faraday rotations produce an average attenuation, but additionally perturb the spectra mode by mode owing to the resultant effective rotation of the polarization state on the sky.
Download figure:
Standard image High-resolution image

**Figure 9.** Left: LST-averaged square magnitude of the simulated Vokes I polarization leakage as a function of frequency ν, e.g., the result of averaging over t in the top right panel of Figure 8. Summing the black “Intrinsic” and pink “Averaged” curves over any of the sub-bands ${ \mathcal B }$ would produce the numerator and denominator for S_k = S₁₀₀ in Equation (46). Right: delay power spectra over the full 100–200 MHz band of the intrinsic Vokes I polarization leakage and a representative sample of the leakage when Faraday rotation for a single day is included. In both ν and τ representations there is a characteristic amplitude and shape, but the ionospheric Faraday rotation for a single day perturbs the spectra. The spectra resulting from averaging the visibilities over 100 days of different Faraday rotations produce an average attenuation, but additionally perturb the spectra mode by mode owing to the resultant effective rotation of the polarization state on the sky.
Download figure:
Standard image High-resolution image

It is worth noting explicitly that by Parseval’s theorem the sums over the delay τ in Equation (46) are equivalent to simply summing the squared (and windowed) visibility amplitude over the frequency sub-band ${ \mathcal B }$ . We write the definition in terms of delay spectra in anticipation of modifications to the definition for use with real data or more realistic simulations in which we have more confidence in the detailed frequency–frequency covariance. In particular, with real data we cannot easily subtract the Stokes I term from our Vokes I data, but we could instead apply a high-pass delay filter by restricting the sum over τ-modes in Equation (46) to $| \tau | \gt | {\tau }_{\mathrm{filter}}|$ since $\widetilde{I}(\tau )$ is compact in delay compared to $\widetilde{Q}(\tau )$ and $\widetilde{U}(\tau )$ . More generally, we reiterate that the notion of “attenuation” is not uniquely defined, and we have made a particular choice here, though a well-motivated one, that is only applicable to simulated data.

In Figure 10 we plot the relative leakage factor ξ_I in each of five frequency bands over the 100 cumulative subsets S_k of S, the set of 100 different days for which visibilities were simulated. The cumulative subsets are ${S}_{1}=\{1\},\,\ldots \,{S}_{N}\,=\{1,2,\,\ldots ,\,N\},\,\ldots \,{S}_{100}=\{1,2,\,\ldots ,\,100\}$ , so in Figure 10 we can write the attenuation ξ_I as a function of N, the number of consecutive sidereal days in the sum. This will be referred to as the “natural attenuation,” since it is “natural” to simply average up all available data after observing for N days. Each panel shows the set of five such discrete attenuation curves for each of the four calendar years and two instrumental response models for which we computed visibilities as discussed in Section 3.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Natural attenuation of the Vokes I polarization leakage as defined in Equation (46) as a function of the number N of consecutive days averaged over, up to 100 days from September 10 for the four years 2009, 2011, 2012, and 2014. Left: from the fiducial simulations with the HERA antenna response model. Right: from the fiducial simulations with the simple Airy dipole model. It is interesting to see that the large change in the magnitude of the RM during the year 2014 (Figure 5) does not produce a correspondingly large change in the attenuation.
Download figure:
Standard image High-resolution image

The smooth nature of the attenuation curves is reminiscent of the correlated and trending walks in the plane considered in Section 2 and reflects the trends of φ(n) shown in Figure 5. We can also recognize the occasional oscillatory behavior as a natural feature of the polarization state undergoing a trending walk in the tangent plane at each point on the sky. The effect of the changing relative orientation of the polarization state to the instrumental response is evident as ξ_I > 1 for some of the curves when N is small.

The lack of uniform ordering by frequency band is due to the interplay between the intrinsic oscillation in the random walk of the polarization state and the frequency dependence of the instrumental coupling over each band ${ \mathcal B }$ . Although the instrumental response changes only slowly over each 20 MHz sub-band, we see that there is significant variation in Figure 2 over the full 100–200 MHz band. This contributes to the deviations from the ordering of attenuation curves by central frequency that we might naively expect from the simple formulae in Section 2. Additionally, there is the changing relative importance of the factor A and the μ-rotation as a function of N. For small N, the amplitude factor A ∼ 1, but for the highest-frequency band the rotation is changing the relative polarization angle between the sky and the instrumental response such that the Faraday-rotated leakage is smaller in magnitude than the intrinsic leakage. This is in contrast to the lower-frequency bands, some of which have the leakage amplified (ξ_I > 1) somewhat by the μ-rotation for small N. As N increases, the amplitude factor A becomes more significant. In particular, the magnitude of A decreases more quickly at lower frequencies. As a result, eventually the attenuation curves of the lower-frequency bands drop into a monotonic, or nearly monotonic, ordering with frequency band. Of course, μ is still a function of N, so its variation should be expected to contribute to eccentricities in the curves as the alignment of the effective polarization state to the instrumental response continues to gradually change.

The attenuation curves for the HERA model show that in 100 days the polarization leakage in our simulations is attenuated by a factor of 10 or more in the three bands in the 100–160 MHz range, and the two higher-frequency bands are not far behind. Comparison with the Hertzian dipole model shows that the details of the coupling between the instrument and the polarization state of the sky do affect the resulting visibility enough to noticeably change the attenuation factor as a function of N. This makes clear that accurate prediction of the attenuation factor is dependent on an accurate beam (and sky) model.

4.2. Sky Model Variance

There is significant variation in the natural attenuation curves with the changing ionospheric Faraday rotations for different years, which reflects the underlying variation of the visibilities. The fact that different rotations of the polarization state of the sky can produce such a change in the power spectrum suggests that if the Faraday rotation field and instrument model are held fixed, different polarized skies could also produce significantly different results.

The attenuation curves we have considered so far are the results of simulations using a single realization of a statistical model of the diffuse polarization. To understand how much our simulations could vary with the choice of sky model, we compute visibilities for 100 realizations of the diffuse polarization generated with CORA.

Figures 11–13 show the resulting natural attenuation curves of 100 different sky model realizations using the HERA instrument model over the same 100-day sequence in each year from 2003 to 2014. For each year and sub-band the geometric mean and geometric variance of the sample of attenuation curves are estimated, and the resulting mean curve and 2σ intervals are also shown.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Distributions of the Vokes I polarization leakage attenuation factor obtained in simulations with different realizations of the sky model. For each year (row) 100 different sky realizations are generated and the attenuation factor is computed in each sub-band (column). Each thin gray curve is the attenuation curve for a realization of the sky model. The thick black points show the geometric mean and geometric 2σ deviation of the distribution of ξ_I(N) for each N.
Download figure:
Standard image High-resolution image

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Same as in Figure 11, but for the years 2007–2010.
Download figure:
Standard image High-resolution image

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Same as in Figure 11, but for the years 2011–2014.
Download figure:
Standard image High-resolution image

There is a significant variance over the different realizations of the sky. Thus, it is not possible to predict with high accuracy the attenuation that occurs for a specific set of measurements without an accurate sky (and instrument) model. Nevertheless, for the purpose of forecasting the likely range of attenuation that might be obtained in HERA data, the variation is constrained enough that we still get a good idea of what to expect. By considering the mean attenuation curves for each year, we see a large variation as a function of the year. This is due to the solar cycle.

4.3. The Solar Cycle

Solar activity—the rate of ionizing flux and charged particle emission from the Sun—waxes and wanes with the ∼11 yr solar cycle and is highly correlated with the rate of sunspot occurrence. Thus, there is a correlation between the average TEC and solar activity (e.g., Sotomayor-Beltran et al. 2013), and we expect a similar correlation with the attenuation factor. Historical sunspot data, as well as future projections, are available from the NOAA Space Weather Prediction Center. In Figure 14 we compare the mean attenuation over the variable sky model simulations at 100 days with the solar cycle as tracked by the number of sunspots observed in each month.

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** Illustration of how the ionospheric attenuation in our simulations follows the solar cycle by overlaying sunspot number data with the inverse of the expected attenuation factors obtained in our sky model variance simulations. Note the two different vertical axes on the left and right. The thin blue and thin black curves show the historical monthly sunspot numbers from the NOAA Space Weather Prediction Center (https://www.swpc.noaa.gov/products/solar-cycle-progression) over solar cycles 23 and 24. The black points show the count for each month, while the blue curve is a smoothed version of the data obtained by averaging the counts over 13 months. Overlaid are the inverses of the mean over sky realizations of the power spectrum attenuation for each band obtained in the sky model variance simulations at N = 100, i.e., the last point on the curves in Figures 11–13.
Download figure:
Standard image High-resolution image

It turns out that the years 2011–2014 approximately span the peak of the current solar cycle 24. In contrast, the year 2008—the year of the previous solar minimum—exhibits little or no attenuation even after 100 days. The HERA array should reach its full complement of antennas and observe for 3 yr from 2021 to 2023. Thus, we expect that the next few seasons will pass through solar minimum and be near solar maximum by the end of data taking.

4.4. Fluctuating Polarized Power in the Vokes Parameters

The fluctuation of the ionospheric Faraday rotation turns the otherwise constant (for each t) polarization leakage term into something resembling a stochastic noise term, which is suppressed similarly to the thermal noise by averaging over multiple days of observation. However, unlike the thermal noise, we should not expect the Faraday-rotated polarization to be optimally suppressed by averaging over all the available visibility samples.

For example, this is obvious in the attenuation curve for the year 2012 in Figure 10, where there is a significant lack of monotonicity in the attenuation curves. In this case, if one had visibilities for only the first 40 days, the optimal attenuation is not obtained by averaging over all 40 days; more attenuation could be obtained by only including the first ∼20 days in the average, i.e., by using only half the available data. The reason for this can be observed in Figure 5, where the value of the RM for the year 2012 between days 20 and 40 can be seen to have a corresponding trend reversal, resulting in more coherent averaging of the polarization leakage over this range of days and thus less attenuation. On the other hand, because of the linear trend over tens of days, it is reasonable to expect that an average over a set of days with more separation between them will produce a greater attenuation, since the difference between each day’s RM tends to increase with the number of intervening days.

Given the set of indices of consecutive sidereal days

it is then likely that there are subsets ${S}_{k}\subset S$ of nonconsecutive days that will produce significantly more attenuation of the polarization terms in the Vokes parameters than the natural attenuation produced by summing over all consecutive days. There will equivalently be subsets that produce significantly less attenuation or even amplification of the polarization terms.

Figure 15 shows some distributions of ${\xi }_{I}({S}_{k},{ \mathcal B })$ obtained from the fiducial visibility simulations over a collection ${\mathfrak{C}}={\{{S}_{k}\}}_{k=1}^{{N}_{s}}$ of ${N}_{s}={10}^{6}$ subsets of S for which the natural attenuation is shown in Figure 10. The number of elements N in each subset S_k is held fixed at N = 50—that is, the elements of each S_k are drawn from S without replacement.

Figure 15. Refer to the following caption and surrounding text. — **Figure 15.** Left: distributions of ${\xi }_{I}({S}_{k},{ \mathcal B })$ . Right: distributions of ${\xi }_{L}({S}_{k},{ \mathcal B })$ . The distributions are obtained by taking 10⁶ random subsets S_k each of length N = 50 days (i.e., elements of S_k are drawn without replacement) from the fiducial HERA simulations in the years 2009 (top), 2011 (middle), and 2014 (bottom). The dashed lines indicate the 100-day natural attenuation value (i.e., S_k = S, N = 100) of ξ_I or ξ_L, not the mean of the distribution. Note the difference in range on the horizontal axes between the left and right panels.
Download figure:
Standard image High-resolution image

The existence of these wide distributions of attenuation factors suggests a null test that is particularly sensitive to polarization leakage in the delay spectrum estimator of the EOR power spectrum. If the thermal noise is sufficiently suppressed for a given subset size N, then the fluctuation of any problematic polarization leakage should dominate the variation of the different power spectrum estimates. Since variation in the ionospheric Faraday rotation of polarization leakage implies that there will be a distribution of delay power spectra similar to those in Figure 15, the absence of such a distribution rules out polarization leakage as the source of a detection above the expected thermal noise. The converse is not necessarily true—there may be other sources of contamination that could also produce such a distribution, so the presence of such a distribution does not imply that the excess power is due to polarization on the sky.

The best-case scenario is that the magnitude of any polarization leakage is below the cosmological signal level. In this case the effect of ionospheric fluctuations would never be observed directly in the Vokes I spectrum. Therefore, as a consistency check we will want to simultaneously observe the ionospheric fluctuation of the Vokes Q and U parameters. Observing fluctuating polarized power in the Vokes Q/U parameters will then show that there is a fluctuation that would have been observed in Vokes I if the magnitude of polarization leakage had been large enough.

While considering the delay spectra of ${{ \mathcal V }}_{Q}$ and ${{ \mathcal V }}_{U}$ individually can be useful for assessing instrumental systematics (Kohn et al. 2016), we can continue to proceed by analogy to the Stokes parameters to define a quantity analogous to L² = Q² + U² that is maximally sensitive to the magnitude of linear polarization on the sky:

For real data this will include a bias due to Stokes I and the instrumental polarization impurity that is unaffected by the changing Faraday rotations (see Equation 44(b)), and at low frequencies, where the polarization fraction on the sky is small, Stokes I will be at least even with Stokes Q/U in contribution to Vokes Q/U. For the purpose of assessing the attenuation of observed polarized power in our simulation, we follow the same approach here as was applied to Vokes I leakage and subtract the Stokes I term from Vokes Q and Vokes U to isolate the terms that are sensitive to the ionospheric Faraday rotation.

We then define a Vokes polarization delay power spectrum

and thus the relative attenuation factor for the Vokes polarization power

Examples of ξ_L as computed from our simulations are also shown in Figure 15—the ξ_L distributions are generated from the exact same collection of subsets S_k that was used for the ξ_I distributions. We can see that the peaks in the distributions over ${\mathfrak{C}}$ and the N = 100 values of ξ_L (dashed vertical lines in the plots) are comparable to the same quantities for ξ_I, but they are not perfectly correlated. It is notable that the variance of the ξ_L distributions is significantly larger than those of ξ_I, because the definitions of ξ_I and ξ_L already divide out an absolute magnitude. Note that this increased variance is not symmetric relative to the peaks of the ξ_I distributions—the distributions of ξ_L are relatively skewed toward smaller values. It appears that it is easier to find a combination of Faraday rotations that can significantly attenuate the Vokes polarization power spectrum than it is for the Vokes I polarization leakage.

Another potential issue with this null test is the computational cost in sampling the collection of possible subsets. The number of possible subsets S_k of any length is ${2}^{{N}_{d}}$ , which is ∼10³⁰ for N_d = 100, while the number of subsets of size N = 50 is $\tfrac{100!}{50!50!}\sim {10}^{29}$ . It is thus not possible to construct the full distributions over all possible subsets even in simulation, much less for actual measurements, which generally require more processing to produce power spectrum estimates. In reality, the number of samples N_s that can be used may only be in the hundreds or thousands.

The distributions in Figure 15 were generated by uniformly sampling the full collection of subsets of length 50. But we are not ignorant of the ionosphere’s behavior, and we would like to be able to use this knowledge to bias the sampling toward the tails of these distributions.

In principle, we might like to use simulations such as the ones done in this paper to inform a selection of subsets with which to estimate the power spectrum. But we have already seen how sensitive the attenuation is to the detailed coupling of the polarization states of the sky to the instrumental response, which casts doubt on our ability to make accurate predictions from simulations given our current levels of knowledge about the relevant functions. Fortunately, for the purpose of this null test we do not need perfect accuracy, only to do a little better than completely random guessing. Additionally, we need not precisely predict the actual magnitude of the attenuation in a particular subset, only its relative place in the distribution.

We attempt to approximate the distribution of ${\xi }_{I}({S}_{k},{ \mathcal B })$ in our fiducial HERA simulations by a simpler functional of the ionospheric RM that is independent of the observed sky and includes only an approximate and generic model of the instrumental response. Define

where ν_* is the central frequency of each sub-band ${ \mathcal B }$ , the function A² is given by Equation (42), and the sums are computed over an nside = 8 HEALPix map, as the integrand does not vary as much on small scales as the functions in our visibility calculation. The Mueller matrix elements used are those of the analytically defined Airy dipole model, computed from the definition of this Jones matrix (Equation (32)) and the formula for the Mueller matrix elements (Equation (20)).

This quantity need not predict the value of the attenuation precisely. We are only interested here in finding subsets that correspond to attenuation factors in the tails of the distributions of ξ_I in Figure 15. Thus, to compare the distributions for ξ_I, ξ_L, and ${{ \mathcal A }}^{2}$ , we compute the z-scores for each variable from the distribution over the chosen collection ${\mathfrak{C}}$ of subsets S_k. The z-score for the variable $X\in \{{\xi }_{I}({S}_{k},{ \mathcal B })$ , ${\xi }_{L}({S}_{k},{ \mathcal B })$ , ${{ \mathcal A }}^{2}({S}_{k},{ \mathcal B })\}$ is

where Mean() and Std() are the mean and standard deviation of X over ${\mathfrak{C}}$ , respectively. We denote the z-scores for each of these variables by ${{ \mathcal Z }}_{I}$ , ${{ \mathcal Z }}_{L}$ , ${{ \mathcal Z }}_{A}$ , respectively.

Figure 16 shows the correlation of ${{ \mathcal Z }}_{A}({S}_{k},{ \mathcal B })$ with ${{ \mathcal Z }}_{I}({S}_{k},{ \mathcal B })$ in the same years and for the same collection of subsets as used in Figure 15. We can see that the subsets that produce values of ${{ \mathcal A }}^{2}$ in the tails of the distribution tend to also find values of ξ_I in the tails of the distribution. The correlation is far from perfect, but as noted, the point is merely to improve the statistical power of the null test—any correlation helps compared to completely uniform sampling. Figure 17 then shows how ${{ \mathcal Z }}_{I}$ and the cuts on ${{ \mathcal Z }}_{A}$ are correlated with ${{ \mathcal Z }}_{L}$ . Additionally, the proxy function ${{ \mathcal A }}^{2}$ is simply an inspired guess based on Equation 44(b). It seems likely that an improved method of sampling these distributions based on the ionospheric RM data could be found; in particular, we have not used the fact that the RM has a significant trend as a function of n.

Figure 16. Refer to the following caption and surrounding text. — **Figure 16.** Correlation of the ionospheric fluctuation tracer ${{ \mathcal A }}^{2}$ with the Vokes I polarization leakage in the fiducial visibility simulation for the years 2009 (top), 2011 (middle), and 2014 (bottom). The collection of 10⁶ subsets used here is the same as the one used to make Figure 15. The red points show a cut on the 500 largest and smallest values of ${{ \mathcal Z }}_{A}$ for each ${ \mathcal B }$ . Such a cut would select the subsets to be used to estimate the power spectrum in our null test.
Download figure:
Standard image High-resolution image

**Figure 16.** Correlation of the ionospheric fluctuation tracer ${{ \mathcal A }}^{2}$ with the Vokes I polarization leakage in the fiducial visibility simulation for the years 2009 (top), 2011 (middle), and 2014 (bottom). The collection of 10⁶ subsets used here is the same as the one used to make Figure 15. The red points show a cut on the 500 largest and smallest values of ${{ \mathcal Z }}_{A}$ for each ${ \mathcal B }$ . Such a cut would select the subsets to be used to estimate the power spectrum in our null test.
Download figure:
Standard image High-resolution image

Figure 17. Refer to the following caption and surrounding text. — **Figure 17.** Correlation of the ionospheric fluctuations in the Vokes polarization band power with the Vokes I polarization leakage band power from the fiducial HERA visibility simulations for the years 2009 (top), 2011 (middle), and 2014 (bottom). The red points here correspond to the subsets from the cut on ${{ \mathcal Z }}_{A}$ , i.e., the red points in Figure 16. This shows that if we see a distribution in the Vokes polarization, we can infer that there exists a distribution in the Vokes I polarization, though we should not necessarily expect to see the *same* distribution.
Download figure:
Standard image High-resolution image

**Figure 17.** Correlation of the ionospheric fluctuations in the Vokes polarization band power with the Vokes I polarization leakage band power from the fiducial HERA visibility simulations for the years 2009 (top), 2011 (middle), and 2014 (bottom). The red points here correspond to the subsets from the cut on ${{ \mathcal Z }}_{A}$ , i.e., the red points in Figure 16. This shows that if we see a distribution in the Vokes polarization, we can infer that there exists a distribution in the Vokes I polarization, though we should not necessarily expect to see the *same* distribution.
Download figure:
Standard image High-resolution image

Figure 18. Refer to the following caption and surrounding text. — **Figure 18.** Example output from `radionopy`: the full-sky TEC content of the ionosphere (in TECU; factors of 10¹⁶ electrons per m²) on a Healpix grid, in this case projected onto a python `Basemap`. This particular snapshot shows the ionosphere at UT 0 hr on 2012 April 11.
Download figure:
Standard image High-resolution image

5. Discussion and Conclusions

1.
Ionospheric attenuation cannot be counted on to suppress polarization leakage in the power spectrum. Given what little is known about the level of polarized power on the sky in the 100–200 MHz frequency band, even at solar maximum it seems as likely as not that this attenuation might suppress polarization leakage to a negligible level. This increases the importance of precise modeling of this systematic, either to show that it will indeed be small relative to the EOR signal, or for the purpose of subtraction.
2.
Our simulations suggest a definitive test for polarization leakage in the power spectrum. This test comprises the following:
- (a)
  From the set S of N_d available sidereal days select a collection ${\mathfrak{C}}$ of subsets ${S}_{k}\subset S$ with the number N < N_d of elements in each S_k held fixed. The number N_d must be large enough to allow significant ionospheric variation over S. Additionally, the fraction N/N_d must be chosen to strike a balance between allowing the ionospheric attenuation to vary significantly between subsets and also ensuring that each subset represents sufficient integration time on the thermal noise.
- (b)
  Compute the power spectra P_I(S_k) and P_L(S_k) for each of the subsets. This produces a distribution of power spectra over ${\mathfrak{C}}$ .
- (c)
  If the Vokes I power spectrum estimator is dominated by Stokes I on the sky, then the changing ionospheric Faraday rotations between different subsets will have no effect and each subset will produce the same spectrum up to an expected distribution owing to the thermal noise.
- (d)
  The distribution of P_L should be significantly and obviously inconsistent with the expected thermal noise distribution.
- (e)
  The null test is passed when both 2(c) and (d) are satisfied, as 2(d) demonstrates that the effective polarized power on the sky has an observable variation over ${\mathfrak{C}}$ , while 2(c) shows that there is no corresponding variation of what is supposed to be Stokes I.
The method by which the elements of ${\mathfrak{C}}$ should be chosen remains open to further investigation. We have shown that a simple proxy function for ionospheric attenuation can reliably bias the sampling toward subsets with relatively high or low attenuation factors. Additional consideration could produce an improved method. The sensitivity of this test as a function of the thermal noise level is explored in a schematic way in Appendix C, but detailed consideration should be the subject of further simulations and analysis that can explore in detail the parameter space of cosmological signal level, thermal noise level, and polarized foreground power level. Additionally, the method of quantifying the consistency of these distributions with an expected thermal noise distribution need not be limited to simply computing the variance. For example, we showed that using our simple proxy function to select subsets can often produce distinctly bimodal distributions. The difference in the means of the high-attenuation collection and the low-attenuation collection could be a useful discriminating statistic. More generally, an advanced subset selection method may go hand in hand with a more robust way of distinguishing the resulting distributions from the expected thermal noise.
3.
The simulations we have used of the polarized sky are intended to be reasonably accurate representations of the expected sky, but their fidelity could certainly be improved. This is necessary for accurate prediction, since we have shown that the level of leakage is sensitively dependent on the correlated structure in the sky model and its alignment with the polarized antenna response, and this does produce large variations in the potential level of leakage. Given this uncertainty, we have purposely avoided considerations of the details of the absolute level of polarization leakage by considering ratios, and we demonstrate that these do show systematic trends independent of the details of the sky model.
4.
Averaging over sidereal days at fixed LST may still be a useful method for suppressing polarized foregrounds even in the situation in which one tries to model and subtract them directly from the visibilities, as the residual (unmodeled) polarization leakage will be attenuated by averaging over many days. This may ease the requirements on the completeness of the polarized model. On the other hand, an increasing level of ionospheric attenuation goes hand in hand with increasing complexity of the ionosphere, and thus increasing complexity of the model that must be constructed in order to perform the subtraction. It remains to be seen whether the global model of the ionospheric Faraday rotation that we have presented here would be adequate for such a task.
5.
The variance in the visibility and resulting power spectrum can be quite large when the polarization angle on the sky is not constrained. While preliminary, the results of our simulations suggest that a statistical foreground model that does not constrain the orientation of the polarization on the sky may be inadequate for predicting polarization leakage levels to the accuracy required for HERA, and possibly other EOR experiments. Determining the extent to which this is true or not through more careful consideration of the parameterization of the sky model and the mapping into the visibility will require further research. Obviously, it is necessary to determine the polarization angle accurately to be able to subtract a model from the visibilities.

This material is based on work supported by the National Science Foundation under grant nos. 1440343 and 1636646, the Gordon and Betty Moore Foundation, and institutional support from the HERA collaboration partners. S.A.K. is supported by a University of Pennsylvania SAS Dissertation Completion Fellowship. J.E.A. acknowledges support from NSF CAREER award no. 1455151.

Appendix A: Comparing Ionospheric RM Outputs

There are now several software packages that interpret CODE ionex files specifically for the use of low-frequency radio interferometers. Two of these are ionFR (Sotomayor-Beltran et al. 2013) and the results shown in Arora et al. (2015). In Figure 19 we show qualitative agreement with both of these works by comparing maps of vertical TEC values over the globe. In Figure 20 we show radionopy and ionFR RM output for a single pointing toward Cassiopeia A (Cas A; R.A. = 23^h23^m279, decl. = +58°48′42 farcs 4) from the LOFAR Core site in the Netherlands, which exhibit quantitative agreement. Slight offsets at the highest RM values that day can be attributed to differences in our interpolation schemes.

Figure 19. Refer to the following caption and surrounding text. — **Figure 19.** Top: vertical TEC from the CODE ionex file for 2011 April 11, overplotted on the globe in a Cartesian projection, as measured in Sotomayor-Beltran et al. (2013) and Arora et al. (2015) (left and right, respectively). Bottom: `radionopy` output for the same times and day. There is qualitative agreement, save for an error resulting in upside-down maps in Sotomayor-Beltran et al. (2013), as pointed-out by Arora et al. (2015).
Download figure:
Standard image High-resolution image

Figure 20. Refer to the following caption and surrounding text. — **Figure 20.** RM of Cas A as viewed from the LOFAR Core site in the Netherlands on 2011 April 11, according to `ionFR` and `radionopy`. The two codes show quantitative agreement; this demonstrates that `radionopy` can be used for single-pointing as well as full-sky RM measurements.
Download figure:
Standard image High-resolution image

Appendix B: The Instrumental Jones Matrix and Basis Transformation

While the instrumental Jones matrix ${\boldsymbol{J}}$ happens to be a 2 × 2 matrix in the case of the antenna with two different feed polarizations, it is better thought of as a list of rank-1 tensor fields ${{\boldsymbol{ \mathcal F }}}_{k}(\nu ,\hat{{\boldsymbol{s}}})$ for the kth feed of N feeds

that correspond to the far-field electric vector fields generated by each feed operated in transmission. Each row of the matrix corresponds to the normalized electric vector field pattern of a single feed of the antenna.

In order to compute Equation 12(a), the instrumental Jones matrix ${\boldsymbol{J}}$ and the coherency matrix ${\boldsymbol{ \mathcal C }}$ must be specified in the same basis. Because of the cylindrical symmetry of the ${\hat{{\boldsymbol{e}}}}_{\alpha }$ , ${\hat{{\boldsymbol{e}}}}_{\delta }$ basis, we can specify the instrumental response in this basis, rather than the alternative of performing a basis transformation on the observed coherency matrix. Observe that the integrand in Equation 12(a) is invariant under a transformation

where ${\boldsymbol{ \mathcal U }}(\hat{{\boldsymbol{s}}})$ is a 2 × 2 unitary matrix field. Any basis transformation (a point-by-point 2 × 2 rotation) is such a unitary matrix. Since the coherency matrix is specified in the ${\hat{{\boldsymbol{e}}}}_{\alpha },{\hat{{\boldsymbol{e}}}}_{\delta }$ basis, we thus require the instrumental response to be specified in this basis.

However, it is generally practical to specify the instrumental response in a basis of spherical coordinates local to the antenna so that the representation is independent of the telescope’s geographic location, and a standard choice of coordinates is the zenith angle θ ∈ (0, π) and local azimuthal angle $\phi \in [0,2\pi )$ . Explicitly, this means that the electric field data generated from an EM simulation of the antenna are specified as the complex coefficient functions of the vector field

which defines the instrumental response as

where ${\hat{{\boldsymbol{s}}}}_{b}$ denotes the direction of the antenna’s boresight. There is an equivalent representation of this vector field ${\boldsymbol{ \mathcal F }}$ in the equatorial basis

The components in the two different bases are related by

which defines a rotation matrix field ${ \mathcal U }(\hat{{\boldsymbol{s}}})$ with elements

For two feeds a and b with the response of each given by the vector fields ${{\boldsymbol{ \mathcal F }}}_{a}$ and ${{\boldsymbol{ \mathcal F }}}_{b}$ , the instrumental Jones matrix is then specified in the equatorial basis as

Appendix C: Effect of Thermal Noise in Polarization Null Test

Since we have not included the effect of thermal noise or an absolute scale for the polarized power in our analysis, we consider a schematic model of how these variables would affect the statistics of the proposed null test. The point is to argue that if polarization leakage were the limiting systematic in the power spectrum, the variance in our null test due to fluctuations in the polarized power will eventually dominate the variance due to thermal noise.

Let ${{ \mathcal P }}_{I}$ be the Stokes I contribution to the power spectrum, ${{ \mathcal P }}_{L}$ the contribution of Stokes Q and U, and ${ \mathcal N }$ the thermal noise with mean $\langle { \mathcal N }\rangle =0$ ; for simplicity of exposition we neglect cross-terms between Stokes parameters. The power spectrum

can then be considered a random variable over the collection ${\mathfrak{C}}$ of subsets of sidereal days, as each subset produces a different realization of the noise, and changing ionospheric attenuation produces a fluctuation in ${{ \mathcal P }}_{L}$ . The ${{ \mathcal P }}_{I}$ term, which represents the cosmological signal, is taken to be constant over the subsets. If ${\widehat{{ \mathcal P }}}_{L}$ is the intrinsic polarized power, then the attenuation factor is

where $\overline{\xi }$ is defined by the mean of ${{ \mathcal P }}_{L}$ over ${\mathfrak{C}}$ ,

The polarized power can also be written as

so we can see that

The mean and variance of ${ \mathcal P }$ are then

If ${{ \mathcal P }}_{I}\gg {\overline{{ \mathcal P }}}_{L}$ , then we detect the cosmological signal with an uncertainty dominated by the thermal noise and any other small systematics. If ${{ \mathcal P }}_{I}\ll {\overline{{ \mathcal P }}}_{L}$ , then we can see that the ionospheric fluctuation of ${{ \mathcal P }}_{L}$ in our null test will dominate the variance due to thermal noise—we will have been thwarted from observing cosmological reionization, but we will not be fooled into thinking otherwise.

In a regime where polarization leakage is comparable to the cosmological signal, we would have ${\overline{{ \mathcal P }}}_{L}\approx {{ \mathcal P }}_{I}$ , and thus the second term is approximately the thermal-noise-to-signal ratio on a detection in the absence of ${{ \mathcal P }}_{L}$ . The HERA experiment is designed to detect the EOR power spectrum at high signal-to-thermal-noise ratio, so even in a regime where ${\overline{{ \mathcal P }}}_{L}$ is slightly smaller than ${{ \mathcal P }}_{I}$ the excess variance in the null test should still be detectable.

Footnotes

1
https://healpix.sourceforge.net
2
https://github.com/UPennEoR/IonRIME. Note that this is intended only to provide documentation for the results, as the code does not constitute a user-friendly software package.
3
https://sourceforge.net/projects/ionfarrot/
4
https://github.com/UPennEoR/radionopy
5
https://github.com/radiocosmology/cora

Please wait… references are loading.