The GMRT-EoR Experiment: A new upper limit on the neutral hydrogen power spectrum at z~8.6

We present a new upper limit to the 21cm power spectrum during the Epoch of Reionization (EoR) which constrains reionization models with an unheated IGM. The GMRT-EoR experiment is an ongoing effort to make a statistical detection of the power spectrum of 21cm neutral hydrogen emission at redshift z~9. Data from this redshift constrain models of the EoR, the end of the Dark Ages arising from the formation of the first bright UV sources, probably stars or mini-quasars. We present results from approximately 50 hours of observations at the Giant Metrewave Radio Telescope in India from December 2007. We describe radio frequency interference (RFI) localisation schemes which allow bright sources on the ground to be identified and physically removed in addition to automated flagging. Singular-value decomposition is used to remove remaining broadband RFI by identifying ground sources with large eigenvalues. Foregrounds are modelled using a piecewise linear filter and the power spectrum is measured using cross-correlations of foreground subtracted images.


INTRODUCTION
Between the recombination epoch at z ∼ 1100 and the first round of star formation at z ∼ 10 the Universe was filled with neutral hydrogen. This neutral gas is thought to have produced 21 cm hyperfine emission with an effective continuum brightness temperature between -500 and 30 mK ) using z = 8.6 and WMAP 7-year cosmological parameters (Komatsu et al. 2010), which have formal error bars of roughly 10 per cent. As the first stars formed, beginning at local peaks in the matter density, the hydrogen gas was locally ionized. This is the start of a period of cosmic evolution called the Epoch (or Era) of Reionization (EoR). Over time the ionized cells grew and overlapped, creating a patchwork of ionized and neutral cells. This general patchy topology is well motivated by theory and simulations ⋆ Email:paciga@astro.utoronto.ca † Email:pen@cita.utoronto.ca (e.g. Furlanetto et al. 2004;Iliev et al. 2006;McQuinn et al. 2007;Zahn et al. 2007; see Trac & Gnedin 2009 for a review) though the exact properties are poorly constrained. Eventually, by z ∼ 6, the ionization of the broadly distributed material was complete, leaving only rare pockets of neutral gas.
In addition to directly constraining the redshift of last scattering z ls ∼ 1100, Cosmic Microwave Background (CMB) polarization data have been used to constrain the redshift to the range 8.0 < z < 12.8 (WMAP7; Komatsu et al. 2010) under the assumption it was instantaneous, although evidence from the CMB indicates that reionization was an extended process (Dunkley et al. 2009). There is also substantial information on the neutral content of the intergalactic medium (IGM) since z ∼ 6 via observations of the Gunn-Peterson trough (Gunn & Peterson 1965) in quasar spectra (Djorgovski et al. 2001;Becker et al. 2001;Fan et al. 2006; Willott et al. 2007). Lyman-alpha (Ly-α) absorption provides a sensitive probe of neutral gas density and can be used to constrain reionization (e.g., Fan et al. 2002). With these two techniques we can observe the start and the end of the neutral era. However, neither CMB nor quasar observations allow detailed examination of the reionization era itself. For this it has long been proposed to use 21 cm fluctuations in the brightness temperature of neutral hydrogen (Sunyaev & Zeldovich 1975;Hogan & Rees 1979). The 21 cm power spectrum, resulting from a combination of the patchy ionized and neutral medium and the underlying mass power spectrum, is generally considered one of the most promising signals (Scott & Rees 1990;Zaldarriaga et al. 2004) and much attention has been paid in the literature toward designing suitable experiments to detect it (e.g., Morales & Hewitt 2004;Morales 2005;Bowman et al. 2006;Harker et al. 2010).
Here the observing frequency is 1420 MHz(1 + z) −1 , in the very high frequency (VHF) range of the radio spectrum near 150 MHz. Several programs are underway to study the EoR in the VHF band including LOFAR 1 (Kassim et al. 2004;Röttgering et al. 2006; for the EoR case see, e.g., Zaroubi & Silk 2005;Harker et al. 2010), MWA 2 (Lonsdale et al. 2009; for the EoR case see, e.g., Bowman et al. 2006;Lidz et al. 2008), PAPER 3 (Parsons et al. 2010), 21CMA 4 (a.k.a. PaST; Peterson et al. 2004), in addition to the current GMRT 5 program, and the first goal of such efforts is to determine the redshift at which roughly half the volume of the universe is ionized. Because the 21 cm emission is an isolated single line, these observations can provide three dimensional information over a large redshift range (Loeb & Zaldarriaga 2004;Furlanetto et al. 2006), allowing a broad search for the half-ionization redshift. In the future 21 cm tomography has the potential to add strong constraints on cosmological parameters (e.g., McQuinn et al. 2006;Cooray et al. 2008;Mao et al. 2008;Furlanetto et al. 2009;Masui et al. 2010).
The brightness temperature of the 21 cm line relative to the CMB is determined by the ionization fraction of hydrogen and the spin temperature of the neutral population, which is in turn governed by the background radiation and the kinetic temperature of the gas (Purcell & Field 1956;Field 1959;Furlanetto et al. 2006). Reionization requires a minimum expenditure of 13.6 eV of energy per hydrogen atom. In contrast, Ly-α photons, when absorbed by a neutral atom, are quickly re-emitted. Ly-α photons are said to undergo 'resonant scattering' and each (multiplyscattered) photon can affect many atoms before cosmic redshifting makes them ineffective (Wouthuysen 1952;Field 1959;Chuzhoy & Shapiro 2006). This means Ly-α pumping of the hyperfine transition requires only about 1 per cent of the UV flux required for ionization. Assuming a gradual increase in UV flux with time, well before flux levels for reionization are reached, Ly-α pumping will couple the spin temperature to the kinetic temperature of the gas. If there was no source of heat at that era other than a weak UV flux, the gas kinetic temperature must have been at its adiabatic expansion value 1.7 K at z = 8.5, and neutral gas will produce a signal due to absorption of CMB photons in excess of stimulated emission and ionized structures would be seen as low-brightness regions on the sky (Chen & Miralda-Escudé 2004). In such a cold-gas model the brightness temperature of the neutral gas against the CMB can be as low as -500 mK.
A small fraction of the mass will have collapsed into minihaloes, which have a temperature higher than this adiabatic temperature . This fraction may even be substantially smaller due to non-perturbative velocity flows (Tseliakhovich & Hirata 2010;Dalal et al. 2010). No study was found on the impact of these non-linear effects on the Ly-α pumped IGM temperature. At linear order, the positive and negative over and under densities cancel, and the non-linear collapse fraction is still small, so we expect the realistic value to be similar to the adiabatic prediction.
Alternatively, X-rays from supernovae or QSOs might have heated the IGM above the CMB temperature of ≈ 30 K before reionization (Madau et al. 1997;Chen & Miralda-Escudé 2004). X-ray heating between 30 K and 10 000 K will result in a largely neutral but warm IGM, which would be seen in emission. In the limit that Ts ≫ 30 K, the volume emissivity becomes independent of temperature. Patchy X-ray heating can also result in large angular scale structure (Alvarez et al. 2010). The sky brightness temperature of the neutral gas has an asymptote at ∼ 30 mK .
The cosmic luminosity of X-rays at z ∼ 9 is not known and difficult to estimate (see e.g., Dijkstra et al. 2004;Salvaterra et al. 2007;Guo et al. 2009). If the rate of core collapse SNe at high redshift matches that of today, the X-ray output at z ∼ 9 would have been sufficient to raise the IGM temperature above the CMB temperature at the onset of reionization. However, the mechanism of and factors affecting core-collapse are not known and most numerical models appear to generically not result in SNe at all (Mezzacappa 2005). It is possible that at high z, core collapse SNe were not as abundant as today, and the IGM was still in absorption during the EoR. We therefore consider two limits: one where the IGM is still cold and T b = −500 mK, and another where it is heated above CMB and T b = 30 mK. A general parametrization of these scenarios was recently proposed by Pritchard & Loeb (2010).
The GMRT-EoR project uses the Giant Metrewave Radio Telescope (GMRT; Ananthakrishnan 1995) near Pune, India, to make measurements of the power spectrum of the neutral hydrogen signal with the hope of characterizing the structure in the range 8.1 < z < 9.2. The FWHM of the GMRT primary beam at 150 MHz is 3.3 • which provides a cylindrical comoving survey volume of (280 h −1 Mpc) 3 , with about equal dimensions in three directions. The primary sensitivity comes from the compact central core which is contained within ≈ 1 km, or 30 dish diametres, which gives images with ∼ 30 2 resolution elements. As recorded, the data have high spectral resolution along the line of sight, but we bin the data for comparable transverse and radial resolution. Each data cube has ∼ 30 3 ≈ 30000 resolution elements so even though the signal to noise ratio of each element is less than one, the spatial power can be measured to 1 per cent. A statistical measurement of the power spectrum provides the best hope of pinning down the EoR half-ionization redshift.
In section 2 we describe the data and analysis, including RFI removal and foreground filters used. In section 3 we extend the analysis to cross-correlations of multiple nights, and a measurement of the power spectrum. We then conclude in section 4. All distances are in comoving coordinates.

