PAPER-64 CONSTRAINTS ON REIONIZATION: THE 21 cm POWER SPECTRUM AT z = 8.4

Redundant calibration has gained attention recently as a particularly powerful way to solve for internal dof in radio interferometric measurements without simultaneously having to solve for the distribution of sky brightness (Wieringa 1992; Liu et al. 2010; Noorishad et al. 2012; Marthi & Chengalur 2014; Zheng et al. 2014; P14). The grid-based configuration of PAPER antennas allows a large number of antenna calibration parameters to be solved for on the basis of redundancy (P12a; P14; Zheng et al. 2014). Multiple baselines of the same length and orientation measure the same sky signal. Differences between redundant baselines result from differences in the signal chain, including amplitude and phase effects attributable to antennas, cables, and receivers. Redundant calibration only constrains the relative complex gains between antennas and is independent of the sky. Since redundant calibration preserves signals common to all redundant baselines, this type of calibration does not result in signal loss.

In practice, redundant calibration often takes on two flavors: log calibration (LOGCAL) and linear calibration (LINCAL) (Liu et al. 2010; Zheng et al. 2014). LOGCAL uses logarithms applied to visibilities,

$$\begin{eqnarray}&&{v}_{{ij}}={g}_{i}^{*}{g}_{j}{y}_{i-j}+{n}_{{ij}}^{\mathrm{res}},\end{eqnarray} \tag{ 1 }$$

where g denotes the complex gain of antennas indexed by i and j, and y represents the "true" visibility measured by the baseline, to give a linearized system of equations

$$\begin{eqnarray}&&\mathrm{log}{v}_{{ij}}=\mathrm{log}{g}_{i}^{*}+\mathrm{log}{g}_{j}+\mathrm{log}{y}_{i-j}.\end{eqnarray} \tag{ 2 }$$

In solving for per-antenna gain parameters with a number of measurements that scales quadratically with antenna number, redundancy gives an over-constrained system of equations that can be solved using traditional linear algebra techniques. While LOGCAL is useful for arriving at a coarse solution from initial estimates that are far from the true value, LOGCAL has the shortcoming of being a biased by the asymmetric behavior of additive noise in the logarithm (Liu et al. 2010).

LINCAL, on the other hand, uses a Taylor expansion of the visibility around initial estimates of the gains and visibilities,

$$\begin{eqnarray}&&{v}_{{ij}}={g}_{i}^{0*}{g}_{j}^{0}{y}_{i-j}^{0}+{g}_{i}^{1*}{g}_{j}^{0}{y}_{i-j}^{0}+{g}_{i}^{0*}{g}_{j}^{1}{y}_{i-j}^{0}+{g}_{i}^{0*}{g}_{j}^{0}{y}_{i-j}^{1},\end{eqnarray} \tag{ 3 }$$

where 0 denotes initial estimates and 1 represents the perturbation to the original estimate and is the solutions we fit for. Using initial estimates taken from LOGCAL, LINCAL constructs an unbiased estimator.

Redundant calibration was performed using OMNICAL¹⁷ —an open-source redundant calibration package that is relatively instrument agnostic (Zheng et al. 2014). This package implements both LOGCAL and LINCAL, solving for a complex gain solution per antenna, frequency, and integration. The solutions are then applied to visibilities and the results are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to solving for gain solutions, OMNICAL also characterizes the quality of the calibration parameters by calculating the χ² for every integration. As defined in Zheng et al. (2014),

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{{ij}}\displaystyle \frac{{| {v}_{{ij}}-{y}_{i-j}{g}_{i}^{*}{g}_{j}| }^{2}}{{\sigma }_{{ij}}^{2}},\end{eqnarray} \tag{ 4 }$$

where σ² is the noise in the visibilities. The χ² measures sum of the deviation of measured visibilities to that of the best fit model derived from the LINCAL relative to a noise model, and gives us a tool to use in order to check the quality of our data. The number of dof, as defined in Zheng et al. (2014), is given by

$$\begin{eqnarray}\mathrm{dof} & = & {N}_{\mathrm{measurements}}-{N}_{\mathrm{parameters}}\\ & = & 2{N}_{\mathrm{baselines}}-2({N}_{\mathrm{antennas}}+{N}_{\mathrm{unique}\ \mathrm{baselines}}),\end{eqnarray} \tag{ 5 }$$

and is effectively the number of visibilities for which χ² is calculated. If the data are noise-dominated, χ²/dof is drawn from a χ² distribution with μ = 1 and σ² = 2/dof. The calculated χ²/dof for every frequency and integration of a fiducial day of observation in this season and for the fiducial power spectrum baselines is shown in Figure 5, demonstrating the stability of the PAPER instrument.

Zoom In Zoom Out Reset image size

We measure a mean χ²/dof of 1.9. This indicates that the redundant calibration solutions, while a substantial improvement over the previous PAPER-32 calibration (P14), do not quite result in residuals that are thermal noise dominated. Possible sources of this excess include instrumental crosstalk and poorly performing signal chains. While the latter will be down-weighted by the inverse of the estimated signal covariance described in Section 5, crosstalk is a defect in the data that must be addressed. Crosstalk caused by the cross-coupling of signals between antennas reveals itself as a static complex bias to a visibility that varies on timescales much longer than typical fringe rates. This effect skews the distribution of the χ² of the residuals away from 1. To minimize crosstalk, we first use OMNICAL to solve for antenna-dependent gains, and then average the residual deviations from redundancy over 10 minute windows before subtracting the average from the original visibilities. This crosstalk removal preserves signals common to redundant baseline groups (such as the 21 cm signal). Unfortunately, it also preserves a term that is the average of the crosstalk of all baselines in the redundant group. This residual crosstalk is removed by a fringe-rate filter later in the analysis.

After solving for the relative complex gains of the antennas using redundant calibration, an overall phase and gain calibration remains unknown. We use the standard self calibration method for radio interferometers to solve for the absolute phase calibration. We used Pictor A, Fornax A, and the Crab Nebula to fit for the overall phase solutions. Figure 6 shows an image of the field with Pictor A (5:19:49.70, −45:46:45.0) and Fornax A (3:22:41.70, −37:12:30.0).

Zoom In Zoom Out Reset image size

We then set our over all flux scale by using Pictor A as our calibrator source with source spectra derived in Jacobs et al. (2013),

$$\begin{eqnarray}&&{S}_{\nu }={S}_{150}\times {\left(\displaystyle \frac{\nu }{150\;\mathrm{MHz}}\right)}^{\alpha },\end{eqnarray} \tag{ 6 }$$

where S₁₅₀ = 381.88 Jy ± 5.36 and α = −0.76 ± 0.01, with 1σ error bars.

To derive the source spectrum from our measurements, we use data that have been LST-averaged prior to the wide-band delay filter described in Section 3.3, for the hour before and after the transit of Pictor A. We image a 30° × 30° field of view for every frequency channel for each 10 minute snapshot and apply uniform weights to the gridded visibilities. We account for the required three-dimensional Fourier transform in wide field imaging by using the w-stacking algorithm implemented in WSclean (Offringa et al. 2014) although we note that the standard w-projection algorithm implemented in CASA¹⁸ gives similar performances as the PAPER array is essentially instantaneously coplanar. A source spectrum is derived for each snapshot by fitting a two-dimensional Gaussian to Pictor A by using the PyBDSM¹⁹ source extractor. Spectra are optimally averaged together by weighting them with the primary beam model evaluated in the direction of Pictor A. To fit our bandpass, we divide the model spectrum by the measured one and fit a 9th order polynomial over the 120–170 MHz frequency range. Figure 7 shows the calibrated Pictor A spectrum and the model spectrum from Jacobs et al. (2013). Also plotted are the 1σ error bars derived from the PyBDSM source extractor and averaged over the multiple snapshots used after being weighted by the beam-squared.

