FIRST SEASON MWA EOR POWER SPECTRUM RESULTS AT REDSHIFT 7

The preprocessing step converts the data from the non-standard memory dump format from the correlator's GPUs to a standard UVFITS format (Greisen 2012), phases the data to the center of the pointing, performs flagging, averages in time and frequency, and writes to the disk. For this step we use the COTTER package (Offringa et al. 2015), which in turn calls the AOFLAGGER ³⁵ to perform RFI flagging (Offringa et al. 2010).

Besides RFI, we manually flag data to remove three other effects. As mentioned in Section 2, the edges of the coarse bands contain aliasing from the filter applied in the digital receivers. We flag two 40 kHz channels on either side of the coarse-band edges (total of four channels per coarse-band edge). In principle, the aliasing could be calibrated out and these channels could be recovered, but this is left for future work. Second, early data from the MWA also contained some instances where antennas were not pointed properly at the beginning of observations due to a beamformer error. This problem has since been resolved, but exists in our data. We therefore flag the first two seconds of each observation to avoid any potentially mis-pointed data. Finally, we flag the central 40 kHz channel of each coarse band. This channel corresponds to the coarse-band DC mode, which has been observed to contain anomalous power, likely due to very low-level rounding errors in the polyphase filter bank of the digital receivers.

In the final step of preprocessing, COTTER averages the data to 2 s, 80 kHz resolution and writes to the disk. At this stage, each snapshot is contained in a single UVFITS file, which serves as the base unit for calibration and imaging.

The next step in our pipeline is to calibrate the visibilities. As shown in Figure 1, a sky model is derived from our imaged data, which feeds back into the calibration. This is naturally an iterative process as the sky model is refined with improved calibration solutions. We will describe the sky model in detail in Section 4, but here we present the production mode of analysis after the model has been determined.

Calibration is accomplished within the FHD package in two steps. First we match the raw data to model visibilities formed using a realistic simulation of the array response to our foreground model. This step results in independent complex gains for each antenna, 80 kHz frequency channel, and polarization. In the second step we impose restrictions on these gains motivated by our understanding of the spectral response of the instrument. By reducing the number of free parameters in this step we increase our signal to noise on calibration solutions and avoid absorbing unmodeled confusion source flux density into our gain solutions which has been shown to contaminate the EoR window (Barry et al. 2016). Recently, Patil et al. (2016) demonstrated excess noise and a suppression of diffuse foregrounds in LOFAR images after calibration and subtraction of bright discrete foreground sources. They suggested a multi-frequency calibration solution as a potential solution to these systematics, similar to the scheme implemented here. Indeed we have not seen evidence of these systematics through simulations and direct propagation of our noise (B. J. Hazelton et al. 2016, in preparation; Barry et al. 2016).

In the first calibration step, we allow our complex gains to account for any direction-independent response on a per antenna, per 80 kHz frequency channel, per polarization basis. To define our gains in this way, we have made two simplifications. First, we push all direction-dependence of the antenna response into the model of the primary beam. FHD is capable of using separate models for each antenna, though this has not been implemented in this work. The second simplification is to ignore any cross terms between our different axes; for example, a mutual coupling between two antennas. In principle, these terms would introduce baseline dependent gains, rather than antenna dependent. We have not yet seen evidence of these terms in MWA data, and so we currently neglect them in our solutions.

Under these assumptions we can express the measured, uncalibrated visibility for an antenna pair (i, j) as

$$\begin{eqnarray}&&{V}_{{ij}}^{\prime }(\nu )\approx {g}_{i}(\nu ){g}_{j}^{* }(\nu ){V}_{{ij}}(\nu ),\end{eqnarray} \tag{ 1 }$$

where ${g}_{i}(\nu )$ and ${g}_{j}(\nu )$ are the complex gains for antennas i and j, respectively, and ${V}_{{ij}}(\nu )$ is the true visibility, which we aim to recover through calibration. Because we are treating the polarizations independently, we allow the antenna subscript to run over both east–west (E–W) and north–south (N–S) polarizations.

FHD utilizes the StEFCal algorithm described in Salvini & Wijnholds (2014) to find the minimum ${\chi }^{2}$ estimate with respect to the complex gains. The result of this operation is estimated gains for every antenna, frequency channel, and polarization. However, with certain known properties of the antenna response we can reduce the number of free parameters in our solution. This takes us to the second step of calibration.

Our initial gain model included an amplitude bandpass common to all antennas, $B(\nu )$, and an antenna-dependent low-order polynomial in frequency. Mathematically, we can express our restricted gain as

$$\begin{eqnarray}&&{\hat{g}}_{i}(\nu )=B(\nu ){P}_{i}(\nu ),\end{eqnarray} \tag{ 2 }$$

where the polynomial term can be further decomposed as

$$\begin{eqnarray}&&{P}_{i}(\nu )=({A}_{i,0}+\nu {A}_{i,1}+{\nu }^{2}{A}_{i,2}){e}^{i({\phi }_{i,0}+\nu {\phi }_{i,1})}.\end{eqnarray} \tag{ 3 }$$

The coefficients $B(\nu )$, ${A}_{i,n}$, and ${\phi }_{i,n}$ are real quantities. This model allows us to capture arbitrary spectral responses due to common antenna factors (such as the polyphase filter shapes or phased array response), as well as slowly varying antenna-dependent deviations. We also found it necessary to include terms dependent on the length of the cable connecting the beamformers to the receivers, which we will discuss in Section 4.1.

Once all calibration terms are found, we complete the calibration by dividing the raw data by the gain estimates to produce calibrated visibilities. The entire calibration process is done for every snapshot of EoR data, providing both a two-minute resolution time dependence of calibration solutions, as well as model visibilities that are used for foreground subtraction and diagnostic purposes.

