An Updated Tomographic Analysis of the Integrated SachsWolfe Effect and Implications for Dark Energy
Abstract
We derive updated constraints on the Integrated SachsWolfe (ISW) effect through crosscorrelation of the cosmic microwave background with galaxy surveys. We improve with respect to similar previous analyses in several ways. First, we use the most recent versions of extragalactic object catalogs: SDSS DR12 photometric redshift (photo) and 2MASS Photo datasets, as well as employed earlier for ISW, SDSS QSO photo and NVSS samples. Second, we use for the first time the WISE SuperCOSMOS catalog, which allows us to perform an allsky analysis of the ISW up to . Third, thanks to the use of photos, we separate each dataset into different redshift bins, deriving the crosscorrelation in each bin. This last step leads to a significant improvement in sensitivity. We remove crosscorrelation between catalogs using masks which mutually exclude common regions of the sky. We use two methods to quantify the significance of the ISW effect. In the first one, we fix the cosmological model, derive linear galaxy biases of the catalogs, and then evaluate the significance of the ISW using a single parameter. In the second approach we perform a global fit of the ISW and of the galaxy biases varying the cosmological model. We find significances of the ISW in the range 4.75.0 thus reaching, for the first time in such an analysis, the threshold of 5 . Without the redshift tomography we find a significance of 4.0 , which shows the importance of the binning method. Finally we use the ISW data to infer constraints on the Dark Energy redshift evolution and equation of state. We find that the redshift range covered by the catalogs is still not optimal to derive strong constraints, although this goal will be likely reached using future datasets such as from Euclid, LSST, and SKA.
I Introduction
We have, at present, strong evidence for Dark Energy (DE) from the large amount of available cosmological data (e.g., Planck Collaboration et al., 2016a). Nonetheless, this evidence is mostly based on precise constraints from the Cosmic Microwave Background (CMB) epoch extrapolated to the present time. Local, or presentday, constraints on DE are, instead, mostly given by SuperNovae (SN) data, which are not yet precise enough for accurately constraining the properties and time evolution of DE (e.g., Betoule et al., 2014).
Thus, it is important to look for alternative local DE probes. In this respect such a DEsensitive mesurement is given by the latetime Integrated SachsWolfe effect (ISW) on the CMB (Sachs and Wolfe, 1967). This effect is imprinted in the angular pattern of the CMB in the case of a nonflat universe, or in the presence of DE for a flat one. The effect is very small and cannot be well measured using the CMB alone since it peaks at large angular scales (small multipoles, ) which are cosmicvariance limited. On the other hand, it was realized that this effect can be more efficiently isolated by crosscorrelating the CMB with tracers of the LargeScale Structure (LSS) of the Universe at low () redshift (Boughn et al., 1998; Boughn and Crittenden, 2002), with most of the signal lying in the range for a standard CDM cosmological model (Afshordi, 2004).
In the past, many ISW analyses were performed using a large variety of tracers at different redshifts (Dupe et al., 2011; Planck Collaboration et al., 2014, 2016b; Afshordi et al., 2004; Boughn and Crittenden, 2004; Fosalba et al., 2003; Cabré et al., 2007, 2006; Giannantonio et al., 2006; Raccanelli et al., 2008; Rassat et al., 2007; Xia et al., 2009; Ferraro et al., 2015; HernandezMonteagudo et al., 2014; Shajib and Wright, 2016; Pietrobon et al., 2006; McEwen et al., 2007; Vielva et al., 2006). In a few cases, global analyses were performed combining different LSS tracers, giving the most stringent constraints and evidence for the ISW effect at the level of (Giannantonio et al., 2008, 2012; Ho et al., 2008). A further idea which was sometimes exploited is to use the redshift information of a given catalog to divide it into different redshift bins, compute the crosscorrelation in each bin, and then combine the information. This tomographic approach was pursued, for example, in the study of 2MASS (Francis and Peacock, 2010) or SDSS galaxies (Scranton et al., 2003; Sawangwit et al., 2010). Typically, the use of tomography does not provide strong improvement over the nobinning case, either because the catalog does not contain a large enough number of objects and splitting them increases the shotnoise, or because the redshift range is not well suited for ISW studies.
Nonetheless, in the recent years, several catalogs with redshift information and with a very large number of objects have become available thanks to the use of photometric redshifts (photos) instead of spectroscopic ones. Although photos are not as accurate as their spectroscopic counterparts, the former are sufficient for performing a tomographic analysis of the ISW with coarse bins. Hence we can exploit these large catalogs, which have the advantage of giving a low shot noise even when divided into subsamples. In this work, we combine for the first time the two above approaches: we use several datasets covering different redshifts ranges, and we bin them into redshift subsamples to perform a global tomography. We show that in this way we are able to improve the significance of the ISW effect from without redshift binning to exploiting the full tomography information. When combining the various catalogs, we take special care to minimize their overlap both in terms of common sources and the same LSS traced, in order not to use the same information many times. This is done by appropriate data cleaning and masking. We then use these improved measurement of the ISW effect to study deviations of DE from the simplest assumption of a cosmological constant.
Finally, the correlation data derived in this work and the associated likelihood will soon be made publicly available, in the next release of the MontePython^{1}^{1}1See http://baudren.github.io/montepython.html package (Audren et al., 2013).
Ii Theory
The expression for the crosscorrelation angular power spectrum (CAPS) between two fields and is given by:
(1) 
where is the presentday power spectrum of matter fluctuations. In the above expression we have assumed an underlying cosmological model, like CDM, in which the evolution of density fluctuations is separable in wavenumber and redshift on linear scales. A different expression applies, for example, in the presence of massive neutrinos (Lesgourgues et al., 2008), where the and evolution is not separable,
For the case of the fluctuation field of a catalog of discrete objects, one has
(2) 
where and represent the redshift distribution and the galaxy bias factor of the sources, respectively, are spherical Bessel functions, is the linear growth factor of density fluctuations and is the comoving distance to redshift .
For the case of crosscorrelation with the temperature fluctuation field obtained from the CMB maps (), the ISW effect in real space is given by (e.g., Nishizawa, 2014)
(3) 
where represents the gravitational potential. In the expression, we neglect a factor of , which introduces an error of the order of , smaller than the typical accuracy achieved in the determination of the ISW itself. Using the Poisson and Friedmann equations,
(4) 
where is the velocity of light, is the cosmological scale factor, is the Hubble parameter today, is the fractional density of matter today, and is the matter fluctuation field in Fourier space, we can write
(5) 
Finally, the equations above can be combined through Eq. (1) to give the CAPS expected for the ISW effect resulting from the correlation between a catalog of extragalactic objects, tracing the underlying mass distribution, and the CMB. Using the Limber approximation (Limber, 1953) the correlation becomes (Ho et al., 2008)
(6) 
The Limber approximation is very accurate at and accurate at the level of at (Limber, 1953), which is sufficient for the present analysis.
In our study, we use the public code class^{2}^{2}2See http://classcode.net (Blas et al., 2011) to compute the linear power spectrum of density fluctuations. As an option, this code can compute internally the spectra and , for arbitrary redshift distribution functions, using either the Limber approximation or a full integral in space. We prefer, nonetheless, to use the Limber approximation since CAPS calculations are significantly faster. Also, to get better performances and more flexibility, we choose to perform these calculations directly inside our python likelihood, reading only from the class output. We checked on a few examples that our spectra do agree with those computed internally by class.
Iii CMB maps
We use CMB maps from the Planck 2015 data release^{3}^{3}3See http://pla.esac.esa.int/pla/#maps (Planck Collaboration et al., 2016a) which have been produced using four different methods of foreground subtraction: Commander, NILC, SEVEM, and SMICA. Each method provides a confidence mask which defines the region of the sky in which the CMB maps can be used. We construct a combined mask as the union of these four confidence masks. This mask is applied on the CMB maps before calculating the crosscorrelation. We will use the SEVEM map as default for the analysis. Nonetheless, we have also tested the other maps to check the robustness of the results. The test is described in more detail in Sec. VIII.
As the ISW effect is achromatic, for further crosschecks we also use CMB maps at different frequencies. In particular we use maps at 100 GHz, 143 GHz, and 217 GHz. The results using these maps are also described in Sec. VIII.
Iv Additional cosmological datasets
In the following we will perform parameter fits using the ISW data obtained with the crosscorrelation. Beside this, in some setups, we will also use other cosmological datasets in conjunction. In particular, we will employ the Planck 2015 public likelihoods^{4}^{4}4See http://pla.esac.esa.int/pla/#cosmology (Planck Collaboration et al., 2016a) and the corresponding MontePython interfaces Planck_highl_lite (for high temperature), Planck_lowl (for low temperature and polarization), and Planck_lensing (CMB lensing reconstruction). Finally we will use BAO data from 6dF (Beutler et al., 2011), SDSS DR7 (Ross et al., 2015) and BOSS DR10&11 (Anderson et al., 2014), which are implemented as bao_boss and bao_boss_aniso in MontePython.
V Catalogs of Discrete Sources
For the crosscorrelation with the CMB, as tracers of matter distribution we use five catalogs of extragalactic sources. As the ISW is a wideangle effect, they were chosen to cover as large angular scales as possible, and two of them are allsky. Furthermore, our study does not require exact, i.e. spectroscopic, redshift information, thus photometric samples are sufficient. Except for one case, the datasets employed here include individual photos for each source, which allows us to perform a tomographic approach by splitting the datasets into redshift bins.
The catalogs we use span a wide redshift range; see Fig. 1 for their individual redshift distributions. Table 1 quantifies their properties (sky coverage, number of sources, mean projected density) as effectively used for the analysis, i.e., after applying both the catalog and CMB masks.
For a plot of the sky maps and masks of the catalogs described below, and for their detailed description, see Cuoco et al. (2017). Below we provide a short summary of the properties of the datasets.
source  sky  number  mean surface 
catalog  coverage  of sources  density [deg] 
NVSS  62.3%  431,724  67.2 
2MPZ  64.2%  661,060  24.9 
WISESCOS  64.5%  17,695,635  665 
SDSS DR12  18.7%  23,907,634  3095 
SDSS DR6 QSO  15.6%  461,093  71.8 
v.1 2mpz
As a tracer of the most local LSS in this study we use the 2MASS Photometric Redshift catalog^{5}^{5}5Available from http://ssa.roe.ac.uk/TWOMPZ.html. (2MPZ, Bilicki et al., 2014). This dataset was built by merging three allsky photometric datasets covering optical, nearinfrared (IR), and midIR passbands: SuperCOSMOS scans of UKST/POSSII photographic plates (Peacock et al., 2016), 2MASS Extended Source Catalog (Jarrett et al., 2000), and Widefield Infrared Survey Explorer (WISE, Wright et al., 2010). Photos were subsequently estimated for all the included sources, by calibrating on overlapping spectroscopic datasets.
2MPZ includes galaxies over almost the full sky. Part of this area is however undersampled due to the Galactic foreground and instrumental artifacts, we thus applied a mask described in Alonso et al. (2015). When combined with the CMB mask, this leaves over 660,000 2MPZ galaxies on of the sky (Table 1).
2MPZ provides the bestconstrained photos among the catalogs used in this paper. They are practically unbiased () and their random errors have RMS scatter , to a good accuracy independent of redshift. We show the 2MPZ redshift distribution in Fig. 1 with the dotdashed green line; the peak is at while the mean . The overall surface density of 2MPZ is sources per square degree.
For the tomographic analysis we split the catalog in three redshift bins: , and . The first two include the bulk of the distribution, approximately divided into two comparable subsamples, while the third bin explores the tail of the where most of the ISW signal is expected.
A precursor of 2MPZ, based on 2MASS and SuperCOSMOS only, was used in a tomographic ISW analysis by Francis and Peacock (2010), while an early application of 2MPZ itself to ISW tomography is presented in Steward (2014). In both cases no significant ISW signal was found, consistent with expectations. Another ISWrelated application of 2MPZ is presented in Planck Collaboration et al. (2016b), where it was applied to reconstruct ISW anisotropies caused by the LSS.
v.2 Wise SuperCOSMOS
The WISE SuperCOSMOS photo catalog^{6}^{6}6Available from http://ssa.roe.ac.uk/WISExSCOS.html. (WISC, Bilicki et al., 2016) is an allsky extension of 2MPZ obtained by crossmatching WISE and SuperCOSMOS samples. WISC reaches roughly 3 times deeper than 2MPZ and has almost 30 times larger surface density. However, it suffers from more severe foreground contamination, and its useful area is of the sky after applying its default mask. This is further reduced to once the Planck mask is also used; the resulting WISC sample includes about 17.5 million galaxies.
WISC photos have overall mean error and distancedependent scatter of . The redshift distribution is shown in Fig. 1 with the dashed orange curve. The peak is at , and the majority of the sources are within . In the tomographic approach, the WISC sample is divided into four redshift bins: , , , and , with approximately equal number of galaxies in each bin.
As far as we are aware, our study employs the WISC dataset for an ISW analysis for the first time. Various studies based on WISE have been performed in the past (Goto et al., 2012; Kovács et al., 2013; Ferraro et al., 2015; Shajib and Wright, 2016). However, the samples used there differed significantly from WISC, and none included individual redshift estimates which would allow for redshift binning.
v.3 SDSS DR12 photometric
Currently there are no allsky photo catalogs available reaching beyond WISC. Therefore, in order to look for the ISW signal at , we used datasets of smaller sky coverage. The first of them, with the largest number density of all employed in this paper, is based on the Sloan Digital Sky Survey Data Release 12 (SDSSDR12) photo sample compiled by Beck et al. (2016); to our knowledge, our study is its first application to an ISW analysis, although earlier versions (DR 6 and DR 8) were used in Giannantonio et al. (2008, 2012) (but without binning).
The parent SDSSDR12 photo dataset includes over 200 million galaxies. Here we however use a subsample described in detail in Cuoco et al. (2017), which was obtained via appropriate cleaning as recommended by Beck et al. (2016), together with our own subsequent purification of problematic sky areas. In particular, as the SDSS galaxies are distributed in two disconnected regions in the Galactic south and north, with most of the area in the northern part, and uneven sampling in the south, we have excluded the latter region from the analysis. After additionally employing the Planck CMB mask, we were left with about 24 million SDSS DR12 sources with mean and mostly within . The resulting sky coverage is and the mean surface density is deg. The redshift distribution is shown in Fig. 1 with the solid blue line.
Thanks to the very large projected density of objects, we were able to split the SDSSDR12 sample into several redshift bins, keeping low shotnoise in each shell. For the tomographic analysis we divided the dataset into six bins: , , , , and . The range is not subdivided further since this redshift range is best covered by WISC, where we already have subbins. The photo accuracy of SDSSDR12 depends on the ‘photo class’ defined by Beck et al. (2016), and each class has an associated error estimate. Our specific preselection detailed in Cuoco et al. (2017) leads to an effective photo scatter of based on the overall error estimates from Beck et al. (2016).
v.4 Sdss Dr6 Qso
As a tracer of high LSS, we use a catalog of photometric quasars (QSOs) compiled by Richards et al. (2009) from the SDSS DR6 dataset (DR6QSO in the following), used previously in ISW studies by e.g. Giannantonio et al. (2008); Xia et al. (2009) and Giannantonio et al. (2012). We apply the same preselections as in Xia et al. (2009), and the resulting sample includes QSOs on of the sky. We exclude from the analysis three narrow stripes present in the south Galactic sky and use only the northern region.
The DR6QSO sources are provided with photos spanning formally but with a relatively peaked and mean (dotted red line in Fig. 1). For tomographic analysis, this QSO dataset will be split into three bins of , , and , selected in a way to have similar number of objects in each bin. We excluded the QSOs in the range in order to minimize the overlap with the other catalogs in this redshift range. Nonetheless, there are very few DR6QSO catalog objects at these redshifts, thus this choice has only a very minor impact on the results. The typical photo accuracy of this dataset is as reported by Richards et al. (2009), and we will use this number for the extended modeling of underlying s per redshift bin in Sec. VIII.
v.5 Nvss
The NRAO VLA Sky Survey (NVSS, Condon et al., 1998) is a catalog of radio sources, most of which are extragalactic. This sample has already been used for multiple ISW studies (e.g. Boughn and Crittenden, 2002; Pietrobon et al., 2006; Vielva et al., 2006; McEwen et al., 2007; Raccanelli et al., 2008). The dataset covers the whole sky available for the VLA instrument; after appropriate cleanup of likely Galactic entries and artifacts, the NVSS sample includes objects fluxlimited to mJy, located at declinations and Galactic latitudes . This is the only of the datasets considered in this work which does not provide even crude redshift information for the individual sources. We thus use it without tomographic binning and, where relevant, assume its to follow the model of de Zotti et al. (2010) (purple shortlongdashed line in Fig. 1). This sample spans the broadest redshift range of all the considered catalogs, namely .
v.6 Masks
In the correlation of the CMB with each catalog we use the CMB mask, described in Sec. III, combined with the specific mask of the given catalog. Beside this, we define specific masks which we use when combining the signal from the different catalogs in order to circumvent including the same information twice, and to avoid the need to take into account the crosscorrelations between various tracers of the same LSS. We proceeded in the following way.

