Detection of Gravitational Lensing in the Cosmic Microwave Background

Detection of Gravitational Lensing in the Cosmic Microwave Background


Gravitational lensing of the cosmic microwave background (CMB), a long-standing prediction of the standard cosmolgical model, is ultimately expected to be an important source of cosmological information, but first detection has not been achieved to date. We report a 3.4 detection, by applying quadratic estimator techniques to all sky maps from the Wilkinson Microwave Anisotropy Probe (WMAP) satellite, and correlating the result with radio galaxy counts from the NRAO VLA Sky Survey (NVSS). We present our methodology including a detailed discussion of potential contaminants. Our error estimates include systematic uncertainties from density gradients in NVSS, beam effects in WMAP, Galactic microwave foregrounds, resolved and unresolved CMB point sources, and the thermal Sunyaev-Zeldovich effect.

I Introduction

Within just two decades, cosmology has progressed from a rather speculative science to one of the most successful fields of physics, driven by an exemplary interplay between experiment and theory. Much of this progress has been owing to the well understood physics underlying the Comic Microwave Background (CMB) anisotropy, seeded by oscillations in the baryon-photon plasma of the early universe.

Measurements of these fluctuations by a number of experiments have given rise to a basic cosmological paradigm, with the tightest current constraints on the cosmological parameter budget coming from combinations of data from the Wilkinson Microwave Anisotropy Probe (WMAP) satellite Bennett et al. (2003a); Spergel et al. (2006) in conjunction with small scale CMB experiments (e.g. Jones et al. (2006); Readhead et al. (2004a); Kuo et al. (2006)), and other rich probes of cosmological clustering and dynamics such as supernovae, galaxy surveys, the Lyman-alpha forest, weak lensing, and others (e.g. Riess et al. (1998); Perlmutter et al. (1999); York et al. (2000); Colless et al. (2001); Seljak et al. (2006); Hoekstra et al. (2002); Astier et al. (2006); Refregier (2003); Tegmark et al. (2006)).

The CMB promises to remain a gold mine for precision cosmology, and two new frontiers lie ahead. First, a polarized component has recently been detected by a number of groups Kovac et al. (2002); Piacentini et al. (2006); Readhead et al. (2004b); Barkats et al. (2005); Montroy et al. (2006); Johnson et al. (2006), offering e.g. the prospects of detecting primordial gravitational waves and constraining recombination physics.

Second, large scale structure between the last scattering surface and us alters the primary CMB anisotropy, through gravitational lensing (for a recent review of the theory see Lewis and Challinor (2006)), through scattering off hot electrons in large scale structure (the Sunyaev-Zel’dovich effects) Sunyaev and Zeldovich (1970, 1980), and through redshifting during the traverse of time-dependent potential fluctuations (the ISW effect) Sachs and Wolfe (1967). A number of specialized instruments will soon begin to study details of these secondary anisotropies Kosowsky (2003); Ruhl et al. (2004).

As important as constraining cosmological and astrophysical parameters, detecting any of these effects is a crucial milestone for cosmological physics. The Sunyaev-Zel’dovich effect has been found by targeting clusters detected in X-ray Birkinshaw et al. (1984); Holzapfel et al. (1997); Dawson et al. (2002); LaRoque et al. (2002), also at high significance level using WMAP Afshordi et al. (2006), and it has been observed in cross-correlation of galaxy surveys with WMAP Hernandez-Monteagudo et al. (2004); Afshordi et al. (2004); Fosalba et al. (2003). The ISW effect has been detected in cross-correlation of WMAP with galaxy surveys and with the hard X-ray background Boughn and Crittenden (2004, 2005); Nolta et al. (2004); Maddox et al. (1990); Fosalba et al. (2003); Scranton et al. (2003); Fosalba and Gaztanaga (2004); Padmanabhan et al. (2005); Cabre et al. (2006); McEwen et al. (2006).

A detection of gravitational lensing in the CMB has so far been outstanding. The main difficulty at millimeter wavelengths is the high angular resolution needed, as typical deflection angles over a cosmological volume are only a few arcminutes. Non-Gaussianity imprinted by lensing into the primordial CMB may allow statistical detection with surveys at lower angular resolution, but the signal-to-noise is currently too low for internal detections. Cross correlation with other tracers of large scale structure offers a way to limit systematics and increase the signal to noise.

A first attempt Hirata et al. (2004) was made by cross correlating the WMAP first year release Bennett et al. (2003a) with data from the Sloan Digital Sky Survey (SDSS) York et al. (2000). These authors used a sample of 503,944 SDSS Luminous Red Galaxies (LRG’s) overlapping with of the sky observed by WMAP. They were not able to find evidence for gravitational lensing within statistical error bounds. While SDSS LRG’s have a well understood redshift distribution, their number density drops rapidly beyond , and has only marginal overlap with the higher redshift range that is geometrically optimal for CMB lensing. Photometric quasars found in SDSS may offer an additional handle.

Here we go a different route, using the 1.9 million radio sources found in the NRAO VLA Sky Survey (NVSS) Condon et al. (1998). The large sky coverage and estimated depth of NVSS make it an excellent candidate for a search for CMB lensing in cross correlation with WMAP. The survey covers 77% of the sky, 58% of which is found to overlap with WMAP, once masks to limit systematics have been applied.

The structure of this paper is as follows.

First, we describe the datasets (§II), theory (§III), and pipeline (§IV) that will be used for detecting CMB lensing by reconstructing the lensing potential from WMAP, and cross-correlating the result to NVSS. The detection is shown, with statistical errors only, in Fig. 5. The rest of the paper is devoted to null tests and assigning systematic errors: NVSS systematics (§V), WMAP beam effects and Galactic foregrounds (§VI), resolved and unresolved point sources (§VII), and Sunyaev-Zeldovich fluctuations (§VIII). We quote our final result including systematic errors (Fig. 19) in §IX, where we also mention future directions.

In our calculations we will assume throughout the cosmological model favored by a combination of WMAP, smaller scale CMB experiments, and other data (the WMAP+ALL analysis, Spergel et al. (2006)): a local expansion rate km/s/Mpc, primordial power spectrum slope , matter and dark energy fractions of and respectively, and amplitude .

Ii Datasets

ii.1 Wmap

Figure 1: Left panel: CMB signal power spectrum, and three-year WMAP noise power spectrum. Right panel: Fiducial NVSS signal power spectrum, and NVSS shot noise ( gal/steradian).

With the goal of producing full sky maps of the CMB with unprecedented accuracy, the WMAP satellite was launched in June 2001. Since then, it has been mapping the sky using 10 Differential Assemblies covering 5 frequency bands centered at 23 (K), 33 (Ka), 41 (Q), 61 (V) and 94 GHz (W). In our analysis we use the 2 Q-band, 2 V-band, and 4 W-band temperature maps produced using 3 years of observations Hinshaw et al. (2006) and made publicly available 1. We will use as a default mask the Kp0 mask, which cuts out the Galactic plane and point sources bright enough to be resolved by WMAP, leaving about 78.46% of the sky Bennett et al. (2003a).

The intrinsic quality of this dataset leaves us with few instrumental systematic effects to worry about Jarosik et al. (2003); Page et al. (2003); Jarosik et al. (2006). Nonetheless, noise inhomogeneities and beam effects could be of particular concern for our lensing statistic. The former will be optimally handled by our estimator. Although the latter are well controlled for the power spectrum estimation Jarosik et al. (2006); Hinshaw et al. (2006), they could potentially affect our lensing estimator as will be discussed below. We will show how the formalism presented in Hinshaw et al. (2006) allows us to control them in our particular context too. Another source of systematic error might come from other astrophysical sources, namely residual galactic foregrounds (synchrotron, free-free and dust), residual point sources and the signature of galaxy clusters via the Sunyaev-Zeldovich effect. These potential contaminants will be discussed in later sections.

ii.2 Nvss