After calibration, we form snapshot image cubes (frequency mapping to line-of-sight direction). The visibilities are gridded using the primary beam as the gridding kernel, placing the data in the "holographic frame," which has been shown to be the optimal weighting scheme for combining images in the sense that it preserves all information for parameter estimation (Bhatnagar et al. 2008; Morales & Matejek 2009). We use a simulation-based primary beam model for MWA antennas developed by Sutinjo et al. (2015), which incorporates an average embedded element pattern for the dipoles, as well as mutual coupling between dipoles within an antenna. At this stage, we also average in frequency, by a factor of two, by gridding pairs of frequency channels to the same (u, v) plane, resulting in a frequency resolution of 160 kHz. This is done to reduce the data volume and has little impact on our results, because our sensitivity to the EoR is extremely low at the corresponding line-of-sight k modes (Beardsley et al. 2013).

The data are gridded into separate cubes for even/odd interleaved time samples. The sum of the even and odd cubes contains power from the sky as well as noise, while the difference between the two contains only noise power. Subtracting the power spectrum of the difference from that of the sum yields an unbiased estimate of the sky power (Jacobs et al. 2016).

The gridded (u, v) data are then Fourier transformed to create sky-frame images. To avoid aliasing, we image out 90° from phase center (gridding resolution of a half wavelength), and crop the image. We found that cropping the image based on beam value resulted in hard edges when integrating images from different snapshots having beams pointed in slightly different directions. To avoid the hard edges, we predetermine a set of HEALPix pixels to interpolate to, and use the same set for all snapshots on a given field. The final cropped field of view for this work is a $21^\circ$ square centered at R.A. = 0^h, decl. = −27°.

In principle the foregrounds could be subtracted from the data immediately after calibration, before gridding and imaging. However, for diagnostic reasons, we have found it beneficial to carry the model through the entire pipeline, and only subtract just before squaring to power spectrum units, allowing us to form power spectra of the dirty, model, and residual data.

FHD provides all inputs needed for the ppsilon package, which include weight cubes (for proper accumulation of data), variance cubes (for error propagation), and even/odd interleaved data cubes (for an unbiased estimator and a direct measurement of the noise). We also retain both E–W and N–S polarizations. At this time, we can remap our coordinates to cosmological units according to the relationships given in Morales & Hewitt (2004), where angular and frequency units $({\theta }_{x},{\theta }_{y},f)$ map to cosmological distance units $({r}_{x},{r}_{y},{r}_{z})$. Often, r_z is referred to as the line-of-sight dimension, and denoted as ${r}_{| | }$. Similarly, $({r}_{x},{r}_{y})$ is the plane perpendicular to the line of sight and is denoted as ${{\boldsymbol{r}}}_{\perp }$.

In the final steps of our pipeline, we integrate the snapshot image cubes, and calculate a power spectrum estimate. Though the integration step is formally an imaging component, we conceptually group it with the power spectrum step. This is because all steps of the pipeline up to this point have been performed on a per-snapshot basis (no communication between snapshots beyond foreground modeling), and the integration step can be used to select subsets of data for diagnosing power-spectrum artifacts. The images produced by FHD are in the holographic frame, which is already properly weighted for combining images, so the integration is simply adding the cubes together and propagating the weights.

The ppsilon package performs a discrete Fourier transform (DFT) on each ${r}_{| | }$ slice of the integrated HEALPix cube, forming a cube in $({{\boldsymbol{k}}}_{\perp },{r}_{| | })$, where ${{\boldsymbol{k}}}_{\perp }$ represents the cosmological wave number in the plane perpendicular to our line of sight, and ${r}_{| | }$ is the distance to the observed redshift slice along the line of sight. The Fourier transform along the ${r}_{| | }$ dimension is treated separately due to incomplete (u, v) sampling and flagged frequency channel; this leads to structure in the frequency sampling along any given k_⊥ pixel. We thus adopt the Lomb–Scargle least-squares method to determine the orthogonal eigenfunctions given our sampling function and estimate the total power in each ${k}_{| | }$ mode (Scargle 1982). In addition, because the spectrally smooth foregrounds contain vastly more power than the expected EoR signal, we apply a Blackman–Harris window function prior to the ${r}_{| | }$ to ${k}_{| | }$ transform. This has the effect of trading lower effective bandwidth for higher dynamic range (e.g., Thyagarajan et al. 2013, 2016).

The pipeline as described up to this point is applied to both the calibrated data and the foreground model. At this point, we subtract the model from the data to form residual data as well. All three sets of data (dirty, model, and residual) are carried through the remaining steps to form corresponding power spectra.

After squaring and dividing by our observation window function (Bowman et al. 2006), we arrive in 3D power spectrum space. We use the even/odd interleaved cubes to form a signal (odd plus even) and noise (odd minus even) power spectrum, which we subtract from one another to form an unbiased estimate of our signal power. This is mathematically equivalent to cross multiplying the even and odd cubes. The 3D power spectrum cube can next be averaged in annuli of constant k_⊥ and ${k}_{| | }$ (i.e., orthogonal to the ${k}_{| | }$ axis) to form 2D power spectra, or spherical shells to form one-dimensional (1D) power spectra where we will ultimately constrain the EoR.

An example 2D residual power spectrum formed from the three-hour golden data set is shown in Figure 2. The bulk of the residual power is in the so-called foreground "wedge," indicated with the diagonal black lines where the solid line corresponds to sky emission at the horizon, and the dashed line corresponds to emission at the edge of the MWA field of view. Above the wedge, we see horizontal lines of contamination due to the periodic frequency sampling function. These lines are often referred to as the coarse-band harmonics. We see vertical streaks at high ${k}_{| | }$ and high k_⊥. These are due to sparse (u, v) sampling beyond the dense core of the MWA (starting $\sim 70\lambda$), which results in non-uniform spectral sampling after gridding. This, in turn, causes the foreground power to mix to high spectral modes (Bowman et al. 2009). The frequency-dependent sampling creates a PSF in the ${k}_{| | }$ for each (u, v) cell. In principle, covariance weighting along the frequency dimension may be able to mitigate this leakage (e.g., Liu & Tegmark 2011), and an ppsilon implementation is under development. We provide axes more closely related to the measurements to assist in connecting these instrumental effects. The right axis shows delay (Fourier dual to frequency), and the top axis shows the baseline length measured in wavelengths. Between the coarse-band lines and to the left of the vertical streaks we see regions that contain both positive and negative (indicated by blue on the color bar) pixels. These regions are where we have successfully isolated the foregrounds and retained a noise-like EoR window in the three-hour integration.

