A NVSS redshift window functions

# Correlation of CMB with large-scale structure: I. ISW Tomography and Cosmological Implications

## Abstract

We cross-correlate large scale structure (LSS) observations from a number of surveys with cosmic microwave background (CMB) anisotropies from the Wilkinson Microwave Anisotropy Probe (WMAP) to investigate the Integrated Sachs-Wolfe (ISW) effect as a function of redshift, covering . Our main goal is to go beyond reporting detections towards developing a reliable likelihood analysis that allows one to determine cosmological constraints from ISW observations. With this in mind we spend a considerable amount of effort in determining the redshift-dependent bias and redshift distribution () of these samples by matching with spectroscopic observations where available, and analyzing auto-power spectra and cross-power spectra between the samples. Due to wide redshift distributions of some of the data sets we do not assume a constant bias model, in contrast to previous work on this subject. We only use the LSS data sets for which we can extract such information reliably and as a result the data sets we use are 2-Micron All Sky Survey (2MASS) samples, Sloan Digital Sky Survey (SDSS) photometric Luminous Red Galaxies, SDSS photometric quasars and NRAO VLA Sky Survey (NVSS) radio sources. We make a joint analysis of all samples constructing a full covariance matrix, which we subsequently use for cosmological parameter fitting. We report a 3.7 detection of ISW combining all the datasets. We do not find significant evidence for an ISW signal at , in agreement with theoretical expectation in CDM model. We combine the ISW likelihood function with weak lensing of CMB (hereafter Paper II Hirata et al. (2008)) and CMB power spectrum to constrain the equation of state of dark energy and the curvature of the Universe. While ISW does not significantly improve the constraints in the simplest 6-parameter flat CDM model, it improves constraints on 7-parameter models with curvature by a factor of 3.2 (relative to WMAP alone) to , and with dark energy equation of state by 15% to [posterior median with “” (16th–84th percentile) range]. A software package for calculating the ISW likelihood function can be downloaded at http://www.astro.princeton.edu/~shirley/ISW_WL.html.

###### pacs:
98.80.Es, 95.36.+x, 98.65.Dx.

## I Introduction

The Cosmic Microwave Background (CMB) has provided us with a wealth of cosmological information. The large-scale anisotropies were first discovered by the Differential Microwave Radiometer (DMR) on Cosmic Background Explorer (COBE) satellite Smoot et al. (1992), and the smaller-scale CMB anisotropies were subsequently measured by various ground-based/balloon-borne experiments. More recently, the Wilkinson Microwave Anisotropy Probe (WMAP) satellite Bennett et al. (2003); Jarosik et al. (2007) produced a cosmic variance limited map of CMB anisotropies down to . The structure of the angular power spectrum when combined with other cosmological probes (such as Sloan Digital Sky Survey, Tegmark et al. (2006), Hubble Key Project Freedman et al. (1994) and 2dF Galaxy Redshift Survey Cole et al. (2005)), allows extremely precise measurements of the cosmological parameters of the CDM model. While most of the fluctuations seen by WMAP and other CMB experiments were generated at the last surface of scattering, structures formed at low redshift also leave imprints on the CMB. These anisotropies, such as the thermal Sunyaev-Zeldovich (tSZ) Sunyaev and Zeldovich (1980a) and kinetic Sunyaev Zeldovich effects (kSZ) Sunyaev and Zeldovich (1980b), the Integrated Sachs-Wolfe (ISW) effect Sachs and Wolfe (1967), and gravitational lensing, contribute only slightly to the CMB power spectrum on scales measured by WMAP, but they can be detected by cross-correlating the CMB with suitable tracers of the large scale structure.

This is the first of two papers that measure the Integrated Sachs-Wolfe effect and gravitational lensing (Paper II) in cross-correlation. In this paper, we focus on large scale galaxy-temperature correlations and their large scale cosmological source, the Integrated Sachs-Wolfe (ISW) effect. The ISW effect results from the red- or blue-shifting of the CMB photons as they propagate through gravitational potential wells. As the potential wells of the Universe (i.e., the spatial metric) evolve, the energy gained by photons falling into the potential well does not cancel out the energy loss as photons climb out of the well. This is important at late times when the Universe is not matter dominated and the gravitational potential is time dependent. It is only significant on large scales, since on small scales the amount of time spent by the photon in each coherence region of the gravitational potential is small and any small scale fluctuations will be smoothed out as the photon go through numerous potential wells along the line of sight.

To measure the above effect, we cross-correlate the CMB temperature anisotropies with maps of galaxies from the Two Micron All Sky Survey (2MASS), luminous red galaxies (LRGs) and quasars from the Sloan Digital Sky Survey, and radio sources from the NRAO VLA Sky Survey (NVSS). This incorporates most of the LSS tracers used by previous efforts Boughn et al. (1998); Fosalba et al. (2003); Scranton et al. (2003); Afshordi et al. (2004); Boughn and Crittenden (2004a); Fosalba and Gaztañaga (2004); Nolta et al. (2004); Padmanabhan et al. (2005a); Gaztañaga et al. (2006); Cabré et al. (2006); Giannantonio et al. (2006); Vielva et al. (2006); Pietrobon et al. (2006); McEwen et al. (2007); Rassat et al. (2007) to detect the ISW effect. Our goal in this work extends this previous literature by going beyond detecting the ISW effect to measuring its redshift evolution and using that to constrain different cosmological models (e.g. the ISW effect due to spatial curvature occurs at significantly higher redshifts than that due to a cosmological constant). We therefore require a large redshift range ( to ) but with sufficient redshift resolution to unambiguously discern any redshift evolution of the signal. In addition, to draw robust cosmological conclusions from an observed redshift evolution, we must constrain both the redshift distribution and evolution of the bias with redshift for each of the samples; the simple assumption of constant bias is in most cases no longer sufficient. These considerations drive our survey selections; we discuss these in more detail in Sec. VIII. Our final product is a likelihood code that can be applied to any cosmological model. In addition to providing complementary constraints on standard cosmological parameters, we expect it can be a strong discriminator of the modified gravity models, which have very distinctive ISW predictions Song et al. (2007).

We review the theory behind the ISW effect in Sec. II. The CMB and LSS data sets used are described in Sec. III; the results of cross-correlating the two are in Sec. IV. Sec. V and  VI constrain the redshift distributions of the samples, and possible systematic contamination of the cross-correlations. Sec. VII presents the cosmological implications of these results, and Sec. VIII summarizes our conclusions. The companion paper (Paper II) uses the same data sets to detect the weak lensing of the CMB. All of the theoretical predictions are made with WMAP 3 year parameters (=, = , , , ) except in Section V or otherwise stated.

## Ii Theory

We briefly review the ISW effect and its cross-correlation with the galaxy density (see also Refs. (Peiris and Spergel, 2000; Cooray, 2002; Afshordi et al., 2004)). The temperature anisotropy due to the ISW effect is expressed as an integral of the time derivative of the gravitational potential over conformal time ,

 ΔTISW(\boldmath^θ)=2∫η0ηrdη∂ϕ∂η, (1)