SDSS catalogs (i.e. SDSS DR6 QSOs and SDSS DR12 galaxies) are used without additional masks. When combining the information with other catalogs we, however, exclude the first SDSS DR12 bin, since the region is best covered by 2MPZ.

To avoid correlations with the SDSS catalogs, when using all the remaining ones (i.e. NVSS, 2MPZ, WISC) we apply a mask which is a complementary of the joint mask of SDSS DR12 galaxies and SDSS DR6 QSOs (in short, SDSS mask in the following).

For 2MPZ and WISC it is not possible to define mutually exclusive masks since both these datasets cover practically the same part of the sky. Nonetheless, we use them together, since WISC was built excluding most of the objects already contained in 2MPZ (Bilicki et al., 2016). The two catalogs, thus, have practically no common sources. In this way the correlation among the two datasets is significantly suppressed, although not totally, since both trace the same underlying LSS in the overlapping redshift ranges. We will, however, not consider the first bin, , of WISC in the combined analysis since in this redshift range 2MPZ has better redshift determination and basically no stellar contamination. Nonetheless, as we will show in Sec. VII, the evidence for ISW in the range , where 2MPZ and WISC have most of the overlap, is very small, so, in practice, this has only a marginal effect on the final ISW significance.

Similarly, also for NVSS, 2MPZ and WISC it is not possible to define a mutually exclusive mask due to the large common area of the sky. In this case, we note that 2MPZ and WISC cover only the low redshift tail of NVSS. Thus, the overlap and correlation among them is minimal.
We will thus use the above setup when reporting combined significances of the ISW from the different catalogs. For simplicity, we will use the same setup also to derive autocorrelations of the single catalogs. In this case the significances could be increased slightly for NVSS, 2MPZ, and WISC if their proper masks were used, but we checked that the improvement is only marginal.
Vi CrossCorrelation Analysis
In the previous section we have presented the catalogs of extragalactic objects that we use in the analysis. Their input format is that of a 2D pixelized map of object counts , where specifies the angular coordinate of the th pixel. For the crosscorrelation analysis we consider maps of normalized counts , where is the mean object density in the unmasked area, and CMB temperature maps, also pixelized with a matching angular resolution.
In our analysis we compute both the angular 2point crosscorrelation function, CCF, , and its harmonic transform, the angular power spectrum , CAPS. However, we restrict the quantitative analysis to the CAPS only. The reason for this choice is that the CAPS has the advantage that different multipoles are almost uncorrelated, especially after binning. Their covariance matrix is therefore close to diagonal, which simplifies the comparison between models and data. Similarly, we compute also the autocorrelation power spectrum of the catalogs (APS) and the related autocorrelation function (ACF).
catalog  z  b  
SDSS  00.1  3.59  3.76  4.11  
0.10.3  1.71  1.68  1.63  
0.30.4  0.64  0.63  0.61  
0.40.5  4.84  4.65  4.99  
0.50.7  6.16  5.86  6.35  
0.71  15.04  14.99  15.16  
WiXSC  00.09  0.46  0.38  0.28  
0.090.21  2.38  2.67  2.62  
0.210.3  10.07  10.14  10.09  
0.30.6  5.62  5.88  5.55  
QSO  01  5.9  5.93  4.97  
0.51  3.09  3.07  3.07  
12  3.61  3.59  3.6  
23  7.08  7.05  7.08  
2MPZ  00.105  4.41  1.30  1.26  
0.1050.195  2.00  2.07  2.17  
0.1950.3  6.54  6.67  6.34  
NVSS  06  3.02  0.64  —  —  
catalog  z  b  
SDSS  01  1.25  1.27  1.59  
WiXSC  00.6  3.15  3.74  3.94  
QSO  03  2.77  2.76  2.4  
2MPZ  00.3  5.19  2.04  1.98 
We use the PolSpice^{7}^{7}7See http://www2.iap.fr/users/hivon/software/PolSpice/ statistical toolkit (Szapudi et al., 2001; Chon et al., 2004; Efstathiou, 2004a; Challinor and Chon, 2005) to estimate the correlation functions and power spectra. PolSpice automatically corrects for the effect of the mask. In this respect, we point out that the effective geometry of the mask used for the correlation analysis is obtained by combining that of the CMB maps with those of each catalog of astrophysical objects. The accuracy of the PolSpice estimator has been assessed in Xia et al. (2015) by comparing the measured CCF with the one computed using the popular LandySzalay method (Landy and Szalay, 1993). The two were found to be in very good agreement. PolSpice also provides the covariance matrix for the angular power spectrum, (Efstathiou, 2004b).
For the case of source catalog APS a further step is required. Contrary to the CAPS, the APS contains shot noise due to the discrete nature of the objects in the map. The shot noise is constant in multipole and can be expressed as , where is the fraction of sky covered by the catalog in the unmasked area and is the number of catalog objects, again in the unmasked area. The above shotnoise has been subtracted from our final estimated APS.
The Planck Point Spread Function and the map pixelization affect in principle the estimate of the CAPS. However, the CAPS contains information on the ISW only up to where these effects are negligible. We will thus not consider them further.
Finally, to reduce the correlation in nearby multipoles induced by the angular mask, we use an binned version of the measured CAPS. The number of bins and the maximum and minimum used in the analysis will be varied to assess the robustness of the results. We indicate the binned CAPS with the same symbol as the unbinned one, . It should be clear from the context which one is used. The in each bin is given by the simple unweighted average of the within the bin. For the binned we build the corresponding covariance matrix as a block average of the unbinned covariance matrix , i.e., , where are the widths of the two multipole bins, and run over the multipoles of the first and the second bin. The binning procedure is very efficient in removing correlation among nearby multipoles, resulting in a block covariance matrix that is, to a good approximation, diagonal. We will use nonetheless the full block covariance matrix in the following, although we have checked that using the diagonal only gives minor differences. When showing CAPS plots, however, we use the diagonal terms to plot the errors on the , , where the sum runs over the multipoles of the bin contributing to .
catalog  z  