Zoom In Zoom Out Reset image size

Figure 3 shows the uncertainty level corresponding to the 2D power spectrum in Figure 2. These levels are found by propagating noise measured from calibrated visibilities through the entire pipeline. The dense core of the MWA results in relatively low noise at low ${k}_{\perp }$. As the (u, v) density of the array decreases at longer baselines, the noise increases. The noise is mostly constant across ${k}_{| | }$ because the thermal noise amplitude is relatively constant over frequency and the Fourier transform distributes the noise power across all modes (Morales & Hewitt 2004). The exception is the second ${k}_{| | }$ bin, which has slightly lower noise than the rest of the ${k}_{| | }$ modes. As explained by B.J. Hazelton et al. (2016, in preparation), this is a result of the Lomb–Scargle periodogram interacting with the Blackman–Harris window function we apply to the data. Most of the bins have the noise equally distributed between cosine and sine terms, with the exception of the ${k}_{| | }=0$ bin, where the cosine term carries all the information. Due to the correlations introduced by the Blackman–Harris window function, the first non-zero ${k}_{| | }$ mode is highly covariant with the bottom bin—therefore, the noise on this mode's cosine term is lower than other higher order cosines.

Zoom In Zoom Out Reset image size

Early in our analysis it became apparent that the bandpass of each MWA antenna requires more free parameters than those described in Equation (2). In particular, differing lengths of cable connecting the beamformers to the receivers lead to different bandpass shapes due to signal attenuation. We therefore allow the amplitude bandpass factor to depend on cable length, ${B}_{\alpha }(\nu )$, where α denotes the cable type.³⁶

The top panel of Figure 4 shows a 2D power spectrum using the gain model described so far. Three of the horizontal lines in the EoR window can be attributed to the coarse-band gaps, as they reside at the harmonics expected for 1.28 MHz periodic sampling. But the sampling cannot account for the fourth line, highlighted by the arrow. The corresponding delay time of the line, $\tau \approx 1.23$ μs, corresponds almost exactly to twice the signal travel time through the 150 m cables (with velocity 0.81c). We therefore introduced a reflection term into our gain model, which allows our restricted gains to account for the interference of the incident signal with a round-trip reflected signal in the cable. A similar approach was taken at lower frequencies in Ewall-Wice et al. (2016). The full gain model is expressed as

$$\begin{eqnarray}&&{\hat{g}}_{i}(\nu )={B}_{\alpha }(\nu )({P}_{i}(\nu )+{R}_{i}(\nu )),\end{eqnarray} \tag{ 4 }$$

where the reflection term can be further decomposed as

$$\begin{eqnarray}&&{R}_{i}(\nu )={R}_{i,0}{e}^{-2\pi i{\tau }_{i}\nu }.\end{eqnarray} \tag{ 5 }$$

The reflection coefficient, ${R}_{i,0}$, is allowed to be complex, while the reflection delay, ${\tau }_{i}$, is real. After fitting for all other parameters, we fit for the reflection mode. We found that if we did not first fit and remove the polynomial terms, our reflection fits would be dominated by the large power in the smooth structure, resulting in solutions that did not match physical reality. The reflection fitting is done for antennas with suspected reflections, the most offensive of which is seen in the antennas with cables of 150 m. Because the cables were not cut at exact lengths (variations on order tens of centimeters), we solve for both the reflection coefficient and delay by performing a direct Fourier transform to a highly over-resolved delay grid and selecting the mode where the reflection amplitude is largest. This fitting produces a single reflection coefficient for each antenna with nominal cable length of 150 m.

Zoom In Zoom Out Reset image size

The resulting 2D power spectrum after fitting for cable reflections is shown in the bottom panel of Figure 4, where the reflection line is suppressed below the noise level. This example demonstrates the power of the 2D power spectrum as an FoM. The power in the reflection line was about five orders of magnitude below the foregrounds, making it very difficult to detect in image-based metrics. However, the power spectrum is specifically designed to be sensitive to low-level spectral structure, and using the two-dimensional spectrum allows us to identify the "shape" of contamination, enabling a precise diagnosis.

A similar, low-level effect that had the potential to contaminate the EoR window is shown in Figure 5. The spectra shown are the model spectra—calculated by propagating the sky model visibilities through the entire pipeline. Despite our foreground model not containing spectral structure, the EoR window in the top panel seems to have a floor at a level $\sim {10}^{7}$ mK² h⁻³ Mpc³, comparable to predicted EoR signals. In addition, there is a faint line in the upper left of the plot that has the same slope as the wedge, but seems to originate far beyond the horizon.

Zoom In Zoom Out Reset image size

The source of this floor and super-horizon line was traced to the resolution at which we formed the gridding kernel while gridding visibilities. In the interest of computational speed and efficiency, the kernel is pre-calculated at a high resolution, then a nearest-neighbor lookup table is used to approximate the beam values at discrete pixels. This is a common practice in most analysis software packages. The kernel resolution used in the top panel of Figure 5 was 0.04 wavelengths—much smaller than the half-wavelength grid corresponding to horizon-to-horizon imaging. However, as a baseline migrates in (u, v) coordinates across frequencies, it undergoes discrete steps between kernel values used. This effectively results in small baseline position errors, which shift periodically in frequency, resulting in power being mixed into the window and a relatively strong harmonic at the frequency of the shifting (the faint line in the upper left of the window). While we illustrate this effect with a model power spectrum, the same problem exists in the data, but is below the noise levels at three hours of integration.

This effect is distinct from more fundamental wedge effects (e.g., Morales et al. 2012; Hazelton et al. 2013) in that it is a result of the computational limitation of the analysis. The primary contribution to the foreground mode mixing is due to the loss of information as each baseline samples a range of (u, v) modes. On the other hand, the effect discussed here is due to small positional offsets in the gridding kernel and only exists after gridding.

