Abstract

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.

1 INTRODUCTION

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.

3 DATA

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 forumla5 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.

Table 1

The apparent brightest 10 point sources in the observing field, listed in descending order of their primary beam attenuated 7C catalogue flux, listed in the F7C-BEAM column; F7C indicates their original flux measured from the 7C Radio Survey at 151 MHz, F74 from the VLSS at 74 MHz, and F365 from the Texas Survey of Radio Sources at 365 MHz. α is the spectral index calculated from the three surveys and Δα are the respective errors on the spectral indices.

Source RA (°) Dec. (°) F7C-BEAM (Jy) F7C (Jy) F74 (Jy) F365 (Jy) α Δα 
3C 200 126.855 29.3128 2.135 13.435 26.78 6.407 0.93 0.07 
5C 7.245 126.487 26.7331 1.511 1.536 2.95 0.731 0.89 0.04 
B2 0825+24 127.172 24.6107 0.992 2.922 10.58 0.671 1.77 0.08 
B2 0819+25 125.560 25.6421 0.951 1.701 3.15 0.845 0.84 0.04 
PKS 0832+26 128.798 26.5764 0.706 2.124 3.43 1.384 0.62 0.10 
B2 0829+28 128.082 27.8792 0.682 1.634 2.64 1.030 0.63 0.08 
5C 7.111 124.785 26.2457 0.670 1.786 3.44 1.086 0.82 0.20 
B2 0832+26A 128.770 26.1785 0.638 1.960 3.48 1.156 0.75 0.11 
B2 0828+27 127.841 26.9659 0.625 0.889 1.63 0.400 0.87 0.03 
B2 0816+26B 124.819 26.7014 0.569 1.417 2.82 0.657 0.94 0.05 
5C 7.223 126.028 26.4672 0.515 0.584 1.060 0.306 0.81 0.06 
Source RA (°) Dec. (°) F7C-BEAM (Jy) F7C (Jy) F74 (Jy) F365 (Jy) α Δα 
3C 200 126.855 29.3128 2.135 13.435 26.78 6.407 0.93 0.07 
5C 7.245 126.487 26.7331 1.511 1.536 2.95 0.731 0.89 0.04 
B2 0825+24 127.172 24.6107 0.992 2.922 10.58 0.671 1.77 0.08 
B2 0819+25 125.560 25.6421 0.951 1.701 3.15 0.845 0.84 0.04 
PKS 0832+26 128.798 26.5764 0.706 2.124 3.43 1.384 0.62 0.10 
B2 0829+28 128.082 27.8792 0.682 1.634 2.64 1.030 0.63 0.08 
5C 7.111 124.785 26.2457 0.670 1.786 3.44 1.086 0.82 0.20 
B2 0832+26A 128.770 26.1785 0.638 1.960 3.48 1.156 0.75 0.11 
B2 0828+27 127.841 26.9659 0.625 0.889 1.63 0.400 0.87 0.03 
B2 0816+26B 124.819 26.7014 0.569 1.417 2.82 0.657 0.94 0.05 
5C 7.223 126.028 26.4672 0.515 0.584 1.060 0.306 0.81 0.06 

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 α,  

1
formula
Each L(ν)T(t) product function is a singular eigenvector of the matrix.

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.

5 CALIBRATION

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.

Figure 1

Apparent pulsar polarization angle. The solid, dotted, dashed, long dashed are the mornings of the 18th, 17th, 9th and 16th, respectively. The vertical axis is twice the polarization position angle, i.e. tan−1(U/Q). The long-term change seen in the data (the dip towards the centre of the plot) is likely due to change in airmass and geomagnetic field angle, while the variation between each point and those immediately following is consistent with pulsar angle variations at the source. The night of the 18th showed a 2 hour period in which the polarization angle departed from the slowly changing behaviour characteristic of the rest of the data set. This is likely an ionospheric anomaly. For this plot, angles are arithmetically averaged over 1 min intervals, then the median angle for each 9 min interval is plotted.

Figure 1

Apparent pulsar polarization angle. The solid, dotted, dashed, long dashed are the mornings of the 18th, 17th, 9th and 16th, respectively. The vertical axis is twice the polarization position angle, i.e. tan−1(U/Q). The long-term change seen in the data (the dip towards the centre of the plot) is likely due to change in airmass and geomagnetic field angle, while the variation between each point and those immediately following is consistent with pulsar angle variations at the source. The night of the 18th showed a 2 hour period in which the polarization angle departed from the slowly changing behaviour characteristic of the rest of the data set. This is likely an ionospheric anomaly. For this plot, angles are arithmetically averaged over 1 min intervals, then the median angle for each 9 min interval is plotted.

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 LR 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.