As a tracer of the large scale density field, we use observations resulting from the NRAO VLA Sky Survey (NVSS) 1.4 GHz continuum survey. This survey covers 82% of the sky with Condon et al. (1998) with a source catalog containing over sources that is 50% complete at 2.5 mJy. It is appropriate for our purpose since most of the bright sources are AGN-powered radio galaxies and quasars whereas the less bright ones correspond to nearby star-forming galaxies. As a consequence, almost all the sources away from the Galactic plane () are extragalactic.

We pixelize the NVSS catalog using HEALPix Gorski et al. (2005) maps with corresponding to around  square pixels 2. As an extra precaution, we removed sources with a flux greater than 1 Jy as well as a 1 degree disk around them. We also mask out pixels at low Galactic latitude () and those unobserved by the survey (). We ended up with sources with an average density gal/steradian.

Iii CMB lensing

Weak lensing by large scale structure remaps the CMB temperature field on the sky; the lensed temperature and unlensed temperature are related by Blanchard and J. (1987)


where is a vector field representing the deflection angles. To first order in perturbation theory, is expected to be a pure gradient:


where the scalar potential is given by the line of sight integral:


where denotes conformal distance along the line of sight in the assumed flat cosmology, is the conformal distance to recombination, and is conformal time today. The integral in Eq. (3) receives contributions from a broad redshift range with median around .

How can CMB lensing be detected in data? At the power spectrum level, lensing slightly smooths the acoustic peaks in the temperature power spectrum and adds power in the damping tail Seljak (1996). However, these effects are too small to be detectable in existing datasets. Going beyond the power spectrum, the effect of CMB lensing on higher-point statistics of the CMB is stronger and requires less instrumental sensitivity to detect Bernardeau (1997).

The theory of CMB lens reconstruction Hu (2001a, b); Okamoto and Hu (2003); Hirata and Seljak (2003) provides a framework for extracting this higher-point signal which we will use throughout this paper. One first defines a quadratic (in the CMB temperature ) estimator for the CMB lensing potential . The simplest higher-point estimator for detecting CMB lensing would be the power spectrum : a quadratic estimator in the reconstruction or a four-point estimator in .

Figure 2: Left panel: Auto power spectrum of the CMB lensing potential, and reconstruction noise power spectrum (Eq. (4)) at three-year WMAP noise levels. Right panel: Cross power spectrum between the CMB lensing potential and NVSS galaxy counts, and the effective noise power spectrum for detecting the cross-correlation. The “boost” in signal-to-noise between the two cases is sufficient to obtain a several-sigma detection of CMB lensing.

However, the three-year WMAP data do not have sufficient sensitivity to detect CMB lensing via the auto power spectrum . This can be seen by considering the statistical “noise” in the reconstruction; in Hu (2001a) it is shown that the reconstruction noise power spectrum is given by:


where is defined by


In Fig. 2 (left panel) we have shown the noise power spectrum for three-year WMAP sensitivity, with the fiducial signal power spectrum shown for comparison. Although the CMB temperature anisotropies are signal-dominated across a wide range of angular scales (Fig. 1), the lens reconstruction is highly noise-dominated. At this level of signal-to-noise, an “internal” (to WMAP) detection of CMB lensing, by measuring the auto power spectrum , is not possible.

It is frequently the case that a signal which is too noisy for internal detection can nonetheless be detected via cross-correlation to a second, less noisy signal. (For example, the first-year WMAP data had poor sensitivity to the polarization signal, but contained a many-sigma detection of CMB polarization via the cross-correlation Kogut et al. (2003)). In this paper, we will detect the lensing signal in WMAP by cross-correlating to radio galaxy counts in NVSS, thus detecting a nonzero cross power spectrum . The galaxy field is much less noisy than (Fig. 1), but the two fields have a significant redshift range in common and so are highly correlated; the correlation in the fiducial model is on angular scales . Therefore, the effective signal-to-noise is higher for the cross-correlation (Fig. 2, right panel). A forecast based on this signal-to-noise ratio, and the assumption of simple scaling, predicts that a sigma detection can be made. If the same forecast is repeated using the parameters from Hirata et al. (2004) (i.e. first-year WMAP sensitivity and Sloan LRG’s over 4000 deg), we find a sigma result, in agreement with previous results.

In addition to the improved statistical errors from higher signal-to-noise, obtaining the detection as a cross-correlation is more robust to systematics, as we will see in detail in §VVIII. Any source of systematic contamination which appears in either WMAP or NVSS, but not both, will not bias our estimates for the cross power spectrum , since it does not correlate the two surveys. At worst, such a contaminant can affect the statistical significance of the detection, by increasing the error bars on each bandpower.

Our estimator for will be defined by cross-correlating the quadratic reconstruction of the lensing potential to the NVSS overdensity field . Thus the estimator is three-point: two-point in the CMB temperature and one-point in the galaxy field. The same three-point estimator can also be derived from the general theory of bispectrum estimation Komatsu et al. (2005); Creminelli et al. (2006); Smith and Zaldarriaga (2006).

The most general three-point correlation between two CMB multipoles and one galaxy multipole which is allowed by rotational and parity invariance is of the form:


