THE MURCHISON WIDEFIELD ARRAY 21 cm POWER SPECTRUM ANALYSIS METHODOLOGY

The MWA is an interferometric array of phased array tiles operating in the 80–300 MHz radio band. Each tile consists of a 4 × 4 grid of dual polarization bow-tie shaped dipoles that are used to form a beam on the sky with a full width of 26°($\lambda /2$ m) at the half power point. Signals from individual antennas are summed by an analog delay-line beamformer which can steer the beam in steps of 68 $\mathrm{cos}({\rm{elevation}})$. The signal is digitized over the entire bandwidth but only 30.72 MHz are available at any one time. This 30 MHz of bandwidth is broken into 1.28 MHz "coarse" bands by a polyphase filter-bank in the field and sent to the correlator (Ord et al. 2015), where it is further channelized to 40 kHz, cross-multiplied and then averaged at 0.5 s intervals. More details on the design and operation of the MWA can be found in Lonsdale et al. (2009) and Tingay et al. (2013).

The MWA reionization observing scheme spans two 30 MHz tunings, 140–170 MHz (9.2 < z < 7.5) and 167–196 MHz (7.5 < z < 6.25) and two primary minimal foreground regions (Field 0 at R.A. 0^h, −27°and Field 1 at 4^h, decl. −27°); both transit the zenith at the MWA's latitude and are near the south galactic pole. A third pointing toward Hydra A is also observed; see Figure 1 for an overview. Here we focus on the low redshift tuning, and the R.A. = 0^h pointing, where the band is chosen for its lower sky temperature and pointing is chosen for its ease of calibration—having fewer bright, resolved sources; see Table 1 for a listing of observing parameters.

Zoom In Zoom Out Reset image size

Table 1.  MWA EoR Observing Parameters

Parameter	Value
Field of View	26°$(\lambda /2){\rm{m}}$ FWHM
Tuning	166–196 MHz redshift range $7.56\lt z\lt 6.25$
Target area	(R.A., decl.) 0h00m, −27°00m
Primary beam pointing quantization	68
Snapshot length	112 s
Time and frequency resolution	0.5 s, 40 kHz
Post-flagging resolution	2 s, 80 kHz
Time	3 hr on 2013 Aug 23, six 30 minute pointings or 96 snapshotsa

Note.

^aThe same data set is used in every pipeline run.

Download table as:  ASCIITypeset image

The analysis presented here is on three hours of data, one of 400 nights which have been collected as part of the observing program; 150 nights are thought to be necessary for a detection of typical models (Beardsley et al. 2013).

During observation, the beam-former was set such that the target region repeatedly drifted through the field of view. With an available beamformer step size of 68; each drift was about 30 minutes long. This was done for a total of 6 pointings in a night, or about 3 hr. The data included here include the two pointings leading up to the target crossing zenith, the zenith pointing, and then three more pointings after the transit crossing. Data were recorded in 112 s units (which we term "snapshots") for a total of 96 snapshots. These snapshots are the basic unit of time on which many operations become independent—e.g., RFI flagging, calibration, and imaging. Though the full set of linear polarization parameters is correlated, and Stokes I images and power spectra are the final product of interest, at this stage of the analysis the instrumental polarizations have been found to be more instructive; with one exception, only the linear east–west polarization is examined here. No significant differences are seen in the north–south data. The same set of snapshots is used in every pipeline run.

Each snapshot is flagged for interference using the AOFlagger (Offringa et al. 2010)²⁸ algorithm, which applies a high pass filter to remove foregrounds, uses the SumThreshold algorithm to search for line-shaped outliers, and then applies a scale invariant rank operator to the resulting mask, which identifies certain morphological features.

As described in Offringa et al. (2015), the interference environment at the Murchison Radio-astronomy Observatory is relatively benign and generally requires flagging of about 1.1% of the data. After flagging the data are averaged to 2 s and 80 kHz, reducing the data volume by a factor of 8.

