All-sky angular power spectra from cleaned WISE$\times$SuperCOSMOS galaxy number counts

TL;DR

This paper develops a method to correct systematic effects in galaxy count data from WISE and SuperCOSMOS, improving angular power spectrum measurements and cosmological parameter estimation on large scales.

Contribution

It introduces a new integrated model for stellar contamination and obscuration effects, and demonstrates its effectiveness on mock and real data for cosmological analysis.

Findings

Improved angular power spectrum measurements at large scales.

Consistent cosmological parameters with Planck $ m{ extLambda CDM}$ model.

Identification of residual systematics in southern hemisphere data.

Abstract

Aiming to extract cosmological information from linear scales of the WISE$\times$SuperCOSMOS photometric redshift catalog, we perform a characterization of the systematic effects associated with stellar content, evidencing the presence of contamination and obscuration. We create an integrated model for these effects (which together we call `usurper contamination'), devise a method to remove both of them simultaneously and show its functionality by applying it to a set of mock catalogs. When administered to WISE$\times$SuperCOSMOS data, our method shows to improve the measurements of angular power spectra on scales $\ell\lesssim15$ and the extraction of cosmological parameters therefrom, even though a significant excess of power remains at these scales. When ignoring scales $\ell<15$, we still find strong indications of systematics, albeit these can be localized in the southern equatorial hemisphere. An independent analysis of the northern hemisphere at $\ell\geq 15$ agrees with a $\Lambda$CDM model with parameters from the Planck satellite and gives $\Omega_{\mathrm{c}}=0.254\pm0.020$ and $\Omega_{\mathrm{b}}<0.065$ at 95% confidence limit when combined with priors on $H_0$, $A_s$ and $n_s$.

Figures

Tables

Table 1. Table 1 : Redshift bin number, photo- z 𝑧 z range and estimated values for the mean spectroscopic redshift and obscuration α 𝛼 \alpha , contamination β 𝛽 \beta and true mean galaxy density n ¯ g subscript ¯ 𝑛 g \bar{n}_{\mathrm{g}} parameters. The uncertainties were estimated from mock realizations described in Sec.  5.1 .

Bin	Photo-$z$ range	Mean spec-$z$	$\alpha$	$\beta$	${\overline{n}}_g$
1		0.1280	0.0902(156)	0.421(39)	3.223(48)
2		0.1753	0.0681(75)	0.259(28)	5.499(41)
3		0.2248	0.0735(56)	0.136(27)	6.601(33)
4		0.2719	0.0818(50)	0.095(25)	5.826(23)
5		0.3180	0.0674(80)	0.051(13)	1.748(08)

Table 2. Table 2 : Maximum multipole (i.e. linear scales) ℓ max subscript ℓ max \ell_{\mathrm{max}} included in our MCMC estimation of the posterior, and the multipole at which the shot noise starts to dominate, ℓ shot subscript ℓ shot \ell_{\mathrm{shot}} , for each C ℓ i ​ j subscript superscript 𝐶 𝑖 𝑗 ℓ C^{ij}_{\ell} .

Bin $i$	Bin $j$	${\ell }_{\max }$	${\ell }_{\mathrm{shot}}$
2	2	70	362
2	3	80	–
2	4	89	–
2	5	103	–
3	3	93	367
3	4	104	–
3	5	117	–
4	4	117	316
4	5	129	–
5	5	140	180

Table 3. Table 3 : Comparison of true obscuration, contamination and mean galaxy density parameters and those recovered from 2500 simulations. The columns are: the redshift bin and the true value x 0 subscript 𝑥 0 x_{0} , the mean value from the simulations x ¯ ¯ 𝑥 \bar{x} , the standard deviation of the simulations σ x subscript 𝜎 𝑥 \sigma_{x} , and the fractional bias ( x ¯ / x 0 − 1 ¯ 𝑥 subscript 𝑥 0 1 \bar{x}/x_{0}-1 ) in percent for each parameter x 𝑥 x .

Bin	${\alpha }_0$	$\overline{\alpha }$	${\sigma }_{\alpha }$	% bias	${\beta }_0$	$\overline{\beta }$	${\sigma }_{\beta }$	% bias	$n_{\mathrm{g0}}$	${\overline{n}}_g$	${\sigma }_{n_g}$	% bias
1	0.0902	0.0934	0.0157	3.4	0.421	0.435	0.039	3.2	3.223	3.22	0.049	-0.1
2	0.0681	0.0705	0.0077	3.5	0.259	0.273	0.029	5.3	5.499	5.498	0.042	-0.0
3	0.0735	0.0755	0.0056	2.6	0.136	0.15	0.028	9.0	6.601	6.601	0.033	-0.0
4	0.0818	0.0836	0.005	2.2	0.095	0.106	0.024	10.5	5.826	5.826	0.023	0.0
5	0.0674	0.0703	0.0081	4.0	0.051	0.056	0.013	8.8	1.748	1.748	0.008	0.0

Table 4. Table 4 : Summary statistics showing the performance of different processing strategies and best-fit models used to fit the angular power spectra of the WSC data. The first column names the strategy adopted. The second displays the best-fit χ 2 superscript 𝜒 2 \chi^{2} , computed with respect to the C ℓ subscript 𝐶 ℓ C_{\ell} s that best-fit that particular processing strategy. The third column contains the number of degrees of freedom, while the last column displays the p 𝑝 p -value. All the approaches shown in the table that are restricted to ℓ ≥ 15 ℓ 15 \ell\geq 15 (or ℓ ≥ 21 ℓ 21 \ell\geq 21 ) were also subjected to the cleaning process.

Approach	${\chi }^2$	d.o.f.	$p$-value
Uncleaned	268	156	$5.87\times {10}^{-8}$
Cleaned	237	156	$2.92\times {10}^{-5}$
$\ell \ge 15$	163	136	$0.058$
$\ell \ge 21$	151	126	$0.063$
North uncleaned	398	156
North cleaned	376	156
North $\ell \ge 15$	148	136	$0.23$
South uncleaned	248	156	$3.77\times {10}^{-6}$
South cleaned	185	156	$0.054$
South $\ell \ge 15$	157	136	$0.10$

Table 5. Table 5 : Best-fit values of the parameters of a Λ Λ \Lambda CDM model and the mean and standard deviation of their marginal posteriors for the full-sky, northern and southern datasets. All datasets shown here were cleaned and fit from ℓ ≥ 15 ℓ 15 \ell\geq 15 to linear scales. Since the posterior of Ω b subscript Ω b \Omega_{\mathrm{b}} reaches the CLASS limit Ω b = 0.0065 subscript Ω b 0.0065 \Omega_{\mathrm{b}}=0.0065 , we quote its 95% confidence upper limit instead of the mean.

	Full-sky $\ell \ge 15$		North $\ell \ge 15$		South $\ell \ge 15$
Parameter	Best-fit	Mean $\pm \sigma$	Best-fit	Mean $\pm \sigma$	Best-fit	Mean $\pm \sigma$
$h$	$0.6703$	$0.6729\left(99\right)$	$0.6758$	$0.6729\left(99\right)$	$0.6758$	$0.6723\left(97\right)$
${\Omega }_b$	$0.0321$		$0.0340$		$0.0076$
${\Omega }_c$	$0.214$	$0.211\left(14\right)$	$0.251$	$0.254\left(20\right)$	$0.178$	$0.184\left(14\right)$
$\ln \left({10}^{10}{A}_s\right)$	$3.103$	$3.088\left(36\right)$	$3.092$	$3.086\left(36\right)$	$3.078$	$3.088\left(36\right)$
$n_s$	$0.9678$	$0.9652\left(64\right)$	$0.9687$	$0.9657\left(63\right)$	$0.9638$	$0.9653\left(65\right)$
$b_2$	$1.432$	$1.424\left(75\right)$	$1.274$	$1.322\left(80\right)$	$1.394$	$1.455\left(82\right)$
$b_3$	$1.284$	$1.274\left(63\right)$	$1.195$	$1.229\left(68\right)$	$1.251$	$1.314\left(67\right)$
$b_4$	$1.149$	$1.144\left(55\right)$	$1.107$	$1.134\left(61\right)$	$1.104$	$1.152\left(58\right)$
$b_5$	$1.264$	$1.266\left(64\right)$	$1.254$	$1.274\left(75\right)$	$1.210$	$1.238\left(70\right)$
${\sigma }_{z,2}$	$0.0575$	$0.0581\left(22\right)$	$0.0587$	$0.0598\left(36\right)$	$0.05591$	$0.0557\left(29\right)$
${\sigma }_{z,3}$	$0.04970$	$0.04969\left(90\right)$	$0.0495$	$0.0498\left(15\right)$	$0.0496$	$0.0500\left(12\right)$
${\sigma }_{z,4}$	$0.04401$	$0.04400\left(73\right)$	$0.0449$	$0.0453\left(11\right)$	$0.0436$	$0.0439\left(10\right)$
${\sigma }_{z,5}$	$0.0353$	$0.0354\left(13\right)$	$0.0384$	$0.0379\left(18\right)$	$0.0344$	$0.0333\left(17\right)$
${\sigma }_8$	$0.7915$	$0.796\left(35\right)$	$0.869$	$0.854\left(43\right)$	$0.792$	$0.763\left(38\right)$

Equations

$$p(z_{\mathrm{s}}|z_{\mathrm{p}})=\frac{\sqrt{2\pi a(z_{\mathrm{p}})}\Gamma \left[a(z_{\mathrm{p}})-\frac{1}{2}\right]}{\Gamma \left[a(z_{\mathrm{p}})\right]}{\left[1+\frac{(z_{\mathrm{s}}-z_{\mathrm{p}}{)}^2}{2a(z_{\mathrm{p}})s^2(z_{\mathrm{p}})}\right]}^{-a(z_{\mathrm{p}})},$$

$$\frac{\mathrm{d}{N}^i}{\mathrm{d}{z}_{\mathrm{s}}}={\int }_{z_{\min }^i}^{z_{\max }^i}p(z_{\mathrm{s}}|z_{\mathrm{p}})\frac{\mathrm{d}N}{\mathrm{d}{z}_{\mathrm{p}}}\mathrm{d}{z}_{\mathrm{p}},$$

$$n_{\mathrm{stars}}(b)=\bar{n}{e}^{-a|b|}+c,$$

$$n_{\mathrm{obs}}^i(\mathbf{\theta })=W(\mathbf{\theta })\{[1-{\alpha }_{i}S(\mathbf{\theta })]n_{\mathrm{g}}^i(\mathbf{\theta })+{\beta }_{i}S(\mathbf{\theta })+{ϵ}_{\mathrm{g},i}(\mathbf{\theta })+{ϵ}_{\mathrm{s},i}(\mathbf{\theta })\}.$$

$${\bar{n}}_{\mathrm{obs}}^i(s)=\frac{1}{A(s)}\int n_{\mathrm{obs}}^i(s,\lambda )\mathrm{d}\lambda =[1-{\alpha }_{i}S(s)]{\bar{n}}_{\mathrm{g},i}(s)+{\beta }_{i}S(s),$$

$$A(s)\equiv \int W(s,\lambda )\mathrm{d}\lambda ,$$

$${\bar{n}}_{\mathrm{g},i}(s)\equiv \frac{1}{A(s)}\int W(s,\lambda )n_{\mathrm{g}}^i(s,\lambda )\mathrm{d}\lambda .$$

$${\sigma }_{\mathrm{obs},i}^2(s)=\frac{1}{A(s)}\int {\delta }_{\mathrm{obs},i}^2(s,\lambda )\mathrm{d}\lambda \simeq [1-{\alpha }_{i}S(s){]}^2{v}_{\mathrm{g}}^i+{\bar{n}}_{\mathrm{obs}}^i(s),$$

$${\sigma }_{\mathrm{obs},i}^2(s)-{\bar{n}}_{\mathrm{obs}}^i(s)=[1-{\alpha }_{i}S(s){]}^2{v}_{\mathrm{g}}^i.$$

$$n_{\mathrm{w}}^i(\mathbf{\theta })\equiv \frac{n_{\mathrm{obs}}^i(\mathbf{\theta })}{1-{\hat{\alpha }}_{i}S(\mathbf{\theta })}\simeq W(\mathbf{\theta })\left[n_{\mathrm{g}}^i(\mathbf{\theta })+{\beta }_i\frac{S(\mathbf{\theta })}{1-{\hat{\alpha }}_{i}S(\mathbf{\theta })}+\frac{{ϵ}_{\mathrm{g},i}(\mathbf{\theta })+{ϵ}_{\mathrm{s},i}(\mathbf{\theta })}{1-{\hat{\alpha }}_{i}S(\mathbf{\theta })}\right].$$

$${\bar{S}}_{\mathrm{w}}^i\equiv \frac{1}{A}\int \frac{S(\mathbf{\theta })}{1-{\hat{\alpha }}_{i}S(\mathbf{\theta })}{\mathrm{d}}^2\theta ,$$

$${\beta }_i\simeq \frac{\int \delta n_{\mathrm{w}}^i(\mathbf{\theta })\delta S_{\mathrm{w}}^i(\mathbf{\theta }){\mathrm{d}}^2\theta }{\int [\delta S_{\mathrm{w}}^i(\mathbf{\theta }){]}^2{\mathrm{d}}^2\theta },$$

$${\hat{n}}_{\mathrm{g},i}=\frac{\int [n_{\mathrm{obs}}^i(\mathbf{\theta })-{\hat{\beta }}_{i}S(\mathbf{\theta })]{\mathrm{d}}^2\theta }{\int W(\mathbf{\theta })[1-{\hat{\alpha }}_{i}S(\mathbf{\theta })]{\mathrm{d}}^2\theta }.$$

$${\hat{\delta }}_{\mathrm{g}}^i(\mathbf{\theta })=W(\mathbf{\theta })\left\{\frac{n_{\mathrm{obs}}^i(\mathbf{\theta })-{\hat{\beta }}_{i}S(\mathbf{\theta })}{{\hat{n}}_{\mathrm{g},i}[1-{\hat{\alpha }}_{i}S(\mathbf{\theta })]}-1\right\}.$$

$$D_{\ell }^{ij}\equiv \frac{1}{2\ell +1}\sum _{m=-\ell }^{\ell }{a}_{\ell m}^i{a}_{\ell m}^{j\ast }.$$

$$D_{\ell }^{ij}=\sum _{{\ell }^{\prime }}{R}_{\ell {\ell }^{\prime }}{F}_{{\ell }^{\prime }}^{ij}.$$

$$C_{\ell }^{ij}=\frac{F_{\ell }^{ij}-N^{ij}}{w_{\ell }^2},$$

$$N^{ij}=\frac{{\delta }_{\mathrm{K}}^{ij}}{{\bar{n}}_{\mathrm{g},i}^2}\int \frac{n_{\mathrm{obs}}^i(\mathbf{\theta })}{[1-{\hat{\alpha }}_{i}S(\mathbf{\theta }){]}^2}{\mathrm{d}}^2\theta$$

Full text

1 Introduction

Galaxy and quasar surveys are key observational tools in cosmology. Through them we can probe the matter distribution in the Universe and, therefore, test predictions of cosmological models and constrain their parameters [1, 2, 3, 4, 5, 6, 7, 8, 9]. A few interesting questions that have been tackled using wide-angle surveys are the homogeneity and isotropy of the Universe [10, 11, 12], the shape of the distribution function of galaxy density [13], the behavior of gravity and the growth of structure [14, 15], the accelerated expansion of the Universe and the properties of dark energy [16, 17], the amount of matter in the Universe [18], the existence of primordial non-Gaussianities [19] and the mass of neutrinos [20, 21].

To improve constraints on the relevant cosmological parameters – especially on the equation of state of dark energy, the mass of neutrinos, and on non-Gaussianities and other inflation signatures – deeper and wider surveys are desirable [22, 23, 24]. Moreover, these surveys are required in order to improve our knowledge in various other interesting topics. For instance, the cross-correlation between galaxy distribution and gravitational lensing of the cosmic microwave background (CMB) – which serves as a consistency check for our cosmological model, besides helping improve constraints on halo bias and on the amplitude of matter perturbations – peaks at large angular scales [25]. This is also the case for the integrated Sachs-Wolfe effect [26], probed by the cross-correlation of CMB temperature maps with galaxy distribution, whose signal comes from the evolution of gravitational potential along the CMB photon path. Another topic that requires almost full-sky galaxy surveys is the measurement of our own peculiar velocity with respect to all other galaxies [27, 28, 29]. Large galaxy surveys can also help us improve our understanding of galaxy bias, that could be enhanced at the largest scales if we assume that galaxies form at density peaks (see e.g. [30]).

This type of studies are among the main aims of the largest forthcoming surveys. The Large Synoptic Survey Telescope222http://www.lsst.org/ (LSST) [31] will cover $18,000$ $\mathrm{deg^{2}}$ up to redshifts $z\lesssim 3$, while the Euclid space telescope333http://www.euclid-ec.org will observe 15,000 $\mathrm{deg^{2}}$ and detect galaxies up to $z\lesssim 2$ [32]. The huge volume probed by these surveys will lead to such a high number of observed galaxies and quasars that spectroscopic measurements will be performed only for a small fraction of them, while for most sources we will have to rely on photometric redshifts (photo-$z$s hereafter) and photometric classification. These techniques are, however, plagued with high contamination levels, large systematic effects and uncertainties, which must be controlled and understood in order to fully mine the relevant datasets.

An existing testbed for such surveys is the WISE$\times$SuperCOSMOS catalog [33] (WSC hereafter), a full-sky photometric galaxy dataset reaching $z\lesssim 0.4$. Our goal in this paper is to estimate cosmological parameters from this sample under the standard cosmological model $\Lambda$CDM. In this process, we characterize the systematic effects associated with stellar density (which, although more prominent in the WSC, will also be present in LSST and Euclid) and devise a method for mitigating such effects, tested on simulations and real data. These efforts were aimed at improving the survey’s reach toward the largest scales ($\ell\lesssim 20$), commonly distrusted in cosmological analyses (e.g. [34, 35]). Finally, we also verified the impact of different assumptions about photo-$z$ properties on the measured parameters.

This work adds to previous WSC analyses that have investigated its observational and cosmological properties, such as hemispherical anisotropies [36], Minkowski functionals [37], cross-correlation with CMB lensing [38, 39] and temperature [40]. It also contributes with a new strategy for mitigating systematics to those previously proposed [41, 42, 43, 44, 45, 46].

The basic technique we employ here – the measurement of angular power spectra in tomographic redshift shells – has been applied before to other photometric galaxy catalogs such as the MegaZ [47] and the 2MASS Photometric Redshift datasets (2MPZ) [34], to photometric quasar samples [42] and spectroscopic galaxy data [35].

The outline of this paper is as follows: in Sec. 2 we describe the datasets used in this work, namely the WSC, the Sloan Digital Sky Survey, and the Gaia Data Release 2 catalogs. In Secs. 3.1 and 3.2 we expose the stellar contamination and obscuration effects existent in WSC. A method of modeling and mitigating these effects is presented in Sec. 3.3. The WSC data had their angular power spectra measured and cosmological parameters inferred using the methodology described in Sec. 4. Using simulations, we validate our treatment of systematic effects and our measurement methodology in Sec. 5. Our results – WSC angular power spectra and inferred cosmological parameters – are shown in Sec. 6. Sec. 7 presents tests for robustness and systematics, and in particular, Sec. 7.7 shows evidence for a systematic effect represented by a tension in the cosmological parameters constrained from the angular clustering measured from different equatorial hemispheres. Finally, we conclude and summarize our work in Sec. 8.

2 The dataset

The main dataset used in this work is the WISE$\times$SuperCOSMOS (WSC) galaxy catalog [33], described in Sec. 2.1. To estimate the stellar contamination and obscuration in WSC, we used the Sloan Digital Sky Survey (SDSS) [48] and Gaia DR2 [49, 50] photometric catalogs as tracers of stars. These are described in Secs. 2.2 and 2.3, respectively.

2.1 WISE$\times$SuperCOSMOS

The WSC catalog is a cross-match between two parent full-sky catalogs: the AllWISE release [51] from the Wide-field Infrared Survey Explorer (WISE) [52], a mid-infrared space survey in four bands $W1$–$W4$ (3.4, 4.6, 12 and 22$\mathrm{\mu m}$); and the SuperCOSMOS Sky Survey [53], a program of automated digitalization of optical photographic plates in the $B$, $R$ and $I$ filters, taken by the United Kingdom Schmidt Telescope (UKST, in the southern hemisphere) and the Palomar Observatory Sky Survey-II (POSS-II, in the northern hemisphere).

To build the WSC catalog, only the bands $W1$, $W2$, $B$ and $R$ were used. Following the magnitude limits of the parent surveys, all WSC sources have magnitudes $W1<17$, $B<21$ and $R<19.5$. To preferentially select galaxies over quasars and stars, WSC only includes extended sources according to SuperCOSMOS morphological classification meanClass$=1$. It should be noted that the SuperCOSMOS morphological classification is not as accurate as most recent optical surveys; its angular resolution is approximately $2^{\prime\prime}$ [54], and the quality of the imaging technique is lower (photographic plates vs. CCDs).

We also point out that the WSC star/galaxy separation is not based only on SuperCOSMOS morphology but also on WISE colours [33]. WSC sources have $W1>13.8$ (a cut aimed at removing bright stars), $R-W2<7.6-4(W1-W2)$ and $W1-W2<0.9$ (two color-cuts aimed at removing quasars). In this work, we further enforced: the requirement $W1-W2>0.2$ for all sources to reduce stellar contamination and maintain a constant galaxy selection function across the sky (in opposition to a cut dependent on Galactic latitude [33]); maximum Galactic extinction, $E(B-V)<0.10$; removal of sources in highly contaminated regions such as near the Galactic plane and bulge, around the Magellanic Clouds and in regions presenting unusually high densities according to a lognormal distribution, all accomplished with the WSC final mask [33]; removal of sources in regions where Gaia stellar density is greater than 7 times its average (this increases the masked region around the bulge; see Sec. 2.3); and removal of a SuperCOSMOS tile with bad photometry.

The photo-$z$s for all the WSC sources were obtained with the artificial neural network code ANN$z$ [56]444http://www.homepages.ucl.ac.uk/~ucapola/annz.html, trained on a WSC cross-match with the complete and deep Galaxy And Mass Assembly (GAMA)-II spectroscopic dataset [57]. The adequacy of GAMA for photo-$z$ training in WSC has been discussed and demonstrated in [33]. As explained there, these photo-$z$s were further corrected for asymmetries between their northern and southern distributions, likely caused by differences between the POSS-II and UKST pass-bands. That paper contains also a comprehensive analysis of photo-$z$ properties based on comparisons with several external spectroscopic datasets. A projection of all sources in the photometric redshift (photo-$z$) range $0.10<z<0.35$ is presented in Fig. 1, and Fig. 2 shows the photo-$z$ distribution of WSC galaxies. In this work we binned the WSC sources in photo-$z$ shells of width $\Delta z=0.05$, covering the range $0.10<z<0.35$. However, for reasons presented in Sec. 3, most of our analysis ignores the first bin ($0.10<z<0.15$).

It was shown in [38] that the WSC photo-$z$ errors can be modeled by a Gaussian with zero mean (i.e. no bias) and standard deviation $\sigma(z_{\mathrm{p}})=0.02+0.08z_{\mathrm{p}}$. Also, that an even better model would be a generalized Lorentzian, where the probability density of a galaxy with photo-$z$ $z_{\mathrm{p}}$ to have true redshift $z_{\mathrm{s}}$ is:

[TABLE]

where $\Gamma(x)$ denotes the Gamma function. The WSC photo-$z$s follows $a(z)=3-4z$ and $s(z)=0.02+0.04z$.

The spectroscopic redshift (spec-$z$) distribution $\mathrm{d}N^{i}/\mathrm{d}z_{\mathrm{s}}$ of the galaxies in our photo-$z$ bin $i$ can be estimated from a convolution of $p(z_{\mathrm{s}}|z_{\mathrm{p}})$ and the photo-$z$ distribution $\mathrm{d}N/\mathrm{d}z_{\mathrm{p}}$:

[TABLE]

where the $i$th bin in photo-$z$ ranges from $z_{\mathrm{min}}^{i}$ to $z_{\mathrm{max}}^{i}$. The resulting distributions in the different redshift bins are shown in Fig. 2. We verified that the expected spec-$z$ distributions in each bin are very well fitted by a Gaussian (generalized Lorentzian) when the photo-$z$ errors are Gaussian (generalized Lorentzian). Also, one can see in Fig. 2 that the mean values of the spec-$z$s in every bin are practically independent of the distribution assumed for the photo-$z$ errors, whereas the FWHM may depend on it. Thus, when computing theoretical angular power spectra $C_{\ell}$s to fit the data, we decided to keep the means of the spec-$z$ distributions fixed and let the widths vary as nuisance parameters.

2.2 Sloan Digital Sky Survey

The Sloan Digital Sky Survey (SDSS, [48]) is one of the largest astronomical surveys in operation. Here we use its $14^{\mathrm{th}}$ Data Release (DR14) [58]. The survey has already covered $\sim 1/3$ of the northern hemisphere, imaging in the five $ugriz$ broad bands. The effective photometric limit of the survey is $r=22.2$ ($95\%$ completeness for point sources). In addition, its spectroscopic counterpart is complete down to $r=17.77$ for the Main Galaxy Sample (MGS) and deeper for other sub-projects, such as the Baryon Oscillation Spectroscopic Survey (BOSS) [16].

Despite containing a wealth of cosmological information, the SDSS data was only used here to verify the existence of stellar contamination and obscuration in the WSC catalog, as described in Sec. 3. The SDSS star-galaxy separation classifies extended objects as galaxies and those point-like as stars. This classification is based on a difference between two types of magnitudes, namely sources with $\mathrm{psfMag-cModelMag}>0.145$ are classified as galaxies and otherwise as stars.555The two types of magnitudes $\mathrm{psfMag}$ and $\mathrm{cModelMag}$ represent magnitudes obtained using the point spread function (PSF) model and the best linear combination of exponential and de Vaucouleurs profiles, respectively. Note that the ‘stars’ according to this definition will include any point sources, also extragalactic ones such as quasars. The PSF of $1.24$ arcsec (full width at half maximum, FWHM) allows for a reliable star-galaxy separation up to $r<21.5$ [59] and accurate astrometric positions ($<0.1$ arcsec per coordinate) at $r<20.5$ [60].

2.3 Gaia

Gaia is a space mission of the European Space Agency (ESA) aimed at measuring the three-dimensional positions and velocities of about $10^{9}$ Milky Way stars [49]. It uses two optical telescopes of $0.7\mathrm{m^{2}}$ collecting area equipped with astrometry, photometry and spectrometry instruments. The astrometry is performed in the range $330-1050$ nm ($G$-band hereafter). Gaia Data Release 2 (DR2, [50]) contains $1.3\times 10^{9}$ point sources with measured parallaxes and proper motions at $\sim 1$ milliarcsecond (mas) and $\sim 1$ mas per year ($\mathrm{mas~{}y^{-1}}$) precision, respectively, up to a $G$ band magnitude of 20.

In order to create a representative full-sky stellar density map using the Gaia DR2 source catalog, we had to ensure minimal observational effects (i.e. anisotropic completeness) and minimal contamination by galaxies and quasars. To achieve the first goal we applied a magnitude cut $G<20$; as the limiting magnitudes in the regions of interest (away from the Galactic plane) are $G\gtrsim 20.8$ [61], we expect uniform sky completeness under this selection criterion. The second goal was achieved by taking advantage of the excellent astrometric measurements of Gaia: we assumed that the majority of observed stars should have either a non-zero parallax or a non-zero total proper motion. Thus, we enforced, for all the selected sources, a 5$\sigma$ detection of at least one of these quantities. We computed the total proper motion as $\mu=\sqrt{\mu_{\delta}^{2}+\mu_{\alpha^{\prime}}^{2}}$, where $\mu_{\alpha^{\prime}}\equiv\mu_{\alpha}\cos\delta$, $\mu_{\alpha}$ is the proper motion in the direction of right ascension $\alpha$ and $\mu_{\delta}$ is the proper motion in the direction of declination $\delta$. The uncertainty on $\mu$ was computed with the usual error propagation, taking into account the covariance between $\mu_{\delta}$ and $\mu_{\alpha^{\prime}}$ measurement errors.

In order to avoid regions of high stellar density, here we only selected sources at Galactic latitudes $|b|>10^{\circ}$. A map of the $282,045,470$ selected Gaia sources, in Galactic coordinates, is presented in Fig. 3.

3 Characterizing and cleaning the WSC catalog

3.1 Stellar contamination

As described in [33], the WSC catalog is affected by considerable stellar contamination, especially near the Galactic plane. In order to demonstrate this contamination, we follow the strategy applied in [33] of cross-matching WSC with SDSS photometric sources, taking advantage of the latter’s better star-galaxy separation. Our final analysis, however, is based on a stellar density template derived from Gaia’s observations (see Sec. 3.3).

We cross-matched the WSC sample and the photometric SDSS DR14, linking the nearest objects in the latter (within 2 arcsec) to objects in the former. In total, our cross-matched sample consists of $4,282,564$ sources which includes both stars and galaxies. The result of this cross-match is shown in Fig. 4. From the bottom panel, it is quite clear that stellar contamination exists and increases towards the Galactic plane and bulge. Also, there are hints that the WSC galaxies are partly obscured by stars, since the top panel shows a deficit of galaxies in the same regions where stars creep in. Overall, 96% percent of the matched sources are galaxies and 4% are stars, according to SDSS.

Fig. 5 shows how the amount of stellar contamination changes with redshift. We can see that, at high redshifts, the contamination (i.e. stellar fraction) increases significantly. This increase (and the one at low redshift) are mostly caused by a drop in the number of observed galaxies and not by an increase in the number of misclassified stars (see the top panel).

The average dependence with Galactic latitude $b$ of the number of stars contaminating WSC (according to SDSS classification) is shown in Fig. 6. We see that both Galactic hemispheres, as well as different redshift ranges, follow very similar trends. For instance, the binned redshift curves differ from the full range one mostly by a constant factor. Despite these similarities, subtle slope changes and hemispherical asymmetries are noticeable.

Figure 6 also shows the result of an exponential fit to the average number of stars per pixel $n_{\mathrm{stars}}$ in the range $0.10<z<0.35$ as a function of Galactic latitude $b$:

[TABLE]

where $\bar{n}=7.14$, $a=0.06$ and $c=0.35$; in the plot, the fit was scaled by 0.71 to ease visualization. We can see that the fit represents the stellar density reasonably well.

3.2 Stellar obscuration

As shown in Fig. 4 and already noted in previous works with SDSS data [41], the density of observed galaxies anti-correlates with the stellar density: regions with high concentration of stars tend to present lower galaxy counts. We call this phenomenon stellar obscuration. This is expected to affect mostly photometric datasets, such as WSC, due to issues such as blending and source misidentification. Spectroscopic surveys can avoid stellar contamination given that spectroscopy can easily separate stars from galaxies, but even they suffer from such obscuration, as shown in [41].

The approach used in [41] to estimate stellar obscuration was to correlate the observed number of galaxies with the density of foreground stars: as galaxies are, on average, homogeneously distributed on the sky, these two should be uncorrelated in the absence of systematic effects. This approach is feasible in the spectroscopic SDSS case because the galaxy sample there is free from stellar contamination. For WSC sources that include a large amount of stars disguised as galaxies, however, this approach is not possible: correlating WSC sources with stellar density would mix contributions from obscuration and contamination.

An alternative approach is to look at the standard deviation of the counts of galaxies inside pixels as a function of stellar density. The basic idea is that stellar obscuration will not only modulate the mean observed density $\bar{n}_{\mathrm{obs}}$ but also the cosmological fluctuations around the mean, since the former multiplies the latter. Let us formalize this statement through a model for the observed projected density of sources $n^{i}_{\mathrm{obs}}(\bm{\theta})$ at angular position $\bm{\theta}$, for a particular redshift bin $i$, that takes into account both contamination and obscuration:

[TABLE]

In the equation above, $W(\bm{\theta})$ is the survey binary window function, $S(\bm{\theta})$ is a template for the expected density of stars, $n^{i}_{\mathrm{g}}(\bm{\theta})$ is the real galaxy density, and $\epsilon_{\mathrm{g},i}(\bm{\theta})$ and $\epsilon_{\mathrm{s},i}(\bm{\theta})$ are Poisson fluctuations in the number of observed galaxies and stars, respectively. The non-negative parameters $\alpha_{i}$ and $\beta_{i}$ control the amount of obscuration and contamination in the $i$th bin, respectively. The noise terms obey the statistical properties $\langle\epsilon_{\mathrm{g},i}(\bm{\theta})\rangle=\langle\epsilon_{\mathrm{s},i}(\bm{\theta})\rangle=0$, $\langle\epsilon_{\mathrm{g},i}(\bm{\theta})\epsilon_{\mathrm{s},i}(\bm{\theta})\rangle=0$, $\langle\epsilon^{2}_{\mathrm{g},i}(\bm{\theta})\rangle=W(\bm{\theta})[1-\alpha_{i}S(\bm{\theta})]\bar{n}_{\mathrm{g},i}$ and $\langle\epsilon_{\mathrm{s},i}^{2}(\bm{\theta})\rangle=W(\bm{\theta})\beta_{i}S(\bm{\theta})$, such that the variance of the total noise is equal to the expected number of observed sources, $\langle[\epsilon_{\mathrm{g},i}(\bm{\theta})+\epsilon_{\mathrm{s},i}(\bm{\theta})]^{2}\rangle=\langle n^{i}_{\mathrm{obs}}(\bm{\theta})\rangle$.

Our model for the observed sources (Eq. 3.2) makes three important simplifications: first, the number of observed galaxies is linearly suppressed by the stellar density; second, both contamination and obscuration depend on the same template; third, the template $S(\bm{\theta})$ is independent of redshift, although the final obscuration and contamination are modulated by bin-dependent parameters $\alpha_{i}$ and $\beta_{i}$. The fact that both obscuration and contamination follow the same template leads to a replacement of galaxies for stars, and thus we call this combined effect “usurper contamination”.

To estimate $\alpha_{i}$ from the observed variance of $n^{i}_{\mathrm{obs}}$, we first write $\bm{\theta}$ as a pair of generalized coordinates $\bm{\theta}=(s,\lambda)$, where $s$ specifies an isocontour of $S(\bm{\theta})$ while by varying $\lambda$ we move along this isocontour. We then compute the mean of $n^{i}_{\mathrm{obs}}$ for a fixed $s$:

[TABLE]

where

[TABLE]

and

[TABLE]

Note that due to our choice of coordinates, the value of $S(\bm{\theta})$ only depends on $s$. Also, we assumed that when averaging the noise $\epsilon_{\mathrm{g},i}+\epsilon_{\mathrm{s},i}$ along $\lambda$, we may approximate the result to zero. Now we can define the difference $\delta n_{\mathrm{obs},i}(s,\lambda)\equiv n^{i}_{\mathrm{obs}}(s,\lambda)-W(s,\lambda)\bar{n}^{i}_{\mathrm{obs}}(s)$ and compute the variance of source number counts for an isocontour $s$ of $S(\bm{\theta})$:

[TABLE]

where $v^{i}_{\mathrm{g}}=\langle[n^{i}_{\mathrm{g}}(s,\lambda)-\bar{n}_{\mathrm{g},i}(s)]^{2}\rangle$ is the true galaxy variance (i.e. it does not include systematic effects nor shot-noise). Eq. (3.6) assumes that the mean source number counts in $s$, $\bar{n}_{\mathrm{obs}}^{i}(s)$, approximates the expected value $\langle n^{i}_{\mathrm{obs}}(s,\lambda)\rangle$, and that $v^{i}_{\mathrm{g}}$ is independent of the angular position. By placing all directly measurable quantities on the left hand side, we get a function of $s$ that might be fitted to extract the obscuration parameter $\alpha_{i}$:

[TABLE]

In the absence of obscuration, the quantity on the left would be the variance due to cosmological fluctuations (and thus we call it $\sigma_{\mathrm{cosmo}}^{2}$) and it should be constant all over the sky. We point out that, in practice, the estimation of $\sigma_{\mathrm{cosmo}}^{2}$ is performed in a discrete way [the integral in Eq. (3.6), for instance, is actually a sum over pixels], and the isocontours are approximated by bands of similar $S(\bm{\theta})$ values.

Fig. 7 shows $\sigma_{\mathrm{cosmo}}$ as a function of Galactic latitude $b$ for each photo-$z$ bin. It is completely evident that the amplitude of the observed fluctuations is not constant and decreases towards the Galactic plane, quantitatively proving the existence of obscuration in the WSC data. This figure also shows that our obscuration model, given by the square root of Eq. (3.7) (depicted as a red line) describes really well the overall behavior of this effect (except for the first bin). Here it is important to emphasize that only $\alpha_{i}$ and $v^{i}_{\mathrm{g}}$ were fitted to $\sigma_{\mathrm{cosmo}}$, while $S(s)$ was held fixed according to Eq. (3.1), estimated from contamination. Thus, Fig. 7 also shows that our assumption that both obscuration and contamination follow the same template $S(s)$ is reasonable.

The error bars in Fig. 7 were computed by error propagation from the uncertainties on $\bar{n}_{\mathrm{obs}}^{i}(s)$ and $\sigma_{\mathrm{obs},i}(s)$, which are given by $\sigma_{\mathrm{obs},i}(s)/\sqrt{N(s)}$ and $\sigma_{\mathrm{obs},i}(s)/\sqrt{2[N(s)-1]}$, respectively [$N(s)$ is the number of pixels taken into account in each isocontour]. We interpret the fact that the scatter of the data points is larger than the one expected from the error bars as an indication that $v^{i}_{\mathrm{g}}$ actually varies from one isocontour to another due to large-scale density fluctuations that are the size of the isocontours. It might also indicate that other factors are in play (e.g. Galactic extinction, seeing and bad photographic plates) or that the template $S(\bm{\theta})$ can be improved.

3.3 Cleaning the WSC projected galaxy distribution

We can note from Fig. 6 that the stellar contamination is not perfectly symmetric across the Galactic plane. Moreover, Fig. 3 shows that the stellar density (and likely contamination) is larger near the Galactic bulge, as expected. Aiming at improving the removal of contamination and obscuration from WSC projected galaxy distribution, we built a stellar density template $S(\bm{\theta})$ from the Gaia DR2 map presented in Fig. 3, which should accurately represent the stellar distribution – and thus contamination and obscuration – all over the sky.

A few unusually bright pixels – with excessive object counts due to the presence of globular clusters and nearby dwarf galaxies – can be seen in Fig. 3. To remove these count peaks, we computed, for all the pixels, the median of the number counts in neighboring pixels; if the central pixel value was more than twice the median, we replaced the former by the latter. Then, to attenuate the Poisson noise (visible in Fig. 6), we applied Gaussian smoothing with standard deviation $\sigma_{\mathrm{G}}=34.4$ arcmin in order to retain all the structure in our template up to the scales of interest. Finally, we applied the same mask as in the data (see Sec. 2.1) and normalized the template by its mean value (the template’s absolute scale does not affect the data treatment). The final result is shown in Fig. 8.

Using the template from Fig. 8, we estimated $\alpha_{i}$ for all redshift bins with Eq. (3.7). Fig. 9 shows the data variance in template isocontours as a function of the template’s value and our fit. Again, our model describes reasonably well the obscuration effect, apart from the first bin.

We next removed the contamination using a slightly modified version of the approach proposed in [62, 46], i.e. we computed the correlation of the (weighted) data with the contamination template. The weighted number counts are defined as:

[TABLE]

In the equation above, $\hat{\alpha}_{i}$ is our estimate for the obscuration parameter, and the first term inside the square brackets shows we assumed $\hat{\alpha}_{i}=\alpha_{i}$. Deviations from this assumption were studied in [44] and they mostly lead to a small re-scaling of the angular power spectrum; we also verified this effect in our simulations (see Sec. 5). We call the term inside the brackets that multiplies $\beta_{i}$ the weighted template.

Now we can compute the average weighted density $\bar{n}_{\mathrm{w}}^{i}$ in the redshift bin $i$ over the survey footprint $W(\bm{\theta})$ and do the same for the weighted template:

[TABLE]

where $A$ is the unmasked area of the survey. Assuming that the Poisson noise and the true galaxy overdensity $\delta n_{\mathrm{g}}^{i}(\bm{\theta})\equiv n_{\mathrm{g}}^{i}(\bm{\theta})-\bar{n}_{\mathrm{g},i}$ do not correlate with our weighted template difference $\delta S_{\mathrm{w}}^{i}(\bm{\theta})\equiv S_{\mathrm{w}}^{i}(\bm{\theta})-\bar{S}_{\mathrm{w}}^{i}$, $\beta_{i}$ can be estimated by:

[TABLE]

where $\delta n_{\mathrm{w}}^{i}(\bm{\theta})\equiv n_{\mathrm{w}}^{i}(\bm{\theta})-\bar{n}_{\mathrm{w}}^{i}$. In practice, we may expect a small random correlation between $\delta n_{\mathrm{g}}^{i}(\bm{\theta})$ and $\delta S_{\mathrm{w}}^{i}(\bm{\theta})$, which was investigated and corrected for in [46]. We nevertheless ignore this bias as it does not affect our analysis significantly (see Sec. 5). Lastly, the true mean galaxy density $\bar{n}_{\mathrm{g},i}$ is estimated from the data as:

[TABLE]

The estimated parameters $\alpha_{i}$, $\beta_{i}$ and $\bar{n}_{\mathrm{g},i}$, are presented in Table 1, where we see that both obscuration and contamination are clearly present in all bins. Once more, the first bin shows higher values of systematics, and the estimated error bars are much larger as well, specially for the obscuration parameter. All these difficulties presented by the first bin indicate that its systematics are not as successfully modeled as for the other bins. On top of that, this redshift range is well covered by data with better constrained photo-$z$s, larger sky coverage, and apparently less affected by systematics – namely the 2MPZ catalog [63, 34]. Thus, we decided to ignore the first redshift bin in the subsequent analysis.

Using the estimated $\alpha_{i}$, $\beta_{i}$ and $\bar{n}_{\mathrm{g},i}$, we can solve Eq. (3.2) for $n_{\mathrm{g}}^{i}(\bm{\theta})$ inside the survey footprint and estimate the true galaxy density contrast $\delta_{\mathrm{g}}^{i}(\bm{\theta})\equiv n_{\mathrm{g}}^{i}(\bm{\theta})/\bar{n}_{\mathrm{g},i}-1$:

[TABLE]

These are shown in Fig. 10. From Eq. (3.12) we see that, even if our template $S(\bm{\theta})$ and our parameter estimates all perfectly match the data, the first term inside the curly brackets might be negative due to Poisson fluctuations in $n^{i}_{\mathrm{obs}}(\bm{\theta})$ [for instance, $n^{i}_{\mathrm{obs}}(\bm{\theta})$ might be zero in a given pixel].

4 Methodology for cosmological analysis

4.1 Angular power spectrum estimation

We measured the auto and cross angular power spectra $C_{\ell}^{ij}$ of the WSC density contrast maps using the NaMaster666http://github.com/LSSTDESC/NaMaster code [64], which is based on a pseudo-$C_{\ell}$ estimator proposed in [65], and includes further improvements [66, 46]. NaMaster first computes the so-called pseudo power spectra $D^{ij}_{\ell}$, which measure the variance of the coefficients $a^{i}_{\ell m}$ used to expand the masked density contrast maps $\delta^{i}_{\mathrm{g}}(\bm{\theta})$ in spherical harmonics $Y_{\ell m}(\bm{\theta})$:

[TABLE]

The coefficients $a_{\ell m}$ represent the spherical harmonic transform of the product of the mask $W(\bm{\theta})$ and the estimate of the underlying galaxy number density contrast, according to Eq. (3.12). Therefore, when transformed to harmonic space, this product translates into a convolution of the underlying (i.e., full sky) power spectrum $F^{ij}_{\ell^{\prime}}$ and a mixing matrix $R_{\ell\ell^{\prime}}$ [66]:

[TABLE]

The mixing matrix – determined from the power spectrum of the mask (see e.g. equation 6 of [34]) – is shown in Fig. 11, where we present on the left hand side a 2D color map of the full $R_{\ell\ell^{\prime}}$, while on the left panel we present particular examples of the $R_{\ell\ell^{\prime}}$ for three values of $\ell^{\prime}$. Solving for $F^{ij}_{\ell^{\prime}}$ by inverting the mixing matrix, our final estimate of the true, full-sky angular power spectra $C^{ij}_{\ell}$ is given by:

[TABLE]

where the denominator removes the effect of the pixelization ($w_{\ell}$ is known as the pixel window function, computed with the Healpix package for our mask resolution) and $N^{ij}$ is the shot-noise contribution, given by:

[TABLE]

with $\delta_{\mathrm{K}}^{ij}$ the Kronecker delta. Table 2 shows that thanks to the large mean number density of galaxies in most redshift bins, the shot-noise term does not dominate on the scales we are interested in. There are two aspects of Eq. (4.4) worth pointing out. First, that these are the full number counts $n_{\mathrm{obs}}^{i}(\bm{\theta})$ (including contamination by stars) that contribute to the shot-noise. Second, since we weight the data by a factor $[1-\hat{\alpha}_{i}S(\bm{\theta})]$, the noise is also amplified.

We will present our measurements of the power spectrum in the form of averages in bins of width $\Delta\ell=6$, and we consider angular scales from $\ell_{0}=3$.

4.2 Estimation of cosmological parameters

In order to set constraints on cosmological parameters from the measured angular power spectra, we used the Markov Chain Monte Carlo (MCMC) code Montepython777http://baudren.github.io/montepython.html [67, 68]. We have adopted a flat $\Lambda$CDM cosmological model with the following free parameters: $h=H_{0}/(100\;\mathrm{km~{}s^{-1}~{}Mpc^{-1}})$, where $H_{0}$ is the Hubble constant; $\Omega_{\mathrm{b}}$, the baryon density parameter; $\Omega_{\mathrm{c}}$, the cold dark matter density parameter; $\ln(10^{10}A_{\mathrm{s}})$, where $A_{\mathrm{s}}$ is the amplitude of the primordial power spectrum; $n_{\mathrm{s}}$, the spectral index of the primordial power spectrum; the galaxy biases $b_{i}$ in each of the four redshift bins $i$; and the widths $\sigma_{\mathrm{z},i}$ of the $i$th bin’s spec-$z$ selection function as nuisance parameters; in total 13 parameters. We also computed the amplitude of matter density fluctuations $\sigma_{8}$ as a derived parameter. As priors we used independent Gaussian distributions for parameters poorly constrained by our data: $h$, $\ln(10^{10}A_{\mathrm{s}})$ and $n_{\mathrm{s}}$, with means and standard deviations according to Planck [69]. We point out that our narrow prior on $\ln(10^{10}A_{\mathrm{s}})$ is translated into prior constraints on $\sigma_{8}$ as well.

The theoretical $C^{ij}_{\ell}$s used in the MCMC were computed by CLASS888http://class-code.net [70, 71], including contributions from redshift space distortions and non-linear structure growth [72, 73]. Together with the Planck priors, we applied the extra restriction $\Omega_{\mathrm{b}}>0.0065$, to avoid a crash in CLASS caused by reionization happening at too late times. To describe the spec-$z$ distribution of galaxies in each photo-$z$ bin, we used their Gaussian approximation; the impact of assuming a generalized Lorentzian distribution for the spec-$z$s is discussed in Sec. 7.4.

To reduce the number of possible sources of systematic deviations between measured and true parameters, we restricted our analysis to linear scales. More specifically, for each $C^{ij}_{\ell}$ we only kept the multipoles for which the contribution from the non-linear part of the matter power spectrum was smaller than 10%. Table 2 shows this maximum $\ell$ for each angular power spectrum.

The binned $C^{ij}_{\ell}$s were combined into a single data vector with 169 entries. We assumed a Gaussian likelihood, with the covariance (see Fig. 12) extracted from 2500 lognormal mock catalogs of the WSC, described in Sec. 5.1. The $C^{ij}_{\ell}$s are measured from these mocks, following the same procedure as used for the data (described in Secs. 3.3 and 4.1). The inverse of the covariance matrix estimated from the simulations has been corrected for the finite sample used to estimate it, following [74].

The sampling of the posterior was performed with 10 Markov chains of 100,000 entries, each one starting from a different point in the parameter space. The convergence of the chains was verified using the Gelman-Rubin diagnostic $R$, and all parameters presented $R<1.15$.

5 Validation of the analysis

To make sure that all the steps of our analysis work properly, we used the Full-sky Lognormal Astro-fields Simulation Kit (FLASK999http://www.astro.iag.usp.br/~flask) [75] to create hundreds of mock WSC catalogs including the obscuration and contamination effects described in Sec. 3. These mocks are described in detail in Sec. 5.1. The results obtained with our methodology are shown in Sec. 5.2.

5.1 Mock WSC maps

To create lognormal realizations that reproduce the basic statistical properties of the data, FLASK may use as input: a set of auto- and cross-$C_{\ell}$s for all the redshift slices being simulated; the mean density of sources $\bar{n}_{\mathrm{g},i}$; and the angular completeness $f(\bm{\theta})$, which includes the mask, obscuration and/or angular variations in density. FLASK output is a set of Healpix maps that include Poisson noise and that reproduce all one- and two-point statistics as required, including cross-correlations between redshift bins. It can also create realizations of stellar densities by assuming $C_{\ell}^{\mathrm{star}}=0$, setting $f(\bm{\theta})$ to our template and Poisson sampling from it.

As the expected densities (angular-position dependent), we used $\bar{n}_{\mathrm{g},i}[1-\alpha_{i}S(\bm{\theta})]$ for galaxies and $\beta_{i}S(\bm{\theta})$ for stars, where $S(\bm{\theta})$ is our template (Fig. 8) and the remaining parameters are set according to Table 1. As input galaxy $C_{\ell}$s we used smooth fits to the measured cleaned $C_{\ell}$s, presented later in Fig. 16, up to $\ell_{\mathrm{max}}=640$. The maps of galaxies and stars, for each bin, were summed into a single map of sources, to which we applied our analysis. Fig. 13 shows an example of such a realization.

5.2 Recovery of mock properties

In this Section we test our methodology to account for the systematic effects and the likelihood machinery on our mock catalogs.

5.2.1 Obscuration and contamination parameters

The first step in our method is to estimate the obscuration $\alpha_{i}$ and contamination $\beta_{i}$ parameters for all bins $i$, along with the mean galaxy density $\bar{n}_{\mathrm{g},i}$. Table 3 shows the mean and standard deviation of these parameters, recovered from 2500 WSC mock realizations. It shows that our method traces quite well the true values in each of the bins, being able to detect obscuration and contamination at no less than 5.7$\sigma$ and 3.9$\sigma$, respectively. It does, however, produce a detectable average bias of $\sim 3\%$ on $\alpha$ and $\sim 7\%$ on $\beta$ (i.e. overcorrecting, on average, both effects). This bias is likely caused by our simplifying assumption that the cosmological variance $\sigma_{\mathrm{cosmo}}^{2}(s)$ for counts-in-pixels, computed from a pixel set, does not depend on the template’s isocontour (Eq. 3.6), which happens to vary both in shape and in area. As density fluctuations are spatially correlated, variances computed in a smaller and compact set will, in general, be smaller. Since, as we show next, this bias does not significantly affect our results, we considered the method’s performance adequate despite of its presence.

5.2.2 Angular power spectra

In order to verify the impact of our cleaning method on the $C_{\ell}$s, we processed 750 simulations in three different ways101010Due to the large amount of time required to measure 10 $C_{\ell}$s in three different ways from each simulation, we did not use all 2500 at this step.: (i) we ignored the obscuration and contamination present in them (equivalent to assuming $\alpha=\beta=0$ even though they are not); (ii) we only corrected for contamination (equivalent to assuming $\alpha=0$); and (iii) we applied the full cleaning method described in Sec. 3.3. We then measured the angular power spectra using our standard procedure (Sec. 4.1). Fig. 14 shows the average of 750 estimated $C_{\ell}$s obtained after applying the three different cleaning processes in comparison to the true $C_{\ell}$ in the case of the cross-correlation between redshift bins 3 and 4, as an example. The general characteristics of the results are the same for all bin pairs.

We see that, in our simulations, the typical effect of contamination is to increase the power on the largest scales ($\ell\lesssim 20$), showing at $\ell\lesssim 10$ deviations with respect to the true spectrum larger than the $1\sigma$ error bars. This is clearly a systematic effect we need to correct for in the data. Similarly, we see that correcting only for contamination gives an estimate of the angular power spectrum which is biased low. This bias, also reported in [44], amounts to $\sim 15\%$ and reaches the size of the error bars on small scales. After correcting for obscuration, this bias is eliminated (although there seems to be a small positive and constant bias of $\sim 1\%$, possibly caused by an overestimated $\alpha$). This shows that when measuring cosmological or astrophysical parameters controlling the amplitude of the angular clustering signal (e.g. galaxy bias and primordial spectrum amplitude), correcting for obscuration is a key step to obtain accurate results. This comparison shows that, assuming our modeling of the systematic is correct, its application significantly improves the accuracy of the estimated power spectrum, leaving only possible residual bias which amounts to no more than $2\%$ of the error bars.

5.2.3 Cosmological parameters

Once we have quantified that our cleaning method properly corrects the measured power spectrum, we use the MCMC methodology to recover the input cosmological and nuisance parameters from cleaned $C_{\ell}$s. The input $C_{\ell}$s for the FLASK code were computed by CAMB sources111111http://camb.info/sources [76], including not only contributions from redshift space distortions and non-linearities, but also from a minimal neutrino configuration and all available effects (e.g. gravitational lensing and general relativistic corrections). Furthermore, the computation of the input $C_{\ell}$s did not use the Limber approximation and adopted higher precision than the $C_{\ell}$s sampled by the MCMC. In this way, we verified that any potential deviations in the recovered parameters are not caused by approximation schemes or by the omission of small physical effects. The input $C_{\ell}$s were computed assuming a flat $\Lambda$CDM model and Gaussian spec-$z$ selection function, and the simulations included the effects of obscuration and contamination. The procedure followed here mimics that applied to the real data.

Figure 15 shows that all marginal posterior distributions for all the parameters measured from two independent simulations are compatible with their input values, thus validating our pipeline. We also point out that, while the covariance matrix used in our likelihood was estimated from simulations based on smooth fits to the measured $C^{ij}_{\ell}$s, the simulated data used in this test was generated from a different set of $C^{ij}_{\ell}$s, thus attesting that our covariance matrix does not bias the results. The reader might notice that posterior distribution for e.g. $\sigma_{\mathrm{z},3}$ from one simulation is displaced from its true value. In this case, we point out that we have 13 independent parameters and, in fact, we should expect a few to deviate from their expected values by chance. The second simulation shows that this is indeed the case. Moreover, we see that the ability to pinpoint $\Omega_{\mathrm{b}}$ is affected by the specific realization of noise and cosmological fluctuations.

6 Results

In Fig. 16 we show the binned angular power spectra measured from the WSC data, both cleaned and uncleaned for contamination and obscuration. We see that the effect of cleaning the WSC is reasonably small on most scales, indicating that galaxy catalogs containing systematics can still provide cosmological information. Despite the smallness of the corrections, we point out two relevant aspects of them: first, applying them significantly reduces the power on the largest scales ($\ell\lesssim 10$) in many cases, and as we showed in Sec. 5.2, this is caused by contamination removal; second, some cleaned $C_{\ell}$s show an almost constant-factor increase in power at all scales, caused by obscuration correction, as we have previously discussed (see Fig. 14).

Figure 16 also shows the best-fit model to the cleaned data (black lines), as well as $C_{\ell}$s obtained using cosmological parameters given by Planck [69] and nuisance parameters fitted “by eye” as a rough reference (green lines). By comparing them we see that WSC $C_{\ell}$s present steeper slopes than the fiducial cosmology on the linear scales, and that these slope differences do not seem to be restricted to the largest scales only ($\ell\lesssim 20$).

A comparison of the best-fit $\chi^{2}$ obtained under different data processing strategies is presented in Table 4. In every case, the best-fit $C_{\ell}$s and $\chi^{2}$ were obtained for that particular dataset by fitting all 13 parameters to it. We see that correcting the WSC count maps for obscuration and contamination improves the fit, but the model remains not very representative of the data. If we ignore multipoles $\ell<15$, the fit improves significantly and reaches a $p$-value level that could be considered adequate (larger than $5\%$). This may indicate that the largest scales are likely still contaminated even after our corrections, possibly due to a mismatch between the real contamination and our template. If we push the multipole cut to higher values, the agreement does not improve significantly, suggesting that the best-fit is not constrained by the largest scales after this point.

Figure 17 shows the 1D and 2D posterior distributions for the cosmological parameters that were not directly restricted by Planck priors, for the first three approaches listed in Table 4, and compares them with the values according to Planck. In any case, the posteriors disagree with Planck values by more than $2\sigma$. The most discrepant parameter is $\Omega_{{\mathrm{c}}}$, which WSC data push to lower values. This is caused by the steep slope of WSC $C_{\ell}$s, as shown in Fig. 16: lower cold matter density leads to less small-scale clustering, tilting the power spectrum and making it redder.

Fig. 17 also shows that the power reduction at the largest scales produced by our cleaning method significantly shifts the posterior closer to the Planck values without increasing its size. This demonstrates, using real data, the importance and advantages of our cleaning method. The removal of the largest scales ($\ell<15$) moves the posterior even closer to the Planck values but increases its width. This suggests that the stellar density template may be improved, or that the largest scales are affected by extra effects besides stellar obscuration and contamination. In either case, it is clear that identifying and taking these effects into account might be a better approach than ignoring the largest scales altogether, as the latter strategy would lead to worse precision.

7 Robustness tests and search for further systematics

To search for other potential sources of systematic effects in the WSC data, we applied a series of modifications to our analysis and checked their impact on the measured $C_{\ell}$s, using the full-sky cleaned $C_{\ell}$s with $\ell\geq 15$ as our starting point. One of these tests, already mentioned in Sec. 6, consisted in ignoring scales $\ell<21$, with the only impact on the results being that of making the posterior distributions wider.

7.1 Different stellar maps and masks

When considering possible residual systematics related to stellar density, we first raised two hypotheses, viz. (i) that problematic regions near the Galactic bulge could be biasing our results, and (ii) that the Gaussian smoothing applied to the Gaia map could have erased Galactic structures at intermediate scales ($\ell\gtrsim 30$) that were not removed by our fiducial template. In order to verify these claims, we created the template $S(\bm{\theta})$ shown in the left panel of Fig. 18, consisting in an unsmoothed map with regions of stellar density three times the mean masked, and applied it in the removal of obscuration and contamination effects.

Another possibility we considered was that the map obtained from Gaia was not representing the contamination in the WSC. Hence, and based on the cross-match with SDSS presented in Sec. 3.1, we created a template with the simple exponential model given by Eq. (3.1) (central panel in Fig. 18). Finally, a third template was created, in which the model fitted to the matched SDSS stars had an exponential disk (with an extra longitudinal dependence) and a Gaussian bulge (Fig. 18, right panel). Note that despite the fact that this approach is based on a direct estimate of the contamination, it relies on an extrapolation from the SDSS footprint to the whole sky.

In all cases, the cleaning process reduced the power on the largest scales and only significantly affected multipoles $\ell\lesssim 20$. The two templates based on Gaia led to very similar outcomes and, for most $C^{ij}_{\ell}$s, improved the results (in terms of large scale power reduction) more than those based on the SDSS cross-match.

7.2 Uncleaned data with $\ell\geq 15$

Although we have shown in Secs. 5.2.2 and 5.2.3 that our cleaning method does not bias the results (and, in particular, does not tilt the $C^{ij}_{\ell}$s), this demonstration assumed that our stellar density template had the correct shape. To test whether our cleaning method could be biasing the results due to a mismatch between the template and the true stellar density distribution, we estimated the posterior distribution of the cosmological parameters from the uncleaned data, while ignoring the largest scales ($\ell<15$).

Once more, the posteriors remained highly compatible with those for the cleaned $\ell\geq 15$ case shown in Fig. 17. This testifies that our cleaning method only affects the large-scale modes $\ell<15$ and is not responsible for the disagreement between WSC and Planck.

7.3 Mixing matrix deconvolution and $H_{0}$ prior

As stated in Sec. 4.1, we estimated the angular power spectra of the unmasked sky by inverting Eq. (4.2), that is, we deconvolved the mixing matrix from the data by multiplying the observed pseudo-$C_{\ell}$s $D_{\ell}^{ij}$ by the inverse $R_{\ell^{\prime}\ell}^{-1}$ of the mixing matrix. As this process combines data from different scales, it could transfer to smaller scales the power coming from systematic effects on the largest scales. An alternative analysis that should not suffer from this effect would be to compare the data pseudo-$C_{\ell}$s (after removing shot-noise and the pixel window function) with theoretical $C_{\ell}$s convolved with the mixing matrix $R_{\ell^{\prime}\ell}$. If no anomalous power transfer happens during the mixing matrix deconvolution in the first approach, both methods should yield the same results. We ran the simple test of comparing the ratio between deconvolved data $C_{\ell}$s and fiducial $C_{\ell}$s with the ratio between data pseudo-$C_{\ell}$s and the same (but convolved) fiducial $C_{\ell}$s. Both ratios agreed to very good extent, demonstrating the robustness of our analysis against forward/backward modeling schemes.

Another analysis choice impact we decided to test was that of the $H_{0}$ prior. It is known that the Hubble constant inferred from Planck ($H_{0}=67.31\pm 0.96\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$) is in tension with direct measurements using type Ia supernovae ($H_{0}=73.24\pm 1.74\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$) [77]. We verified if the Planck $H_{0}$ prior could be responsible for the observed tension between $\Omega_{\mathrm{c}}$ estimated from WSC and Planck data by replacing it with direct measurements prior. Although this change displaced the WSC posteriors, it actually increased the tension on $\Omega_{\mathrm{c}}$ slightly.

7.4 Redshift distributions

As mentioned in Sec. 2.1, the WSC photo-$z$ errors can also be described by a generalized Lorentzian distribution, which results (from Eq. 2.2) in a spec-$z$ distribution in each photo-$z$ bin that can also be well modeled by the same distribution. To test whether the Gaussianity assumption used in our analysis could bias the results, we employed generalized Lorentzian functions as spec-$z$ distributions, and used it to generate a model of $C_{\ell}$s and sample the posterior distribution of the constrained parameters.

Generalized Lorentzian functions have an extra parameter $a$ that controls the importance of its tails (see Eq. 2.1). We ran the MCMC under two configurations: (i) we let $a_{i}$ for each $i$ redshift bin vary freely between 1.5 and 4.4; and (ii) we fixed $a_{i}$ to the values estimated in [38]. In both cases, the posterior on the cosmological parameters remained quite similar to the one in Fig. 17.

We also investigated the scenario in which ignoring the photo-$z$ correction due to the north-south hemisphere asymmetry (used in our main analysis; see [33] for details) could bias our results: using the photo-$z$ computed directly by ANN$z$, we performed again the redshift binning, the cleaning and the measurement of the $C_{\ell}$. As a result, no significant change was observed.

7.5 Ignoring cross-$C_{\ell}$s

Our main analysis extracted cosmological information not only from the power spectra in each redshift bin but also from the cross-correlation between different bins. These cross-$C_{\ell}$s are more sensitive to the spec-$z$ distribution of the galaxies in each bin than the auto-$C_{\ell}$s, since the overlap in the true redshift distributions is one of the main sources of signal in their case. To test if the cross-$C_{\ell}$s could be biasing our results, we repeated the procedure to estimate cosmological parameters using only the auto-$C_{\ell}$s. Although we noticed no significant displacement in the mean values, the width of the 1D posterior distributions are larger, showing that the cross-$C_{\ell}$s do contain extra information.

7.6 Different color and magnitude cuts

The WSC catalog originally had a variable color cut dependent on the distance from the Galactic center, which we removed by applying a more stringent constant cut all over the sky (see Sec. 2.1). We computed the WSC $C_{\ell}$s using this original color cut and verified that it results in more power on the largest scales than when using our cut $W1-W2>0.2$, without much change on smaller scales. This effect is expected since a position-dependent color cut creates a variable galaxy selection function, likely making their expected number and bias also position-dependent. In other words, the homogeneity in the number counts discussed by [33] is achieved by the original WSC color cut not by removing stars at low latitudes but by removing both stars and galaxies.

Another potential source of systematics is angular variations of the survey depth, caused by changes in observing conditions (in SuperCOSMOS) or total exposure times (number of observations in WISE). To test for the presence of such features, we applied more conservative magnitude cuts, going from the WSC fiducial limits $W1<17$, $B<21$ and $R<19.5$ to $W1<16.5$, $B<20.7$ and $R<19.1$. Apart from a constant factor increase in power in all $C_{\ell}$s, expected as a consequence of the larger bias usually associated with brighter galaxies, no significant change was seen in the angular power spectra.

7.7 Analysis by hemispheres

In the full-sky analysis above we have found that, although our corrections for stellar systematics do shift the WSC posterior in the direction of Planck (with respect to the no-correction case), there is still considerable discrepancy between WSC and Planck parameters. As there is no significant evidence from other LSS studies at similar redshifts that this should be physical [78, 34, 79, 35], here we analyze a possible observational effect related to the WSC catalog. Namely, the underlying SuperCOSMOS data are based on measurements from two telescopes which had different effective passbands [53, 80]. The division between the two ‘hemispheres’ is at $\delta_{1950}=2.5^{\circ}$. As noted in [33], despite color-dependent corrections applied in SuperCOSMOS to match north and south [80], within the WSC flux limits (Sec. 2.1) there is still residual difference in mean galaxy surface density, the north one being up to 4% larger than in the south. As this might be a sign of a genuine systematic in the data, we decided to perform the same analysis as for the full sky also for the two hemispheres. Similar tests were done for 2MPZ [34], which also incorporates SuperCOSMOS as one of the input samples, but no significant differences between resulting cosmological parameters were found. We note that the second input sample included in WSC – WISE – does present non-uniformities in the data but not of hemispherical nature, but rather related to the satellite’s polar orbit and Moon avoidance maneuvers121212See http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec3_4.html. Since these non-uniformities are much more involved, we did not investigate their possible impact on our results.

Both hemispheres present similar behavior as the full-sky analysis with respect to the impact of our cleaning process and removal of multipoles $\ell<15$. When including all scales with $\ell\geq 3$ and ignoring contamination, both hemispheres display a bad $\chi^{2}$ and a clear tension with the results from Planck. The cleaning process improves the fit and reduces the tension with Planck by lowering the power on the largest scales ($\ell<9$). Fig. 19 shows the northern $C_{\ell}$s before and after cleaning. The southern hemisphere – which already presented less excess of power on these scales – is better improved than the northern hemisphere and actually reaches an acceptable $\chi^{2}$, although both remain inconsistent with Planck. When we remove the scales with $\ell<15$, the $\chi^{2}$ becomes acceptable for both hemispheres, but only the northern hemisphere agrees with Planck. This may indicate that different systematics dominate each hemisphere. The $\chi^{2}$ values are shown in Table 4 and the posterior distributions for each hemisphere under the $\ell\geq 15$ cut are shown and compared with the full-sky posteriors in Fig. 20. The values of all 13 parameters associated with these posteriors are presented in Table 5.

8 Conclusions and summary

In this work we have tested a method for estimating and correcting for stellar obscuration and contamination in galaxy redshift catalogs, aiming at extracting cosmological information therefrom. In particular, we have implemented it on the WISE$\times$SuperCOSMOS (WSC) photo-$z$ dataset [33]. We first exposed in Sec. 3 the presence of both systematic effects in WSC, a feature that will be common – albeit with a smaller amplitude – in every photometric survey, including LSST, J-PAS and Euclid. In that Section we also proposed a model for this feature that ties both systematic effects together with a single stellar density template (Eq. 3.2), and showed that such an assumption is supported by the data (Figs. 7 and 9).

The fact that both effects – obscuration and contamination – are linked by the template demands a strategy other than simple correlation [62, 46] to quantify them. We proposed the usage of variance measurements in bands of constant stellar density as a way of determining obscuration (Eq. 3.7), which later can be used to estimate the contamination and real galaxy density (Eqs. 3.10 and 3.11). We verified that this strategy is able to detect and quantify all these parameters with less than 10% bias (Table 3). Although Figs. 14 and 15 show that the effect of any residual bias is negligible in the rest of our analysis, we propose that this cleaning method can be improved if we take into account that galaxy densities are spatially correlated.

It is worth pointing out that one may disentangle the effects of obscuration and contamination directly on the power spectra – since they result in different signatures (see Sec. 5.2.2) – by treating them as nuisance parameters and fitting them together with the cosmology. In this process, it is relevant to not only analyse the data $C_{\ell}$s but also the cross power spectra between data and stellar density template. The contamination level may be better constrained by it, whereas the measurement of obscuration might benefit from a cross-$C_{\ell}$ between the template and the observed density difference squared, $\delta n_{\mathrm{obs}}^{2}$. Testing this approach should be part of a future work.

All of our methodology was validated by applying it to simulations created by the FLASK code [75] and comparing the measured power spectrum and the constrained parameters with their true input values (see Sec. 5). That analysis demonstrated that applying our cleaning method to data containing obscuration and contamination significantly improved the cosmological parameter inference. This method was also applied to real WSC data, and Table 4 and Fig. 17 show that it improved the power spectra fit at $\ell<15$ and increased the agreement between the cosmology derived from WSC power spectra and from the Planck satellite.

Despite the improvement obtained both in the full-sky and in separate hemispheres by removing contamination and obscuration, the constraints derived from the angular clustering of WSC still disagree with those derived from the Planck mission if we include the largest scales. In fact, the minimum $\chi^{2}$ shows that the standard model is a bad fit to the data in the full-sky and northern hemisphere case (see Table 4) even for parameter values significantly different from Planck. The fit is significantly improved if we ignore multipoles $\ell<15$ (reaching a $p$-value of $\sim$6% for the full-sky case, which could be considered adequate), proving that much of the deviation between data and model happens at the largest scales. The reasons for this excess of power is still unclear and should be further investigated. Possible explanations are that (i) stellar density templates, used to remove contamination, must be improved to better represent WSC true distribution; and that (ii) important systematic effects not related to stars – possibly Galactic extinction and photometric/observational errors – need to be taken into account. We do point out, however, that we have tested three other templates (Fig. reffig:alt-templates) and they have led to similar (or more biased) results.

As the scales $\ell\sim 15$ are close to the power spectrum turn-over at redshifts probed by WSC, we could claim that the excess of power at scales larger than this comes from physical effects not considered by our cosmological model, e.g.: by non-Gaussian initial conditions caused by certain inflationary models [81]; by the fact that galaxy forms at matter density peaks, which could lead to a scale-dependent bias [30]; or by alternative dark energy models that produce different redshift space distortions [82]. These effects, however, are expected to be small and probably should wait for the investigation of other possible systematics sources. It is worth pointing out that an excess of power is constantly seen on the largest scales of galaxy surveys [47, 34, 35], and due to its correlation with stellar density and other reasons, these modes are mostly discarded. This is unfortunate as the physical effects mentioned above remain untested.

We have tested our data for several potential sources of systematics using the full-sky data limited to $\ell\geq 15$ as a reference. In Sec. 7 we verified that our results are not affected by photo-$z$ error distributions, cleaning strategy, color and magnitude cuts, photo-$z$ estimations, mixing matrix deconvolution, the use of cross-$C_{\ell}$s and the restriction to $\ell\geq 21$. We did find, however, a hemispherical pattern in the data. Fig. 20 shows that the independent analyses of the two equatorial hemispheres lead to incompatible cosmological parameters: while the analysis restricted to $\ell\geq 15$ and to the northern hemisphere agrees with Planck (with $\Omega_{c}=0.254^{+0.018}_{-0.022}$, $\sigma_{8}=0.854^{+0.042}_{-0.045}$ and $\Omega_{\mathrm{b}}<0.065$ at 95% confidence limit when combined with Planck priors on $H_{0}$, $A_{s}$ and $n_{s}$), the same analysis in the southern hemisphere does not. This difference in extracted cosmological parameters is traced back to different slopes in the measured power spectra, with the southern and full-sky $C_{\ell}$s being steeper than the northern one. Considering it is unlikely that such a large difference has a cosmological origin – especially since cosmological parameters have been extracted from the southern hemisphere before and shown to agree with Planck [34, 79] – this indicates that systematics affecting scales $\ell\gtrsim 15$ are located in the southern part of WSC data.

Investigating what is the actual source of these systematics is beyond the scope of the present work. One possible explanation is residual miscalibration of SuperCOSMOS data in the south. Namely, the input photographic material for that sample was based on observations from two telescopes (POSS-II in the north and UKST in the south) with effectively different passbands. Despite the efforts to match the passbands using specific color terms, anchored to SDSS as the calibrator [80], differences may still exist. As the WSC sample is selected using maximum apparent magnitude cuts, in particular in the SuperCOSMOS $B$ and $R$ bands, this may propagate into inconsistent sample selections in the two hemispheres. Indeed, already in [33], a difference in effective source densities between north and south was reported.

However, the most plausible explanation for the north-south difference in large-scale power comes from residual uncertainties in photometric zero points. Where direct SDSS calibration was lacking, plate zero points were given an initial estimate from mean optical-2MASS color, and these values were then perturbed in an attempt to optimize agreement in plate overlap zones. It seems likely that this process will yield zero points that are imperfect at the level of a few hundredths of a magnitude [83], imprinting a pattern of $5\times 5\,\mathrm{deg^{2}}$ tiles on the galaxy distribution with an rms variation in surface density of several percent. This is the right order of magnitude to cause the offsets that we measure (J. Peacock, private communication). The second input dataset of WSC – WISE – is unlikely to be at the origin of the hemispherical differences as any non-uniformity in that survey is related mostly to telescope’s polar orbit and Moon avoidance strategy. Thus, our present findings seem to suggest that the observed asymmetry is due to residual systematics in the UKST data.

Acknowledgments

The authors thank Vinicius Placco and Laura Sampedro for advice on the use of Gaia data, Filipe Abdalla for comments related to power spectrum measurements and cosmological parameter extraction and John Peacock for valuable information on SuperCOSMOS and useful advice on data analysis and interpretation.

We thank the Wide Field Astronomy Unit (WFAU) at the Institute for Astronomy, Edinburgh, for storing the WISE$\times$SuperCOSMOS data, available for download from http://ssa.roe.ac.uk/WISExSCOS. We have also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is http://www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work has made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by the Brazilian agency FAPESP (grant 2009/54006- 4) and the INCT-A. HSX acknowledges FAPESP Brazilian funding agency for the financial support. ABA acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the Severo Ochoa program SEV-2015-0548. MB is supported by the Netherlands Organization for Scientific Research, NWO, through grant number 614.001.451 and by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12.