We resolved this issue by increasing the resolution at which we form the kernel. While the most accurate answer is to model at infinite resolution, it is not computationally feasible to do so. Instead, we chose a resolution at which the effect no longer materially impacts the power spectrum. With experimentation, we found that a kernel resolution of 0.007 wavelengths was sufficient. The resulting model power spectrum is shown in the bottom panel of Figure 5. The contamination within the EoR window now drops significantly lower. Of course, without knowing the exact level of the cosmological signal, we cannot know if this level is sufficiently low. This effect may need to be revisited if the EoR power spectrum is lower than predicted. The improvement in beam model resolution came at the cost of 20% increase in memory usage for our imaging pipeline.

Our pipeline uses two modes of FHD—the full deconvolution mode to identify point sources and build a catalog, and a production mode, where we calibrate using the sky model and subtract it from the data without further fitting. Because full deconvolution on every observation from the MWA is computationally not feasible, we currently restrict ourselves to the golden data set to build our model. The process of building the catalog is presented in Carroll et al. (2016). They applied machine learning classification methods to select reliable detections from the full deconvolution FHD mode, culminating in the KGS catalog.³⁷ The work here uses an early iteration of the KGS catalog that was available at the time of analysis.

It has been shown that subtracting a foreground model strictly within the main lobe of the primary beam will not be sufficient to suppress the power spectrum wedge and unlock the EoR window (Thyagarajan et al. 2015a, 2015b; Pober et al. 2016). We are in the process of repeating the model building described above using additional MWA observations pointed away from our field of interest. Until that work is complete, we supplement the extent of our point-source model using additional catalogs. We accomplish this through a hierarchical catalog pulling from our early KGS catalog, the MWA Commissioning Survey³⁸ (MWACS, Hurley-Walker et al. 2014), the Culgoora catalog (Slee 1977), and the Molonglo Reference Catalog (MRC, Large et al. 1981). Source flux densities are prioritized in this order based on our confidence in their predicted flux density at 182 MHz. We first cluster the source lists to avoid redundant sources, using a 3.5 arcmin neighborhood radius, we then select the flux density from the highest priority catalog. Spectral indices were obtained from the MWACS and Culgoora catalogs when available, otherwise a two-point spectral index was estimated for Culgoora–MRC matches or MRC–SUMSS matches. All other sources were given a spectral index of −0.8, the previously determined median spectral index of sources below 1.4 GHz (Oort et al. 1988; Hunstead 1991; De Breuck et al. 2000; Mauch et al. 2003). The spectral index was used to extrapolate the flux density from the catalog frequency to 182 MHz, but a uniform spectral index of −0.8 was used within our band when forming model visibilities for calibration and foreground subtraction.

The resulting hierarchical catalog is shown in Figure 6. The EoR0 field corresponds to the red KGS source patch. All sources above the horizon and within a primary beam value greater than 1% maximum (including sidelobes) are included in both our calibration and foreground subtraction models.

Zoom In Zoom Out Reset image size

To show the effect our hierarchical catalog has on the power spectrum FoM, we ran our entire analysis pipeline on the golden data set with two foreground models—first with only the MWACS catalog, and again with the full hierarchical catalog. We then compare the resulting 2D power spectra to inspect whether the new catalog results in more accurate calibration and foreground subtraction.

Because we use the input foreground model to calibrate our visibilities, the gain estimates for our two runs differed. In particular, we saw the overall flux scale of the calibrated visibilities was higher when using the full hierarchical catalog, owing to a more complete sky model requiring lower amplitude gains to describe the data. This difference is on the order of a part in 10³ in mK² units, but is largely amplified when observing the difference in the power spectra. In order to put both residual 2D power spectra onto the same scale for comparison, we first divide each $({k}_{\perp },{k}_{| | })$ cell by the corresponding dirty power spectrum pixel, then subtract to arrive at a ratio difference,

$$\begin{eqnarray}&&{\mathtt{ratio}}\ {\mathtt{difference}}=\displaystyle \frac{{{\mathtt{residual}}}_{1}}{{{\mathtt{dirty}}}_{1}}-\displaystyle \frac{{{\mathtt{residual}}}_{2}}{{{\mathtt{dirty}}}_{2}}.\end{eqnarray} \tag{ 6 }$$

Here we use subscript 1 to represent the MWACS catalog spectra, and 2 for the hierarchical catalog spectra. The ratio difference is shown in Figure 7. The entire wedge being positive (blue) demonstrates that using the new catalog resulted in successfully subtracting a higher fractional power. Not only does this indicate that we subtract more power, but it also confirms that our calibrated data is more closely matched to the model, meaning that our calibration solutions, in general, are more accurate.

Zoom In Zoom Out Reset image size

Line et al. (2016) performed a similar analysis, but with an emphasis on the positional accuracy of catalog entries. By using simulations of visibilities, they were able to isolate the effect of using perfect versus offset source positions. In their simulations, source position offsets of 14% of the synthesized beam-width had a significant impact on foreground subtraction, confirming the importance of a complete and accurate point-source catalog in order to retain a clean EoR window.

In addition to a point-source foreground model, we also introduce a diffuse emission model within the primary field of view of the MWA. Though the EoR0 field was chosen to be relatively devoid of Galactic emission, we find this diffuse structure still highly contaminating due to the very high sensitivity of the MWA at large scales. For computational purposes, we have divided the diffuse emission into two regimes—the faint clouds in our main lobe, and the bright plane and other structures in the sidelobes. While the latter has been studied by many people and a Global Sky Model (GSM) is readily available (de Oliveira-Costa et al. 2008), the computational obstacle of simulating the instrument response to a full-sky diffuse model is yet prohibitive, though under active pursuit (Thyagarajan et al. 2015b). We instead focus here on the diffuse structure within the primary field and leave the full-sky model to future iterations of analysis.