This equation defines the bispectrum . (More properly, with the prefactor included, we have defined the “reduced bispectrum” in Eq. (8); with this prefactor reduces to the flat sky bispectrum in the limit of large Komatsu and Spergel (2001). Whenever we write bispectra in this paper, are understood to denote CMB multipoles and denotes a galaxy multipole.

From this perspective, the CMB lensing signal simply gives a contribution to the bispectrum which we want to measure. The lensing bispectrum is proportional to :


One can think of this as a single bispectrum which is estimated to give an overall detection, or a linear combination of independent bispectra corresponding to bandpowers in .

In Appendix B, we show that the lens reconstruction and bispectrum formalisms are equivalent, so that it is a matter of convenience which to use. In this paper, we have generally used the lens reconstruction formalism, but will occasionally refer to the bispectrum formalism when it provides additional perspective.

Figure 3: Mean contribution to the squared total detection significance per NVSS galaxy multipole (left panel), and per unit increase in maximum CMB multipole (right panel). Most of the statistical weight comes from galaxy multipoles near , and CMB multipoles near .

One issue which is clearer from the bispectrum perspective is the distribution of statistical weight. Suppose we consider the total squared detection significance , rather than splitting the signal into bandpowers. Starting from the bispectrum in Eq. (9), one can write as a sum over multipoles . In Fig. 3, we have split up this sum to show the contribution per multipole. (Since there are two CMB multipoles, we show the contribution per unit increase in the maximum multipole .) It is seen that the greatest statistical weight comes from galaxy multipoles near , and CMB multipoles near corresponding to an acoustic trough in the primary CMB. In bispectrum language, most of the signal is in “squeezed” triangles where the galaxy wavenumber is much smaller than the two CMB wavenumbers. This corresponds to the intuitive statement that lens reconstruction estimates degree-scale lenses indirectly through their effect on smaller-scale hot and cold spots in the CMB.

Iv Pipeline

(1)Gaussian fields {, , }(2)(4)Lensed CMB (3)NVSS data(7)Filtered galaxyfield WMAP data(5)Filtered CMB (6)Reconstructedpotential (8)Lensing estimator
Figure 4: Simulation + analysis pipeline used in this paper; the stages (1)-(7) are described in detail in §IV.

In this section, we describe our simulation and analysis pipeline for estimating the cross power spectrum from the WMAP and NVSS datasets, and present results with statistical errors. (Systematics will be treated in §VVIII.)

iv.1 Pipeline description

Our pipeline is shown in Fig. 4. Steps (1)-(4) represent the simulation direction and produce simulated WMAP and NVSS datasets with CMB lensing. Steps (5)-(8) are the analysis direction and produce power spectrum estimates in bands , starting from the WMAP and NVSS datasets. We now describe each step in detail.

The first step (1) is simulating Gaussian fields: the unlensed CMB temperature, lensing potential, and (shot noise free) radio galaxy field . We use the power spectra in the fiducial model. The last two are computed using the Limber approximation (e.g. Bartelmann and Schneider (2001)) and a simple constant galaxy bias model: we take the galaxy overdensity to be given by the line of sight integral


using a fiducial redshift distribution and galaxy bias that will be discussed in the next section.

In step (2), we compute the lensed CMB from the lensing potential and unlensed CMB. The lensing operation


is performed directly in position space (rather than relying on an approximation to Eq. (11) such as the gradient approximation). The right-hand side of Eq. (11) is evaluated using cubic interpolation on a high resolution ( arcmin) map.

In step (3), we simulate the eight Q, V, and W-band channels of WMAP. The maps are simulated at Healpix resolution and downgraded to to minimize pixelization artifacts. To simulate each map, we first convolve with the beam and pixel window in harmonic space:


where is the beam transfer function (distinct for each channel) and is the pixel window function. We then take the spherical transform and add Gaussian noise to each pixel. The noise RMS is pixel-dependent but the noise is assumed uncorrelated between pixels.

As the last step in the simulation direction, in step (4) we simulate NVSS, including clustering which is controlled by the Gaussian field , by generating a Poisson galaxy count in each pixel whose mean is given by


where is the mean number of galaxies per pixel over the survey. We simulate NVSS at and downgrade to .

Step (5) is the first step in the analysis direction: we start with the pixel-space maps corresponding to the eight Q, V, and W-band WMAP channels, and compute a single harmonic-space map representing the inverse signal + noise filtered temperature . This reduction step is a common ingredient in many types of optimal estimators Bond et al. (1998); Oh et al. (1999); Creminelli et al. (2006); Smith and Zaldarriaga (2006); Jewell et al. (2004); Wandelt et al. (2004). The general principle is that the filtering operation completely incorporates the sky cut and noise model, so that optimal estimators can be constructed by simple subsequent operations directly in harmonic space. For example, the optimal power spectrum estimator is obtained by straightforwardly computing the power spectrum .

Here and throughout the body of the paper, we will defer technical details of the estimators to Appendices AB and concentrate on conveying intuition. In this case, the idea is that the filter simply weights each mode of the data by the inverse of its total variance, so that poorly measured modes are filtered out. For example, the sky cut is incorporated into the noise covariance by assigning infinite noise variance to pixels which are masked (in implementation, we use rather than and set the relevant matrix entries to zero). Data outside the sky cut is then completely filtered out: the map is independent of the map values in masked pixels, and everything “downstream” in the analysis pipeline will be blind to the masked data. As a similar example, we marginalize the CMB monopole and dipole modes by assigning them infinite variance. Finally, the beam transfer functions (Eq. (12)) are kept distinct in the filtering operation, so that optimal frequency weighting is performed: the filtered map will receive contributions from all frequencies at low , but will depend mainly on the highest-frequency channels (i.e., the channels with narrow beams) at high . The filtered map can also be thought of as the least-squares estimate of the signal, given data from all channels.

In step (6), we perform lens reconstruction. Given the filtered CMB temperature from step (5), we compute the reconstructed potential , defined by the equation:


where and are defined by


As explained in Hu (2001a), is a noisy reconstruction of the CMB lensing potential (or more precisely, the inverse noise weighted potential , where is the noise covariance of the reconstruction) which is quadratic in the CMB temperature. Note that both and are defined in harmonic space, but Eq. (14) involves multiplication and derivative operations in real space; in Appendix B, we explain in detail how is computed.

In step (7), we perform inverse signal + noise filtering on the NVSS data: given pixel-space galaxy counts, we compute the harmonic-space map where the noise covariance represents shot noise. This is analagous to the WMAP filtering operation in step (5), but there is one new ingredient. In addition to marginalizing data outside the sky cut, and the monopole and dipole, we marginalize any mode which is independent of the angular coordinate in equatorial coordinates. (In harmonic space, this is equivalent to marginalizing modes with .) This is needed to remove a systematic effect in NVSS which we will discuss in detail in §V; for now we remark in advance that all results in this paper include this marginalization.

Finally, in step (8), we compute the bandpower estimator by cross-correlating the fields and from steps (6) and (7). There is one wrinkle here: as we show in Appendix B, to obtain the optimal estimator, we must include an extra term which subtracts the Monte Carlo average taken over unlensed simulations of WMAP:


where is a normalization constant to be discussed shortly. (We have included the factor since we estimate bandpowers assuming that is flat in each band.) Note that the Monte Carlo average vanishes for symmetry reasons in the case of full sky coverage and isotropic noise, but sky cuts or noise inhomogeneities will give rise to a nonzero average. The extra term in Eq. (17) simply improves the variance of the estimator by subtracting the spurious cross-correlation between this average and the galaxy field .

We determine the estimator normalization by end-to-end Monte Carlo simulations of the pipeline, including a nonzero in the simulations for calibration. (Strictly speaking, the normalization should be a matrix which couples bands , but we have neglected the off-diagonal terms, which are small for our case of large sky coverage and wide bands.) As we will see in Appendix B, the normalization is proportional to a cut-sky Fisher matrix element, which must be computed by Monte Carlo unless an approximation is made such as simple scaling. In addition, Monte Carlo simulations are also needed to compute the one-point term in Eq. (17).

This concludes our description of the pipeline. We have not motivated the details in the construction of our lensing estimator , but in Appendix B we show that the estimator is optimal, by proving that it achieves statistical lower limits on the estimator variance, so that the best possible power spectrum uncertainties are obtained. This justifies the combination of ingredients presented here: inverse signal + noise filtering (steps 5 and 7), keeping the lensing potential in harmonic space (step 6), and including the one-point term in the cross-correlation (step 8); and shows that no further improvements are possible.

iv.2 Results

Figure 5: Detection of CMB lensing via the cross power spectrum between the reconstructed potential and galaxy counts. The three 1 error bars on each bandpower represent different Monte Carlo methods: WMAP simulations vs NVSS simulations (left/black), WMAP data vs NVSS simulations (middle/blue), and WMAP simulations vs NVSS data (right/red). These error bars represent statistical errors only; the result with systematic errors included will be shown in Fig. 19.

The result of applying this analysis pipeline to the WMAP and NVSS datasets is shown in Fig. 5. We emphasize that the uncertainties are purely statistical. Systematic errors will be studied in §VVIII, and an updated version of the result shown in §IX, where we also show that the detection significance with systematic errors included is 3.4.

Our error bars were obtained by Monte Carlo, cross-correlating simulations of WMAP and NVSS. As a consistency check, Fig. 5 shows that nearly identical error bars are obtained if WMAP simulations are cross-correlated to the real NVSS data, or vice versa. This is an important check; if it failed, then we would know that our simulations were failing to capture a feature of the datasets which contributes significant uncertainty to the lensing estimator. In addition, it shows that the uncertainties only depend on correctness of one of the simulation pipelines. Suppose, for example, that the NVSS dataset contains unknown catastrophic systematics which invalidate our simulations. Because the same result is obtained by treating NVSS as a black box to be cross-correlated to WMAP simulations, it is still valid (provided that WMAP contains no “catastrophic” systematics!)

Figure 6: CMB lensing detection obtained by analyzing Q-band (left/black error bar in each triple), V-band (middle/blue), and W-band (right/red) data from WMAP separately, showing consistency of the result between CMB frequencies.

As another consistency check, in Fig. 6 we show the detection that is obtained if each frequency in WMAP is analyzed separately. No signs of inconsistency are seen, although we have not attempted to quantify this precisely: the results obtained from different frequencies are correlated even though the CMB noise realizations are independent, because NVSS is identical and so is the underlying CMB realization. For the same reason, we caution the reader that the three sets of error bars in Fig. 6 cannot be combined in a straightforward way to obtain an overall result. The best possible way of combining the data is already shown in Fig. 5: the maps from the three frequencies are combined into a single CMB map which is cross-correlated to NVSS.

iv.3 Curl null test

Our lensing estimator detects a gradient component in the deflection field via cross-correlation to radio galaxy counts. If we instead decompose the deflection field into gradient and curl:


then one can similarly devise an estimator to detect the curl component. Since the curl component is expected to be absent cosmologically, this is a null test Cooray et al. (2005). Note that we have parameterized the curl component by a pseudoscalar potential , for notational uniformity with the gradient component which is parameterized by its scalar potential .

In Appendix B, we show that the optimal estimator is constructed as follows. First, we define a reconstructed potential which is quadratic in the CMB temperature:


with as in Eqs. (15), (16). Second, we define a power spectrum estimator by cross-correlating to galaxy counts, subtracting the one-point term:


This construction is identical to our construction (Eqs. (14), (17))) of the lensing estimator , except that a rotation has been included (via the antisymmetric tensor ) in Eq. (19).

