MEASUREMENT OF 21 cm BRIGHTNESS FLUCTUATIONS AT z ∼ 0.8 IN CROSS-CORRELATION

Measurements of neutral hydrogen are essential to our understanding of the universe. Following cosmological reionization at z ∼ 6, the majority of hydrogen outside of galaxies is ionized. Within galaxies, it must pass through its neutral phase (H i) as it cools and collapses to form stars. The quantity and distribution of neutral hydrogen is therefore intimately connected with the evolution of stars and galaxies, and observations of neutral hydrogen can give insight into these processes.

Above redshift z = 2.2, the Lyα line redshifts into optical wavelengths and H i can be observed, typically in absorption against distant quasars (Prochaska & Wolfe 2009). Below redshift z = 0.1, H i has been studied using 21 cm emission from its hyperfine splitting (Zwaan et al. 2005; Martin et al. 2010). There, the abundance and large-scale distribution of neutral hydrogen are inferred from large catalogs of discrete galactic emitters. Between z = 0.1 and z = 2.2 there are fewer constraints on neutral hydrogen, and those that do exist (Meiring et al. 2011; Lah et al. 2007; Rao et al. 2006) have large uncertainties.

While the 21 cm line is too faint to observe individual galaxies in this redshift range, one can nonetheless pursue three-dimensional (3D) intensity mapping (Chang et al. 2008; Loeb & Wyithe 2008; Ansari et al. 2012; Mao et al. 2008; Seo et al. 2010; Mao 2012). Instead of cataloging many individual galaxies, one can study the large-scale structure (LSS) directly by detecting the aggregate emission from many galaxies that occupy large ∼1000 Mpc³ voxels. The use of such large voxels allows telescopes such as the Green Bank Telescope (GBT) to reach z ∼ 1, conducting a rapid survey of a large volume.

Aside from being used to measure the hydrogen content of galaxies, intensity mapping promises to be an efficient way to study the large-scale structure of the universe. In particular, the method could be used to measure the baryon acoustic oscillations to high accuracy and constrain dark energy (Chang et al. 2008). However, intensity mapping is a new technique that is still being pioneered. Ongoing observational efforts such as the one presented here are essential for developing this technique as a powerful probe of cosmology.

Synchrotron foregrounds are the primary challenge to this method because they are three orders of magnitude brighter than the 21 cm signal. However, the physical process of synchrotron emission is known to produce spectrally smooth radiation (Oh & Mack 2003; Seo et al. 2010). If the calibration, spectral response and beam width of the instrument are well controlled and characterized, the subtraction of foregrounds should be possible because the foregrounds have fewer degrees of freedom than the cosmological signal. We find that this allows the foregrounds to be cleaned to the level of the expected signal. The auto-correlation of intensity maps is biased by residual foregrounds, and minimizing and constraining these residuals is an active area of work. However, because residual foregrounds should be uncorrelated with the cosmological signal, they only boost the noise in a cross-correlation with existing surveys. This makes the cross-correlation a robust indication of neutral hydrogen density fluctuations in the 21 cm intensity maps (Chang et al. 2010; Vujanovic et al. 2012).

The first detection of the cross-correlation between LSS and 21 cm intensity maps at z ∼ 1 was reported in Chang et al. (2010), based on data from GBT and the DEEP2 galaxy survey. Here we improve on these measurements by cross correlating new intensity mapping data with the WiggleZ Dark Energy Survey (Drinkwater et al. 2010). Our measurement improves on the statistical precision and range of scales of the previous result, which was based on 15 hr of GBT integration time over 2 deg. sq.

Throughout, we use cosmological parameters from Komatsu et al. (2009), in accord with Blake et al. (2011).

The observations presented here were conducted with the 680–920 MHz prime-focus receiver at the GBT. The unblocked aperture of GBT's 100 m offset paraboloid design results in well-controlled sidelobes and ground spill, advantageous to minimizing radio-frequency contamination and overall system temperature (∼25 K). The receiver is sampled from 700 MHz (z = 1) to 900 MHz (z = 0.58) by the Green Bank Ultimate Pulsar Processing Instrument (GUPPI) pulsar back-end systems (DuPlain et al. 2008).