SDSS  00.1  0.07  1.224  1.219  0.005  
0.10.3  0.87  3.89  3.12  0.76  
0.30.4  1.57  4.47  2.01  2.45  
0.40.5  2.03  6.57  2.45  4.12  
0.50.7  2.28  9.28  4.06  5.22  
0.71  0.37  6.76  6.62  0.13  
WiXSC  00.09  1.08  2.84  1.68  1.16  
0.090.21  0.33  4.63  4.52  0.11  
0.210.3  1.1  3.62  2.4  1.21  
0.30.6  1.41  4.91  2.92  1.99  
QSO  01  1.52  5.95  3.64  2.31  
0.51  1.45  7.46  5.34  2.11  
12  1.52  3.99  1.68  2.31  
23  0.38  3.11  2.96  0.14  
2MPZ  00.105  0.36  1.26  1.13  0.13  
0.1050.195  0.3  1.12  1.03  0.09  
0.1950.3  0.71  1.66  1.16  0.5  
NVSS  06  2.97  14.9  6.11  8.79 
catalog  

SDSS  3.29  30.96  20.11  8.46  
WixSC  1.67  13.16  10.39  2.76  
Quasars  2.13  14.55  10.01  2.99  
2MPZ  0.81  4.04  3.38  0.65  
SDSS+WixSC  3.49  44.12  31.94  11.21  
SDSS+Quasars  3.9  45.51  30.28  11.45  
SDSS+WiXSC+Quasars  4  58.67  42.66  14.2  
SDSS+WiXSC+Quasars+NVSS  4.97  73.57  48.85  21.52  
SDSS+WiXSC+Quasars+NVSS+2MPZ  5  77.61  52.61  22.16 
catalog  

