CHIPS: THE COSMOLOGICAL H i POWER SPECTRUM ESTIMATOR

The power spectral density (power spectrum) measures the spatial covariance of a signal, integrated over a spatial volume, and corresponds to the Fourier transform of the two-point correlation function (i.e., the autocorrelation function), $\xi (r)$:

$$\begin{eqnarray}P({\boldsymbol{k}}) & = & {\displaystyle \int }_{V}\xi ({\boldsymbol{r}})\mathrm{exp}(-2\pi i{\boldsymbol{k}}\cdot {\boldsymbol{r}})d{\boldsymbol{r}},\\ & = & {\displaystyle \int }_{V}\langle T({{\boldsymbol{r}}}_{1})T({{\boldsymbol{r}}}_{1}+{\boldsymbol{r}}){\rangle }_{{{\boldsymbol{r}}}_{1}}\mathrm{exp}(-2\pi i{\boldsymbol{k}}\cdot {\boldsymbol{r}})d{\boldsymbol{r}},\\ & = & \displaystyle \frac{1}{V}{\displaystyle \int }_{V}T({{\boldsymbol{r}}}_{1})T({{\boldsymbol{r}}}_{1}+{\boldsymbol{r}})\mathrm{exp}(-2\pi i{\boldsymbol{k}}\cdot {\boldsymbol{r}})d{\boldsymbol{r}}d{{\boldsymbol{r}}}_{1},\\ & = & \displaystyle \frac{1}{V}\tilde{T}({\boldsymbol{k}}){\tilde{T}}^{*}({\boldsymbol{k}}),\end{eqnarray} \tag{ 2 }$$

where $T({\boldsymbol{r}})$ is the temperature fluctuation (relative to the mean) at vector position ${\boldsymbol{r}}$, and $\langle \rangle$ denotes an average over positions in the volume, V, as a proxy for an ensemble sampling of the distribution. The power spectrum has the units of temperature squared times volume (${{\rm{K}}}^{2}\;{{\rm{Mpc}}}^{3}$) and describes the integrated power on a given spatial scale, ${\boldsymbol{k}}$, averaged over the volume, V.

The power spectrum can be computed directly from an image cube, using a three-dimensional Fourier transform of image intensities to find $\tilde{S}({\boldsymbol{k}})$, and then squaring and normalizing by the cube volume, where $S({\boldsymbol{x}})$ is the measured brightness temperature in units of Jy beam⁻¹. The conversion from brightness temperature to flux density is given by

$$\begin{eqnarray}&&{S}_{\mathrm{Jy}}={10}^{26}\displaystyle \frac{2{{kT}}_{K}}{{\lambda }^{2}}\;{\rm{\Omega }}\;{\rm{Jy}},\end{eqnarray} \tag{ 3 }$$

which is the product of the specific intensity and the observation solid angle, Ω, and k is the Boltzmann constant.

Alternatively, the power spectrum can be computed directly from observed interferometric visibilities. The visibility is the mutual coherence of the two electrical signals detected by the antennas forming a baseline (Thompson et al. 1986). In the flat-sky approximation, where the field of view (FOV) of the instrument is small such that sky curvature can be ignored, the visibility is identically the Fourier transform of the product of the sky signal and the beam response. For wide FOV instruments, such as the MWA, this approximation breaks down, and there is a curvature convolving kernel in the measured visibility, in addition to the sky signal and beam response. Thyagarajan et al. (2015a, 2015b) studied this effect with MWA data, but current estimators, in general, ignore this effect. The visibility is the natural measurement space of the instrument (the radiometric noise is uncorrelated between visibilities) and will be used in CHIPS to compute the power spectrum directly (without tranforming to and from image space). Curvature sky terms are handled explicitly in CHIPS, and this design choice offers a natural departure from other, image-based power spectrum estimators (Liu & Tegmark 2011; Patil et al. 2014; Paul et al. 2014; Dillon et al. 2015; B. Hazelton et al. 2015, in preparation), but has been used for the angular power spectrum (Choudhuri et al. 2014). A related estimator, the "delay spectrum" (Parsons et al. 2012, 2014, and references therein), directly Fourier transforms along each visibility's frequency channels, matching a temporal delay with the position of given sky emission relative to the phase center. This visibility-based estimator approximates the power spectrum on large angular scales, but breaks down for the longer baselines.

The spherically-averaged (1D) power spectrum computes the power on three-dimensional spatial scales, $k=| {\boldsymbol{k}}|$, under the assumption that the statistics of the signal are isotropic and translationally invariant (the latter being a fundamental assumption of the power spectrum, according to the Wiener–Khinchin theorem). The 1D power spectrum has the advantage of integrating over the largest parameter space, potentially increasing signal detectability. The dimensionless 1D power spectrum, which integrates the total power on a given spatial scale over the volume, is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}(k)=\displaystyle \frac{{k}^{3}}{2{\pi }^{2}}P(k).\end{eqnarray} \tag{ 4 }$$

PAPER has published the most competitive 1D limits to date, relative to theoretical expectations, yielding (22.4 mK)² at z = 8.4 and $0.15\;h$ Mpc${}^{-1}\lt k\lt 0.5\;h$ Mpc⁻¹ (Ali et al. 2015). Patil et al. (2014) have demonstrated the power of the LOFAR variance statistic, which computes the overall variance in each spectral channel, with simulated data.

It is advantageous, however, to compute the 2D power spectrum as a first step, to discriminate between continuum contaminating foreground sources and the cosmological spectral-line signal. The structure of foregrounds in 2D space is described in Section 5.1; however, we comment here that it is substantially different from the cosmological signal. The 2D power spectrum has arguments (${k}_{\perp },{k}_{\parallel }$), where ${k}_{\perp }=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}$ resides in the plane of the sky (angular scales) and ${k}_{\parallel }\propto \eta$ is the line-of-sight component. Here η is the Fourier dual of frequency, where we use the fact that the observed frequency of the neutral hydrogen spectral line is linearly proportional to distance. Following Morales & Hewitt (2004) and McQuinn et al. (2006), the transformations between Fourier dimensions and cosmological coordinates are

$$\begin{eqnarray}&&{k}_{\perp }=\displaystyle \frac{2\pi | {\boldsymbol{u}}| }{{D}_{M}(z)},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{k}_{\parallel }=\displaystyle \frac{2\pi {H}_{0}{f}_{21}E(z)}{c{(1+z)}^{2}}\eta ,\end{eqnarray} \tag{ 6 }$$