The MWA RTS (Mitchell et al. 2008; Ord et al. 2010) was initially designed to make wide-field images in real time from the MWA 512-tile system (Mitchell et al. 2008). On the de-scoped 128 element array, it has been implemented as an offline system, where it has been adjusted to compensate for the lower filling factor (Ord et al. 2010). The RTS incorporates algorithms intended to address a number of known challenges inherent to processing MWA data, including; wide-field imaging effects, direction-dependent (DD) antenna gains and polarization response, and ionospheric refraction of low-frequency radio waves. Each MWA observation (112 s) is processed by the RTS separately, in series. The RTS is also parallelized over frequency so that each coarse channel (1.28 MHz broken into 40 kHz channels) is processed largely independent of the other coarse channels, with only information about peeled source offsets communicated between processing nodes.

The RTS calibration strategy is based upon the "peeling" technique proposed by Noordam (2004) and a foreground model using a cross-matching of heritage southern sky catalogs²⁹ with the MWA Commissioning Survey.³⁰ The brightest apparent calibrators in the field of view are sequentially and iteratively processed through a Calibrator Measurement Loop (CML). During each pass through the CML, (i) the expected (model) visibilities of known catalog sources are subtracted from the observed visibilities. For the data processed in this work, 1000 sources are subtracted for each observation. (ii) The model visibilites for the targeted source are added back in and phased to the catalog source location. Any ionospheric offset of the source can now be measured by fitting a phase ramp to the phased visibilities. (iii) The strongest sources are now used to update the DD antenna gain terms, while weaker sources are only corrected for ionospheric offsets. For this work, 5 sources are used as full DD calibrators and 1000 sources are set as ionospheric calibrators. The CML is repeated until the gain and ionospheric fits converge to stable values. A single bandpass for each tile is found by fitting a 2nd order polynomial to each 1.28 MHz wide coarse channel. The ∼1000 strongest sources are then subtracted from the calibrated visibilities. Calibration and model subtraction parameters are summarized in Table 3. Model subtracted visibilities are passed to the RTS imager and to the CHIPS power spectrum estimator.

The RTS imager uses a snapshot imaging approach to mitigate wide-field and DD polarization effects. Following calibration, the residual visibilities are first gridded to form instrumental polarization images which are co-planar with the array. These images are then regridded into the HEALpix (Górski et al. 2005) frame with wide-field corrections. Weighted instrument polarization images are stored, along with weight images containing the Mueller matrix terms, so that further integration can be done outside of the RTS. It is also possible to use the fitted ionospheric calibrator offsets to apply a correction for ionospheric effects across the field during the regridding step or subtraction of catalog sources, but in this work this correction has not been applied. These snapshot data and weight cubes are then integrated in time to produce a single HEALpix cube. This cube, averaged over the spectrum, is shown, with and without foregrounds, in Figure 3.

Zoom In Zoom Out Reset image size

FHD (Sullivan et al. 2012) is an imaging algorithm designed for very wide field of view interferometers with direction- and antenna-dependent beam patterns. Simulated beam patterns are used to grid visibilities to the uv plane and the reverse operation of turning gridded cubes into time-ordered visibilities. This simulation of weights provides a necessary accounting of information loss caused by the inherent size of the dipole element. As it is designed from the ground up to account for the widefield effects inherent in the MWA it has grown to include a full range of tasks such as calibration and simulation. In this pipeline the framework is used to both calibrate and image the data.

The FHD calibration pipeline generates a model data set, computes a calibration solution which minimizes the difference with the data, smooths the calibration solution to minimize the number of free parameters, and outputs the residual. The calibration model is formed from sources found by deconvolving, in broadband images, about 75 of the 96 snapshots included here and retains those which are common to all snapshots and pass other consistency checks (P. Carroll et al. 2016, in preparation). In each snapshot sources are included in the model if they are at or above 1% of the peak primary beam, this amounts to about 7000 sources and a flux limit of about 80 mJy with slight variations snapshot to snapshot. Most sources in this catalog have spectral indexes between 0 and −2, with the majority near a mean of −0.8, which corresponds to a 13% difference across the 30 MHz. Spectral index is not directly modeled; spectra are simulated to be flat; however, during catalog subtraction, the data are multiplied by a positive spectral index of 0.8 such that most sources will appear to be flat. Full spectral modeling of catalog sources will be included in future analyses.