SDSS  1.49  5.3  3.09  2.21  
WixSC  1.02  5.28  4.24  0.65  
Quasars  2.03  5.55  1.41  3.94  
2MPZ  0.39  0.87  0.72  0.15  
NVSS  2.97  14.9  6.11  8.79  
SDSS+WixSC  2.23  18.47  13.48  4.96  
SDSS+Quasars  2.35  19.85  14.33  5.2  
SDSS+WiXSC+Quasars  2.84  33.02  24.97  7.95  
SDSS+WiXSC+Quasars+NVSS  4.02  47.91  31.76  15.27  
SDSS+WiXSC+Quasars+NVSS+2MPZ  4.08  51.95  35.28  15.92 
Parameter 
68% limits 

Vii Derivation of the ISW significance
In this section we illustrate the two methods we use to quantify the significance of the ISW. We will assume for the first method a flat CDM model with cosmological parameters , , , , at Mpc, and , in accordance with the most recent Planck results (Planck Collaboration et al., 2016a).
vii.1 Method 1
This is the usual method employed in previous publications to study the significance of the ISW. In this case we fix the cosmological model to the bestfit one measured by Planck, and we derive with class the matter power spectrum , which is used to calculate the expected autocorrelation for each catalog for the appropriate redshift bin. The measured autocorrelation is then used to fit the linear bias, as a proportionality constant in the predicted . An example of this fit is shown in the left panel of Fig. 2. A simple over the bins of the autocorrelation is used for the fit:
(7) 
where and represent the model and the measured CAPS, and the sum is over all bins.
As mentioned in Sec. VI we tested that the use of the full covariance matrix with respect to the diagonal expression for the above does not give appreciable differences. Table 2 summarizes the various measured biases, and the default binning used for the autocorrelations. We tested the robustness of the fitted biases changing the number of bins from 4 to 6 and the maximum from 4080, and we found stable results, with variations of the order of 10%. A maximum of 4080 is chosen since above this range typically nonlinear effects become significant. As the default case, we use 4 bins in the range 1060.
As a further test we checked the impact of using nonlinear corrections to the matter power spectrum to model the autocorrelation of the catalogs. The nonlinear corrections were implemented through the version of Halofit (Takahashi et al., 2012) implemented in class v2.6.1. The last 2 columns of Table 2 show the bias and the bestfit obtained using the nonlinear model. It can be seen that the biases obtained with and without nonlinear corrections are fully compatible. The only exception is the first redshift bin of 2MPZ where the bestfit bias changes at the level. More importantly, the fit shows a visible improvement from to . This is expected, since at these low redshifts even the small s correspond mostly to small, nonlinear, physical scales. As we show below, however, 2MPZ presents little or no imprint of the ISW effect, so we conclude that the use of the linear has a negligible impact on the study of the ISW effect in this analysis.
In the second step, all the galaxy biases are fixed to bestfit values previously derived, and only the measured crosscorrelations are used. At this point only a single parameter is fitted using as data either a single measured crosscorrelation or a combination of them, with the statistics:
(8) 
where and represent the model (for the standard CDM cosmological model considered) and the measured catalog – CMB temperature crosscorrelation for a given redshift bin, respectively, the sum is over all the bins, and over different catalogs and different redshift bins. The linear parameter quantifies the agreement with the above standard model expectation. In the denominator we use the error provided by Polspice discussed in the previous section. In principle, however, one should use an error where the model is taken into account. For the case of binned data, however, this is a small effect (see for example discussion in Fornasa et al. (2016)).
An example of measured crosscorrelation and fit to the model is shown in the right panel of Fig. 2. Table 3 summarizes the results of the fit for each single bin of each catalog, for each catalog combining the different bins, and for different combinations of the catalogs, where, again, for each catalog binning has been used. For the default case we use four multipole bins between of 4 and 100, but, again, we have verified that the results are stable when changing the number of bins from 4 to 6 and the maximum from 60 to 100, which is expected, since the ISW effect is rapidly decreasing as a function of , and not much signal is expected beyond .
To quantify the significance of the measurement we use as teststatistic the quantity
(9) 
where is the minimum , and is the of the null hypothesis of no ISW effect, i.e. of the case . TS is expected to behave asymptotically as a distribution with a number of degrees of freedom equal to the number of fitted parameters, allowing us to derive the significance level of a measurement based on the measured TS. In this case, since there is only one fitted parameter, the significance in sigma is just given by . From Table 3 one can see that the maximum significance achieved with Method 1 when using all the catalogs in combination is . From the different results it can also be seen that the main contribution is given by NVSS and SDSS DR12 galaxies. We remind that the crosscorrelation with NVSS is calculated masking the area of the sky used to calculate the correlation with SDSS. The two are, thus, completely independent. A smaller, and comparable, contribution, is given by WISC and SDSSQSO. 2MPZ instead show basically no sign of ISW, which is expected given the very low range.
We also show in Table 4 the result of the fit when no redshift binning is used. It is clear that without such binning the significance of the ISW is significantly reduced, especially for SDSSDR12 and WISC , while the significance of SDSSQSO is almost unchanged. Overall, when no redshift binning is used, the significance of the ISW effect combining all the catalog is 4.0 , which is significantly reduced with respect to the 4.7 achieved with the redshift binning.
Parameter  AC+CC  CC  PL+AC+CC 