GMRT and Data Description
The GMRT is a radio interferometer consisting of 30 antennas, each with a diametre of 45 m. Fourteen of these are arranged in a dense central core within 1 km which allows the high brightness sensitivity required to search for the dim EoR signal (Pen et al. 2009). The longest separation between antennas is about 25 km. For this experiment the telescope was operated at 150 MHz with a 16.7 MHz bandwidth. For 21 cm (1421 MHz) emission, this frequency range probes a redshift range of 8.1 < z < 9.2. The highest angular resolution at this observing frequency at GMRT is about 20 arcseconds. Each antenna provides a pair of left and right circularly polarized signals, which are passed through a new signal processing system, built in part for this project, called the GMRT Software Backend (GSB; Roy et al. 2010). The final recorded visibilities have a resolution of 7.8 kHz and 1/4 s in frequency and time, respectively.
The field is centred on the pulsar B0823+26. This pulsar has a period of 0.53 seconds and an average flux of 350 mJy at 150 MHz (Hobbs et al. 2004). It is situated in a relatively cold part of the sky at galactic latitude ≈ 30 • with few nearby bright sources. The on-pulse flux is about 6 Jy, brighter than all other sources in the field, making it a good calibrator. The calibration of this data was described in Pen et al. (2009) and is summarized below. The positions of the brightest sources in the field are shown in Fig. 1.
The data are folded in time into 16 'gates' such that the pulse from the pulsar is contained within one gate. The period of the pulsar is much shorter than the time-scale at which ionospheric fluctuations dominate, so by comparing the 'on' gate with the neighbouring 'off' gates, everything in the field that is constant over the period of the pulsar can be removed. This includes sky sources and most radio frequency interference (RFI). By comparing the pulsar signal received by each antenna, the relative system gain of each antenna can be calibrated. This technique allows calibration of both phase and polarization in real time, with 100 per cent observing efficiency, since no additional time is needed for phase or flux calibration. Errors in clock synchronization sometimes cause the pulsar pulse to straddle two gates. In these cases it is necessary to include the flanking gates when identifying the pulsar signal. This is not ideal since the pulsar signal is diluted over multiple gates, but the correction is limited to only the portion of observations that require it.
Since the pulsar amplitude varies from pulse to pulse, the absolute system gain needs to be measured separately. This is done using a noise injection system, which the GMRT software backend decodes to calibrate the absolute gain of the system, while a sky radio source is used to transfer the noise source calibration on to the sky. The primary flux calibrator used was 5C 7.245, a radio galaxy at z ≈ 1.6 (Willott et al. 2001), located within the field of view (see B2 0816+26B 5C 7.223 B 0823+26 Figure 1. The 12 brightest sources in the field used for these observations. Pulsar PSR B0823+26 is at the centre of the field. The circle denotes the half-power diametre of the main beam. Fig. 1). However, since this radio galaxy is an extended object, it can only be used as a calibrator for short baselines where its structure is not resolved. For antennas in longer baselines we determine the relative gain calibration using the pulsar.