The result of the curl null test is shown in Fig. 7. The for the null test is 12.1 with 8 degrees of freedom, so the null test passes.

Figure 7: Result of the curl null test (). As in Fig. 5, the three error bars on each bandpower represent different Monte Carlo methods: WMAP simulations vs NVSS simulations (left/black), WMAP data vs NVSS simulations (middle/blue), and WMAP simulations vs NVSS data (right/red).

How strong is the null test obtained by demanding that be consistent with zero? One might hope that astrophysical contaminants, such as point sources or the Sunyaev-Zeldovich effect, would contribute both gradient and curl components to the reconstructed deflections, and thus be monitored by the null test. However, parity invariance requires even when . Since astrophysical contaminants are expected to obey parity invariant statistics, they will not bias on average. Our null test therefore only monitors contaminants which can violate parity invariance, such as Galactic foregrounds or instrumental systematics. This is analgous to the null test in CMB polarization experiments: it is not sensitive to all sources of contamination, but is nevertheless an important sanity check.

We remark that for a detection of CMB lensing which is internal to the CMB (detecting lensing via the auto power spectrum , rather than the cross spectrum considered here), one would have one null test () which can monitor parity-invariant contaminants, and one null test () which cannot.

V NVSS systematics

In the previous section, we obtained a statistical detection of CMB lensing (Fig. 5) by cross-correlating WMAP and NVSS, and showed that two consistency checks were satisfied: frequency independence (Fig. 6) and a curl null test (Fig. 7). The rest of the paper is devoted to studying potential instrumental and astrophysical contaminants of the lensing signal, to show that the observed lensing cross-correlation is not due to systematic contamination. In this section, we will consider NVSS systematics.

Figure 8: Maximum likelihood NVSS galaxy power spectrum, calculated without (top panel) and with (bottom panel) marginalization of modes in equatorial coordinates. In the bottom panel, fiducial spectra are shown (both for ) from the model for by Dunlop and Peacock (1990) (dotted line) and our fit in Eq. 21 (dashed line).

If a maximum likelihood galaxy power spectrum is calculated from NVSS using the sky cut described in §II, the power spectrum shown in the top panel of Fig. 8 is obtained. The very high bandpower in the lowest band is a clear sign of systematic contamination. If the low modes are isolated by low-pass filtering the NVSS galaxy counts to , the resulting map shows azimuthal “striping” when plotted in equatorial coordinates (Fig. 9). This is a known systematic effect in NVSS Blake and Wall (2002): due to calibration problems at low flux densities, the galaxy density has a systematic dependence on declination, which can mimic long-wavelength modes in the galaxy field.

Figure 9: NVSS galaxy overdensity field in equatorial coordinates, low-pass filtered to multipoles , showing visible azimuthal striping.

To remove this contaminant, we analyze NVSS in equatorial coordinates, and marginalize any modes in the data which are constant in the azimuthal coordinate . The marginalization is performed by modifying the NVSS noise model so that all such modes are assigned infinite variance, as described in Appendix A. Thus any signal which is constant in is completely filtered out in the inverse signal+noise weighted map which appears in our estimators. Note that treating the marginalization as part of the noise model means that the loss in sensitivity due to marginalizing modes is already included in the statistical errors; it is not necessary to assign systematic errors separately.

Figure 10: Change in maximum likelihood galaxy power spectrum, when NVSS is analyzed with marginalization vs no margnialization (top panel) or marginalization vs marginalization (bottom panel) in equatorial coordinates. The error bars represent the RMS shift obtained when Monte Carlo simulations are analyzed in the same way.

After including this marginalization in the analysis, the NVSS galaxy power spectrum shown in the bottom panel of Fig. 8 is obtained, showing reasonable agreement with our fiducial . Marginalizing modes produces a large shift in the lowest bandpower and a much smaller shift in higher bands. In Fig. 10 (top panel), we show the shift in each bandpower when modes are marginalized, relative to an error bar which shows the RMS shift obtained when the same marginalization is performed in NVSS simulations. It is seen that the shift is statistically significant not only in the lowest band, but all the way to . We conclude that declination gradients in NVSS are an important systematic on a range of scales and should always be marginalized in cosmological studies.

Has marginalizing completely removed the systematic? To answer this, we tried marginalizing the Fourier mode in the azimuthal coordinate , in addition to the mode. In this case, we find (Fig. 10, bottom panel) that the shift in bandpowers is consistent with simulations. (There is a possible glitch at , but this is outside the range of angular scales which contribute to the lensing detection.) Therefore, we believe that marginalizing all modes with in equatorial coordinates completely removes the systematic; there is no evidence that the contamination extends to higher .

In addition to declination gradients, there is another NVSS systematic which has been relevant for cosmological studies: multicomponent sources Blake and Wall (2002); Blake et al. (2004). Radio galaxies whose angular size is sufficiently large to be resolved by the 45-arcsec NVSS beam will appear as multiple objects in the NVSS catalog. This can contribute extra power to the auto spectrum , at a level which is a few percent of the shot noise. At worst, this could increase the variance of our cross-correlation estimator by a few percent without biasing the estimator. Furthermore, as can be seen in Fig. 8 (bottom panel), we see no evidence for galaxy power in excess of fiducial in the highest band, which is most sensitive to this systematic. We conclude that multicomponent sources are a negligible source of systematic error for CMB lensing.

Next we consider uncertainties in the NVSS redshift distribution and galaxy bias . These uncertainties affect our fiducial power spectra in a given cosmology, and would need to be understood in detail if we wanted to constrain cosmological parameters from our lensing detection. However, since we are merely measuring the cross spectrum , there is only one effect to consider: the Monte Carlo error bars we assign depend on the fiducial galaxy spectrum used in the simulations. (We verified in simulations that the fiducial cross spectrum does not significantly affect the error bars.) If we use a fiducial with too little power, we will underestimate our errors. Therefore, it is important to check that our fiducial agrees with the galaxy power spectrum obtained from the data.

Estimates for the radio luminosity function inspired by optical and infrared observations were given in Dunlop and Peacock (1990). Using their mean-z, model 1 for average sources, Boughn and Crittenden (2002) were able to reproduce the NVSS auto-correlation function rather well. However the dotted curve in Fig. 8 shows the galaxy power spectrum , calculated using a mean bias of (in agreement with the values in Boughn and Crittenden (2002); Blake et al. (2004)) and the same model for . For our fiducial value of , the model power spectrum is deficient relative to the observed power spectrum.

Therefore, we search for a NVSS redshift distribution that better reproduces our angular power spectrum measurement. We find that for , a near Gaussian which is lopsided toward low redshift and centered at :


results in a good fit. This match to the NVSS angular power spectrum is shown in the dashed curve of Fig. 8. We have used this fiducial in all simulations in this paper.