The data used in this analysis were collected between 2011 February and November as part of a 400 hr allocation over four fields. This allocation was specifically to corroborate previous cross-correlation measurements (Chang et al. 2010) over a larger survey area, and to search for auto-power of diffuse 21 cm emission. The analysis here is based on a 105 hr integration of a 45 × 24 "15 hr deep field" centered at $14^{\rm h}31^{\rm m}28\buildrel{\mathrm{s}}\over{.}5$ right ascension, 2°0' declination and an 84 hr integration on a 70 × 43 "1 hr shallow" field centered at 0^h52^m0^s right ascension, 2°9' declination. The beam FWHM at 700 MHz is 0314 and at 900 MHz it is 025. At band-center, the beam width corresponds to a comoving length of 9.6 h⁻¹ Mpc. Both fields have nearly complete angular overlap and good redshift coverage with WiggleZ.

Our observing strategy consists of sets of azimuthal scans at constant elevation to control ground spill. We start the set at the low right ascension (right hand) side of the field and allow the region to drift through. We then re-point the telescope to the right side of the field and repeat the process. For the 15 hr field, this set of scans consists of 8 one-minute scans each with a stroke of 4°. For the 1 hr field, a set of scans consists of 10 two-minute scans, each 8° in length. Note that since we observe over a range of local sidereal times, our scan directions cover a range of angles with respect to the sky. This range of crossing angles makes the noise more isotropic, and allows us to ignore the directional dependence of the noise in the 3D power spectrum. The survey regions have most coverage in the middle due to the largest number of intersecting scans. Observations were conducted at night to minimize radio-frequency interference (RFI).

The optical data are part of the WiggleZ Dark Energy Survey (Drinkwater et al. 2010), a large-scale spectroscopic survey of emission-line galaxies selected from UV and optical imaging. It spans redshifts 0.2 < z < 1.0 across 1000 deg². The selection function (Blake et al. 2010) has angular dependence determined primarily by the UV selection, and redshift coverage which favors the z = 0.6 end of the radio band. The galaxies are binned into volumes with the same pixelization as the radio maps and divided by the selection function, so that we consider the cross-power with respect to optical overdensity.

Here we describe our analysis pipeline, which converts the raw data into 3D intensity maps, then correlates these maps with the WiggleZ galaxies.¹¹

3.1. From Data to Maps

The first stage of our data analysis is a rough cut to mitigate contamination by terrestrial sources of RFI. Our data natively have fine spectral resolution with 4096 channels across 200 MHz of bandwidth. This facilitates the identification and flagging of RFI. In each scan, individual channels are flagged based on their variance. Any RFI not sufficiently prominent to be flagged in this stage is detected as increased noise later in the pipeline and subsequently down-weighted during map-making. Additional RFI is detected as frequency–frequency covariance in the foreground cleaning and subtracted in the map domain. While RFI is prominent in the raw data, after these steps, it was not found to be the primary limitation of our analysis.

In addition to RFI, we also eliminate channels within 6 MHz of the band edges (where aliasing is a concern) and channels in the 800 MHz receiver's two resonances at roughly 798 MHz and 817 MHz. Before mapping, the data are re-binned to 0.78 MHz-wide bands (corresponding to roughly 3.8 h⁻¹ Mpc at band-center).

For a time-transfer calibration standard, we inject power from a noise diode into the antenna. The noise diode raises the system temperature by roughly 2 K and we switch it at 16 Hz so that the noise power can be cleanly isolated. Calibration is performed by first dividing by the noise diode power (averaged over a scan) in each channel, and then converting to flux using dedicated observations of 3C286 and 3C48. The gain for X and Y polarizations may differentially drift and so these are calibrated independently. Our absolute calibration uncertainty is dominated by the calibration of the reference flux scale (5%; Kellermann et al. 1969), measurements of the calibration sources with respect to this reference (5%; see also Scaife & Heald 2012), and uncertainty of our measurement of these fluxes (5%). Receiver nonlinearity, uncertainty in the beam shape, and variations in the diffuse galactic emission in the on- and off-source measurements are estimated to contribute of the order of 1% each. These are all assumed to be uncorrelated errors and give 9% total calibration systematic error.