The function of FHD's calibration is to minimize the number of free parameters, with the twin goals of minimizing potential signal loss and building a deep understanding of instrumental systematics. Here we give a brief description of the instrument model; a listing of all the parameters mentioned here is also given in Table 3 along with a rough accounting for the total number of fitting parameters. Initial complex gain solutions are computed using the Alternating Direction Implicit technique described in Salvini & Wijnholds (2014) for each antenna, channel, and polarization. This generates a gain and phase for every channel on every tile, for each 112 s snapshot. Most antennas have similar solutions with the main features corresponding to the exact type and length of analog cable feed of which there are 6 different types; per-tile solutions are further averaged into per-cable-type and averaged over the entire 3 hr observation. After these solutions are divided out, the residual per-tile solutions are further fit for a second order amplitude spectral polynomial and a 1st order phase slope. This is done on every snapshot to account for temperature-driven amplifier gain changes. One systematic easily visible in the power spectra is a small reflection corresponding to the 150 m cables. This is fit and removed as a phase delay with a ∼0.1 dB amplitude in the time averaged per-tile bandpass solutions.

The residual time-ordered visibilites are then passed to CHIPS and to FHD imaging for formation of spectral cubes.

FHD Imaging. The imager portion of the FHD framework is the FHD algorithm, which is based on loss-less optimal map-making developed for the CMB Tegmark (1997). It is a variant on the class of wide-field imaging algorithms such as A-projection (Bhatnagar et al. 2013), peeling (Mitchell et al. 2008), and forward modeling (Bernardi et al. 2011; Pindor et al. 2011) which improve upon standard radio interferometric algorithms by including the primary beam in the calculation of the instrumental response. It is holographic in the sense that it uses the known primary beam of the antenna when performing simulation and gridding operations, and fast due to its use of a sparse matrix representation of the instrument transfer function. FHD is described in detail by Sullivan et al. (2012). Deconvolution is not used in this pipeline paper, however, FHD's relevant focus on highly accurate primary beam modeling is used to form optimally weighted images. FHD does not currently grid w terms separately and is therefore limited to planar arrays and relatively short periods of rotation synthesis where long baselines do not rotate significantly with the Earth.

As deployed in the pipeline described here, the FHD imager is used to grid the time ordered visibilities into a uvf cube weighted according to the Fourier transform of the primary beam using an electromagnetic model described by Sutinjo et al. (2015). This weighting scheme is similar in effect to "natural" weighting scheme, which provides the highest possible supernova remnant by weighting uv samples according to the inverse variance. Natural weighting is traditionally not favored in interferometric imaging with sparse arrays, as it can dramatically impact the point-spread function by downweighting baselines for which there are few samples. We choose it here because it maximizes the sensitivity of our densely sampled uv plane core.

Deep mosaicks are made by using the warped snapshot method (Cornwell et al. 2012). The data are split into even and odd samplings at the 8 s cadence and then imaged at a 2 min cadence with a uv resolution chosen to result in a 90 degree field of view. The even odd split is carried through to the final power spectrum analysis where they are cross-multiplied. In each uv-frequency (uvf) voxel we accumulate weighted data, weights, and the square of the weights according to

$$\begin{eqnarray}{D}_{{uvf}} & = & \displaystyle \sum _{i}\;{B}_{i;{uvf}}{V}_{i;{uvf}}\\ {W}_{{uvf}} & = & \displaystyle \sum _{i}\;{B}_{i;{uvf}}\\ {{\rm{Var}}}_{{uvf}} & = & \displaystyle \sum _{i}\;{B}_{i;{uvf}}{B}_{i;{uvf}}^{* },\end{eqnarray} \tag{ 1 }$$

where B_i is the Fourier transform of the cross power beam evaluated in that uvf voxel at snapshot time i,V_i is visibility at time i.³¹ Essentially D is numerator of the mean performed in the mosaicking step with W the normalizing denominator of the mean; the same is done for the error cube Var/W².³² The cubes are Fourier transformed, corrected for the w projection coordinate warping, and gridded into the HEALpix frame.

Mosaicking. The 112 snapshot HEALpix cubes are summed in time, keeping pixels with a beam weight of 1% or more, a cut which effectively limits the field of view to ∼20°. The resulting mosaick is handed on to εppsilon (Section 3.4) and EmpCov (Section 3.6). This image, averaged over the spectral dimension, is shown, with and without foregrounds, in Figure 3. This image retains the full weighting proportional to the number of samples in each uvf cell and is therefore very similar to natural weighting. Though beam weighting theoretically gives an optimal inverse variance weighting for each snapshot it does not capture the change in variance due to changing system temperature. The mosaicking as performed here weights each snapshot equally, under the assumption that the system temperature, which is dominated by the temperature of the sky in the direction of the phased array pointing, stays roughly constant through the three hour tracked integration.

