Bright Source Subtraction Requirements For Redshifted 21 cm Measurements
Abstract
The Hi 21 cm transition line is expected to be an important probe into the cosmic dark ages and epoch of reionization. Foreground source removal is one of the principal challenges for the detection of this signal. This paper investigates the extragalactic point source contamination and how accurately bright sources ( Jy) must be removed in order to detect 21 cm emission with upcoming radio telescopes such as the Murchison Widefield Array (MWA). We consider the residual contamination in 21 cm maps and power spectra due to position errors in the skymodel for bright sources, as well as frequency independent calibration errors. We find that a source position accuracy of arcsec will suffice for detection of the Hi power spectrum. For calibration errors, accuracy in antenna gain amplitude is required in order to detect the cosmic signal. Both sources of subtraction error produce residuals that are localized to small angular scales, Mpc, in the twodimensional power spectrum.
Subject headings:
Cosmology: Early Universe, Galaxies: Intergalactic Medium, Radio Lines: General, Techniques: Interferometric, Methods: Data Analysis1. Introduction
The cosmological epoch of reionization (EoR) is a key milestone in the history of structure formation, marking the transition from a fully neutral to a highly ionized intergalactic medium (IGM) due to the ultraviolet and Xray radiation of early stars, galaxies, and black holes. Recent observations of the GunnPeterson effect, i.e., Ly absorption by the neutral IGM, toward the most distant quasars (), and the large scale polarization of the CMB, corresponding to Thompson scattering during reionization, have set the first constraints on the reionization process. These results suggest significant variance in both space and time, starting perhaps as far back as (Komatsu et al. (2010); WMAP 7 year data) and extending to (Fan et al., 2006). Previous WMAP five year data indicates the 5 detection of the Emode of polarization which rules out any instantaneous reionization at at 3.5 level. In case of the GunnPeterson effect, the IGM becomes optically thick to Ly absorption for a neutral fraction as small as . In order to overcome these limitations, it has been widely recognized that mapping the redshifted Hi 21 cm line has great potential for direct studies of the neutral IGM during reionization (Barkana & Loeb, 2001; Fan et al., 2006; Furlanetto et al., 2006; Morales & Wyithe, 2009).
There are number of upcoming lowfrequency arrays with key science goals to detect the Hi 21 cm signal from the EoR. This includes the Murchison Widefield Array [MWA] (Mitchell et al., 2008; Lonsdale et al., 2009; Bowman et al., 2006), Precision Array to Probe Epoch of Reionization [PAPER] (Backer et al., 2007; Parsons et al., 2009), Low Frequency Array [LOFAR] (Harker et al., 2010; Jelić et al., 2008; Labropoulos et al., 2009) and Giant Meterwave Radio Telescope [GMRT] (Pen et al., 2009). One of the major challenges for all of these upcoming arrays will be the removal of the continuum foreground sources in order to detect the faint Hi signal from the EoR.
A variety of continuum foregrounds complicate redshifted 21 cm measurements of the EoR (Shaver et al., 1999). Diffuse Galactic synchrotron emission dominates the lowfrequency radio sky and is approximately four orders of magnitude brighter than the mK 21 cm signal at the frequency relevant to reionization ( MHz). In addition, Galactic and extragalactic freefree emission contribute additional flux to the diffuse foreground. Radio point sources from AGN, radio galaxies, and local Galactic sources are numerous and particularly challenging. The brightest of these sources have fluxes well above Jy and are seven or eight orders of magnitude above the EoR signal in lowfrequency radio maps. The distribution of point sources also extends to very faint levels such that the brightness temperature due to confused sources in upcoming arrays will be K, or three orders of magnitude brighter the than the 21 cm background.
In this paper we discuss how the radio interferometric imaging techniques are going to affect the foreground source modeling and subsequent removal from the dataset in order to search for the EoR signal. Recently, there has been extensive research on foreground source modeling at these low frequencies (Di Matteo et al., 2002; Jelić et al., 2008; Thomas et al., 2009). Similar effort has also been made in exploring different techniques to remove the foregrounds from the EoR dataset by Morales et al. (2006a, b); Bowman et al. (2009); Liu et al. (2009); Parsons et al. (2010); Gleser et al. (2008); Harker et al. (2009b, a). Since attempting to observe a signal below the confusion limit of foreground sources is a novel aspect of 21 cm experiments, most of these works primarily focus on the removal of faint and confused sources that fall below a specified cutoff flux limit, ( Jy). They do not consider the foreground sources brighter than and how accurately they need to be removed. Indeed, most of these analysis implicitly assume that the bright foreground sources above have been removed perfectly. But in reality imperfect instrument calibration or any errors in the subtracted foreground model will introduce artifacts and leave residual contamination in the data after bright source removal, even by traditional techniques such as “peeling”. These residuals may interfere with either the subsequent faint source subtraction or the ultimate detection and characterization of the redshifted 21 cm signal.
In Datta et al. (2009) we dealt with the bright point sources above and the limitations that will be caused due to imperfect removal of such sources in the image plane. In this paper we extend the initial analysis in order to estimate the residual contamination in the power spectral domain of improper bright source subtraction. The objective of this paper is to demonstrate how the accuracy in the foreground removal affects the detection of Hi 21 cm power spectra with the MWA. In Section 2, we discuss our choice of sky model and outline the simulation parameters, including the array specifications and data reduction procedure, and describe the two categories of corruption terms that we will consider: source position errors and residual calibration errors. The results obtained for the residual angular power spectrum, spherically averaged threedimensional power spectrum, and twodimensional power spectrum are presented in Section 3. Finally, in the last section we summarize the implications of the results from our simulations.
2. The Simulation
2.1. Sky Model
Our main aim is to explore the level of accuracy needed in instrument calibration and foreground modeling in order to ensure 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 – distribution of sources and is termed the “Global Sky Model” (GSM) from now onwards. 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):
(1) 
For a fieldofview of the total number of sources above 1 Jy is , following the above powerlaw distribution. The entire flux range, between 1 Jy and Jy, has been subdivided into several bins (in 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 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 fieldofview, 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 ().
Figure 1 shows a simulated image of our bright source sky model that has been produced using the procedure described below. A widefield variant of the wellknown ClarkCLEAN algorithm that utilizes a wprojection algorithm for 3dimensional 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.
2.2. Array Specifications
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.
Parameters  Values 

No. of Tiles  512 
Central Frequency  158 MHz () 
Field of View  15 at 158 MHz. () 
Synthesized beam  4.5’ at 158 MHz. () 
Effective Area per Tile  17 
Maximum Baseline  1.5 km 
Total Bandwidth  32 MHz 
250 K  
Channel Width  32 kHz 
Thermal Noise  7.55 mK 
(5000 hours MHz) 
Note. – Array parameters have been influenced by the MWA specifications as mentioned in Mitchell et al. (2008) and Bowman et al. (2009). Original MWA Fieldofview is at 150 MHz.
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 hours in Right Ascension and 26 degree in Declination. Most of the upcoming lowfrequency telescopes, including the MWA, will only be able to observe a field around its transit. We have used 6 hours of integrations for all the simulations, assuming that the telescopes will observe a field between 3 hours in Hour Angle from transit.
2.3. Data Reduction Procedure
The fieldofview will include bright sources ( Jy). The individual flux densities of these foregrounds are 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 datasets. 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 realtime 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 level from the datasets in these iterative rounds of realtime calibration and imaging procedure. As a result the residual imagecubes 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 :

The raw datasets from the correlator will go through realtime calibration and subsequent removal of the bright sources based on some Global Sky Model (GSM), down to level, in the UV domain.

The residual datasets will be imaged and stored as a cube for the future processing and removal of sources which are below .
The simulated data reduction pathway that we follow in this paper is:
(i) First, the observed visibilities () are simulated for a 6 hour observation ( hours in Hourangle) using the GSM and the array configuration from Section 2.2. In the above notation, denotes the Fourier conjugate of the sky coordinate () and is the frequency of observations.
(ii) Next, we generate the foreground model () 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:
(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 realtime.
For the residual calibration errors, the model visibilities are given by:
(3) 
where are the antennadependent complex gains. The parameters and denote small amplitude and phase deviations, respectively, and are each drawn from there own Gaussian distribution with standard deviation or (Datta et al., 2009).
The MWA will produce calibration solutions in realtime with an second 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 realtime processing will be largely independent or highly correlated between individual 8second 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 8second 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 bestcase 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 8second calibration solutions are perfectly correlated for the duration of one 6hour 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 daytoday. Hence, in our simulation, and are used to draw a calibration error value ( and ) from a Gaussian distribution only once per antenna per night and that specific error is applied to all the simulated 8second solutions for the given antenna throughout the 6hour 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.
(iii) Now we are ready to calculate the residual visibilities by subtracting the foreground model produced in step (ii) from the simulated observation of step (i) according to:
(4) 
In the sections below, we refer to this step as “UVSUB” since 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:
(5) 
(iv) 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 sourcesubtracted 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 loworder 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, . 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, . We refer to this step as “IMLIN” for the remainder of this paper because it was implemented with 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 IMLIN are also shown in Datta et al. (2009, panels (b) of their Figures 5, 7, and 9).
Using higher order polynomials in the IMLIN 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 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 lowestorder polynomial likely to be sufficient for removing the faint continuum sources given their powerlaw 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 fieldofview, depending on the channel width in frequency (Cornwell et al., 1992), and was found to be inadequate for our purposes.
(v) 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 for a narrow frequency channel, the spherically averaged threedimensional power spectrum from the entire data cube, and the twodimensional power spectrum found by averaging over transverse modes in the full threedimensional 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.
In order to simulate the observed visibilities (), we have used the simulator tool in the CASA software ^{1}^{1}1http://casa.nrao.edu/. We have also used CASA to perform the imaging and the subsequent IMLIN step. The rest of the operations are performed using separately written python ^{2}^{2}2http://www.python.org/ scripts.
3. Results
We begin our discussion of the results of the residual power spectrum determination by reviewing our initial findings from Datta et al. (2009). In that work, we explored the source position and calibration accuracy needed to allow direct imaging of Stromgren spheres with very deep integrations by the MWA. Our simulations demonstrated that knowledge of the true positions of the bright foreground sources in an MWA target field is required to within arcsec, assuming Gaussian errors, in order for the residual contamination following subtraction to be below the 21 cm signal from Stromgren spheres in image maps that could be acquired by the MWA with 5000 hours of integration. Similarly, in Datta et al. (2009) we found that, for the case of calibration errors corrupting the measurements under the same conservative assumptions outlined in step (ii), a calibration accuracy of in gain amplitude (or degree in phase) is needed for the residual contamination to be below the thermal noise in a part of the image map far from any bright sources for a long integration by the MWA.
Here, we focus our attention on the residual contamination that can be tolerated in measurements of the power spectrum of a target field observed by the MWA. One of the primary motivations for seeking to first detect and characterize the 21 cm power spectrum, rather than immediately attempt to image the background, is that the MWA will only require hours of observing to have sufficient sensitivity to detect the sphericallyaveraged 21 cm power spectrum at , assuming the IGM is not fully ionized at that time. Because the power spectrum measurement differs significantly from a direct imaging observation, there are several key questions that we seek to address: 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 in step (v) (section 2.3) that the MWA is likely to study: angular, spherically averaged, and twodimensional. 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 0.01, 0.1, and 1 in gain amplitudes, which also translates to 0.01, 0.1, and 1 in phase (Datta et al., 2009).
3.1. Angular Power Spectrum
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):
(6) 
where and under flatfield approximation. Here, is the unweighted visibilities from the residual images and is the number of visibilities entering each cell.
In Figures 5 and 6, we have plotted 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 UVSUB residuals of step (iii) in our analysis procedure. The bottom panels (b) show the angular power spectra from the final IMLIN residuals following step (iv). The vertical lines in panel (b) of both figures corresponds to . The IMLIN 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 (corresponding to ). All of the plots have been restricted to to match the size of the MWA synthesized beam (4.5 arcmin).
The total thermal noise power is much stronger than the angular power spectra of the Hi 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 crosscorrelating the data cubes from the two epochs (Bowman et al., 2009). This approach preserves the persistent Hi 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 crosscorrelation), leaving only the thermal uncertainty. Hence, the relevant noise figure for the angular power spectra measurement is given by:
(7) 
where and 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 hours of integration assuming a system temperature of K, channel width of 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 , or approximately ten bins per decade. For the reference Hi 21 cm signal, we show the power spectrum for a fully neutral IGM at (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 Hi power spectrum. For simplification, we have used this single realistic model to compare with our residual power spectrum. The specific conclusions regarding the scalesize 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 UVSUB images are well above the thermal uncertainty power spectrum, as well as the model Hi 21 cm signal power spectrum. In Figure 5(b), it is evident that the residual angular power in the IMLIN image is greatly reduced; and for two of our fiducial source position error levels ( and arcsec), the angular power spectra intercepts the Hi signal power spectrum at . This shows that the IMLIN 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 arcsec would allow the detection of Hi 21 cm signal at scales.
Similar features are seen in Figure 6 for the case of calibration errors, where the residual angular power spectrum from the UVSUB image only intercepts the thermal noise power spectrum near and only the crosses below the model Hi 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 IMLIN image is much lower, particularly below the 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 should allow the detection of the Hi 21 cm signal.
3.2. 1D SphericallyAveraged Power Spectrum
The sphericallyaveraged threedimensional 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 lowfrequency, widefield 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 threedimensional measurement cube to a one dimensional power spectrum.
The first step is to transform our residual image cubes into by performing a three dimensional Fourier transform denoted by the operator . It should be noted here that before performing the Fourier transform, we have changed the units of the residual images from flux unit () to brightness temperature unit (mK). Hence, we get:
(8) 
where . Then, we transform the measurement coordinates into the cosmological coordinates .
(9)  
(10) 
where denotes the Jacobian of the coordinate transformation from (in units of and ) to (in units of ). 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 three dimensional residual visibility cube in the Fourier conjugate coordinates of comoving Mpc.
Assuming isotropy of space and ignoring redshiftspace distortions inherent in converting our observed data cube to cosmological coordinates, the power spectrum can be taken as approximately spherically symmetric in cosmological coordinates. Hence, the power spectrum can be approximated to the square of the , averaged over spherical shells:
(11) 
Thus, we obtain the one dimensional total power spectrum (Morales & Hewitt, 2004), or the more common dimensionless power spectrum given by .
While deriving the 1D power spectrum, we have weighted the individual measurements 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:
(12) 
where denotes the total number of visibilities contributing per cell. Here, we should explicitly mention that the used in the above equation are the unweighted visibilities obtained from the residual images.
The total thermal noise power is much stronger than the 1D spherically averaged power spectra of the Hi 21 cm signal. Hence, similar to the angular power spectrum case, we compare our results with the thermal noise uncertainty given by:
(13) 
where and 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 Hi 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 sphericallyaveraged thermal noise uncertainty power spectrum from 300 hours of observation with the MWA, as mentioned in Equation 13. The thermal uncertainty spectrum is shown assuming the sphericallyaveraged power spectrum has been binned in logarithmic shells of width , or approximately five bins per decade. As discussed in Lidz et al. (2008), the MWA512 will be sensitive primarily to scales Mpc. The vertical lines on Figures 7(b) and 8(b) are at Mpc and indicate the scales below which the IMLIN polynomial fitting step removes significant power. These lines correspond to the threshold in the angular power spectra plots. The higher end of the 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 imagecube. The maximum value of is attained along the axis only, and hence few or no transverse (angular) modes contribute to the power spectrum at small scales above Mpc in Figures 7 and 8. The angular resolution of the MWA of arcmin (synthesized beam) corresponds to Mpc.
Figure 7 shows the 1D sphericallyaveraged power spectra from the residual images. These are the residual images after the foreground subtraction in presence of source position errors. Figure 7(a) shows that the 1D power spectra from the UVSUB image is well below the thermal uncertainty power spectrum and the Hi signal power spectrum. In Figure 7(b), it is evident that the residual 1D sphericallyaveraged power spectra with and arcsec from the final IMLIN image are below the thermal uncertainty power spectrum and the Hi signal power spectrum. Hence, a source position accuracy of arcsec would allow the detection of Hi 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 UVSUB image for calibration error of is well below the thermal uncertainty power spectrum and the Hi signal power spectrum. In Figure 8(b), it is evident that the 1D sphericallyaveraged power spectra with and from the IMLIN images are below the thermal uncertainty power spectrum and the theoretical Hi signal power spectrum. Hence, the residual calibration accuracy of would allow the detection of Hi 21 cm signal.
In comparison to the angular power spectrum, we can infer that the 1D sphericallyaveraged power spectra has 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 sphericallyaveraged spectrum is produced with the total bandwidth of 32 MHz.
3.3. Twodimensional Power Spectrum
In the previous section, we showed the analysis of the 1D sphericallyaveraged power spectrum. However, this formulation mixes the contribution from the and directions. It is useful, therefore, to break the averaging from the three dimensional space to the one dimensional space into two steps since both the foregrounds and a full treatment of the predicted redshifted 21 cm signal that includes redshiftedspace distortions have aspherical structure in the Fourier domain. Following McQuinn et al. (2006), we first average over the transverse (angular) direction in the full threedimension power spectrum to obtain , . This is conceptually similar to averaging over the values and keeping the values in a CMB analysis, except we still have the lineofsight dimension. Next, we obtain the 2D power spectrum based on the maximum likelihood formalism following the same approach as used for the sphericallyaveraged power spectrum in Equations 11 and 12.
Figures 9 through 11 illustrate the results of the simulation for the 2D power spectrum. In all panels of these figures, the color scale is held constant to facilitate comparison between the plots. We show in Figure 9(a) a realization of the 2D thermal noise uncertainty after 300 hours of integration with the MWA on our target field. Figure 9(b) shows the 2D Hi signal power spectrum, , of the Hi 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. It should be noted that these plots are in different units than Figure 7 and 8. 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 arcsec and for residual calibration errors, we use .
Following our convention, Figure 10(a) shows that the 2D power spectrum from the UVSUB image, while Figure 10(b) shows the same for the IMLIN image. From Figure 10(a), we can conclude that the source position errors are more localized towards higher values. Based on the Hi signal shown in Figure 9(b), we see that the signal dominates over the residuals at Mpc. For the IMLIN image in Figure 10(b), the power spectrum shows that Hi signal dominates at Mpc. The thermal uncertainty map (Figure 9(a)) shows that thermal uncertainty dominates over the source position errors at lower values, as mentioned above. There are also regions with high (at Mpc) where the thermal uncertainty dominates over source position errors.
Figure 11 shows a very comparable pattern arising from the residual calibration errors. Hence, we find that the GSM position accuracy of arcsec and calibration accuracy of is sufficient to detect the Hi 21 cm signal in 2D power spectrum. The advantage of the 2D spectrum over the 1D spherically averaged power spectrum is that one can even search for Hi signal at scales around Mpc along the axes, which is fairly clean at lower values. However, in the 1D power spectrum similar scales are dominated by the residual errors due to combined contribution from and .
4. Conclusion
With the results from Section 3, we can address our three key questions for bright source subtraction residuals in the power spectral domain. First, we find that the level of the source subtraction accuracy required for power spectral detection of the 21 cm signal is roughly comparable to the accuracy required for direct imaging of the Hi signal (Datta et al., 2009). The power spectrum tolerance does suffer, however, compared to the direct imaging case in one regard. For direct imaging, it is possible to find areas in the final image map that are far from bright sources and have very low residual contamination. Whereas, for the power spectrum analysis, we use the entire image map for the calculation, mixing both the good and the bad areas in the image cube. Because of this difference, the power spectrum analysis requires a more stringent calibration accuracy of compared to for direct imaging. This problem is most pronounced for the angular power spectrum analysis since the residual contamination is dominated by angular power. It is mitigated partially in the 1D sphericallyaveraged power spectrum by the inclusion of spectral information–with its lower residual contamination power–in the calculation, and also in the 2D power spectrum. However, in the case of the 1D sphericallyaveraged power spectrum, it should be noted that the calibration accuracy is determined on the basis of a hrs of integration by MWA512. If we increase the amount of observing time, the requirement on the accuracy of calibration can be lowered because our model allows the calibration error to average toward zero with additional nights of observations.
The 2D power spectrum, in particular, addresses our second key question, showing clear advantages for separating the residual contamination from the desired signal through distinct localization of the respective contributions in the the and plane. The results for both source position errors and residual calibration errors indicate that at Mpc, we are able to probe most of the scales where the Hi signal is dominant over the residual errors. In the 1D power spectrum we see dominant contribution from the residual errors around Mpc, which can be probed in the 2D power spectra along the axes.
Finally, our third key question was whether we might expect to build a template library of residual contamination errors in the power spectrum domain in order to facilitate interpretation of the final power spectrum. The results of this work indicate that it will indeed be possible. Further, the significant similarities between the 2D power spectra for the source position error case and the gain calibration error case suggest that there may be common and easily identifiable properties of bright source residual contamination that are largely independent of the specific error causing the contamination. This will be a valuable tool for upcoming experiments.
We have found that the IMLIN polynomial subtraction step is crucial not only for faint point source and diffuse foreground subtraction as studied in other works, but also for the success of the bright source foreground subtraction that we explored in this paper.
For the simulations included in this paper, we have performed the foreground subtraction of bright sources from a dataset of a minimum of 6 hours of observation (extrapolated to 300 and 5000 hours) in order to have the full effect of earth rotation synthesis. However, the MWA may perform much of its bright source removal over much shorter timescales ( minutes or less) as part of its realtime calibration pipeline. The major implication for shorter timescale removal of the foregrounds would be to break our assumption when modeling source position errors that the position errors are constant for the entire observation.
We also made the assumption in this paper that each antenna’s calibration errors are perfectly correlated for an entire 6hour observation night, but uncorrelated between observing nights. If the residual calibration error where instead perfectly random between every 8second cycle of the realtime calibration, then we estimate it could be possible to achieve the desired residual contamination noise level and detect the redshifted 21 cm Hi signal from reionization with a significantly larger calibration error of . We have not performed our detailed simulation under this assumption, however, nor have we used the exact parameters that will be employed for the realtime calibration pipeline of the MWA.
The results from the angular power spectrum puts more stringent contraints on the accuracies in source position and calibration. Using a larger chunk ()of frequency width might have reduced the constraints. However, we conclude that if a wider bandwidth is available, it is more advantageous to perform a 1D spherically averaged power spectrum than the angular power spectrum. This is due to the inclusion of the axes contribution to the power spectrum.
Lastly, we would like to emphasis that similar constraints can also be derived for other upcoming arrays, such as LOFAR and PAPER, as well as for future arrays like the Square Kilometer Array or a lunar array. But detailed simulations with the unique array specifications for each instrument would be required, which is beyond the scope of this paper. We expect to build on our present analysis in future work by exploring other arrays, addressing the modified scenarios described above, and including additional calibration issues such as widefield gain calibration of the primary beam and ionosphere.
References
 Backer et al. (2007) Backer, D. C., Parsons, A., Bradley, R., Parashare, C., Gugliucci, N., Mastrantonio, E., Herne, D., Lynch, M., Wright, M., Werhimer, D., Carilli, C., Datta, A., & Aguirre, J. 2007, in Bulletin of the American Astronomical Society, Vol. 38, Bulletin of the American Astronomical Society, 967–+
 Barkana & Loeb (2001) Barkana, R. & Loeb, A. 2001, Phys. Rep., 349, 125
 Blake & Wall (2002) Blake, C. & Wall, J. 2002, MNRAS, 329, L37
 Bowman et al. (2006) Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2006, ApJ, 638, 20
 Bowman et al. (2009) —. 2009, ApJ, 695, 183
 Cornwell et al. (1992) Cornwell, T. J., Uson, J. M., & Haddad, N. 1992, A&A, 258, 583
 Datta et al. (2009) Datta, A., Bhatnagar, S., & Carilli, C. L. 2009, ApJ, 703, 1851
 Datta et al. (2007) Datta, K. K., Choudhury, T. R., & Bharadwaj, S. 2007, MNRAS, 378, 119
 Di Matteo et al. (2002) Di Matteo, T., Perna, R., Abel, T., & Rees, M. J. 2002, ApJ, 564, 576
 Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415
 Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181
 Gleser et al. (2008) Gleser, L., Nusser, A., & Benson, A. J. 2008, MNRAS, 391, 383
 Hales et al. (1988) Hales, S. E. G., Baldwin, J. E., & Warner, P. J. 1988, MNRAS, 234, 919
 Harker et al. (2010) Harker, G., Zaroubi, S., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Ciardi, B., Jelić, V., Koopmans, L. V. E., Labropoulos, P., Mellema, G., Offringa, A., Pandey, V. N., Pawlik, A. H., Schaye, J., Thomas, R. M., & Yatawatta, S. 2010, MNRAS, 612
 Harker et al. (2009a) Harker, G., Zaroubi, S., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Ciardi, B., Jelić, V., Koopmans, L. V. E., Labropoulos, P., Mellema, G., Offringa, A., Pandey, V. N., Schaye, J., Thomas, R. M., & Yatawatta, S. 2009a, MNRAS, 397, 1138
 Harker et al. (2009b) Harker, G. J. A., Zaroubi, S., Thomas, R. M., Jelić, V., Labropoulos, P., Mellema, G., Iliev, I. T., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Ciardi, B., Koopmans, L. V. E., Pandey, V. N., Pawlik, A. H., Schaye, J., & Yatawatta, S. 2009b, MNRAS, 393, 1449
 Hobson & Maisinger (2002) Hobson, M. P. & Maisinger, K. 2002, MNRAS, 334, 569
 Jelić et al. (2008) Jelić, V., Zaroubi, S., Labropoulos, P., Thomas, R. M., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Ciardi, B., Harker, G., Koopmans, L. V. E., Pandey, V. N., Schaye, J., & Yatawatta, S. 2008, MNRAS, 389, 1319
 Komatsu et al. (2010) Komatsu, E., Smith, K. M., Dunkley, J., Bennett, C. L., Gold, B., Hinshaw, G., Jarosik, N., Larson, D., Nolta, M. R., Page, L., Spergel, D. N., Halpern, M., Hill, R. S., Kogut, A., Limon, M., Meyer, S. S., Odegard, N., Tucker, G. S., Weiland, J. L., Wollack, E., & Wright, E. L. 2010, ArXiv eprints
 Labropoulos et al. (2009) Labropoulos, P., Koopmans, L. V. E., Jelic, V., Yatawatta, S., Thomas, R. M., Bernardi, G., Brentjens, M., de Bruyn, G., Ciardi, B., Harker, G., Offringa, A., Pandey, V. N., Schaye, J., & Zaroubi, S. 2009, ArXiv eprints
 Lidz et al. (2008) Lidz, A., Zahn, O., McQuinn, M., Zaldarriaga, M., & Hernquist, L. 2008, ApJ, 680, 962
 Liu et al. (2009) Liu, A., Tegmark, M., & Zaldarriaga, M. 2009, MNRAS, 394, 1575
 Lonsdale et al. (2009) Lonsdale, C. J., Cappallo, R. J., Morales, M. F., Briggs, F. H., Benkevitch, L., Bowman, J. D., Bunton, J. D., Burns, S., Corey, B. E., Desouza, L., Doeleman, S. S., Derome, M., Deshpande, A., Gopala, M. R., Greenhill, L. J., Herne, D. E., Hewitt, J. N., Kamini, P. A., Kasper, J. C., Kincaid, B. B., Kocz, J., Kowald, E., Kratzenberg, E., Kumar, D., Lynch, M. J., Madhavi, S., Matejek, M., Mitchell, D. A., Morgan, E., Oberoi, D., Ord, S., Pathikulangara, J., Prabu, T., Rogers, A., Roshi, A., Salah, J. E., Sault, R. J., Shankar, N. U., Srivani, K. S., Stevens, J., Tingay, S., Vaccarella, A., Waterson, M., Wayth, R. B., Webster, R. L., Whitney, A. R., Williams, A., & Williams, C. 2009, IEEE Proceedings, 97, 1497
 McQuinn et al. (2006) McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815
 Mitchell et al. (2008) Mitchell, D. A., Greenhill, L. J., Wayth, R. B., Sault, R. J., Lonsdale, C. J., Cappallo, R. J., Morales, M. F., & Ord, S. M. 2008, IEEE Journal of Selected Topics in Signal Processing, vol. 2, issue 5, pp. 707717, 2, 707
 Morales et al. (2006a) Morales, M. F., Bowman, J. D., Cappallo, R., Hewitt, J. N., & Lonsdale, C. J. 2006a, New Astronomy Review, 50, 173
 Morales et al. (2006b) Morales, M. F., Bowman, J. D., & Hewitt, J. N. 2006b, ApJ, 648, 767
 Morales & Hewitt (2004) Morales, M. F. & Hewitt, J. 2004, ApJ, 615, 7
 Morales & Wyithe (2009) Morales, M. F. & Wyithe, J. S. B. 2009, ArXiv eprints
 Myers et al. (2003) Myers, S. T., Contaldi, C. R., Bond, J. R., Pen, U., Pogosyan, D., Prunet, S., Sievers, J. L., Mason, B. S., Pearson, T. J., Readhead, A. C. S., & Shepherd, M. C. 2003, ApJ, 591, 575
 Parsons et al. (2009) Parsons, A. R., Backer, D. C., Bradley, R. F., Aguirre, J. E., Benoit, E. E., Carilli, C. L., Foster, G. S., Gugliucci, N. E., Herne, D., Jacobs, D. C., Lynch, M. J., Manley, J. R., Parashare, C. R., Werthimer, D. J., & Wright, M. C. H. 2009, ArXiv eprints
 Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., Wright, M. C. H., Bradley, R. F., Gugliucci, N. E., Parashare, C. R., Benoit, E. E., Aguirre, J. E., Jacobs, D. C., Carilli, C. L., Herne, D., Lynch, M. J., Manley, J. R., & Werthimer, D. J. 2010, AJ, 139, 1468
 Peebles (1993) Peebles, P. J. E. 1993, Principles of physical cosmology
 Pen et al. (2009) Pen, U., Chang, T., Hirata, C. M., Peterson, J. B., Roy, J., Gupta, Y., Odegova, J., & Sigurdson, K. 2009, MNRAS, 399, 181
 Shaver et al. (1999) Shaver, P. A., Windhorst, R. A., Madau, P., & de Bruyn, A. G. 1999, A&A, 345, 380
 Thomas et al. (2009) Thomas, R. M., Zaroubi, S., Ciardi, B., Pawlik, A. H., Labropoulos, P., Jelić, V., Bernardi, G., Brentjens, M. A., de Bruyn, A. G., Harker, G. J. A., Koopmans, L. V. E., Mellema, G., Pandey, V. N., Schaye, J., & Yatawatta, S. 2009, MNRAS, 393, 32
 White et al. (1999) White, M., Carlstrom, J. E., Dragovan, M., & Holzapfel, W. L. 1999, ApJ, 514, 12
 Zaldarriaga et al. (2004) Zaldarriaga, M., Furlanetto, S. R., & Hernquist, L. 2004, ApJ, 608, 622