Bright Source Subtraction Requirements For Redshifted 21 cm Measurements

Bright Source Subtraction Requirements For Redshifted 21 cm Measurements

A. Datta11affiliation: New Mexico Tech, Socorro, NM 87801, USA 22affiliation: National Radio Astronomy Observatory, Socorro, NM 87801, USA , J.D. Bowman 33affiliation: California Institute of Technology, Pasadena, CA 91125, USA 44affiliation: Hubble Fellow and C.L. Carilli 22affiliation: National Radio Astronomy Observatory, Socorro, NM 87801, USA , adatta@nrao.edu
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 sky-model 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 two-dimensional power spectrum.

Subject headings:
Cosmology: Early Universe, Galaxies: Intergalactic Medium, Radio Lines: General, Techniques: Interferometric, Methods: Data Analysis

1. 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 ultra-violet and X-ray radiation of early stars, galaxies, and black holes. Recent observations of the Gunn-Peterson 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 E-mode of polarization which rules out any instantaneous reionization at at 3.5  level. In case of the Gunn-Peterson 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 red-shifted 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 low-frequency 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 low-frequency 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 free-free 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 low-frequency 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 data-set 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 data-set 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 three-dimensional power spectrum, and two-dimensional 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 field-of-view of the total number of sources above 1 Jy is , following the above power-law 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 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 ().

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 3-dimensional 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.

Figure 1.— Simulation of the sky model centered on RA=4h and DEC= as would be observed by the MWA. Clark-CLEAN has been applied to this image using w-projection (256 planes) and natural weighting (Datta et al., 2009).

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 Field-of-view is at 150 MHz.

Table 1Array Specifications
Figure 2.— Array layout for the 512 elements with maximum baseline of 1.5 km.

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 low-frequency 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.

Figure 3.— (a) The spectral profile along two lines of sight in the the final residual image cube following the IMLIN step. In this case, the two pixels were chosen to be next to the positions of sources in the input model. (b) Same as panel (a), but here, the two pixels were chosen to be far from any sources in the input model. Synthesized beam area of 4.5 arcmin 4.5 arcmin is used to convert the flux densities to surface brightness.

2.3. Data Reduction Procedure

The field-of-view 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 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 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 :

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

  • The residual data-sets will be imaged and stored as a cube for the future processing and removal of sources which are below .

Figure 4.— 1D spherically averaged power spectrum of the input Global Sky Model showing the total power of the bright sources in the sky model before any foreground removal has been applied. The thermal noise uncertainty for a 300 hour observation by the MWA is also shown, along with the Hi 21 cm signal power spectrum for a fully neutral IGM () at (Furlanetto et al., 2006).

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 Hour-angle) 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 real-time.

For the residual calibration errors, the model visibilities are given by:

(3)

where are the antenna-dependent 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 real-time 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 real-time processing will be largely independent or highly correlated between individual 8-second 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-second 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-second calibration solutions are perfectly correlated for the duration of one 6-hour 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, 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 8-second solutions for the given antenna throughout the 6-hour 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 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, . 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 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.

(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 three-dimensional power spectrum from the entire data cube, and the two-dimensional power spectrum found by averaging over transverse modes in the full three-dimensional 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 111http://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 222http://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 spherically-averaged 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 two-dimensional. 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

Figure 5.— (a) Angular power spectrum of the UVSUB residual image made after subtraction of a foreground model with source position errors of 0.01, 0.1, and 1 arcsec. (b) Same as panel (a) but for the IMLIN residual image, , that is produced after polynomial fitting and subtraction has been applied along each sight-line. The vertical line in panel (b) corresponds to , below which the polynomial fitting is expected to remove most of the structure. Both panels include the thermal noise uncertainty power spectrum assuming 5000 hours of observation with the MWA and the Hi 21 cm signal power spectrum for a fully neutral IGM (, ; Furlanetto et al. (2006)). These angular power spectra are what would be expected from the MWA if it integrated deep enough to directly image a typical cosmic Stromgren Sphere.
Figure 6.— Same as Figure 5, but for the residuals due to calibration errors. In this case, the three residual angular power spectra are for errors of 0.01, 0.1, and 1.

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 flat-field approximation. Here, is the un-weighted 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 cross-correlating 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 cross-correlation), 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 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 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 Spherically-Averaged Power Spectrum