Gridding the data from the time ordered data to a map is done in two stages. We follow cosmic microwave background (CMB) map-making conventions as described in Tegmark (1997). The map maker treats the noise to be uncorrelated except for deweighting the mean and slope along the time axis for each scan. Each frequency channel is treated independently. In the first round of map-making, the noise is estimated from the variance of the scan. This is inaccurate because the foregrounds dominate the noise. This yields a suboptimal map which nonetheless has high a signal-to-noise ratio on the foregrounds. This map is used to estimate the expected foreground signal in the time ordered data and to subtract this expected signal, leaving time ordered data that are dominated by noise. After flagging anomalous data points at the 4σ level, we re-estimate the noise and use this estimate for a second round of map-making, yielding a map which is much closer to optimal. In reality, it is a bad assumption that the noise is uncorrelated. We have observed correlations at finite time lag and between separate frequency channels in our data. Exploiting these correlations to improve the optimality of our maps is an area of active research. For all map-making, we use square pixels with widths of 00627, which corresponds to a quarter of the beam's FWHM at the high frequency edge of our band. Figure 1 shows the 15 hr field map.

Zoom In Zoom Out Reset image size

In addition to the observed maps, we develop signal-only simulations based on Gaussian realizations of the nonlinear, redshift-space power spectrum using the empirical-non-linear (NL) model described by Blake et al. (2011).

3.2. From Maps to Power Spectra

The approach to 21 cm foreground subtraction in literature has been dominated by the notion of fitting and subtracting smooth, orthogonal polynomials along each line of sight. This is motivated by the eigenvectors of smooth synchrotron foregrounds (Liu & Tegmark 2011, 2012). In practice, instrumental factors such as the spectral calibration (and its stability) and polarization response translate into foregrounds that have more complex structure. One way to quantify this structure is to use the map itself to build the foreground model. To do this, we find the frequency–frequency covariance across the sample of angular pixels in the map, using a noise inverse weight. We then find the principal components along the frequency direction, order these by their singular value, and subtract a fixed number of modes of the largest covariance from each line of sight. Because the foregrounds dominate the real map, they also dominate the largest modes of the covariance.

There is an optimum in the number of foreground modes to remove. For too few modes, the errors are large due to residual foreground variance. For too many modes, 21 cm signal is lost, and so after compensating based on simulated signal loss (see below), the errors increase modestly. We find that removing 20 modes in both the 15 hr and 1 hr field maximizes the signal. Figure 1 shows the foreground-cleaned 15 hr field map.

We estimate the cross-power spectrum using the inverse noise variance of the maps and the WiggleZ selection function as the weight for the radio and optical survey data, respectively. The variance is estimated in the mapping step and represents noise and survey coverage. The foreground cleaning process also removes some 21 cm signal. We compensate for signal loss using a transfer function based on 300 simulations where we add signal simulations to the observed maps (which are dominated by foregrounds), clean the combination, and find the cross-power with the input simulation. Because the foreground subtraction is anisotropic in k_⊥ and k_∥, we estimate and apply this transfer function in two dimensions (2D). The GBT beam acts strictly in k_⊥, and again we develop a 2D beam transfer function using signal simulations with the beam.

The foreground filter is built from the real map which has a limited number of independent angular elements. This causes the transfer function to have components in both the angular and frequency direction (Nityananda 2010), with the angular part dominating. This is accounted for in our transfer function. Subtleties of the cleaning method will be described in a future methods paper.