—  —  
—  —  
vii.2 Method 2
The first method is, in principle, not fully selfconsistent, because the autocorrelations, and hence the biases, are sensitive to the underlying matter power spectrum. We fixed the matter power spectrum to the Planck CDM bestfitting model, but this may not be the best fit to the autocorrelation data. The induced error should be negligible when CMB and BAO are also used, since they impose to be very close to the fiducial model. But more importantly, the crosscorrelation determines a given amount of ISW, and this has in principle an effect on cosmology, since a different ISW means a different Dark Energy model and thus also a different . For these reasons it is more consistent to fit to the data at the same time as the bias parameters, the cosmological parameters, and the parameter used to assess the detection significance.
We perform such a fit using the MontePython environment. The fit typically involves many parameters () which can present degeneracies which are not known in advance. To scan efficiently this parameter space we run MontePython in the Multinest mode (Feroz et al., 2009) . In this way we can robustly explore the posterior with typically likelihood evaluations, and efficiencies of the order of 10%. We consider two cases.
In the first case, we only use crosscorrelation and autocorrelation measurements. We call this dataset AC+CC, and we fit a total of 22 parameters, i.e, 15 biases, , and the six CDM parameters (, , , , , ). When Planck data are used, we also include the nuisance parameter (Planck Collaboration et al., 2016a). For all cosmological parameters except , we use Gaussian priors derived from a fit of Planck+BAO summarized in Table 5, which are consistent with those published in Planck Collaboration et al. (2016a). The error bars from Planck+BAO are so small that we find essentially the same result for when fixing these five parameters to their best fit values instead of marginalizing over them with Gaussian priors. Our results for this fit are shown in the first column of Table 6. As expected, the constraint on coming from the AC+CC data is weaker than that from Planck+BAO data, by about a factor 6. Also, the bestfit of the AC+CC analysis is lower than the Planck+BAO fit, by about 2. The fitted galaxy biases are typically compatible with those of Method 1, although in several cases they are 1020% larger, which can be understood as a consequence of the lower , resulting in a lower normalization. Indeed, the measured autocorrelations basically fix the product of the squared biases and of the overall amplitude. Comparing the case with free to the one with , we find TS== 22, giving a significance of 4.7 , identical to the one found in Method 1.
In the second case, we fit the same parameters to the data, but we now include the full Planck+BAO likelihoods instead of Gaussian priors on five parameters. Formally, we use the Planck and BAO likelihoods combined with the from the AC+CC data:
(10) 
It should be noted that the use of other data besides AC+CC does not affect the ability to derive the significance of the ISW detection, which is only encoded in the parameter entering the AC+CC likelihood. Results of this fit are shown in the third column of Table 6. The main difference with respect to the previous fit is the value of , now driven back to the Planck bestfit. This upward shift in results, again, in a global downward shift of the biases, by about 1020%, giving now a better compatibility with the results of Method 1.
In general, apart from the small degeneracy with resolved by the inclusion of Planck+BAO data, the biases are well constrained by the fit. This means that the subspace of biases is approximately orthogonal to the rest of the global parameter space, which simplifies the fit and speeds up its convergence. To measure the significance, in this case we define the test statistic as TS= , which shares the same properties of the TS defined in terms of the . Comparing the case with free to the one with , we now get TS= = 24.9, which gives a significance of 5.0 . Since the cosmology is basically fixed by the Planck+BAO data to a point in parameter space very close to the fiducial model of Method 1, this improvement in significance comes, apparently, from fitting jointly the biases and (while in Method 1 the biases were kept fixed using the results of the first step of the method). The joint fit explores the correlations which exist between the biases and . This results in a better global fit, and also in a slightly enhanced significance, reaching the 5 threshold.
Viii Robustness tests
In this section we describe some further tests performed to verify the robustness of the results.
As mentioned in Sec. III, several CMB maps are available from Planck, resulting from different foreground cleaning methods. In Fig. 3 we show the results of the crosscorrelation using four CMB maps cleaned with four different methods. We pick up as an example the crosscorrelation with the full 2MPZ catalog, without subdivision in redshift bins. It clearly appears that the use of different maps has no appreciable impact on the result.
Another important aspect is the possible frequency dependence of the correlation. In particular, while the ISW effect is expected to be achromatic, some secondary effects, like a correlation due to a SunyaevZel’dovich (Sunyaev and Zeldovich, 1972) or ReesSciama (Rees and Sciama, 1968) imprint in the CMB map, are expected to be frequency dependent. To test this possibility, we use available Planck CMB maps at 100 GHz, 143 GHz and 217 GHz. Again the the full 2MPZ catalog is used as example, since these effects are expected to peak at low redshift. Fig. 4 shows the result of the correlation at different frequencies. We observe a very small trend of the CAPS with frequency, especially for the first bin, but this effect is negligible with respect to the error bars of the data points. Results are similar for the other catalogs, showing no frequency dependence.
Finally, we tested the effect of photo errors. In the basic setup, the theoretical predictions for the auto and crosscorrelation functions per redshift bin are modeled by assuming that the true redshift distribution is well approximated by the photo one, i.e. . In reality, sharp cuts in will correspond to more extended tails in because the photos are smeared out in the radial direction. However, we can easily take photo errors into account if we know their statistical properties. In the case of 2MPZ, the photo error is basically constant in and has roughly Gaussian scatter of centered at , while for WISC the scatter is with also approximately zero mean in . For SDSS QSOs it is also approximately constant in and equal to 0.24. Finally, for SDSS DR12 the error is (see Sec. V). We thus derive the effective true redshift distribution of a given bin by convolving the measured photo selection function in that bin with a dependent Gaussian of width . The resulting true distribution is a smoothed version of the photo distribution, presenting tails outside the edges of the bin. We then use this distribution to fit again the auto and crosscorrelations data. The results are shown in the last column of Table 2. We find that the effect of photo errors has some impact on the determination of the biases. The effect is most important in the high tails of various catalogs, and, in particular, WISC and SDSS DR12. This is not surprising since, in these cases, the photo errors increase with redshift and are largest at high. The effect is at the level of 1020%. This corresponds to a decrease in of the same amount in these bins. Nonetheless, since the above bins only have a limited weight on the combined fit, the impact on the final determined from the global fit of all bins and catalogs is basically negligible.
Parameter  68% limits 