Table 2

Sources in the observing field that are used to estimate the overall calibration flux scale.

Source RA (°) Dec. (°) F7C (Jy) F151 (Jy) Error 
5C 7.245 126.487 26.7331 1.536 1.74 0.13 
B2 0819+25 125.560 25.6421 1.701 1.51 0.11 
B2 0828+27 127.841 26.9659 0.889 0.93 0.05 
5C 7.223 126.028 26.4672 0.584 0.53 0.09 
Source RA (°) Dec. (°) F7C (Jy) F151 (Jy) Error 
5C 7.245 126.487 26.7331 1.536 1.74 0.13 
B2 0819+25 125.560 25.6421 1.701 1.51 0.11 
B2 0828+27 127.841 26.9659 0.889 0.93 0.05 
5C 7.223 126.028 26.4672 0.584 0.53 0.09 

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 forumla3. This profile is a good fit to our holographically measured primary beam power out to about 2° radius.

7 MODELLING

After the primary calibration, the dominant polarized source is the pulsar. The second brightest apparent polarized source is 3C 200, forumla5 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.

Figure 2

Polarized point-source map at the Faraday depth φ= 0 rad m−2. The radius (mean of semimajor and semiminor axes) of each ellipse is proportionate to apparent flux. Ellipticity is proportionate to polarization. Map is averaged over 16 MHz. The asterisk in the centre is the pulsar. 5C 7.245 has been omitted for clarity. The ellipticity is plotted eight times larger than the apparent polarization. The global polarization angle was rotated to align the pattern radially. The feed support structures could lead to such a polarization pattern (see the text). The dashed ring is the half power point of the primary beam. Sources outside that radius appear more strongly polarized.

Figure 2

Polarized point-source map at the Faraday depth φ= 0 rad m−2. The radius (mean of semimajor and semiminor axes) of each ellipse is proportionate to apparent flux. Ellipticity is proportionate to polarization. Map is averaged over 16 MHz. The asterisk in the centre is the pulsar. 5C 7.245 has been omitted for clarity. The ellipticity is plotted eight times larger than the apparent polarization. The global polarization angle was rotated to align the pattern radially. The feed support structures could lead to such a polarization pattern (see the text). The dashed ring is the half power point of the primary beam. Sources outside that radius appear more strongly polarized.

Figure 3

Apparent polarization of sources at the Faraday depth φ= 0 rad m−2. The small squares are the apparent polarizations. The asterisks are multiplied by the primary beam, which is an estimate of their impact to the total polarization in the image. We see that most of the apparent polarized flux contribution is from the periphery of the beam.

Figure 3

Apparent polarization of sources at the Faraday depth φ= 0 rad m−2. The small squares are the apparent polarizations. The asterisks are multiplied by the primary beam, which is an estimate of their impact to the total polarization in the image. We see that most of the apparent polarized flux contribution is from the periphery of the beam.

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 uv−λ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.

9 RESULTS

We compute the angular power spectrum at different rotation measures. The visibilities were directly binned into a power spectrum assuming a forumla3 Gaussian beam. In this case,  

2
formula
where V is the visibility in jansky for a baseline u measured in wavelengths and θb is the FWHM of a Gaussian beam. The left-hand side is measured in K2. It follows that l= 2πu.

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 uv 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 uv 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 ν,  

3
formula
where VLR=Q+ i U is the complex LR visibility which contains both Q and U parameters. This is a conservative overestimate of true polarization in the sky, since this statistic includes thermal noise, polarization leakage, RFI, etc. By construction, this power is positive-definite. We measured in frequency bins Δν= 1 MHz, so the Faraday rotation within a bin is not expected to play an important role.

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 forumla 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 uv 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.

Figure 4

Power spectrum of unpolarized (I2, long dash) and total polarized emission (P2, dotted) over the central square degree, after subtracting the 31 brightest sources in the primary beam. This polarized signal includes leakage and thermal noise. The upper bound symbols are a 2σ estimate of polarization bound using baseline cross-correlation. They are averages of the polarization power in 1 MHz bands. The lower solid and dot–dashed lines are the total power polarization measured at φ=−14 and 14 over a width Δφ= 3.5, as determined from the map power. Using the cross-correlation technique described in the text, we plot a less conservative estimate with circles for φ=−14 and squares for φ=+14. Open symbols denote negative power, and filled are positive values. For details on the unit conversions, see the appendix.

