NEW LIMITS ON 21 cm EPOCH OF REIONIZATION FROM PAPER-32 CONSISTENT WITH AN X-RAY HEATED INTERGALACTIC MEDIUM AT z = 7.7

In preprocessing, we use delay/delay-rate (DDR) filters (Parsons & Backer 2009) to help identify radio-frequency interference (RFI) events and as part of a data compression technique that reduces data volume by over a factor of 40. Further details are supplied in Appendix A. The result is that the 2048 original frequency channels become 203, and the 60 original time samples per 10 minute file become 14, for an effective integration time of 43 s. Our chosen output integration time is slightly shorter than the 45.7 s interval dictated by the maximum fringe rate of a 300 m baseline. These data contain all emission that rotates with the sky and falls in the delay range |τ| < 1000 ns, corresponding to −0.5 < k_∥ < 0.5 h Mpc⁻¹ at z = 7.7 (164 MHz) but compressed in data volume by over a factor of 40.

Phase calibration is achieved using standard interferometric self-calibration with the addition of constraints derived from the substantial instantaneous redundancy of the array (Wieringa 1992; Liu et al. 2010). The GPS surveying and total-station laser ranging located all antennas to within ±10 cm from a perfect grid, based on agreement between the two surveying methods. These antenna coordinates are used to calibrate the electrical delay of each antenna signal path on the basis of redundancy. Because of the high degree of instantaneous redundancy in this array, we are able to derive per-antenna delays for each signal path relative to only two parameters that could not be solved for on the basis of redundancy: the relative electrical delay between the antennas in fiducial east–west and north–south baselines (0–16 and 0–1, respectively). We find per-antenna delays averaged over each day of observation to agree within ±0.1 ns between days. Because this residual delay is a factor of 1000 smaller than the 100 ns width of the delay bins used in the power spectrum analysis in Section 3.6, the signal attenuation that would result from averaging over the phase variation associated with this residual delay is very small (of order 10⁻⁶). On this basis, we adopt a single per-antenna delay solution for the entire 92 day observation. These per-antenna delay solutions for all baselines are incorporated into a global self-calibration step that derives the remaining two unknown phase parameters using the central components of Centaurus A, Fornax A, and Pictor A.

The use of redundancy in calibration allows antenna-based gain and phase parameters to be derived from the ratios of baselines within instantaneously redundant sets. Because these instantaneously redundant baselines measure the same sky signal, the signal with each set of redundant baselines essentially becomes a single degree of freedom that, by taking ratios, is projected out prior to calibration. Hence, there is no possibility of signal loss associated with the redundant aspect of calibration.

3.3.1. An Absolute Flux Scale from Pictor A

The first challenge we encounter in flux calibration is establishing an accurate spectrum of a calibration reference. Pictor A is chosen as a calibrator source of necessity; given the grating lobes of the synthesized beam generated by our maximum-redundancy array configuration, side lobes of Pictor A and Fornax A dominate nearly all beam pointings. Of these two sources, Pictor A is the least resolved.

Unfortunately, the spectrum of Pictor A is particularly poorly characterized in the 100–200 MHz band where PAPER operates. Previous measurements imply that the spectrum below 400 MHz deviates substantially from the power law (with spectral index α = −0.85) that holds at higher frequencies (Perley et al. 1997). There are also serious inconsistencies in the measured flux of Pictor A below 1 GHz (e.g., Slee 1995; Kuehr et al. 1981; Large et al. 1981; Burgess & Hunstead 2006) and a dearth of measurements at nearby frequencies in order to establish the spectral slope of Pictor A in this band.

This shortcoming in the literature has recently been addressed with observations from a 64 antenna, single-polarization deployment of the PAPER array in a minimum-redundancy configuration more suited for imaging. The details of these observations and their calibration are given in Jacobs et al. (2013), along with the resulting spectrum of Pictor A that is shown to be consistent with a single power law that is best fit by

$$\begin{equation} S_\nu = S_{150} \left(\frac{\nu }{150\, {\rm MHz}}\right)^\alpha \end{equation} \tag{ 1 }$$

with S₁₅₀ = 382 Jy ± 5.4, with α = −0.76 ± 0.01. Error ranges indicate 76% confidence intervals.

3.3.2. Gain Calibration Augmented by Redundancy

For these maximum-redundancy observations, gain calibration is derived on the basis of both redundancy and self-calibration. First, we correct for the known dependence of amplifier gain and cable attenuation on ambient temperature using measured temperature data and temperature coefficients that were characterized in laboratory measurements (Parashare & Bradley 2009) and confirmed in field measurements (Pober et al. 2012). These corrections correspond to approximately a 5% variation in amplitudes over the course of the observations. This is the only time-dependent parameter used in the gain calibration.

Next, redundant self-calibration is used to derive the amplitude of each antenna signal path relative to a single overall bandpass function for each (XX,YY) polarization. As mentioned in the context of phase calibration, per-antenna gains are derived from the ratios of baselines within instantaneously redundant sets, and because the sky signal cancels in these ratios, there is no possibility of signal loss associated with redundant gain calibration.

This overall bandpass is then calibrated to the flux scale established in Jacobs et al. (2013) by phasing to the position of Pictor A and averaging over all intercolumn baselines. Because Pictor A is resolved at the 20% level on the longest (300 m) east–west PAPER baselines, we use per-baseline estimates of resolution weights to down-weight and correct for resolution effects. These resolution weights are estimated from the 90 cm map of Pictor A presented in Perley et al. (1997). Averaged over all of the baselines in the array, the inclusion of resolution effects change the measured flux by less than 5%. As in Jacobs et al. (2013), a fringe-rate filter is applied to the resulting beam to suppress side lobes that introduce variations deviating more than ±0.1 mHz from the fringe rate of the source in question.

To derive a source spectrum from a drift-scan source profile, we average time samples using weights derived from a model of PAPER's primary beam response as an effective approximation of inverse-variance weighting (Pober et al. 2012) and average over 92 days of drift-scan observations. A ninth-order polynomial is fit over a 100 MHz bandwidth to bring each measured spectrum in accordance with the unpolarized model of the Pictor A spectrum described in Section 3.3.1. The resulting Pictor A spectrum, plotted for each polarization in Figure 3 along with the reference spectrum derived from minimum-redundancy PAPER observations (Jacobs et al. 2013), characterizes the flux scale used in the subsequent power-spectrum analysis.

Zoom In Zoom Out Reset image size

This polynomial bandpass model uses a limited number of terms to prevent the coupling of beam side lobes, which are the most likely cause of the residual gain variation in Figure 3, into the overall gain calibration. While such side lobes are difficult to distinguish from intrinsic frequency structure, we err on the side of caution (using a limited number of degrees of freedom) to avoid introducing any frequency-dependent gain terms not in the data. In this cautious approach, it is possible that higher-order frequency structure in the instrumental bandpass may not be calibrated out. To address this possibility, we make the argument that evidence of such structure must be present in the measured Pictor A spectrum after the application of the bandpass calibration, and to the extent that some higher-order structure is demonstrably present, this structure can only be ascribed to uncalibrated structure in the instrumental bandpass if it modulates foreground structure to higher-order modes in the measured power spectrum. The fact that, at the upper limits we report in Section 4.2, we do not see this structure argues for its absence. In any case, this conservative modeling of the PAPER bandpass with as few degrees of freedom as necessary does not compromise the validity of the upper limits that are reported.

The restricted order of the polynomial used in calibrating the bandpass also serves to insulate against the signal loss that would result from inadvertently including 21 cm EoR fluctuations in the bandpass model. We use simulations with many realizations of low-level white noise added to the polynomial bandpass model we fit using Pictor A to quantify the expected signal loss as a function of k_∥ in the power spectra measured in Section 4.2. In each case, we compare the power spectrum derived from the original injected noise with the assumed polynomial bandpass model to the power spectrum of the injected noise after dividing by a new fit of a ninth-order polynomial to the bandpass. We find signal loss in the range −0.06 < k < 0.06 h Mpc⁻¹ to be less than 8% and signal loss outside of this range to be below 0.2% and to decrease rapidly with |k_∥|. This simulation calculates signal loss strictly on the basis of fitting a ninth-order polynomial bandpass model to a spectrum containing the 21 cm EoR signal and ignores the use of multiple baseline lengths and Earth-rotation synthesis in measuring the Pictor A spectrum. This averaging of independent modes attenuates the coupling of the 21 cm EoR signal into the Pictor A spectrum that is used in the bandpass calibration in proportion to the square root of the number of modes averaged. Hence, this simulation establishes an upper bound on signal loss that results from bandpass calibration and shows it to be a small effect for the k_∥ modes that we use to set our upper limits in Section 4.2. Because polynomial fitting is a linear process, the fractional signal loss we calculate using a white-noise signal model is independent of the amplitude of the model signal.