This first iteration of modeling the diffuse structure within our primary field was done by simply using output point-source subtracted residual images of FHD from the three-hour golden data set. We combine the images from three hours to leverage the rotation of the Earth to improve the PSF of the instrument. This integrating of images is done with uniform weighted data to minimize the effect of double instrument convolution (once when the data are taken, and again when we use the output as a model). We then form a pseudo Stokes I image by adding beam-weighted E–W and N–S polarizations. We also integrate in frequency to form a single continuum image for our model. By forcing our model to contain no spectral structure, we mitigate the risk of subtracting a cosmological signal.

Future iterations of this model will contain a spectral index and multiple polarization components. Lenc et al. (2016) demonstrated the presence of strong polarized large-scale structures in the EoR0 field at 154 MHz with varying Faraday depths. They observed rotations in the Stokes Q–U plane due to ionospheric conditions, which will need to be taken into account in our future diffuse models.

Ultimately, we need model visibilities to subtract from our data. However, rather than storing a model for all observations (different pointings and phase centers), we find it simplest to treat the model as a single HEALPix image that can be used to create model visibilities for each independent observation. FHD treats each pixel in the image as a point source at the pixel location with the total flux density equal to the surface brightness of the diffuse image times the area of the pixel, and it creates model visibilities for each snapshot in the same way that it imports a catalog of point sources. This is similar to the strategy employed by Thyagarajan et al. (2015b) to model diffuse structure.

The diffuse model used for this work is shown in Figure 8. While the actual model image is uniform weighted with resolution $\sim 6$ arcmin, we show it smoothed to degree scales to emphasize the large-scale structure and to approximate what the MWA instrument observes with a natural weighting.

Zoom In Zoom Out Reset image size

For this iteration of analysis, we only use the diffuse model for subtraction, and omit it for calibration. At the time of writing, short baselines were not producing reliable results in the calibration loop. Instead, we chose to omit the diffuse model and mask baselines shorter than 50 wavelengths for calibration purposes. However, we add the diffuse model and unmask short baselines when performing foreground subtraction. We compare the total residual power in the images by squaring and summing the residual image cube from before and after the diffuse subtraction. This includes all scales measured by the instrument, but, because the image cubes are naturally weighted, the ${\boldsymbol{k}}$ space weighting is the inverse of Figure 3, meaning that ${k}_{\perp }\lesssim 0.05\,h$ Mpc⁻¹ dominates the average. We saw a 70% reduction in residual power when including our diffuse model—a strong indication that our diffuse model improves our subtraction. Future iterations will incorporate the model in calibration and gradually improve the model itself.

We demonstrate the impact of our diffuse model again by running the golden data set through our entire analysis pipeline with and without using the diffuse model. Because the foreground models used for calibration in the two runs are identical (diffuse is omitted for calibration), the calibration solutions were identical and therefore the dirty power spectra were identical as well. We compare the model and residual power spectra by direct subtraction in Figure 9. The top panel shows our point-source foreground model power spectrum minus the model power spectrum when including diffuse. Because our diffuse model added power to the full foreground model, the entire plot is negative (red). Even the EoR window is red (albeit at a much lower level than the wedge) because of power leakage from non-uniform spectral sampling and from the dynamic range limit of the Blackman–Harris window function. The bottom panel is the difference in residual power. The wedge is completely positive (blue), indicating that the diffuse model successfully subtracted from the dirty visibilities. The differences in the window are positive and negative, indicating that the differences in the foreground model in this region are below the three-hour noise level.

Zoom In Zoom Out Reset image size

All 2780 snapshots in the data set are preprocessed, calibrated, and imaged using the pipeline described in the previous sections of this article. Similar to tests on the golden data set, we rely heavily on the 2D power spectrum as a diagnostic tool, this time to identify and excise poor-quality data.

We use the jackknife method to filter out bad data and detect patterns. This involves dividing the data into many subgroups and forming power spectra. The goal is to find observational parameters that affect the quality of the data, such as the day or time of night. A powerful grouping is to divide the data into the observation day and "pointing" (which direction the antennas were pointed while tracking the EoR0 field). We show an example of one day's worth of per-pointing power spectra in Figure 10. The pointings are labeled sequentially with −5 corresponding to five pointing steps prior to zenith transit, 0 corresponding to zenith, and +4 corresponding to four steps after zenith. This numbering scheme is shown graphically in Figure 11. Early in the night (−5 through −3) the bulk of the Galactic plane was still above the horizon, and, despite being very far from our primary field of view, highly contaminated our observations through the sidelobes of the instrument. We then have relatively well-behaved pointings, with the EoR window dominated by noise, until the final pointing of the night, when we can see evidence of the Galaxy rising again indicated by strong lines of power at the edge of the foreground wedge. We saw the same contamination on all days of observation and decided to cut all −5, −4, −3, and +4 pointings from our data set. In principle, it may be possible to model the Galactic plane well enough to account for its presence, but we leave this to future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Motivated by the manual classification done with power-spectrum jackknives, we developed another metric to quickly predict the power spectrum quality for each snapshot in our data set. We do this by forming a delay spectrum (Parsons et al. 2012) from the raw, uncalibrated visibilities. We then calculate an estimated total EoR window power by integrating the power above the horizon line and below the first coarse-band line. While we would ideally use calibrated visibilities for this metric, we have found that the uncalibrated power is strongly correlated to the calibrated power, and the cuts we make below are independent of calibration. This is encouraging because future analysis could use this metric before processing the snapshots, saving valuable computing resources.³⁹

We plot the window power for each snapshot in our data set in Figure 12(a). We first note that our conclusion from the jackknife tests is confirmed here—early pointings contain strong contamination and have high window power. We also see many outlier snapshots with excess power. Inspecting power spectra from these individual snapshots confirmed poor data quality, even after calibration. The cause of these contaminants is yet unknown, but could easily be attributed to low-level RFI that was missed by the AOFLAGGER, or intermittent hardware failures. To remove the poor data, we made the cut shown with the black box, keeping only snapshots inside.

Zoom In Zoom Out Reset image size

