The Giant Metrewave Radio Telescope (GMRT) reionization effort aims to map out the large-scale structure of the Universe during the epoch of reionization (EoR). Removal of polarized Galactic emission is a difficult part of any 21 cm EoR programme, and we present new upper limits to diffuse polarized foregrounds at 150 MHz. We find no high-significance evidence of polarized emission in our observed field at mid-galactic latitude (J2000 08h26m+26). We find an upper limit on the two-dimensional angular power spectrum of diffuse polarized foregrounds of (l2Cl/2π)1/2 < 3 K in frequency bins of width δν= 1 MHz at 300 < l < 1000. The three-dimensional power spectrum of polarized emission, which is most directly relevant to EoR observations, is [k3PP(k)/2π2]1/2 < 2 K at k⊥ > 0.03 h Mpc−1, k < 0.1 h Mpc−1. This can be compared to the expected EoR signal in total intensity of [k3P(k)/2π2]1/2∼ 10 mK. We find that polarized structure is substantially weaker than suggested by extrapolation from higher frequency observations, so the new low upper limits reported here reduce the anticipated impact of these foregrounds on EoR experiments. We discuss the Faraday beam and depth depolarization models and compare predictions of these models to our data. We report on a new technique for polarization calibration using pulsars, as well as a new technique to remove broad-band radio frequency interference. Our data indicate that, on the edges of the main beam at the GMRT, polarization squint creates ∼3 per cent leakage of unpolarized power into polarized maps at zero rotation measure. Ionospheric rotation was largely stable during these solar minimum nighttime observations.
A current frontier in observational and theoretical cosmology is the search for and understanding of large-scale structure during the epoch of reionization (EoR). Detection of reionization at z∼ 10 would allow study of the emergence of the first abundant luminous objects. Such observations hold the promise of examination of astrophysical and cosmological processes at epochs as early as a few hundred million years after the start of the expansion of the Universe.
The Wilkinson Microwave Anisotropy Probe (WMAP) satellite has measured polarization in the cosmic microwave background (CMB) at large angular scales. This polarization is believed to arise from the Thomson scattering of the CMB photons near the EoR (Komatsu et al. 2009; Nolta et al. 2009). The observed optical depth τ∼ 0.089 ± 0.016 corresponds to an instantaneous reionization redshift of zreion= 10.8 ± 1.4.
One way to study the reionization transition in detail is by imaging redshifted 21 cm emission. At redshifts above the EoR transition, the gas is neutral and is predicted to glow with about 25 mK sky brightness temperature. After reionization is complete, this glow is absent. At redshifts close to the transition, a patchy sky is expected. Simulations (Iliev et al. 2008) suggest that with existing telescopes a measurement near 150 MHz may allow for statistical detection of ∼20 Mpc patchiness in the neutral hydrogen. This detection would pin down the reionization redshift and begin the process of more detailed study of the transition.
The 21 cm EoR signal, in the 10 mK range, is much smaller than the expected ∼10 K patchiness due to Galactic synchrotron emission and extragalactic foregrounds. One of the most difficult challenges of 21 cm EoR astronomy is the removal of these bright foregrounds. The proposed strategy makes use of the very different spatial and spectral character of the two components of sky structure. Since the received frequency of the 21 cm line emission encodes the distance of the source, the EoR signal consists of a series of source screens that are largely independent. For example, if one detects a patchy layer of cosmic material at 150 MHz due to the EoR transition, another patchy layer must be present at 151 MHz, but these two images come from distinct radial shells and the bright and dim patches in the two images are not expected to be aligned in the sky. The characteristic size of the ionization structures near the epoch of ionized-cell overlap is approximately 10 h−1 Mpc, which corresponds to 5 arcmin in angular size, or about 1 MHz in radial distance for the redshifted 21cm line at z∼ 9. In contrast, the foreground emission is highly correlated from one frequency to the next, since most radio-bright structures emit with a very smooth synchrotron spectrum. To remove the foregrounds, some type of frequency-differencing technique will be needed.
At frequencies above 300 MHz, many radio sources with a synchrotron spectrum have been shown to be linearly polarized, and this raises a concern for the frequency-differencing scheme. In propagation through the interstellar medium, the polarization rotates in the sky because of the Faraday effect. In the polarization basis of the telescope, this means that each linear polarization component can oscillate with frequency. If the oscillation period falls close to the frequency separation used for the foreground subtraction, incomplete removal of Galactic synchrotron emission may leave a residual that masks the EoR signal.
One way to avoid the Faraday rotation problem is to measure total intensity, described by the Stokes I parameter, using an instrument that rejects the linearly polarized Stokes components Q and U. However, all real instruments have various sources of leakage between these components. To design these instruments and plan future EoR observations, one needs information on the polarized sky brightness in the EoR band near 150 MHz.
Recent studies at 350 MHz (de Bruyn et al. 2006) report that the polarized component has much more structure on arcminute angular scales than the total intensity. Polarized structure amounting to several kelvins has been reported. Scaling using a synchrotron spectral index of 2.6 (de Oliveira-Costa et al. 2008), the polarized structure could be tens of kelvins in the EoR band. This means that instrumental leakage from Q/U to I could severely limit the ability to remove Galactic emission.
This paper places new constraints on polarized sky brightness in the EoR range of frequencies near 150 MHz. We find that the level of polarized sky brightness is well below that expected by extrapolation from higher frequencies. We discuss possible mechanisms which would lead to this strong depolarization, and suggest future techniques to differentiate between them.
2 The GMRT EoR project
Our group has initiated an effort to search for 21 cm structures at the EoR using the Giant Metrewave Radio Telescope (GMRT) in India.
The telescope consists of 30 antennas of 45 m diameter. 14 are designated as ‘central core’ antennas and lie within a 1 km area. The remaining 16 are along three arms of length up to 10 km, designated as east, west and south. The dense layout of the core allows high brightness sensitivity, which is needed for the search for reionization. For this experiment, the 150 MHz feeds are used. These consist of orthogonal pairs of folded dipoles, backed by a ground plane. These antennas couple to X and Y linear polarizations, but the two signals pass through a hybrid coupler before entering the amplifiers. For each antenna, this results in a pair of right and left circularly polarized signals which are amplified, upconverted and transmitted optically to the receiver room. Later processing allows measurement of the full set of Stokes parameters.
To attempt the EoR experiment, and for other work at the GMRT, we built a new signal processing system for the telescope, the GMRT software backend (Roy et al., in preparation). This consists of 16 commercial analogue to digital (AD) sampling boards installed in an array of off-the-shelf computers. The AD boards have four input channels, and are connected to a common clock and trigger signal, allowing synchronous sampling. The sampling system is capable of 32 MHz bandwidth, but for this experiment only 16 MHz was used. The passband was defined by an Intermediate Frequency (IF) filter. All 60 signals were sampled with 8 bit precision at 33 MSample/s. Each AD board transfers the digitized data into its individual host computer. These streams are Fourier transformed in blocks of length 4096 samples at 16 bit precision. The 2048 complex Fourier coefficients are then rescaled to 4 bits, and sent over a gigabit network to correlation nodes. Each block of 128 frequencies is sent to one of the 16 software correlation nodes. These products, which we call visibilities, are accumulated for 1/4 s in 16 distinct gates. Thus, the initial visibilities have a frequency and time resolution of 7.8 kHz and 1/4 s (without accounting for the gates), respectively.
The ‘gates’ are essential to our polarization calibration technique. These gates in time are synchronized so each covers one of the 16 segments of the pulsation cycle of a pulsar in the field. The pulsar is a known source of polarized emission. By comparing pulsar-on to pulsar-off, we can measure the system gain directly using a sky source. The raw 1/4 s averaged visibilities were stored on disk, as well as 16 fringe-stopped gated visibilities integrated for 16 s, and averaged over frequency into 128 frequency channels.
All these signal processing calculations occur simultaneously, and are structured as individual asynchronous pipeline processes. The processor and network speeds are sufficient that each calculation is completed in real time, and there is less than 10 per cent data loss.
Our target field was centred on pulsar PSR B0823+26 (J2000: 08 26 51.38+26 37 23.79). The field was chosen to contain a bright pulsar, and minimal other bright sources. The galactic latitude is b= 30°, and happens to be a relatively cold spot in the sky.
The brightest radio source within the primary beam, and our primary flux reference, is 5C 7.245, located 12 arcmin from the pulsar. Its flux is estimated to be 1.5 Jy at 151 MHz from the 7C Lowdec Survey (Waldram et al. 1996). It is a radio galaxy at z= 1.61 (Willott, Rawlings & Blundell 2001). The source has two components separated by 16 arcsec, and is resolved by the longest baselines at the GMRT.
A bright source outside the main beam, 3C200, is 5 to the north. It is attenuated by the primary beam profile, but shows up as the strongest source in the raw polarized maps due to I to P leakage in the side lobes.
In Table 1, we list the 10 brightest sources in the field after the primary beam power attenuation effect from the 7C Lowdec Survey. Also listed are their fluxes from the VLA (Very Large Array) Low-frequency Sky Survey (VLSS) at 74 MHz (Cohen et al. 2007) and the Texas Survey of radio sources at 365 MHz (Douglas et al. 1996). Based on these three surveys, we calculate the spectral indices and the errors, as listed in Table 1.
|Source||RA (°)||Dec. (°)||F7C-BEAM (Jy)||F7C (Jy)||F74 (Jy)||F365 (Jy)||α||Δα|
|Source||RA (°)||Dec. (°)||F7C-BEAM (Jy)||F7C (Jy)||F74 (Jy)||F365 (Jy)||α||Δα|
This region was mapped at 150 MHz (Landecker & Wielebinski 1970) and found to have brightness temperature 170K. de Oliveira-Costa et al. (2008) have reconstructed a variety of sky maps and used them to estimate spectral indices. From these reconstructions, we estimate a spectral index of 2.6 for this region, so our sky brightness temperature varies from 200 to 150 K over our spectral band (140–156 MHz).
Data were taken from 2007 December 7–18 starting at 10 pm local time, going until 7 am. This period has minimal radio frequency interference (RFI). The data were recorded in 1 hour ‘scans’. The analysis in this paper is primarily based on the data taken during the night of December 9, which was fully reduced through the described pipeline at the time of writing of this paper.
4 INTERFERENCE REMOVAL
Terrestrial RFI was removed in two stages. First, spectral-line RFI was flagged. In each frequency subset of 128 channels, the distribution of intensities is calculated. The upper and lower quartile boundaries are used to estimate the amplitude of the noise. Any outliers further than 3σ are masked. Each mask is 8 kHz wide and 4 s long. Approximately 1 per cent of the data are flagged by this process for deletion.
A bigger problem at the GMRT has been broad-band RFI. This is particularly problematic for short baselines. The rationale for our approach to broad-band RFI segregation is that sources in the sky produce visibilities that oscillate as the Earth turns. Meanwhile, sources on the ground, although their intensity typically varies with time, remain at fixed lag for each baseline. We used a singular value decomposition (SVD) to sort among the terrestrial and celestial classes.
We assemble the visibilities in a two-dimensional (2D) rectangular matrix, Vi(ν, t), which is a function of baseline i, frequency ν and time t. We find that Earth-fixed broad-band sources tend to flicker synchronously in all baselines at all frequencies, and are therefore factorizable as VRFI=Li(ν) T(t). We call Li(ν) the visibility template and Ti(t) the temporal template. The RFI visibility can thus be written as a sum of individual RFI sources α,
In contrast, objects on the celestial sphere exhibit fringe rotation at a unique rate for each baseline, and thus do not factor. Celestial sources produce very small eigenvalues under an SVD.
A decomposition of the form (1) can be accomplished via an SVD algorithm as long as the eigenvectors are orthogonal. A pair of visibility vectors due to two physically separate sources will tend to have little overlap, since the dimensionality of this space is huge. There are of the order of 302= 900 baselines between antenna pairs, each with four polarization combinations. Two distinct sources may accidentally fall at the same lag for one pair but will be separated from each other when other pairings are considered. The time templates might not be as orthogonal. In particular, electric arcing from AC transmission lines, which is likely the cause of a substantial subset of the broad-band RFI, might occur more often at the AC waveform peak. This timing can be common to many sources. For this work, only about 100 of the largest SVD eigenvalues are flagged as RFI, in a matrix with about 109 entries, so even if the eigenvectors are not perfectly orthogonal, the very sparse space spanned by the first 100 eigenmodes is still given by an SVD decomposition.
To carry out the procedure, we begin by organizing the data in a matrix. Each hour scan has 14 400 time records, which we call the rows of our matrix, the column number 2048 × 60 × 61 entries which correspond to the number of frequency channels and baselines. Each matrix entry is a complex visibility. This 7495 680 × 14 400 matrix is then SVD decomposed. The first 100 right eigenvectors are tagged as ‘RFI’, and removed. This removes 0.7 per cent of the matrix entries.
We have also used the SVD procedure to physically locate RFI sources. We selected several of the brightest eigenmodes for further study. Near-field imaging using these eigenvectors has allowed us to localize the sources on the Earth's surface with 100 m precision. Visiting these locations with radio-direction-finding equipment, we have localized candidate emitters in several cases. The brightest two appear to originate from extraneous thin wires hanging vertically from high-tension power lines. Other sources appear related to intermittent radiation from transformers. So far, we have visited only a few sources.
The brightest SVD eigenmodes do appear to represent RFI, not celestial sources, however there are some limitations to the technique. The antennas slowly turn as they track objects in the sky, and this rotation leads to phase changes for terrestrial RFI sources. The telescope feed moves by half a wavelength every 12 min. For an RFI source located between two antennas, the relative delay changes by a full wavelength. We believe that this sometimes causes individual sources to break up, occupying several RFI modes. However, as long as the source is bright, it is still identifiable, even if split into components.
Another limitation: the SVD analysis tends to misidentify, as RFI, celestial sources with fringes that rotate slowly. A baseline that is oriented perfectly north–south does not produce any fringe rotation, so baselines with small v are more prone to erroneous projection of celestial sources on to the RFI modes. We find that setting the RFI cut at 100 eigenmodes has a notable impact only on baselines with |v| < 10λ. Increasing the number of cut eigenmodes widens the impacted v range, while substantially decreasing the number of cut eigenmodes leaves significant residual RFI power.
Empirically, a cut at 100 eigenmodes works very well, most visible broad-band RFI had been removed from the data after this process.
Polarization calibration is particularly challenging at low frequencies. The primary beam is wide, telescope gain is low, interference is frequent, the ionosphere is variable and the sky offers few polarized point sources. Normally, an iterative procedure is used, where one starts with a guess for the sky model, determines the system gain and phases, makes a new map and uses this as a model. Errors in the sky model are entangled with errors in the telescope system calibration.
We used a novel approach by observing a field containing a bright pulsar. By subtracting the off-phase from the on-phase, one can lock into the modulated emission. Of the 16 gates, we took the gate containing the pulsar-on state, and subtracted the average of the two flanking gates. After this subtraction, the measured visibilities contain only the pulsar. Other celestial or terrestrial signals remain only if they happen to be modulated in sync with the pulsar. While the pulsar flux may vary from pulse to pulse, each antenna sees the same flux, and the pulsar serves as a common source allowing the relative gain of the 30 antennas to be measured. The pulsar also serves as a guide star. Pulsar fluxes can be variable in time, but their positions are very stable.
This pulsar is strongly polarized for individual pulses but the polarization angle varies from pulse to pulse. The polarization fraction for individual pulses has been measured to be 40 per cent (Hankins & Rankin 2008) at 430 MHz. Many pulsars have polarization fractions that decrease with frequency so the fraction may be higher at 150 MHz. Our data were stored as 16 s averages, which corresponds to 30 pulses. We found that the average polarization fraction for these averages was 20 per cent, with an rms scatter of 8°. The polarization fraction decreases as longer averages are assembled. This is consistent with a highly polarized pulse whose angle varies from pulse to pulse, with only a slight alignment preference. This random angle behaviour is also observed at 430 MHz. Fig. 1 shows the 9 min averaged angle, after correcting for parallactic angle rotation.
We also used the GMRT noise injection system to measure the gain of the electronic signal chains. At each prime focus, a noise waveform is injected into the two polarization channels, modulated at 0.5 Hz. We used the noise source on a setting which injects 40 K, common between the L and R channels. We measured the strength of the noise source in the L–R correlation of each antenna and used this as a gain reference.
The goal of the calibration process is to determine two complex gains for each of the 60 signals. One is the response to its own polarity, and the other is the leakage from the orthogonally polarized mode. The pulsar model appears to work well, with typical rms phase closure errors around 1° averaged over 20 min. This model reduces the 1830 complex visibilities into 120 complex gains for each frequency channel. The visibilities are reproduced from this reduced degree of freedom model. This is called the ‘bandpass calibration’. In addition, the model allows for 1792 complex gains, one for each time interval of approximately 16 s, corresponding to 8 hours of integration. The model assumes that the gain factor as a function of time times a function of frequency. The 1792 × 1830 numbers at each frequency are reduced to 1792+120. Polarization calibration was achieved by using the parallactic angle rotation to decompose the pulsar into three components: the ‘unpolarized’ source which is invariant under parallactic angle rotation, and two rotating components. More details can be found in Pen et al. (in preparation).
The system is fully calibrated using only the pulsar and noise sources, except for an overall scale described below. We find that the unpolarized pulsar flux varies by about 2 per cent when averaged over 1 hour, relative to the average of all noise injection sources. The noise source and system gains appear stable for most antennas, with a few exceptions. The pulsar flux can be readily subtracted from the time-averaged data.
This calibration procedure results in a polarization map referenced to the pulsar polarization angle at each frequency. We then correct this with the pulsar's known RM = 5.9 rad m−2 (Manchester 1974; Manchester et al. 2005) for the rest of this paper.
6 FLUX SCALE
We use the pulsar for relative antenna calibration, but not to establish the overall gain, since the pulsar flux varies with time. To accomplish the overall flux calibration, we use four well-measured bright continuum radio sources in the field. We estimate the flux of the four brightest sources that lie within the full width half-maximum (FWHM) of the primary beam using the estimated spectral indices of these sources from the catalogues listed in Table 1 and correcting for the primary beam power attenuation. We calculate the source flux values from our data at 151 MHz and match the average flux to that of the 7C fluxes. The estimated flux and respective errors are listed in Table 2. The resulting fractional errors on these four source fluxes are less than 13 per cent, and the error on the flux scale in the calibration procedure is 6 per cent.
|Source||RA (°)||Dec. (°)||F7C (Jy)||F151 (Jy)||Error|
|Source||RA (°)||Dec. (°)||F7C (Jy)||F151 (Jy)||Error|
Note.F7C is again the flux from the 7C catalogue at 151 MHz, F151 is the flux estimated from our data at 151 MHz and error is the respective fractional error on the estimated flux.
The calibration of the power spectrum amplitude depends on the point-source flux calibration described above and also on the primary beam width and profile. Throughout, we approximate the primary beam as a round Gaussian beam with a FWHM of 3. This profile is a good fit to our holographically measured primary beam power out to about 2° radius.
After the primary calibration, the dominant polarized source is the pulsar. The second brightest apparent polarized source is 3C 200, 5 to the north, with an intrinsic flux of 13 Jy at 150 MHz. The source is probably not significantly intrinsically polarized, but rather because the side lobes are strongly polarized, it leaks into the polarized maps.
Fig. 2 shows the apparent luminosities and polarizations of the brightest point sources in the field. Objects appear more polarized as they are located away from the image centre, suggesting the primary beam as the origin of this apparent polarization. Fig. 3 shows the radial apparent polarization.
Precise cleaning of the polarized image shown in Fig. 2 will require detailed knowledge of the off-axis leakage at the GMRT, which we have not yet measured. Despite this imperfect knowledge, we have been able to make some progress using a modified version of a standard ‘clean’ procedure. We fringe stop the data iteratively to sources in descending order of flux as predicted from the 7C catalogue convolved by the Gaussian primary beam model. At each source position, we measure the equivalent point-source flux by summing all baselines. The point source is modelled as two complex numbers, one for the intensity and one for polarization. Separate values are fit every 4 MHz and every 20 min. A point source with this flux is then subtracted from the data. We repeat this for the 31 brightest sources in the field. For the short baselines, sources are confused, so we downweighted the central square antennas by a factor of 10. A baseline containing one central square antenna gets 10x less weight, and a baseline with two such antennas is downweighted by 100x. Similarly, we downweight the last two antennas on each arm, since some of the sources are resolved by them. This effectively uses the intermediate baselines to solve the model, and uses the cross-correlations from arm antenna to central square to extrapolate the flux on the shortest and longest baselines.
Each of the 31 point sources is fitted with 80 parameters: four frequencies by 20 temporal fluxes (one per 20 min, total of 400 min), for a total of 248 degrees of freedom. To quantify the amount of spurious signal removal, we added a 30 mJy trace in the form of a Gaussian footprint centred at u= 100, v=−100 with an rms of σ= 20 in the Fourier space. Applying the same clean on this data and subtracting the original clean data recover the original footprint to better than 99 per cent, while leakage into originally empty baselines acquires up to 6 per cent of the peak flux. We conclude that this process has at most a small impact on the flux at |u| ∼ 140, i.e. l∼ 900.
8 ROTATION MEASURE SYNTHESIS
At low frequencies, the Faraday rotation can have a large effect. The pulsar has a rotation measure of RM = 5.9 rad m−2. The apparent polarization angle rotates as χ= RMλ2+constant. We use units of rad m−2 throughout this paper. At 2 m wavelength, the angle changes 5.3 rad over the 16 MHz observed band. Since polarization angles are periodic over 180°, this results in almost two polarization rotation periods. After pulsar referenced polarization calibration, all angles are defined relative to the pulsar. A direct 2D map would strongly suppress sources whose RM is more than ΔRM > 3.5 away from the pulsar's, since the rotating polarization angle causes strong cancellation of the source flux.
To make maps at various Faraday depths, we apply a method known as rotation measure synthesis (Brentjens & de Bruyn 2006). The frequency channels are mapped into a three-dimensional (3D) grid of u−v−λ2. We chose a 2048 × 2048 × 128 grid, where all our frequency channels are scaled to fit in the first 64λ2 channels. Padding with 64 zeros increases the sampling in the Faraday depth and results in several reconstructed slices per Rotation Measure Transfer Function (RMTF) width (although note that this does not affect the width of the RMTF). A 3D FFT then makes a 3D map with axes (x, y, φ) The frequencies are binned to the nearest grid point, resulting in a 2/π= 64 per cent bandwidth depolarization at the Nyquist limit of φ= 56 rad m−2.
At the Faraday depth φ= 0, we can clearly see the beam squint: point sources in the side lobes appear radially polarized. Our calibration procedure equalizes the polarized response at the pulsar position, but the leakage will vary across the beam in a quadrupolar pattern, producing the radial ellipses in Fig. 2. This squint, a difference between the width of the E- and H-plane antenna patterns, can arise from a variety of reasons. For example, the prime focus instruments are mounted on four support legs, which are each narrower than our observing wavelength (2 m). They carry reflecting currents more efficiently along the legs than transverse. The illumination pattern in one linear polarization has a horizontal line blocked, while the orthogonal polarization pattern would have the opposite line blocked. To estimate the size of such a squint term, we assume that wavelength of the order of 1 is removed from a 45 m dish. The effective area of the dish is about 1000 m2, of which 10 per cent is blocked, and this blockage is different for the two polarizations. This could make the E and H beams in the sky different. The blockage attenuates the main lobe of the beam, while generating an extended pedestal in the pattern. In such a model, one expects unpolarized objects within the first null of the primary beam to appear radially polarized (in terms of the electric field). Beyond the null, objects would be tangentially polarized. Another source of beam squint occurs because the illumination of the dish by the feed is not exactly round.
In principle, one could measure the primary beam for both polarizations and correct the leaked maps to produce a clean Stokes I map. This would require a non-local calculation, both in the sky and in the Fourier space, and would require the use of full matrix inversion. Codes like cbi-gridr (Myers et al. 2003) could be adapted for this purpose.
In order to minimize the impact of non-coplanar w terms, we only used baselines with |w| < 200. We convolved each visibility with a Gaussian of rms σ= 23 in the u–v space to reduce the field of view to the central square degree, which further reduces the impact of non-coplanarity, and the beam squint. The phase errors from non-coplanar w terms are Δφ∼wΔθ2/2, and we expect phase errors up to a few degrees at the edge of the reduced field.
A given point in u–v space is typically covered by several baselines, at different frequencies. The RM synthesis maps are a sum of these points, weighted by the phase factor that accounts for the oscillation due to the chosen rotation measure. Power spectra can be used to summarize the statistics of these maps (Appendix A). There are several types of power spectra that can be computed, both 2D and 3D, containing different types of information. Perhaps the most straightforward is the square of each gridded visibility in a frequency interval ν,
Since different baselines see the sky at different parallactic angles, one can also measure the cross-correlation between distinct visibilities. To reduce the impact of leakage, RFI and thermal noise, we can discard the self-terms in the power spectrum. In equation (3), we note that the visibility location u=u(b, t) is a function of the label of baseline b and the time of observation t. Each gridded visibility can contain contributions from different baselines taken at different times. Instead of summing and squaring, we can square and drop all contributions where b=b′.
Computationally, it is expensive to store all baselines contributing to a u–v grid point, so we randomly separate baselines into two bins, grid the visibilities from each bin and compute their cross-correlation. We repeat this many times to include all possible pairs of baselines contributing to a product, in this particular case 129. This results in an average power 〈Cl〉, as well as a standard deviation δCl. The latter is a rough estimate of the standard deviation on 〈Cl〉. The variance on the mean is smaller than the variance on each subsample. If the baselines were densely distributed, the variance is overestimated by a factor of 2. The actual baselines are not dense, so sometimes one of the elements of the pair will be empty. This causes the variance to be further overestimated. The average 〈Cl〉 appears consistent with zero, so we roughly interpret δCl as an effective 2σ upper bound, shown in Fig. 4.
If one fixes a rotation measure, all frequencies can be summed into a single u–v grid point, which further improves sensitivities, and localizes the RM of the power spectrum. Before subtracting the point sources, excess power is visible at φ= 0 when imaging a 4° field. After point-source subtraction, we saw no rotation measure at which there is a statistically significant signal. Restricting the field to 1° using a primary beam convolution in u–v space should further suppress any potential contamination. We apply θb= 1° in equation (3). We show two representative lines in Fig. 4 at φ=−14 and 14 rad m−2, where we have plotted the absolute value of the power. Because we discarded autocorrelation terms, noise can lead to positive or negative correlations. The effective width of the window in the Faraday measure space is determined by the 16 MHz bandwidth, which corresponds to bins of width Δφ∼ 3.5 rad m−2.
[Note: there are several ways to measure the width of the RMTF. For a top-hat range of frequencies, |RMTF|(φ) = |sinc(φΔλ2)| (the phase depends on the wavelength to which one de-rotates). The width relevant for power spectrum analysis (Appendix A) is since if a slice is suppressed by the RMTF then its power spectrum is suppressed by rmtf2. One could alternatively compute the FWHM of |RMTF|, which for the sinc function is 3.79/Δ(λ2).]
The observed apparent polarized power could either be polarized structures in the sky, or residual leakage from beam squint. Towards the edge of the beam, point sources appear in the polarized maps at up to 10 per cent of their unpolarized flux, with a quadrupolar pattern characteristic of beam squint: objects on two diametric sides of the fields appear as positive sources, and at 90° appear negative.
We show a map of the point-source-subtracted sky in Stokes Q and U in Fig. 5 at φ= 3 rad m−2, half way to the pulsar RM. Excess power is still discernible in the region of the primary beam. The point-source subtraction reduces the impact of this leakage, but it could still leave residuals. Our apparent measured polarized signal might not be truly polarized in the sky, but rather leakage from the unpolarized sky brightness or residual man-made interference. We treat the measured power in this image as an upper limit to the polarized sky brightness.
9.1 3D power spectrum
EoR structure is expected to be 3D, so we need to compare the 3D power spectrum of the polarized foreground to that of the expected signal. Theorists have traditionally written the EoR power spectrum in terms of the temperature variance per logarithmic range in k (Bharadwaj & Ali 2004; Barkana & Loeb 2005):.
It is appropriate to evaluate EoR foregrounds by projecting the Galactic polarization structure on to the cosmological volume using the usual mapping ν→z= 1420 MHz/ν− 1, and estimating the polarized 3D power spectrum that one would infer. Details of this power spectrum definition are given in Appendix A. At our central frequency ν= 148 MHz, and with the WMAP determination of Ωm= 0.26 (Dunkley et al. 2009) the constants described in the appendix are: the redshift z= 8.59; comoving radial distance D= 6720 h−1 Mpc (where H0= 100 h km s−1 Mpc−1); H/c= 5.04 × 10−3h Mpc−1 and C= 12.9 h−1 Mpc MHz−1. The mapping of equation (A38) between the Faraday depth and the apparent radial wavenumber is k∥= 0.0086 h Mpc−1(rad m−2)−1 φ. The 3D polarized power spectrum PP(k, μ) is simply related to the 2D power spectrum CPl of the Faraday depth slice in the RM synthesis map at this value of φ. One can then obtain
The 3D power spectra Δ2P(k) are shown in Fig. 6; the power per ln k is at most ∼2 K at k≤ 0.1h Mpc−1. This rotation measure corresponds to polarized emission at a line-of-sight wavenumber of k∥= 0.12 h Mpc−1. The implied 3D power spectra are shown in Fig. 6.
If extended polarized emission is bright at 150 MHz, leakage from Stokes Q or U into Stokes I may make removal of this foreground from 21 cm EoR images quite difficult. Removal of the synchrotron intensity foreground from EoR images requires projecting out structure that varies slowly with frequency (i.e. has k∥≈ 0). This also removes any component of the EoR signal that varies slowly along the line of sight, i.e. it suppresses modes with |k∥| less than a few times (CΔν)−1. The anticipated EoR structure is at scales of the order of k∼ 0.1 h Mpc−1; this range of scales that will be contaminated by Q, U→I leakage of polarized structure at the Faraday depth of the order of |φ| ∼ 10 rad m−2 (see equation A38). This is the Faraday depth range for which polarization leakage causes the most difficulty for an EoR search, although if the polarization leakage varies rapidly with frequency other values of φ are also of concern (see Appendix A4).
Previous work at 350 MHz (de Bruyn et al. 2006) indicated polarized emission with spatial variations of several kelvins on angular scales of tens of arcminutes. This polarized structure is reported to be brighter than the unpolarized component at the same angular scale. Extrapolating this polarized structure to 150 MHz with a spectral index of 2.6, this corresponds to several tens of kelvins in the EoR band. To remove this foreground would require that instrumental Q or U to I leakage in EoR observations be controlled to ≪10−3.
Our observed upper limit to polarized sky structure in frequency channels of width 1 MHz is in the range of 1–3 K at 300 < l < 1000, with weaker upper limits at smaller scales. This is shown in Fig. 4. For the field we studied, the polarized sky structure is an order of magnitude smaller than one would have predicted by scaling from the published 350 MHz observations.
The results reported in de Bruyn et al. (2006) are at different locations in the sky, so a direct comparison is not straightforward. The (de Bruyn et al. 2006) field is centred at 13h14m+45 (J2000), while our field is at 08h26m+26. However, one way to attempt a comparison is to look at the WMAP measurement of polarization emission on the two fields (Hinshaw et al. 2009).
Smoothing the WMAP data to 7 resolution, the 13 hour field has a polarized temperature of 13 μK, while our 08h field is 17 μK. The noise is comparable to the difference, so we estimate that our field has a similar net polarization at high frequency. In the Haslam 408 MHz intensity maps (Haslam et al. 1982), the 13 hour field has a temperature of 20 K which is similar to the 08h field at 18 K. Using the available data, the fields seem quite similar.
de Bruyn et al. (2006) explain the observed polarization structures as due to a magnetic field that is smooth across their field, combined with an electron density that varies. This causes structure in the Faraday depth allowing the polarized sky structure to exceed the total intensity variations. We see no reason to think our field would be different in this regard.
Vinyaikin et al. (1996) provide a model of diffuse Galactic polarization towards the North Celestial Pole. They find polarization of 3.5 ± 1.0 K at 88 MHz (with a 24° FWHM beam) and 2.15 ± 0.25 K at 200 MHz (8° FWHM).1 Their result is difficult to compare to ours, for several reasons. First, there is no overlap in (u, v) coverage between their single-dish (or phased antenna array) measurements and our interferometer results (which have no zero-length baselines). Secondly, one must interpolate between frequencies to compare to our 150 MHz results; Vinyaikin et al. (1996) can fit physical models with polarization nulls between 88 and 200 MHz, but the number and location of such nulls are not uniquely determined. Finally, the observations are in different parts of the sky, and one would not expect geometry-dependent properties such as the frequencies of polarization nulls to be exactly the same.
10.1 Depolarization models
Several depolarizing effects have been predicted at low frequency (Burn 1966; Sokoloff et al. 1998; Vinyajkin 2004). These include beam depolarization due to field gradients, or sub-beam structure, and depth depolarization due to colocation of the rotating and emitting regions. We now consider each effect in turn as possible explanations for the low level of polarization at 150 MHz.
First, we consider beam depolarization (Burn 1966). de Bruyn et al. (2006) reported RM gradients of ∇RM = 2–4 rad m−2 deg−1. At 150 MHz (λ= 2 m), a gradient of this size would result in a periodic modulation of the polarization vector field at
Alternatively, beam depolarization could arise if there was a small-scale random magnetic field. An equipartition magnetic field in the galaxy could be several μG, but this is a factor of several larger than the coherent field. The coherent field is observed to contribute an RM = 5.9 rad m−2 to the pulsar, which results in a 23 rad rotation of the polarization vector. The dispersion measure DM = 19.4751 cm−2 pc (Manchester 1974; Manchester et al. 2005), results in a mean electron-weighted magnetic field along the line-of-sight field of 〈B∥〉= 1.232 RM/DM = 0.37 μG. This coherent part of the field is lower than the equipartition value making plausible the idea that there is fine structure in the field.
Substantial beam depolarization will occur if the small-scale field structure produces variation of rotation of the order of unity across the synthesized beam. Given the current sensitivity at the GMRT, structures at l≫ 104 are unconstrained. A statistically isotropic field on scales smaller than that would have a correlation length 104 shorter than the line of sight. Random rotation measures add stochastically, so the net angle would be equivalent to a 100 times weaker, but coherent field. So a random magnetic field which has five times the field strength of the coherent field and a correlation length 104 times shorter than the rotation screen would lead to beam depolarization. Taking an electron scaleheight of 1 kpc (Cordes & Lazio 2002), which at b= 30° is 2 kpc in projection, would correspond to a fine-structure scale of ∼0.2 pc. A magnetic field of 3 μG on these scales would suffice to depolarize, which is in the plausible range of galactic magnetic fields. The spatial scales are a bit smaller than the outer scales of 1 pc proposed in turbulence models (Haverkorn et al. 2008), but may still be plausible given the general nature of these arguments. At our galactic latitudes, the outer scale may also be substantially larger.
We now take the analysis of fine-scale beam depolarization one step further by assuming a turbulence power spectrum. We consider a Faraday thin emitting region with a rotating medium in the foreground. We model the RM with a Kolmogorov-like power spectrum on small angular scales, i.e. ΔRM2(θ) ∝θ5/3 (Haverkorn et al. 2008). In such a model, the polarization direction is coherent over patches of size θc where ΔRM(θc)λ2∼ 1. For the Kolmogorov spectrum, this implies that θc∝λ−12/5. The peak of the observed polarized power spectrum then scales as lpeak∼θ−1c∝λ12/5. Observations at 350 MHz (de Bruyn et al. 2006) typically find polarized emission at scales of the order of 10 arcmin. There do not appear to be published power spectra of these polarization maps, but if the Q and U Stokes parameters change sign on this scale, as one would expect from the abundant depolarization canals (Haverkorn, Katgert & de Bruyn 2004), then these structures contribute mainly at l∼ 103. Scaling to 150 MHz would move this power up to smaller scales by a factor of (350/150)12/5∼ 8, or lpeak∼ 104. This reduces the tension with our observations, which have upper limits on (l2Cl/2π)1/2 of a few tens of K on these scales (see Fig. 4). This is close to the expected signal from scaling from higher frequencies; the additional presence of depth depolarization or different turbulence parameters in our field would eliminate any tension.
Next we consider depth depolarization, which was already discussed by Burn (1966). This occurs because the emission and rotating media are mixed, and some rotation occurs within the source. This effect has been reported at 1.4 GHz for observations looking right into the Galactic plane, where rotation measures can exceed 1000 rad m−2 (Uyaniker et al. 2003). Our regime has a similar product of RMλ2, so it is tempting to draw analogies. For this part of the discussion, we assume a uniform magnetic field. The synchrotron emission arises from relativistic electrons, which make up a small fraction of the total electron count. However, the Faraday rotation is caused by all free electrons.
The strength of the depolarization depends on the profile of the emitting region. An example of a rapidly depolarizing model is a Gaussian emission region embedded in a medium with uniform ne and B. This might be a reasonable model of an emission region due to cosmic rays diffusing across magnetic field lines. Some prominent structures in the sky, such as the Cygnus loop, may be examples of such sources. This results in a Faraday dispersion function that is a Gaussian
However, it is unlikely that the electron density follows the precise Gaussian profile needed to achieve the deep cancellation of the polarization from different regions. We now consider an alternative scenario, which would lead to much less depth depolarization. Assume that the emission spatial structure traces the ionized gas structure perfectly. Two scenarios could lead to depolarization: (i) the ionized and relativistic populations are uniformly distributed, but the magnetic field is confined to flux tubes; or (ii) the ionized and relativistic components trace each other (e.g. if they were carried by the same magnetic field). In both cases, the Faraday dispersion function is a top-hat, i.e. F(φ) =P0/|φ0| between 0 and φ0. This leads to a net polarization
The emission region could also be only partially overlapping with the rotation region. Many complex permutations are plausible. The simple models presented above are meant to illustrate that strong depth depolarization is possible, given the known parameters.
To further study the nature of the depolarization, multiwavelength polarized observations of the same field are needed. Repeated observations during daytime would give ionospheric rotation measure variations that could separate beam squint from intrinsic polarization. Much more sensitive polarized maps then could be made. Such observation would likely concentrate on fields that show substantial polarized flux.
Several steps could also be implemented in the future. Just as parallactic angle rotation allows us to decompose the polarization components of the pulsar, one could use time variable ionospheric rotation measure to unleak the sources. The ionospheric Faraday rotation changes through the day–night cycle, as we see in Fig. 1. After a few months, the same point in the sky will be viewed through a substantially thicker ionosphere, and the Faraday depth of all sources will be modulated by an ionospheric contribution Δφ. The apparent angle of polarization will change, while unpolarized light stays unchanged. This allows a separation of polarized and unpolarized components which does not require modelling of beam squint. Specifically, in circularly polarized baselines, a two-epoch observation at 1,2 measures
To test the scenario of small-angle beam depolarization, one could compare the rotation measure of pulsars in globular clusters. A search of the Australia Telescope National Facility (ATNF) pulsar data base (Manchester et al. 2005) did not turn up any close pairs. M15 and 47 Tuc have multiple pulsars separated by arcminutes which would be predicted to have fractional variations in RM of the order of 10 per cent in a beam depolarization scenario.
As far as the search for high-redshift reionization is concerned, one wants fields with low Galactic polarized flux. The upper limit to polarization we report indicates that suitable fields are available.
Other steps to reduce the polarization leakage could include mosaicking a wide field, where the angular patterns would cancel, as was done in the work of de Bruyn et al. (2006). In addition, a raster of observation of the pulsar would allow measurement of the polarization calibration into the skirts of the main beam.
Given the low level of polarized sky structure we report, along with the estimates of polarization leakage at the GMRT presented in Fig. 3, we can estimate that even a modest polarization calibration will allow measurement of EoR sky structure. The 1σ upper limit to sky structure measured as a 3D power spectrum [k3PP(k)/2π2]1/2 at scales relevant to EoR (k⊥∼k∥∼ 0.1h Mpc−1) is around 2 K, and the I to Q, U leakage at the edge of the main beam is 0.03. Assuming that the Q, U to I leakage is also around 0.03, the polarized sky structure on this field could create apparent sky structure that would interfere with detection of the EoR signal at the 60 mK in 3D. This is the level of polarization artefact that may occur at the edge of the main beam. Leakage is smaller near the centre. A factor of 10 primary beam leakage correction would be sufficient to bring potential polarized leakage levels below the ≲10 mK expected EoR signal. Polarization calibration made by rastering the pulsar across the main beam, along with mosaic observations, may allow the beam squint terms to be measured and corrected for to achieve this additional factor of 10 reduction in the polarized emission artefacts, allowing a measurement of large-scale structure during the EoR. The required calibration should be attainable using the pulsar technique. There are caveats however; for example, if there is much stronger polarization in neighbouring fields than polarization leakage through the side lobes or diffraction spikes may be an issue. As we have seen in this paper, bright unpolarized sources outside the beam can appear polarized, but at zero Faraday depth.
The GMRT EoR data have resulted in the strongest upper limits to polarized emission at 150 MHz to date. The observed polarization is an order of magnitude smaller than had been expected from higher frequency observations. We note that comparison was not made on the same fields. We estimated the expected beam and depth polarization effects, but these are model-dependent and poorly constrained by the data. Within the uncertainties, either mechanism could account for the lack of observed polarization. If beam depolarization were the cause, it would imply a significant small-scale magnetic field with structure at l > 10 000. Others have reported depth depolarization at similar RMλ2, so depth depolarization may play a role here as well. When compared to the sensitivities required for EoR observations, the upper bounds we report are still two orders of magnitude larger than the sensitivities required. The measured polarization leakage at the GMRT is less than 10 per cent, so a modest leakage correction, accurate to 10 per cent, should make EoR observations feasible. The low polarized emission we report is promising for the search for reionization.
We acknowledge helpful discussions with Ethan Vishniac, Michel Brentjens, Marc Antoine Miville-Deschenes and Judd Bowman. We acknowledge financial support by CIfAR, NSERC and NSF. The GMRT is operated by the National Centre for Radio Astrophysics (NCRA). We thank the NCRA for their extended support in this project.
APPENDIX A: POWER SPECTRUM CONVENTIONS
We report our final results in this paper in terms of (limits on) polarized power spectra. This is appropriate because the theoretical predictions for 21 cm fluctuations from the reionization epoch are usually presented as power spectra, either 2D (Cl) or 3D [P(k)]. An additional advantage is that unlike many other statistical measures, power spectra can be measured unambiguously over some range of l (or k) from data that are missing zero-length baselines or contain significant instrument noise. The purpose of this appendix is to summarize the definition of the power spectrum for both total intensity and polarization, and show how they can be derived from visibilities.
A1 2D power spectra: total intensityis the mean intensity of the field. equation (A1),
The power spectrum can also describe the rms fluctuations of a sky map smoothed by a beam. For example, convolving I(n) with a Gaussian beam of FWHM θb to generate a smoothed map Ism(n) is equivalent to multiplying its Fourier transform by the Fourier transform of a Gaussian:times the power spectrum of the original map, so the rms fluctuation of the smoothed map is
We finally relate the power spectrum to the variance of visibilities observed on an interferometer. We first consider the visibility V(u) observed on an interferometer with baseline u (separation transverse to line of sight in units of the wavelength λ), primary beam FWHM θb. For our angular coordinate n, we take the origin to be at the centre of the primary beam. A point source at position n with flux F (in Jy) will give a visibilityequation (A1), and noting that the rapid modulation of the fringe phase implies that the mean intensity gives no contribution, we find equation (A2):
In our case, if the visibilities are calibrated in Jy, then the intensity I(I) has units of Jy sr−1 and the power spectrum is determined in Jy2 sr−2. However, we want the power spectrum in units of K2, which means that the visibilities have to be converted to K sr units. The conversion factor isequation (A12) yields or 0.058 rad, this becomes
Unlike, for example, the microwave background or extragalactic large-scale structure, Galactic foregrounds are not statistically isotropic. Nevertheless, one may still consider the power spectrum in a particular region of sky using equation (A14).
A2 2D power spectra: polarization
The same Fourier decomposition used for the total intensity can also be applied to linear polarization. One may decompose the linear polarization in analogy to equation (A1) as
In studies of the microwave background, it is common to decompose the polarization into so-called E- and B-mode components, in which one rotates the reference axes of the Stokes parameters to be aligned with the Fourier wavevector l:Stebbins 1996; Kamionkowski, Kosowsky & Stebbins 1997; Seljak 1997; Zaldarriaga & Seljak 1997).2 One can then define power spectra CEl and CBl in analogy to equation (A2). This is useful because in the case of the microwave background the source of polarization is the Thomson scattering, and hence the polarization direction is correlated with spatial morphology, leading to different power spectra in the E and B modes. For low-frequency Galactic emission, however, the Faraday rotation makes it unlikely that the observed polarization angles would have anything to do with the spatial morphology. Therefore, we instead work in terms of the sum of the power spectra, CPl≡CEl+CBl, which incorporates both Stokes parameters. Then, we have equation (A4), is
A3 3D power spectra: total intensity
The predicted 21 cm signal is 3D in the sense of having structure in both the angular and radial (frequency) directions. The statistic most commonly computed by theorists is power spectrum, P(k, μ), and the power per logarithmic range in wavenumber, k3P(k, μ)/(2π2) (Bharadwaj & Ali 2004; Barkana & Loeb 2005). To define this statistic, we begin by decomposing the intensity into its Fourier modes in analogy to equation (A21):
Even in the case of a cosmological signal, the power spectrum may depend on the direction of k relative to the line-of-sight vector which we take to be . This is because the observed frequency of emission depends not just on the location of the emitting hydrogen gas, but also on its radial velocity.3 We therefore introduce
In order to connect P(k, μ) to observations, we need to find the conversion from the comoving coordinates r to the observables, angular position n and frequency ν. In the direction transverse to the line of sight, the relation is r⊥=Dn, where D is the comoving angular diameter distance to redshift z=λ/(21 cm) − 1. In the direction along the line of sight, the relation is more complicated:
We can now construct the 2D power spectrum, which can be measured via equation (A12), in terms of the 3D power spectrum. We construct a 2D image I2D(n) with bandwidth Δν,equation (A1), we can find the Fourier transform of the 2D map. For l≠ 0, equation (A28), yielding and approximate equation (A14), this allows one to estimate the 3D power spectrum from frequency-stacked 2D visibilities.
A4 3D power spectra: polarization
The machinery used to describe intensity fluctuations can also be used to describe polarization fluctuations. Polarized emission from EoR is expected to be intrinsically very small and if present would be strongly affected by the Galactic Faraday rotation (although see Babich & Loeb 2005). Polarized emission from the Galaxy is a potentially severe foreground, and we characterize it by describing its power spectrum using the same conventions used for the EoR signal. To do this, we decompose the polarized flux density into Fourier modes just as was done for the intensity:
The 3D polarized power spectrum can be obtained by stacking in analogy to equation (A29). We define stacked 2D polarization maps,equation (A31) gives, in the polarization case, that the 2D power spectra of the complex maps P2D(n, τ) are
In the case of the cosmological signal in the unpolarized intensity, we identified k∥ with the line-of-sight wavenumber of a mode in large-scale structure. In the case of polarized Galactic synchrotron radiation, the rapid variation of the signal with frequency is instead due to the Faraday rotation. In fact, equation (A35) is (at least over a limited range of frequency) an implementation of RM synthesis, with the exponent replaced by e2πiτν. That is, it is the RM synthesis map at the Faraday depth
We thus conclude that the mathematical procedure to extract the EoR power spectrum from a foreground-cleaned data cube is (in principle!) similar to that of computing a 2D power spectrum on RM synthesis maps, except that it is applied to the real total intensity I (via equation A29) rather than complex polarization Q+ iU. In particular, a simple type of polarization-to-intensity leakage such asequation (A38) if (γQ, γU) is frequency-independent (or varies slowly with frequency), since it results in an apparent Stokes I signal with a sinusoidal variation in ν. In this case, the observed power spectrum would be equation (A38).