Following calibration to a consistent flux scale, XX and YY polarizations for each baseline are directly summed to form a coarse estimate of the Stokes I parameter. This simplistic construction of Stokes I neglects differences in the beam responses of each linear polarization and, as a result, contains contributions from Stokes Q away from beam center (Jelić et al. 2010; Moore et al. 2013). Nonetheless, for analysis of a single baseline type, correcting for direction-dependent gains is highly nontrivial at best, and we attempt no such correction here. On the basis of the asymmetry of the PAPER beam model under 90° rotation, we estimate that this construction of Stokes I will contain approximately 4% of Stokes Q as well (Moore et al. 2013). For a more thorough treatment of polarization in the context of drift-scanning arrays with delay-spectrum analysis, we refer the reader to Moore et al. (2013).

The next major step in our data reduction pipeline is to remove a best estimate of foreground emission arising from the galactic synchrotron and extragalactic point sources. This is an essential step, as foreground emission exceeds the expected level of EoR emission by nine orders of magnitude (Pober et al. 2013a) and can conceal low-level RFI and crosstalk systematics. Moreover, as described in P12b, smooth-spectrum foreground emission can corrupt higher-order k modes in power-spectrum measurements as a result of side lobes that arise from RFI flagging and the finite (windowed) bandwidth used in the line-of-sight Fourier transform. Although power-spectrum measurements that aim to constrain 21 cm reionization must be limited to ∼10 MHz to avoid significant signal evolution within the band, foreground emission is generally coherent over the entire 100 MHz band of PAPER observations.

A delay-domain filter applied over a wide bandwidth can achieve a much sharper degree of foreground isolation than is possible in any sub-band, and this can be used to effectively model and remove a smooth foreground model from each sub-band of interest. Following the procedure outlined in P12b and implemented in Pober et al. (2013a), we perform a delay transform of each integration over the entire observing band, according to the equation

$$\begin{equation} \tilde{V}(\tau)=\int {W(\nu) S(\nu) V(\nu)e^{-2\pi i\nu \tau } d\nu }, \end{equation} \tag{ 2 }$$

where V(ν) is the measured visibility as a function of spectral frequency ν, τ is delay (the Fourier complement to ν), W(ν) is a Blackman–Harris windowing function (Harris 1978) used for minimizing side lobes and band-edge effects, S(ν) encodes the frequency-dependent sample weights that result from RFI flagging, and ∼ denotes the delay transformation. Following the procedure outlined in Parsons & Backer (2009), we deconvolve the delay-domain kernel that results from the product W(ν)S(ν) using a CLEAN-likeiterative forward modeling deconvolution algorithm (Högbom 1974) in the delay domain, restricting model components to fall inside of 15 ns beyond the horizon limit of each baseline. We chose a 15 ns standoff from the horizon to suppress the first level of suprahorizon emission that results from the inherent width of $\tilde{S}(\tau)$. Because the sampling function is known exactly, the restriction of model components to low-order delay modes ensures that the deconvolution process cannot introduce power at higher-order delay modes beyond what was scattered there by side lobes of the sampling function. Moreover, because of the least-square metric used in evaluating CLEAN models for divergence, any delay-domain residuals associated with the sampling function will be unbiased.

The resultant CLEAN model is subtracted from $\tilde{V}(\tau)$ to effect a delay-domain filter suppressing smooth-spectrum emission. On the 30 m, nearly east–west baselines used in this analysis, in a band centered at 164 MHz, this corresponds to a cutoff at k = ±0.06 h Mpc⁻¹. This delay filter is solely responsible for −60 dB of suppression (in mK² power) of the foreground signal seen between the unfiltered signal (top black curve) and the final residual (lowest black curve) shown in Figure 4. More importantly, in the narrow-band delay transformation that forms the basis of power-spectrum measurements in Section 3.6, this filter heavily suppresses the coupling of delay modes beyond the horizon to the foreground-dominated modes within the horizon. There is a small amount of signal loss associated with this filtering. Using the simulation methodology described in Section 3.7, we compute this signal loss to be 4.8% for the two modes on either side of the horizon limit, 1.3% for the next modes outside of that, and <0.015% elsewhere.¹⁷ The cost of this signal loss is far outweighed by the benefits of suppressing couplings to foreground-dominated modes where, despite the low coupling strength, the sheer magnitude of the foreground signal creates a dominant bias.

Zoom In Zoom Out Reset image size

After the removal of the bulk of the foregrounds, the signal remaining for each baseline is visibly noise-dominated. At this point, we apply a final pass of RFI excision, flagging 3σ outliers, and proceed to crosstalk removal to remove systematics that were previously undetected beneath the bright, smooth-spectrum emission. Crosstalk removal proceeds by subtracting the one hr time average of the visibilities for each baseline from each integration. This process, described in Parsons et al. (2010), distinguishes oscillating fringes associated with celestial emission from the static phase bias associated with crosstalk. This crosstalk-removal filter essentially constitutes a high-pass fringe-rate filter, as described in Section 3.5. The width of the stop-band of the crosstalk filter is much narrower than low-pass fringe-rate filters described in that section.

Night time data are then averaged in LST over the 92 day observation using 43 s time bins that match the integration interval of the data after the compression step described in Section 3.1. In the LST binning process, sporadic non-Gaussian events that evade the RFI flagging process can easily skew the average value in an LST bin away from the median value. To avoid bias from non-Gaussian tails, we flag the 10% of data in each LST bin, for each frequency channel, that deviate most from the median value, and exclude them from the computed average for that LST bin. This last flagging step suppresses non-Gaussian outlying events without biasing the computed average in an LST bin. We simulate this median filtering process for samples in a normal distribution around an underlying random signal assigned to each LST bin and verify that the computed variance of the underlying signal is not biased by the median filter, regardless of the level of the signal with respect to the noise. We also use this simulation to verify that bootstrapping correctly reproduces the range of variation in the final estimator as expected for Gaussian statistics and that the bootstrap error bars are unaffected by the median filter. These simulations indicate that the LST binning process described above has no associated bias or signal loss for recovering the cosmological EoR signal.

In the last step prior to forming power spectra, we apply a fringe-rate filter to effect time-domain integration, using the effective time interval that a baseline measures a single k mode to integrate coherently (with noise decreasing as $\sqrt{t}$, in units of mK), before measurements at different times represent independent modes that must be squared before further integration (with noise now decreasing as $\sqrt{t}$, in units of mK²).

One way of handling this additional integration is via gridding in the uv plane. Each measured visibility in a wide-field interferometer represents the integral over a kernel in the uv plane that reflects the primary beam of the elements (Bhatnagar et al. 2008; Morales & Matejek 2009) and the w component of the baseline (Cornwell et al. 2003). As noted in Sullivan et al. (2012) and Morales & Matejek (2009), in order to optimally account for the mode-mixing introduced by these kernels, gridding kernels must be used that correctly distribute each measurement among the sampled uv modes, such that, in the ensemble average over many measurements by many baselines, each uv mode becomes an optimally weighted estimator of the actual value given the set of measurements.

However, this approach has a major shortcoming when applied to maximum-redundancy array configurations. In order to maximize sensitivity, such configurations are set up to deliberately sample identical visibilities that reflect the same combinations of modes in the uv plane with few nearby measurements (P12a). As a result, such array configurations tend to lack enough measurements of different combinations of uv modes to permit the ensemble average to converge on the true value. Said differently, maximum-redundancy array configurations tend to produce measurement sets that, when expressed as linear combinations of uv modes of interest, are singular.