where ${D}_{M}(z),{H}_{0},{f}_{21},z$ are the transverse comoving distance, Hubble constant, rest frequency of the hyperfine transition (∼1420 MHz), and observation redshift, respectively. E(z) is (Hogg 2000)

$$\begin{eqnarray}&&E(z)\overset{\underset{\mathrm{def}}{}}{=}\sqrt{{{\rm{\Omega }}}_{M}{(1+z)}^{2}+{{\rm{\Omega }}}_{k}{(1+z)}^{2}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}.\end{eqnarray} \tag{ 7 }$$

The high-level design philosophy of CHIPS is to use the data in such a way as to allow for computationally efficient and robust parameter estimation and error covariance estimation, while retaining maximum information. As such, the following design features have been employed:

• 1.  
  Calibrated visibilities will be the primary data input—visibility space is the natural measurement space of interferometric data, where radiometric (stochastic) noise is uncorrelated and Gaussian distributed with color dependent on the frequency-dependent system temperature. In addition, correlations due to point-spread function sidelobes inherent in image space are naturally accounted for in visibility space, where the measured information is clearly defined.
• 2.  
  A least-squares spectral analysis (LSSA) method will be used to compute the line-of-sight transform from frequency to spectral space, not the Fourier transform (Kay 1993)—the MWA bandpass has regular missing channels, owing to the coarse bandpass edges, and intermittent missing channels owing to flagged radio frequency interference (RFI; Offringa et al. 2015). The Fourier transform cannot be used for irregular and missing data. LSSA can compute the optimal spectral representation of the data using all of the available information.
• 3.  
  An ML estimate of the cosmological signal (a quadratic estimator; Kay 1993; Liu & Tegmark 2011; Dillon et al. 2013) will be computed (an inverse-covariance weighting)—the ML solution asymptotes to the optimal solution for large quantities of data and is appropriate for the data set acquired here. It retains the full information contained within the data and allows for computation of the full uncertainty covariance matrix. (Note that this is not the only estimator that can be used efficiently to approach this problem.)
• 4.  
  Foregrounds will be modeled a priori and included in the estimator ("foreground suppression")—as discussed in Section 5.1, foregrounds can be approached in different ways, depending on the degree of knowledge one assumes about them (Liu & Tegmark 2011; Dillon et al. 2013; Chapman et al. 2014; Bonaldi & Brown 2015). Here we take a two-tier approach, where known foregrounds (sky model of point and extended sources) are subtracted from the data coherently, and the remaining signal is treated statistically. CHIPS therefore uses the full Fourier space for estimation (contrast with "foreground avoidance" techniques, where the contaminated regions are excised).

We aim to take the large number of measured, calibrated visibilities and grid them onto a regular grid representing the uv-plane, in each frequency channel. The true sky signal is shaped by the antenna primary beam, and the measured visibilities sample a range of modes. Each datum represents information from a region of the uv-plane, according to the beam, and we grid a measurement across this region. For an instrument with linear crossed-polarization feeds, xx and yy, we model the correlator output as

$$\begin{eqnarray}&&\left[\begin{array}{c}{\tilde{V}}_{{xx}}\\ {\tilde{V}}_{{xy}}\\ {\tilde{V}}_{{yx}}\\ {\tilde{V}}_{{yy}}\end{array}\right]=\left[\begin{array}{c}{\tilde{B}}_{{xx}}\\ {\tilde{B}}_{{xy}}\\ {\tilde{B}}_{{yx}}\\ {\tilde{B}}_{{yy}}\end{array}\right]\\ \quad \quad \ast \left(\begin{array}{cccc}{g}_{{xx}} & {g}_{{xx}}\mathrm{cos}2\psi & {g}_{{xx}}\mathrm{sin}2\psi & 0\\ {g}_{{xy}}{{\rm{\Delta }}}_{{xy}} & -{g}_{{xy}}\mathrm{sin}2\psi & {g}_{{xy}}\mathrm{cos}2\psi & i\\ {g}_{{yx}}{{\rm{\Delta }}}_{{yx}} & -{g}_{{yx}}\mathrm{sin}2\psi & {g}_{{yx}}\mathrm{cos}2\psi & -i\\ {g}_{{yy}} & -{g}_{{yy}}\mathrm{cos}2\psi & -{g}_{{yy}}\mathrm{sin}2\psi & 0\end{array}\right)\left[\begin{array}{c}\tilde{I}\\ \tilde{Q}\\ \tilde{U}\\ \tilde{V}\end{array}\right],\end{eqnarray} \tag{ 8 }$$

which maps the Stokes visibilities in sky coordinates, $(\tilde{I},\tilde{Q},\tilde{U},\tilde{V})$, to the measurement set in instrumental coordinates (Hamaker et al. 1996). Here we are explicitly describing the relationship in the flat-sky approximation, although the curvature terms are introduced below. The $4\times 4$ matrix encodes projection effects from the parallactic angle, ψ, polarization leakage, $({{\rm{\Delta }}}_{{xy}},{{\rm{\Delta }}}_{{yx}})$, and the direction-independent, antenna-based complex gains, $({g}_{{xx}},{g}_{{xy}},{g}_{{yx}},{g}_{{yy}})$. The convolution accounts for the primary beam shape, $({\tilde{B}}_{{xx}},{\tilde{B}}_{{xy}},{\tilde{B}}_{{yx}},{\tilde{B}}_{{yy}})$. To extract the total intensity, $\tilde{I}$, we could use all four polarization outputs from the correlator. Instead, to reduce computational load, we use only ${V}_{{xx}}$ and ${V}_{{yy}}$. This choice is made at the expense of some signal when the zenith angle is large. Failure to treat individual linear feeds will lead to polarization leakage if the visibilities are combined without consideration for beam shape differences (Moore et al. 2013).

Absorbing the complex gains into the beams, ${\boldsymbol{B}}$, we write the xx and yy visibilities as

$$\begin{eqnarray}{\tilde{V}}_{{xx}}(u,v,w) & = & (\tilde{I}+\tilde{Q}\mathrm{cos}2\psi +\tilde{U}\mathrm{sin}2\psi )\\ & & \ast \;{\tilde{B}}_{{xx}}\ast \;\displaystyle \int \displaystyle \frac{\mathrm{exp}[-2\pi i({\boldsymbol{D}}\cdot {\boldsymbol{s}})]}{\sqrt{1-{l}^{2}-{m}^{2}}}{dldm}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}{\tilde{V}}_{{yy}}(u,v,w) & = & (\tilde{I}-\tilde{Q}\mathrm{cos}2\psi -\tilde{U}\mathrm{sin}2\psi )\\ & & \ast {\tilde{B}}_{{yy}}\ast \displaystyle \int \displaystyle \frac{\mathrm{exp}[-2\pi i({\boldsymbol{D}}\cdot {\boldsymbol{s}})]}{\sqrt{1-{l}^{2}-{m}^{2}}}{dldm},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{D}}\cdot {\boldsymbol{s}}={ul}+{vm}+w(\sqrt{1-{l}^{2}-{m}^{2}}-1))\end{eqnarray} \tag{ 11 }$$

is the projection of the baseline vector onto the sky coordinates, and we have now explicitly included the w-terms in the equations through the phase integral (final term).

Using the convolution theorem, the mapping from the underlying sky Fourier representation and the measured visibilities can be written as

$$\begin{eqnarray}&&{\tilde{V}}_{{xx}}={\tilde{G}}_{w}\ast {\tilde{B}}_{{xx}}\ast (\tilde{I}+\tilde{Q}\mathrm{cos}2\psi +\tilde{U}\mathrm{sin}2\psi )\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\tilde{V}}_{{yy}}={\tilde{G}}_{w}\ast {\tilde{B}}_{{yy}}\ast (\tilde{I}-\tilde{Q}\mathrm{cos}2\psi -\tilde{U}\mathrm{sin}2\psi ),\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}{\tilde{G}}_{w} & \overset{\underset{\mathrm{def}}{}}{=} & \displaystyle \int \displaystyle \frac{\mathrm{exp}[-2\pi i(w(\sqrt{1-{l}^{2}-{m}^{2}}-1))]}{\sqrt{1-{l}^{2}-{m}^{2}}}\\ & & \times \mathrm{exp}[-2\pi i({ul}+{vm})]{dldm}\end{eqnarray} \tag{ 14 }$$

explicitly highlights the additional convolution due to curvature terms. These functions encode the deviation from a strict 2D Fourier transform between sky and visibilities and the spectral leakage due to the primary beam.

We now expand from considering a single visibility to a set and describe these as a vector. In doing so, we extend the convolutions to describe the transform from the underlying sky, $\tilde{I}$, to the full measurement set of data. We can write the discrete convolutions as matrix operations, to encode this transformation, and find the measurements for a set of baselines:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{V}}}}_{{xx}}={{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{xx}}(\tilde{{\boldsymbol{I}}}+\tilde{{\boldsymbol{Q}}}\mathrm{cos}2\psi +\tilde{{\boldsymbol{U}}}\mathrm{sin}2\psi )\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{V}}}}_{{yy}}={{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{yy}}(\tilde{{\boldsymbol{I}}}-\tilde{{\boldsymbol{Q}}}\mathrm{cos}2\psi -\tilde{{\boldsymbol{U}}}\mathrm{sin}2\psi ).\end{eqnarray} \tag{ 16 }$$

We are interested in determining the underlying Stokes I sky distribution at a given angular scale, ${u}_{\alpha }$, given the measured visibilities, $\tilde{{\boldsymbol{V}}}(u)$. A vector of ($N\times 1$) measured visibilities, $\tilde{{\boldsymbol{V}}}$, is generated from an underlying ($M\times 1$) sky distribution, $\tilde{{\boldsymbol{I}}}$, via an ($N\times M$) matrix transform, ${G}_{w}B$. We can rewrite and combine Equations (15) and (16) such that

$$\begin{eqnarray}\tilde{{\boldsymbol{I}}}(u,v,\nu ) & = & \displaystyle \frac{1}{2}{({{\boldsymbol{B}}}_{{xx}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{xx}})}^{-1}{{\boldsymbol{B}}}_{{xx}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{\tilde{{\boldsymbol{V}}}}_{{xx}}\\ & & +\displaystyle \frac{1}{2}{({{\boldsymbol{B}}}_{{yy}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{yy}})}^{-1}{{\boldsymbol{B}}}_{{yy}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{\tilde{{\boldsymbol{V}}}}_{{yy}},\end{eqnarray} \tag{ 17 }$$

where the square matrix inversion involves inverting an $(M\times M)$ Hermitian complex-valued matrix, which will be (almost) diagonal for independent modes, ${\boldsymbol{u}}$. Equation (17) effectively unwraps the effects of w-terms and the polarized beams, before combining xx and yy information to remove the polarized component of the sky signal ($\tilde{Q},\tilde{U}$). The ${{\boldsymbol{G}}}^{\dagger }$ matrix serves to distribute the footprint of the visibilities across the uv-plane, as dictated by the amplitude of the w-term; for large w, the Fourier-space beam is effectively broader and the G matrix captures this. In practice, as discussed by Dillon et al. (2013), the uv-plane needs to be fully sampled (the matrix full rank) for this matrix inverse to exist (i.e., there needs to be independent information in all modes). For the MWA and using rotation synthesis, this requirement is met for parts of the uv-plane, but not in general.

There are multiple approaches to handling this, including computation of the pseudoinverse and simplification of the matrix by considering only diagonal components. The former inverts the matrix in regions where there is sufficient information. The latter effectively ignores the correlations introduced by the finite extent of the beam across the uv-plane and reduces the matrix inverse step to a simple weighting of the data according to the measurements (averaging). It does, however, still encode the distribution of power across multiple modes in the uv-plane. The approximation is reasonable given that the matrix is already highly diagonal. In the testing phases of CHIPS, as described here, we will implement the latter:

$$\begin{eqnarray}\tilde{{\boldsymbol{S}}}({\boldsymbol{u}},{\boldsymbol{v}},\nu ) & = & \displaystyle \frac{1}{2}{({{\boldsymbol{B}}}_{{xx}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{xx}})}_{\mathrm{diag}}^{-1}{{\boldsymbol{B}}}_{{xx}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{\tilde{{\boldsymbol{V}}}}_{{xx}}\\ & & +\displaystyle \frac{1}{2}{({{\boldsymbol{B}}}_{{yy}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{{\boldsymbol{G}}}_{w}{{\boldsymbol{B}}}_{{yy}})}_{\mathrm{diag}}^{-1}{{\boldsymbol{B}}}_{{yy}}^{\dagger }{{\boldsymbol{G}}}_{w}^{\dagger }{\tilde{{\boldsymbol{V}}}}_{{yy}}.\end{eqnarray} \tag{ 18 }$$

This introduces an error, for each polarization, with the XX polarization expression:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{S}}}}_{{XX}}({\boldsymbol{u}},\nu )={{\boldsymbol{R}}}_{{XX}}{\tilde{{\boldsymbol{I}}}}_{{XX}}({\boldsymbol{u}},\nu ).\end{eqnarray} \tag{ 19 }$$

Here the error matrix, ${{\boldsymbol{R}}}_{{XX}}$, is

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{{XX}}\overset{\underset{\mathrm{def}}{}}{=}{({\rm{diag}}({{\boldsymbol{B}}}^{\dagger }{\boldsymbol{B}}))}_{{XX}}^{-1}{({{\boldsymbol{B}}}^{\dagger }{\boldsymbol{B}})}_{{XX}}.\end{eqnarray} \tag{ 20 }$$

The same applies to the YY polarization, where the error matrix, in general, is different, depending on the shape of the beams for each telescope pointing. In general, this can invalidate the canceling of the Q and U linear polarization signals in Equation (17) because there is now the potential for mismatch between the two polarizations. In practice, the XX and YY polarizations are not combined in CHIPS at this point (the results described in this work remain in the two instrumental poarizations) because the current MWA beam model is known to be imperfect. In the future, this mismatch will need to be further studied for leakage effects.

Upon gridding the data onto the uv-plane for each frequency channel, the next step is to estimate the line-of-sight wavenumber information, by transforming the data at each grid point from frequency, ν, to η. For a complete and regularly sampled set of N complex-valued discrete data points embedded within white Gaussian noise, $\tilde{{\boldsymbol{S}}}(\nu )$, the information contained within a spectral mode, η, can be obtained using the DFT:

$$\begin{eqnarray}&&{ \mathcal S }(u,v,\eta )=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \sum _{j=0}^{N-1}\tilde{S}(u,v,{\nu }_{j})\mathrm{exp}[-2\pi {ij}\eta /N],\end{eqnarray} \tag{ 21 }$$

where $\eta \in [0,N-1]$. This can be written in matrix formalism, where the Fourier transform has convenient properties as a Vandermonde matrix:

$$\begin{eqnarray}&&\tilde{{\boldsymbol{S}}}(\eta )={{ \mathcal F }}^{\dagger }\tilde{{\boldsymbol{S}}}(\nu ),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{ \mathcal F }}^{\dagger }{ \mathcal F }={ \mathcal I }.\end{eqnarray} \tag{ 23 }$$

Here $\tilde{{\boldsymbol{S}}}(\eta )$ is the complex-valued estimate of the spectral information in mode η, and ${ \mathcal F }$ is a square, Hermitian matrix containing the trigonometric kernel. We have also explicitly dropped the parameterization by $u,v$ because the transform is performed for each point on the uv-plane, with the understanding that these parameters are carried along for the computations presented in this section.

When the data have generalized covariance (correlated samples, unequal weightings), the signal in some spectral mode can be estimated using a generalized ML Fourier transform. If the data are distributed such that

$$\begin{eqnarray}&&\tilde{{\boldsymbol{S}}}(\nu )\sim { \mathcal C }{ \mathcal N }(\bar{S}(\nu ),{\boldsymbol{C}}),\end{eqnarray} \tag{ 24 }$$

the optimal estimate is given by (i.e., the inverse-covariance weighting estimator; Liu & Tegmark 2011; Dillon et al. 2015)

$$\begin{eqnarray}&&\tilde{{\boldsymbol{S}}}(\eta )\sim { \mathcal C }{ \mathcal N }({{ \mathcal F }}^{\dagger }{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}(\nu ),{{ \mathcal F }}^{\dagger }{\boldsymbol{C}}{ \mathcal F }),\end{eqnarray} \tag{ 25 }$$

where we have used the fact that the Fourier transform is unitary (${{ \mathcal F }}^{\dagger }={{ \mathcal F }}^{-1}$). Effectively, this prewhitens the signal (by suppressing data with a lot of noise and removing correlation) and computes the Fourier transform. The result has noise properties consistent with a prewhitening operation and with the reordering of the data according to the summing of phased data points.

Finally, when the data have incomplete sampling, the above formalism naturally generalizes and is described by the LSSA. When there are only M sampled points at locations ${\nu }_{j}$, within N frequency channels, the generalized least-squares estimate is

$$\begin{eqnarray}&&\tilde{{\boldsymbol{S}}}(\eta )={ \mathcal H }\tilde{{\boldsymbol{S}}}(\nu )\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\overset{\underset{\mathrm{def}}{}}{=}\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \sum _{j=0}^{M-1}\tilde{S}({\nu }_{j})\mathrm{exp}[-2\pi i{\nu }_{j}\eta /M],\end{eqnarray} \tag{ 27 }$$

where $\eta =[0,\ldots ,\;{N}_{\mathrm{freq}}]$ are evenly spaced and ${N}_{\mathrm{freq}}\lt N$. The optimal estimator is

$$\begin{eqnarray}&&\tilde{{\boldsymbol{S}}}(\eta )\sim { \mathcal C }{ \mathcal N }({({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })}^{-1}{{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}(\nu ),{({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })}^{-1}),\end{eqnarray} \tag{ 28 }$$

where the LSSA operator has been labeled ${ \mathcal H }$. Note that now $M\lt N$ and the ${ \mathcal H }$ matrix is not square. In general, this method leads to correlated information between different modes, but this can be reduced by estimating ${N}_{\mathrm{freq}}\lesssim M/2$ modes (Vio et al. 2013).

Now that we have described the computation of the coherently averaged (gridded and LSSA) data from the underlying measured visibilities, we can write the likelihood function of the data to compute the power estimates.

We describe the data with a generalized multivariate Gaussian, with zero mean, and general covariance matrix, ${\boldsymbol{C}}$. This form describes data that are drawn from a statistical distribution of possible universes, with temperature fluctuations that are Gaussian distributed. This covariance matrix contains all of the terms contributing to the signal in a visibility. The joint likelihood for the vector of complex-valued coherently averaged data, ${\boldsymbol{S}}(u,v,\eta )$, is given by

$$\begin{eqnarray}&&L(\tilde{{\boldsymbol{S}}};{\boldsymbol{C}})=\displaystyle \frac{1}{{\pi }^{N}{\rm{det}}({\boldsymbol{C}})}\mathrm{exp}[-{\tilde{{\boldsymbol{S}}}}^{\dagger }{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}],\end{eqnarray} \tag{ 29 }$$

where the complex-valued covariance matrix is

$$\begin{eqnarray}&&{\boldsymbol{C}}\overset{\underset{\mathrm{def}}{}}{=}\langle \tilde{{\boldsymbol{S}}}{\tilde{{\boldsymbol{S}}}}^{\dagger }\rangle \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&=\;{{\boldsymbol{C}}}_{\mathrm{FG}}+{{\boldsymbol{C}}}_{{\rm{N}}}+{{\boldsymbol{C}}}_{{\rm{P}}},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}{{\boldsymbol{C}}}_{{\rm{P}}}=\left(\begin{array}{cccc}{p}_{11} & {p}_{12} & \cdots & {p}_{1N}\\ {p}_{12} & {p}_{22} & \cdots & {p}_{2N}\\ \vdots & \vdots & \ddots & \vdots \\ {p}_{1N} & {p}_{2N} & \cdots & {p}_{{NN}}\end{array}\right)\end{eqnarray} \tag{ 32 }$$

describes the parameters of interest (the mode powers, ${p}_{\alpha \alpha }$) and the correlations between powers imprinted by the instrument and the experiment ("window functions"). Strictly, the cosmological information forms a diagonal matrix, with uncorrelated power between modes, but the imperfect data correlate contiguous modes. In this implementation of CHIPS, we assume a diagonal cosmological signal matrix, but do encode the correlations between modes in the data by capturing and propagating the covariances introduced by the line-of-sight transform (LSSA). These correlations are important when combining modes together to perform the binning to 1D (see Section 4.3.1). Note that here we are explicitly identifying the power in which we are interested with the likelihood function of the coherently averaged data. The other terms in Equation (31) describe the statistical contribution to mode ($u,v,\eta$) from foregrounds (${{\boldsymbol{C}}}_{\mathrm{FG}}$) and measurement (radiometric) noise (${{\boldsymbol{C}}}_{{\rm{N}}}$). In uv space, measurement noise is uncorrelated and Gaussian distributed (colored Gaussian noise [CGN]), where the "color" for a coherently averaged datum describes the weighting of the noise due to the number of visibilities contributing to that cell. The statistical foreground contribution will be treated in Section 5.1.

The ML estimate, obtained by finding the value of a given parameter that maximizes the likelihood function (or minimizes the negative log likelihood), is asymptotically efficient (achieves optimal estimation precision for large data sets). We aim to determine an expression for the parameters of interest (mode powers) in terms of the data and data covariance matrix. The log likelihood function is minimized by differentiating with respect to the parameter, setting to zero, and solving for that parameter. The derivative of the log likelihood for a zero-mean generalized Gaussian is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{ln}}L}{\partial {p}_{\alpha }}={\rm{tr}}\left({{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)-{\tilde{{\boldsymbol{S}}}}^{\dagger }(u,v,\eta ){{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}(u,v,\eta ),\end{eqnarray} \tag{ 33 }$$

where "tr" denotes the trace of the matrix. For clarity, we now drop the dependence of the data on ($u,v,\eta$), with the understanding that we are working with the data in wavenumber space. Setting Equation (33) to zero, and then using in the first term the identity ${\mathbb{I}}={{\boldsymbol{C}}}^{-1}{\boldsymbol{C}}$ and expanding the covariance matrix into its constituent parts, we find

$$\begin{eqnarray}&&{\hat{p}}_{\alpha }{\rm{tr}}\left({{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)\overset{\underset{\mathrm{def}}{}}{=}{\rm{tr}}\left({{\boldsymbol{C}}}_{{\rm{P}}}{{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&=\;{\tilde{{\boldsymbol{S}}}}^{\dagger }{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}-{\rm{tr}}\left(({{\boldsymbol{C}}}_{\mathrm{FG}}+{{\boldsymbol{C}}}_{{\rm{N}}}){{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right),\end{eqnarray} \tag{ 35 }$$

where ${\hat{p}}_{\alpha }={\hat{p}}_{\alpha \alpha }$ denotes our estimate of the power in mode α, as described in the diagonal elements of the cosmological power matrix, Equation (32). Note here that the first term on the right-hand side of Equation (35) is squaring the coherently averaged data and therefore producing the desired power-like quantity.

Therefore, the estimate of the power in mode α is given by computation of

$$\begin{eqnarray}{\hat{p}}_{\alpha } & = & \displaystyle \frac{1}{{\rm{tr}}\left({{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)}\left({\tilde{{\boldsymbol{S}}}}^{\dagger }{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}{{\boldsymbol{C}}}^{-1}\tilde{{\boldsymbol{S}}}\right.\\ & & -\left.{\rm{tr}}\left(({{\boldsymbol{C}}}_{\mathrm{FG}}+{{\boldsymbol{C}}}_{{\rm{N}}}){{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)\right),\end{eqnarray} \tag{ 36 }$$

where we decouple dependence of the estimator on the mode powers by making the approximation

$$\begin{eqnarray}&&{\boldsymbol{C}}\approx {{\boldsymbol{C}}}_{\mathrm{FG}}+{{\boldsymbol{C}}}_{{\rm{N}}},\end{eqnarray} \tag{ 37 }$$

under the assumption that ${{\boldsymbol{C}}}_{\mathrm{FG}}+{{\boldsymbol{C}}}_{{\rm{N}}}\gg {{\boldsymbol{C}}}_{{\rm{P}}}$.

In practice, the data can be split into two and cross-correlated. This removes thermal noise power from the final power estimate, at the expense of a factor of $\sqrt{2}$ in the final signal-to-noise ratio (S/N) (in temperature; a factor of 2 in power). Although noise power has been removed from the estimate, thermal noise uncertainty remains in the final covariance of the power estimates. Typically, alternate correlator output data sets are processed into alternate coherent data vectors, to ensure uniform uv coverage between the two (this is important for applying the same data covariance matrix to both data sets). Using this scheme, the cross-correlation power formed from coherent data vectors, ${}_{1}{\boldsymbol{S}}$ and ${}_{2}{\boldsymbol{S}}$, is

$$\begin{eqnarray}{\hat{p}}_{\alpha } & = & \displaystyle \frac{1}{{\rm{tr}}\left({{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)}\left({{}_{1}{\boldsymbol{S}}}^{\dagger }{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}{{\boldsymbol{C}}}_{}^{-1}{}_{2}{{\boldsymbol{S}}}_{\mathrm{}}\right.\\ & & -\left.{\rm{tr}}\left({{\boldsymbol{C}}}_{\mathrm{FG}}{{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}^{-1}\displaystyle \frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}\right)\right).\end{eqnarray} \tag{ 38 }$$

The expected values and covariances are derived in the Appendix and demonstrate that the expected value of the ML estimate yields the cosmological power.

The second term in Equation (38) is the foreground power bias. As discussed in a number of recent papers (Liu & Tegmark 2011; Dillon et al. 2013, 2015), subtraction of this term requires an accurate model for the foregrounds. Inaccurate subtraction leads to bias in the power. Given the relative power of the foreground and cosmological signal, mis-subtraction will occur when the foreground power is incorrect at the fraction of a percent level. To avoid this, we retain the foreground power bias and allow the noise covariance to indicate which modes are contaminated (i.e., the second term is set to zero).

4.3.1. Explicitly Including LSSA

To connect the LSSA transform with the ML estimate, we can take the expression for the optimal power, Equation (38), which is a function of ($u,v,\eta$), and insert the LSSA estimate of $\tilde{{\boldsymbol{S}}}(u,v,\eta )$ from the underlying gridded data, $\tilde{{\boldsymbol{S}}}(u,v,\nu )$. To do so, we combine Equations (28) and (38) to compute the optimal (ML) power estimate:

$$\begin{eqnarray}&&{\hat{p}}_{\alpha }=\displaystyle \frac{{({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}_{}^{-1}{}_{1}\bar{S}(\nu ))}^{\dagger }\frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}_{}^{-1}{}_{2}\bar{S}(\nu ))}{{\rm{tr}}\left(({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })\frac{\partial {\boldsymbol{C}}}{\partial {p}_{\alpha }}({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })\right)}.\end{eqnarray} \tag{ 39 }$$

The denominator (the normalization) reduces to

$$\begin{eqnarray}&&\displaystyle \sum _{i}| {A}_{i\alpha }{| }^{2},\end{eqnarray} \tag{ 40 }$$

where A_ij are the elements of the inverse-covariance matrix,

$$\begin{eqnarray}&&{\boldsymbol{A}}\overset{\underset{\mathrm{def}}{}}{=}{{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H }.\end{eqnarray} \tag{ 41 }$$

This then yields the sum over the square of the matrix elements for the normalization. Finally, the vector of power estimates is ${\chi }^{2}$-distributed:

$$\begin{eqnarray}&&\hat{{\boldsymbol{p}}}(\eta )\sim {\chi }^{2}({\boldsymbol{p}},2{(({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H }))}^{-1}),\end{eqnarray} \tag{ 42 }$$

where the (identity) diagonal cosmological signal matrix, $d{\boldsymbol{C}}/{dp}$, has been dropped.

To compute the spherically-averaged (1D) power spectrum from the 2D output, we average in elongated cylindrical shells, allowing for differences in the line of sight and angular modes for the cosmological signal in redshift space in the conversion from $u,v,\eta$ to k (Mpc⁻¹). The ML estimate of the 1D power is given by weighting the individual 2D mode powers contributing to a cell, k, with the inverse of their power covariance matrix. This is the inverse of the covariance term from Equation (42), with

$$\begin{eqnarray}&&{ \mathcal A }\overset{\underset{\mathrm{def}}{}}{=}({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H })({{ \mathcal H }}^{\dagger }{{\boldsymbol{C}}}^{-1}{ \mathcal H }).\end{eqnarray} \tag{ 43 }$$

Thus,

$$\begin{eqnarray}&&{P}_{k}=\displaystyle \frac{{\displaystyle \sum }_{i\in k}{{ \mathcal A }}_{,i}{\boldsymbol{p}}}{{\rm{tr}}({D}^{\dagger }{{ \mathcal A }}_{,i}D)},\end{eqnarray} \tag{ 44 }$$

where the covariance matrix of data is used quadratically to reflect that we are now working to optimally average power, and D is a binning matrix that reduces the full 2D space to the subspace spanned by $i\in k$. The numerator terms here come from performing the matrix multiplication from Equation (42) and including all of the covariances between 2D modes, i, contributing to binned mode k. The denominator provides the variances and covariances between the 1D k modes, encoding the degree of uncertainty on each parameter and the relationship between estimates. The syntax, ${}_{,i}$, denotes the matrix, ${ \mathcal A }$, as formed from projecting onto the bin, including parameter covariances.

Dealing with foreground contamination can be approached in different ways, all of which rely on the basic distinction that foreground sources are continuum emitters and therefore have a smooth spectrum over a small bandwidth, and the cosmological signal is a spectral line with structure reflecting the brightness temperature at different spatial depths. This distinction fundamentally dissociates angular modes from line-of-sight modes and is the reason for operating with the 2D power spectrum. In 2D, the power from flat-spectrum foregrounds that could be observed through a perfect instrument (no frequency-dependent instrumental response, no spectral incompleteness, and complete sampling of the Fourier plane) would be contained within the DC ($\eta =0$) mode of ${k}_{\parallel }$ because the flat-spectrum foregrounds add a simple amplitude across frequency. For a real instrument and nonflat sources, foregrounds occupy a broader region of parameter space, colloquially termed the "wedge" because of its broad wedge-like structure in (${k}_{\perp },{k}_{\parallel }$) space. The wedge has been well studied (Datta et al. 2010; Morales et al. 2012; Parsons et al. 2012; Trott et al. 2012; Vedantham et al. 2012; Hazelton et al. 2013; Thyagarajan et al. 2013). It results from the incomplete sampling of the Fourier modes and/or spectral modes by an interferometer and the migration of the angular mode sampled by a given baseline as a function of observation frequency ("mode mixing"). A mathematical derivation of the expected structure of the wedge for the MWA is included in this section.

Broadly, there are two primary foreground treatment design options employed in the literature: (1) "avoidance," where it is attempted to contain foregrounds to a region of parameter space and ignore this region (e.g., PAPER's approach; Parsons et al. 2012; Jacobs et al. 2015); and (2) "suppression," where the foreground contribution to the signal is modeled (using an a priori source model, nonparametric methods, or a data-driven model) and the estimator uses this knowledge to suppress contaminated information. The latter option includes the nonparametric fitting of low-order polynomials (Bowman et al. 2009) and the more sophisticated Component Analysis techniques discussed by Chapman et al. (2012, 2014), Bonaldi & Brown (2015), and references therein. It also includes model- and data-driven (parametric) statistical descriptions of the sky itself and incorporation of that information (Liu & Tegmark 2011; Dillon et al. 2013, 2015). These approaches all have their own advantages and disadvantages. Broadly, nonparametric approaches are simple to implement but, requiring no input physical knowledge, can destroy cosmological information and retain foreground signal. Typically, the metric for an adequate solution is not well defined. Conversely, parametric models are designed to include as much physical information about the foregrounds and cosmological signal as is available. However, they can be computationally difficult to implement, are limited by our knowledge of the sky at low radio frequencies, and can also destroy information if the models are incorrect. Surveys such as the MWA's GLEAM (Wayth et al. 2015) and LOFAR's MSSS (Heald et al. 2015) will be crucial for improving our understanding.

Ultimately, the key to a robust estimator is to have an understanding of the limitations and biases of one's method and incorporate that knowledge into the methodology and results. CHIPS aims to do this and handles foregrounds using a two-tiered approach: (1) known foregrounds (sky model of point and extended sources, with updated calibration solutions and ionospheric corrections) are subtracted from the data coherently (visibility data) using the Real-Time System (RTS) calibration (see Section 7); (2) remaining signal is handled statistically, using an a priori model for the expected distribution and spectral structure of sources (CHIPS). We then use a perturbative calculation to ascertain our degree of confidence in the model (Liu et al. 2015).

We initially employ a two-component statistical foreground model, consisting of extragalactic point sources and Galactic synchrotron, described here. The component describing the Galactic plane is omitted. The Galactic plane is more difficult to describe statistically, with definite sky position and skewed statistics. We leave this as an open problem for future work.

Point-source Covariance Matrix

Galactic Synchrotron Covariance Matrix

Galactic synchrotron emission, from electrons spiraling along our Galaxy's magnetic field lines, produces emission on large scales. We use the parameters of Jelić et al. (2008, 2010) to motivate our model, which also includes the effects of the MWA primary beam and instrumental chromaticity.

The intrinsic temperature power spectrum is modeled as

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{T}_{\mathrm{GS}}^{2}\rangle (u,\nu )={(\eta {T}_{B})}^{2}{\left(\displaystyle \frac{u}{{u}_{0}}\right)}^{-2.7}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-2.55}\;{{\rm{K}}}^{2},\end{eqnarray} \tag{ 55 }$$

where $\eta =0.01$ is the fluctuation level relative to the uniform brightness temperature, ${T}_{B}=253\;{\rm{K}}$, ${u}_{0}=10$ wavelengths, and ${\nu }_{0}=100\;{\rm{MHz}}$ is the reference frequency. The covariance matrix is a function of the parameter u, to reflect that we expect the statistics to be rotationally invariant. The apparent steep spectral index in temperature units is flattened once converting to flux density (integrated) units for a given instrument, such that the intrinsic power is

$$\begin{eqnarray}&&{P}_{\mathrm{GS}}(u,\nu )={\left(\displaystyle \frac{2k}{{\lambda }^{2}}\right)}^{2}{\rm{\Omega }}{(\eta {T}_{B})}^{2}{\left(\displaystyle \frac{u}{{u}_{0}}\right)}^{-2.7}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-2.55}\;{\mathrm{Jy}}^{2}\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&=\;{\left(\displaystyle \frac{2k}{\lambda }\right)}^{2}\displaystyle \frac{1}{{A}_{\mathrm{eff}}}{(\eta {T}_{B})}^{2}{\left(\displaystyle \frac{u}{{u}_{0}}\right)}^{-2.7}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-2.55}\;{\mathrm{Jy}}^{2}\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&=\;\left(\displaystyle \frac{{(2k)}^{2}}{{A}_{\mathrm{eff}}}\right){(\eta {T}_{B})}^{2}{\left(\displaystyle \frac{u}{{u}_{0}}\right)}^{-2.7}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{-0.55}\;{\mathrm{Jy}}^{2}.\end{eqnarray} \tag{ 58 }$$

The truncation of the spectrum due to the primary beam can be obtained via convolution (via a Fourier transform), and Figure 1(a) shows the intrinsic and convolved power spectra. Here the Fourier transform is used (rather than LSSA) because the plot uses complete and regularly sampled spacing for display purposes. The actual covariance contains the same incompleteness as the data, when implemented in CHIPS. Note that the primary beam convolution is frequency dependent, generating complex structure at small wavenumbers.

Zoom In Zoom Out Reset image size

Finally, the instrumental chromaticity is included in a similar manner as for the point-source model, yielding the following model for the spectral covariance matrix:

$$\begin{eqnarray}{{\boldsymbol{C}}}_{\mathrm{GS}}(\nu ,\nu ^{\prime} ;u) & = & 2\pi \sqrt{{P}_{\mathrm{GS}}^{\prime }(u,\nu ){P}_{\mathrm{GS}}^{\prime }(u,\nu ^{\prime} )}\\ & & \times \displaystyle \int {J}_{0}(2\pi ({ul})f(\nu )){ldl}\ {\mathrm{Jy}}^{2},\end{eqnarray} \tag{ 59 }$$

where $P^{\prime}$ denotes the beam-convolved power.

5.1.1. Total Covariance Matrix

The final component of the data covariance corresponds to the thermal (measurement, radiometric) noise. This is the measurement uncertainty on each visibility due to the finite number of data samples (information) that contribute to it. We are effectively estimating the sky signal with a fixed amount of information, which depends on the signal strength and the number of samples. The former is set by the system temperature, which is dominated by sky temperature at low frequencies, and the latter is set by the bandwidth and sampling time for each visibility. For a single polarization,

$$\begin{eqnarray}&&\sigma ={10}^{26}\displaystyle \frac{2{{kT}}_{\mathrm{sys}}}{{A}_{\mathrm{eff}}}\displaystyle \frac{1}{\sqrt{{\rm{\Delta }}\nu {\rm{\Delta }}t}}\;{\rm{Jy}}.\end{eqnarray} \tag{ 60 }$$

When gridding the visibilities onto the uv-plane, as described in Section 4.1, the noise decreases coherently (with square-root improvement). This noise reduction therefore follows the evolution of sampled points in the uv-plane for a nightly track of the EoR field. Identical observations on subsequent nights yield the same coherent reduction in power. The thermal noise contribution to the power spectrum therefore maps the distribution of points in the uv-plane for a nightly track, thereby reflecting the array configuration.

The set of high-band EoR data were taken from a single night of observations during the first semester of MWA EoR observing (2013 August 23). These data were chosen to test, refine, and compare different calibration and analysis methodologies, as described in Jacobs et al. (2015). The original data set consists of 160 112 s observations. These data were refined to a final set of 94 observations, removing pointings that were $\gt 25$° from zenith and a pointing heavily affected by Galactic emission in the sidelobes. As described in Section 7, the data were calibrated and provided as 8 s visibilities, yielding 14 time steps per observation and seven per interleaved data set per observation (temporally interleaved data are used to compute the cross-power spectrum). After calibration, we performed a uv cut at a maximum distance of 300 wavelengths at the lowest frequency, within which the EoR signal is expected to fall, and used the full bandwidth data set for the initial analysis (30.72 MHz).

The cross-power spectra are produced with CHIPS both with and without the foreground covariance matrix included in the estimator. Without foreground covariance, the data are weighted purely by the system temperature for that observation, the relative weight of the visibility (determined from the number of 10 kHz channels contributing), and the uv-sampling of the instrument and correspond to the most straightforward application of a power spectrum analysis. With foreground covariance, the data are downweighted according to the modeled and measured noise and signal contamination within the data.

Figure 4 shows the instrumental north–south and east–west cross-power polarizations for the full bandwidth. The solid and dashed diagonal lines correspond roughly to the expected locations of foregrounds within the main lobe of the MWA primary beam and the horizon, respectively. There is clear foreground contribution within the expected region. Note that this wedge emission is consistent with the input foreground model, as shown in Figure 1. There is also the clear imprint of the MWA's uv-sampling function, which has few very short and few very long ($k\gt 0.08$ h Mpc⁻¹) baselines. The copies of the wedge are visible at regular intervals in ${k}_{\parallel }$, corresponding to the comb-like bandpass sampling function of the MWA. The errors and expected noise and the measured signal-to-error ratio can also be computed (Figures 5 and 6). The error is computed using Equation (42) for the thermal-noise-only data covariance matrix, as considered for the estimator (foreground data covariance contribution is omitted here). The expected thermal noise is computed by gridding visibilities with a thermal noise contribution given an average system temperature of ${T}_{\mathrm{sys}}=240$ K and propagating through the same estimator. This system temperature was estimated by matching the observed noise (obtained through the difference in the even and odd data sets) to unity-valued gridded visibilities. Note that this is only the thermal noise contribution, and the signal-to-error plot shows all of the high ratio detections expected in the foreground-dominated regions. The signal-to-error plot shows behavior consistent with thermal noise away from the contaminated wedge and bandpass harmonics regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the foreground contribution is included in the data covariance matrix, power is effectively removed from the contaminated regions, and the errors reflect those parts of the parameter space (Figure 7(b)). The errors at the locations of the coarse bandpass features are also elevated, but this is not obvious from the plots as shown. The small signal-to-error ratio in the wedge indicates that these modes are highly contaminated and should be downweighted in the final binning from 2D to 1D. The ratio is close to unity in the wedge, suggesting that the foreground covariance is capturing the contamination appropriately. These data can then be averaged in cylindrical bins and normalized to obtain the spherically-averaged dimensionless power spectrum. Instead of treating the full bandwidth, which corresponds to a redshift range of 6.2–7.5 and therefore is highly likely to contain signal evolution, we follow Dillon et al. (2015) and split the coherent data into three contiguous 10.24 MHz bins, corresponding to redshifts z = [6.2–6.6], [6.6–7.0], and [7.0–7.5], and compute the 1D power spectra for those (Figure 8) for the E-W polarization (which shows reduced foreground leakage from the Galaxy). This is effectively using the top and bottom plots from Figure 7 as the power and errors, but with the reduced bandwidth (it is not exactly this process, because the inverse-covariance estimator also uses all of the off-diagonal error terms, which are not represented in the 2D error plots). Unlike Dillon et al. (2015), however, we do not exclude any wedge contribution here and instead allow the estimator to "work within the wedge" and downweight contaminated modes. The only data cut performed is to remove the ${k}_{\perp }=0$ and ${k}_{\parallel }=0$ bins, which show large contamination and a poor foreground model response. Both the thermal-noise-only and the full thermal noise+foreground data covariance matrices are shown. The error bars reflect 95% confidence regions in both dimensions. Inclusion of the foreground model increases the uncertainties on those modes, corresponding to the substantial contamination in those regions. The full data covariance model increases the uncertainties such that the results are consistent with thermal noise+foregrounds across most of parameter space. The "detections" in the ${k}_{\perp }\quad =\quad$[0.03–0.10] h Mpc⁻¹ region are due to power bias from unmodeled foregrounds, recalling that the power bias term from Equation (36) is omitted from the estimator.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From these results, we can set 2σ upper limits on the EoR power spectrum at the point where the 1D power is at a minimum and we achieve a "detection," corresponding to k = 0.05 h Mpc⁻¹. (We omit the inner detection because the power associated with this bin includes power on the scale of the primary beam size). Table 1 lists the upper limits at this wavenumber for the 3 hr of data in the E-W polarization. While not competitive with the results from deeper studies with other telescopes, it is consistent with previously published values, with a best detection value of ${{\rm{\Delta }}}_{k}^{2}$ = (275 mK)² in the redshift range z = [6.2–6.6].

Table 1.  95% Confidence (2σ) Upper Limits on the EoR Power in Three Redshift Bins at $k=0.05\;{\rm{h}}$ Mpc⁻¹, Using the Full CHIPS Estimator and Including All Data

Redshift	${{\rm{\Delta }}}_{k}^{2}$ (mK2)
z = [6.2–6.6]	7.6 × 104
z = [6.6–7.0]	8.8 × 104
z = [7.0–7.5]	16.5 × 104

Download table as:  ASCIITypeset image