Next we compare the power in our two instrumental polarizations. When conducting our jackknife tests, we saw contamination could exist in one polarization and not the other—potentially due to RFI, or an instrumental failure. We plot the ratio of N–S power to E–W power in Figure 12(b), after applying the previous cut. We see that the ratio is generally flat, with the exception of some outliers and almost all of pointing +3. We suspect that the excess N–S power in the later pointing is due to the Galaxy leaking back into the sidelobe of the N–S polarization, but not in the other. Again, we make the cut shown with a black box.

Our final data cut was made when inspecting the snapshot residual continuum images output from FHD. We estimate a proxy for the residual flux density from foreground sources by calculating the rms of the fractional residual flux densities for all subtracted sources greater than 0.5 Jy within half-beam power. The fractional residual flux density is measured as the ratio of the residual image pixel value (Jy) at the source position to the integrated flux density of the source. We found that most snapshots had residual flux density rms $\lt 10 \%$ with some outliers, which we cut from the data. This cut largely overlapped with previous cuts, but removed an additional 95 observations. The cause of the high-residual flux density rms is yet unknown and is not localized by observation day, LST, pointing, or any other jackknife we have performed.

After the cuts described above, we are left with 1029 snapshots, or just over 32 hours of data. This is an aggressive cut, with the aim of removing all poor-quality data and retaining a data set as clean as we can determine from our data quality checks. As mentioned earlier, more sophisticated analysis could allow us to recover data cut from this study in the future (e.g., modeling and removing the Galaxy in observations pointing far off zenith).

In this section, we discuss the results after processing the remaining 32 hours of data through the ppsilon pipeline, and place an upper limit on the EoR power spectrum.

The bandwidth of the MWA and our data set is 30.72 MHz, but in order to avoid the effects of cosmic evolution over the span of our measured redshift range, we limit observations to about 8 MHz, or ${\rm{\Delta }}z\sim 0.3$ at our frequency. We do this by dividing the band into three overlapping sub-bands of 15.36 MHz each, which, after applying the Blackman–Harris window function, will result in effective bandwidths of 7.68 MHz. We label these sub-bands as "low" ($z\approx 7.1$), "mid" ($z\approx 6.8$), and "high" ($z\approx 6.5$).

The 32-hour integrated residual power spectra for our three sub-bands and both instrumental polarizations are shown in Figure 13. A number of features can be seen in these spectra. In all cases the foreground wedge is prominent, with extra power at large scales (low k_⊥). This is an indication that our diffuse model needs to improve for deeper subtraction.

Zoom In Zoom Out Reset image size

The lowest region of the EoR window, above the horizon line and around ${k}_{| | }\approx 0.25\,h$ Mpc⁻¹, contains purely positive bins, indicating non-noise-like power. The cause for this seeming leakage is yet unknown, and it limits our integration in this analysis. One suspected origin of this contamination was individual snapshots with high contamination that slipped through our cuts. However, dividing the data into random subsets and comparing results indicated that there is not a small number of offending observations, instead, the leakage appears to exist at a low level in all of the data.

Another potential cause of the foreground leakage into the EoR window is insufficient calibration quality—particularly the spectral shape of the instrumental bandpass. Investigations are underway to combine snapshot observations for higher signal-to-noise calibration solutions, as well as more sophisticated parameterization of the gain model in Equation (4) to account for more physical effects, such as dependence on ambient temperature variation.

Moving up in the window, between the first two coarse-band harmonic lines we can see an additional faint line. This is the re-emergence of the 150 m cable reflections. While the method described in Section 3.2 sufficiently calibrated out this reflection line for the three-hour golden set integration, the lower noise level in this deeper integration shows it is not completely removed. Future analysis will require higher signal to noise on the reflection fitting, which may require combining snapshots for calibration solutions.

Between the coarse-band harmonic lines and the reflection line, we do see regions where our spectra appear noise-like (positive and negative values). This is encouraging, despite the leakage at lower ${k}_{| | }$, and motivates us to be selective when binning to 1D power spectra.

Figure 14 shows slices of power that help isolate the contamination. The two-dimensional power spectrum is the N–S, mid-band spectrum from Figure 13. We have drawn a horizontal solid blue line to indicate the slice used to plot power as a function of k_⊥ for a fixed ${k}_{| | }$ (top right), and a vertical dashed line to show the slice used to plot power as a function of ${k}_{| | }$ for a fixed k_⊥ (bottom right). In the k_⊥ power plot, we can see leakage at high k_⊥ from the residual of the vertical streaks due to poor sampling at large (u, v). While the increased noise would mean that we will down weight the bins at large k_⊥, the contamination is even larger, meaning it would bias our estimate even though those bins are down weighted. The low end of k_⊥ contains exceptionally bright foregrounds at large scales, which, due to the width of the coarse-band harmonics, contaminate most of the ${k}_{| | }$ modes. With both of these regimes in consideration, we exclude bins with $u\lt 10\lambda$ and $u\gt 70\lambda$, shown with gray boxes. Between these two regimes are many bins consistent with the noise level, which we include in our 1D averages.

Zoom In Zoom Out Reset image size

Similarly, we inspect the ${k}_{| | }$ power in the bottom right plot of Figure 14. Here, we see the huge contamination from foregrounds at low ${k}_{| | }$. While the leakage drops significantly, our measurements do not reach the noise level until after the first coarse-band harmonic. However, we chose to include bins below the first coarse band with ${k}_{| | }\gt 0.15\,h$ Mpc⁻¹ because the cosmological signal is expected to contain substantially higher power at large scales, and ultimately these low ${k}_{| | }$ bins provide our most competitive limits, despite being systematic limited. We exclude the coarse-band harmonics by masking out the harmonic bin itself and two bins on either side (total of five bins per harmonic).

The final mask we apply is the wedge. We found that a small buffer beyond the horizon was necessary to completely mask out the wedge, consistent with Dillon et al. (2015a). We implement the buffer by increasing the slope of the horizon line by 14%. This line is shown as the dashed diagonal line in the left panel of Figure 14.

The various masking described above is summarized with the black contours shown on the two-dimensional power spectra in Figure 13. The contours show the cuts in two-dimensional space, but the masking is actually performed directly in the 3D power spectrum cube and averaged directly to one-dimension.