We estimate the errors and their covariance in our cross-power spectrum by calculating the cross-power of the cleaned GBT maps with 100 random catalogs drawn from the WiggleZ selection function (Blake et al. 2010). The mean of these cross powers is consistent with zero, as expected. The variance accounts for shot noise in the galaxy catalog and variance in the radio map either from real signal (sample variance), residual foregrounds or noise. Estimating the errors in this way requires many independent modes to enter each spectral cross-power bin. This fails at the lowest k values and so these scales are discarded. In going from the 2D power to the 1D powers presented here, we weight each 2D k-cell by the inverse variance of the 2D cross-power across the set of mock galaxy catalogs. The 2D–1D binning weight is multiplied by the square of the beam and foreground cleaning transfer functions. Figure 2 shows the resulting galaxy–H i cross-power spectra.

Zoom In Zoom Out Reset image size

To relate the measured spectra with theory, we start with the mean 21 cm emission brightness temperature (Chang et al. 2010),

$$\begin{equation} T_b=0.29\frac{\Omega _{\rm {H\,{\scriptsize I}}}}{10^{-3}}\left(\frac{\Omega _m+ (1+z)^{-3}\Omega _\Lambda }{0.37}\right)^{-\frac{1}{2}}\left(\frac{1+z}{1.8} \right)^{\frac{1}{2}} {\rm mK}. \end{equation} \tag{ 1 }$$