Through the parallel-but-convergent development of these imagers have emerged two very similar systems; however, some differences remain in the analysis captured here. The two primary differences are in the treatment of calibration and in the subtracted catalogs.

In both pipelines the calibration is a two-step process. First, calibration solutions for each channel, and antenna are computed by solving for the least-squares difference with a model data set. Next, those solutions are fit to a model of the array; for example, fitting a polynomial to the bandpass. FHD and RTS take different approaches to this step, a fact reflected in the number of free parameters in this fit. A smaller number of parameters minimizes the possibility of cosmological signal loss or the spurious incorporation of unmodeled foreground emission in the calibration solutions (Barry et al. 2016); more free parameters can absorb physics missing from the instrument model. As tabulated in Table 3 the RTS fits for ∼6420 free parameters while FHD fits for ∼880.

In practice, some parameters will be averaged over more than a single night, which will further reduce the number of free parameters per observation, though hundreds to thousands of free parameters is still typical. This is a large number but it is considerably smaller than the 180 million data points typically recorded in single snapshot obdservation.

As has been noted, there is nearly an order of magnitude difference in the number of free parameters between the two pipelines, which is worth considering. The primary difference is in the treatment of the passband. There are a number effects which show up in the passband calibration. The edges of the 1.28 MHz bands are known to be subject to aliasing from adjacent coarse channels as well as under-sampling when cast to 4 bit integers by the correlator (van Vleck corrections) and so are flagged. This flagging creates a regular sampling function which shows up as the characteristic horizontal lines in a 2D power spectrum (see Section 4.1). Added to this is a small amount of interference flagging. Additionally, reflections at analog cable junctions show up as additional spectral ripple corresponding to the length of the cables.

The RTS fits for a low order polynomial on every 1.28 MHz chunk on every antenna, while FHD averages each channel over all antennas to get a common passband for all and then fits a low order polynomial to get any tile to tile variation. This significantly reduces the number of free parameters and the likelihood of signal loss, though leaving open the possibility of additional un-modeled instrumental effects.

The construction of the foreground subtraction model is also a point of difference between the two pipelines. As noted in Table 3, foreground/calibration models contain different numbers of sources which have been derived by different means. The RTS catalog cross-matches multiple heritage southern sky catalogs with the MWA Commissioning Survey using the Bayesian cross-matcher PUMA (J. Line et al. 2016, in Preparation). The FHD subtraction model contains sources found in a deep deconvolution of this same data set. Both catalogs have the goal of producing a reliable set of sources that minimizes false positives and accurately reflects resolved components, though they go about it in different ways. The FHD catalog focuses on the reliability aspect by performing a deconvolution on every snapshot used in the observation and selecting sources which appear in most observations (P. Carroll et al. 2016, in preparation). The RTS catalog has used the somewhat less precise MWA commissioning catalog but by cross-matching these sources against many other catalogs of known sources and fitting improved positions and fluxes, the accuracy is seen to increase.

εppsilon calculates a power spectrum estimate from image cubes and directly propagates errors through the full analysis, see B. Hazelton et al. (2016, in preparation) for a full description. The design criteria for this method is to make a relatively quick and uncomplicated estimate of the power spectrum to provide a quick turnaround diagnostic.

The accumulated data, weight, and variance cubes are Fourier transformed along the two spatial dimensions into uvf space, where the spatial covariance matrix is assumed to be diagonal. This is approximately true if the uv pixel size is well matched to the primary beam size, so the εppsilon direct Fourier transform grid size is restricted to being equal to the width of the Fourier transformed primary beam; i.e., 1/(field of view). The uvf data cubes are then divided by the weight cubes to arrive at a uniformly weighted 3D Fourier cube. The variance cubes are similarly divided by the square of the weights. Next, the sums and differences of the even and odd cubes are computed with variances given by adding the reciprocal of the even and odd variances in quadrature. The difference cube then contains only noise and the sum cube contains both sky signal and noise.