We make no claim that our fiducial is a more accurate model for the real NVSS redshift distribution than the previously considered model. It is just a device for generating simulations with the same power spectrum as the data, so that we do not underestimate our error bars. As a check, in Section IV we compared Monte Carlo based error estimates for WMAP data versus NVSS data on one hand, and WMAP data versus NVSS simulations on the other, and obtained agreement (Fig. 5). Using the dotted line in Fig. (8) would underestimate the power spectrum errors by due to the disagreement with the power spectrum seen in the data. We have not investigated the reason for the disagreement in detail since it is somewhat peripheral to the primary purpose of this paper. However, the redshift distribution and galaxy bias assumed in the modeling would be critical if we were to infer constraints on cosmological parameters (such as the normalization of matter fluctuations or the total matter density ) from our measurements of the NVSS angular power spectrum and the cross correlation . We return to this issue in §IX.

Vi WMAP systematics

Because our lensing estimator receives contributions from CMB anisotropies on small angular scales (Fig. 3), the WMAP systematics most likely to affect the detection are point sources and beam effects. In our pipeline, beam effects are incorporated by convolving the CMB with an isotropic beam (Eq. (12)) which is different for each DA. This is approximate in two ways: first, the real WMAP beams are not perfectly isotropic, but contain asymmetries which also convolve small-scale modes of the CMB by a sky varying kernel defined by the the details of the scanning strategy. Second, the isotropic part of each beam is not known perfectly; uncertainty in the beam transfer function acts as a source of systematic error in our lensing detection. We study these two effects in §VI.1VI.2.

In §VI.3, we consider Galactic microwave foregrounds and show that their effect on the lensing detection is small. Point sources and thermal SZ will be treated separately in §VII, §VIII. The ISW effect Sachs and Wolfe (1967) does not affect our lensing estimator, since the signal is negligible on CMB angular scales () which contribute. The Rees-Sciama effect Rees and Sciama (1968) would give a small contribution on these scales, but we will ignore it since it is negligible compared to the SZ signal.

vi.1 Beam asymmetry

The WMAP beams are asymmetric due to: 1) the feeds not being at the primary focus, and 2) substructure caused by 0.02 cm rms deformations in the primary mirror Page et al. (2003). The Q-band beams are elliptical with minor/major axis ratio of . The V and W-band beams show significant substructure at the to dB level, leading to distortions in the inferred power spectrum Hinshaw et al. (2006).

Although deviations from azimuthal symmetry of the beams have a small effect when estimating the WMAP temperature power spectrum, it is unclear whether the same is true when estimating lensing. At an intuitive level, CMB lens reconstruction recovers degree-scale modes of the lensing potential indirectly, through their distorting effect on smaller-scale hot and cold spots in the CMB. Beam asymmetries which convolve the small-scale CMB modes have a qualitatively similar effect and may be degenerate with lensing. For example, a beam quadrupole imparts an overall ellipticity or shear to the hot and cold spots.

Figure 11: Result of convolving a single noiseless CMB realization with the WMAP V1 beam, including beam asymmetry. We have shown the output map separated into contributions from different beam multipoles: (isotropic component, top left), (top right), (bottom left), and (bottom right). Each map has been scaled independently for visibility; the RMS temperature in the maps is 88, 0.4, 1.0, 0.04 K. The convolution with the multipoles is scan dependent and shows alignments with the ecliptic poles reflecting the WMAP scan strategy.

To incorporate beam asymmetry into our pipeline, we expand the beam profile in spherical harmonics . The multipoles of the beam represent the azimuthally averaged beam and are already incorporated in both the analysis and simulation directions of our pipeline. The higher- multipoles have been estimated by the WMAP team and represent corrections to the azimuthally symmetric approximation. In Appendix D, we show how to incorporate the higher multipoles into the simulation direction of the pipeline, generalizing the convolution in Eq. (12). In contrast to the multipoles, convolving with the higher multipoles depends on the scan strategy; our method incorporates the details of the WMAP scan based on full timestream pointing. In Fig. 11, we illustrate our simulation procedure for a single noiseless realization in -band, showing the contribution of the multipoles to the beam-convolved map.

It would be very difficult to incorporate asymmetric beams into the analysis direction of the pipeline, so our approach is to treat beam asymmetry as a source of systematic error. We assign each lensing bandpower a systematic error given by the Monte Carlo RMS change in the bandpower when the same WMAP + NVSS simulation is analyzed with and without including beam asymmetry in the simulation pipeline. We find that the systematic error in each band is small compared to the statistical error. The result is shown, as part of a larger systematic error budget, in the “Beam asymmetry” column of Tab. 1 in §IX.

vi.2 Beam uncertainty

We have shown that systematic errors from beam asymmetry are small, so that the beam may be treated as the simple convolution in Eq. (12) to a good approximation. This leaves only one remaining beam-related source of systematic error: measurement uncertainty in the beam transfer function .

We model the beam transfer function uncertainty following (Hinshaw et al., 2006, §A.2). The beam covariance matrix is dominated by a small number of modes. We SVD decompose the matrix for each DA and keep only the 10 most significant modes. Then we construct realizations of the beam transfer function using


where is the standard beam transfer function, are unit-variance normal random deviates, and are the beam covariance modes.

Armed with this simulation procedure, we assign systematic errors by computing the RMS change in each bandpower when the same simulation is analyzed with and without simulated beam uncertainty. We find that the systematic errors are extremely small.

vi.3 Galactic foregrounds

Figure 12: Foreground templates used in this paper, shown with Kp0 mask (§II) applied. Left panel: Dust template, based on Finkbeiner et al. (1999) with frequency dependence given by Eq. (23). Right panel: Free-free template, based on Finkbeiner (2003); Bennett et al. (2003b) with frequency dependence given by Eq. (24). The masked RMS of the templates in V-band is 6.4 K and 4.8 K respectively.

In addition to the CMB, the sky at microwave frequencies contains other foreground signals which must be considered as a source of systematic error in lensing. We will find that the most important of these are point sources and the thermal Sunyaev-Zeldovich effect, which will be discussed in §VII and §VIII respectively. The other relevant microwave foregrounds are Galactic in origin: dust, free-free emission, and synchrotron radiation. For descriptions of the foreground components, we refer the reader to Bennett et al. (2003b).

Following Hinshaw et al. (2006), we will model dust contamination by adding a template derived from “Model 8” from Finkbeiner et al Finkbeiner et al. (1999), evaluated at 94 GHz and scaling to frequency by:


where denotes antenna temperature. The dust template is shown in Fig. 12, left panel.

When we cross-correlate simulations of WMAP and NVSS, we find that including the dust template in the WMAP simulation results in a very small change in the estimated lensing signal. We take the Monte Carlo RMS average of the change in each bandpower when the same pair of simulations is analyzed with and without the template as a systematic error estimate, shown in the “Dust” column of Tab. 1 in §IX.

One might worry that this way of assigning systematic errors, based entirely on simulations, is too optimistic because it fails to account for unknown correlations between the templates and the real datasets. As a check, we obtain consistent results if we cross-correlate an ensemble of WMAP simulations against the real NVSS data, or the real WMAP data (with and without template subtraction) against an ensemble of NVSS simulations. Finally, when the real WMAP and NVSS datasets are cross-correlated with and without template subtraction, the change in each bandpower is consistent with our systematic error estimates, and no evidence for an overall bias is seen.

We treat free-free emission similarly; in this case we use the full-sky H map from Finkbeiner (2003), with the correction for dust extinction from Bennett et al. (2003b), and frequency dependence:


where K/Rayleigh and denotes the H intensity. Again we find consistent systematic errors in the simulation-simulation, simulation-data, and data-data cases described in the previous paragraph. The results are shown in the “Free-free” column of Tab. 1 in §IX; the systematic errors from free-free are slightly higher than dust, but still small.

Finally, we turn to Galactic synchrotron emission. The WMAP team has derived synchrotron templates both from the Haslam 408 MHz survey Bennett et al. (2003b); Haslam et al. (1981), and internally by differencing the K and Ka band WMAP channels Hinshaw et al. (2006). However, both of these templates are intended for use at degree scales, and do not have sufficient resolution to measure the sychrotron signal on the angular scales () which contribute to our lensing estimator. Therefore, it would not be meaningful to assign systematic errors from synchrotron emission by using either of these templates.