Our alternative approach avoids this and many of the difficulties outlined in Hazelton et al. (2013) by applying a carefully tailored fringe-rate filter to each time series of visibility spectra. As outlined in Appendix A in the context of data compression, we take the Fourier transform of the time series in each channel and apply a low-pass filter that preserves fringe rates that geometrically correspond to sources rotating on the celestial sphere. For a planar array with transit observations, fringe rates vary according to declination, with fringe rates reaching a maximum (f_max) at δ = 0°, decreasing to 0 at δ = −90°, and for an array such as PAPER deployed near −30°S latitude, reaching a minimum of ≈ − f_max/2 at δ = −60° on the far side of the south celestial pole. In order to avoid introducing undesirable frequency structure, we apply the same filter, tuned to the width set by the highest frequency of the sub-band used in the power spectrum analysis described in Section 3.6, to each channel, even though maximum fringe rates are generally frequency-dependent. In a future paper, we will explore the idea of employing fringe-rate filters that purposely down-weight fringe-rate modes on the sky according to the expected signal-to-noise ratio in each mode. Such filters would essentially correspond to a one-dimensional implementation of the inverse primary beam uv gridding discussed in Morales & Matejek (2009), and have many features in common with m-mode synthesis described in Shaw et al. (2014).

Since thermal noise scatters equally into all fringe rate bins, applying a filter that passes only fringe rates corresponding to the celestial emission has the effect of denoising the data. We apply such a filter to the data, choosing the bounds of the filter to match the geometric bounds set by a 30 m east–west baseline, according to the equation

$$\begin{equation} f_{\rm max} = \frac{|{\bf b}_{\rm eq}|}{c} \omega _\oplus \nu, \end{equation} \tag{ 3 }$$

where f_max is the maximum fringe rate, b_eq is the baseline vector projected parallel to the equatorial plane, c is the speed of light, ω_⊕ is the angular frequency of the Earth's rotation, and ν is the spectral frequency. At 174 MHz (the highest frequency in a 20 MHz window centered on 164 MHz that is used in Section 3.6), f_max = 1.27 mHz, corresponding to a fringe period of 788 s. Hence, the fringe-rate filter that is applied passes fringe rates in the range −0.63 < f < 1.27 mHz. The width of this filter corresponds in sensitivity to an effective integration time of 525 s. We note that this filtering could have been applied during the data compression described in Section 3.1 but was implemented separately to enable the compression to work uniformly over all baselines in the array without additional information about antenna location.

After applying this filter, we transform the data back to time domain in preparation for forming power spectra via the delay transform. It should be noted that, in time domain, the data are now heavily oversampled; adjacent samples are no longer statistically independent. Hence, when averaging power spectra versus time, noise will not beat down according to the strict number of samples but rather according to the actual number of statistically independent samples underlying the time series.

In the final stage of our core analysis, we apply the delay transform in Equation (2) to a 20 MHz band centered at 164 MHz, again using a Blackmann–Harris window that yields an effective (noise-equivalent) bandwidth of 10 MHz, centered at redshift z = 7.7. In this paper, we limit our analysis to a single cosmological bin where, out of the entire PAPER observing bandwidth, the noise and foreground residual in Figure 4 is lowest. As discussed in P12b, the τ modes that result from the delay transform have a peaked response in k_∥ that allows each (u, v, τ) mode to be interpreted as sampling the mode k = 2π(u/X, v/X, τ/Y), with X and Y representing conversion factors from angle and frequency to comoving distance, respectively. As such, we use k_∥ and τ modes interchangeably in the following discussion.

In contrast to the wide-band delay filtering described in Section 3.4, where a CLEAN-like deconvolution with a restricted range of model components was used to suppress the side lobes of smooth-spectrum emission resulting from unsmooth sampling functions, we do not deconvolve covariances in $\tilde{V}(\tau)$ introduced by windowing and sampling functions. First, in the band chosen for this analysis, the sampling function S(ν) in Equation (2) is nearly unity at all frequencies, thanks to the exquisite RFI environment at the Karoo site in South Africa. Secondly, the covariances between nearby τ modes introduced by the windowing function go hand-in-hand with the suppression of covariances between modes at larger separations. Hence, we defer the removal of covariances introduced by the windowing and sampling functions at this stage until the covariance diagonalization process described in Section 3.7.

From these delay-transformed visibilities, we construct unbiased estimators of the power spectrum, $\widehat{P}({\bf k})$, by cross-multiplying delay spectra between the seven baseline groups constructed in Section 3.4 and summing over all cross-multiples according to the equation

$$\begin{equation} \widehat{P}({\bf k}_{t\tau }) = \left(\frac{\lambda ^2}{2k_{\rm B}}\right)^2\frac{X^2Y}{\Omega B} \left\langle {\tilde{V}_i(\tau,t) \tilde{V}_j^*(\tau,t)}\right\rangle _{i<j}, \end{equation} \tag{ 4 }$$

which follows from Equation (12) of P12a, with λ being the observing wavelength, k_B is Boltzmann's constant, X²Y is a cosmological scalar with units of (h⁻³ Mpc³/sr Hz), and Ω is the angular area,¹⁸B is the bandwidth, 〈...〉_{i < j} indicates the ensemble average over instantaneously redundant baseline measurements indexed by i, j, and $\tilde{V}(\tau,t)$ is the delay-transformed visibility, expressed in terms of delay τ and time t. We use t as a subscript on k to denote the different modes sampled by a baseline as the sky rotates and τ to indicate the dependence of k on the delay mode in question.

Equation (4) represents the diagonal-covariance limit of an optimal quadratic estimator (Dillon et al. 2013; Liu & Tegmark 2011), which we review in detail and adapt to delay-transformed visibilities in Appendix C. The diagonal limit generates optimal power spectrum estimators in the simple case of ideal signal covariance among statistically independent τ modes. In theory, the translation-invariance of the cosmological signal (along with the interchangeability of τ and k) ensures that the signal covariance is diagonal. Over narrow frequency bands such as the one used for these measurements, a similar invariance (and thus, diagonality) along the frequency direction should hold for the noise and foregrounds (Dillon et al. 2013). In practice, we encounter substantial systematics whose covariance deviates from this ideal. In Section 3.7, we examine how off-diagonal covariances arising from systematics can be suppressed in a straightforward way and incorporated into Equation (4) to produce power-spectrum estimators that are, to leading order, optimal. Because this later analysis follows the same trajectory as the analysis presented, we will first bring the current analysis to completion before examining how it is modified to incorporate off-diagonal covariances.

For N redundant samples of a k mode, N(N − 1)/2 cross-multiples are constructed. At this stage, we do not use i > j pairings, choosing instead to preserve the imaginary component of the power spectrum estimator, even though P₂₁(k) is real valued, as a diagnostic. In the final power-spectrum estimate, we drop the imaginary component, which is equivalent to including the i > j baseline pairings. Hence, we ultimately have N² − N measurements of P(k), with noise decreasing proportional to N in mK² units, as would be expected for coherently integrating redundant samples before squaring them.

Next, we average over measurements of independent k modes that statistically reflect the same underlying power spectrum, to produce our best power-spectrum estimate,

$$\begin{equation} \widehat{\Delta }^2_{21}(k) = \frac{k^3}{2\pi ^2} \langle \widehat{P}({\bf k}_{t\tau })\rangle _{|{\bf k}_{t\tau }|=k}, \end{equation} \tag{ 5 }$$

where the three-dimensional (3D) symmetry of the power spectrum is invoked to average over all independent measurements of modes in a shell of |k| = k, with independent measurements indexed here by t. As described in Section 3.5, the number of independent modes that are averaged (with noise decreasing with number of modes, M, as $\sqrt{M}$ in mK² units; see P12a) is determined by overall observing window and the number of fringe-rate bins that are preserved in the fringe-rate filtering process. Noise levels are estimated using bootstrapping, as described in Section 3.8. Since we have not decimated the number of integrations to the critical sampling rate corresponding to the width of the applied fringe-rate filter, M is not the number of integrations. However, we are free to average the power-spectrum estimates for each integration, even though nearby samples do not have statistically independent noise, understanding that noise will decrease according to the number of underlying independent samples.