Zoom In Zoom Out Reset image size

Fitting a polynomial to the bandpass has the potential for signal loss which would include suppressing modes that may contain the cosmological signal. In order to quantify the signal loss associated with fitting a ninth degree polynomial to the bandpass, we run a Monte Carlo simulation of the effect the bandpass has on a model 21 cm reionization signal. We construct a model baseline visibility as a Gaussian random signal multiplied by the derived bandpass for every independent mode measured. We calculate the total number of independent modes by counting the number of independent uv-modes sampled for the different baseline types over the two hour time interval used to measure the bandpass. We average each mode together and fit a 9th degree polynomial. Using this as our measured bandpass for this simulated signal, we finally compare the power spectrum from the output of the simulated signal to the input power spectrum as a function fo k-mode. We find that between $-0.06\lt k\lt 0.06$, the width of our wideband delay filter described below, the signal loss is less than 3% and at the mode right outside the above limit is 2 × 10⁻⁷%. We apply the latter correction factor for all modes outside the width of the delay filter to the final power spectrum.

Before implementing our foreground removal techniques, we combine the two linear polarizations for an estimate of Stokes I as per Moore et al. (2013). Namely, Stokes I can be estimated as

$$\begin{eqnarray}&&{V}_{{\rm{I}}}=\displaystyle \frac{1}{2}({V}_{\mathrm{XX}}+{V}_{\mathrm{YY}}),\end{eqnarray} \tag{ 7 }$$

where ${V}_{\mathrm{XX}}$ and V_YY are the visibilities of the two linear polarizations measured by the interferometer. There are some important caveats to the estimate of Stokes I provided by Equation (7). One important caveat is that it neglects the beam asymmetry between the two linear polarization states. This mismatch can cause polarization leakage from Stokes Q into Stokes I, thus contaminating our measurement of the power spectrum with any polarized emission from the sky. This effect for PAPER, as shown in Moore et al. (2013), leaks 4% of Q in to I in amplitude (2.2 × 10⁻³ in the respective power spectra). We take the conservative approach and do not correct for this effect, noting that the leakage of Q in to I will result in positive power, increasing our limits.

Foreground removal techniques discussed in the literature include spectral polynomial fitting (Wang et al. 2006; Bowman et al. 2009a; Liu et al. 2009a), principal component analysis (Liu & Tegmark 2011; Paciga et al. 2011; Masui et al. 2013; Paciga et al. 2013), non-parametric subtractions (Harker et al. 2009; Chapman et al. 2013), and inverse covariance weighting (Liu & Tegmark 2011; Dillon et al. 2013, 2014; Liu et al. 2014a, 2014b), Fourier-mode filtering Petrovic & Oh (2011), and per-baseline delay filtering described in P12b. This delay-spectrum filtering technique is well-suited to the maximum redundancy PAPER configuration, which is not optimized for the other approaches where high fidelity imaging is a prerequisite. The delay-spectrum foreground filtering method is described in detail by P14; its application is unchanged here. In summary; we Fourier transform each baseline spectrum into the delay domain

$$\begin{eqnarray}{\tilde{V}}_{\tau } & = & \displaystyle \int {W}_{\nu }{A}_{\nu }{I}_{\nu }{e}^{-2\pi i{\tau }_{g}}\cdot {e}^{2\pi i\tau \nu }d\nu \\ & = & {\tilde{W}}_{\tau }\ast {\tilde{A}}_{\tau }\ast {\tilde{I}}_{\tau }\ast \delta ({\tau }_{g}-\tau ),\end{eqnarray} \tag{ 8 }$$

where A_ν is the frequency dependent antenna response, W_ν is a sampling function that includes RFI flagging and a Blackman-Harris tapering function that minimizes delay-domain scattering from RFI flagging, and I_ν is the source spectrum. In the delay domain, a point source appears as a δ-function at delay τ_g, convolved by the Fourier transforms of the source spectrum, the antenna response, and the sampling function. We note that the antenna response effectively determines a finite bandpass, which imposes a lower bound of 1/B ≈ 10 ns on the width of any delay-domain convolving kernel. As per Parsons & Backer (2009) and P14, we deconvolve the kernel resulting from W(τ) using an iterative CLEAN-like procedure (Högbom 1974) restricting CLEAN components to fall within the horizon plus a 15 ns buffer that includes the bulk of the kernels convolving the δ-function in Equation (8). To remove the smooth spectrum foreground emission we subtract the CLEAN components from the original visibility.

Applying the delay filter to fiducial baselines used in the power spectrum analysis, foregrounds are suppressed by ∼4 orders of magnitude in power, or −40 dB of foreground suppression, as seen in Figure 8. As discussed in P14, there is a small amount of signal loss associated with this filter. For the baselines and filter parameters used, the loss was found to be 4.8% for the first mode outside of the horizon, 1.3% for the next mode out, and less than 0.0015% for the higher modes.

Zoom In Zoom Out Reset image size

After the wideband delay filter, we remove a second layer of RFI that was overshadowed by the foreground signal. RFI are excised with a filter which flags values 3σ above the median using a variance calculated in a localized time and frequency window.

We then average the entire season in LST with 43 s bin widths, matching the cadence of the compressed data. The full season was 135 days long; of these, 124 days were included in the average. We make two separate LST-binned data sets, averaging every other Julian day together to obtain an "even" and "odd" data set. The use of these two data sets allows us to construct an unbiased power spectrum estimate.

Sporadic RFI events result in measurements that, in any individual LST bin, deviate from the Gaussian distribution characteristic of thermal noise. To catch these events, we compute the median of a LST bin for each frequency and flag values 3σ above the median, before averaging. Since we are narrowing the distribution of visibilities about the median, the measured thermal noise variance is not preserved under this filter. However, since the central value is preserved, the expectation value of the measured visibility in each LST bin is unchanged, and there is no associated signal loss for power spectrum measurements. Moreover, because errors are estimated empirically through bootstrapping (see Section 5.4), the slight increase in measurement error associated with truncating the tails of the Gaussian distribution are naturally accounted for.

By averaging visibilities in time, we aim to maximize sensitivity by coherently combining repeated measurements of k-modes before squaring these measurements and averaging over independent k-modes to estimate the power spectrum amplitude. This is mathematically similar to the more traditional process of gridding in the uv plane, but applied to a single baseline. However, rather than applying a traditional box-car average, we can apply a kernel—a so-called "fringe-rate" filter—that weights different temporal rates by the antenna beam corresponding to the parts of the sky moving at the same rate.

For a given baseline and frequency, different parts of the sky exhibit different fringe-rates. Maximum fringe rates are found along the equatorial plane, where the rotation rate of the sky is highest, and zero fringe rates are found at the poles, where the sky does not rotate and hence sources do not move through the fringes of a baseline (Parsons & Backer 2009). Fringe rates are not constant as a function of latitude. Bins of constant fringe rate correspond to rings in R.A. and decl., where the east–west projection of a baseline projected toward a patch of the sky is constant. We use this fact in conjunction with the rms beam response for each contour of constant fringe rate to construct a time average kernel or "fringe-rate filter."

As examined in Parsons et al. (2015), it is possible to tailor fringe-rate filters to optimally combine time-ordered data for power-spectrum analysis. Fringe-rate filters can be chosen that up-weight points of the sky where our instrument is more sensitive and down-weight those points farther down in the primary beam, which are less sensitive. For white noise, all fringe-rate bins will contain the same amount of noise, but the amount of signal in each bin is determined by the primary beam response on the sky. By weighting fringe-rate bins by the rms of the beam response, we can get a net increase in sensitivity.