where and are the conformal time at recombination and today, respectively, and we ignored the effect of Thomson scattering suppression, which is negligible for the redshift range of interest here. For scales sufficiently within the horizon, the gravitational potential is related to the mass fluctuation in Fourier space by the Poisson equation:

 ϕ(k,z)=−32H20c2Ωm(1+z)δ(k,z)k2, (2)

where is the ratio of the matter density to the critical density today, is the Hubble constant today, is the speed of light, is the redshift, and is the comoving wave number. On large scales where the mass fluctuation , the perturbations grow according to linear theory = .

We are interested in cross–correlating the temperature anisotropies, , with the observed projected galaxy overdensity . The intrinsic angular galaxy fluctuations are given by:

 g(\boldmath^θ)=∫dzb(z)Π(z)δ(χ(z)% \boldmath^θ,z), (3)

where is an assumed scale-independent bias factor relating the galaxy overdensity to the mass overdensity, i.e. , is the normalized selection function, and is the comoving distance to redshift . We focus on the cross-spectrum of the galaxies with the CMB temperature fluctuation:

 CgTℓ=2π∫k2dkP(k)[g]ℓ(k)[T]ℓ(k) (4)

where is the matter power spectrum today as a function of the wave number , and the functions and are

 [g]ℓ(k)=∫dzbi(z)Π(z)D(z)jℓ(kχ(z)) (5)

and

 [T]ℓ(k) = 3H20c2ΩmTCMB (6) ×∫dzddz[D(z)(1+z)]jℓ(kχ(z))k2.

The Limber approximation, which is quite accurate when is not too small (), can be obtained from Eq. (4) by setting and using the asymptotic formula that (when ). We find that the substitution is a better approximation to the exact expressions than . This gives

 CgTℓ = 3ΩmH20TCMBc21(ℓ+1/2)2 (7) ×∫dzb(z)Π(z)H(z)cD(z)ddz[D(z)(1+z)] ×P(ℓ+1/2χ).

The above discussion ignores the effects of gravitational lensing, which alters the expected signal through two competing effects – changing the flux limit of the survey as well as the observed galaxy density. Both of these effects can be thought of as altering the redshift distribution of the tracers, and so we defer the discussion to Sec. V.

## Iii Data

We describe the CMB and galaxy data sets used in our analysis below; these are summarized in Table 1. The data sets not used in this paper are discussed further in the Sec. VIII, where we provide detailed explanations for the choices made. All large scale structure data were pixelized in the HEALPix system with the resolution and sky coverage shown in Table 1.

### iii.1 CMB temperature from WMAP

The WMAP mission Bennett et al. (2003); Jarosik et al. (2007) measured the all-sky maps of the CMB at multipoles up to several hundred. We use the second public data release of the WMAP data with the first three years of observations. The all-sky CMB maps are constructed in the following bands: K (23 GHz), Ka (33 GHz), Q (41 GHz), V (61 GHz) and W (94 GHz). These maps are pixelized in the HEALPix Górski et al. (2005) resolution 9 format with pixels, each 47.2 sq. arcmin in area. These maps are not beam-deconvolved and this, with the scan strategy of WMAP, results in nearly uncorrelated Gaussian uncertainties on the temperature in each pixel Jarosik et al. (2007). We limit our analysis to Ka through W band as the K-band is heavily contaminated by the Galactic emission. We trim all masks with the WMAP Kp0 mask and point source mask to remove regions contaminated by Galactic emission and point sources, leaving 76.8% ( resolution 9 HEALPix pixels) of the sky for the ISW analysis. We choose not to use either the WMAP “Internal Linear Combination” (ILC) map or the foreground cleaned map to avoid a number of practical difficulties as these maps lose frequency dependence of the original maps and have complicated pixel-pixel noise correlations.

### iii.2 Two Micron All Sky Survey (2MASS)

We use galaxies from the Two Micron All Sky Survey (2MASS) Extended Source Catalog (XSC) Skrutskie et al. (1997); Jarrett et al. (2000); Skrutskie et al. (2006) as mass tracers of the low redshift Universe. The median redshift of these objects is . We use , the -band isophotal magnitude measured inside a circular isophote with surface brightness of 20 mag arcsec, as our default flux measure. We extinction correct the magnitudes from the catalog using the reddening maps Schlegel et al. (1998):

 K20=K20,raw−AK, (8)

where (Afshordi et al., 2004). Note that we ignore changes to the isophotal radius due to extinction. We remove regions with in the dataset as the galaxy density starts to drop drastically. We visually inspects how the galaxy density changes with and decide to cut with as there is a drastic drop. There are galaxies in the 2MASS XSC after removing known artifacts and sources in close proximity to a large galaxy (’a’ and ’z’) and requiring (which rejects duplicate observations of the same part of the sky). The 2MASS XSC can miss objects near bright stars or overlapping artifacts, and so we used the XSC coverage map Jarrett et al. (2000) and masked out pixels with coverage, thus of the sky.

We divided the 2MASS sample into 4 flux bins: , , , . Note that the redshift distribution of these 4 bins actually overlap significantly. Our sample selection for 2MASS is similar to Afshordi et al. (2004) except the pixelization.

### iii.3 Data from Sloan Digital Sky Survey (SDSS)

The Sloan Digital Sky Survey has taken CCD images of deg of the high-latitude sky (York et al., 2000). A dedicated 2.5m telescope (Gunn et al., 1998, 2006) at Apache Point Observatory images the sky in photometric conditions (Hogg et al., 2001) in five bands () (Fukugita et al., 1996; Smith et al., 2002) using a drift-scanning, mosaic CCD camera (Gunn et al., 1998). All the data processing are done by completely automated pipelines, including astrometry, source identification, photometry (Lupton et al., 2001; Pier et al., 2003), calibration (Tucker et al., 2006; Padmanabhan et al., 2007a), spectroscopic target selection (Eisenstein et al., 2001; Strauss et al., 2002; Richards et al., 2002), and spectroscopic fiber placement (Blanton et al., 2003). The SDSS is well underway, and has produced seven major releases (Stoughton et al., 2002; Abazajian et al., 2003, 2004, 2005; Adelman-McCarthy et al., 2006, 2007; Adelman-McCarthy and for the SDSS Collaboration, 2007).

In addition to constructing LRG and quasar maps, we constructed three additional maps that we use to reject region sheavily affected by poor seeing or stellar contamination. These include (i) a map of the full width at half-maximum (FWHM) of the point-spread function (PSF) in band; (ii) a map of stellar density ( stars, smoothed with a 2 degree FMHM Gaussian); and (iii) a similar map using only the red stars ().

All SDSS magnitudes used here are extinction-corrected using the maps of Ref. Schlegel et al. (1998). We use SDSS model magnitudes for the LRGs, and PSF magnitudes for the quasars and stars.

#### Luminous Red Galaxies

We use the photometric Luminous Red Galaxies (LRGs) from Sloan Digital Sky Survey (SDSS) constructed as described in Padmanabhan et al. (2005b). The LRGs have been very useful as a cosmological probe since they are typically the most luminous galaxies in the Universe, thus probing a larger volume than most other tracers. On top of this, they also have very regular spectral energy distributions and a prominent 4000Å break, making photo- acquisition much easier than the other galaxies. We will not be repeat our selection criteria here as it is thoroughly described in Padmanabhan et al. (2005b). We only accept sky regions with (almost identical to as in Padmanabhan et al. (2005b)) and an band FWHM arcsec.

Furthermore, there are a few regions in SDSS that have 60% more red stars than typical for their galactic latitude; we suspect photometric problems and rejected these regions. The red star cut removed deg in assorted parts of the sky.

We slice our LRG sample into two redshift bins for the ISW analysis: and .

#### Photometric quasars

We select quasars photometrically from the Sloan Digital Sky Survey by first generating a candidate quasar catalog consisting of UVX objects Richards et al. (2004). These are point sources with excess UV flux (i.e. ) observed magnitudes fainter than 14.5 (to avoid saturation problems), extinction corrected magnitudes brighter than , and u-band error less than mag ( detection in ). We call this the ALL-UVX catalog. We also have the public catalog of photometric quasars from Data Release 3 (DR3) generated by Ref. Richards et al. (2006), which we will call DR3-QSO objects. We also construct a UVX object list from only DR3 data, denoted DR3-UVX. This catalog is used to extend the selection and photometric redshifts from the DR3 region to the ALL region. Ideally the catalog would have been based on running the algorithm of Ref. Richards et al. (2006) on the ALL region but this option was not available at the time we constructed the quasar catalog.

We first match the DR3-UVX objects to the DR3-QSO objects and then assign the photometric redshifts from the DR3-QSO objects to the matched DR3-UVX object. For objects that are in DR3-UVX catalog, but not in the DR3-QSO catalog, we mark them as rejects. We now have a DR3-UVX catalog with every object either assigned a redshift or marked as a reject. The reject rate for DR3-UVX (ALL-UVX) is (). Then, we lay down the DR3-UVX catalog in color (,,,) space, and then for each ALL-UVX object, we find its nearest neighbor in this color space, then assigning it the same “redshift” as its matched DR3-UVX neighbor. If the DR3-UVX object has a redshift (not a reject), then the ALL-UVX object is classified as a quasar with the same redshift (photo- only), otherwise it is rejected. This procedure generates a photometric catalog of quasars in the full survey area, based on the matching against DR3 quasars in color space. However, this catalog only has the photometric redshifts, but not the actual redshift distribution. The actual redshift distribution will be discussed in Sec. VI. The average color offsets of the quasar candidate to its match for , , and are , , and , while the typical errors on the colors of the candidates are (), (), () and (). As the color differences between the match and the candidate are well within the error of the colors, we conclude that the quasar candidates are matched with high accuracy.

We then cut the catalog according to and FWHM arcsec. These cuts are determined when we look at the variation of the quasar number overdensity over a range of extinction and seeing. Also, since quasars are more sensitive than LRGs to extinction (as a result of the importance of the filter in selecting quasars), we cut the catalog at a lower . We also imposed a cut rejecting regions with more than twice average stellar density, i.e. we require stars/deg.

We further divide the sample into two redshift (photo-z) bins: (low-) and (high-). This division of sample is due to the fact that there are strong emission lines (e.g. Mgii) that redshift from one filter into the next around the redshifts of , and , causing these two redshift bins to be relatively free of cross-contamination. However, as we will see, they do contain significant contamination from redshifts below and above . We therefore constrain their redshift distribution by cross-correlating these with auxiliary data sets; we discuss this further in Sec. V.

The construction of the full sample using the DR3 catalog as described above introduces one potentially worrying systematic, namely the possibility that regions of the sky observed after DR3 would have a different density of sources than DR3 regions as a result of the nearest-neighbor method misbehaving in low-density regions of color space. This would provide a spurious feature in the quasar maps that resembles the DR3 coverage map. In order to check for this problem, we look for correlations between observing dates (if the ALL sample is misbehaving, it will be different from DR3 sample) with galaxy overdensity, and we do not find any significant correlations (Fig. 2). We also look at the correlation between quasar overdensity and the stellar number density to see if there is significant stellar contamination, we do not find any either (Fig. 2).

### iii.4 NRAO VLA Sky Survey (NVSS)

The NRAO VLA Sky Survey (NVSS) is a 1.4 GHz continuum survey covering the entire sky north of declination using the compact D and DnC configurations of the Very Large Array (VLA) Condon et al. (1998). The images all have 45 arcsec FWHM resolution and nearly uniform sensitivity and yield a catalog of almost discrete sources stronger than mJy.

This survey has several potentially major artifacts: Galactic synchrotron emission, spurious power from bright sources and a declination-dependent striping problem. All of these have to be treated properly before one can claim that the power coming from the cross-/auto-correlation is not due to some spurious issues. The Galactic synchrotron emission can in principle be an issue because it contributes significantly to the noise temperature of the VLA, and for realistic number counts, increased noise temperature could change the number of sources with measured flux above some threshold. (As an interferometer the VLA is not directly sensitive to the diffuse synchrotron foreground.) This issue is treated by incorporating a template – the Haslam map Haslam et al. (1981) – in the cross-correlation analysis and projecting out the power that are correlated to this template. Even though the Haslam map is at 408 MHz, the frequency dependence of the galactic synchrotron emission is fairly flat, allowing us to use it as a template of the Galactic synchroton radiation. The bright sources are problematic since the VLA has a finite dynamic range ( in snapshot mode with limited -plane coverage) and thus the identification of faint sources in fields with a bright source is unreliable. This issue is mitigated by masking out all the bright sources. Striping is a known systematic effect in NVSS (Blake and Wall, 2002): the galaxy density has a systematic dependence on declination, which can mimic long-wavelength modes in the galaxy field. To deal with the above potential problems, we first impose a flux limit of 2.5 mJy (where NVSS is 50% complete), mask out a 0.6 degree radius around all the bright sources ( Jy). Then to reduce striping, we also include templates to project out the synchrotron and declination-striping modes. The implementation of this projection of spurious power will be further discussed in Sec. IV.

## Iv Cross-correlation power spectrum analysis

### iv.1 Methodology

We start by organizing the temperature fluctuations and the galaxy overdensities into a single data vector,

 x=(xB,T,xg), (9)

where is a vector with the measured CMB temperature (with the monopole and dipole subtracted) in band at every HEALPix pixel; analogously, is the tracer number overdensity. The vector has a total length where and are the number of accepted pixels for the CMB and LSS maps respectively. We suppress the band subscript for simplicity, with the implicit understanding that we always refer to the cross correlation of a single WMAP band with the tracer overdensity. The covariance matrix of is,

 C=Cdiag+(0CgT†CgT0), (10)

where is given by,

 Cdiag=(CTT+NTT00Cgg+Ngg), (11)

where is the noise matrix. The submatrices , and are defined by

 Cabij=∑lmCablY∗lm(^nai)Ylm(^nbj), (12)

where is the position (on the sky) of the point of the vector . The temperature-temperature, galaxy-galaxy and galaxy-temperature angular power spectra are denoted by and respectively.

The galaxy power spectrum is first estimated using a pseudo- estimator Efstathiou (2004), and fit by the non-linear power spectrum of Smith et al. (2003), multiplied by a constant linear bias. We project out the monopole and dipole of both these power spectra by setting the power in the modes to a value () much greater than the true power spectrum.

We parametrize as a sum of bandpowers, , with amplitudes to be estimated,

 CgTl=∑ici~Pi,l. (13)

We consider “flat” bandpowers given by

 ~Pi,l={B(l)li,min≤l

where is the product of the beam transfer function Page et al. (2003), and the HEALPix pixel transfer functions at WMAP and LSS resolution. This parametrizes the power spectrum as a sum of step functions and is useful when the shape of the power spectrum is unknown.

We estimate the by forming quadratic combinations of the data Tegmark (1997); Seljak (1998),

 qi=12xtC−1diag∂C∂ciC−1diagx. (15)

These are related to the estimated by the response matrix ,

 ^ci=∑j(F−1)ijqj, (16)

where

 Fij=12tr[C−1diag∂C∂ciC−1diag∂C∂cj]. (17)

If , then the are good approximations to the maximum likelihood estimates of the . The covariance matrix of the is the inverse of the response matrix, if the fiducial power spectra and noise used to compute correctly describe the data (in this case is the Fisher matrix, hence the notation). The matrix determines the weighting and is often called a “prior” in quadratic estimation theory. Note that this usage has nothing to do with Bayesian priors – in particular, Eq. (16) is unbiased regardless of the choice of prior (though for bad choices the estimator is not minimum variance). Implementing the above algorithm is complicated by the sizes of the datasets; the implementation we use is in (Padmanabhan et al., 2003; Hirata et al., 2004; Padmanabhan et al., 2005a), and we refer to the reader to the discussion there.

[In addition to the cross-power spectra in Eq. (14), in quadratic estimator theory one usually tries to estimate the CMB and galaxy auto-power spectra as well. Because our prior is diagonal, however, these decouple, i.e. the entries in that couple the auto-powers and cross-powers are zero. For this reason we can leave the auto-powers out of the quadratic estimator.]

As mentioned earlier, the NVSS dataset has issues that require additional processing. Assume a systematic that we characterize as follows:

 xobs=xtrue+λE. (18)

If estimate , even if is the true covariance, we will still have a biased answer. However, the substitution

 C=Ctrue+ζEEt (19)

yields an unbiased estimate of when . One can add as many systematic templates (i.e. modes to project out of the map) as desired. To immunize the NVSS correlations from possible systematics, we break the NVSS map into 74 declination rings, and for each ring include a template map consisting of either (for pixels within the declination ring) or (for all other pixels). This removes the declination-dependent stripes. We also put in the 408 MHz Haslam map Haslam et al. (1981) (technically K) as a template for the Galactic synchrotron radiation. We experimented with the values of and found that the cross-spectra are converged with the choice for the declination rings and K for the synchrotron map.

### iv.2 Priors

To generate the priors for the cross-correlation power spectrum analysis, we need the approximate autopower spectrum of the galaxies. The auto-correlation is done using the same methodology as described in Sec. IV.1. The resulting autopower spectra must be smoothed, before being used as priors. This avoids statistical fluctuations in over- or under-weighting the corresponding monopoles in the cross-correlation, which could result in underestimation of signal since we would artifically down-weight multipoles that had accidentally high power in galaxies and place more weight on multipoles that had little power. We did the smoothing in two different ways. For the cases where the redshift distribution was available early enough in the analysis (2MASS or LRG), we fit the auto-power spectrum to the non-linear matter power spectrum Smith et al. (2003) to get the linear bias. In other cases (quasars, NVSS) we did not have the redshift distribution at the time the priors were created; we created the priors by using a smoothed, splined auto-power spectrum of the sample as the prior.

In the cases where we did a fit using the nonlinear matter power spectrum, the fit biases are 1.15, 1.18, 1.20, and 1.22 (2MASS, brightest to faintest); 1.92 (LRG low-); and 1.86 (LRG high-). After generating the priors, we made several modifications to the analysis, including the inclusion of redshift-dependent bias in 2MASS. Thus while the priors were not updated since they give a good fit to the observed autopower spectrum, it should be noted that these bias values are not used in the cosmological analysis (i.e. for ISW prediction purposes).

To generate priors for the CMB, we generate the priors using the theoretical s from WMAP and take into the account of the effect of pixelization and beams by convolving with the pixel and beam window functions.

### iv.3 Results of cross-correlation

Figs. 3, 4, and 5 plot the cross-correlation between WMAP and the 2MASS, SDSS and NVSS samples respectively; the four different symbols in each of these plots correspond to the four WMAP bands we use. The observed achromatic nature of the signal is consistent with it being ISW, and is an important check for frequency dependent systematics. The two quasar samples are at the highest redshifts we can probe, so if there is an ISW cross-correlation at –2, it would mean that there is significant gravitational potential change at these redshifts. This is not expected in simplest CDM cosmology, but could be present either in models where dark energy equation of state is rapidly changing with redshift or in models where curvature plays a role. The observed lack of a signal for these redshifts therefore strongly constrains such models. Note however that the NVSS cross-correlation cannot be automatically interpreted as a detection of high redshift ISW, as (see below) it covers a wide redshift range.

## V Redshift Distributions

The basic problem is to determine for each galaxy sample and each cosmological model the function that relates the matter density to the two-dimensional galaxy overdensity :

 gi(^n)=∫∞0fi(z)δ[^n,χ(z)]dz. (20)

Eq. (20) is understood to be valid on scales where the galaxies trace the matter distribution. In the absence of magnification bias, the function is simply the product of the bias and the redshift distribution: , where is the probability distribution for the galaxy redshift. In the presence of magnification bias, which is important for the SDSS quasars and possibly the NVSS radio sources, takes on the more complicated form

 fi(z)=bi(z)Πi(z)+∫∞zW(z,z′)[α(z′)−1]Πi(z′)dz′, (21)

where is the slope of the number counts of the galaxy density as a function of flux: . Here is the lensing window function:

 W(z,z′) = 32ΩmH201+zcH(z)sin2Kχ(z) (22) ×[cotKχ(z)−cotKχ(z′)],

where is the radial comoving distance, is the sine like function (equal to in a flat Universe), and is the cotangent like function (equal to in a flat Universe).

It is in fact the function that is required if one is to predict the ISW effect in a given cosmology. It is this same function that is required to predict the linear-regime angular power spectrum of the galaxies. This section describes the method by which is obtained for each of the samples. The methods are quite different due to the different types of information available for each sample. In particular there are very few spectroscopic redshifts available for NVSS. Note however that all methods include galaxy clustering data, as this is needed to determine the bias even if the redshift probability distribution is known perfectly.

All of the numbers and plots in this section only that depend on cosmology are computed using the original WMAP third-year flat 6-parameter CDM cosmology (, , , , and ), i.e. from the first release of Ref. Spergel et al. (2007). However in the Markov chain, the functon is re-computed for each cosmological model and used to predict the ISW signal.

### v.1 2mass

The 2MASS samples go down to a limiting magnitude of . At this relatively bright magnitude, almost all objects (97.9%, after correcting for the fiber collisions) have SDSS spectra, provided of course that they lie within the spectroscopic mask. In practice there are two subtleties that can occur. One is that the bias cannot be obtained to high accuracy from linear theory because even the moderate multipoles () are nonlinear, especially for the nearest 2MASS slice, and the lowest multipoles suffer from cosmic variance. The other is that the bias varies with redshift: even though the 2MASS galaxies cover a narrow range in redshift during which the Universe expands by only %, the use of apparent magnitude to define the samples means that the typical luminosity of a galaxy varies by several magnitudes across the redshift range of interest. more biased, this effect shifts the peak of the effective redshift distribution to higher redshifts than the actual distribution .

We match the 2MASS galaxies with the SDSS MAIN galaxy sample by first defining the 2MASS sample as discussed in III.2, then we select 2MASS galaxies only within mask that is more than 90% complete. We then try to match all the 2MASS galaxies with the SDSS MAIN galaxies that are within and found that almost all of the objects from 2MASS sample have SDSS spectra. We thus use the spectroscopic redshifts of the matched SDSS galaxies to identify the redshifts of the 2MASS galaxies. The redshift distribution is binned with . The redshift distribution for each of the four slices is shown in Fig. 6.

The problem of nonlinear evolution is generally very complicated, however for ISW work we only need a solution accurate to a few tens of percent. Therefore we have used the -model Cole et al. (2005), which relates the galaxy power spectrum to the linear power spectrum via

 Pgal(k)=b21+Qk21+AkPlin(k), (23)

where is the linear bias appearing in Eq. (21). Cole et al. Cole et al. (2005) found in simulations that this function fits the galaxy power spectrum in simulations for Mpc, while the required value of varies depending on the sample. Our method is to compute the theoretical angular galaxy power spectrum via the Limber integral, and fit this to the measured treating and as free parameters. This procedure can be done either assuming is constant with redshift, or (better) taking into account the redshift-dependent bias,

 Pgal(k,z)=b2⋆b2rel(z)1+Qk21+AkPlin(k,z), (24)

where is known and is a free parameter. While there is very little evolution in the 2MASS redshift range, the nearby and distant galaxies can have very different biases because they correspond to different luminosity ranges. The results for each are shown in Table 2. is based on taking the r-band luminosities of the galaxies and using from Tegmark et al. (2006). Note that the prominent peak of redshift distribution at is a supercluster known as the Sloan Great Wall. (In principle can depend on redshift as well, so one should be careful about interpreting the fit value and indeed one can see from Table 2 that fit in this way is not stable. However the changes in seen in the table when we restrict to much lower suggest that this is not a large effect on the bias.)

The -model fits for the 2MASS sample (and the LRGs) are shown in Fig. 7.

### v.2 SDSS LRGs

Next we consider the photometric LRG sample from SDSS. The sample is faint enough that spectroscopic redshifts are unavailable for most of the objects. Fortunately, precise photometric redshifts are available for LRGs since they have very uniform spectra whose main broadband feature is a break at 400nm. This break passes through the SDSS and filters in the interesting redshift range, so the and colors of an LRG correlate very strongly with its redshift Padmanabhan et al. (2005b). The error distribution of the photometric redshifts has been calibrated using spectro-s from the 2SLAQ surveyCannon et al. (2006); this procedure, and an inversion method used to determine the actual redshift distribution given the photo- distribution, are described in Padmanabhan et al. Padmanabhan et al. (2005b). These methods were applied to determine the redshift probability distribution for the LRGs used in this sample. The redshift distributions so obtained are shown in Fig. 8.

The bias is determined by the same -model fitting procedure as we used for 2MASS. The maximum values of considered are 240 for the low- slice and 400 for the high- slice, which correspond to roughly Mpc at the typical redshifts of these samples. For the fiduciual cosmology, the low- LRG slice gives a bias of and ; the high- slice gives and . In order to reduce the possible impact of the nonlinear regime on our results, we also did fits where the maximum value of was reduced by a factor of 2 or 4. The results are shown in Table 3 and the bias estimates are seen to be consistent with each other. In what follows we have used the original () fits for the LRG bias, noting that the remaining uncertainty in is small compared to the uncertainty (change in number of sigma detection is: 0.0043 (0.0388) for low-z LRG (high-z LRG)) resulting from statistical error in the ISW signal. However we note that it is not clear how well the -model works for LRGs at small scales, and we recommend more detailed analysis before taking the very small statistical error in at face value. The -model fits are shown in Fig. 7.

For the LRGs – unlike the 2MASS galaxies – each of the two photo- slices covers a narrow redshift range and the threshold luminosity varies slowly across that range, so we expect the bias to not vary significantly across the redshift range. This expectation has been confirmed in previous angular clustering studies which found % variation from to Padmanabhan et al. (2007b), and also by our own bias analysis which finds no significant difference between the two bins. Thus we conclude that for the purposes of ISW work (where we have a sigma signal for low-z LRG (high-z LRG) correlation), variation of the LRG bias within an individual photo- bin (0.2–0.4 or 0.4–0.6) can be neglected.

We calculate the possible contribution from magnification bias given the redshift distribution of the LRGs and also an assumed cosmology. We find that the possible contribution from magnification bias is times (depending on the scale) smaller than the actual signal. Therefore magnification bias is not contributing significantly to our signal.

### v.3 SDSS quasars

The function for the quasars is more uncertain than for the LRGs. This is in part due to the limited spectroscopic coverage available, but also the difficulty of constructing quasar photo-s and the lower clustering amplitude, which leads to noisier estimates of bias parameters. The basic procedure for obtaining is thus to find a region of sky with as high spectroscopic completeness as possible while still retaining a large area; use this to obtain a preliminary estimate ; and then fit for the bias parameters using clustering data, of which several are needed if is multimodal. The remainder of this section describes the details of the determination and what possible errors can be introduced by spectroscopic incompleteness, stellar contamination, redshift-dependent bias, and cosmic magnification.

In order to determine the redshift probability distribution, we began by constructing a set of five rectangles that lie within the coverage area of the SDSS, 2QZ Croom et al. (2004), 6QZ Croom et al. (2004), and 2SLAQ Richards et al. (2005) surveys. These rectangles lie along the equator (the declination range is to ) and cover the five RA ranges 137–143, 150–168, 185–193, 197–214, and 218–230. There is a significant amount of area with coverage from all surveys that is rejected as it was found to have lower completeness in 2SLAQ because there is less plate overlap. Spectra in SDSS were required to have high confidence (zConf) Stoughton et al. (2002) and those in 2QZ, 6QZ, and 2SLAQ were required to be of high quality (quality) Croom et al. (2004).

Our coverage rectangles contained a total of 1410 low-redshift and 1269 high-redshift photo- quasars; these numbers are lower than the product of the spectroscopic coverage area and the number density of photo- quasars because some parts of the latter catalogue were rejected by our stellar density cuts. Of the low-redshift photo- quasars, we found that 257 (18%) had no spectroscopic redshift determination or low quality ones, 58 (4%) were identified as stars, and the remaining 1095 (78%) are extragalactic. For the high-redshift sample these numbers are 208 (16%), 13 (1%), and 1048 (83%) respectively. From this data we construct a preliminary redshift probability distribution for each of the photo- slices using a kernel density estimator,

 Πprelim(z)=1NexNex∑k=11√2πσe−(z−zk)2/2σ2, (25)

where is the number of matches to extragalactic objects, is the redshift of the th object, and is the kernel width. The estimator is consistent in the limit that the number of objects and at fixed . In practice, must be chosen to be small compared to the width of any real features in the redshift distribution (otherwise these are artificially smoothed out), and large enough to smooth out shot noise (and redshift-clustering noise, if significant). We have used (using changes the fit bias by only 5%). This preliminary distribution is shown in the top panel of Fig. 9. The redshift distributions in the two photo- quasar slices are multimodal due to the nature of the photo- error distribution: the quasar spectra redward of Lyman- are usually characterized by a roughly power-law continuum with superposed emission lines. This means that quasar colors oscillate as emission lines redshift into and out of the SDSS filters, resulting in an (approximately) self-intersecting locus in color space and many degeneracies in the photo- solution.

If the quasar bias were constant and magnification bias negligible, then we would have simply , with the proportionality constant being the product of the bias and the probability for a photo- quasar to actually be extragalactic. This constant could then be determined by fitting the amplitude of the quasar autocorrelation function, as has been done in most past ISW studies. However, in the real Universe quasars are known to have an evolving bias, which is potentially significant across the redshift range considered, and at redshifts lensing magnification can become significant. The magnification can be calculated from the slope of the quasar counts near the magnitude limit, which gives for the low- sample and for the high- sample. In principle the cut on the -band magnitude error () could have an additional effect since magnification will reduce ; however this is not an issue for us since at the threshold, for UVX objects we will have where the typical magnitude error is even accounting for extinction (). Since for these samples is small, we compute the magnification bias using in place of the true distribution . That is, we replace Eq. (21) with

 fi(z) ≈ bi(z)Πi(z) (26) +∫∞zW(z,z′)[α(z′)−1]Πi,prelim(z′)dz′.

This leaves only the problem of constraining the product using the clustering data, i.e. the quasar power spectrum and quasar-LRG cross-power. Unfortunately the data is not capable of constraining a full model-independent distribution, so instead we write

 bi(z)Πi(z)D(z)=A(z)Πi,prelim(z), (27)

where is the growth factor, and is a piecewise constant function of . This is equivalent to assuming that the clustering amplitude (divided by spectroscopic completeness) of the quasars is constant in redshift slices, which has been found to be a better approximation than constant bias in most quasar surveys Myers et al. (2006). For comparison, the empirical “Model 3” of Ref. Croom et al. (2001) predicts to change by only 5% from to 1.45, and by 13% from to 2.00. For the more recent model, Eq. (15) of Ref. Croom et al. (2005), these numbers are 24% and 15% respectively. At higher redshifts () there is a sharp increase in Shen et al. (2007) but UVX-selected samples do not contain objects from this redshift range.

We constrain in as many redshift slices as can be constrained using the data. In particular since the quasar redshift distributions are multimodal, we would like to be able to fit a different clustering amplitude in each peak. The treatment of the two quasar samples is slightly different due to the availability of different information in their redshift ranges, so we now discuss their redshift distributions separately. In each case, the autopower spectra were fit to linear theory up to (Mpc at ) and the quasar-LRG cross-spectra were fit up to (Mpc at ).

#### Low-z sample: 0.65<zphoto<1.45

For the low- quasar sample, we can only constrain one redshift slice. An examination of Fig. 9 shows that the distribution is actually trimodal, with peaks at , , and . A fit assuming a constant yields , with dof (). Almost all of the weight for this comes from the central () peak. We also ran two-slice fits to determine whether the clustering data constrain the amplitudes of the low- and high-redshift peaks. The first such fit is of the form

 A(z)={A1z<0.52A2z≥0.52, (28)

which allows the low-redshift slice to vary ( is the local minimum of ). This fit gives and , with dof. We also tried a two-parameter fit in which the high-redshift slice is allowed to vary:

 A(z)={A′1z<1.83A′2z≥1.83 (29)

(the local minimum of between the main and high-redshift peaks is at ). This fit gives and (), with dof. The errors on are highly asymmetric in this case because the constraint comes mainly from the quasar autopower; and are then degenerate because one only knows the total power, not how much comes from each redshift slice. The shape of the power spectrum breaks this degeneracy in principle, however in practice it is far too noisy. The fact that the high-redshift slice cannot give negative power accounts for the “hard” upper limit on .

From this exercise we conclude that the clustering data cannot independently measure the bias in either the low- or high-redshift peak. The reasons are different in each case. The low-redshift peak contained only 1.7% of the spectroscopic identifications, and thus almost certainly contains only a very small fraction of our quasars. This peak lies at the same redshift as the low- SDSS LRGs, and the quasar-LRG cross-correlation is the major constraint on . Unfortunately this cross-correlation is drowned out by the enormous Poisson noise contributed by the quasars in the other two peaks, and is detected at only . On the other hand, the LRGs oversample the cosmic density field on linear scales and cover the same region of sky as the quasars. One would thus expect that since the LRG-quasar correlation is only seen at this low significance, and the ISW effect from this redshift range contributes only a small fraction of the power in the CMB, the contribution of the low-redshift peak to the quasar-ISW correlation would be statistically insignificant. We find that the predicted peak of the quasar-ISW for only the low-redshift peak quasars is lower than the entire sample (high-z QSO) by , which is significantly smaller than the error on the cross correlation. This is run using a WMAP-3yr parameters.

The high-redshift peak contains 10% of the quasars. Its amplitude must be measured in autocorrelation due to the lack of other samples at that redshift, which is a serious drawback since only 1% of quasar pairs come from the high-redshift peak. An alternative approach to constraining its amplitude would be cross-correlation against the spectroscopic quasar sample at , but we did not pursue this approach here.

#### High-z sample: 1.45<zphoto<2.00

The high- photometric quasar sample also has a trimodal distribution: there is one peak at , a second at , and a third at . In this case however, it is the highest-redshift peak that contains most of the objects, with the middle peak in second place and only a few objects in the lowest-redshift peak. This situation makes it both possible and necessary to fit separate amplitudes for the peaks; in this case we will find that two amplitudes can be constrained, one for the two low-redshift peaks and one for the main (high-redshift) peak.

As a first step, we attempt to fit all three of the peaks with separate amplitudes,

 A(z)=⎧⎪⎨⎪⎩A1z<0.33A20.33≤z<1.18A3z≥1.18. (30)

This leads to the results , , and (1), with dof. The large error bar on indicates that this parameter cannot be constrained from the data, so we instead try a two-slice fit in which we fix . This fit gives the tighter constraints and , with dof, and it is what we use for the rest of the paper.

#### Redshift Distribution Summary

The quasar autopower spectra and quasar-LRG cross-spectra, along with the model fits, are shown in Fig. 10. For the QSO0 sample, there is excess power ( above the prediction) in the lowest- bin, corresponding to a % RMS fluctuation in the number density on scales of degrees. The two most obvious sources of such power are stellar contamination and photometric calibration errors. Given that % of the photometric “quasars” are actually stars Richards et al. (2004) and that the relative photometric calibration across the sky in SDSS is estimated to be % in the band (the worst band, but one very important for quasar work) Padmanabhan et al. (2007b), either of these seems plausible. In any case, these very low multipoles were not used in fitting the redshift distribution in either auto- or cross-power.

It is essential to test the robustness of the quasar fits, in particular against the possibility of nonlinear clustering affecting the range of multipoles used in the fits. The first way we do this is by repeating our analysis using the nonlinear matter power spectrum of Smith et al. Smith et al. (2003) in place of the linear power spectrum. In the analysis with the nonlinear spectrum, the amplitude for the low- quasar slice increases by , and the amplitudes for the and parts of the high- quasar slice increase by and , respectively. If we restrict our attention to the lowest multipoles (instead of cutting at 140 or 160), these changes are , , and . In each case the change is very small compared with the error bars. Thus we do not believe that nonlinear clustering is affecting our estimates.

### v.4 Nvss

The function for NVSS is the hardest to obtain because there are no spectroscopic samples of NVSS objects that have sufficiently high completeness to obtain the redshift distribution. Past ISW analyses Boughn and Crittenden (2004a); Nolta et al. (2004) with the NVSS have been based on the radio luminosity function of Dunlop & Peacock Dunlop and Peacock (1990), which itself was fit to a combination of source counts, redshifts for some of the brightest sources, and the local luminosity function. A constant bias was then assumed. The redshift distribution so obtained is reasonable, however it has three major drawbacks: (i) the redshift probability distribution for the faint sources (which make up most of the sample) is constrained only by the functional form used for the luminosity function and not by the data; (ii) it does not give the redshift dependence of the bias, which could be very important since the redshift range is broad, and the typical luminosity of the sources varies with redshift; and (iii) the absolute bias is constrained using the NVSS autopower spectrum, which is known to contain power of instrumental origin and hence is probably a less reliable constraint than the cross-correlation against other surveys. The alternative method to measure is by cross-correlation against the other samples whose redshift distributions are known. This method is adopted here, since it does not have any of the aforementioned problems. Its main drawback is that the other samples only probe the range out to , and little data is available to constrain above that.

#### Procedure

In order to measure the effective redshift distribution of NVSS, we must first obtain the cross-correlation of NVSS with each of the eight other samples (the four 2MASS samples, and two samples each of LRGs and quasars). This is done by using the same angular cross-spectrum estimation method as was used for the ISW analysis, and the cross-spectra are shown in Fig. 11. The main subtlety that arises is that the cross-spectrum (where and are LSS samples) can actually contain Poisson noise if there are objects that are in both samples. The Poisson noise term is of the form

 Cijℓ=Cijℓ(LSS)+¯nij¯ni¯nj, (31)

where is the number of sources per steradian in catalog , and is the number of sources per steradian that appear in both catalogs. In order to measure we must match the NVSS to each of the other samples. Note that the positional errors in NVSS are typically several arc seconds, and consequently there will always be some false matches. Therefore we estimate the fraction of matches as

 ¯ni,NVSS¯ni=NmatchNi−πθ2max¯nNVSS, (32)

where is the number of matches within some radius , and is the number of sources in catalog in the NVSS mask. This was estimated for radii of 40 and 20 arcsec, and the results are shown in Table 4.

We next computed the cross-power spectra between NVSS and each of the other samples. These spectra (after subtraction of the Poisson term) are shown in Fig. 11. The redshift distribution was then fit to the cross-power spectra. In this fit the minimum multipole used is (below which there is a large amount of spurious power in the NVSS map) and the highest- bin used was determined by the formula , where Mpc is the smallest scale to be fit and is the distance corresponding to the 20th percentile of the window function for that sample as defined in Appendix A. We have fit with a -distribution,

 fNVSS(z)=αα+1zα+1⋆Γ(α)beffzαe−αz/z⋆. (33)

This function has three free parameters, , , and . Of these the normalization may be viewed as an effective bias in the sense that ; in the absence of cosmic magnification this would be the bias averaged over the redshift distribution. The peak of the distribution is at , and controls the width of the distribution. The parameter fit gives , , and .

#### High-redshift tail

The above analysis of the NVSS distribution involved cross-correlations against several samples at . (The QSO0 sample has a small number of objects at , however they have no significant impact on the fitting of the QSO0NVSS cross-spectrum.) Thus it leaves open the issue of whether there is a tail of objects at high redshift, . Since is a product of bias times redshift probability distribution, it need not be normalized – can have any value – so there is no way to tell from the cross-correlation analysis alone whether a portion of the sample is missing. If we also use the NVSS autopower spectrum then in principle one can determine whether an additional source of angular fluctuations is necessary. However the angular clustering at fixed angular scale is much stronger at low than high redshift, and the NVSS autopower spectrum is of low signal-to-noise ratio and possibly contaminated by systematics, so we have not chosen this strategy.

An alternative approach to the high- tail is to directly match against optical/ NIR catalogs. One can then use the relation or (if multiband imaging is available) photometric redshifts. There are always some radio sources without optical identifications, however this method enables one to set an upper limit to the number of NVSS sources that can be at high redshift. For our analysis, we have matched against the COSMOS field, which has a modest solid angle (2 deg), multiband imaging allowing good photometric redshifts, and deep high-resolution coverage with the VLA. Area is required due to the low density of NVSS sources (40 deg), and high-resolution radio images are required to uniquely identify an NVSS source with an optical counterpart due to the large positional uncertainty in the NVSS ( arcsec for faint sources) (Condon et al., 1998).

The COSMOS field contains 87 NVSS sources that pass our cuts. We began by matching these to the VLA-COSMOS observations, which are much deeper and have typical positional uncertainties of arcsec Schinnerer et al. (2007). Of the NVSS sources, 79 have a match within 30 arcsec (we take the nearest source in the event of multiple matches). The 79 VLA-COSMOS sources that match to NVSS are then matched to the optical catalog Capak et al. (2007); there are 64 successful matches within 1 arcsec. This represents 74% of the original NVSS catalog. It is of course possible that there are some false matches. By adding up for each NVSS source, where is the density of VLA-COSMOS sources and is the distance to the nearest VLA-COSMOS source (or 30 arcsec if the NVSS source had no match), we estimate that there are false NVSS/VLA-COSMOS matches. A similar argument suggests that false matches of VLA-COSMOS to the COSMOS optical/NIR catalog. Thus we expect that 58.5 of the matches are correct, corresponding to 67% of the initial NVSS catalog.

We show the photometric redshift distribution of the matches (according to Mobasher et al. Mobasher et al. (2007)) in Fig. 12. Our best-fit (with the distribution) has 24% of the bias-weighted source distribution at and 8% at ; if the source bias increases with redshift, as usually found for optical quasars, this number would be lower. From Fig. 12 we see that only 2 out of 64 matches fall at , i.e. the high-redshift tail of the distribution can only exist in reality if (i) most of the 26% of the sources with failed matches to COSMOS optical/NIR data are actually at , or (ii) the sources at have a large bias. Both (i) and (ii) are physically plausible but we have no direct evidence for them.

The conservative solution in this case is to consider two limiting cases for the redshift distribution of the sources at . One case, which gives the minimal lensing signal for all cosmologies, and the minimal (maximal) ISW signal for CDM (closed) cosmologies, is to set at . In the opposite limiting case, we have assumed that all failed and incorrect NVSS matches, and all sources with (i.e. a total of 35%) are at , and have four times the clustering amplitude measured for the optical quasars (QSO1 sample), e.g. (where is the growth factor) for the fiducial cosmology; the shape of at was left unchanged from the -distribution fit. In order to understand the change of ISW and CMB-lensing signals due to changes of our assumption of the high-z end of the redshift distribution of NVSS, we look at two different redshift distributions, one with nothing at (minimal model) and the other with a ”maximal” number of sources (assuming clustering strength 4 times of the optical quasars and all the failed optical IDs are at ). We find that the signals for both ISW (average: ) and CMB-Lensing change by less than 10%, therefore, one won’t expect the unidentified high-z tail of the NVSS sources be a problem in our analysis.

#### Constraints, robustness, and alternatives

While the fit parameters are formally determined by the , it is useful to graphically display the constraints in order to show what parts of the distribution are constrained by which data. This we have done in Fig. 13. For each of the eight samples, we have plotted on the vertical axis the constant value that provides the best fit to cross-correlation with that sample and its error bar. The horizontal position is determined by the following procedure. We show in Appendix A that the estimated constant is actually given by an integral over some window function,

 ⟨^fNVSS⟩=∫∞0W(z)fNVSS(z)dz, (34)

where the window function integrates to unity. The horizontal position of the data points in Fig. 13 is the median of the window function, i.e. the redshift where . The error bars extend from the 20th to the 80th percentile of the window function.

Finally we wish to compare the redshift distribution we have obtained to that used in previous ISW studies. The previous results were based on the radio luminosity function of Dunlop & Peacock Dunlop and Peacock (1990). In each case, it appears that the authors used the luminosity function and -correction based on the spectral index to infer the redshift distribution, assumed constant bias and negligible magnification, and determined the one free parameter (the bias) by fitting to the autopower spectrum. If we do this using the fiducial WMAP cosmology and our autopower spectrum we find , and the function obtained is shown as the dashed line in Fig. 13. This curve, while roughly consistent with the NVSS-quasar and NVSS-LRG correlations, badly overpredicts the NVSS-2MASS correlation. Note that the problem cannot be fixed by changing the single bias parameter: if were reduced by a factor of to fit the 2MASS data, then the LRG and quasar data would be discrepant.

There are several possible explanations for this:

• The shape of is being modified by magnification bias.

• The extrapolation of the luminosity function to faint sources at high redshift by Dunlop & Peacock is in error.

• It is possible that the Dunlop & Peacock redshift distribution accurately describes the NVSS sources, but the bias increases with redshift so as to produce the shape seen in Fig. 13.

• The cut imposed by us (and by other ISW groups) that requires NVSS sources to be unresolved is selecting against nearby objects, and hence pulling down the low- part of the curve.

Of these, possibility number 1 is easy to rule out. Application of Eq. (21) implies that has a maximum change due to magnification bias of (), and a smaller change at lower redshift ( at ), where is the source count slope. The NVSS point source counts suggest a slope of between 2.5 and 5.0 mJy, and between 5 and 10 mJy, which suggests that the effect of magnification bias on is at most of order . In order to accommodate the discrepancy of between our result and the Dunlop & Peacock distribution of at , we would need an absurd slope, .

Distinguishing among the remaining three possibilities is harder. We believe possibility number 2 is unlikely because the discrepancy between Dunlop & Peacock and our work occurs at low redshift where their luminosity function should be most reliable: this regime is constrained by the local source counts rather than by extrapolation. Redshift-dependent bias (possibility number 3) exists for most samples of objects and there is no reason to expect it to be absent for NVSS. However, based on the Dunlop & Peacock and our , the bias would have to change from at to at . Such a large variation, combined with the unusually low value of the bias at , suggests that this is not the full explanation. The final possibility (4) is the removal of extended sources. This is hard to assess because of the low density of extended NVSS sources above our flux cut (deg). Of the 20 such sources in the COSMOS field, 19 match to VLA-COSMOS and 13 of these matches are found in the COSMOS optical/NIR catalog. It is worth noting that 8 of these (62%) have , versus 30/64 (30%) for the unresolved NVSS sources. This appears to go in the right direction, however it is difficult to make quantitative statements about whether the extended sources actually resolve the discrepant redshift distributions because of the unknown (but probably large, especially for the low- part of the distribution) sampling variance error bars.

In summary, while the full explanation for the difference between our and that of Dunlop & Peacock remains unknown, it seems likely (based on process of elimination) that a combination of redshift-dependent bias and our rejection of the unresolved NVSS sources plays a role. Magnification bias is ruled out as the explanation, and the discrepancy occurs in a regime where the extrapolations used in Dunlop & Peacock probably do not matter.

## Vi Systematics

We investigate various systematic effects in our correlations utilizing a specific multipole range. We choose these multipole bins based on two criteria. First, they should not be affected by non-linearities. Second, they should not be affected by any of the systematic effects in a significant way. We therefore only utilize the multipoles corresponding to Mpc and we also discard the first -bin for all samples since it is affected by the galactic foreground contamination. The specific -bins that are utilized are tabulated in Table 5.

### vi.1 Dust Extinction

Since it is possible that incorrect dust extinction systematically adds signals to our ISW cross correlation, we cross correlate the reddening maps Schlegel et al. (1998) in the same manner as we cross correlate each of our sample to the cosmic microwave background. If there is a systematic effect contributed via dust extinction, it will show up as a correlation, we can then estimate the effect and correct it from our tracer-cmb correlation.

In order to the verify that dust extinction does not affect our results, we constructed a vector of the estimated spurious cross-spectra . The spurious cross-spectra were computed by taking the cross-power spectrum of the CMB with the reddening map and multiplying by an estimate of . Note that has an entry for each -bin for each sample, so it has a total length of 42. We then compute the quantity (the derivation of this quantity and its relevance to understand contamination from extinction is detailed in Appendix C):