We also apply this same process for constructing power-spectrum estimators to data that include smooth-spectrum foreground emission before it was suppressed with the wide-band delay filter described in Section 3.4. While these data are not expected to yield accurate measurements of P(k) at higher k owing to side lobes of data flagging and low-level systematics, as mentioned previously, they do yield effective measurements of the bright foreground emission that falls within the limits of the wide-band delay filter (top right panel, Figure 5). At the baseline length used in this analysis, this emission is dominated by bright point sources such as Pictor A and Fornax A (Parsons et al. 2010; Jacobs et al. 2011). In order to present a complete picture of the power spectrum of emission on the sky that is valid at all k modes, we stitch together the final power spectrum using measurements derived from foreground data at k modes where emission is suppressed by the wide-band delay filter, and using measurements from the foreground-suppressed data elsewhere.

Zoom In Zoom Out Reset image size

The result of the analysis up to this point is the power spectrum illustrated by the cyan curve in Figure 6, which exhibits a residual off-diagonal covariance structure that is illustrated in the top right panel of Figure 5.

Zoom In Zoom Out Reset image size

The off-diagonal covariances shown in the upper-right panel of Figure 5, which are estimated empirically from measured data, are indicative of systematics corrupting what should nominally be nearly statistically independent estimators of k modes of reionization. This diagnosis is reinforced by the substantial variation of the covariance structure between different baselines; we expect the sky signal, whether from foregrounds or reionization, to be common to all baseline cross-products (as in the top left panel of Figure 5). Hence, we now turn to investigating techniques for removing these systematics while retaining, to the greatest extent possible, the reionization signal we are after. As we describe briefly here and in detail in Appendix C, this can be achieved by using the quadratic estimator formalism (Dillon et al. 2013; Liu & Tegmark 2011).

As derived in Appendix C, the special case of applying optimal quadratic estimators to delay-domain visibilities can, to leading order, be approximated as a modification of the visibility. In brief, knowledge of the covariances allows off-diagonal leakages between delay bins to be predicted and removed. If we take the output of the off-diagonal covariance removal process described in Equation (C20) as a perturbed visibility, $\tilde{V}^\prime (\tau)$ from Equation (C21), such that then for this modified visibility, Equation (4) encodes the remainder of the optimal estimator formalism.

There are, however two subtle issues to address in applying the off-diagonal covariance removal process described in Appendix C. The first problem pertains to our limited ability to empirically deduce the true covariances present in our data. Because we have a limited number of independent samples from which to estimate covariances in the data, we measure residual noise and signal covariances that, in an infinite sample set, would not exist. If we naively subtract these residual covariances, we overfit the noise and suppress the 21 cm signal.

There is a simple prescription for avoiding this problem, given that, as the top right panel of Figure 5 illustrates, the dominant systematics we see in the covariance of k modes exhibit a baseline-dependent structure. Because celestial emission, whether from foregrounds (top left panel, Figure 5) or 21 cm reionization, should have identical structure among redundant baselines, we may subtract the average covariance among many baseline cross-products (the "panels" in each subplot of Figure 5) from each individual cross-product to remove the residual signal covariance and leave behind only the covariance of noise and systematics from which off-diagonal covariances can be removed even if they are residuals that would not be present in the true covariance matrix without attenuating the desired celestial signal. This lossless removal of systematics is another substantial benefit of the highly redundant configurations presented in P12a.

We apply this technique for each baseline cross-product, using all baselines except the two being cross-multiplied for estimating the average covariance in order to avoid, to the greatest extent possible, coupling baseline-dependent systematics into each estimate of the common-mode covariance. Iterating this process a modest number of times (two or three), produces an improvement in the results as initial baseline-specific systematics are first removed from the baseline, allowing an improved estimate of the average covariance to be subtracted, such that the remainder of the baseline-specific systematics are removed. Further iteration is not necessary, as the process rapidly converges to the result shown in the bottom-left panel of Figure 5.

This process, though quite effective, introduces a second issue that must be addressed. In Equation (4), we were careful to exclude products between measurements involving the same baseline in order to avoid incurring the noise bias that would result. The benefits of avoiding this bias far outweigh the slight improvement in sensitivity that including such "auto-products" would produce. Unfortunately, by using the average of many baselines to estimate and subtract an average covariance from each baseline cross-product and then subtracting off-diagonal terms in the residual, we couple the baseline-averaged noise into the data from each baseline. The result is a low-level residual noise bias approximately equal to the noise in the power-spectrum estimate.

The most straightforward approach we have found for eliminating this residual noise bias, other than direct subtraction (which introduces additional complications), is to divide baselines into four separate groups of approximately equal size. Within each group, we apply the off-diagonal covariance removal process, including the subtraction of an average covariance, using data only from baselines within that group. Then, to avoid incurring a noise bias, we use only cross-products of baselines between these four groups to estimate P(k). By excluding intragroup cross-products, this approach sacrifices a factor of approximately 15% in sensitivity (in mK²), but as before, we find the benefits of avoiding noise bias to outweigh the loss in sensitivity.

The last step in the process of suppressing off-diagonal covariances in our data is, for select modes, to relax the constraint of not subtracting a baseline-averaged covariance for selected modes, even if it results in overfitting the noise and the signal suppression that is associated with it. Particularly, we note the interior seven k modes (the five inside of the horizon limit shown in Figure 6 as well as the first modes beyond this limit on either side) that we measure are more than an order of magnitude brighter than other modes and are so corrupted by smooth-spectrum foreground emission that, barring a heroic effort aimed at modeling and removing these foregrounds, they are unlikely to be useful for constraining high-redshift 21 cm emission. We find it advantageous to remove all off-diagonal covariances associated with these modes, even if it means overfitting the noise. The result of this process is shown in the bottom right panel of Figure 5.

Because we have overfit the noise in this final step, it now becomes necessary to investigate how we have affected the 21 cm signal that we aim to measure because failure to account for signal attenuation can lead to erroneous constraints. We use Monte Carlo simulations (see, e.g., Masui et al. 2013; Paciga et al. 2013; and Switzer et al. 2013) to estimate the expected signal attenuation through the analysis chain from Section 3.6 onward, as illustrated by the red data-flow path in Figure 2. Because our analysis does not model out modes of arbitrary shape (such as occurs in principal component analysis), and because delay filtering and covariance removal operates in the same space in which the power spectrum is measured, the simulations needed to characterize signal loss for our analysis are relatively straightforward. In our simulations, we apply the analysis pipeline used on the data, including off-diagonal covariance suppression, to a battery of simulated data sets that contain both a randomly generated sky signal and thermal noise. In each data set, the model sky signal consists of a spatially random, flat-spectrum signal common to all baselines that rotates with the celestial sphere, modeling the output of the fringe-rate filtering step that precedes the narrow-band delay transform in Figure 2. The noise that is added to this sky model is a Gaussian random signal unique to each baseline, to which a fringe-rate filter is applied.¹⁹

We examine two cases in our simulations. In the first, we apply the off-diagonal covariance removal process using empirically determined covariances internal to the simulated data. This case effectively treats the simulated signal as if it were the data and investigates the effects of overfitting noise in the data used to produce the covariance matrices. In the second case, we apply the exact diagonalization matrices used to remove covariance terms in the measured data to the simulated signal. This case examines the effect that the operations applied to the data would have on a simulated signal that is uncorrelated with the data, as would be the case for a faint signal that is subdominant to foregrounds and instrumental systematics in the data. These two cases serve to characterize any combination of signals that are correlated or uncorrelated with the data, and hence they accurately bound the signal loss to which a potential 21 cm reionization signal would be subject in this analysis.

In both simulation cases, the amplitude of each k mode in the output power spectrum of the recovered signal is compared to the known input amplitude of the simulated signal, averaged over 1000 independent simulations. In the first case, results indicate that removing the off-diagonal covariances associated with the central seven modes results in a ∼5% reduction in signal amplitude (in mK²) for other modes and a ∼30% reduction for the seven modes in question. We compensate by dividing by these signal attenuations when reporting the power spectrum limits in Figure 6. In the second case, which tests the effects of our exact analysis on the noise statistics of an uncorrelated signal, we find that noise that is uncorrelated with the data is not significantly affected by this analysis.

We note that none of the analysis described above constitutes formal model subtraction. Because we estimate covariances from data, it is possible to overfit noise and subtract signal, but we take steps to avoid this and to compensate for signal attenuation where subtraction does happen by estimating signal attenuation via Monte Carlo simulations. If the full covariances were known a priori, this analysis would be strictly lossless (Tegmark 1997) and, as shown in Appendix C, optimal to leading order.