In the absence of a template for synchrotron, the best we can do is to make the assumption that the synchrotron contamination at is comparable to the other Galactic foregrounds. In V-band, synchrotron, free-free, and dust emission all contaminate the CMB at roughly similar levels Bennett et al. (2003b). In addition, synchrotron and dust appear to have similar spatial distributions (Hinshaw et al., 2006, Fig. 5), so the dust template should give us a reasonable estimate of possible synchroton contamination. However, a direct test of this assumption will have to await future higher-resolution measurements of synchrotron emission.

These results and the consistency of our measurement between frequencies (Fig. 6) lead us to conclude that our lensing detection is not contaminated by significant residual foregrounds. However, we quantify it by assigning each lensing bandpower a total systematic error from foregrounds by adding the systematic errors from the dust and free-free templates (treating the two as correlated) and then doubling each RMS error to account for a synchrotron contribution with the same order of magnitude. The result is shown in the “Total Galactic” column of Tab. 1 in §IX.

Vii Point source contamination

Point sources which are bright enough to be resolved by WMAP are excluded by the Kp0 mask (§II), but unresolved point sources act as a contaminating signal in the CMB. If the unresolved CMB point source signal were uncorrelated to NVSS, we would not expect point sources to affect our lensing estimator significantly. However, NVSS radio galaxies will contribute some nonzero flux at microwave frequencies and so appear directly as part of the point source contribution to the CMB. In addition, CMB point sources which do not actually appear as objects in NVSS may be correlated to NVSS objects in some way, e.g. if both are tracers of the same large-scale potential. Therefore, point sources are a possible contaminant of our lensing detection.

In this section, we will place limits on the level of point source contamination and assign systematic errors. Point sources will turn out to be our dominant source of systematic error, and so we will devote considerable effort to constructing reliable error estimates.

vii.1 Point source estimator

It is difficult if not impossible to construct a realistic model which would allow the level of point source contamination to be reliably estimated from general principles. At radio and microwave frequencies, several populations of point sources have been identified Toffolatti et al. (1998); Lin and Mohr (2006); Coble et al. (2006); Giommi et al. (2007) with significant uncertainties in spectral index and clustering properties.

Therefore, our approach will be to estimate the level of point source contamination directly from the data. In this subsection, we will motivate and construct an estimator which is optimized for detecting point sources instead of CMB lensing, to use as a monitor for point source contamination. The first candidate for the point source estimator is simply the cross power spectrum .

However, consider the following toy model for point sources: suppose that there are distinct populations of unclustered Poisson point sources which appear as objects in the NVSS catalog, and the -th population has number density and constant flux per source at CMB frequencies. In this model, the cross power spectrum is


whereas the bias to the lensing estimator is proportional to


Because the right-hand sides of Eqs. (25), (26) are not related in any model-independent way, one cannot translate a value of the cross spectrum to an estimate of the point source contamination in the lensing estimator, without making implicit assumptions about the point source model.

For this reason, we next consider a different candidate for the point source estimator: the three-point estimator optimized to detect the “Poisson” bispectrum


where, following Eq. (8), denote CMB multipoles and denotes a galaxy multipole. (We will construct the estimator shortly; for now we “define” the point source estimator by writing down the bispectrum which we want to detect.)

To motivate this form, we note that the bispectrum in our toy model is


Comparing to Eq. (26), we see that each point source population makes contributions to the Poisson bispectrum and lensing estimator which are proportional. Therefore, an estimate of the Poisson bispectrum will directly translate to a systematic error estimate for the lensing estimator.

This aspect of our toy model illustrates a general point: a statistical contaminant, such as unresolved point sources, affects the lensing detection by making a contribution to the bispectrum which may be coupled to the lensing bispectrum (Eq. (9)) which is measured by our estimator. Therefore, when trying to understand point source contamination, one should first ask: what bispectrum do point sources contribute?

We will actually consider a more general point source bispectrum than the Poisson form in Eq. (27), which relaxes two assumptions of the toy model. First, we have assumed that point sources do not cluster (i.e., are purely Poisson distributed). Furthermore, we have assumed that each CMB point source appears as an object in NVSS; there is a second case to consider in which the point sources do not actually appear as objects, but are merely clustered in a way which is correlated to NVSS.

Consider a population of clustered point sources which are tracers of a Gaussian field . (We assume that the bias is absorbed into the definition of , so that the probability of a point source at position is .) For our second case, where the point sources do not appear as NVSS objects, a short calculation shows that the point source bispectrum is:


In the first case, where the sources do appear as NVSS objects, the bispectrum is given by:


where is the average temperature at CMB frequencies, is the number density of the point source population, and is the number density of NVSS.

In Eq. (30), the first term represents contributions from Poisson statistics, the second represents point source clustering on the galaxy angular scales () which contribute to the lensing detection, and the third represents clustering on CMB angular scales (). We will assume that the third term is small compared with the first two and can be neglected. This is a critical assumption for our methodology and so we justify it carefully, giving two arguments.

The first argument is that a realistic point source clustering power spectrum will be rapidly decreasing with and so the factors in the third term (with ) will be small compared with the factor in the second term (with ).

The second argument is more formal and shows that the third term in Eq. (30) is small compared to the first term. The ratio of the third and first terms is given by


In the second line, we have used the fact that the contribution to the NVSS galaxy power spectrum from the point source population alone is given by . In the third line, we have used our measurement of (Fig. 8), which shows that for . The intuition behind this formal argument is that if point source clustering were important on small angular scales, we would see this signal in the NVSS power spectrum.

We have now shown that the most general point source bispectrum is a combination of Eq. (29), and Eq. (30) with the third term neglected. This motivates our final choice of point source estimator: we will use the three-point estimator optimized to detect any bispectrum of the form


where is arbitrary (our estimator will estimate in bands). This generalizes the Poisson bispectrum considered previously (Eq. (27)).

We have shown that Eq. (32) is a sufficiently general form of the point source bispectrum to allow an arbitrary clustering power spectrum between point sources, an arbitrary cross-correlation to the NVSS overdensity field , and applies whether the CMB point sources actually appear as objects in NVSS, or are merely correlated to NVSS. Indeed, by putting an arbitrary dependence in Eq. (32), we have been conservative by allowing a very general point source contribution. However, there is one caveat: we have assumed that point sources are biased tracers of Gaussian fields. Non-Gaussian contributions from nonlinear evolution have not been included. In halo model language Cooray and Sheth (2002), we have incorporated one-halo and two-halo terms in the bispectrum but not the three-halo term.

Now that we have determined the most general bispectrum contributed by point source contamination (Eq. (32)), how do we construct the point source estimator? In Appendix B, we show that the optimal estimator for this bispectrum is constructed in a way which is analagous to the lensing estimator (or the curl null test ). First, we define a field which is quadratic in the CMB:


where was defined previously in Eq. (15). Then we cross-correlate to galaxy counts, subtracting the one-point term as usual:


This defines the optimal estimator for the point source bispectrum (Eq. (32)), with the galaxy multipole binned into a bandpower .

Intuitively, the field can be thought of as a “quadratic reconstruction” of CMB point source power, in the same sense that is a quadratic reconstruction of the CMB lensing potential. Our estimator is obtained by cross-correlating to the filtered galaxy field : we are only interested in point source power which is correlated to NVSS. By using to directly estimate the bispectrum due to point sources from data, we can assign systematic errors to the lensing bandpower which do not depend on the details of the point source model, as we will now see.

vii.2 Results

Figure 13: Point source estimator applied to the WMAP and NVSS datasets, showing no evidence for CMB point source power which is correlated to NVSS. The error bars were obtained from Monte Carlo WMAP+NVSS simulations without point sources.