Ix Dark energy fit
In this section we investigate the power of the crosscorrelation data to constrain DE, in a similar framework as presented in Corasaniti et al. (2005); Pogosian et al. (2005) and Pogosian (2005). For this purpose, we do not use the parameter employed in Sec. VII, since it is only an artificial quantity necessary to evaluate the ISW significance from the crosscorrelation data. However, as shown in Sec. VII, there is indication that the bestfit value of is above 1 at slightly more than 1. This suggests (although with low statistical significance) that DE could differ from a simple cosmological constant. To investigate this more in detail, we perform a fit with Method 2 of Sec. VII, but with , and with extra parameters accounting for dynamical Dark Energy. For simplicity, we use the empirical parametrization (Linder, 2003) and the parameterized postFriedmann framework of Hu and Sawicki (2007) and Fang et al. (2008), which are implemented in class, to study models with . We test several different fit setups. In particular, since the AC dataset is a cosmological probe with its own sensitivity to the cosmological parameters, we test various combinations in which the AC and CC data are used separately. A further reason to study the AC data separately from the CC ones is that the APS of extragalactic objects are typically difficult to model accurately, even at small , due to the nonlinearity and possible stochasticity of the galaxy bias with respect to matter. Separate fits to the AC and CC data could then reveal inconsistencies that might be associated to our minimal assumption that the bias is linear and scaleindependent.
We perform the following fits: (a) Planck+BAO, (b) Planck+BAO+CC+AC, (c) Planck+BAO+CC, (d) CC only, (e) AC only, (f) AC+CC. Case (a) has the standard 6 CDM parameters, plus , , and one Planck nuisance parameter, , required for the evaluation of the Planck likelihood (Planck Collaboration et al., 2016a), thus 9 parameters in total. The results of this baseline fit are shown in Table 7. Case (b) includes CC and AC datasets, and uses additionally 15 bias parameters (24 parameters in total). Case (c) is similar to (b) but without AC data. Since the biases are still needed for the CC fitting, they are still included in the fit, but with a Gaussian prior coming from fit (b). We verified that just fixing the biases to the best fit (b), instead of including them in the fit with Gaussian priors, does not actually change the results. Similarly, the result does not change if the biases are taken from another fit than (b), like (e) or (f). For fit (d), featuring only CC data, all cosmological parameters except (, , ) and all bias parameters are either fixed or marginalized with Gaussian priors. For fit (e), featuring only AC data, all cosmological parameters except (, , ) are fixed to bestfit values, while the biases are left free, since they are constrained by the AC data. Finally fit (f) combines AC and CC data, and uses the same setup as fit (e).
Parameter  AC+CC  CC+bias priors  AC  PL+AC+CC  PL+CC+bias priors 