One consequence of the process for removing off-diagonal covariances described above is that baseline-specific covariances, including residual covariances from thermal noise, are heavily suppressed. As a result, it is not possible to estimate noise residuals simply from the variation between baseline cross-products at the conclusion of this process. Instead, we use bootstrapping of baseline samples at the input of the off-diagonal covariance removal process to estimate error in the output power-spectrum estimate.

There are, however, a few subtleties in applying bootstrap resampling in this context. One complication is that, as described in Section 3.7, in order to avoid incurring a noise bias, it is vital that the removal of off-diagonal covariances proceed on independent groups of baselines. To reflect the fact that our power spectrum estimator is unbiased, it is important to restrict the resampling with replacement to avoid including data from the same baseline into the four independent baseline groups. This is done by first randomly assigning baselines uniquely to one of the four groups, and then applying sampling with replacement only within each of these groups. At each iteration of resampling, the assignment of baselines to groups is randomized.

Another issue in applying bootstrapping to this case pertains to how a second bias can be mistakenly included if we are not careful to keep track of repeated entries of the same baseline data within a group. In this case, the problem results from averaging used to suppress residual signal covariances in the case that such cross-products may actually be the product of a baseline's data with itself. The noise bias in these auto-products is coupled into the average of the covariances and biases the result of the off-diagonal covariance removal in each baseline group. In the cross-multiplication of baselines between groups, these biases (which are positive) couple into the result.

The solution in this case is to exclude auto-products from the determination of the average covariance among baseline cross-products. Since only a few independent cross-products are necessary to distinguish baseline-specific features from a residual signal covariance, this restriction is relatively harmless. However, it does require that there be at least two independent baselines in each resampled baseline group. By imposing this restriction, as well as the restriction to not repeat baselines between two independent groups, we limit the sample space that bootstrapping explores. It is therefore important to check that there remain enough data permutations to adequately estimate errors. Even with these restrictions on the resampling of 28 redundant baselines, the number of distinct permutations is much larger than the number of independent samples we have, so the limiting factor in estimating errors will be the number of samples, not the number of bootstrap resamplings.

The mean and 2σ errors shown in Figure 6 are derived from 100 bootstrap resamplings, with error bars calculated to enclose 95% of the bootstrap samples. The errors we report do not include covariances between k modes, which are not treated in bootstrapping. However, as implied in Dillon et al. (2014), if these modes are not added together in binning, then the variances derived from bootstrapping are sufficient to characterize the errors.

Figure 7 shows images of a single frequency channel using Earth-rotation synthesis, including all 496 baselines of PAPER's 32 antenna maximum-redundancy array. As a reminder, the power spectra we report in Section 4.2 is estimated directly from a restricted set of visibilities. These images are for diagnostic purposes. As indicated by the marked reduction of emission between the upper right and lower left panels, the application of the delay filter suppresses peak emission by at least three orders of magnitude. The further decrease in residual emission with the application of the fringe-rate filter between the lower left and lower right panels indicates that residuals in the lower left panel are noise-dominated, with the lower right panel showing that peak emission from the raw image has been suppressed by 3.3 orders of magnitude.

Zoom In Zoom Out Reset image size

This result is paralleled in Figure 4, which shows the rms temperature of visibilities as a function of frequency in LST bins between 2:00 and 5:00, before and after the application of the wide-band delay filter, as well as how noise in the filtered data behaves through the various integration steps described in Section 3. As can be seen by the ratio of the unfiltered data to the residual, the wide-band delay filter is responsible for suppressing foreground emission in the middle of the band by nearly three orders of magnitude (in mK), demonstrating the efficacy of this approach for isolating smooth-spectrum foregrounds from k modes that might be used to detect 21 cm reionization. This filter substantially outperforms time-differencing or frequency-differencing the visibilities for foreground suppression because of the sinusoidal nature of the fringes.

The gap in the lower set of curves of Figure 4 between the raw and differenced visibilities indicates the presence of an interfering signal that varies relatively smoothly in time and frequency. Such characteristics are generally indicative of foreground emission and suggest that near the limits of the sensitivity of these observations, foregrounds are not completely suppressed. Examining the set of visibilities in Figure 4 from which foregrounds have been removed via derivatives in time or frequency, we find a ratio of approximately 5 between the single-baseline and baseline-averaged curves in the 140–180 MHz frequency range. This ratio is roughly consistent with the value of 5.3 that would be expected for integrating 28 baselines with independent noise. Below 130 MHz, this ratio decreases substantially, indicating the presence of a correlated signal between different baselines. This correlated signal is shown to be suppressed by time-differencing and frequency-differencing filters, indicating that it is likely to be the result of unsuppressed foregrounds near the minimum delay threshold imposed by the wide-band delay filter. This interpretation is consistent with the power spectrum reported in Figure 6.

Scaling the noise-dominated single-baseline curve by the square root of the number of measurements contributing to each temperature residual, we determine the system temperature at 160 MHz to be 560 K. This value is in excess of the expected level of sky noise arising from galactic synchrotron emission, but it represents a substantial improvement over the system temperature in Parsons et al. (2010).

Next, we examine the set of visibilities derived from data after the application of a fringe-rate filter, which, although weighted for a sharp cutoff in fringe-rate domain, effectively corresponds to a ten-fold increase in integration time over the 43 s LST bins used in the upper plot. The ratio between the curves with and without fringe-rate filtering is approximately 3 in the 140–170 MHz frequency range. This ratio is roughly consistent with the value of 3.2 that would be expected for the ten-fold increase in integration time that the application of the fringe-rate filter represents. The ratio between the filtered and unfiltered curves is substantially lower on either end of this band, indicating a departure from noise-like behavior. As in other cases, the suppression of the systematic with frequency-differencing suggests that this reflects residual foregrounds that are relatively smooth. Nonetheless, the increasing systematics near the edge of the band form the basis of one of the arguments for observing with a wide instantaneous bandwidth. It also provides the motivation for limiting our power spectrum analysis in this initial investigation to the portion of the spectrum that appears consistent with noise.

In Figure 6, we show the spherically averaged 3D power spectrum derived from our measurements. In this figure, the central five data points in the range −0.06 ⩽ k_∥ ⩽ 0.06 h Mpc⁻¹ fall within the horizon limit of the delay spectrum and come from data without the application of the wide-band delay filter. We also illustrate the data before (cyan) and after (black) the application of the off-diagonal covariance removal described in Section 3.7. Averaging over ±k_∥, we place a 2σ upper limit on $\Delta _{21}^2(k)$ of 1660 mK² for k ≈ 0.27 h Mpc⁻¹. This limit is more than an order of magnitude more stringent than the previous best upper limit, which was recently revised to 61,500 mK² at k = 0.5 h Mpc⁻¹ at z = 8.6 (Paciga et al. 2013).

As this plot illustrates, the residual foreground signal that we detect in the frequency domain (see Figure 4) is concentrated at low |k_∥| modes in delay space, with the largest excess in a region just beyond the horizon limit described in P12b. In the range −0.22 ⩽ k_∥ ⩽ 0.22 h Mpc⁻¹, bins show evidence of a systematic bias that makes them inconsistent with zero with high significance. The source of this biasing signal is likely to be foreground emission, as predicted in P12b, but instrumental systematics are another possibility. We observe a second feature at moderate (3σ) significance at k_∥ ≈ ±0.35 h Mpc⁻¹ that appears relatively symmetrical with respect to positive and negative k_∥. We presume that this feature is caused by instrumental systematics or foreground emission, such as leakage from a polarized signal with high rotation measure (Moore et al. 2013). The underlying cause is still under investigation and will require substantial follow-up observation and careful analysis to validate.

As a diagnostic, we apply our analysis to random noise injected into the analysis before the narrow-band delay transform (see Figure 2) at the level expected for these observations with a system temperature of 560 K, which was measured at 160 MHz using frequency-domain visibilities after the application of a wide-band delay filter (Figure 4). These noise simulations agree with Equation (27) of Parsons et al. (2012a), which is the basis of the dashed cyan curve in Figure 6. We find good agreement between this expected noise level, the estimated power spectrum, and the errors derived from bootstrapping.