In Fig. 13, we show the result of applying the point source estimator , constructed in the previous section, to the WMAP and NVSS datasets. The to zero is 11.7 with 12 degrees of freedom. Therefore, no evidence for point source contamination is seen. This lets us put strong constraints on the systematic error in lensing due to point sources: the point source contribution must be small enough to be hidden in Fig. 13, even though the estimator is optimized for point sources. The rest of this subsection is devoted to assigning systematic errors based on this observation.

We find that for distinct bands , the point source and lensing estimators in band are uncorrelated to the estimators in band . This is unsurprising; it follows from the definitions that the bands are independent for all-sky coverage and homogeneous noise, so that the only correlation is due to inhomogeneities. Since we have large sky coverage and wide bands, the correlations should be small. We will treat each band independently, for consistency with our point source model, which allows an arbitrary dependence in the point source amplitude (Eq. (32)). We will illustrate our method in detail for the band .

Figure 14: Histogrammed 1.4 GHz flux distribution in NVSS, with the fitting function in Eq. (35) shown for comparison.

First, we use simulations to study the effect of point sources on the estimators , using the following fiducial point source model. (We will show shortly that the final result does not depend on the details of the point source model.) Each simulated NVSS galaxy is assigned a randomly generated flux between 2 mJy and 1 Jy, drawn from the distribution


This distribution was obtained empirically from the flux distribution seen in the real NVSS data (Fig. 14). We then assign the flux


at each WMAP frequency , where is a constant which will be varied to simulate different overall levels of point source contamination. Following Bennett et al. (2003b), we take spectral index in our fiducial point source model.

Figure 15: Ensemble of simulations in the fiducial point source model (Eqs. (35), (36)) with varying point source amplitude . For each realization, we show the observed point source level in the band and the change in the lensing estimator due to the point source contribution. The dotted vertical line shows the point source level in this band estimated from the real WMAP + NVSS data; the smaller vertical error bar shows the mean and RMS among simulations whose observed point source level matches the measured value.

In Fig. 15, we show the values of the point source estimator obtained in an ensemble of simulations with varying point source amplitude , and the change in the lensing bandpower which is due to the point source contribution. (Note that we do not show the true point source amplitude for each simulation; we show the observed point source level , estimated the same way as in the data.)

We find that the results can be fit by treating as a Gaussian variable with mean and variance which depend on :


where K, , K.

Based on this picture, how can we assign systematic errors due to point sources? Consider the distribution of values obtained by considering only realizations whose observed point source level agrees with the value ( K) observed in the data (indicated by the dotted vertical line in Fig. 15.) Note that this distribution includes realizations with a range of values for the true point source amplitude ; we are effectively averaging over point source levels allowed by the observed value of (i.e. the posterior distribution). By Eq. (37), we get a Gaussian distribution with parameters:


indicated by the vertical error bar in Fig. 15.

We have now arrived at an distribution (Eq. (38)) for the change in which is due to the point source contribution. The central value of this distribution is nonzero; point source contamination makes a negative contribution on average, as can be seen in Fig. 15. To be conservative, we will not shift our estimate for in the positive direction by the central value (this would allow point sources to “help” the lensing detection), but will include the shift as part of the systematic error. Thus we would quote the systematic error in as: .

As we have described it, this procedure appears to depend on the fiducial point source model (Eqs. (35), (36)). However, we find that the final systematic error estimate in each band is relatively robust even under drastic changes to the model. We tried the following extreme cases: assigning constant flux to each source rather than using Eq. (35), taking spectral index in Eq. (36) rather than , and finally simulating point sources which are merely correlated to NVSS rather than appearing as NVSS objects. All of these models give similar results to within a factor . (Note that our point source estimator in Eq. (34) is actually optimized for point sources with a blackbody spectral distribution, but these results show that we obtain robust systematic error constraints across a reasonable range of spectral indices.)

Repeating this procedure for every band, we obtain a systematic error estimate for each lensing bandpower . Since we have considered several point source models, we assign the systematic error for each band using the model which gives the largest error in that band. The results are shown in the “Resolved point source” column in Tab. 1 in §IX. We find a systematic error which is smaller than the statistical error in all bands, but is the largest overall source of systematic error.

The relative robustness of our error estimate to the point source model is consistent with our discussion in the previous subsection: regardless of the details of the model, the contamination to the lensing estimator is proportional to the level of the point source bispectrum (Eq. (32)) contributed by point sources. By directly estimating the bispectrum, we can obtain a relatively model-independent constraint on the systematic error due to point sources. This would not be possible if a simpler statistic were used, such as the cross power spectrum .

The procedure we have described is similar to the Fisher matrix based method that is frequently used to marginalize point sources when estimating primordial non-Gaussianity from the CMB bispectrum Komatsu et al. (2005), but differs in several details. First, we use a general form of the point source bispectrum (Eq. (32)) which allows point source clustering, and also allows CMB point sources to appear or not appear as NVSS objects. Second, we do not shift the lensing estimator by the central value of the posterior distribution in Eq. (38), but treat the shift as part of the systematic error. Finally, the Fisher matrix formalism would not predict the increased variance in in the presence of point sources (Eq. (37)). This is included in the Monte Carlo based procedure presented here. The Fisher matrix does predict the overall negative slope in Fig. 15, which is a property of the point source and CMB lensing bispectra. As a check, if we directly compute the Fisher matrix (see Eq. (43) below), we find a weak negative correlation ( in each band) between the lensing and point source shapes.

vii.3 Resolved point sources

Figure 16: Illustration of our procedure (Eq. (39)) for simulating correlations between the NVSS galaxy field (left) and source mask (right). For visibility, we have bandlimited the galaxy field to , and used 100 sources with masking radius rather than the mask parameters of the datasets (§II).

Now that we have analyzed systematic errors in lensing from unresolved CMB sources, we consider resolved sources. Resolved CMB point sources have been treated in the pipeline by simply masking each source (§II). If the sources are correlated to radio galaxies, so that the WMAP mask is correlated to NVSS, one may wonder whether the masking procedure can bias the lensing detection.

We can prove the following general result (Appendix C): in the absence of CMB lensing, correlations between the mask and galaxy field cannot fake the lensing signal, i.e. the expectation value is zero even if the mask is correlated. Interestingly, our proof depends on the presence of the one-point term in the estimator (Eq. (17)) and does not rule out the possibility of bias if this term is omitted.

Given this general result, the lowest-order effect that might be expected from mask-galaxy correlations is a bias proportional to the lensing signal, i.e. a calibration error. We looked for a calibration error in simulations, by randomly generating a point source mask by assigning a point source to pixel with probability


(This is an extreme case, corresponding to a linear bias model in the maximally biased limit .) An example of this simulation procedure is shown in Fig. 16.

With the source mask density of the real datasets (§II), we see no evidence for a calibration error after 1024 Monte Carlo simulations of the full pipeline. The same result was obtained replacing the NVSS overdensity by the lensing potential on the right-hand side of Eq. (39), or bandlimiting the right-hand side for several choices of band.

Since we do not have a general proof that the calibration error is small, we can only conclude that it is smaller than the statistical limit from our Monte Carlo sample. In Tab. 1, we have assigned each bandpower a 3% systematic calibration error in the “Resolved point sources” column, but we see no evidence for the effect and it may be much smaller.

Viii Sunyaev-Zeldovich fluctuations

A further source of possible contamination of the WMAP-NVSS correlation comes from re-scattering of the primordial microwave background off hot electrons inside the large scale structure field that also underlies the distribution of NVSS sources. The largest effect is the thermal Sunyaev-Zel’dovich (SZ) effect Sunyaev and Zeldovich (1970, 1980), due to inverse Compton scattering which shifts photons away from their originally black-body spectrum. The kinetic Sunyaev-Zel’dovich (kSZ) effect, due to Doppler scattering of CMB photons by large scale structure moving along the line of sight, is expected to be a concern for lensing reconstruction with future CMB experiments that are able to frequency clean the thermal effect Amblard et al. (2004); Huffenberger and Seljak (2005). On the angular scales relevant for WMAP, the kinetic effect is much smaller and more Gaussian than the thermal effect, and we neglect it in the following analysis.