The next step is to Fourier transform in the frequency direction. Here we choose to weight by a Blackman–Harris window function, which heavily downweights the outer half of the band and decreases the leakage of bright foreground modes into other power spectrum modes as described in Thyagarajan et al. (2013), Parsons et al. (2012), and Vedantham et al. (2012), among others. The transform of the spectral dimension is done using the Lomb & Scargle periodogram to minimize the effects of regular gaps in the spectrum, which occur every 1.28 MHz. The sky signal power is estimated by the square of the sum cube minus the square of the difference cube, which is mathematically identical to the even/odd cross power if the even and odd variances are identical

$$\begin{eqnarray}&&{({(E+O)(E+O)}^{\dagger }-(E-O)(E-O)}^{\dagger })\displaystyle \frac{1}{4}\\ &&=({{EE}}^{\dagger }+{{OO}}^{\dagger }+{{EO}}^{\dagger }+{{OE}}^{\dagger }\\ &&\quad -({{EE}}^{\dagger }+{{OO}}^{\dagger }-{{EO}}^{\dagger }-{{OE}}^{\dagger }))\displaystyle \frac{1}{4}\\ &&=(2{{EO}}^{\dagger }+2{{OE}}^{\dagger })\displaystyle \frac{1}{4}\\ &&={\mathfrak{R}}({{EO}}^{\dagger }),\end{eqnarray} \tag{ 2 }$$

while the square of the difference cube provides a realization of the noise power spectrum.

Diagnostic power spectra are generated by averaging cylindrically to a two-dimensional ${k}_{\parallel },{k}_{\perp }$ power spectrum. The diagnostic power spectra are computed over the full 30 MHz bandwidth to provide the highest possible ${k}_{\parallel }$ resolution of the foreground and systematic effects. These are shown in the left column of Figure 4. One-dimensional power spectra (shown in Figure 7), are calculated by a similar process but only using 10 MHz³³ of bandwidth, which corresponds to a redshift range of 0.3 (at 182 MHz) and averaging along shells of constant k, masking points within the foreground wedge³⁴ and weighting by inverse variance.

Zoom In Zoom Out Reset image size

The error bars are estimated as the sum of the beam weights accumulated in the mapping process (e.g., Equation (1)) assuming a constant noise figure for each data point

$$\begin{eqnarray}&&{\sigma }_{{uvf}}^{2}=\displaystyle \sum _{i}\;{{Var}}_{{uvf}}^{i}{/({W}_{{uvf}}^{i})}^{2}{\sigma }_{{if}}^{2},\end{eqnarray} \tag{ 3 }$$

where the noise ${\sigma }_{{if}}^{2}$ is an average noise spectrum, per snapshot i, calculated by differencing even–odd data sets and computing an rms over all baselines for each snapshot i. These errors are then propagated into the 2D and 1D power spectra by quadratic sum.

As noted above, this noise model assumes that the apparent sky and receiver temperatures are roughly constant through the 3 hr tracked observing period. Another way of estimating error is to form power spectra from the even–odd difference. Comparing these two methods, we find that they agree to within a few percent in both 2D and 1D power spectra.

The CHIPS power spectrum estimation method computes the maximum likelihood estimate of the 21 cm power spectrum using an optimal quadratic estimator formalism and is more completely described in Trott et al. (2016). The design criteria for this method were to fully account for instrumental and foreground induced covariance in the estimation of the power spectrum. The approach is similar to that used by Liu & Tegmark (2011), but with the key difference of being performed entirely in uvw-space, where the data covariance matrix is simpler (block diagonal), and feasible to invert. This approach also allows straightforward estimation of the variances and covariances between sky modes by direct propagation of errors. CHIPS takes as input calibrated and foreground subtracted time-ordered visibilities. Tapping into the pipeline post-calibration but before imaging, CHIPS uses its own internal instrument model to estimate and propagate uncertainty.