Here, Ω_H I is the comoving H i density (in units of today's critical density), and Ω_m and Ω_Λ are evaluated at the present epoch. We observe the brightness contrast, $\delta T = T_b \delta _{\rm {H\,{\scriptsize I}}}$, from fluctuations in the local H i overdensity δ_H I. On large scales, it is assumed that neutral hydrogen and optically selected galaxies are biased tracers of the dark matter, so that $\delta _{\rm {H\,{\scriptsize I}}} = b_{\rm {H\,{\scriptsize I}}} \delta$, and δ_opt = b_optδ. In practice, both tracers may contain a stochastic component, so we include a galaxy–H i correlation coefficient r. This quantity is scale-dependent because of the k-dependent ratio of shot noise to large-scale structure, but should approach unity on large scales. The cross-power spectrum is then given by $P_{\rm {H\,{\scriptsize I}},\rm {opt}}(k)=T_b b_{\rm {H\,{\scriptsize I}}}b_{\rm {opt}}rP_{\delta \delta }(k)$ where P_δδ(k) is the matter power spectrum.

The large-scale matter power spectrum is well known from CMB measurements (Komatsu et al. 2011) and the bias of the optical galaxy population is measured to be b²_opt = 1.48 ± 0.08 at the central redshift of our survey (Blake et al. 2011). Simulations including nonlinear scales (as in Section 3.1) are run through the same pipeline as the data. We fit the unknown prefactor Ω_H Ib_H Ir of the theory to the measured cross-powers shown in Figure 2, and determine Ω_H Ib_H Ir = [0.44 ± 0.10(stat.) ± 0.04(sys.)] × 10⁻³ for the 15 hr field data, and Ω_H Ib_H Ir = [0.41 ± 0.11(stat.) ± 0.04(sys.)] × 10⁻³ for the 1 hr field data. The systematic term represents the 9% absolute calibration uncertainty from Section 3.1. It does not include current uncertainties in the cosmological parameters or in the WiggleZ bias, but these are subdominant. Combining the two fields yields Ω_H Ib_H Ir = [0.43 ±  0.07(stat.) ±  0.04(sys.)] ×  10⁻³. These fits are based on the range 0.075 h Mpc⁻¹ < k < 0.3 h Mpc⁻¹ over which we believe that errors are well estimated (failing toward larger scales where there are too few k modes in the volume) and under the assumption that nonlinearities and the beam/pixelization (failing toward smaller scales) are well understood. A less conservative approach is to fit for 0.05 h Mpc⁻¹ < k < 0.8 h Mpc⁻¹ where the beam, model of nonlinearity and error estimates are less robust, but which shows the full statistical power of the measurement, at 7.4σ combined. Here, Ω_H Ib_H Ir = [0.40 ± 0.05(stat.) ± 0.04(sys.)] × 10⁻³ for the combined, Ω_H Ib_H Ir = [0.46 ± 0.08] × 10⁻³ for the 15 hr field and Ω_H Ib_H Ir = [0.34 ± 0.07] × 10⁻³ for the 1 hr field.

To compare with the result in Chang et al. (2010), Ω_H Ib_relr = [0.55 ± 0.15(stat.)] × 10⁻³, we must multiply their relative bias (between the GBT intensity map and DEEP2) by the DEEP2 bias b = 1.2 (Coil et al. 2004) to obtain an expression with respect to b_H I. This becomes Ω_H Ib_H Ir = [0.66 ± 0.18(stat.)] × 10⁻³, and is consistent with our result.

The absolute abundance and clustering of H i are of great interest in studies of galaxy and star formation. Our measurement is an integral constraint on the H i luminosity function, which can be directly compared to simulations. The quantity Ω_H Ib_H I also determines the amplitude of 21 cm temperature fluctuations. This is required for forecasts of the sensitivity of future 21 cm intensity mapping experiments. Since r < 1 we have put a lower limit on Ω_H Ib_H I.

To determine Ω_H I alone from our cross-correlation requires external estimates of the H i bias and stochasticity. The linear bias of H i is expected to be ∼0.65 to ∼1 at these redshifts (Marín et al. 2010; Khandai et al. 2011). Simulations to interpret Chang et al. (2010) find values for r between 0.9 and 0.95 (Khandai et al. 2011), albeit for a different optical galaxy population. Measurements of the correlation coefficient between WiggleZ galaxies and the total matter field are consistent with unity in this k-range (with r_{m, opt} ≳ 0.8) (Blake et al. 2011). These suggest that our cross-correlation can be interpreted as Ω_H I between 0.45 × 10⁻³ and 0.75 × 10⁻³.

Measurements with the Sloan Digital Sky Survey (Prochaska & Wolfe 2009) suggest that before z = 2, Ω_H I may have already reached ∼0.4 × 10⁻³. At low redshift, 21 cm measurements give Ω_H I(z ∼ 0) = (0.43 ± 0.03) × 10⁻³ (Martin et al. 2010). Intermediate redshifts are more difficult to measure, and estimates based on Mg ii lines in DLA systems observed with the Hubble Space Telescope find Ω_H I(z ∼ 1) ≈ (0.97 ±  0.36) × 10⁻³ (Rao et al. 2006), in rough agreement with z ≈ 0.2 DLA measurements (Meiring et al. 2011) and 21 cm stacking (Lah et al. 2007). This is in some tension with a model where Ω_H I falls monotonically from the era of maximum star formation rate (Duffy et al. 2012). Under the assumption that b_H I = 0.8, r = 1, the cross-correlation measurement here suggests Ω_H I ∼ 0.5 × 10⁻³, in better agreement, but clearly better measurements of b_H I and r are needed. Redshift space distortions can be exploited to break the degeneracy between Ω_H I and bias to measure these quantities independently of simulations (Wyithe 2008; Masui et al. 2010). This will be the subject of future work.

Our measurement is limited by both the number of galaxies in the WiggleZ fields and by the noise in our radio observations. Simulations indicate that the variance observed in our radio maps after foreground subtraction is roughly consistent with the expected levels from thermal noise. This is perhaps not surprising, our survey being relatively wide and shallow compared to an optimal LSS survey, however, this is nonetheless encouraging.

We thank John Ford, Anish Roshi, and the rest of the GBT staff for their support; Paul Demorest and Willem van-Straten for help with pulsar instruments and calibration; and K. Vanderlinde for helpful conversations.

K.W.M. is supported by NSERC Canada. E.R.S. acknowledges support by NSF Physics Frontier Center grant PHY-0114422 to the Kavli Institute of Cosmological Physics. J.B.P. and T.C.V. acknowledge support under NSF grant AST-1009615. X.C. acknowledges the Ministry of Science and Technology Project 863 (under grant 2012AA121701), the John Templeton Foundation and NAOC Beyond the Horizon program, and the NSFC grant 11073024. A.N. acknowledges financial support from the Bruce and Astrid McWilliams Center for Cosmology.

Computations were performed on the GPC supercomputer at the SciNet HPC Consortium.