The induced temperature change of the thermal SZ compared to the CMB, , is proportional to the line of sight integral over the cluster gas pressure, (the Compton-y parameter), where is the free electron density, the Boltzmann constant, the electron temperature, the electron mass, and the Thomson scattering cross section. It also has a characteristic frequency dependence, given in terms of the dimensionless frequency by


This frequency dependence causes changes in the expected amplitude of the SZ between the WMAP V (Q) and W channels. These differences are smaller than the statistical error of our WMAP-NVSS cross correlation measurement, making it impossible to distinguish the SZ effect from lensing on frequency basis alone.

We therefore rely on angular separation. Our preferred way to describe the SZ effect and assign systematic errors would be to use full hydrodynamical simulations of the effect (e.g. Springel et al. (2000); White et al. (2002); Seljak et al. (2001)). Unfortunately these have to date only been performed on scales of comoving Megaparsec, allowing modeling of secondary anisotropies on scales of only a few square degrees. Our lensing estimator on the other hand receives contributions from , requiring simulation on scales substantially larger than square degrees. A somewhat less computationally expensive route would be to establish halo catalogues based on perturbation theory schemes (e.g. Scoccimarro and Sheth (2002); Monaco et al. (2002)) that are then decorated with semi-analytic gas pressure profiles. Even these procedures are however very costy for our purposes of covering 40,000 square degrees on the sky at a depth of about 4 comoving Gigaparsec, under the necessary requirement of resolving halos down to solar masses in order to reliably model SZ fluctuations below Komatsu and Seljak (2002).

As we will argue in this section however, on the scales relevant to a lensing detection using WMAP, SZ contamination can be treated as part of the point source contribution which has been studied in the previous section.

To begin with, notice that although at WMAP frequencies SZ clusters contribute a temperature decrement to the CMB, their contribution to the point source estimator is positive, because the estimator is quadratic in the CMB. Therefore our “point source” estimator will be able to serve as a monitor for the sum of point source and SZ contamination. This is yet another advantage of the three-point estimator over the cross spectrum discussed in §VII.1: because point sources make a positive contribution to the cross spectrum but the SZ contribution is negative, the cross spectrum cannot constrain both contaminants.

Next consider the spatial distribution of SZ. The vast majority of the thermal SZ signal stems from collapsed regions with a gas density contrast of hundreds of times the mean density of the universe (see e.g. White et al. (2002)). If cluster profiles could be approximated as -functions, then they could be treated as biased tracers of large scale structure that is correlated to NVSS galaxies. Since our point source model (Eq. (32)) allows clustering and cross-correlation to NVSS, this would allow us to treat the SZ contribution as part of the point source contribution.

To quantify the deviation from pointlike profiles, in Fig. 17, we show galaxy cluster profiles in angular multipole space, calculated with the universal gas-pressure profile model of Komatsu and Seljak (2001), at z=0.1 and z=1.0. This redshift range is chosen to span roughly the range where the SZ might be correlated to NVSS sources. It can be seen that many of the relevant clusters fall below the angular scale () where our lensing reconstruction gathers most of its information, but some large nearby SZ clusters have profiles as extended as tens of arcminutes, and show some slope at the relevant angular scales.

Figure 17: The Compton-y profile for three different cluster masses at z=1 (thick lines) and z=0.1 (thin lines). The profiles have been normalized to 1 at l=0 to facilitate comparison. According to this panel, at high redshift it may be possible to approximate even rare and massive clusters as point sources on the scale where our lensing estimator gathers most of its information, .

To determine whether this slope is important at WMAP resolution, we consider the angular power spectrum , which is an average over redshift and mass of all clusters. In cross correlation with NVSS, this integral would be modulated by the source redshift number density. Since the NVSS redshift distribution is not very well understood, here we apply uniform weight to all objects to obtain an estimate for the scale dependence of the power spectrum. We calculate the power spectrum including both the Poisson (1-halo) and clustering (2-halo) contributions, following the formalism of Komatsu and Kitayama (1999); Komatsu and Seljak (2002). The results are shown (for the low frequency (Rayleigh-Jeans) limit in which ) in Figure 18.

Figure 18: The thermal Sunyaev-Zel’dovich angular power spectrum contributions (in the Rayleigh-Jeans limit) from Poisson and clustering terms. On the scale of most interest for our lensing reconstruction, , the SZ Poisson term dominates by an order of magnitude over the clustering part. The angular power spectra were calculated using the gas pressure profile model by Komatsu and Seljak (2001, 2002).

It is seen that the SZ power spectrum is not flat at , owing to the contribution of the most massive and nearby clusters, but has the rough scaling


over the relevant range of angular scales.

We can incorporate this scale dependence into the analysis by considering a bispectrum of the form


To quantify the effect of scale dependence on the lensing estimator, we compute the correlation between this shape and the point source shape (Eq. (32)), using the Fisher matrix formalism Heavens (1998). According to this, the Fisher matrix element between two bispectra is defined by


To a good approximation, when bispectra are estimated from data, the covariance matrix is given by:


When we compute the Fisher matrix for the point source (Eq. (32)) and scale-dependent (Eq. (42)) shapes at WMAP and NVSS noise levels, we find a correlation coefficient . At this level of correlation, the point source shape and SZ shape can not be distinguished to , unless a 6 detection of the point source shape can also be made. Since we do not find any evidence for point source contamination in the data (Fig. 15), we conclude that the difference between the point source and SZ bispectra should be negligible in the context of the WMAP and NVSS data sets.

As an additional check, we tried modifying our point source simulations by giving each point source an profile, and SZ frequency dependence (Eq. (40)), including the negative sign. This crude procedure is of course not an accurate method for simulating SZ in detail, but does incorporate two qualitative features which distinguish SZ from point sources at WMAP resolution: the scale dependence (Eq. 41) and frequency dependence (Eq. 40). We find that the systematic errors in lensing (obtained from Monte Carlo simulations as described in §VII) are within the range of point source models previously considered, showing that neither of these deviations from pure point source behavior significantly affects our method.

Finally, there is one assumption in our point source model which we can check explicitly for the case of SZ: that clustering is unimportant on scales of (see Eq. 30). This can be seen directly from Fig. 18; the clustering term is dominated by the Poisson term by an order of magnitude.

Ix Final result and discussion

Beam Galactic Point source + SZ
Statistical Asymmetry Uncertainty Total Dust Free-free Total Unresolved Resolved Total Stat + systematic
Table 1: Final estimated bandpowers, together with statistical uncertainties and systematic errors from point sources. All entries in the table are in multiples of .

Figure 19: Final result from Tab. 1, showing statistical errors alone (blue/inner error bars) and statistical + systematic errors (red/outer).

In Tab. 1 and Fig. 19, we show our final result: an estimated value of in bandpowers, together with statistical and systematic uncertainties. Our procedure for combining errors is as follows. We combine the errors from beam asymmetry (§VI.1) and beam uncertainty (§VI.2) into a “total beam” error assuming that the two are completely correlated. We obtain a “total Galactic” error from Galactic CMB foregrounds by combining the dust and free-free systematic errors (§VI.3) assuming correlated errors, and double the result to account for synchrotron (where no template is available on the relevant angular scales). We obtain a “total point source” error by combining the errors from unresolved and resolved sources, assuming that the two are correlated. (As we have shown in §VIII, the “point source” errors apply to the total systematic error from CMB point sources and the thermal SZ effect.) We then obtain our final result by combining the statistical, total beam, total Galactic, and point source errors, assuming that the four are uncorrelated.

What is the total statistical significance of our detection? To assess this, we combine our bandpower estimates into a single estimator , giving each bandpower a weight proportional to its fiducial expectation value (not the measured value in Tab. 1) and inversely proportional to its total (statistical + systematic) variance: