BRIGHT SOURCE SUBTRACTION REQUIREMENTS FOR REDSHIFTED 21 cm MEASUREMENTS

Our main aim is to explore the level of accuracy needed in instrument calibration and foreground modeling in order to ensure that the residual errors from bright foreground source removal do not obscure the detection of the signal from cosmic reionization. With this goal in mind, we use a simple sky model for our simulations that only includes bright radio point sources. No diffuse emission from the Galaxy is included as a part of the sky model, and the 21 cm signal and thermal noise are also omitted. Our sky model is derived from the log N–log S distribution of sources and is termed the "Global Sky Model" (GSM) from now onward. Since the GSM only includes sources above 1 Jy, we follow the source counts from the 6C survey at 151 MHz (Hales et al. 1988):

$$\begin{equation} dN/dS = 3600 S_{{\rm Jy}}^{-2.5} {\rm \;Jy}^{-1} {\rm \;str}^{-1}. \end{equation} \tag{ 1 }$$

For a field of view of 15° the total number of sources above 1 Jy is ∼170, following the above power-law distribution. The entire flux range, between 1 Jy and 10³ Jy, has been subdivided into several bins (in the logarithmic scale) and populated with the number of sources that corresponds to the flux range of each bin (according to Equation (1)). Inside each bin, we have assigned each source a flux density following a normal distribution. The strongest source in our GSM is ∼200 Jy. The observed distribution of radio sources shows evidence for only very weak angular clustering and the brightest extragalactic sources in the sky are not clustered at all (Blake & Wall 2002). Therefore, in order to assign a position to each of these sources within the field of view, a uniform random number generator has been used which predicted the offset from the field center for respective sources. In the GSM all the foreground sources are flat spectrum, i.e., with zero spectral index (α = 0).

Figure 1 shows a simulated image of our bright source sky model that has been produced using the procedure described below. A wide-field variant of the well-known Clark-CLEAN algorithm that utilizes a w-projection algorithm for 3D imaging was applied to the simulated image. The apparent angular size of the sources in Figure 1 reflects the size of the synthesized beam. The input sources in the model are treated as ideal point sources. This input sky model is used for all the simulations presented in this paper.

Zoom In Zoom Out Reset image size

Table 1 outlines the instrumental parameters that we have assumed for this analysis. Most of these parameters reflect the current specifications for the MWA, but we note that the array is presently under development and some properties may be subject to change. In addition, we have intentionally reduced the simulated field of view compared to the actual MWA in order to reduce the computational overhead of the simulation. Figure 2 shows the array layout for the 512 element array with maximum baseline of 1.5 km.

Zoom In Zoom Out Reset image size

Table 1. Array Specifications

Parameters	Values
No. of tiles	512
Central frequency	158 MHz (z ∼ 8)
Field of view	∼15° at 158 MHz (∝λ)
Synthesized beam	∼45 at 158 MHz (∝λ)
Effective area per tile	∼17 m2
Maximum baseline	∼1.5 km
Total bandwidth	32 MHz
Tsys	∼250 K
Channel width	∼32 kHz
Thermal noise	∼7.55 mK
	(5000 hr and 2.5 MHz)

Notes. Array parameters have been influenced by the MWA specifications as mentioned in Mitchell et al. (2008) and Bowman et al. (2009). The actual MWA field of view is ∼25° at 150 MHz.

Download table as:  ASCIITypeset image

For the purposes of modeling earth rotation synthesis in the instrumental response, the center of the target field is chosen such that it coincides with one of the cold spots in the foreground Galactic synchrotron emission visible from the southern hemisphere location of the MWA. The exact field center used for the GSM is 4 hr in right ascension and −26° in declination. Most of the upcoming low-frequency telescopes, including the MWA, will only be able to observe a field around its transit. We have used 6 hr of integrations for all the simulations, assuming that the telescopes will observe a field between ±3 hr in Hour Angle from transit.