The limits shown in Figure 6 are not yet at the level of the brighter predictions favored by numerical and seminumerical simulations (McQuinn et al. 2007; Zahn et al. 2007; Iliev et al. 2008; Santos et al. 2008; Mesinger et al. 2011) that include, among other things, the predicted effects of X-ray heating on the intergalactic medium (IGM). Nonetheless, we can, for the first time,²⁰ constrain the 21 cm brightness temperature, 〈T_b〉, of the neutral gas at mean density for several reionization models and ionization fractions.

We recall that the brightness temperature contrast of the 21 cm transition relative to the cosmic microwave background (CMB; neglecting velocity gradient terms and assuming low optical depth, τ₂₁) can be written in the form

$$\begin{eqnarray} \Delta T &\approx \frac{T_S - T_{{\rm CMB}}}{1 +z} \tau _{21} \nonumber \\ &\approx T_0(z) \left(\frac{T_S - T_{{\rm CMB}}}{T_S} \right) x_{{\rm H}\,\scriptsize{I}} (1 + \delta), \end{eqnarray} \tag{ 6 }$$

where T_S is the spin temperature, $x_{{\rm H}\,\scriptsize{I}}$ is the neutral fraction, and δ is the density contrast, all of which are spatially dependent. The pre-factor is then

$$\begin{eqnarray} T_0(z) & = 23\, {\rm mK} \left(\frac{\Omega _b h^2}{0.02} \right) \left[ \left(\frac{0.15}{\Omega _m h^2} \right) \left(\frac{1+z}{10} \right) \right]^{1/2}\nonumber \\ & = 25\, {\rm mK} \left(\frac{1+z}{10} \right)^{1/2} \end{eqnarray} \tag{ 7 }$$

(see, e.g., Zaldarriaga et al. 2004; Furlanetto et al. 2006), where we have used the Planck parameters (Planck Collaboration et al. 2013).

Most reionization models assume an early period of UV-photon production from stars that, while insufficient to ionize the IGM, serves to drive the spin temperature toward the gas temperature via the Wouthuysen–Field effect. This is followed by heating of the neutral gas, primarily through X-rays associated with star formation and black hole accretion. In these scenarios, heating drives T_S ≫ T_CMB during reionization, so ξ ≡ (T_S − T_CMB)/T_S ≈ 1. The brightness contrast can be significantly larger, however, if T_S < T_CMB. In this case, the 21 cm line is seen in absorption against the CMB.

The maximum value of this contrast corresponds physically to a situation in which the spin temperature remains tightly coupled to the gas temperature via the Wouthuysen–Field effect²¹ (Wouthuysen 1952; Field 1958; Hirata 2006), but in which heating of the neutral gas by any mechanism (X-rays, Lyα, or shocks) is negligible. While a full characterization of this "no-heating" scenario would require a detailed simulation, we note that a reasonable rescaling of power spectra for which ξ was assumed to be one (as in Lidz et al. 2008) can be obtained in the following way. Because reionization tends to rapidly produce either fully neutral or fully ionized regions, $x_{{\rm H}\,\scriptsize{I}}$ can be locally represented as either 0 (for an ionized region) or 1 (for a neutral one). Because we have assumed no heating of the neutral medium, the spatial dependence of ξ is negligible. This turns the product T₀ξ into an overall multiplicative factor on the temperature fluctuations in Equation (6) and thus rescales the dimensionless ionization power spectrum, $\Delta ^2_i(k)$ (calculated from analytic or numerical models), as

$$\begin{equation} \Delta ^2_{21}(k) = T_0^2(z) \xi ^2(z) \Delta ^2_i(k) \equiv \langle T_b\rangle ^2 \Delta ^2_i(k). \end{equation} \tag{ 8 }$$

We can estimate the maximum 〈T_b〉 produced in the no-heating scenario by noting that the CMB and 21 cm spin temperature must stay in equilibrium due to residual Compton scattering until about z ∼ 150 (Furlanetto et al. 2006). Maximum contrast occurs if T_S tracks the gas temperature thereafter, so we can parameterize the maximum value of ξ as

$$\begin{equation} \xi _{\rm max} = 1 - \frac{1+z_{\rm dec}}{1 + z} \approx - \frac{150}{1+z}, \end{equation} \tag{ 9 }$$

where z_dec is the redshift at which the kinetic temperature of the IGM gas and the CMB drop out of equilibrium. Combining Equations (7) and (9) leads to

$$\begin{equation} \max [T_0(z) \xi (z)] \approx 370\, {\rm mK} \left(\frac{10}{1+z} \right)^{1/2}. \end{equation} \tag{ 10 }$$

To obtain a bound on 〈T_b〉, we treat ξ as a free parameter ranging from ξ_max ⩽ ξ ⩽ 1 that parameterizes the relative degree of heating of the IGM. We first note that brightness temperature fluctuations in excess of Equation (10) are ruled out a priori at a given redshift, as indicated by the maximum color scale corresponding to 〈T_b〉 > 400 mK in Figure 8. We then examine three different models for $\Delta ^2_i(k)$ at various ionization states and determine for each the maximum 〈T_b〉, beyond which the predicted value for $\Delta ^2_{21}(k)$ is inconsistent at the 2σ level with the measurements shown in Figure 6.

Zoom In Zoom Out Reset image size

We examine two models presented in Lidz et al. (2008) based on simulations by McQuinn et al. (2007). Although the ionization power spectra are associated with suggested redshifts in that paper, they apply to any redshift for the indicated ionization fraction, x_i. The first model that we examine is the fiducial "S1" model. For the x_i = 0.5 and x_i = 0.75 ionization power spectra, we obtain limits of 〈T_b〉 < 275 mK and 〈T_b〉 < 291 mK, respectively. These limits are illustrated in the upper two panels in Figure 9.

Zoom In Zoom Out Reset image size

In the second model that we examine, "S3," reionization is dominated by emission from massive halos, and the 21 cm power spectrum peaks at lower k modes in the final stages of reionization because of the increased spacing between the dominant sources of ionizing photons. For this model, we obtain limits of 〈T_b〉 < 217 mK and 〈T_b〉 < 264 for x_i = 0.5 and 0.75, respectively. These limits are illustrated in the lower two panels in Figure 9.

The third model that we consider is a toy "patchy" reionization model described by Equation (21) in P12a, which approximates $\Delta ^2_{21}(k)$ as being flat between the bounds of a minimum (k_min) and maximum (k_max) cutoff in k, scaled appropriately so that total power is conserved:

$$\begin{equation} \Delta ^2_i(k) = \left(x_{{\rm H}\,\scriptsize{I}}-x_{{\rm H}\,\scriptsize{I}}^2 \right) / \ln (k_{\rm max}/k_{\rm min}), \end{equation} \tag{ 11 }$$

where $x_{{\rm H}\,\scriptsize{I}}=1-x_i$ is the neutral hydrogen fraction. Using this toy model, we tune k_min and k_max and determine a maximum mean temperature above which Equation (8) becomes inconsistent with our measurements, as shown in Figure 8. While our measurements are not yet able to exclude 〈T_b〉 ≈ 30 mK predicted for an X-ray-heated IGM under the reasonable assumption that k_max/k_min > 10, we are able rule out 〈T_b〉 ⩾ 370 mK for k_min < 0.1 h Mpc⁻¹ for nearly all physically motivated values of k_max in both the x_i = 0.3 and x_i = 0.5 models. For x_i = 0.1, this result becomes dependent on the value of k_max.