Applying this filter effectively weights the data by another factor of the beam area, changing the effective primary beam response,²⁰ A(l, m) (Parsons et al. 2015). By utilizing prior knowledge about the beam area, we are selectively down-weighting areas on the sky contributing little signal. This will result in a net improvement in sensitivity depending on the shape of the beam and the decl. of the array. For PAPER, this filter roughly doubles the sensitivity of our measurements.

Generally, a fringe-rate filter integrates visibilities in time. For a fringe-rate filter, f_fr, the effective integration time can be calculated by comparing the variance statistic before and after filtering:

$$\begin{eqnarray}&&{t}_{\mathrm{int},\mathrm{after}}={t}_{\mathrm{int},\mathrm{before}}\displaystyle \frac{\int {\sigma }_{f}^{2}{df}}{\int {\sigma }_{f}^{2}{f}_{\mathrm{fr}}^{2}{df}},\end{eqnarray} \tag{ 9 }$$

where ${t}_{\mathrm{int},\mathrm{before}}$ is the integration time before filtering, σ_f denotes the noise variance in fringe rate space and the integral is taken over all possible fringe rates for a given baseline and frequency. As discussed in Parsons et al. (2015), the signal re-weighting associated with this fringe-rate filter can be interpreted as a modification to the shape of the primary beam.

For the fiducial baseline at 151 MHz, the integration time, as given in Equation (9), associated with an optimal fringe rate filter is 3430 s. The number of statistically independent samples on the sky decreases from 83 to 1 sample per hour. As discussed in Section 5.3, empirically estimating a covariance matrix with a small number of independent samples can lead to signal loss in the OQE. In order to counteract the signal loss, we degrade the optimal fringe-rate filter, as shown in Figure 9, to have an effective integration time of 1886 s, increasing the number of independent modes to 2 per hour. The fringe rate filter is now sub-optimal, but is still an improvement on the boxcar weighting as used in P14. As documented in Table 1, the correction factor for the associated signal loss of the filter we have chosen is 1.39.

Zoom In Zoom Out Reset image size

We implement the modified filter on a per baseline basis by weighting the fringe-rate bins on the sky by the rms of the beam at that same location. In order to obtain a smooth filter in the fringe-rate domain, we fit a Gaussian with a hyperbolic tangent tail to this filter. In addition, we multiply this response with another hyperbolic tangent function that effectively zeros out fringe rates below 0.2 mHz. This removes the slowly varying signals that we model as crosstalk. We convolve the time-domain visibilities with the Fourier transform of the resulting fringe-rate filter, shown in Figure 9, to produce an averaged visibility. The effect on the data can be seen in Figure 8.

In order to build sensitivity to the 21 cm reionization signal, it is critical that PAPER be able to integrate coherently measurements made with different baselines on different days. Figure 10 shows the visibility repeatability between baselines and nights as a function of LST. Specifically, we histogram the real part of the visibilities for all redundant fiducial baselines in a given LST bin for foreground contained data. We see that for a given LST bin, the spread in values over all the baselines is ∼50 Jy which corresponds with our observed T_sys ∼ 500 K. We get more samples per LST bin in the range of 2–10 hr due to our observing season, therefore the density of points in this LST region is greater, as shown by the color scale. This density plot shows that redundant baselines agree very well with one another; OMNICAL has leveled the antenna gains to within the noise.

Zoom In Zoom Out Reset image size

Delving in a little deeper, we also examine the stability in time for measurements in a particular LST bin. In order to quantify the stability in time we extract one channel for a given baseline for every observation day and LST bin. We then Fourier transform along the time direction for every LST bin and compute the power spectrum. As shown in Figure 11, for timescales greater than one day, we see that signal variance drops by almost four orders of magnitude, with the exception of an excess on two-day timescales caused by the changing alignment of the 42.9 s integration timescale relative to a sidereal day. The implication of this measurement is that, after calibration, PAPER measurements are sufficiently stable to be integrated coherently over the entire length of a 135 day observation. This implies day-to-day stability of better than 1%, contributing negligibly to the uncertainties in the data.

Zoom In Zoom Out Reset image size

During the LST binning step, the variance of the visibilities that are averaged together for a given frequency and LST bin are recorded. Using these variances, we calculate the system temperature as a function of LST, averaging over each LST hour.

$$\begin{eqnarray}&&{T}_{\mathrm{rms}}={T}_{\mathrm{sys}}/\sqrt{2{\rm{\Delta }}\nu t},\end{eqnarray} \tag{ 10 }$$

where Δν is the bandwidth, t is the integration time, and T_rms is the rms temperature, or the variance statistic described above. Figure 12 shows the results of this calculation. In this observing season, the system temperature drops just below previous estimates as in P14 and Jacobs et al. (2015) of ${T}_{\mathrm{sys}}=560\;{\rm{K}}$, at ${T}_{\mathrm{sys}}=500\;{\rm{K}}$ at 160 MHz. However, this estimate is more consistent with the results derived in (Moore et al. 2015), where ${T}_{\mathrm{sys}}=505\;{\rm{K}}$ at 164 MHz. The change in the system temperature can be attributed to the reduced range of LST used in the calculation. We note that at 7:00 LST, there is an increase in the system temperature due to the rising of the galactic plane as seen in Figure 3.

Zoom In Zoom Out Reset image size

When calculating the system temperature using the variance in the visibilities for a given LST and frequency, we take into account the fact that we flag 3σ outliers from the median. To calculate an effective correction factor to account for the filtering, we assume the visibilities follow a Gaussian distribution which would require a correction factor of 1.34 for the removal of data points that are 3σ above the median. In other words, we are accounting for the wings of the Gaussian that would contribute to the variance in the visibility.

Previous estimates of the system temperature (P14; Jacobs et al. 2015) relied on differencing and averaging baselines, time samples, and/or frequency channels. The relative agreement between these various methods of estimating the system temperature provides a robust measure of the system temperature of the PAPER instrument. Agreement between the instantaneous measurements of the system temperature, the LST repetition variance, and the predicted power spectrum noise level (see below) indicates a robustly stable system with no significant long term instability contributing appreciable noise.

We use the OQE method to estimate our power spectrum as done in Liu & Tegmark (2011), Dillon et al. (2013), Liu et al. (2014a), Liu et al. (2014b), and Trott et al. (2012). Here we briefly review the OQE formalism with an emphasis on our application to data, which draws strongly from the aforementioned works, but also relies on empirical techniques similar to those used in P14. The end goal of this analysis is to estimate the 21 cm power spectrum, ${P}_{21}({\boldsymbol{k}})$, defined such that

$$\begin{eqnarray}&&\langle {\tilde{T}}_{b}({\boldsymbol{k}}){\tilde{T}}_{b}^{*}({{\boldsymbol{k}}}^{\prime })\rangle ={(2\pi )}^{3}{\delta }^{{\rm{D}}}({\boldsymbol{k}}-{{\boldsymbol{k}}}^{\prime }){P}_{21}({\boldsymbol{k}}),\end{eqnarray} \tag{ 11 }$$

where ${\tilde{T}}_{b}({\boldsymbol{k}})$ is the spatial Fourier transform of the brightness temperature distribution on the sky, $\langle \rangle$ denotes an ensemble average, and δ^D is the Dirac delta function.

In order to make an estimate of the power spectrum in the OQE formalism, one begins with a data vector ${\boldsymbol{x}}$. This vector could, for example, consist of a list of brightness temperatures on the sky for an imaging-based data analysis, or (in our case) a list of measured visibilities. We form the intermediate quantity,