The 15° field of view will include ∼170 bright sources (>1 Jy). The individual flux densities of these foregrounds are ∼10⁵–10⁷ times higher than the signal from cosmic reionization that these instruments are aiming to detect. So the challenge lies in calibration and subsequent removal of such bright sources from the raw data sets. The data rate of 19 GB s⁻¹ (Mitchell et al. 2008) will not allow the MWA to store the raw visibilities produced by the correlator. Hence real-time calibration and imaging needs to be done in order to reduce the data volume and store the final product in the form of image cubes (Mitchell et al. 2008). The critical steps include removal of the bright sources above the S_cut level from the data sets in these iterative rounds of real-time calibration and imaging procedure. As a result the residual image cubes will not be dominated by these bright sources and the rest of the foregrounds can be removed in the image domain.

However, the accuracy of the foreground source removal strategies are strongly dependent on the data reduction procedure. The likely data reduction procedure which will be followed by the upcoming telescopes can be broadly outlined as

• 1.  
  the raw data sets from the correlator will go through real-time calibration and subsequent removal of the bright sources based on some GSM, down to S_cut level, in the UV domain.
• 2.  
  The residual data sets will be imaged and stored as a cube for the future processing and removal of sources which are below S_cut.

The simulated data reduction pathway that we follow in this paper is

• 1.  
  first, the observed visibilities ($V_{ij}^{{\rm Obs}}(\mathbf {u},\nu) \equiv V_{ij}^{{\rm GSM}_{{\rm perfect}}}(\mathbf {u},\nu)$) are simulated for a 6 hr observation (±3 hr in Hour Angle) using the GSM and the array configuration from Section 2.2. In the above notation, u ≡ (u, v) denotes the Fourier conjugate of the sky coordinate (θ_x, θ_y) and ν is the frequency of observations.
• 2.  
  Next, we generate the foreground model (V^mod_ij(u, ν)) that will be subtracted from the observation. In this case, the model is corrupted to either simulate errors in the assumed positions of the sources or to simulate calibration errors. For the source position errors, the model visibilities are given simply by

  $$\begin{equation} V_{ij}^{{\rm mod}}(\mathbf {u},\nu) = V_{ij}^{{\rm GSM}_{{\rm imperfect}}}(\mathbf {u},\nu), \end{equation} \tag{ 2 }$$

  where the position of each source has been slightly moved from its original location by a distance drawn from a Gaussian distribution with standard deviation σ_θ. We assume that the source position errors are constant throughout the entire duration of the observation, as would be the case for a foreground model constructed from either an outside catalog or from the data itself at the conclusion of the observation. This is an idealization that may be broken in practice if sources are "peeled" in real time.For the residual calibration errors, the model visibilities are given by

  $$\begin{equation} V_{ij}^{{\rm mod}}(\mathbf {u},\nu) = g_i(t)g_j^{*}(t)V_{ij}^{{\rm GSM}_{{\rm perfect}}}(\mathbf {u},\nu), \end{equation} \tag{ 3 }$$

  where $g_i(t) \approx (1+a_i)e^{i\phi _i}$ are the antenna-dependent complex gains. The parameters a_i and ϕ_i denote small amplitude and phase deviations, respectively, and are each drawn from their own Gaussian distribution with standard deviation σ_a or σ_ϕ (Datta et al. 2009).The MWA will produce calibration solutions in real time with an ∼8 s cadence. This rapid pace is planned in order to simultaneously calibrate both the instrument hardware properties and the ionospheric phase screen. It is not known, yet, if in practice the residual calibration errors from the real-time processing will be largely independent or highly correlated between individual 8 s solutions. This is an important experimental property to consider in our simulation, because the degree of correlation greatly affects the level of accuracy needed in individual calibration solutions. If the individual errors are largely independent, then each 8 s sample can be modeled as coming from a Gaussian distribution and the accuracy tolerance will be relatively loose since many samples will be available and tend to average toward zero. Such a situation would be the best-case scenario. On the other hand, if the calibration errors are highly correlated, then each calibration solution must meet a much more stringent accuracy level to achieve the same residual contamination at the end of the full integration.For our simulation procedure, we assume a relatively conservative scenario that the residual errors in a given antenna's 8 s calibration solutions are perfectly correlated for the duration of one 6 hr observing night, but perfectly uncorrelated between successive observing nights. We further assume that the residual errors between antennas are perfectly uncorrelated at all times. This choice is somewhat arbitrary given the current level of knowledge, but we believe it is a plausible fiducial case since both the overall ionospheric properties and the ambient conditions may change significantly from day-to-day. Hence, in our simulation, σ_a and σ_ϕ are used to draw a calibration error value (a_i and ϕ_i) from a Gaussian distribution only once per antenna per night and that specific error is applied to all the simulated 8 s solutions for the given antenna throughout the 6 hr period of rotation synthesis. When the next night's observing block commences, a new error is drawn from the distribution for each antenna, and so on.
• 3.  
  Now we are ready to calculate the residual visibilities by subtracting the foreground model produced in step (2) from the simulated observation of step (1) according to

  $$\begin{equation} V_{ij}^{{\rm res}}(\mathbf {u},\nu)=V_{ij}^{{\rm obs}}(\mathbf {u},\nu)-V_{ij}^{{\rm mod}}(\mathbf {u},\nu). \end{equation} \tag{ 4 }$$

  In the sections below, we refer to this step as GSM subtraction and it was implemented using the UVSUB algorithm (Cornwell et al. 1992). For the residual calibration errors, we can reduce Equation (4) by substituting in with Equation (3) and simplifying to obtain

  $$\begin{equation} V_{ij}^{{\rm res}}(\mathbf {u},\nu) = (1 - g_i g_j^*) V_{ij}^{{\rm GSM}_{{\rm perfect}}}. \end{equation} \tag{ 5 }$$
• 4.  
  At this point, we have completed the subtraction of bright sources from our simulated observation, leaving only residual contamination due to the differences between our simulated observation and the corrupted model. Example images of the residual contamination at this stage are shown in Datta et al. (2009, panels (a) of their Figures 5, 7, and 9). In practice, this bright source-subtracted data cube will be the starting point for the second stage of redshifted 21 cm foreground subtraction that aims to remove faint and confused sources by fitting and subtracting a low-order polynomial along the frequency axis for each line of sight in the data cube. We want to understand how this additional process affects the end result of the bright source removal, so we approximate the faint source polynomial fitting here by applying a Fourier transform to the UV map generated from the residual visibilities in order to produce a residual dirty image cube, I^res. In this dirty image cube, we fit a third-order polynomial in frequency along each line of sight and subtract it. Thus, we obtain the final residual image, $I^{{\rm res}}_{{\rm polysub}}(\vec{\theta }, \nu)$. We refer to this step as GSM+polynomial subtraction for the remainder of this paper and it was implemented using the IMLIN algorithm (Cornwell et al. 1992). To illustrate this final result, Figure 3 shows residual spectral profiles along four lines of sight in the dirty image cube after polynomial fitting and subtraction. Example images of the final residual contamination after GSM+polynomial subtraction are also shown in Datta et al. (2009, panels (b) of their Figures 5, 7, and 9).Using higher order polynomials in the GSM+polynomial subtraction step removes structures at increasingly smaller scales (McQuinn et al. 2006). This improves the foreground cleaning, but since the 21 cm reionization signal has significant structures on scales that correspond to ∼2.5 MHz or approximately 10% of the bandwidth over which the polynomial is fit, it also has the potential to remove much of the 21 cm signal. We have restricted our attention to a third-order polynomial in this work because it is the lowest-order polynomial likely to be sufficient for removing the faint continuum sources given their power-law spectral shapes.We also explored using the UVLIN algorithm (Cornwell et al. 1992), which fits and subtracts polynomials in the UV domain instead of the image domain, eliminating the need to convert our residual data sets into image cubes. However, UVLIN works perfectly only within a small field of view, depending on the channel width in frequency (Cornwell et al. 1992), and was found to be inadequate for our purposes.
• 5.  
  The final step in our procedure is to calculate power spectra from the residual image cubes and compare these residual foreground power spectra to the theoretically predicted 21 cm power spectrum and expected thermal noise power spectrum for the MWA. We calculate three forms of the residual power spectra from our final data cubes: the derived angular power spectrum C_ℓ for a narrow frequency channel, the spherically averaged 3D power spectrum P(k) from the entire data cube, and the 2D power spectrum P(k_⊥, k_||) found by averaging over transverse modes in the full 3D power spectrum. Each of these cases is discussed in more detail in Section 3. As a reference, we show in Figure 4 the spherically averaged power spectrum for our input GSM before any source removal has been applied.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to simulate the observed visibilities (V^Obs_ij), we have used the simulator tool in the CASA software.⁵ We have also used CASA to perform the imaging, GSM subtraction, and the subsequent GSM+polynomial subtraction step. The rest of the operations are performed using separately written python⁶ scripts.

In this work, we focus our attention on the residual contamination that can be tolerated in measurements of 21 cm power spectra, rather than direct imaging of the 21 cm background. The 21 cm field from the reionization epoch, as observed by a radio interferometer, is a fully 3D data set, where the two angular dimensions translate to transverse distances on the plane of the sky and the spectral dimension corresponds to a line-of-sight distance within the target field through the redshift of the 21 cm line (modulo small deviations due to peculiar velocities). The full treatment of the power spectrum of such an observation is also 3D.

The primary motivation for seeking to first detect and characterize 21 cm power spectra, rather than immediately attempt to image the background, is that the signal-to-noise ratio (S/N) in a properly chosen power spectrum measurement can be significantly higher than the per-pixel S/N in a direct imaging observation. This is because the information in the 21 cm field is expected to be effectively compressed in a power spectrum by averaging over many Fourier modes.

As we briefly mentioned in step (5) above (Section 2.3), we consider three different power spectra that can be calculated from MWA observations. The first is the angular power spectrum, found by taking the 2D Fourier transform of an image map at any specific frequency in the observed data cube. This type of power spectrum is widely used for characterizing CMB anisotropies and was the first to be proposed in connection to redshifted 21 cm measurements. For 21 cm measurements, however, it has become recognized over the last several years that angular power spectra will be highly susceptible to contamination by diffuse foregrounds (Di Matteo et al. 2002, 2004) and will likely have limited utility for characterizing 21 cm fluctuations. Nevertheless, most existing experiments, including PAPER (Parsons et al. 2010) and the GMRT EOR project (Paciga et al. 2010), use angular power spectra as effective diagnostics of their instrumental performance and foreground subtraction routines. We begin our investigation with angular power spectra in Section 3.2.

More recently, the preferred method of deriving a power spectrum for characterizing redshifted 21 cm fluctuations has been to average over the full 3D 21 cm power spectrum, P(k), in logarithmic shells of constant radius, k, thus yielding P(k) (Morales & Hewitt 2004). This method is possible because the 21 cm field is homogeneous and isotropic, and hence, will have spherical symmetry in the 3D power spectrum (in practice, redshift-space effects due to peculiar velocities distort this symmetry, however; albeit in a reasonably predictable manner). The spherically averaged power spectrum has the advantage of greatly increasing the S/N over the angular power spectrum by effectively reducing the entire observed data set into ∼10 independent measurements (as opposed to reducing only a single image plane in the case of the angular power spectrum). For the MWA, this means that only ∼300 hr of observing are required to have sufficient sensitivity to detect the spherically averaged 21 cm power spectrum at z ≈ 8, assuming the IGM is not fully ionized at that time, whereas ∼5000 hr would be necessary to directly image the 21 cm background. We present results for the spherically averaged power spectrum of an entire simulated MWA data cube in Section 3.3.

Finally, in Section 3.4, we consider an intermediate level between the spherically averaged power spectrum and the full 3D power spectrum. In this case, we hope to facilitate an approach that retains the S/N benefits of the spherically averaged power spectrum, while anticipating some of the possible problems. It has been shown (Bowman et al. 2009) that different regions in the full 3D power spectrum will be corrupted by diffuse foregrounds to varying levels. But in our somewhat "naive" spherically averaged treatment described above, we assume that the entire 3D power spectrum derived from MWA observations contains equally valid measurements. It would be better if we could assess which regions of the full 3D power spectrum to use before performing the average and then discard the regions anticipated to be most corrupted by foreground contamination when averaging. This should be possible with a 2D power spectrum that is calculated in a manner similar to the spherically averaged power spectrum, but, instead of averaging over P(k) in shells of constant k, we average over transverse annuli of constant $k_{\bot }\equiv \sqrt{\vphantom{A^A}\smash{\hbox{${k_x^2 + k_y^2}$}}}$ for each line-of-sight plane k_|| ≡ k_z in the 3D power spectrum. This yields the 2D power spectrum, P(k_⊥, k_||). Because both the diffuse and bright source foregrounds are generally isotropic in the sky, they should exhibit axial symmetry around the line-of-sight axis in the 3D power spectrum, as opposed to the (approximate) spherical symmetry of the 21 cm signal (Morales et al. 2006b). Hence, it should be possible to retain sufficient information to largely isolate the foreground and signal contributions in P(k_⊥, k_||). We anticipate that the 2D power spectrum will likely be the default technique applied to most 21 cm measurements as experiments progress, hence the results in Section 3.4 should be the most pertinent to many readers.