Figure 4

Power spectrum of unpolarized (I2, long dash) and total polarized emission (P2, dotted) over the central square degree, after subtracting the 31 brightest sources in the primary beam. This polarized signal includes leakage and thermal noise. The upper bound symbols are a 2σ estimate of polarization bound using baseline cross-correlation. They are averages of the polarization power in 1 MHz bands. The lower solid and dot–dashed lines are the total power polarization measured at φ=−14 and 14 over a width Δφ= 3.5, as determined from the map power. Using the cross-correlation technique described in the text, we plot a less conservative estimate with circles for φ=−14 and squares for φ=+14. Open symbols denote negative power, and filled are positive values. For details on the unit conversions, see the appendix.

If one fixes a rotation measure, all frequencies can be summed into a single uv 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 uv 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 forumla 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.

Figure 5

Polarized sky image at φ= 3 rad m−2. Left-hand panel is the Stokes Q intensity, while right-hand panel is the Stokes U. All baselines with |w| < 200 are included. The field of view is forumla78, which is three primary beams. Real polarized sky structure as well as leakage I into Q,U will have higher contrast within the main beam region. Sources that generally have the same contrast inside and outside the main beam area include RFI and the aggregate extended image artefacts of multiple bright sources far outside the main beam. These images do show a weak enhancement of contrast in the main beam region, but the enhancement is too weak to allow confidence in detection of sky structure. For this reason, we interpret the measured power spectra from these images as upper limits to the polarized sky structure. Artefacts from one known bright source 3C 200 are faintly visible as concentric rings centred on the source location. 3C 200 is in the upper centre of the image just outside the main beam.

Figure 5

Polarized sky image at φ= 3 rad m−2. Left-hand panel is the Stokes Q intensity, while right-hand panel is the Stokes U. All baselines with |w| < 200 are included. The field of view is forumla78, which is three primary beams. Real polarized sky structure as well as leakage I into Q,U will have higher contrast within the main beam region. Sources that generally have the same contrast inside and outside the main beam area include RFI and the aggregate extended image artefacts of multiple bright sources far outside the main beam. These images do show a weak enhancement of contrast in the main beam region, but the enhancement is too weak to allow confidence in detection of sky structure. For this reason, we interpret the measured power spectra from these images as upper limits to the polarized sky structure. Artefacts from one known bright source 3C 200 are faintly visible as concentric rings centred on the source location. 3C 200 is in the upper centre of the image just outside the main beam.

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):  

4
formula
where μ is the cosine of the angle between the wavevector and the line of sight, and PT is the 3D temperature power spectrum. It is also common to work in terms of the radial and transverse components of the wavenumber, k=kμ and forumla.

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  

5
formula

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.

Figure 6

The 3D polarization power spectra based on our observed visibilities. These are estimated from the power of maps at fixed rotation measures. The transverse scale has been converted from angular to comoving units at z= 8.5. The parallel scale is converted from φ according to equation (A38). The contours are of Δ2P (see the text), and start at 0.0069 K2 in the left middle near k∼ 0, increasing logarithmically by 10 per contour level. On large scales, k < 0.3 h Mpc−1 our observations set an upper limit equivalent to a few kelvins per ln k in 3D space.

Figure 6

The 3D polarization power spectra based on our observed visibilities. These are estimated from the power of maps at fixed rotation measures. The transverse scale has been converted from angular to comoving units at z= 8.5. The parallel scale is converted from φ according to equation (A38). The contours are of Δ2P (see the text), and start at 0.0069 K2 in the left middle near k∼ 0, increasing logarithmically by 10 per contour level. On large scales, k < 0.3 h Mpc−1 our observations set an upper limit equivalent to a few kelvins per ln k in 3D space.

10 DISCUSSION

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, UI 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 forumla7 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 forumla 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  