Figure 7.— (a) 1D spherically-averaged power spectrum of the UVSUB residual image made after subtraction of a foreground model with source position errors of 0.01, 0.1, and 1 arcsec. (b) Same as panel (a) but for the IMLIN residual image, , that is produced after polynomial fitting and subtraction has been applied along each sight-line. The vertical line in panel (b) corresponds to  Mpc, below which the polynomial fitting is expected to remove much of the structure. Both panels include the thermal noise uncertainty spectrum assuming 300 hours of observation with the MWA and binning into logarithmic spherical shells of width , or approximately five bins per decade. The Hi 21 cm signal power spectrum for a fully neutral IGM (, ; Furlanetto et al. (2006)) is also shown. Detecting the spherically-averaged 21 cm power spectrum is the primary goal of the MWA.
Figure 8.— Same as Figure 7, but for the residuals due to calibration errors. In this case, the three residual spherically-averaged power spectra are for errors of 0.01, 0.1, and 1.

The spherically-averaged three-dimensional 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 three-dimensional 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 co-moving 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 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 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 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 spherically-averaged 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 spherically-averaged power spectrum has been binned in logarithmic shells of width , or approximately five bins per decade. As discussed in Lidz et al. (2008), the MWA-512 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 image-cube. 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 spherically-averaged 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 spherically-averaged 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 spherically-averaged 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 spherically-averaged 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 spherically-averaged spectrum is produced with the total bandwidth of 32 MHz.

Figure 9.— (a) Realization of the 2D thermal uncertainty power spectrum after cross-correlating simulated thermal noise maps (300 hours of total integration) from two different epochs (Bowman et al., 2009). (b) 2D power spectrum of the Hi 21 cm signal (Furlanetto et al., 2006) given by , where . Note that the quantity plotted here and the following Figures is in units of mK Mpc. The color scale is shown in .
Figure 10.— (a) 2D power spectrum of the UVSUB residual image made after subtraction of a foreground model with source position errors of  arcsec. (b) Same as panel (a) but for the IMLIN residual image, , that is produced after polynomial fitting and subtraction has been applied along each sight-line. The color scale is shown in .
Figure 11.— Same as Figure 10, but for the residuals due to calibration errors. In this case, the residual 2D power spectrum is shown for a fiducial calibration error level of .

3.3. Two-dimensional Power Spectrum

In the previous section, we showed the analysis of the 1D spherically-averaged 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 redshifted-space distortions have aspherical structure in the Fourier domain. Following McQuinn et al. (2006), we first average over the transverse (angular) direction in the full three-dimension 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 line-of-sight dimension. Next, 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 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 spherically-averaged 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 spherically-averaged power spectrum, it should be noted that the calibration accuracy is determined on the basis of a  hrs of integration by MWA-512. 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 2-D 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 data-set 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 real-time calibration pipeline. The major implication for shorter time-scale 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 6-hour observation night, but uncorrelated between observing nights. If the residual calibration error where instead perfectly random between every 8-second cycle of the real-time 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 real-time 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 wide-field gain calibration of the primary beam and ionosphere.

AD and CC are grateful for support from the Max-Planck Society and the Alexander von Humboldt Foundation through the Max Planck Forshungspreise 2005. JDB is supported by NASA through Hubble Fellowship grant HF-01205.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. The National Radio Astronomy Observatory is operated by Associated Universities, Inc, under cooperative agreement with the National Science Foundation.

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 e-prints
  • 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 e-prints
  • 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. 707-717, 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 e-prints
  • 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 e-prints
  • 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
60278
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description