Because power spectra measurements differ significantly from direct imaging, there are several key questions that we seek to address in the remainder of this section: (1) are the tolerances on the source position and calibration errors greater (or lesser) than in the direct imaging case, (2) is a particular region of the power spectrum likely to be more affected by residual contamination than another, and (3) can we hope to build a library of template models for foreground contamination that could be used to marginalize out some of the contamination during the analysis of the power spectrum?

We will address these questions for each of the three classes of power spectra listed above. For each class of power spectrum, we will present residual power spectra for both corruption models: source position errors and calibration errors. And for each corruption model, we will use three fiducial levels of error in our investigation: for source position errors, our fiducial cases are σ_θ = {0.01, 0.1, and 1 arcsec} while for the calibration errors our fiducial levels are σ_a = {0.01%, 0.1%, and 1%} in gain amplitudes, which also translates to σ_p = {001, 01, and 1°} in phase (Datta et al. 2009).

White et al. (1999) describe the technique to derive the angular power spectrum from radio interferometric data. Using the flat-field approximation (Datta et al. 2007)

$$\begin{equation} C_{\ell }=\frac{\sum _{2\pi |\mathbf {u}|=\ell }W(\mathbf {u})|V(\mathbf {u})|^2}{\sum _{2\pi |\mathbf {u}|=\ell }W(\mathbf {u})}, \end{equation} \tag{ 6 }$$

where $|\mathbf {u}| = \sqrt{u^2 + v^2}$ and ℓ ≃ 2π|u| under flat-field approximation. Here, V(u, ν) is the un-weighted visibilities from the residual images and W(u) is the number of visibilities entering each u cell.

In Figures 5 and 6, we have plotted ℓ(ℓ + 1)C_ℓ/(2π) calculated for one frequency bin of width 125 KHz from our simulated residual image maps. Figure 5 shows the angular power spectrum resulting from using the foreground model that is corrupted by source position errors. Figure 6 illustrates the same result for the case of residual calibration errors. The top panels (a) of both figures show the angular power spectra derived from the GSM-subtracted residuals of step (3) in our analysis procedure. The bottom panels (b) show the angular power spectra from the final GSM+polynomial-subtracted residuals following step (4). The shaded region in panel (b) of both figures corresponds to ℓ ⩽ 250. The GSM+polynomial subtraction step, which involves fitting a third-order polynomial over a total bandwidth of 32 MHz, is expected to remove most of the significant structures for scales larger than this. All of the plots have been restricted to ℓ ≲ 5000 to match the size of the MWA synthesized beam (4.5 arcmin).

Zoom In Zoom Out Reset image size

The total thermal noise power is much stronger than the angular power spectra of the ${\rm H\,\mathsc{i}}$ 21 cm signal. Hence, we have assumed that the final power spectrum from the real data will be generated by dividing the observation into different epochs of equal duration and then cross-correlating the data cubes from the two epochs (Bowman et al. 2009). This approach preserves the persistent ${\rm H\,\mathsc{i}}$ 21 cm signal and eliminates the thermal noise power (which will be independent between the two observing epochs and, therefore, average to zero during the cross-correlation), leaving only the thermal uncertainty. Hence, the relevant noise figure for the angular power spectra measurement is given by

$$\begin{equation} C_{\ell }^N = \left\langle\frac{ \sum _{\ell } \left| N_1(\ell)^{\ast } N_2(\ell) + N_2(\ell)^{\ast } N_1(\ell) \right|}{2} \right\rangle, \end{equation} \tag{ 7 }$$