The constraints shown in Figures 8 and 9 indicate that, for several reionization models and ionization fractions, our measurements are inconsistent with the 21 cm brightness temperature from neutral regions that would result from the gas in the IGM cooling strictly according to adiabatic expansion after recombination without additional heating. With the caveats that (1) the models we use may not correspond to the actual ionization power spectrum of reionization and (2) we do not know whether the ionization fraction at z ≈ 8 is near to the ionization fractions we investigate (although models such as Zahn et al. 2007 and Lidz et al. 2008 do predict x_i = 0.5 near z ≈ 8), the implication of our measurements is that the IGM may have been heated with respect to adiabatic cooling. X-ray heating from inverse Compton scattering from supernovae, high-mass X-ray binaries, or miniquasars is the most likely culprit (Madau et al. 2004; Ricotti & Ostriker 2004; Mirabel et al. 2011; Tanaka et al. 2012), although shock heating, which is currently disfavored (McQuinn & O'Leary 2012), has not been completely ruled out as a possibility.

We have used the PAPER results to explore some basic models for reionization. Under the simplifying, although perhaps not unreasonable, assumption that 〈T_b〉 is roughly factorable from the growth of structure (e.g., ionization bubbles and/or cosmic density), we constrain 〈T_b〉 in Figures 8 and 9 using limits to brightness temperature fluctuations observed by PAPER. For example, assuming a typical minimum k scale of ionization fluctuations of 0.1 h Mpc⁻¹ and a maximum of 10 h Mpc⁻¹, then as shown in Figure 8, the PAPER results imply upper limits on |〈T_b〉| of 191 and 175 mK for neutral fractions of 30% and 50%, respectively, at z = 7.7. For currently favored reionization models dominated by emission from massive halos (McQuinn & O'Leary 2012), PAPER results imply upper limits between 217 and 264 mK.

In summary, for massive-halo-dominated reionization models and for a range of simple patchy reionization models, PAPER data imply that the value of 〈T_b〉 is less than 400 mK. The value of 400 mK is important because it reflects the expected 〈T_b〉 in the case where the spin temperature of the 21 cm line is coupled to the gas kinetic temperature, but the neutral gas is never heated on large scales due to, e.g., X-rays from high-mass X-ray binaries or miniquasars. Physically, this could correspond to a situation where the Lyα photons from the first galaxies are sufficient to couple the spin and kinetic temperatures throughout the neutral IGM via the Wouthuysen–Field effect, but they do not substantially heat the gas (Chen & Miralda-Escudé 2004), nor do supernovae, X-ray binaries, or miniquasars produce enough X-rays to heat the large-scale neutral IGM. We point out that while the X-ray heating constraints we present require a coupling between the spin and kinetic temperatures, we know that this coupling must be occurring at z = 7.7. Observed star formation rates at z = 8 (Bouwens et al. 2010; Schenker et al. 2013) exceed by almost an order of magnitude the threshold of 10⁻³ M_☉ yr Mpc³ required for UV photons to couple these temperatures (McQuinn & O'Leary 2012).

While this no-heating model is considered physically unlikely (Furlanetto et al. 2006), the argument can be turned around: the fact that we constrain fluctuations at the 400 mK level provides, for the first time, suggestive evidence for large-scale heating of the neutral IGM during reionization at z > 7. Such constraints on X-ray heating may soon have repercussions on the global signal (Pritchard & Loeb 2012) and on the morphology of ionization (Mesinger et al. 2013). Note that our measurements pertain to the thermal evolution of the neutral gas during reionization and not to the subsequent thermal evolution of the ionized gas after reionization (Hui & Haiman 2003).

We begin by reviewing the quadratic estimator formalism. In this formalism, the value of the power spectrum p_α in the αth k bin (i.e., P(k_α), also known as the "bandpower") can be optimally estimated by the quantity²⁴

$$\begin{equation} \widehat{p}_\alpha = \frac{1}{2} \sum _\beta M_{\alpha \beta } {\bf x}^{\dag} {\bf C}^{-1} {\bf Q}^\beta {\bf C}^{-1} {\bf x}, \end{equation} \tag{ C1 }$$

where the hat signifies that $\widehat{p}_\alpha$ is an estimate of the true bandpower p_α, x is a vector containing binned data with ${\bf x}_\alpha \equiv \tilde{V}(\tau _\alpha)$ from Equation (2), C ≡ 〈xx†〉 is the covariance matrix of the data, M_αβ (discussed below) normalizes the estimator, and Q^β is the response of the covariance to the βth bandpower. The family of Q matrices (one for each value of β) maps the vector space of bandpowers to the vector space of data and is constructed to satisfy the relation

$$\begin{equation} {\bf C}^{\rm sky} = \sum _\beta p_\beta {\bf Q}^\beta, \end{equation} \tag{ C2 }$$

where C^sky is the contribution of sky emission to the covariance. This contribution includes both the foregrounds and cosmological signal. The foregrounds are included with the signal because we have conservatively avoided the direct modeling (and subtraction) of foregrounds from the data (Dillon et al. 2014). The total covariance is given by

$$\begin{equation} {\bf C}= {\bf N}+ {\bf C}^{{\rm junk}} + {\bf C}^{{\rm sky}}, \end{equation} \tag{ C3 }$$

where N is the instrumental noise covariance, and C^junk is the contribution of all contaminants except for foregrounds and noise.

Thinking of M_αβ as being an element of a matrix M that lives in a vector space of bandpowers, one is free to choose between several different choices of M, all of which are lossless (because M is generally taken to be an invertible matrix and all steps prior to its application are also lossless) but have different error properties. Consider for example the window function matrix W, which relates the expectation value of our bandpower estimators to the true bandpowers:

$$\begin{equation} \langle {\bf \widehat{p}} \rangle = {\bf W} {\bf p}, \end{equation} \tag{ C4 }$$

where we have grouped the bandpowers into vectors. The window function matrix can be shown to take the form

$$\begin{equation} {\bf W} = {\bf M}{\bf F}, \end{equation} \tag{ C5 }$$

where F is the Fisher information matrix, given by

$$\begin{equation} F_{\alpha \beta } = \frac{1}{2} {\rm tr} \left[ {\bf C}^{-1} {\bf Q}^\alpha {\bf C}^{-1} {\bf Q}^\beta \right]. \end{equation} \tag{ C6 }$$

Different choices for M give window functions (rows of W) with different properties, and one is subject only to the restriction that M must be normalized in a way that ensures that each window function integrates to unity to conserve power. We pick M = F⁻¹ without loss of generality, a fact that we will explain below (see Dillon et al. 2014 for other contexts where one might entertain other choices for M, such as taking M to be diagonal).

It is important to note that Equation (C1) estimates power spectra from single-baseline data. The length of x is therefore equal to the number of delay bins, and the covariance denoted here relates different delay bins from the same baseline only. The covariance matrix shown in Figure 5, in contrast, includes correlations between delay bins of different baselines as well. However, its analysis in Section 3.7 serves only to provide a better estimate for C by taking advantage of redundant baselines. In other words, our final result in Appendix C.2 will yield a procedure for removing off-diagonal systematics that operates on a per-baseline basis once C has been estimated, with information from different baselines mixing only in the cross-multiplication to form the power spectrum.

We now turn to applying the quadratic estimator formalism to delay-transformed visibilities. As derived in P12b, τ modes are a good approximation to k_∥, and their mutual orthogonality makes the simplified analysis in Section 3.6 a good approximation to optimal analysis. As a result, we only apply the full quadratic formalism to leading order in order to safeguard against the mismodeling of covariances, which are estimated empirically from a finite data set and are therefore imperfect representations of the true covariances. In this case, the leading-order optimal estimator reduces to the simple prescription given in Equation (4), provided that we used modified visibilities (as shown below) in place of $\tilde{V}(\tau,t)$.

Working in delay space, the Q matrices become particularly simple because delay modes are to an excellent approximation the same as Fourier modes along the line of sight (P12b). In this approximation, C_sky is diagonal, with elements equal to the bandpowers, i.e.,

$$\begin{equation} {\bf C}^{{\rm sky}}_{AB} = \delta _{AB} p_A, \end{equation} \tag{ C7 }$$

where A, B are index delay modes, and δ_AB is the Kronecker delta function. It follows that in this basis, the Q matrices simply pick out the correct delay/Fourier component:

$$\begin{equation} {\bf Q}^\beta _{AB} = \frac{\Omega B}{X^2Y}\left(\frac{2k_{\rm B}}{\lambda ^2}\right)^2 \delta _{AB} \delta _{A\beta }. \end{equation} \tag{ C8 }$$

Like the sky signal, the instrumental noise covariance can be modeled as a diagonal matrix, so that

$$\begin{equation} {\bf N}_{AB} = \delta _{AB} \sigma _A^2, \end{equation} \tag{ C9 }$$

where σ_A is the rms noise fluctuation in the Ath delay bin.

With ideal data, where C_junk = 0, the preceding discussion implies that C should be diagonal. However, as we can see from Figure 5, the actual data covariance contains nonideal off-diagonal entries. We model the real covariance matrix as

$$\begin{equation} {\bf C}= \left(\begin{array}{@{}ccccccccc@{}}r_1^2 & \varepsilon _{12} r_1 r_2 & \varepsilon _{13} r_1 r_3 & \\ \varepsilon _{21} r_2 r_1 & r_2^2 & \varepsilon _{23} r_2 r_3 & \cdots \\ \varepsilon _{31} r_3 r_1 & \varepsilon _{32} r_3 r_2 & r_3^2 & \\ & \vdots & & \ddots \\ \end{array} \right), \end{equation} \tag{ C10 }$$

where $r_A \equiv \sqrt{p_A + \sigma _A^2}$ is the total rms power (sky plus noise) in the Ath delay bin, and ε_AB is the coupling constant between the Ath and Bth bins. Since C is Hermitian, we note that $\varepsilon _{AB} = \varepsilon _{BA}^*$. Note that so far, our parameterization of C is still completely general because we have yet to specify the ε parameters; the fact that the off-diagonal entries are scaled to the diagonals is simply a convenient convention driven by the observation that our empirically derived covariances seem to roughly follow this form. By comparing Equations (C3), (C7), and (C9), we see that²⁵

$$\begin{equation} {\bf C}_{{\rm junk}} = \left(\begin{array}{@{}ccccccccc@{}}0 & \varepsilon _{12} r_1 r_2 & \varepsilon _{13} r_1 r_3 & \\ \varepsilon _{21} r_2 r_1 & 0 & \varepsilon _{23} r_2 r_3 & \cdots \\ \varepsilon _{31} r_3 r_1 & \varepsilon _{32} r_3 r_2 & 0 & \\ & \vdots & & \ddots \\ \end{array} \right). \end{equation} \tag{ C11 }$$

We will now proceed to compute the various quantities required for our power spectrum estimator, Equation (C1). First, we compute C⁻¹, which will be our primary tool for mitigating the influence of the junk contribution. Instead of performing a brute-force inversion of C, we can make the observation that it is equal to²⁶

$$\begin{eqnarray} {\bf C}^{-1} &\equiv & ({\bf N}+ {\bf C}^{{\rm sky}}+ {\bf C}^{{\rm junk}})^{-1} = ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} \nonumber\\ && - ({\bf N}+ {\bf C}^{{\rm sky}}+ {\bf C}^{{\rm junk}})^{-1} {\bf C}^{{\rm junk}} ({\bf N}+ {\bf C}^{{\rm sky}})^{-1}. \quad\quad \end{eqnarray} \tag{ C12 }$$

Taking advantage of the fact that (N + C^sky + C^junk)⁻¹ appears on the right-hand side, we can iterate this expression into a series:

$$\begin{eqnarray} {\bf C}^{-1} &=& ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} - ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} {\bf C}^{{\rm junk}} ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} \nonumber \\ &&+ ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} {\bf C}^{{\rm junk}} ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} {\bf C}^{{\rm junk}} ({\bf N}+ {\bf C}^{{\rm sky}})^{-1} \nonumber\\ && + \dots \end{eqnarray} \tag{ C13 }$$