The resulting 1D power spectra are shown in Figure 15, where the measured power is shown with solid blue, the 1σ noise level is shown with thin red, and the 2σ upper limit for each bin is shown with magenta. Where our unbiased estimator is negative, the absolute value is shown, but with a dotted blue line. For consistency with the literature, we plot our 1D power spectra as ${{\rm{\Delta }}}^{2}(k)={k}^{3}{P}_{21}(k)/(2{\pi }^{2})$, which has units of mK².

Zoom In Zoom Out Reset image size

In all bands and polarizations, we are heavily signal dominated in our most sensitive region (low k). This is not surprising based on the 2D spectra we examined earlier. The gaps in the data are due to the excised coarse-band harmonics. We can see the cable reflection line between the first two gaps in all polarizations and sub-bands. Between the coarse-band harmonic lines, especially at larger k, are bins that approach the noise level and are consistent with zero (notably the N–S polarization for the low and mid bands). The bin just above the first coarse-band harmonic is also consistent with zero in the low band. These regions are encouraging because they have the potential to continue integrating down with more data without further analysis improvements. Further suppression of the coarse-band harmonic and cable reflection lines will also improve the limits presented here.

Our best upper limits for a cosmological signal are at low k, despite the strong leakage. Because any foreground leakage should not correlate with the EoR signal, we can assume that the power of the sum of leakage and cosmological signal is greater than the power of the EoR signal alone, allowing us to set an upper limit. We quote the best 2σ upper limits, ${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$, for each polarization and band in Table 1. These limits and their context in the field are further discussed in Section 6.

Table 1.  Upper Limits on the EoR Power Spectrum for Our Three Sub-bands and Two Polarizations

			FHD+ppsilon		RTS+CHIPS
Sub-band	z0	Pol	k	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$	k	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$
Low	7.1	E–W	0.231	3.67 × 104
Low	7.1	N–S	0.27	2.70 × 104	0.16	3.2 $\times {10}^{4}$
Mid	6.8	E–W	0.24	3.56 × 104
Mid	6.8	N–S	0.24	3.02 × 104	0.14	2.6 × 104
High	6.5	E–W	0.20	4.70 × 104
High	6.5	N–S	0.24	3.22 × 104	0.14	2.5 × 104

Note. Upper limits, ${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$, are at the 97.7% confidence level. Cosmological wave numbers, k, are in Units of h Mpc⁻¹ and upper limits are in units of mK². The last two columns, RTS+CHIPS, are produced from the reference pipeline, which is discussed in Section 5.3.

Download table as:  ASCIITypeset image

As demonstrated in Jacobs et al. (2016), comparison between independent analysis pipelines is essential to evaluate the validity of algorithms under active development. Such comparison provides insight into the strengths and weaknesses of the pipelines, as well as confirmation of results. The analysis discussed so far has been based on the FHD-to-ppsilon pipeline, and here we compare with the Real Time System (RTS, Mitchell et al. 2008; Ord et al. 2010) to Cosmological H I Power Spectrum (CHIPS, Trott et al. 2016) pipeline. Below, we present a brief overview of the reference pipeline and highlight the lessons learned from the comparison.

The RTS uses a fundamentally different calibration and foreground subtraction approach than FHD, as well as a different input point-source model. The primary differences lie in the use of ionospheric corrections, the bandpass calibration, and the lack of a diffuse model. RTS accomplishes the interferometric calibration in two steps: a direction-independent calibration to a single compound calibrator created by combining the 1000 apparently brightest sources in the sky, followed by an individual calibration and subtraction of those sources inside the Calibration Measurement Loop (CML, Mitchell et al. 2008). The CML estimates individual ionospheric corrections for all sources and performs a full direction-dependent calibration for the 5 very brightest ones. These tasks are distributed over independent parallel threads for each 1.28 MHz coarse channel.

The bandpass calibration is also fundamentally different to FHD. The bandpass for each antenna is fitted to a third-order polynomial within each 1.28 MHz coarse band and each direction and polarization individually. The delay of the 150 m cable reflection discussed earlier is comparable to the scale of the polynomial fit, so we expect it to naturally be fit out of the gain solutions. Indeed, the RTS has not yet seen evidence of the reflection in the processed data (Jacobs et al. 2016).

The source catalog used with RTS was a combined cross-match created with the Positional Update and Matching Algorithm (PUMA, Line et al. 2016). Full details of PUMA and this cross-matched catalog can be found in Line et al. (2016), but we include a brief description here. PUMA is designed to cross-match low radio-frequency ($\lesssim 1\,$ GHz) catalogs. It utilizes a Bayesian probabilistic positional cross-match approach based on Budavári & Szalay (2008), combined with criteria based on catalog resolution and fitting the spectral energy distribution to a power-law. These criteria allow for surveys with differing resolutions, which can introduce confusing sources into the matched process. Primarily, the source catalog was based on the MWACS (Hurley-Walker et al. 2014), which was cross-matched to the 74 MHz Very Large Array Low Frequency Sky Survey redux (VLSSr, Lane et al. 2012), the MRC, (Large et al. 1981), the $843\,$ MHz Sydney University Molonglo Sky Survey (SUMSS, Mauch et al. 2003), and the $1.4\,$ GHz NRAO VLA Sky Survey (NVSS, Condon et al. 1998). MWACS covers approximately $20\buildrel{\rm{h}}\over{.} 5\lt {\rm{R}}.{\rm{A}}.\,\lt \,8\buildrel{\rm{h}}\over{.} 5$, $-58^\circ \lt \delta \lt -14^\circ$, so to complete the sky coverage needed, MRC was used as a base catalog outside of the MWACS coverage and was matched to VLSSr, SUMSS and NVSS.

Since RTS has the added step of individual source peeling, there is a greater possibility for numerical calibration failure. Before passing to CHIPS, a simple quality assurance step is considered, which is based upon the variance of residual visibilities in each 1.28 MHz coarse band. Observations with divergent or unstable calibrations were readily identified by variances far outside the otherwise observed distribution of values. Of the 1029 observations calibrated and used by the FHD pipeline, the RTS+CHIPS pipeline successfully calibrated and included 1003. The remainder failed to calibrate successfully, potentially due to adverse ionospheric conditions.

The CHIPS power spectrum estimator applies a full maximum-likelihood estimator to the data. In its full form, the data are weighted by the thermal noise, and a model for the residual point-source signal, within a full frequency-dependent description of the instrument. For this comparison, the foreground weighting is not used, in order for a direct comparison with the weighting scheme used by ppsilon. The data were divided into the same three redshift bins, and the 2D and 1D power spectra formed.

Table 1 displays comparison 2σ upper limits, demonstrating broad consistency between the pipelines. The RTS+CHIPS best limits appear at a lower k bin than FHD+ppsilon, owing to the different structure of the foreground leakage for the two pipelines. Figure 14, which shows the resulting FHD+ppsilon slices in ${k}_{\parallel }$ and ${k}_{\perp }$, also displays the equivalent power profiles for RTS+CHIPS. The comparison pipeline shows significantly more power at low ${k}_{\perp }$, consistent with no diffuse model being subtracted from the data. Aside from this area, the two pipelines show consistent results. Figure 16 compares the 1D power from both pipelines for the three central redshifts studied in this work. The magenta line reproduces the FHD+ppsilon results from Figure 15, while the gray line plots the equivalent RTS+CHIPS results (with coarse-band harmonics included). Despite the differences in the two pipelines, they produce consistent power across a range of scales, supporting the robustness of the results.

Zoom In Zoom Out Reset image size

Through the reduction of the first season of MWA high-band EoR observations, we have identified several opportunities to improve our analysis. Here, we compile a list of future improvements and their potential impact on the power-spectrum results. Broadly speaking, these improvements emphasize a need to better understand the foregrounds and the instrument, and are likely to be relevant to imaging based analyses in general.

Extended point-source model. The hierarchical catalog used here was a definite improvement over the MWACS catalog alone. This point-source model will continue to improve with the release of the GLEAM Survey, as well as dedicated MWA observations to target sources in the sidelobes of the EoR fields. Improving the foreground catalog in our sidelobes will be crucial in controlling the foreground leakage into the EoR window (Pober et al. 2016). Incorporating these improved catalogs into our analysis will improve both the foreground subtraction and the calibration.

Diffuse emission model. The diffuse model here contained no spectral information and represented only the Stokes I polarization. Nevertheless, subtracting this model removed 70% of the residual power in the foreground wedge. Adding a spectral index and multiple polarization components to the model will improve this subtraction even further. As Lenc et al. (2016) demonstrated, substantial polarized diffuse structure exists in the EoR0 field. The contributions to Stokes Q and U depend on the ionospheric conditions, which will be accounted for in future analyses to appropriately model this emission. In addition, the diffuse model will be extended to an all-sky model. It is evident from the jackknives shown in Figure 10 that strong emission near the horizon is leaking into our power spectra during off-zenith observations. By modeling and subtracting this emission we may be able to recover these observations in future analyses.

Primary beam model. In order to properly model the instrument response to the sky we also need an accurate primary beam model. This is especially important for EoR power spectrum measurements because the (often bright) near-horizon emission resides at the edge of the foreground wedge (Thyagarajan et al. 2015a, 2015b; Pober et al. 2016), but the low elevation beam can be difficult to model at the necessary precision (Thyagarajan et al. 2016). Asad et al. (2016) demonstrated the need for a polarized primary beam model in order to mitigate the risk of diffuse polarized emission leaking into Stokes I, conflating with a potential cosmological measurement. The MWA primary beam model is continuously improving to incorporate mutual coupling between the phased dipoles, and embedded element patterns for each of the 32 dipoles on each antenna (Sutinjo et al. 2015). Other methods are also being pursued to directly measure the full beam response in situ using the ORBCOMM satellite constellation (Neben et al. 2015), and drone-flown calibrator sources (D.C. Jacobs et al. 2016, in preparation). The latter method is similar to the strategy Virone et al. (2014) used to verify SKA beam patterns.

Calibration. The most limiting factor in our analysis is likely to be insufficiently precise calibration solutions, which can cause foreground power to leak into the EoR window, as we see in the integrated 2D power spectra (Figure 13). Because our calibration is based on sky and instrument models, the improvements described above will improve our calibration. We also plan to combine many snapshot observations when estimating gains and achieving higher signal to noise, which in turn will enable more sophisticated parameterizations of the antenna bandpasses. This will require an overhaul of the data flow depicted in Figure 1 because snapshots will need to be combined in the calibration step. Preliminary investigations have shown that the gains are stable enough to be fit across several hours of observation.

Between 2013 and the time of writing, the MWA has collected over 2,000 hours of data targeting the EoR fields. These observations include the three designated sky fields and two frequency bands. In addition, observations are underway to further characterize the sources in the instrument sidelobes. The second phase of the MWA will be commissioned in the latter half of 2016, and will include two sets of highly redundant cores of antennas. This hybrid configuration will enable a unique opportunity to study the calibration techniques used in imaging and delay spectrum analyses, providing another view of the systematics limiting our current results.

While the analysis presented here has not reached full potential, it represents the deepest power spectrum integration to date produced by an imaging pipeline. Imaging involves many difficulties (e.g., efficient gridding and mapmaking, foreground modeling), but if systematics can be overcome, it has the potential to compete with other more targeted experiments and analysis styles. More broadly, imaging analyses will be necessary to perform cross-correlation studies with complimentary probes such as galaxy surveys or other intensity mapping experiments (e.g., Lidz et al. 2009; Doré et al. 2014; Beardsley et al. 2015; Silva et al. 2015; DeBoer et al. 2016; Vrbanec et al. 2016), which will ultimately unlock the full potential of 21 cm observations. Here, we have identified a crucial region of power spectrum space that is currently contaminated, with suggestions for improvement in future analysis. With improved calibration techniques, primary beam models, and understanding of foregrounds, the MWA and other imaging analyses will be able to quickly approach the most competitive upper limits in the field.