$$\begin{eqnarray}&&{\hat{q}}_{\alpha }=\displaystyle \frac{1}{2}{{\boldsymbol{x}}}^{\dagger }{{\boldsymbol{C}}}^{-1}{{\boldsymbol{Q}}}_{\alpha }{{\boldsymbol{C}}}^{-1}{\boldsymbol{x}}-{b}_{\alpha },\end{eqnarray} \tag{ 12 }$$

which will be needed to form the OQE of our power spectrum. Here, ${\boldsymbol{C}}\equiv \langle {{\boldsymbol{xx}}}^{\dagger }\rangle$ is the true covariance matrix of the data vector ${\boldsymbol{x}},$ ${{\boldsymbol{Q}}}_{\alpha }$ is the operator that takes visibilities into power spectrum k-space and bins into the αth bin, and b_α is the bias to the estimate that needs to be subtracted off. In general, ${{\boldsymbol{Q}}}_{\alpha }$ represents a family of matrices, one for each k bin indexed by α. Each matrix is defined as ${{\boldsymbol{Q}}}_{\alpha }\equiv \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}$, i.e., the derivative of the covariance matrix with respect to the band power p_α. The bandpower ${p}_{\alpha }$ can be intuitively thought of as the value of the power spectrum in the αth k bin. Therefore, ${{\boldsymbol{Q}}}_{\alpha }$ encodes the response of the data covariance matrix to the αth bin of the power spectrum.

The bias term b_α in Equation (12) will include contributions from both instrumental noise and residual foregrounds. Their presence in the data is simply due to the fact that both contributions have positive power. One approach to dealing with these biases is to model them and to subtract them off, as is suggested by Equation (12). An alternate approach is to compute a cross-power spectrum between two data sets that are known to have the same sky signal but independent instrumental noise realizations. Labeling these two data sets as ${{\boldsymbol{x}}}_{1}$ and ${{\boldsymbol{x}}}_{2}$ and computing

$$\begin{eqnarray}&&{\hat{q}}_{\alpha }=\displaystyle \frac{1}{2}{{\boldsymbol{x}}}_{1}^{\dagger }{{\boldsymbol{C}}}^{-1}{{\boldsymbol{Q}}}_{\alpha }{{\boldsymbol{C}}}^{-1}{{\boldsymbol{x}}}_{2},\end{eqnarray} \tag{ 13 }$$

one arrives at a cross-power spectrum that by construction has no noise bias. There is thus no need to explicitly model and subtract any noise bias, although any residual foreground bias will remain, since it is a contribution that is sourced by signals on the sky, and therefore must exist in all our data sets.

The set of ${\hat{q}}_{\alpha }{\rm{s}}$ do not yet constitute a properly normalized estimate of the power spectrum (as evidenced, for example, by the extra factors of ${{\boldsymbol{C}}}^{-1}$). To normalize our results, we group the unnormalized bandpowers into a vector $\hat{{\boldsymbol{q}}}$ and apply a matrix ${\boldsymbol{M}}$ (whose exact form we specify later), so that

$$\begin{eqnarray}&&\hat{{\boldsymbol{p}}}={\boldsymbol{M}}\hat{{\boldsymbol{q}}}\end{eqnarray} \tag{ 14 }$$

is a normalized estimate $\hat{{\boldsymbol{p}}}$ of the true power spectrum ${\boldsymbol{p}}$. We emphasize that the vector space that contains $\hat{{\boldsymbol{q}}}$ and $\hat{{\boldsymbol{p}}}$ is an "output" vector space over different k-bins, which is separate from the "input" vector space of the measurements, in which ${\boldsymbol{x}}$ and ${\boldsymbol{C}}$ reside.

To select an ${\boldsymbol{M}}$ matrix that properly normalizes the power spectrum, we must compute the window function matrix ${\boldsymbol{W}}$ for our estimator. The window matrix is defined such that the true bandpowers ${\boldsymbol{p}}$ and our estimates $\hat{{\boldsymbol{p}}}$ of them are related by

$$\begin{eqnarray}&&\hat{{\boldsymbol{p}}}={\boldsymbol{Wp}},\end{eqnarray} \tag{ 15 }$$

so that each row gives the linear combination of the true power that is probed by our estimate. With a little algebra, one can show that

$$\begin{eqnarray}&&{\boldsymbol{W}}={\boldsymbol{MF}},\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\alpha \beta }=\displaystyle \frac{1}{2}{\rm{tr}}({{\boldsymbol{C}}}^{-1}{{\boldsymbol{Q}}}_{\alpha }{{\boldsymbol{C}}}^{-1}{{\boldsymbol{Q}}}_{\beta }),\end{eqnarray} \tag{ 17 }$$

which we have suggestively denoted with the symbol ${\boldsymbol{F}}$ to highlight the fact that this turns out to be the Fisher information matrix of the bandpowers. In order to interpret each bandpower as the weighted average of the true bandpowers, we require each row of the window function matrix to sum to unity. As long as ${\boldsymbol{M}}$ is chosen in such a way that ${\boldsymbol{W}}$ satisfies this criterion, the resulting bandpower estimates $\hat{{\boldsymbol{p}}}$ will be properly normalized.

Beyond the normalization criterion, a data analyst has some freedom over the precise form of ${\boldsymbol{M}}$, which effectively also re-bins the bandpower estimates. One popular choice is ${\boldsymbol{M}}={{\boldsymbol{F}}}^{-1}$, which implies that ${\boldsymbol{W}}={\boldsymbol{I}}$. Each window function is then a delta function, such that bandpowers do not contain leakage from other bins, and contain power from only that bin. However, the disadvantage of this becomes apparent if one also computes the error bars on the bandpower estimates. The error bars are obtained by taking the square root of the diagonal of the covariance matrix, which is defined as

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}=\mathrm{Cov}(\hat{{\boldsymbol{p}}})=\langle \hat{{\boldsymbol{p}}}{\hat{{\boldsymbol{p}}}}^{\dagger }\rangle -\langle \hat{{\boldsymbol{p}}}\rangle {\langle \hat{{\boldsymbol{p}}}\rangle }^{\dagger }.\end{eqnarray} \tag{ 18 }$$

Since $\hat{{\boldsymbol{p}}}={\boldsymbol{M}}\hat{{\boldsymbol{q}}}$, it is easily shown that

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={{\boldsymbol{MFM}}}^{\dagger }.\end{eqnarray} \tag{ 19 }$$

The choice of ${\boldsymbol{M}}={{\boldsymbol{F}}}^{-1}$ tends to give rather large error bars. At the other extreme, picking ${{\boldsymbol{M}}}_{\alpha \beta }\propto {\delta }_{\alpha \beta }/{{\boldsymbol{F}}}_{\alpha \alpha }$ (with the proportionality constant fixed by our normalization criterion) leads to the smallest possible error bars (Tegmark 1997), at the expense of broader window functions. In our application of OQEs in the following sections, we will pick an intermediate choice for ${\boldsymbol{M}}$, one that is carefully tailored to avoid the leakage of foreground power from low k modes to high k modes.

Here we describe the specifics of our application of the OQE formalism to measure the power spectrum. Doing so requires defining various quantities such as ${\boldsymbol{x}},$ ${\boldsymbol{C}}$, ${{\boldsymbol{Q}}}_{\alpha }$ for our analysis pipeline.

First, we consider ${\boldsymbol{x}},$ which represents the data in our experiment. Our data set consists of visibilities as a function of frequency and time for each baseline in the array. In our analysis, we group the baselines into three groups of redundant baselines (described in Section 2), in the sense that within each group there are multiple copies of the same baseline. In the description that follows, we first estimate the power spectrum separately for each group. Power spectrum estimates obtained from the different redundant groups are then combined in a set of averaging and bootstrapping steps described in Section 5.4. Note that because our data have been fringe-rate filtered in the manner described in Section 3.5, we may reap all the benefits of coherently integrating in time simply by estimating the power spectrum for every instant in the LST-binned data before averaging over the time-steps within the LST-binned day (Parsons et al. 2015).

For the next portion of our discussion, consider only the data within a single redundant group. Within each group there are not only multiple identical copies of the same baseline, but in addition (as discussed in Section 3.3), our pipeline also constructs two LST-binned data sets, one from binning all even-numbered days in our observations, and the other from all odd-numbered days. Thus, we have not a single data vector, but a whole family of them, indexed by baseline (i) and odd versus even days (r). Separating the data out into independent subgroups allows one to estimate cross-power spectra rather than auto-power spectra in order to avoid the noise bias, as discussed in the previous section. The data vectors take the form

$$\begin{eqnarray}&&{{\boldsymbol{x}}}^{{ri}}(t)=\left(\begin{array}{c}{V}^{{ri}}({\nu }_{1},t)\\ {V}^{{ri}}({\nu }_{2},t)\\ \vdots \end{array}\right),\end{eqnarray} \tag{ 20 }$$

where ${V}^{{ri}}(\nu ,t)$ is the visibility at frequency ν at time t. Each data vector is 20 elements long, being comprised of 20 channels of a visibility spectrum spanning 10 MHz of bandwidth centered on 151.5 MHz.

Having formed the data vectors, the next step in Equation (12) is to weight the data by their inverse covariance. To do so, we of course require the covariance matrix ${\boldsymbol{C}}$, which by definition, is the ensemble average of ${{\boldsymbol{xx}}}^{\dagger }$, namely ${\boldsymbol{C}}=\langle {{\boldsymbol{xx}}}^{\dagger }\rangle$. Unfortunately, in our case the covariance is difficult to model from first principles, and we must resort to an empirically estimated ${\boldsymbol{C}}$. We make this estimation by taking the time average of the quantity ${{\boldsymbol{xx}}}^{\dagger }$ over 8.5 hr of LST, estimating a different covariance matrix for each baseline and for odd versus even days. While an empirical determination of the covariance is advantageous in that it captures features that are difficult to model from first principles, it carries the risk of cosmological signal loss (Switzer & Liu 2014). We will discuss and quantify this signal loss in Section 5.3.

To gain some intuition for the action of ${{\boldsymbol{C}}}^{-1}$ on our data, let us examine the combination

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{{ri}}={({{\boldsymbol{C}}}^{{ri}})}^{-1}{{\boldsymbol{x}}}^{{ri}}\end{eqnarray} \tag{ 21 }$$

for select baselines. This is a crucial step in the analysis since it suppresses coherent frequency structures (such as those that might arise from residual foregrounds). Note that the inverse covariance weighting employed here differs from that in P14, in that P14 modeled and included covariances between different baselines, whereas in our current treatment we only consider covariances between different frequency channels. Figure 13 compares the effects of applying the inverse covariance matrix to a data vector that contains foregrounds (and thus contains highly correlated frequency structures) to one in which foregrounds have been suppressed by the wideband delay filter described in Section 3.3. In the figure, the top row corresponds to the data vector ${{\boldsymbol{x}}}^{{ri}}$ for three selected baselines in the form of a waterfall plot of visibilities, with frequency on the horizontal axis and time on the vertical axis. The middle section shows the empirical estimate of the covariance by taking the outer product of ${\boldsymbol{x}}$ with itself and averaging over the time axis. Finally, the last row shows the results of inverse covariance weighting the data, namely ${{\boldsymbol{z}}}^{{ri}}$. In every row, the foreground-dominated data are shown in the left half of the figure, while the foreground-suppressed data are shown in the right half.

Zoom In Zoom Out Reset image size

Consider the foreground-dominated ${{\boldsymbol{x}}}^{{ri}}$ in Figure 13, and their corresponding covariance matrices. The strongest modes that are present in the data are the eigenmodes of the covariance matrix with the largest eigenvalues. Figure 14 shows the full eigenvalue spectrum and the four strongest eigenmodes. For the foreground-dominated data, one sees that the eigenvalue spectrum is dominated by the first few modes, and the corresponding eigenmodes are rather smooth, highly suggestive of smooth spectrum foreground sources. The application of the inverse covariance weighting down-weights these eigenmodes, revealing waterfall plots in the bottom row of Figure 13 that look more noise-dominated. With the foreground-suppressed portion (right half) of Figure 13, the initial ${{\boldsymbol{x}}}^{{ri}}$ vectors already appear noise dominated (which is corroborated by the relatively noisy form of the eigenvalue spectra in Figure 14). The final ${{\boldsymbol{z}}}^{{ri}}$ vectors remain noise-like, although some smooth structure (perhaps from residual foregrounds) has still been removed, and finer scale noise has been up-weighted.

Zoom In Zoom Out Reset image size

With intuition established for the behavior of ${{\boldsymbol{C}}}^{-1}$, we may group our identical baselines into five different sets and average together ${{\boldsymbol{z}}}^{{ri}}$ vectors for baselines within the same set. That is, we form

$$\begin{eqnarray}&&{{\boldsymbol{z}}}_{A}^{r}=\displaystyle \sum _{i\in A}{({{\boldsymbol{C}}}^{{ri}})}^{-1}{{\boldsymbol{x}}}^{{ri}},\end{eqnarray} \tag{ 22 }$$

where A ranges from 1 to 5 and indexes the baseline set. At this point, we have 10 weighted data vectors ${\boldsymbol{z}}$ (5 baseline sets, each of which has an even day and odd day version) for every LST-binned time-step. As discussed in the previous section, instrumental noise bias may be avoided by forming cross-power spectra rather than auto-power spectra. Generalizing Equation (13) to our present case where we have 10 different data vectors, we have

$$\begin{eqnarray}&&{\hat{{\boldsymbol{q}}}}_{\alpha }=\displaystyle \sum _{\begin{array}{c}A,B,r,s\\ r\ne s,A\ne B\end{array}}{{\boldsymbol{z}}}_{A}^{r\dagger }{{\boldsymbol{Q}}}_{\alpha }{{\boldsymbol{z}}}_{B}^{s},\end{eqnarray} \tag{ 23 }$$

so that auto-power contributions from identical baseline groups or identical even/odd indices never appear. Residual foreground bias will remain in Equation (23), but in order to avoid possible signal loss from an overly aggressive foreground bias removal scheme, we conservatively allow the foreground bias to remain. Since foreground power will necessarily be positive, residual foregrounds will only serve to raise our final upper limits.

In order to implement Equation (23), it is necessary to derive a form for ${{\boldsymbol{Q}}}_{\alpha }\equiv \partial {\boldsymbol{C}}/\partial {p}_{\alpha }$. To do so, we follow the delay spectrum technique of P12a, where it was shown that

$$\begin{eqnarray}&&P({{\boldsymbol{k}}}_{t\tau })\approx {\left(\displaystyle \frac{{\lambda }^{2}}{2{k}_{{\rm{B}}}}\right)}^{2}\displaystyle \frac{{X}^{2}Y}{{\rm{\Omega }}B}\langle {\tilde{V}}_{i}(t,\tau ){\tilde{V}}_{j}^{*}(t,\tau )\rangle ,\end{eqnarray} \tag{ 24 }$$

where ${V}_{i}(t,\tau )$ is the delay transform of baseline visibilities given by Equation (8), X and Y are the constants that convert from angles and frequency to the co-moving coordinate, respectively, Ω is the power squared beam (see Appendix B of P14), B is the bandwidth, λ is the spectral wavelength, and k_B is Boltzmann's constant. This suggests that in order to estimate the power spectrum from visibilities, one only needs to Fourier transform along the frequency axis (converting the spectrum into a delay spectrum) before squaring and multiplying by a scalar. Thus, the role of ${{\boldsymbol{Q}}}_{\alpha }$ in Equation (23) is to perform a frequency Fourier transform on each copy of ${\boldsymbol{z}}$. It is therefore a separable matrix of the form ${{\boldsymbol{Q}}}_{\alpha }={{\boldsymbol{m}}}_{\alpha }{{\boldsymbol{m}}}_{\alpha }^{\dagger }$, where ${{\boldsymbol{m}}}_{\alpha }$ is a complex sinusoid of a specific frequency corresponding to delay mode α. We may thus write

$$\begin{eqnarray}&&{\hat{{\boldsymbol{q}}}}_{\alpha }=\displaystyle \sum _{\begin{array}{c}A,B,r,s\\ r\ne s,A\ne B\end{array}}{{\boldsymbol{z}}}_{A}^{r\dagger }{{\boldsymbol{m}}}_{\alpha }{{\boldsymbol{m}}}_{\alpha }^{\dagger }{{\boldsymbol{z}}}_{B}^{s}.\end{eqnarray} \tag{ 25 }$$

With an explicit form for ${{\boldsymbol{Q}}}_{\alpha }$, one now also has the necessary ingredients to compute the Fisher matrix using Equation (17).

Having computed the ${\hat{{\boldsymbol{q}}}}_{\alpha }{\rm{s}}$, we group our results into a vector $\hat{{\boldsymbol{q}}}$. This vector of unnormalized bandpowers is then normalized to form our final estimates of the power spectrum ${\boldsymbol{p}}$. As noted above, the normalization occurs by the ${\boldsymbol{M}}$ matrix in Equation (14), and can be any matrix of our desire. Even though the choices of the normalization matrices described above have certain good properties, e.g., small error bars or no leakage, we opt for a different choice of window function, as follows. We first reorder the elements in $\hat{{\boldsymbol{q}}}$ (and therefore in ${\boldsymbol{F}}$, ${\boldsymbol{M}}$, and $\hat{{\boldsymbol{p}}}$ for consistency) so that the k-modes are listed in ascending order, from low k to high k, with the exception that we place the highest k bin third after the lowest two k bins. (The reason for this exception will be made apparent shortly). We then take the Cholesky decomposition of the Fisher matrix, such that ${\boldsymbol{F}}={{\boldsymbol{LL}}}^{\dagger }$, where ${\boldsymbol{L}}$ is a lower triangular matrix. Following that, we pick ${\boldsymbol{M}}={{\boldsymbol{DL}}}^{-1}$, where ${\boldsymbol{D}}$ is a diagonal matrix chosen to adhere to the normalization constraint that ${\boldsymbol{W}}={\boldsymbol{MF}}$ has rows that sum to unity. In this case, the window function matrix becomes, ${\boldsymbol{W}}={{\boldsymbol{DL}}}^{\dagger }$. This means that ${\boldsymbol{W}}$ is upper triangular, and with our ordering scheme, has the consequence of allowing power to leak from high to low k, but not vice versa. Since our k axis is (to a good approximation) proportional to the delay axis, foregrounds preferentially appear at low k because their spectra are smooth. Reducing leakage from low k to high k thus mitigates leakage of foregrounds into the cleaner, more noise-dominated regions. Additionally, our placement of the highest k bin as the third element in our reordering of $\hat{{\boldsymbol{p}}}$ prevents leakage from this edge bin that will contain aliased power. Figure 15 shows the resulting window functions.

Zoom In Zoom Out Reset image size

Our choice of normalization matrix also has the attractive property of eliminating error correlations between bandpower estimates. Using Equation (19), we have that

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={{\boldsymbol{DL}}}^{-1}{{\boldsymbol{LL}}}^{\dagger }{{\boldsymbol{L}}}^{-\dagger }{\boldsymbol{D}}={{\boldsymbol{D}}}^{2}.\end{eqnarray} \tag{ 26 }$$

The error covariance matrix on the bandpowers is thus diagonal, which implies that our final data points are uncorrelated with one another. This stands in contrast to the power-spectrum estimator used in P14, where the Blackmann–Harris taper function induced correlated errors between neighboring data points.

We now discuss some of the subtleties associated with empirically estimating the covariance matrix from the data. Again, the covariance matrix is defined as the ensemble average of the outer product of a vector with itself, i.e.,

$$\begin{eqnarray}&&{\boldsymbol{C}}=\langle {{\boldsymbol{xx}}}^{\dagger }\rangle ,\end{eqnarray} \tag{ 27 }$$

where ${\boldsymbol{x}}$ is the data (column) vector used in the analysis. In our analysis, we do not have a priori knowledge of the covariance matrix. and thus we must resort to empirical estimates (Dillon et al. 2015). As we have alluded to above, we replace the ensemble average with a time average that runs from 0 to 8:30 LST hours.

Since the OQE method for power spectrum estimation requires the inversion of ${\boldsymbol{C}},$ it is crucial that our empirically estimated covariance be a full rank matrix. With our data consisting of visibilities over 20 frequency channels, the covariance matrix is a 20 × 20 matrix. Thus, a necessary condition for our estimate to be full rank is for there to be at least 20 independent time samples in our average. As noted in Section 3.5 the fringe-rate filter used corresponds to averaging time samples for 31 minutes. Over the LST range used in this analysis, this corresponds to roughly 20 statistically independent modes in our data after fringe-rate filtering. We therefore have just enough samples for our empirical estimate, and in practice, our covariance matrices are invertible and allow OQE techniques to be implemented.

Another potential problem that occurs from empirically estimating covariances is that it leads to models of the covariance matrix that over-fit the noise. In this scenario, the covariance matrix tells us that there may be modes in the data that should be down-weighted, for example, but if the empirical covariance estimates are dominated by noise, these may just be random fluctuations that need not be down-weighted. Said differently, the weighting of the data by the inverse covariance is heavily influenced by the noise in the estimate of the covariance matrix and thus has the ability to down-weight valid high-variance samples. Over-fitting the noise in this manner carries with it the possibility of cosmological signal loss. This seems to contradict the conventionally recognized feature of OQEs as lossless estimators of the power spectrum (Tegmark 1997). However, the standard proofs of this property assume that statistics such as ${\boldsymbol{C}}$ are known a priori, which is an assumption that we are violating with our empirical estimates.

In order to deal with possible signal loss, we perform simulations of our analysis pipeline, deriving correction factors that must be applied to our final constraints. We simulate visibilities for Gaussian temperature field with a flat amplitude in P(k) that rotates with the sky, which is fringe-rate filtered in the same way as the data for our fiducial baselines. This signal is processed through our pipeline, and the output power spectrum compared to the input power spectrum, for various levels of input signal amplitude. We repeat this for 40 sky realizations at each signal level. Figure 16 shows the resultant signal loss associated with estimating the covariance matrix from the data. Error bars were obtained through bootstrapping.

Zoom In Zoom Out Reset image size

As a function of the increasing input amplitude of the simulated power spectra, we find that the ratio of output power to input power decreases, which we interpret as signal loss through the use of our empirical OQE of the power spectrum. However, since the transfer function through this analysis is an invertible function, we can correct for the transfer by using the output value to infer a signal loss that is then divided out to obtain the original input signal level. In Figure 16, we see that deviations from unity signal transfer begin at power spectrum amplitudes of ${10}^{7}\;{\mathrm{mK}}^{2}{({h}^{-1}\;\mathrm{Mpc})}^{3}$. For the range of output power spectrum amplitudes in our final estimate of the 21 cm power spectrum (Figure 18), we show that signal loss is <2% at 95% confidence.

As shown in Table 1, the signal loss we characterize for quadratic estimation of the power spectrum band powers is tabulated along with the signal loss associated with each other potentially lossy analysis stage (see Figure 2). We correct for the signal loss in each stage by multiplying the final power spectrum results by the typical loss for each stage, except for modes within the horizon limit and immediately adjacent to the horizon limit, where we apply the maximum signal loss correction to be conservative.

When estimating our power spectra via OQEs, we generate multiple samples of the power spectrum in order to apply the bootstrap method to calculate our error bars. In detail, the power spectrum estimation scheme proposed above requires averaging at several points in the pipeline.

• 1.  
  Visibilities are averaged into five baseline groups after inverse covariance weighting (see Equation (22)).
• 2.  
  Power spectrum estimates from each of the three redundant baseline types (described in Section 2) are averaged together.
• 3.  
  Power spectrum estimates from each LST are averaged together.

With the bootstrapping technique, we do not directly perform these averages. Instead, one draws random samples within the three-dimensional parameter space specified above, with replacement, until one has as many random samples as there are total number of parameter space points. These random samples are then propagated through the power spectrum pipeline and averaged together as though they were the original data. This forms a single estimate (a "bootstrap") of $P({\boldsymbol{k}})$. Repeating random draws allows one to quantify the inherent scatter—and hence the error bars—in our estimate of $P({\boldsymbol{k}})$. When plotting ${{\rm{\Delta }}}^{2}(k)\equiv {k}^{3}P(k)/2{\pi }^{2}$ instead of $P({\boldsymbol{k}})$, we bin power falling in +k and −k, and so we additionally randomize the inclusion of positive and negative k bins.

We compute a total of 400 bootstraps. In combining independent samples for our final power spectrum estimate, we elect to use the median, rather than the mean, of the samples. One can see the behavior of both statistics in Figure 17, where we show how the absolute value of Δ²(k) integrates down as more independent samples are included in the mean and median. In this plot, one can see modes integrating down consistent with a noise-dominated power spectrum until they bottom out on a signal. In the noise-dominated regime, the mean and the median behave similarly. However, we see that the median routinely continues to integrate down as noise for longer. This is an indication that the mean is skewed by outlier modes, suggesting variations beyond thermal noise. The magnitude of the difference is also not consistent with the Rayleigh distribution expected of a cosmological power spectrum limited by cosmic variance. For a Rayleigh distribution, the median is $\mathrm{ln}2\sim 0.69$ times the mean. Instead, we interpret the discrepancy as a sign of contributions from foregrounds, which are neither isotropic nor Gaussian distributed. Since median provides better rejection of outliers in the distribution that might arise from residual foreground power, we choose to use the median statistic to combine measurements across multiple modes. As listed in Table 1, we apply a $1/\mathrm{ln}2$ correction factor to our power spectrum estimates to infer the mean from the median of a Rayleigh distribution.

Zoom In Zoom Out Reset image size

To summarize the previous section, we follow the power spectrum analysis procedure outlined in Section 5.2, we incoherently combine independent power spectrum measurements made at different times and with different baseline groups using the median statistic. As described in Section 5.4, we bootstrap over all of these independent measurements, as well as over the selection of baselines included in the power spectrum analysis for each baseline group, in order to estimate the error bars on the spherically averaged power spectrum P(k), where positive and negative ${k}_{\parallel }$ measurements are kept separate for diagnostic purposes. In the estimation of the dimensionless power spectrum ${{\rm{\Delta }}}^{2}(k)\equiv {k}^{3}P(k)/2{\pi }^{2}$, the folding of $\pm {k}_{\parallel }$ is handled along with the rest of the bootstrapping over independent modes. Finally, the measured values for P(k) and Δ²(k) are corrected for signal loss through all stages of analysis, as summarized in Table 1.

The final results are plotted in Figure 18. For the first two modes outside of the horizon where Δ²(k) is measured, we have clear detections. We attribute these to foreground leakage from inside the horizon related to the convolution kernels in Equation (8) (either from the chromaticity of the antenna response, or from the inherent spectrum of the foregrounds themselves). Somewhat more difficult to interpret are the 2.4σ excess at $k\approx 0.30h\;{\mathrm{Mpc}}^{-1}$ and the 2.9σ excess at $k\approx 0.44h\;{\mathrm{Mpc}}^{-1}$. Having two such outliers is unlikely to be chance.

Zoom In Zoom Out Reset image size

In examining the effects on the power spectrum of omitting various stages of analysis (see Figure 19), we see a pronounced excess in the green curve corresponding to the omission of crosstalk removal in fringe-rate filtering. While the signal is heavily attenuated in the filtering step, it remains a possibility that the remaining detections are associated with instrumental crosstalk. We do note, however, that the qualitative shape of the excess in the crosstalk-removed data does not appear to match that of the crosstalk-containing data.

Zoom In Zoom Out Reset image size

Another likely possibility is that the signal might be associated with foregrounds. Foregrounds, which are not generally isotropically distributed on the sky, are likely to be affected by the spatial filtering associated with fringe-rate filtering, whereas a statistically isotropic signal is not. Indeed, we see that excesses in many modes measured using the P14 stype time-domain filtering (blue in Figure 19) decrease significantly using the improved fringe-rate filter. As discussed in Parsons et al. (2015), the normalization applied to Ω_eff for fringe-rate filtering correctly compensates for the effect of this filtering on power-spectral measurements of a statistically isotropic Gaussian sky signal. We can surmise from any significant change in amplitude of the excess under fringe-rate filtering that it arises from emission that violates these assumptions. We conclude, therefore, that this excess is unlikely to be cosmic reionization, and is more likely the result of non-Gaussian foregrounds. As discussed earlier, one possible culprit is polarization leakage (Jelić et al. 2010, 2014; Moore et al. 2013), although further work will be necessary to confirm this. The interpretation of the signal as polarization leakage is, however, rather high to be consistent with recent measurements in Stokes Q presented in Moore et al. (2015), where the leakage is constrained to be <100 mK² for all k.

That the excesses at k ≈ 0.30 and $0.44h\;{\mathrm{Mpc}}^{-1}$ are relatively unaffected by the filtering could be an indication that they are more isotropically distributed, but more likely, it may mean that the simply arise closer to the center of the primary beam where they are down-weighted less. Both excesses appear to be significantly affected by omitting OMNICAL calibration (orange in Figure 19). This could be interpreted as indicating the excess is a modulation induced by frequency structure in the calibration solution. However, OMNICAL is constrained to prohibit structure common to all baselines, so a more likely interpretation is that this faint feature decorrelates without the precision of redundant calibration. To determine the nature of these particular excesses, further work will be necessary.

In order to aggregate the information presented in the power spectrum into a single upper limit, we fit a flat Δ²(k) model to measurements in the range $0.15\lt k\lt 0.5h\;{\mathrm{Mpc}}^{-1}$. We use a uniform prior of amplitudes between −5000 and 5000 mK², and assume measurement errors are Gaussian. Figure 20 shows the posterior distribution of the fit. From this distribution, we determine a mean of (18.9 mK)² and a 2σ upper limit of (22.4 mK)². The measured mean is inconsistent with zero at the 4.7σ level, indicating that we are detecting a clear power spectrum excess at $k\gt 0.15h\;{\mathrm{Mpc}}^{-1}$.