6
formula
We could detect such a periodic variation, but do not see such a feature. If the spatial gradients were substantially larger in our field, the periodic structures could occur on scales smaller than we constrain. However, this would require an RM gradient larger than 20 rad m−2 deg−1. Beam depolarization due to presence of a smooth RM gradient seems an unlikely explanation for the low polarization.

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(θc2∼ 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  

7
formula
The observed polarization vector is  
8
formula
This model has an exponential suppression with argument ∝λ4, so a factor of ∼2 suppression at 350 MHz would be a factor of 109 suppression at 150 MHz. Using the dispersion and rotation measures to the pulsar, an emission structure at a distance similar to that of the pulsar, with a Gaussian width 10 per cent of the extent of the rotating medium would produce an order of magnitude attenuation of the net polarization between 350 and 150 MHz.

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  

9
formula
This leads to an oscillating depolarization with a ∼λ−2 envelope, which at 150 MHz gives an additional factor of only 5 suppression compared to the 350 MHz data, less than the increase due to the spectral index. In this model, one would expect 5–10 K polarized emission, which is not observed at l < 103.

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.

10.2 Prospects

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  

10
formula
We observe the two visibilities V1,2 at the same azimuth and elevation, so knowing Δφ allows us to solve for the sky and leakage visibilities.

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 (kk∼ 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.

11 CONCLUSIONS

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.

1
Note that the numbers in table 1 of Vinyaikin et al. (1996) are actually forumla due to their normalization convention for polarized brightness.
2
Some early papers denoted the rotated polarizations by G and C, or ⊕ and ⊗, instead of E and B, or used different normalizations. Also note that this decomposition of polarization also applies to other spin 2 fields, such as gravitational lensing shear.
3
An analogous situation occurs in galaxy redshift surveys (Kaiser 1987).

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.

REFERENCES

Babich
D.
Loeb
A.
,
2005
,
ApJ
 ,
635
,
1
Barkana
R.
Loeb
A.
,
2005
,
ApJ
 ,
624
,
L65
Bharadwaj
S.
Ali
S. S.
,
2004
,
MNRAS
 ,
352
,
142
Brentjens
M. A.
De Bruyn
A. G.
,
2006
,
Astron. Nachr.
 ,
327
,
545
Burn
B. J.
,
1966
,
MNRAS
 ,
133
,
67
Cohen
A. S.
Lane
W. M.
Cotton
W. D.
Kassim
N. E.
Lazio
T. J. W.
Perley
R. A.
Condon
J. J.
Erickson
W. C.
,
2007
,
AJ
 ,
134
,
1245
Cordes
J. M.
Lazio
T. J. W.
,
2002
, preprint (0207156)
De Bruyn
A. G.
Katgert
P.
Haverkorn
M.
Schnitzeler
D. H. F. M.
,
2006
,
Astron. Nachr.
 ,
327
,
487
De Oliveira-Costa
A.
Tegmark
M.
Gaensler
B. M.
Jonas
J.
Landecker
T. L.
Reich
P.
,
2008
,
MNRAS
 ,
388
,
247
Douglas
J. N.
Bash
F. N.
Bozyan
F. A.
Torrence
G. W.
Wolfe
C.
,
1996
,
AJ
 ,
111
,
1945
Dunkley
J.
et al.,
2009
,
ApJS
 ,
180
,
225
Hankins
T. H.
Rankin
J. M.
,
2008
, preprint (arXiv:0802.1202)
Haslam
C. G. T.
Salter
C. J.
Stoffel
H.
Wilson
W. E.
,
1982
,
A&AS
 ,
47
,
1
Haverkorn
M.
Katgert
P.
De Bruyn
A. G.
,
2004
,
A&A
 ,
427
,
549
Haverkorn
M.
Brown
J. C.
Gaensler
B. M.
McClure-Griffiths
N. M.
,
2008
,
ApJ
 ,
680
,
362
Hinshaw
G.
et al.,
2009
,
ApJS
 ,
180
,
225
Iliev
I. T.
Mellema
G.
Pen
U.-L.
Bond
J. R.
Shapiro
P. R.
,
2008
,
MNRAS
 ,
384, 863
Kaiser
N.
,
1987
,
MNRAS
 ,
227
,
1
Kamionkowski
M.
Kosowsky
A.
Stebbins
A.
,
1997
,
Phys. Rev. D
 ,
55
,
7368
Komatsu
E.
et al.,
2009
,
ApJS
 ,
180
,
330
Landecker
T. L.
Wielebinski
R.
,
1970
,
Aust. J. Phys. Astrophys. Suppl.
 ,
16
,
1
Manchester
R. N.
,
1974
,
ApJ
 ,
188
,
637
Manchester
R. N.
Hobbs
G. B.
Teoh
A.
Hobbs
M.
,
2005
,
AJ
 ,
129
,
1993
Myers
S. T.
et al.,
2003
,
ApJ
 ,
591
,
575
Nolta
M. R.
et al.,
2009
,
ApJS
 ,
180
,
296
Seljak
U.
,
1997
,
ApJ
 ,
482
,
6
Sokoloff
D. D.
Bykov
A. A.
Shukurov
A.
Berkhuijsen
E. M.
Beck
R.
Poezd
A. D.
,
1998
,
MNRAS
 ,
299
,
189
Stebbins
A.
,
1996
, preprint (9609149)
Uyaniker
B.
Landecker
T. L.
Gray
A. D.
Kothes
R.
,
2003
,
ApJ
 ,
585
,
785
Vinyaikin
E. N.
Kuznetsova
I. P.
Paseka
A. M.
Razin
V. A.
Teplykh
A. I.
,
1996
,
Astron. Lett.
 ,
22
,
582
Vinyajkin
E. N.
,
2004
, in
Uyaniker
B.
,
Reich
W.
,
Wielebinski
R.
, eds,
The Magnetized Interstellar Medium. Copernicus GmbH, Katlenburg-Linday
 . p.
93
Waldram
E. M.
Yates
J. A.
Riley
J. M.
Warner
P. J.
,
1996
,
MNRAS
 ,
282
,
779
Willott
C. J.
Rawlings
S.
Blundell
K. M.
,
2001
,
MNRAS
 ,
324
,
1
Zaldarriaga
M.
Seljak
U.
,
1997
,
Phys. Rev. D
 ,
55
,
1830

Appendix

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 intensity

We first consider the 2D case, e.g. where we consider diffuse emission at only one frequency. The intensity I(n) is then a function of angular position n. It can be Fourier-decomposed as  

(A1)
formula
where l is the Fourier wavenumber traditionally used in cosmology; it is related to u= (u, v) (which has units of waves per radian) by l= 2πu. Finally, forumla is the mean intensity of the field.

The power spectrum Cl of a diffuse field is a statistical measure defined by  

(A2)
formula
where δ(2) is the 2D Dirac delta function. It is related to the correlation function of the field by use of equation (A1),  
(A3)
formula
In particular, the rms intensity fluctuations in the map can be written as  
(A4)
formula
where in the last equality we have integrated over the direction of l and noted that the area element d2l in the l plane is 2πl dl. This integral is the reason why l2Cl/(2π) is a commonly used quantity: it describes the contribution to the intensity variance coming from a particular range of Fourier modes (specifically, per ln l).

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:  

(A5)
formula
The power spectrum of the smoothed map is then forumla times the power spectrum of the original map, so the rms fluctuation of the smoothed map is  
(A6)
formula

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 visibility  

(A7)
formula
where the first factor gives the fringe phase and the second is the primary beam suppression. In a diffuse field, the observed visibility will be the integral of the intensity over the field of view,  
(A8)
formula
Using equation (A1), and noting that the rapid modulation of the fringe phase implies that the mean intensity forumla gives no contribution, we find  
(A9)
formula
The variance of this visibility can be found from equation (A2):  
(A10)
formula
Thus, the variance of the visibilities is directly related to the power spectrum at values of l within ∼θ−1b of 2πu. In the practical case that the fringe spacing is much smaller than the primary beam (θb|u| ≫ 1) and the power spectrum varies slowly with scale, we may approximate Cl by its value at 2πu and write  
(A11)
formula
This implies that the power spectrum can be determined in terms of the visibilities:  
(A12)
formula

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 is  

(A13)
formula
Inserting this conversion factor into V(u) in equation (A12) yields  
(A14)
formula
If we refer the beam to the intrinsic FWHM of the GMRT primary beam, forumla or 0.058 rad, this becomes  
(A15)
formula

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  

(A16)
formula
and similarly for U.

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:  

(A17)
formula
where αl= tan−1(ly/lx) (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, CPlCEl+CBl, which incorporates both Stokes parameters. Then, we have  
(A18)
formula
and the variance of the polarization, analogous to equation (A4), is  
(A19)
formula
Thus, l2CPl/(2π) can be viewed as the contribution to the polarization variance per ln l.

The same methodology for estimating the intensity power spectrum from the variance of the visibilities can be applied to the polarization. We find that  

(A20)
formula

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):  

(A21)
formula
where r is a 3D comoving position and k is the wave vector. Then, a power spectrum is defined via  
(A22)
formula

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 forumla. 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  

(A23)
formula
which is between −1 and +1. Then, the total variance of fluctuations is  
(A24)
formula
where we have used that the volume element d3k is 2πk2 dμ dk after integrating out the azimuthal angle. For an isotropic power spectrum, P(k, μ) =P(k), we can thus interpret k3P(k)/(2π2) as the variance in intensity per logarithmic range in wavenumber. Often one denotes Δ2(k, μ) =k3P(k, μ)/(2π2). The power spectrum is also sometimes written in components transverse to and radial along the line of sight, P(k, k).

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:  

(A25)
formula
where c is the speed of light, H is the Hubble constant at redshift z and ν is the observed frequency. We can then define the radial conversion factor C≡ (1 +z)/(cHν). It follows that the comoving coordinates are related to observed coordinates via (r, r) ↔[Dn, −C(ν−ν0)], where ν0 is a reference frequency (the zero-point of r is arbitrary).

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 Δν,  

(A26)
formula
Using equation (A1), we can find the Fourier transform of the 2D map. For l≠ 0,  
(A27)
formula
where sinc x≡ sin x/x. In the last equality, the d2k integral is trivial due to the δ-function, and the r integral gives the sinc function. By taking the variance, we extract the 2D power spectrum,  
(A28)
formula
Note that the amount of power in the 2D map depends on the bandwidth Δν. This is a direct consequence of the 3D nature of the EoR signal: fluctuations with |k| > (CΔν)−1 are cancelled by the combination of frequencies.

It is possible to measure the full 3D power spectrum by stacking maps at several neighbouring frequencies with a phase shift to avoid cancellation of frequencies, i.e.  

(A29)
formula
where τ is the phase shift in units of cycles per MHz. This is now a complex map but this causes no difficulty. The power spectrum of I2D(n, τ) can be obtained by the same methodology as equation (A28), yielding  
(A30)
formula
Thus, this stacked map is sensitive to fluctuations with radial structure and wavenumbers in the range k≈ 2πτ/C, within a width δk∼ (CΔν)−1. If the power spectrum is slowly varying as a function of k over this range, one may use the fact that forumla and approximate  
(A31)
formula
By combining with 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:  

(A32)
formula
and similarly for U. One may define a polarized power spectrum via  
(A33)
formula
and similarly for PQQ, PUQ and PUU. Since Q and U are real, it follows that PQQ and PUU are real and non-negative and satisfy PQQ(k) =PQQ(−k) and PUU(k) =PUU(−k); also PQU(k) =P*UQ(−k). A decomposition of the polarization tensor into E and B modes is possible, just as in the 2D case, but is probably not interesting if the radiation has propagated through a Faraday thick medium. We therefore focus on measuring the power spectrum PP(k) =PQQ+PUU:  
(A34)
formula

The 3D polarized power spectrum can be obtained by stacking in analogy to equation (A29). We define stacked 2D polarization maps,  

(A35)
formula
Repeating the steps that lead to equation (A31) gives, in the polarization case, that the 2D power spectra of the complex maps P2D(n, τ) are  
(A36)
formula
where ℑ denotes the imaginary part. Note that for a signal that respects statistical parity invariance, the second term here must vanish, but this may not be true for the Galactic emission in a particular patch of sky. This term may be removed by averaging signals from positive and negative τ.

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 forumla replaced by e2πiτν. That is, it is the RM synthesis map at the Faraday depth  

(A37)
formula
Using the conversions from τ to k (k= 2πτ/C), this becomes  
(A38)
formula

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 as  

(A39)
formula
will leak polarized structure at the Faraday depth φ into the intensity map at wavenumber k given by equation (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  
(A40)
formula
More complex types of frequency-dependent polarization leakage can occur in real EoR observations, especially if observations at the same (u, v) but different frequencies are taken non-simultaneously [as is necessary at the GMRT, since a given point in the (u, v) plane corresponds to a different physical baseline length at a different frequency]. One example would be polarized side lobes that rotate relative to the sky due to the altazimuth mount. Another would be changes in the ionospheric Faraday rotation between the epochs of observation. Finally, the beam squint will vary across the bandpass. These frequency-dependent leakages are a concern because polarization at any Faraday depth can be injected into the Stokes I map at wavenumbers that do not obey equation (A38).