—  —  —  
—  —  —  
Figs. 56 show the results for and (marginalized over all the remaining parameters) for some of these fits. Table 8 gives the confidence intervals on all the parameters for all our fits. The most evident result is that the AConly fit selects a region of parameter space significantly in tension with the Planck+BAO constraints, basically excluding the standard case at more than 3 . This is very likely a consequence of the linear bias model not being accurate enough to provide reliable cosmological constraints. This might be particularly true for the autocorrelation of the catalogs in the highest redshift bins, which are the most sensitive to deviations from a standard cosmological constant, but also the ones lying in the tail of the redshift distribution of the catalog, where different population of galaxies are probably selected, which requires more accurate modeling. More sophisticated approach to the modeling of the catalog autocorrelations might be thus required to address properly this issue. Various bias models have been proposed beyond linear bias, like for instance models based on the halo occupation distribution of the catalog objects (see for instance Ando et al. (2017)). We leave a systematic study of this subject for future work. Another possible issue can be related to intrinsic artifacts in the catalog, like nonuniformity in the sky coverage, or large errors in the photo determination. These problems can become more evident especially in the tails of the redshift distribution. Indeed, the largest for AC fits from Table 2 are for the bins in the tail of the distribution, especially for SDSS DR12 and QSOs, indicating a poor match between the model and the data. This can be seen more explicitly also in the related plots in Appendix A.
Hence, in deriving DE constraints it is more conservative to discard information from AC and focus on CC only. We see that the constraints from the CC data are compatible with Planck+BAO results. However, given the relatively low significance of the ISW effect, the former are about three times weaker than the latter for each parameter. The direction of the degeneracy between and is approximately the same in the two fits, which was not obvious a priori, since the two data sets are sensitive to Dark Energy through different physical effects. It appears that the valley of wellfitting models with always corresponds to crossing in the range , but with very different derivatives . Even when is very large, all models in this valley do feature accelerated expansion of the Universe in the recent past, but not necessarily today. In fact, when increases while decreases simultaneously, the stage of accelerated expansion is preserved but translated backward in time.
Since the CC data are less sensitive than Planck and do not feature a different direction of degeneracy, the joint constraints from Planck+BAO+CC are basically unchanged with respect to Planck+BAO only.
X Discussion and Conclusions
We derived an updated measurement of the ISW effect through crosscorrelations of the cosmic microwave background with several galaxy surveys, namely, 2MASS Photometric Redshift catalog (2MPZ), NVSS, SDSS QSOs, SDSS DR12 photometric redshift dataset, and WISE SuperCOSMOS; the two latter are here used for the first time for an ISW analysis. We also improved with respect to previous analyses performing tomography within each catalog, i.e., exploiting the photometric redshifts and dividing each catalog into redshift bins. We found that the current crosscorrelation data provide strong evidence for the ISW effect and thus for Dark Energy, at the 5 level.
However, current catalogs are still not optimal to derive cosmological constraints from the ISW, for two main reasons. First, the clustering of objects requires complicated modeling, probably beyond the simple linear bias assumption. On this last point, improvements are possible using more sophisticated modeling, but at a price of introducing more nuisance parameters. Also, the tails of the redshift distributions of the objects might be more strongly affected by catalog systematics such as uneven sampling or large photo errors.
Second, the data used in this paper are sensitive mostly to the redshift range , while the ISW effect is expected to be important for . Several planned or forthcoming wideangle galaxy surveys will cover this redshift range and should thus bring (major) improvement for ISW detection via crosscorrelation with CMB. For the Euclid satellite, the predicted significance of such a signal is (Amendola et al., 2016), and one should expect similar figures from the Large Synoptic Survey Telescope (LSST Science Collaboration et al., 2009), and the SquareKilometer Array (Jarvis et al., 2015). The very high S/N of ISW from these deep and wide future catalogs will not only allow for much stronger constraints on dark energy than we obtained here, but even on some modified gravity models which often predict very different ISW signatures than CDM (e.g. Renk et al., 2017).
acknowledgments
Simulations were performed with computing resources granted by RWTH Aachen University under project thes0263.
MB is supported by the Netherlands Organization for Scientific Research, NWO, through grant number 614.001.451, and by the Polish National Science Center under contract #UMO2012/07/D/ST9/02785.
Some of the results in this paper have been derived using the HEALPix package^{8}^{8}8http://healpix.sourceforge.net/ (Górski et al., 2005).
This research has made use of data obtained from the SuperCOSMOS Science Archive, prepared and hosted by the Wide Field Astronomy Unit, Institute for Astronomy, University of Edinburgh, which is funded by the UK Science and Technology Facilities Council.
Funding for SDSSIII has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSSIII web site is http://www.sdss3.org/. SDSSIII is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSSIII Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.
Some of the results in this paper have been derived using the GetDist package^{9}^{9}9https://github.com/cmbant/getdist.
References
 Planck Collaboration et al. (2016a) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, N. Bartolo, E. Battaner, R. Battye, K. Benabed, A. Benoît, et al., A&A 594, A13 (2016a), arXiv:1502.01589 .
 Betoule et al. (2014) M. Betoule, R. Kessler, J. Guy, J. Mosher, D. Hardin, R. Biswas, P. Astier, P. ElHage, M. Konig, S. Kuhlmann, J. Marriner, R. Pain, N. Regnault, C. Balland, B. A. Bassett, P. J. Brown, H. Campbell, R. G. Carlberg, F. CellierHolzem, D. Cinabro, A. Conley, C. B. D’Andrea, et al., A&A 568, A22 (2014), arXiv:1401.4064 .
 Sachs and Wolfe (1967) R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
 Boughn et al. (1998) S. P. Boughn, R. G. Crittenden, and N. G. Turok, New Astron. 3, 275 (1998), arXiv:astroph/9704043 .
 Boughn and Crittenden (2002) S. P. Boughn and R. G. Crittenden, Phys. Rev. Lett. 88, 021302 (2002), arXiv:astroph/0111281 [astroph] .
 Afshordi (2004) N. Afshordi, Phys. Rev. D D70, 083536 (2004), arXiv:astroph/0401166 [astroph] .
 Dupe et al. (2011) F.X. Dupe, A. Rassat, J.L. Starck, and M. J. Fadili, A&A 534, A51 (2011), arXiv:1010.2192 [astroph.CO] .
 Planck Collaboration et al. (2014) Planck Collaboration, P. A. R. Ade, N. Aghanim, C. ArmitageCaplan, M. Arnaud, M. Ashdown, F. AtrioBarandela, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, N. Bartolo, E. Battaner, and et al., A&A 571, A19 (2014), arXiv:1303.5079 .
 Planck Collaboration et al. (2016b) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, E. Battaner, K. Benabed, A. Benoît, A. BenoitLévy, J.P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, A. Bonaldi, and et al., A&A 594, A21 (2016b), arXiv:1502.01595 .
 Afshordi et al. (2004) N. Afshordi, Y. Loh, and M. A. Strauss, Phys. Rev. D 69, 083524 (2004), arXiv:astroph/0308260 .
 Boughn and Crittenden (2004) S. Boughn and R. Crittenden, Nature (London) 427, 45 (2004), arXiv:astroph/0305001 .
 Fosalba et al. (2003) P. Fosalba, E. Gaztanaga, and F. Castander, Astrophys. J. 597, L89 (2003), arXiv:astroph/0307249 [astroph] .
 Cabré et al. (2007) A. Cabré, P. Fosalba, E. Gaztañaga, and M. Manera, MNRAS 381, 1347 (2007), arXiv:astroph/0701393 .
 Cabré et al. (2006) A. Cabré, E. Gaztañaga, M. Manera, P. Fosalba, and F. Castander, MNRAS 372, L23 (2006), arXiv:astroph/0603690 .
 Giannantonio et al. (2006) T. Giannantonio, R. G. Crittenden, R. C. Nichol, R. Scranton, G. T. Richards, A. D. Myers, R. J. Brunner, A. G. Gray, A. J. Connolly, and D. P. Schneider, Phys. Rev. D 74, 063520 (2006), arXiv:astroph/0607572 .
 Raccanelli et al. (2008) A. Raccanelli, A. Bonaldi, M. Negrello, S. Matarrese, G. Tormen, and G. de Zotti, MNRAS 386, 2161 (2008), arXiv:0802.0084 .
 Rassat et al. (2007) A. Rassat, K. Land, O. Lahav, and F. B. Abdalla, MNRAS 377, 1085 (2007), arXiv:astroph/0610911 .
 Xia et al. (2009) J. Xia, M. Viel, C. Baccigalupi, and S. Matarrese, J. Cosmology Astropart. Phys. 9, 3 (2009), arXiv:0907.4753 [astroph.CO] .
 Ferraro et al. (2015) S. Ferraro, B. D. Sherwin, and D. N. Spergel, Phys. Rev. D D91, 083533 (2015), arXiv:1401.1193 [astroph.CO] .
 HernandezMonteagudo et al. (2014) C. HernandezMonteagudo et al., MNRAS 438, 1724 (2014), arXiv:1303.4302 [astroph.CO] .
 Shajib and Wright (2016) A. J. Shajib and E. L. Wright, Astrophys. J. 827, 116 (2016), arXiv:1604.03939 [astroph.CO] .
 Pietrobon et al. (2006) D. Pietrobon, A. Balbi, and D. Marinucci, Phys. Rev. D D74, 043524 (2006), arXiv:astroph/0606475 [astroph] .
 McEwen et al. (2007) J. D. McEwen, P. Vielva, M. P. Hobson, E. MartinezGonzalez, and A. N. Lasenby, MNRAS 376, 1211 (2007), arXiv:astroph/0602398 [astroph] .
 Vielva et al. (2006) P. Vielva, E. MartinezGonzalez, and M. Tucci, MNRAS 365, 891 (2006), arXiv:astroph/0408252 [astroph] .
 Giannantonio et al. (2008) T. Giannantonio, R. Scranton, R. G. Crittenden, R. C. Nichol, S. P. Boughn, A. D. Myers, and G. T. Richards, Phys. Rev. D D77, 123520 (2008), arXiv:0801.4380 [astroph] .
 Giannantonio et al. (2012) T. Giannantonio, R. Crittenden, R. Nichol, and A. J. Ross, MNRAS 426, 2581 (2012), arXiv:1209.2125 [astroph.CO] .
 Ho et al. (2008) S. Ho, C. Hirata, N. Padmanabhan, U. Seljak, and N. Bahcall, Phys. Rev. D 78, 043519 (2008), arXiv:0801.0642 .
 Francis and Peacock (2010) C. L. Francis and J. A. Peacock, MNRAS 406, 2 (2010), arXiv:0909.2494 [astroph.CO] .
 Scranton et al. (2003) R. Scranton, A. J. Connolly, R. C. Nichol, A. Stebbins, I. Szapudi, D. J. Eisenstein, N. Afshordi, T. Budavari, I. Csabai, J. A. Frieman, J. E. Gunn, D. Johnston, Y. Loh, R. H. Lupton, C. J. Miller, E. S. Sheldon, R. S. Sheth, A. S. Szalay, M. Tegmark, and Y. Xu, ArXiv eprints (2003), arXiv:astroph/0307335 .
 Sawangwit et al. (2010) U. Sawangwit, T. Shanks, R. D. Cannon, S. M. Croom, N. P. Ross, and D. A. Wake, MNRAS 402, 2228 (2010), arXiv:0911.1352 [astroph.CO] .
 Audren et al. (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, J. Cosmology Astropart. Phys. 1302, 001 (2013), arXiv:1210.7183 [astroph.CO] .
 Lesgourgues et al. (2008) J. Lesgourgues, W. Valkenburg, and E. Gaztanaga, Phys. Rev. D D77, 063505 (2008), arXiv:0710.5525 [astroph] .
 Nishizawa (2014) A. J. Nishizawa, PTEP 2014, 06B110 (2014), arXiv:1404.5102 [astroph.CO] .
 Limber (1953) D. N. Limber, Astrophys. J. 117, 134 (1953).
 Blas et al. (2011) D. Blas, J. Lesgourgues, and T. Tram, J. Cosmology Astropart. Phys. 1107, 034 (2011), arXiv:1104.2933 [astroph.CO] .
 Beutler et al. (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. StaveleySmith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, MNRAS 416, 3017 (2011), arXiv:1106.3366 [astroph.CO] .
 Ross et al. (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, MNRAS 449, 835 (2015), arXiv:1409.3242 [astroph.CO] .
 Anderson et al. (2014) L. Anderson et al. (BOSS), MNRAS 441, 24 (2014), arXiv:1312.4877 [astroph.CO] .
 Cuoco et al. (2017) A. Cuoco, M. Bilicki, J.Q. Xia, and E. Branchini, ApJS 232, 10 (2017), arXiv:1709.01940 [astroph.HE] .
 Bilicki et al. (2014) M. Bilicki, T. H. Jarrett, J. A. Peacock, M. E. Cluver, and L. Steward, ApJS 210, 9 (2014), arXiv:1311.5246 [astroph.CO] .
 Peacock et al. (2016) J. A. Peacock, N. C. Hambly, M. Bilicki, H. T. MacGillivray, L. Miller, M. A. Read, and S. B. Tritton, MNRAS 462, 2085 (2016), arXiv:1607.01189 .
 Jarrett et al. (2000) T. H. Jarrett, T. Chester, R. Cutri, S. Schneider, M. Skrutskie, and J. P. Huchra, AJ 119, 2498 (2000), astroph/0004318 .
 Wright et al. (2010) E. L. Wright, P. R. M. Eisenhardt, A. K. Mainzer, et al., AJ 140, 18681881 (2010), arXiv:1008.0031 [astroph.IM] .
 Alonso et al. (2015) D. Alonso, A. I. Salvador, F. J. Sánchez, M. Bilicki, J. GarcíaBellido, and E. Sánchez, MNRAS 449, 670 (2015), arXiv:1412.5151 .
 Steward (2014) L. Steward, Master’s thesis, University of Cape Town (2014).
 Bilicki et al. (2016) M. Bilicki, J. A. Peacock, T. H. Jarrett, M. E. Cluver, N. Maddox, M. J. I. Brown, E. N. Taylor, N. C. Hambly, A. Solarz, B. W. Holwerda, I. Baldry, J. Loveday, A. Moffett, A. M. Hopkins, S. P. Driver, M. Alpaslan, and J. BlandHawthorn, ApJS 225, 5 (2016), arXiv:1607.01182 .
 Goto et al. (2012) T. Goto, I. Szapudi, and B. R. Granett, MNRAS 422, L77 (2012), arXiv:1202.5306 [astroph.CO] .
 Kovács et al. (2013) A. Kovács, I. Szapudi, B. R. Granett, and Z. Frei, MNRAS 431, L28 (2013), arXiv:1301.0475 .
 Beck et al. (2016) R. Beck, L. Dobos, T. Budavári, A. S. Szalay, and I. Csabai, MNRAS 460, 1371 (2016), arXiv:1603.09708 .
 Richards et al. (2009) G. T. Richards, A. D. Myers, A. G. Gray, R. N. Riegel, R. C. Nichol, R. J. Brunner, A. S. Szalay, D. P. Schneider, and S. F. Anderson, ApJS 180, 67 (2009), arXiv:0809.3952 .
 Condon et al. (1998) J. J. Condon, W. D. Cotton, E. W. Greisen, Q. F. Yin, R. A. Perley, G. B. Taylor, and J. J. Broderick, AJ 115, 1693 (1998).
 de Zotti et al. (2010) G. de Zotti, M. Massardi, M. Negrello, and J. Wall, A&ARv 18, 1 (2010), arXiv:0908.1896 [astroph.CO] .
 Szapudi et al. (2001) I. Szapudi, S. Prunet, and S. Colombi, ApJ 561, L11 (2001).
 Chon et al. (2004) G. Chon, A. Challinor, S. Prunet, E. Hivon, and I. Szapudi, MNRAS 350, 914 (2004), astroph/0303414 .
 Efstathiou (2004a) G. Efstathiou, MNRAS 348, 885 (2004a), astroph/0310207 .
 Challinor and Chon (2005) A. Challinor and G. Chon, MNRAS 360, 509 (2005), astroph/0410097 .
 Xia et al. (2015) J.Q. Xia, A. Cuoco, E. Branchini, and M. Viel, ApJS 217, 15 (2015), arXiv:1503.05918 .
 Landy and Szalay (1993) S. D. Landy and A. S. Szalay, Astrophys. J. 412, 64 (1993).
 Efstathiou (2004b) G. Efstathiou, MNRAS 349, 603 (2004b), astroph/0307515 .
 Takahashi et al. (2012) R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Astrophys. J. 761, 152 (2012), arXiv:1208.2701 [astroph.CO] .
 Fornasa et al. (2016) M. Fornasa et al., Phys. Rev. D D94, 123005 (2016), arXiv:1608.07289 [astroph.HE] .
 Feroz et al. (2009) F. Feroz, M. P. Hobson, and M. Bridges, MNRAS 398, 1601 (2009), arXiv:0809.3437 [astroph] .
 Sunyaev and Zeldovich (1972) R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4, 173 (1972).
 Rees and Sciama (1968) M. J. Rees and D. W. Sciama, Nature (London) 217, 511 (1968).
 Corasaniti et al. (2005) P.S. Corasaniti, T. Giannantonio, and A. Melchiorri, Phys. Rev. D D71, 123521 (2005), arXiv:astroph/0504115 [astroph] .
 Pogosian et al. (2005) L. Pogosian, P. S. Corasaniti, C. StephanOtto, R. Crittenden, and R. Nichol, Phys. Rev. D D72, 103519 (2005), arXiv:astroph/0506396 [astroph] .
 Pogosian (2005) L. Pogosian, J. Cosmology Astropart. Phys. 0504, 015 (2005), arXiv:astroph/0409059 [astroph] .
 Linder (2003) E. V. Linder, Phys. Rev. Lett. 90, 091301 (2003), arXiv:astroph/0208512 [astroph] .
 Hu and Sawicki (2007) W. Hu and I. Sawicki, Phys. Rev. D 76, 104043 (2007), arXiv:0708.1190 .
 Fang et al. (2008) W. Fang, W. Hu, and A. Lewis, Phys. Rev. D 78, 087303 (2008), arXiv:0808.3125 .
 Ando et al. (2017) S. Ando, A. BenoitLévy, and E. Komatsu, ArXiv eprints (2017), arXiv:1706.05422 .
 Amendola et al. (2016) L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker, M. Baldi, N. Bartolo, A. Blanchard, C. Bonvin, S. Borgani, E. Branchini, C. Burrage, S. Camera, C. Carbone, L. Casarini, M. Cropper, C. de Rham, J. P. Dietrich, C. Di Porto, R. Durrer, A. Ealet, P. G. Ferreira, F. Finelli, J. GarciaBellido, T. Giannantonio, L. Guzzo, A. Heavens, L. Heisenberg, C. Heymans, et al., ArXiv eprints (2016), arXiv:1606.00180 .
 LSST Science Collaboration et al. (2009) LSST Science Collaboration, P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Arnett, S. J. Asztalos, T. S. Axelrod, S. Bailey, D. R. Ballantyne, J. R. Bankert, W. A. Barkhouse, J. D. Barr, L. F. Barrientos, A. J. Barth, J. G. Bartlett, A. C. Becker, J. Becla, and et al., ArXiv eprints (2009), arXiv:0912.0201 [astroph.IM] .
 Jarvis et al. (2015) M. Jarvis, D. Bacon, C. Blake, M. Brown, S. Lindsay, A. Raccanelli, M. Santos, and D. J. Schwarz, Advancing Astrophysics with the Square Kilometre Array (AASKA14) , 18 (2015), arXiv:1501.03825 .
 Renk et al. (2017) J. Renk, M. Zumalacárregui, F. Montanari, and A. Barreira, ArXiv eprints (2017), arXiv:1707.02263 .
 Górski et al. (2005) K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, Astrophys. J. 622, 759 (2005), arXiv:astroph/0409513 .
Appendix A Autocorrelation results
In this appendix we show the measured APS and the related bestfit model for all the catalogs and bins considered in the analysis. Dots refer to the measured single multipoles, while data points with error bars refer to binned measurements.
Appendix B Crosscorrelation results
In this appendix we show the measured CAPS and the related bestfit model for all the catalogs and bins considered in the analysis. Dots refer to the measured single multipoles, while data points with error bars refer to binned measurements.