Writing C⁻¹ in this form makes our results more robust to "noise" that arises from our empirical determination of the covariance matrix. For example, our series expansion does not require the inversion of C_junk, which may be unstable under inversion thanks to our imperfect estimation of the covariance. One only needs to invert the combination N + C^sky, which is trivial because both N and C^sky are diagonal. With a systematic series expansion, we can also henceforth work only to leading order in the ε parameters. We make this choice based on the intuition that uncertainties in our covariance estimation mean that higher order terms cannot be reliably estimated. A leading-order approximation may therefore be preferable to using an overly "noisy" covariance to determine C⁻¹, as that would result in an overfitting of the noise in the data (Tegmark et al. 1998).

Taking only the first two terms of Equation (C13), one obtains

$$\begin{equation} {\bf C}^{-1} \approx \left(\begin{array}{@{}ccccccccc@{}}1/r_1^{2} & -\varepsilon _{12}/r_1 r_2 & -\varepsilon _{13}/r_1 r_3 & \\ -\varepsilon _{21}/r_2 r_1 & 1/r_2^{2} & -\varepsilon _{23}/r_2 r_3 & \cdots \\ -\varepsilon _{31} r_3 r_1 & -\varepsilon _{32}/r_3 r_2 & 1/r_3^{2} & \\ & \vdots & & \ddots \\ \end{array} \right). \end{equation} \tag{ C14 }$$

Inserting this into Equation (C6), it is straightforward to show that the effects of C_junk affect the Fisher matrix only to second order in the ε parameters and may therefore be neglected:

$$\begin{equation} F_{\alpha \beta } = \left(\frac{\Omega B}{X^2Y}\right)^2\left(\frac{2k_{\rm B}}{\lambda ^2}\right)^4 \frac{\delta _{\alpha \beta }}{2 r_i^4} + \mathcal {O}(\varepsilon ^2). \end{equation} \tag{ C15 }$$

The final ingredient that we require before we can estimate power spectra is the normalization matrix M. As we remarked above, we will pick M = F⁻¹. If we had instead picked some popular alternatives such as taking M to be diagonal, it would have made no difference because we have just shown that the Fisher matrix is diagonal to leading order.

With a diagonal Fisher matrix, we can further simplify our power spectrum estimator. We first define a family of matrices

$$\begin{equation} {\bf G}_\alpha \equiv \frac{1}{2} \sum _\beta M_{\alpha \beta } {\bf Q}^\beta \end{equation} \tag{ C16 }$$

so that our estimator can be written compactly as

$$\begin{equation} \widehat{p}_\alpha = | {\bf G}_\alpha ^{1/2} {\bf C}^{-1} {\bf x}|^2. \end{equation} \tag{ C17 }$$

Inserting M = F⁻¹ as well as Equations (C8) and (C15) into our definition of G_α yields

$$\begin{equation} {\bf G}^\alpha _{AB} = \frac{X^2Y}{\Omega B}\left(\frac{\lambda ^2}{2k_{\rm B}}\right)^2 r_\alpha ^4 \delta _{AB} \delta _{A\alpha }. \end{equation} \tag{ C18 }$$

Substituting this and Equation (C14) into Equation (C17), we obtain

$$\begin{equation} \widehat{p}_\alpha = \frac{X^2Y}{\Omega B}\left(\frac{\lambda ^2}{2k_{\rm B}}\right)^2 {|} x_\alpha - \sum _{\beta \ne \alpha } \varepsilon _{\alpha \beta } \frac{r_\alpha }{r_\beta } x_\beta {|}^2. \end{equation} \tag{ C19 }$$

This, then, is our leading-order optimal estimator for the bandpower p_α in the presence of a "junk" covariance that couples different delay bins. Intuitively, our estimator instructs us to subtract from each delay bin a predicted leakage from all other bins, which is determined by a coupling constant and the ratio of rms power levels between the bins.²⁷ The resulting quantity is then squared to yield a final power estimate.

Written in this way, it is also clear how one would generalize the result to cross-power spectra. The optimal estimator for cross-power spectra between data from baseline i and baseline j, for example, would be

$$\begin{eqnarray} \widehat{p}_\alpha &=& \frac{X^2Y}{\Omega B}\left(\frac{\lambda ^2}{2k_{\rm B}}\right)^2 \left(x_\alpha - \sum _{\beta \ne \alpha } \varepsilon _{\alpha \beta } \frac{r_\alpha }{r_\beta } x_\beta \right)_i^{\dag} \nonumber\\ && \times \left(x_\alpha - \sum _{\gamma \ne \alpha } \varepsilon _{\alpha \gamma } \frac{r_\alpha }{r_\gamma } x_\gamma \right)_j. \end{eqnarray} \tag{ C20 }$$

In summary, if we define

$$\begin{equation} \tilde{V}^\prime _i(\tau _\alpha)\equiv \left(x_\alpha - \sum _{\beta \ne \alpha } \varepsilon _{\alpha \beta } \frac{r_\alpha }{r_\beta } x_\beta \right)_i, \end{equation} \tag{ C21 }$$

then Equation (C20) reveals that Equation (4) is to leading order an optimal estimator, provided we use $\tilde{V}^\prime$ in lieu of $\tilde{V}$.