The method involves four major steps: (1) Grid and weight time-ordered visibility channels onto a uvw-cube using the primary beam model, (2) compute the least squares spectral transform along the frequency dimension to obtain the best estimate of the line-of-sight spatial sky modes (this technique is comparable to that used by εppsilon), (3) compute the maximum-likelihood estimate of the power spectrum, incorporating foregrounds and radiometric noise, averaging k_x and k_y modes into annular modes on the sky, ${k}_{\perp };$ (4) compute the uncertainties and covariances between power estimates. The first step is the most computationally intensive, requiring processing of all the measured data. The main departure point for CHIPS from εppsilon is in the much finer resolution of the uv grid. Using an instrument model, CHIPS calculates the covariance between uv samples as a function of frequency. Since the beam and uv sampling function are both highly chromatic, extra precision in this inversion is thought to be highly beneficial. After a line of sight transform similar to that used by εppsilon, this covariance information is inverted to find the Fisher Information, the maximum likelihood power spectrum, and covariances between measurements. The maximum likelihood estimate of the power in each ${k}_{\perp },{k}_{\parallel }$ mode is shown in the right column of Figure 4 and averaged in spherical bins in Figure 7.

Before this last averaging step one can optionally include an additional weighting by the known power spectrum of a confused foreground in a process described in more detail for these data by Trott et al. (2016). A 2D foreground weighted power spectrum is shown in Figure 5. The power spectra shown in Figure 7 show an excess of power in excess of the expected noise. This excess is notably similar between both calibration/foreground subtraction pipelines. The amount of power in the excess, as compared with the error bars, also depends rather dramatically on the range of k bins included in the final averaging to the 1D. These are discussed in more detail in Section 4.1.

Zoom In Zoom Out Reset image size

The EmpCov power-spectrum estimation method computes both 2D and 1D power spectra using the quadratic estimator formalism. The method and its application to this data are described in more detail by Dillon et al. (2015a).

The quadratic estimator method of Liu & Tegmark (2011) treats foreground residuals in maps as a form of correlated noise and simultaneously downweights both noisy and foreground-dominated modes, keeping track of the extra variance they introduce into power spectrum estimates. This technique can be computationally demanding but using acceleration techniques described by Dillon et al. (2013), has been applied to the previous MWA 32T results of Dillon et al. (2014), while a very similar technique, working on visibilities rather than maps but also using the data itself to estimate covariance, was used for the recent PAPER 64 results of Ali et al. (2015). Dillon et al. (2015a) build on these methods to mitigate errors introduced by imperfect mapmaking and instrument modeling through empirical covariance estimation, assuming all data covariance is sourced by foregrounds.

EmpCov takes as input FHD calibrated images with foregrounds subtracted as well as possible, split into even and odd time-slices and averaged over many observations. The differences between the even and odd time-slices, which are assumed to be pure noise, are used to calibrate the system temperature in a noise model derived from observation time in uv cells. In order to avoid directly propagating instrumental chromaticity into the foreground residual covariance models (Dillon et al. 2015b), EmpCov uses the data itself to properly downweight residual foregrounds as seen by the instrument. It does this by estimating the frequency–frequency foreground residual covariance in annuli in uvf space, assuming that different uv cells independently sample the foreground residuals. This assumption, similar to that made by CHIPS, allows the combined foreground and noise covariance to be inverted directly and used to downweight the cubes when binning into 2D (and eventually 1D) band powers. As part of the quadratic estimator formalism, EmpCovalso calculates error covariances and window functions (i.e., horizontal error bars). The resulting power spectrum is shown in Figures 5 and 7.

One benefit from having multiple pipelines is the freedom to investigate different optimization axes. The design of the εppsilon power spectrum estimator emphasizes speed and relative simplicity, choices motivated by the need to understand the effect, on the power spectrum, of processing decisions such as observation protocol, flagging, and calibration. Using εppsilon we have discovered and corrected multiple systematic effects, primarily those of a spectral nature which were not obvious in imaging but quite apparent in the 2D power spectrum. With the ability to quickly form power spectra on different sets of data, εppsilon has been an important tool for selecting sets of high quality data.

In contrast, CHIPS starts from time-ordered data and in its calculations emphasizes a more full accounting of instrumental and residual foreground covariance. Not only does this higher resolution covariance calculation provide a more accurate accounting of the instrumental window function on the power spectrum, but it also allows for more precise weighting schemes based on knowledge of the statistical properties of the residual foregrounds. This is useful when making 1D power spectra where foreground-like modes can be downweighted in the average.

Somewhere in the middle of these two is EmpCov, which, like εppsilon, uses image cubes and associated weighting variances but performs a more formal quadratic estimator in which additional covariances can be downweighted and the effects of the instrument window function be factored into the calculation of the power spectrum bins and error bars. It also demonstrates the impact of inverse covariance weighting by forming a measure of covariance from the data. This measure encapsulates both the residual foregrounds modeled by CHIPS as well as any other residuals resulting from mis-calibration.

We apply both εppsilon and the unweighted version of CHIPS to both of our calibration and foreground subtraction pipes to produce a total of four different power spectra (Figures 4 and 6). Each power spectrum estimator has been developed to target the output from a "primary" calibration and foreground subtraction process—the diagonal panels of Figure 4—and have been highly optimized to that up stream source of data. The off-diagonal power spectra were created using auxiliary links which import the data and the metadata produced by the foreground subtraction step. Since they are less highly optimized, lacking as they do the advantage of a close working relationship, these pathways represent an upper limit on the variance to be expected from small analysis differences but allow us to look for effects common to foreground subtraction or to the power spectrum method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Properties shared by all are the large amount of power at low k_∥ roughly at an amplitude of 10¹⁵ mK²/(Mpc/h)³. This emission is approximately flat over most of k_⊥ but rises steeply rise at low ${k}_{\perp }$. The amplitude agreement is particularly apparent in Figure 6 where we plot a slice of the 2D power spectrum at ${k}_{\parallel }=0$ where most foreground power is expected to lie. A model of smooth galactic emission has not been subtracted, which likely contributes to this steep rise. The "wedge" shaped linear dependence on baseline length in the 2D power spectra is due to the inherently chromatic response of a wide field instrument to smooth spectrum foregrounds; sources entering far from the phase center appear as bright pixels at higher ${k}_{\parallel }$ with sources on the horizon at the edge indicated by Figure 4's solid black line. The solid and dotted lines in the figure indicate the upper boundaries of power from sources at the horizon and at the beam half power point, respectively. With the exception of some instrumental features foreground power is well isolated within this expected boundary. This emission is also visible in the image cubes as side-lobes extending from outside of the imaged area which move as function of frequency. Observations recorded when the Galactic plane is near the horizon have a much larger wedge component and have been excluded from this analysis. See Thyagarajan et al. (2015a, 2015b) and Pober et al. (2016) for a detailed discussion of the foreground contributions to the power spectra in this data.

The two main instrumental systematics are horizontal striping due to missing or poorly calibrated data at the edges of regular coarse passbands and vertical striping due to spectral variation near uneven uvf sampling. As described in Section 2, the MWA reads out spectra which are divided into 1.28 MHz wide "coarse bands." These bands have small known aliasing and gain instabilities in the edge channels, so during initial flagging we flag the edge-most 85 kHz. This regular gap in the spectra corresponds to a poor sampling in Fourier space at integer values of $\eta =781$ ns or ${k}_{\parallel }=0.45$ h Mpc⁻¹. In the 2D power spectrum this manifests as horizontal stripes of high power which are in fact sidelobes of the foreground wedge, and higher error bars reflecting a lack of information about these modes. These sidelobes can be minimized by working to improve the accuracy of the passband calibration and thus having fewer flagged channels. They can also be downweighted by accounting for the covariance between modes as is done by EmpCov.

A similar issue also arises from gaps in the uv coverage. In places where coverage is not uniform—such as at longer baseline lengths where baseline density drops as $1/{r}^{2}$—beam weights can vary dramatically as a function of frequency leading to a vertical striping effect. It is most prominent above $| u| \gt 100\lambda$ or ${k}_{\perp }\gt 0.1$ h Mpc⁻¹, which corresponds to when uv sampling begins to drop below unity beyond the densest part of the MWA core. This effect can be ameliorated by using an accurate model for the beam and array positions to grid these samples according to the optimal mapmaking procedure, the approach taken by the FHD imager, or the CHIPS approach of calculating and inverting the full instrumental covariance.

The most noticeable difference between the different pipeline paths is in the power level in the window above the horizon and below the first coarse passband line (between 0.1 and 0.3 ${k}_{\parallel }$ and 0.01 and 0.05 ${k}_{\perp }$). FHD to εppsilon displays a noise-like window in the 2D space, with a number of points dipping below zero while the other methods are noise-like only at at much higher ks. One commonality between all power spectra with this positive bias is a relatively higher amplitude of the coarse passband lines.