where N₁ and N₂ are simulated noise measurements from two different epochs (Bowman et al. 2009).

The residual angular power spectra in the figures can be compared to the thermal noise uncertainty in the observations (Equation 7) and a fiducial 21 cm signal. We plot the expected thermal noise uncertainty angular power spectrum of the MWA after 5000 hr of integration assuming a system temperature of T_sys = 250 K, channel width of 125 KHz, and the observing strategy described in Section 2.2. The thermal uncertainty spectrum is shown assuming the angular power spectrum has been binned in logarithmic intervals of width Δℓ = 0.1 or approximately 10 bins per decade. For the reference ${\rm H\,\mathsc{i}}$ 21 cm signal, we show the power spectrum for a fully neutral IGM at z = 8 (Furlanetto et al. 2006). Modeling the 21 cm signal using a fully neutral IGM provides a reasonable fiducial expectation since recent reionization simulations (Lidz et al. 2008) show that the amplitude of the power spectrum over the scales probed by the MWA is likely to be even larger than the fully neutral level when the universe if roughly 50% ionized. It should be noted that different models predict different amplitudes for the ${\rm H\,\mathsc{i}}$ power spectrum. For simplification, we have used this single realistic model to compare with our residual power spectrum. The specific conclusions regarding the scale-size dependence of where residual power will dominate the 21 cm signal will change depending on the reionization model.

Figure 5(a) shows that the angular power spectrum from the GSM-subtracted images are well above the thermal uncertainty power spectrum, as well as the model ${\rm H\,\mathsc{i}}$ 21 cm signal power spectrum. In Figure 5(b), it is evident that the residual angular power in the GSM+polynomial-subtracted image is greatly reduced; and for two of our fiducial source position error levels (σ_θ = 0.1 and 0.01 arcsec), the angular power spectra intercepts the ${\rm H\,\mathsc{i}}$ signal power spectrum at ℓ ≲ 700. This shows that the GSM+polynomial-subtracted step is very crucial not only for removing faint and confused continuum foreground sources, but also for removing residual power left over after subtracting the bright foreground sources. A source position accuracy of ≲0.1 arcsec would allow the detection of ${\rm H\,\mathsc{i}}$ 21 cm signal at 250 ≲ ℓ ≲ 600 scales.

Similar features are seen in Figure 6 for the case of calibration errors, where the residual angular power spectrum from the GSM-subtracted image only intercepts the thermal noise power spectrum near ℓ ∼ 2000 and only the σ_a = 0.01 % crosses below the model ${\rm H\,\mathsc{i}}$ 21 cm signal power spectrum and the thermal uncertainty spectrum. Again, from Figure 6(b), it is evident that the residual angular power spectrum from the GSM+polynomial-subtracted image is much lower, particularly below the ℓ = 250 threshold. We have not investigated in detail how scales larger than this threshold will be affected by the polynomial subtraction, but it is likely that some of the signal will be removed, as well. A calibration accuracy of σ_a ≲ 0.05% should allow the detection of the ${\rm H\,\mathsc{i}}$ 21 cm signal.

Zoom In Zoom Out Reset image size

The spherically averaged 3D 21 cm power spectrum is the primary reionization observable targeted by the MWA. There has been extensive research on the statistical EoR power spectrum measurement of the brightness temperature fluctuations in low-frequency, wide-field radio observations. Detailed formulation has been developed in the literature by Morales & Hewitt (2004) and Zaldarriaga et al. (2004). The approaches described in these efforts are inspired by the techniques that have been employed successfully for interferometric measurements of CMB anisotropies (White et al. 1999; Hobson & Maisinger 2002; Myers et al. 2003). The primary approach is to convert the full 3D measurement cube to a one-dimensional (1D) power spectrum.

The first step is to transform our residual image cubes $I(\vec{\theta }, \nu)$ into V(u, η) by performing a 3D Fourier transform denoted by the operator $\mathbf {F}(\lbrace \mathbf {u},\eta \rbrace,\lbrace \vec{\theta }, \nu \rbrace)$. It should be noted here that before performing the Fourier transform, we have changed the units of the residual images from flux unit (Jy beam⁻¹) to brightness temperature unit (mK). Hence, we get

$$\begin{equation} V(\mathbf {u},\eta)=\mathbf {F}(\lbrace \mathbf {u},\eta \rbrace,\lbrace \vec{\theta }, \nu \rbrace)I(\vec{\theta }, \nu), \end{equation} \tag{ 8 }$$

where u, η ≡ (u, v, η). Then, we transform the measurement coordinates u, ν into the cosmological coordinates k:

$$\begin{eqnarray} V(\mathbf {k})&=& \mathbf {J}(\mathbf {k},\lbrace \mathbf {u},\eta \rbrace)V(\mathbf {u},\eta) \\ &=& \mathbf {J}(\mathbf {k},\lbrace \mathbf {u},\eta \rbrace)\mathbf {F}(\lbrace \mathbf {u},\eta \rbrace,\lbrace \vec{\theta }, \nu \rbrace)I(\vec{\theta }, \nu), \end{eqnarray} \tag{ 9 }$$

where J(k, {u, η}) denotes the Jacobian of the coordinate transformation from u, η (in units of λ and Hz⁻¹) to k (in units of Mpc⁻¹). We have mainly followed the definition in Peebles (1993) and the formulation detailed in Morales & Hewitt (2004). Hence, we transformed a residual image cube (in sky coordinates) to a 3D residual visibility cube in the Fourier conjugate coordinates of comoving Mpc.

Assuming isotropy of space and ignoring redshift-space distortions inherent in converting our observed data cube to cosmological coordinates, the power spectrum can be taken as approximately spherically symmetric in cosmological k ≡ (k_x, k_y, k_z) coordinates. Hence, the power spectrum can be approximated to the square of the V(k), averaged over spherical shells:

$$\begin{equation} P(k) = \langle | V(\mathbf {k}) |^2 \rangle_{|\mathbf {k}|=k}. \end{equation} \tag{ 11 }$$

Thus, we obtain the 1D total power spectrum (Morales & Hewitt 2004) or the more common dimensionless power spectrum given by Δ² = k³P(k)/(2π²).

While deriving the 1D power spectrum, we have weighted the individual measurements |V(k)|² by the per cell visibility contributions. This scheme is similar to the natural weighting scheme which is applied to the raw visibilities before imaging, and follows the form

$$\begin{equation} P(k) = \frac{\sum _{|\mathbf {k}|=k} W_u(\mathbf {k})\left| V(\mathbf {k}) \right|^2 }{\sum _{|\mathbf {k}|=k} W_u(\mathbf {k})}, \end{equation} \tag{ 12 }$$

where W_u(k) denotes the total number of visibilities contributing per k cell. Here, we should explicitly mention that the V(k) used in the above equation are the un-weighted visibilities obtained from the residual images.

The total thermal noise power is much stronger than the 1D spherically averaged power spectra of the ${\rm H\,\mathsc{i}}$ 21 cm signal. Hence, similar to the angular power spectrum case, we compare our results with the thermal noise uncertainty given by

$$\begin{equation} P^N(k) = \left<\frac{ \sum _{|\mathbf {k}|=k} \left| N_1(\mathbf {k})^{\ast } N_2(\mathbf {k}) + N_2(\mathbf {k})^{\ast } N_1(\mathbf {k}) \right|}{2} \right>, \end{equation} \tag{ 13 }$$

where N₁ and N₂ are simulated noise measurements from two different epochs (Bowman et al. 2009).

Figures 7 and 8 show the 1D spherically averaged power spectrum from the residual images. As with the angular power spectrum, these figures also show theoretical ${\rm H\,\mathsc{i}}$ 21 cm power spectrum. However, here, instead of using a total thermal noise power spectrum as we did for the angular power spectrum plots, we show the spherically averaged thermal noise uncertainty power spectrum from 300 hr of observation with the MWA, as mentioned in Equation 13. The thermal uncertainty spectrum is shown assuming that the spherically averaged power spectrum has been binned in logarithmic shells of width Δk/k = 0.5 or approximately five bins per decade. As discussed in Lidz et al. (2008), the MWA-512 will be sensitive primarily to scales 0.1 ≲ k ≲ 1 Mpc⁻¹. The shaded regions in Figures 7(b) and 8(b) are at k ⩽ 0.03 Mpc⁻¹ and indicate the scales where the GSM+polynomial subtraction step removes significant power. These regions correspond to the ℓ ⩽ 250 threshold in the angular power spectra plots. The higher end of the k value for the residual power spectrum is restricted due to the cell size in the image domain and frequency resolution of the channels in the residual image cube. The maximum value of k is attained along the k_z axis only, and hence few or no transverse (angular) modes contribute to the power spectrum at small scales above k ≳ 0.6 Mpc⁻¹ in Figures 7 and 8. The angular resolution of the MWA of ∼4.5 arcmin (synthesized beam) corresponds to k ≈ 0.6 Mpc⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the 1D spherically averaged power spectra from the residual images. These are the residual images after the foreground subtraction in the presence of source position errors. Figure 7(a) shows that the 1D power spectra from the GSM-subtracted image with σ_θ = 0.01 arcsec are only below the thermal uncertainty power spectrum and the ${\rm H\,\mathsc{i}}$ signal power spectrum. In Figure 7(b), it is evident that the residual 1D spherically averaged power spectra with σ_θ = 0.01 and 0.1 arcsec from the final GSM+polynomial-subtracted image are below the thermal uncertainty power spectrum and the ${\rm H\,\mathsc{i}}$ signal power spectrum. Hence, a source position accuracy of σ_θ ≲ 0.1 arcsec would allow the detection of ${\rm H\,\mathsc{i}}$ 21 cm signal with the MWA.

Turning to the case of the calibration errors, Figure 8(a) shows that only the 1D power spectrum from the GSM-subtracted image for calibration error of 0.01% is well below the thermal uncertainty power spectrum and the ${\rm H\,\mathsc{i}}$ signal power spectrum. In Figure 8(b), it is evident that the 1D spherically averaged power spectra with σ_a = 0.01% and 0.1% from the GSM+polynomial-subtracted images are below the thermal uncertainty power spectrum and the theoretical ${\rm H\,\mathsc{i}}$ signal power spectrum. Hence, the residual calibration accuracy of σ_a ≲ 0.05% would allow the detection of ${\rm H\,\mathsc{i}}$ 21 cm signal.

In comparison to the angular power spectrum, we can infer that the 1D spherically averaged power spectra have a better tolerance for both the source position and residual calibration errors. This also reflects the fact that the angular power spectrum has been produced using a single channel map of 125 KHz, whereas the 1D spherically averaged spectrum is produced with the total bandwidth of 32 MHz.

In the previous section, we showed the analysis of the 1D spherically averaged power spectrum. However, this formulation mixes the contribution from the k_⊥ and k_|| directions. It is useful, therefore, to break the averaging from the 3D k-space to the 1D k-space into two steps since both the foregrounds and a full treatment of the predicted redshifted 21 cm signal that includes redshifted-space distortions have aspherical structure in the Fourier domain (Morales et al. 2006b). Following McQuinn et al. (2006), we average over the transverse (angular) direction in the full three-dimension power spectrum to obtain P(k_⊥, k_||). We obtain the 2D power spectrum based on the maximum likelihood formalism following the same approach as used for the spherically averaged power spectrum in Equations (11) and (12).

Figures 9–12 illustrate the results of the simulation for the 2D residual power spectra. We show in Figure 9(a) the estimated 2D thermal noise uncertainty after 300 hr of integration with the MWA on our target field. Figure 9(b) shows a theoretical 2D ${\rm H\,\mathsc{i}}$ signal power spectrum, P(k_⊥, k_||), of the ${\rm H\,\mathsc{i}}$ signal in units of mK² Mpc⁻³. Figures 10 and 11 show the 2D power spectra of the residual image cubes for source position errors and calibration errors, respectively. In this Section, we have analyzed residual images for only one of our fiducial error levels for each type of model corruption. For source position errors, we use σ_θ = 0.1 arcsec and for residual calibration errors, we use σ_a = 0.1%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4.1. The Wedge

3.4.2. Comparison to Expected Signal