Data Analysis
Data from 2007 December 10, 11, 14, 16, 17, and 18 are included in the current analysis. A software pipeline has been developed to automate the calibration and interference removal steps described below. The current configuration takes approximately 11 hours to process one hour of data on the CITA Sunnyvale computing cluster, or 20 hours on the National Centre for Radio Astronomy (NCRA) HP computing cluster in Pune, India. The NCRA facility can process four hours in parallel, completing one night of observations in slightly less than two days. While the Sunnyvale cluster has the capacity to process as much as 12 hours of data at once, storage capacity, shared usage, and other bottlenecks reduce this significantly.
One of the limiting factors for measurement of the EoR signal is broadband radio frequency interference (RFI), which dominates the signal at 150 MHz. In addition to the standard procedure of flagging bad data, we also use a singular value decomposition (SVD) to remove broadband RFI from the data, and also to identify and remove interference sources. We believe this approach is unique among other RFI mitigation strategies in the radio community.
Narrow line interference was removed by masking points in each frequency bin with an intensity above some threshold. This was done twice in the data reduction pipeline. Input data are initially flagged with a threshold of 8σ on a Gaussian scale, and then again at 3σ after broadband in- In each panel, the horizontal axis is frequency, covering approximately 1 MHz, and the vertical is time, increasing downward, and covering approximately one hour. The grey-scale is the crosspower spectral density. The top row is the C0-C8 baseline and the bottom row the C0-W4 baseline, which are approximately 560 m and 9400 m respectively. Large bright patches indicate broadband RFI, and vertical lines indicate line RFI. The first panel in each row is the initial input data. The second is after the initial 8σ mask, with most interference still visible. The third panel is after removing the largest eigenvalues in the SVD. The broadband interference is no longer visible. The last panel shows the final 3σ mask removing the line RFI, leaving a nearly uniform image.
terference was removed. The first mask can not be too aggressive or the techniques to remove broadband interference, discussed below, will fail. Singular value decomposition (SVD) is used to separate broadband radio sources on the ground from those in the sky. Ground-based sources contribute most to the largest eigenvalues since they do not move as a function of time with respect to the array, while sky sources rotate. One hour of data has 14396 time records, and each record has about 7.5 million entries corresponding to the number of frequency channels and baselines between the 60 antennas. This is treated as a matrix with each time record as one row. The 50 largest eigenvalues are identified through an SVD and flagged as noise to be removed. A sample of the data at a few intermediate stages showing the successful removal of both line and broadband RFI can be seen in Fig. 2. The RFI patterns in (u, v) space both before and after the SVD are illustrated in Fig. 3.