Zoom In Zoom Out Reset image size

We suspect that the excess in our measured power spectrum is likely caused by crosstalk and foregrounds. We therefore suggest ignoring the lower bound on the power spectrum amplitude as not being of relevance for the cosmological signal. On the other hand, since foreground power is necessarily positive, the 2σ upper limit of (22.4 mK)² at z = 8.4, continues to serve as a conservative upper limit. This significantly improves over the previous best upper limit of (41 mK)² at z = 7.7 reported in P14. As we show below and in greater detail in Pober et al. (2015), this limit begins to have implications for the heating of the IGM prior to the completion of reionization.

In this section, we examine the implication of the measured upper limits on 21 cm emission in Figure 18 on the spin temperature of the 21 cm line at z = 8.4. In a forthcoming paper (Pober et al. 2015), we conduct a thorough analysis of the constraints that can be put on the IGM using a simulation-based framework. As a complement to that more thorough analysis, we focus here on a simpler parameterization of the shape of the 21 cm power spectrum signal.

The brightness temperature of the 21 cm signal, $\delta {T}_{b}$, arising from the contrast between the cosmic microwave background, ${T}_{\gamma },$ and the spin temperature, ${T}_{{\rm{s}}},$ is given by

$$\begin{eqnarray}&&\delta {T}_{b}=\displaystyle \frac{{T}_{{\rm{s}}}-{T}_{\gamma }}{1+z}(1-{e}^{-\tau })\approx \displaystyle \frac{{T}_{{\rm{s}}}-{T}_{\gamma }}{1+z}\tau ,\end{eqnarray} \tag{ 28 }$$

where temperatures are implicitly a function of redshift z, and the approximation holds for low optical depth, τ. The optical depth is given by (Zaldarriaga et al. 2004)

$$\begin{eqnarray}&&\tau =\displaystyle \frac{3{c}^{3}{\hslash }{A}_{10}{n}_{{\rm{H}}\;{\rm{I}}}}{16k{\nu }_{0}^{2}{T}_{{\rm{s}}}H(z)}\end{eqnarray} \tag{ 29 }$$

where A₁₀ is the Einstein A coefficient for the 21 cm transition, n_{H i} is the density of the neutral hydrogen, H(z) is the Hubble constant, x_{H i} is the neutral fraction of hydrogen, δ is the local baryon overdensity, ν₀ is the rest frequency of the 21 cm transition, and the remainder are the usual constants. Plugging in the cosmological parameters from Planck Collaboration et al. (2015), we get

$$\begin{eqnarray}&&\delta {T}_{b}\approx {T}_{0}\;{x}_{{\rm{H}}\;{\rm{I}}}(1+\delta )\xi ,\end{eqnarray} \tag{ 30 }$$

where $\xi \equiv 1-{T}_{\gamma }/{T}_{{\rm{s}}}$ and ${T}_{0}\equiv 26.7\;\mathrm{mK}\sqrt{(1+z)/10}$.

If the spin temperature is larger than ${T}_{\gamma },$ we get the 21 cm signal in emission with respect to the CMB, and $\xi \sim 1$. However, if ${T}_{{\rm{s}}}$ is less than ${T}_{\gamma },$ δT_b is negative and ξ can potentially become large.

As in P14, we consider a "weak heating" scenario in which ${T}_{{\rm{s}}}$ is coupled to the gas temperature via the Wouthuysen-Field effect (Wouthuysen 1952; Field 1958; Hirata 2006), but little heating has taken place prior to reionization, so that ${T}_{{\rm{s}}}\lt {T}_{\gamma }$. In this scenario, because we have assumed little heating, we can approximate ξ as having negligible spatial dependence, and therefore ${T}_{0}^{2}{\xi }^{2}$ becomes a simple multiplicative scalar to the 21 cm power spectrum:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{21}^{2}(k)={T}_{0}^{2}{\xi }^{2}(z){{\rm{\Delta }}}_{i}^{2}(k),\end{eqnarray} \tag{ 31 }$$

where ${{\rm{\Delta }}}_{i}^{2}(k)$ is the dimensionless H i power spectrum.

As shown in P14, the maximum value of the prefactor in Equation (31) is given by a no-heating scenario where the spin temperature follows the kinetic gas temperature, which is held in equilibrium with the CMB via Compton scattering until z_dec ≈ 150 (Furlanetto et al. 2006) and then cools adiabatically as ${(1+z)}^{2}$. In this case, ξ is given by

$$\begin{eqnarray}&&\xi =1-\displaystyle \frac{1+{z}_{\mathrm{dec}}}{1+z}\approx -\displaystyle \frac{150}{1+z}.\end{eqnarray} \tag{ 32 }$$

At z = 8.4, this corresponds to a minimum bound on the spin temperature of ${T}_{{\rm{s}}}\gt 1.5\;{\rm{K}}$.

We can now flip this argument around and, for a measured upper bound on ${{\rm{\Delta }}}_{21}^{2}(k)$, we can use models for ${{\rm{\Delta }}}_{i}^{2}(k)$ in Equation (31) to place a bound on ${T}_{{\rm{s}}}$. We consider a class of "patchy" reionization models (P12a; P14) which approximates the ionization power spectrum as flat between minimum and maximum bubble sizes, ${k}_{\mathrm{min}}$ and ${k}_{\mathrm{max}},$ respectively:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}^{2}(k)=({x}_{{\rm{H}}\;{\rm{I}}}-{x}_{{\rm{H}}\;{\rm{I}}}^{2})/\mathrm{ln}({k}_{\mathrm{max}}/{k}_{\mathrm{min}}).\end{eqnarray} \tag{ 33 }$$

For combinations of ${k}_{\mathrm{min}}$ and ${k}_{\mathrm{max}},$ we determine the minimum spin temperature implied by the 2σ 21 cm power spectrum upper limits shown in Figure 18. Figure 21 shows the results of these bounds for neutral fractions of x_{H i} = 0.1, 0.3, 0.5, 0.7, and 0.9. In almost all cases (excepting x_{H i} = 0.1, 0.9 for ${k}_{\mathrm{min}}\lt 0.1h\;{\mathrm{Mpc}}^{-1}$), we find that ${T}_{{\rm{s}}}\gtrsim 3\;{\rm{K}}$, indicating that our measurements are inconsistent with the spin temperature being coupled to a kinetic temperature governed strictly by adiabatic expansion.

Zoom In Zoom Out Reset image size

Our results become more interesting in the range of ${k}_{\mathrm{min}}\sim 0.1$ and ${k}_{\mathrm{max}}\sim 30$ representative of fiducial simulations (Zahn et al. 2007; Lidz et al. 2008). For neutral fractions of 0.3, 0.5, and 0.7, we find that ${T}_{{\rm{s}}}\gtrsim 4\;{\rm{K}}$. Pober et al. (2015) improves on these results by using a simulation-based framework, rather than relying on coarse parametrizations of the power spectrum shape. They compare the limits they find to the amount of heating possible given the currently observed star formation rates in high-redshift galaxy populations (Bouwens et al. 2014; McLeod et al. 2015) and assumptions about the relationship between star formation rates and X-ray luminosities (Furlanetto et al. 2006; Pritchard & Loeb 2008; Fialkov et al. 2014). Assuming the midpoint of reionization lies close to z = 8.4 (a reasonable assumption given that Planck Collaboration et al. 2015 suggests a midpoint of z = 8.8), both the bounds found in this paper and Pober et al. (2015) show evidence for heating that places constraints on the possible values for the star formation rate/X-ray luminosity correlation given certain models of the star formation rate density redshift evolution. We refer the reader to Pober et al. (2015) for a detailed examination of these results.