The relative amplitude of the vertical striping is probably the largest difference between the four power spectrum methods. FHD-εppsilon sees vertical striping largely consistent with noise, and the other methods see the striping at varying levels with both CHIPS spectra showing the largest. As discussed in detail by Morales et al. (2012) and shown in data by Pober et al. (2016), this vertical striping is very sensitive to the accuracy of the weights used to average multiple samples together. εppsilon relies on the imager (FHD or RTS) to simulate the instrument and generate optimally weighted maps, while CHIPS uses an internal instrument model to calculate covariance. FHD used the second generation Sutinjo et al. (2015) beam model, which takes into account cross-coupling within a tile, while the rest used an analytic short-dipole approximation.

We also compare the results of CHIPS and EmpCov using analogous foreground downweighting schemes. A quantitative comparison of these power spectra is difficult, since the quadratic estimator's downweighting scheme does not preserve foreground power. However, the results in Figure 5 are largely similar, showing the familiar wedge structure and the brightest foreground contamination at low ${k}_{\perp }$ where galactic foregrounds dominate. EmpCov excludes long baselines where coverage and sensitivity is poor and as such does not probe to the same range in ${k}_{\perp }$ as CHIPS. EmpCov appears to more successfully remove foreground contamination near the wedge, which likely means that the foreground models employed in CHIPS have room for improvement. Likewise, EmpCov can successfully remove the lines in constant ${k}_{\parallel }$ that arise from flagged channels due to the MWA's coarse band structure, but is still contaminated by the 90 m cable reflection at ${k}_{\parallel }\approx 0.45$ h Mpc⁻¹ (Dillon et al. 2015a).

The final analysis step is to average into 1D power spectra along shells of constant k. These are shown in Figure 7 for three of the four analysis tracks³⁵ shown in Figure 4 with the addition of the Dillon et al. (2015a) points and a theoretical sensitivity curve calculated using the 21CMSENSE sensitivity code³⁶ by Pober et al. (2014).

The positive biases visible in the 2D power spectra are also apparent here. Only a few points are fully consistent with zero at 2σ; however, most are very close to the theoretical sensitivity curve and have errors matching those predicted for noise. The power spectra fall into two groups, those calculated from input image cubes (εppsilon and EmpCov) and those calculated directly from visibilities using CHIPS. The image-based points are somewhat deeper at low k, as noted in the 2D plots. Points from CHIPS are biased more strongly at low k but the slope is flatter and converges with the other pipelines at higher k.

Part of the CHIPS bias is due to the calculation of weightings. Default CHIPS analysis uses a statistical model of confused foregrounds to downweight biased modes, particularly those correlated with the wedge power. For this reason it is desirable to include the wedge modes in that 1D average. However, it significantly changes the interpretation of error bars; points in which a significant amount of power have been downweighted will have error bars much larger than thermal. In the interest of comparing with the other methods, the power spectra in Figures 4, 6, and 7 have been calculated using points only lying outside the wedge horizon. This limits the amount of wedge-to-window covariance CHIPS can remove and contributes to the larger bias.

Including the full wedge in the CHIPS covariance calculation offers foreground suppression, but also introduces a foreground component into the error bars. Compare in Figure 8 the RTS→CHIPS power spectrum in Figure 7 with that given by Trott et al. (2016), which used the same data shown here, though only 1/3 of the 30 MHz band. In both, CHIPS has downweighted by a model of foreground covariance formed by propagating a statistical model of confused sources. The only difference is that black excludes the wedge but red does not. When the wedge is included the modeled foregrounds in those voxels dominate the covariance weights. Applying these weights essentially moves the foreground bias into the error bars and asserts that, given our best model of foregrounds, the power spectra are completely consistent with noise and foregrounds.

Zoom In Zoom Out Reset image size

The power spectra in Figure 7 show the range of results possible given the same input data. Though they do not all agree, they do paint a consistent picture. Differences partly come from the definition of error bars but also indicate the relative difficulty of methods. Methods which rely on an imager seem to perform somewhat better. This is perhaps unsurprising. CHIPS computes the instrumental correlation matrix in visibility space using beam, bandpass, etc. As the CHIPS analysis exists entirely in the visibility space, errors in modeling the instrument are perhaps more difficult to detect than they are in the image space. However, we do not suggest that visibility-based calculations like CHIPS are doomed to failure; rather the opposite. The instrument models will continue to improve, and this improvement will be easily validated by comparison with the other pipelines.