Physical Removal of Radio Interference Sources
While some RFI can be removed in the data analysis, this risks removing sky signal as well. Ideally one would prevent RFI from occurring at all by identifying and correcting the physical sources. We can take advantage of the fact that as an interferometer GMRT is able to make images of the near- field, and create maps of bright RFI sources near the GMRT itself. Candidate sources detected by a single baseline appear as a hyperbola of equal light arrival time in near-field images of a single SVD mode. When a source is detected by many baselines, the corresponding hyperbolas intersect at a single point in the image. These near-field images, an example of which is shown in Fig. 4, become the 'RFI maps' which are used to isolate interference sources. For sources with a high enough duty cycle, the GPS position given by the RFI maps are accurate enough that a handheld yagi antenna and portable radio receiver can be used to find the precise source. With the handheld antenna direction to a source can be determined by ear to better than 30 degrees, and thereafter localised by successive triangulation.
By comparing the position of a noise transmitter relative to the candidate interference source both physically and in the RFI map, we can confirm whether we have identified the correct source. Occasional misidentifications are likely due to the intermittent nature of the sources; there is no guarantee that the radiating source in the images will still be active when an attempt at identification is made. Since the beginning of this effort, calibration has been improved to locate sources within about 100 m. Sources could in principle be located using only the array with a precision 10 m, but the accuracy of available maps, GPS equipment, and other factors limit the real-world precision.
Candidate RFI sources identified so far include transformers, power line junctions, and loose wires in contact with power lines. Table 1 lists some candidate sources recorded on 2008 December 6, and descriptions. In February 2009, GMRT began a collaboration effort with MSEDCL, which controls the power transmission lines in the area around GMRT, to assist us in removing local interference sources. A wire hanging over 132 kV high tension power lines (#1 in Table 1) which was by far the brightest interference source within 20 km, appearing in 30 of the 50 largest eigenmodes, was removed in February 2009. Other bright sources have also been removed since then. Although these sources are still present in the raw data from 2007, this procedure has shown that the SVD algorithm applied in the present work is successful in identifying real RFI sources during analysis.

Single-night Images
After RFI removal has been completed in each one hour scan, polarization calibration is used to combine these into images of an entire night, typically of about 8 hours, as can be seen in Fig. 5. This step accounts for leakage between the left and right polarization signals, as well as effects caused by the relative rotation of the array and sky source by taking into account the changing parallactic angle with time.
Although the array elements very nearly lie in a single plane, the rotation of the Earth turns this plane over time with respect to the line of sight. To correctly stack many hours together, one must consider the frequency and directional dependence of the measurement. One effect to consider is the change of the primary beam with frequency. For a given frequency, the primary beam is well defined and independent of (u, v). Deconvolution is difficult, but one can pick the frequency with the smallest beam and restrict the field of view at other frequencies to match. Computationally this can be done by convolving with the ratio of the beams, which is possible because of the commutativity of the primary beam operators.
However, this is a small effect compared to the effects of a w term. Typically relative antenna positions are considered to be on a two dimensional plane with coordinates parametrized by (u, v) in units of the observing wavelength. Since baselines become non-coplanar over the course of the night, a third dimension, w, must be included, without which the image becomes blurred.
For a thorough explanation of w term issues, see for example Cornwell et al. (2008) and, for a more general overview, Thompson et al. (2001). As the observing frequency changes across the band, a single baseline samples different points in (u, v) space. At any given frequency, the (u, v) plane is typically sparsely sampled, but with many frequency channels there will often be data at a single (u, v) point for at least a few channels. One can subtract the average from these, and under the assumption that foreground sources have the same spectral index as the calibrator, this removes the foregrounds very well. Unfortunately, when considering the full (u, v, w) cube, the data is too sparse for this to work. Since the w coordinate is dependent on the position in (u, v) space, and can become very large, the effects of the w term can not be corrected at different frequencies in the same way as the beam.
Strategies for correcting this effect for GMRT are being developed, which build on those used in CMB studies (Myers et al. 2003;Hobson & Maisinger 2002), but we have not yet applied these to our data. In the interim we restrict ourselves to the short baselines for which w is small and such corrections are not needed. The short baselines are the ones most sensitive to the EoR signal, so we retain most of the EoR sensitivity of the array.
The flux scale of the field is set by calibrating relative to the pulsar in the centre of the field for every baseline, receiver, and frequency. Since the pulsar flux is known to be variable, to set the absolute flux scale we identify the source with the highest flux in the sky image, and set this peak value equal to the known flux of that source. The exact value found depends on the angular resolution, determined by the maximum baseline length used. For all days analysed in this work, the brightest apparent source is the radio galaxy 5C 7.245. The flux of this source was measured as 215.1 ± 12.6 mJy at 1415 MHz by Willis et al. (1976) at the Westerbork Radio Telescope, and as 685 ± 47 mJy at 408 MHz in the 5C catalogue (Pearson & Kus 1978) for a spectral index of 0.93 ± 0.07. Following the same argument as Pen et al. (2009), to achieve the best match with other bright sources in the field, we adopt a value of 1.6 Jy for 5C 7.245. Although the flux of this source changes by as much as 10 per cent across the band, this change is less than that due to the uncertainty in the spectral index. Additionally, this galaxy has two components separated by 12 arcsec (Willott et al. 2001). GMRT is capable of resolving this structure with |b| 2600, where b = (u, v) is the distance between the two antennas in the baseline in units of the observing wavelength. To calibrate the flux we require a point source so that the flux is not spread over multiple image pixels. This is achieved by limiting the maximum baseline length used while calibrating to |b| 2000. At much smaller |b|, we become similarly limited by confusion.
Removing the foregrounds adequately is essential for detecting the EoR signal (see, e.g., Wang et al. 2006) and much work has been done on simulating foregrounds (e.g.   Table 1, identified in (a). Morales & Wyithe (2010) and references therein. To remove foreground sources, we take advantage of the fact that the flux of such sources do not vary greatly with frequency. A signal originating from the EoR will appear as additional variation on top of the foreground signal. To model spectrally smooth sources, for every timestamp and baseline we subtract a piecewise linear fit between the median fluxes in frequency bins of 8, 2, or 0.5 MHz. This is simple to implement and has the advantage over a polynomial fit of a having straightforward and local window function which is simple to interpret. Foregrounds which do not vary over this range are removed, while features with variability at higher frequencies remain. This method results in an upper limit to the EoR signal since the measured power may still include residual foreground variation. By subtracting the mean flux over 2 MHz to remove galactic foregrounds, noise levels can be lowered to about 2 mJy for maps of approximately 10 degrees. Bright point sources may play a significant role (Datta et al. 2010;Di Matteo et al. 2004), but are not treated separately here.

2009, and others). For a review see
To first order, this is the same as a boxcar average 6 which we can write as where D(ν) is the input data,D(ν) is the data after filtering, and the window function is Here, δ(x) is the Dirac delta function and H(x, ∆ν) is a step function centred at x with width ∆ν corresponding to our chosen filter, and of unit area. In Fourier space, this window function is  Figure 6. Fourier transform of our chosen window function for ∆ν = 0.5 MHz (solid line). Dashed lines indicate where w(k ) = 0.5, which sets the minimum k to which we are sensitive. This value scales inversely with ∆ν. The peak of w(k ) 2 /k 3 (the dotdashed line) indicates the k to which we are most sensitive under the assumption that the power spectrum is proportional to k −3 .
shown in Fig. 6, where we have denoted the wave number as kν to make explicit that it is in units of inverse frequency. This can be rewritten in terms of k using the fact that in our redshift range 1 MHz ≈ 11.6 h −1 Mpc. When converted to units of hMpc −1 , we write the wave number as k to emphasize that the window function acts on structure along the line of sight only.
After filtering in this way, we denote the kν for which w(kν) = 0.5 as corresponding to the minimum k along the line of sight to which we are sensitive, while smaller k are removed by the filter and will not contribute as strongly to the power spectrum. The three filters of 8 MHz, 2 MHz, and 0.5 MHz correspond to minimum k of 0.04, 0.16, and 0.65 hMpc −1 respectively. Since the power spectrum is thought to be proportional to k −3 (Iliev et al. 2008), the line w(k ) 2 /k 3 indicates the k to which we are most sensitive.
The 0.5 MHz foreground filter reduces the peak and RMS values by a factor of 50, as can be seen in Fig. 7. The filter is most effective within the primary beam. Fig. 8 demonstrates the dominance of RFI in the filtered dirty maps, and the improvement that the SVD step provides. Without the SVD, the maps in Fig. 7 would have been a factor of 4 noisier.

Differences Between Nights
Since RFI and foreground signals are so much larger than the EoR signal, small errors in the subtraction of foregrounds could easily result in a spurious signal. To avoid this we use cross-correlations between multiple days to make a statistical measurement.
To gauge how successful a cross-correlation measurement might be, we are interested in the relative similarity of the different nights, which we can gauge by taking the difference of visibilities. Days which subtract well will show mostly noise, while days which subtract poorly will still show Figure 7. Data from Dec 10. The top row is the sky image with a 11.4 degree field of view, with |b| < 200 and |u| > 60 binned in the (u, v) plane by ∆u = 5. The bottom row is the visibilities in the same range with ∆u = 0.4 to show structure. The left column is before any foreground subtraction. In this image the dominant source is B2 0825+24 (or 4C 24.17) just south of the FWHM of the primary beam, with a peak value of 2.2 Jy. RMS within the beam is 343 mJy. The right column is after a 0.5 MHz subtraction. The peak value of this image is 47 mJy with an RMS of 6.2 mJy, lower by a factor of about 50. If put on the same grey-scale as the image without foreground subtraction no features would be visible. The dominant source after this filter is 3C 200, well outside the beam. Figure 8. Dirty sky images before the SVD RFI removal step (left) and after the RFI removal (right) with |b| < 200, after the 0.5 MHz foreground subtraction. The field of view is 11.4 degrees. RFI clearly dominates the map on the left, with peak value of 73 mJy and an RMS of 18 mJy. After RFI removal, the peak drops to 32 mJy with an RMS of 4.2 mJy. evidence of bright sources. Visibilities are only used when data exists at the same position in both days. This serves as a test of the important of ionospheric fluctuations. Though the ionosphere is calibrated for along the light of sight to the pulsar in the centre of the field, the ionosphere could change across the field of view. Pairs of days for which this change Table 2. Peak flux and RMS of the difference of each two day pair available, in mJy, with a maximum (u, v) distance of 600, which corresponds to a maximum baseline of 1.2 km. All values are in mJy. To remove foregrounds, a 2 MHz linear filter was applied.

Unfiltered
With 2  away from the field of view is different will not subtract well.
The results of the subtractions are shown in Table 2. It can be seen in Table 2 that December 11 gives consistently poor results in both the unfiltered and filtered images, and was thus excluded from all subsequent analysis. From these maps we can also conclude that RFI is not dominating the differences from day to day. If this had been the case, since RFI is typically isolated in (u, v) space, it would be visible across the whole sky image, including far outside the primary beam. However, it can be seen in Fig. 7 that structure decreases rapidly away from the centre of the field, indicating an astronomical source. This is true of all subtraction pairs.

Cross-correlations and Power Spectra
The power spectrum of sky structure can be determined directly from the visibilities . To find the cross-power of two days, we take the product of the visibilities after gridding with a cell spacing equal to the size of the beam (∆u = 20). Calculating the power spectrum from the visibilities instead of from the correlations in the sky image takes advantage of the fact that the visibilities have a nearly diagonal correlation matrix in the noise (White et al. 1999). We then find the weighted average of visibilities in annuli of (u, v) space which gives the power in units of Jy 2 . Since the amount of data at large |b| decreases, we increase the width of each successive annuli by 60 per cent with increasing |b|. The smallest bin width is equal to the size of the beam. This prevents large artificial variability in power due to sparse sampling. Visibilities are weighted 'naturally', by the inverse of the noise. Since we expect our sensitivity to the EoR signal to diminish rapidly with increasing baseline length, we look only at the first few points, averaged over all possible cross-correlations. Additionally, it is known that the SVD will introduce a loss of power at low |u|. To avoid this, we impose a limit of |u| > 60 when taking the cross-correlations, determined by requiring that the power It can be seen that power is lost in the SVD. Similarly, the solid and dashed lines are before and after the SVD step, respectively, this time with a |u| > 60 limit imposed. In this case the two lines are almost identical, meaning the SVD had little effect on the total power. spectrum before and after the SVD differ by less than 1σ as shown in Fig. 9. The part of the (u, v) plane that is lost with this cut can be seen in Fig. 3. The power spectrum of the cross-correlation can be converted to units of K 2 using l 2 2π C l l=2π|b| = |b| 11.9 is the visibility coordinate in units of wavelength, C l is the power measured in K 2 , θ b is the primary beam size, and ν is the wavelength (Pen et al. 2009). The quantity in the angle-brackets on the right is equal to the power in Jy 2 found above. This conversion is written in terms of the GMRT observations with a primary beam of θ b = 3.3 • and ν = 150 MHz. If gridding in the (u, v) plane is too fine, the data become noisy, while very coarse gridding requires the assumption that the data are constant across the whole cell. As mentioned, we use a gridding equal to the size of the beam. A fully optimized estimate would require a maximum likelihood code which is being worked on in the context of the Cosmic Background Imager (CBI) gridder (Myers et al. 2003). Fig. 10 shows the weighted-average of all cross-correlation pairs, excluding December 11. Bernardi et al. (2009) reported a power spectrum of foregrounds without subtraction in the galactic plane at a level comparable to our measurement both here and in Pen et al. (2009).
We have plotted the power spectra with 2σ bootstrap errors (Efron 1979) which were derived as follows: Using five days of data, there are ten possible cross-correlations. From these, ten are randomly sampled, with replacement, resulting in a slightly different power spectrum. This is repeated 10 4 times, and the variance on this set of power spectra is calculated to give the error on the original. Formally this quantifies the error when taking independent samples of a statistical distribution. In our case this is not a rigorous er- Each point is shown with a 2σ upper limit derived from a bootstrap error analysis, which is in most cases smaller than the size of the point. The points are logarithmically spaced as described in the text, from left to right covering the ranges 377 < l < 578, 578 < l < 899, and 899 < l < 1414. Triangles are the power before subtracting foregrounds, diamonds are after 8 MHz mean subtraction, squares are after 2 MHz mean subtraction, and circles are after 0.5 MHz subtraction. The curved solid line is the theoretical EoR signal from Jelić et al. (2008), and the dashed line is the theoretical EoR signal with a cold absorbing IGM as described in the text.
ror, but a suitable straightforward estimate given the main complications discussed earlier.

Comparison to Models
Fig. 10 can be compared to simulated results from the Low Frequency Array (LOFAR) EoR project in Jelić et al. (2008), which assumes Ts ≫ TCMB. At low l, their simulated EoR signal is approximately (10 mK) 2 , while our lowest point with a similar 0.5 MHz bandwidth filter is (50 mK) 2 with a 2σ upper limit of (70 mK) 2 . These results are comparable to the sensitivities LOFAR expects after 400 hours of the EoR project. We have also considered the case where reheating of the IGM does not occur, so the spin temperature remains coupled to the kinetic temperature of the gas (Ciardi & Madau 2003). In this case, the IGM cools adiabatically after decoupling from the CMB at z ≈ 150. The temperature fluctuations scale with (1 + TCMB/Ts), and the power scales with the same factor squared. Using T k = TCMB(1 + z)/150 at z = 8.6, the power becomes approximately 275 times larger. This line is shown in Fig. 10, and is comparable to the data. Scaling up the warm IGM power spectrum from Iliev et al. (2008) or Jelić et al. (2008) in this way is a reasonable approximation to what the expected signal in a cold IGM. For a more detailed studies of the signal in such an absorption regime, we point the reader to Baek et al. (2009) and Baek et al. (2010). The power spectrum of reionization is intrinsically three dimensional (Morales & Hewitt 2004). The strongest constraints on the 3D power spectrum for the ∆ν = 2 and 0.5 MHz foreground filter case are shown in Fig. 11. This uses k 2 = k 2 + k 2 ⊥ , where k is given by the windowing func- . 3D power spectrum for the same data shown in Fig. 10, using k 2 = k 2 + k 2 ⊥ . This is dominated by k , so the bin width in k ⊥ does not influence the horizontal position of the limits. The strongest constraints from the 2 MHz and 0.5 MHz filters are shown (square and circle, respectively). Upper limits are 2σ bootstrap errors. Three possible signals are shown. The dashed line is the prediction from Iliev et al. (2008) and the double-dotdashed line is the same for a cold IGM. The solid line comes from the single-scale bubble model as described in the text for a cold IGM, using k = 2.5/R to show the maximum power at all k. For the two points shown, the bubble diametres which achieve this maximum power are 27 and 7.4 h −1 Mpc respectively. Only the 0.5 MHz point imposes a limit on the diametre. For a warm IGM case, this signal would be reduced by the same factor as in the two dashed lines. tion of the filter and k ⊥ = l/6 h −1 Gpc. When comparing this to the prediction from Iliev et al. (2008), one should note that our windowing function will also reduce the predicted signal by at most a factor of two.
We also consider an idealized case in which the ionized bubbles during reionization are of uniform scale and nonoverlapping. Then for a given k there will be a characteristic bubble radius R at which the power is maximized. By taking the 3D Fourier transform of a perfectly ionized bubble, and requiring that the universe is 50 per cent ionized, one can show the power is given by and is maximized when kR ≈ 2.5. In this case, k 3 /(2π 2 )P (k) ≈ T 2 b /5, where T b is the brightness temperature ≈ 30 mK in an X-ray heated IGM or almost -500 mK in a cold absorbing IGM. This signal would be more than an order of magnitude larger than the predictions by e.g. Iliev et al. (2008) or Jelić et al. (2008). In Fig. 11 we have included the power spectrum from this model with k chosen to maximize the power in the range of interest. The data currently impose a limit on the size of bubbles in this single-scale model. Our upper limits with a 0.5 MHz foreground subtraction rule out bubbles with diametres from 2.2 to 12.4 h −1 Mpc in the redshift range 8.1 < z < 9.2. The cold IGM constraint is applicable even in the case of simulations, since only UV photons are included which themselves do not heat the IGM.

CONCLUSION
The data analysis has been completed on six days from December 2007 with a noise level of approximately 2 mJy on most nights. The SVD removal strategies for broadband RFI used lower noise by a factor of 4 in temperature, or 16 in power, which flagging alone cannot achieve. We have also tested for ionospheric variations and found that our pulsar calibration is sufficient for dealing with these effects. After RFI removal and foreground subtraction, we have measured a power spectrum which represents a new upper limit on the 21 cm brightness temperature fluctuations during the epoch of reionization. These results can be used to constrain assumptions about the state of the IGM at these times, particularly in the case of a Ly-α pumped, but cold, IGM.
The previous best limit on 21 cm signal at comparable redshift was by Bebbington (1986), who reported no features down to 5 K at z = 8.4. Parsons et al. (2010) have reported a similar limit of about 5 K using PAPER. The upper limit we present here is approximately 70 mK on the variance in 21 cm brightness temperature at z = 8.6, almost two orders of magnitude better than these previous limits. Residual foreground contamination and RFI may still be contributing to this power, but the EoR signal can not be larger than this.
The analysis of additional observations of the B0823+26 field made since 2007, containing approximately 2 times more data than has been treated in this work, is ongoing. We expect to continue to improve these results, and the planned addition of a second field will strengthen this statistical measurement. The dominant uncertainty source remains RFI and foreground modelling.
These results provide a first-look at the progress made at GMRT in detecting EoR. The GMRT-EoR team has been continuing observations, most recently with an additional 150 hours allocated in the summer of 2010. We continue to improve both the system temperature of GMRT antennas and the RFI monitoring and removal strategies and expect to improve upon these results as the analysis of the newer data progresses.

ACKNOWLEDGMENTS
We thank Chris Hirata for his contributions to the analysis pipeline and a reading of the paper, and the anonymous reviewer for many useful comments. We also thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. The computations were performed on CITA's Sunnyvale clusters which are funded by the Canada Foundation for Innovation, the Ontario Innovation Trust, and the Ontario Research Fund. The work of UP is supported by the National Science and Engineering Research Council of Canada. The map image in Fig. 4 is copyright by OpenStreetMap 7 contributors and used under license 8 .