Cosmological constraints from a joint analysis of cosmic microwave background and spectroscopic tracers of the large-scale structure

TL;DR

This paper demonstrates that combining CMB data with spectroscopic galaxy and quasar surveys through a joint analysis significantly improves constraints on cosmological parameters, including neutrino mass and dark energy properties.

Contribution

The study introduces a multi-probe joint analysis method using pseudo-spectra that enhances parameter constraints beyond individual datasets.

Findings

Joint analysis constrains neutrino mass and dark energy equation of state.

Combining datasets improves parameter constraints compared to separate analyses.

Method demonstrates potential for future cosmological studies with multi-probe data.

Abstract

The standard model of cosmology, {\Lambda}CDM, is the simplest model that matches the current observations, but it relies on two hypothetical components, to wit, dark matter and dark energy. Future galaxy surveys and cosmic microwave background (CMB) experiments will independently shed light on these components, but a joint analysis that includes cross-correlations will be necessary to extract as much information as possible from the observations. In this paper, we carry out a multi-probe analysis based on pseudo-spectra and test it on publicly available data sets. We use CMB temperature anisotropies and CMB lensing observations from Planck as well as the spectroscopic galaxy and quasar samples of SDSS-III/BOSS, taking advantage of the large areas covered by these surveys. We build a likelihood to simultaneously analyse the auto and cross spectra of CMB lensing and tracer overdensity maps before running Monte-Carlo Markov Chains (MCMC) to assess the constraining power of the combined analysis. We then add the CMB temperature anisotropies likelihood and obtain constraints on cosmological parameters ($H_0$, $\omega_b$, $\omega_c$, ${\ln10^{10}A_s}$, $n_s$ and $z_{re}$) and galaxy biases. We demonstrate that the joint analysis can additionally constrain the total mass of neutrinos ${\Sigma m_{\nu}}$ as well as the dark energy equation of state $w$ at once (for a total of eight cosmological parameters), which is impossible with either of the data sets considered separately. Finally, we discuss limitations of the analysis related to, e.g., the theoretical precision of the models, particularly in the non-linear regime.

Tables

Table 1. Table 1: Cosmological constraints from the full joint analysis of CMB temperature, CMB lensing and BOSS spectroscopic tracers with the optimistic cut (see section   3.3 ). We test four different cosmological models and use the same flat priors in all cases. Constraints on the tracers’ biases are also given. The last three rows give constraints on the derived parameters Ω m subscript Ω m \Omega_{\mathrm{m}} , σ 8 subscript 𝜎 8 \sigma_{8} and τ 𝜏 \tau . This table gives the median and asymmetric 68% error bars.

Parameter [unit]	$\Lambda$CDM	$\Lambda$CDM + $\Sigma m_{\nu }$	$w$CDM	$w$CDM + $\Sigma m_{\nu }$	Prior (flat)
$H_0$ [$\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}$]	${68.77}_{-0.85}^{+0.87}$	${67.6}_{-1.4}^{+1.3}$	${69.7}_{-6.3}^{+10.2}$	${65.6}_{-4.8}^{+6.7}$	$[50,100]$
${\omega }_b$	$\left({225.1}_{-1.9}^{+1.8}\right)\times {10}^{-4}$	$\left(224.6\pm 1.8\right)\times {10}^{-4}$	$\left(224.8\pm 1.9\right)\times {10}^{-4}$	$\left(224.8\pm 1.9\right)\times {10}^{-4}$	$[0.01,0.03]$
${\omega }_c$	$\left(116.7\pm 1.9\right)\times {10}^{-3}$	$\left(117.2\pm 1.9\right)\times {10}^{-3}$	$\left(117.0\pm 2.2\right)\times {10}^{-3}$	$\left({116.8}_{-2.1}^{+2.0}\right)\times {10}^{-3}$	$[0.05,0.2]$
$\ln {10}^{10}{A}_s$	${3.092}_{-0.034}^{+0.036}$	${3.104}_{-0.036}^{+0.038}$	${3.085}_{-0.049}^{+0.050}$	${3.119}_{-0.046}^{+0.047}$	$[2.5,3.5]$
$n_s$	$\left({976.2}_{-5.1}^{+5.0}\right)\times {10}^{-3}$	$\left({975.0}_{-5.2}^{+4.8}\right)\times {10}^{-3}$	$\left({975.3}_{-5.6}^{+5.8}\right)\times {10}^{-3}$	$\left({976.0}_{-5.5}^{+5.6}\right)\times {10}^{-3}$	$[0.8,1.1]$
$z_{\mathrm{re}}$	${10.2}_{-1.7}^{+1.6}$	$10.8\pm 1.7$	${9.9}_{-2.6}^{+2.3}$	${11.5}_{-2.2}^{+2.0}$	$[0,20]$
$\Sigma m_{\nu }$ [eV]	–	${0.16}_{-0.10}^{+0.12}$	–	$0.17\pm 0.11$	$[0,10]$
$w$	–	–	$-{1.04}_{-0.32}^{+0.21}$	$-{0.93}_{-0.22}^{+0.16}$	$[-2,0]$
$b_{\text{LOWZ}}$	${1.837}_{-0.033}^{+0.034}$	${1.880}_{-0.051}^{+0.059}$	${1.826}_{-0.110}^{+0.100}$	${1.91}_{-0.11}^{+0.10}$	$[0,10]$
$b_{\text{CMASS}}$	$2.086\pm 0.032$	${2.130}_{-0.050}^{+0.058}$	${2.083}_{-0.055}^{+0.059}$	${2.146}_{-0.072}^{+0.077}$	$[0,10]$
$b_{\text{QSO}}$	${2.19}_{-0.44}^{+0.45}$	${2.24}_{-0.45}^{+0.46}$	${2.21}_{-0.43}^{+0.44}$	$2.24\pm 0.46$	$[0,10]$
${\Omega }_m$	$0.296\pm 0.011$	${0.309}_{-0.015}^{+0.018}$	${0.288}_{-0.068}^{+0.058}$	${0.327}_{-0.058}^{+0.054}$	–
${\sigma }_8$	${0.822}_{-0.010}^{+0.011}$	${0.805}_{-0.023}^{+0.019}$	${0.830}_{-0.055}^{+0.082}$	${0.787}_{-0.046}^{+0.060}$	–
$\tau$	${0.082}_{-0.019}^{+0.020}$	$0.088\pm 0.020$	${0.078}_{-0.027}^{+0.028}$	${0.096}_{-0.025}^{+0.026}$	–

Equations

$$A(\hat{\mathbf{n}})={\int }_0^{\infty }dz{W}^A(z)\delta (\chi (z)\hat{\mathbf{n}},z),$$

$$\text{expectationvalue}A(\hat{\mathbf{n}})B^{\ast }({\hat{\mathbf{n}}}^{\prime })=\sum _{\ell =0}^{\infty }\frac{(2l+1)}{4\pi }{P}_{\ell }(\hat{\mathbf{n}}\cdot {\hat{\mathbf{n}}}^{\prime })C_{\ell }^{AB}$$

$$C_{\ell }^{AB}$$

$$\times \int dk\frac{2}{\pi }{k}^{2}P(k,z,z^{\prime })j_{\ell }(k\chi (z))j_{\ell }(k\chi (z^{\prime })),$$

$$C_{\ell }^{AB}={\int }_0^{z_{\ast }}dz\frac{H(z)}{c{\chi }^2(z)}{W}^A(z)W^B(z)P\left(k=\frac{\ell +1/2}{\chi (z)},z\right)$$

$$\mathbf{\alpha }=-2{\int }_0^{z_{\ast }}dz\frac{c}{H(z)}\frac{\chi (z_{\ast })-\chi (z)}{\chi (z)\chi (z_{\ast })}{\nabla }_{\hat{\mathbf{n}}}\Psi (\chi \hat{\mathbf{n}},z).$$

$$W^{{\kappa }_{\mathrm{CMB}}}(z)=\frac{3}{2}\frac{{\Omega }_{\mathrm{m}}{H}_0^2}{c}\frac{(1+z)}{H(z)}\chi (z)\frac{\chi (z_{\ast })-\chi (z)}{\chi (z_{\ast })},$$

$$g(\hat{\mathbf{n}})={\int }_0^{\infty }dz{W}^g(z)\delta (\chi \hat{\mathbf{n}},z).$$

$$W^g(z)=b(z)\frac{dn}{dz}+\frac{3{\Omega }_{\mathrm{m}}}{2c}\frac{H_0^2}{H(z)}(1+z)\chi (z)(5s-2)g(z)$$

$$g(z)={\int }_z^{z_{\ast }}d{z}^{\prime }\left(1-\frac{\chi (z)}{\chi (z^{\prime })}\right)\frac{dn}{d{z}^{\prime }}.$$

$$s={\frac{d\log N(<m)}{dm}|}_{m=m_{\max }}.$$

$$⟨{\tilde{C}}_{\ell }^{AB}⟩=\sum _{{\ell }^{\prime }}{M}_{\ell {\ell }^{\prime }}^{AB}{C}_{{\ell }^{\prime }}^{AB},$$

$$M_{\ell\ell^{\prime}}^{AB}=\frac{2\ell+1}{4\pi}\sum_{\ell^{\prime\prime}}(2\ell^{\prime\prime}+1)\mathcal{W}_{\ell^{\prime\prime}}^{\,AB}\left(\begin{array}[]{ccc}\ell&\ell^{\prime}&\ell^{\prime\prime}\\ 0&0&0\end{array}\right)^{2}.$$

$${\mathcal{W}}_{{\ell }^{\prime \prime }}^{AB}=\frac{1}{2{\ell }^{\prime \prime }+1}\sum _m{w}_{{\ell }^{\prime \prime }m}^A(w_{{\ell }^{\prime \prime }m}^B{)}^{\ast },$$

$$w_{\ell m}^A=\int d\hat{\mathbf{n}}{\mathcal{W}}^A(\hat{\mathbf{n}})Y_{\ell m}^{\ast }.$$

$$C_i=\frac{N_{\mathrm{obs},\mathrm{i}}}{N_{\mathrm{targ},\mathrm{i}}}.$$

$$w_{\mathrm{tot}}=w_{\mathrm{star}}{w}_{\mathrm{seeing}}(w_{\mathrm{cp}}+w_{\mathrm{noz}}-1),$$

$${\delta }_p=\frac{N_{w(p)}}{\overline{N}}-1.$$

$${\delta }_p=\frac{\text{ifrac}{N}_p{C}_{s(p)}}{\overline{N}}-1,$$

$${\hat{\tilde{C}}}_{\ell }^{AB}=\frac{1}{2\ell +1}\sum _{m=-\ell }^{m=+\ell }{\tilde{A}}_{\ell m}{\tilde{B}}_{\ell m}^{\ast }.$$

$$⟨{\hat{\tilde{C}}}_{\ell }^{AB}⟩=\sum _{{\ell }^{\prime }}{M}_{\ell {\ell }^{\prime }}^{AB}{C}_{{\ell }^{\prime }}^{AB}+{\delta }_{AB}{\tilde{N}}_{\ell }^A$$

$${\tilde{N}}_{\ell }^A=\sum _{{\ell }^{\prime }}{M}_{\ell {\ell }^{\prime }}^{AA}{N}_{{\ell }^{\prime }}^A.$$

$${\hat{\tilde{C}}}_{\ell }^{AB}=\frac{1}{2\ell +1}\sum _{m=-\ell }^{m=+\ell }{\tilde{A}}_{\ell m}{\tilde{B}}_{\ell m}^{\ast }-{\delta }_{AB}\frac{1}{N_{\mathrm{sim}}}\sum _{i=1}^{N_{\mathrm{sim}}}{\hat{\tilde{N}}}_{\ell }^{A,i},$$

$$\mathrm{Cov}\left({\tilde{C}}_{\ell }^{AB},{\tilde{C}}_{{\ell }^{\prime }}^{CD}\right)=$$

$$+\sqrt{D_{\ell }^{AC}{D}_{{\ell }^{\prime }}^{AC}{D}_{\ell }^{BD}{D}_{{\ell }^{\prime }}^{BD}}{\text{Y}}_{\ell {\ell }^{\prime }}^{ABCD}$$

$$D_{\ell}^{AB}=\left\{\begin{array}[]{ll}C_{\ell}^{AB}&\text{if }A\neq B\\ C_{\ell}^{AA}+N_{\ell}^{A}&\text{if }A=B\\ \end{array}\right.,$$

$$\mathcal{L}\left({\widetilde{\text{mn@boldsymbol}C}}_{\ell }^{\mathrm{obs}}|b_g,{\Theta }_{\mathrm{cosmo}}\right)=\frac{1}{(2\pi {)}^{n/2}{|\text{Cov}|}^{1/2}}{e}^{-{\chi }^2/2},$$

$${\chi }^2={\left({\widetilde{\text{mn@boldsymbol}C}}^{\mathrm{obs}}-{\widetilde{\text{mn@boldsymbol}C}}^{\mathrm{th}}\right)}^{\mathrm{T}}{\left[\text{Cov}\right]}^{-1}\left({\widetilde{\text{mn@boldsymbol}C}}^{\mathrm{obs}}-{\widetilde{\text{mn@boldsymbol}C}}^{\mathrm{th}}\right),$$

$$n_p\sim \mathrm{Poisson}\left({\lambda }_p\right)\text{with}{\lambda }_p=\overline{N}\left(1+{\delta }_p\right),$$

$${\tilde{\rho }}_{\ell }\equiv \frac{{\tilde{C}}_{\ell }^{\mathrm{syst}\times {\delta }_g}}{\sqrt{{\tilde{D}}_{\ell }^{{\delta }_g\times {\delta }_g}{\tilde{D}}_{\ell }^{\mathrm{syst}\times \mathrm{syst}}}},$$

$$\mathrm{Cov}\left({\tilde{\rho }}_{\ell },{\tilde{\rho }}_{\ell }\right)\approx \frac{1}{f_{\mathrm{sky}}^2(2{\ell }^{\prime }+1)}{M}_{\ell {\ell }^{\prime }}^{gg}$$

Full text

††pubyear: 2017††pagerange: Cosmological constraints from a joint analysis of cosmic microwave background and spectroscopic tracers of the large-scale structure–E

1 Introduction

The large amount of cosmological data collected in the last few decades has been shedding light on the content of the Universe. Assuming General Relativity (GR) and the cosmological principle, the combination of different cosmological probes, such as type Ia supernovæ, primary anisotropies of the cosmic microwave background (CMB), and large-scale structure (LSS) information, among others, indicates that the universe is almost flat, is dominated today by a dark energy (DE) component driving the current accelerated expansion phase of the Universe, and has some form of cold dark matter in addition to baryons and radiation (Planck Collaboration et al., 2016e). The flat $\Lambda$CDM model is currently the simplest model compatible with the data of these combined probes.

We are reaching a precision era in cosmology and we may be able, in the near future, to distinguish between various cosmological models and achieve a better understanding of the fundamental nature of the DE and DM components. Upcoming photometric and spectroscopic galaxy surveys such as the Large Synoptic Survey Telescope (LSST, LSST Science Collaborations and LSST Project et al., 2009), Euclid (Refregier et al., 2010), the Wide Field Infrared Survey Telescope (WFIRST, Spergel et al., 2013) and the Dark Energy Spectroscopic Instrument (DESI, Levi et al., 2013; Schlegel et al., 2009), aim at shedding light on those questions by probing the matter density field with ground-breaking precision. They will provide the data necessary for a deeper investigation of $\Lambda$CDM and its competitors, hopefully allowing us to distinguish them. Additionally, secondary anisotropies of the CMB due to gravitational lensing, the thermal (tSZ) and kinetic (kSZ) Sunyaev-Z’eldovich effects and the integrated Sachs-Wolfe (ISW) effect encode much information about dark matter and dark energy (Peiris & Spergel, 2000). Therefore, future CMB experiments, such as the Simons Observatory (Suzuki et al., 2016) and the Stage-IV CMB experiment (CMB-S4, Abazajian et al., 2016), will provide valuable complementary observations.

While various observations from multiple telescopes will provide exquisite and hopefully complementary data sets (Jain et al., 2015; Rhodes et al., 2017), they will all observe the same sky, i.e. the same underlying matter density field. Therefore, the observables they will measure are potentially statistically correlated. In this context, the cross-correlation between cosmological probes of different experiments yields new information, that is less prone to biases since different experiments are assumed to have uncorrelated noise and independent systematic effects. This correlation needs to be taken into account in the joint statistical analysis of multiple data sets in order to properly extract as much information as possible from it, without underestimating error bars on cosmological parameters (Rhodes et al., 2015). If this makes the analysis more demanding, the outcome is expected to provide stronger constraints on e.g. dark energy, dark matter, the total mass of neutrinos (Pearson & Zahn, 2014) or primordial non-gaussianities (Takeuchi et al., 2012).

Initially, some of the best-explored cross-correlation information was that from CMB and galaxy surveys in order to measure the ISW signal (Crittenden & Turok, 1996; Boughn et al., 1998; Fosalba et al., 2003; Cabré et al., 2006; Moura-Santos et al., 2016). But over the last decade, many different cross-correlation signals have been detected, combining various probes: the CMB anisotropies themselves, the CMB lensing potential, galaxy clustering, cosmic shear from the observations of galaxy weak lensing, etc. In particular, correlations of the gravitational lensing of the CMB with positions of galaxies (Hirata et al., 2008; Ho et al., 2008; Bleem et al., 2012; Sherwin et al., 2012; Geach et al., 2013; Bianchini et al., 2015; Fornengo et al., 2015; Allison et al., 2015; Baxter et al., 2016; Giannantonio et al., 2016; Harnois-Déraps et al., 2016; Bianchini et al., 2016) and lensing of galaxies (Hand et al., 2015; Kirk et al., 2016) have been measured for various surveys. These measurements can provide unbiased estimates of galaxy biases, which encode the link between baryonic and dark matter, or the shear multiplicative bias (Vallinotto, 2012; Liu et al., 2016; Schaan et al., 2017). Finally, they have also been used to detect new signals, e.g. the first detection of CMB lensing by cross-correlation with the NRAO VLA sky survey (Smith et al., 2007), the tSZ effect (Hajian et al., 2013; Hill & Spergel, 2014), the kSZ effect (Hand et al., 2012; Schaan et al., 2016), and the position-dependent Lyman-$\alpha$ power spectrum (Doux et al., 2016).

Since cross-correlation signals are reaching high signal-to-noise ratio, joint analyses, i.e. multi-probe analyses that exploit cross-correlations, are rapidly developing. In particular, joint analyses based on real-space correlation functions or power spectra have been used to test the consistency of cosmological constraints derived from different observations and to cross-calibrate nuisance parameters. Joachimi & Bridle (2010) forecasted that the joint analysis of galaxy weak lensing and galaxy density, including cosmic shear, galaxy clustering and the galaxy-galaxy lensing cross-correlation, could, at once, self-calibrate intrinsic alignments and constrain parameters of the cosmological model (see Krause et al. (2016) as well). Since, several analyses were published that make use of the correlation between galaxy lensing and either galaxy density or CMB lensing on available data. Baxter et al. (2016) used SPT and DES-SV lensing data together with tracers of the large-scale structure and took advantage of the low systematic level of these angular cross-correlations functions to infer cosmological constraints, demonstrating their consistency across data sets. Kwan et al. (2017) used the galaxy clustering and galaxy-galaxy lensing signals of the DES-SV data to obtain robust constraints on $\sigma_{8}$ and $\Omega_{\mathrm{m}}$. Singh et al. (2017) combined the lensing and clustering of SDSS galaxies with CMB lensing from Planck to constrain clustering and shear biases and measure distance ratios, found to be consistent with Planck predictions. Very recent studies combined these three correlation functions in a single analysis: van Uitert et al. (2017) with KiDS111The Kilo-Degree Survey, de Jong et al. (2013). and GAMA222Galaxies And Mass Assembly (Driver et al., 2009), a database of low-redshift surveys spanning the electromagnetic spectrum from radio waves to the ultraviolet domain. galaxies, Joudaki et al. (2017) with KiDS and BOSS galaxies (using the quadrupole of the power spectrum as well), DES Collaboration et al. (2017) with the first year of DES data (forecasts were obtained in Park et al. (2016) and the robustness of the likelihood analysis pipeline was tested in Krause et al. (2017)), and Miyatake et al. (2015); More et al. (2015) used the overlapping area between CFHT and BOSS galaxies. Taking a step further, Nicola et al. (2016, 2017) used, in a fully joint analysis, information from CMB temperature, CMB lensing, photometric surveys (both galaxy positions and lensing) and distance measurements from supernovae and direct $H_{0}$ measurements. Finally, Baxter et al. (2018) recently laid out the methodology for a joint analysis of five DES and SPT+Planck data two-point auto- and cross-correlation functions.

In this work, we aim at performing a joint analysis of Planck CMB data with the spectroscopic LSS tracers of the Baryon Oscillation Spectroscopic Survey (SDSS-III/BOSS, Dawson et al., 2013; Reid et al., 2016; Pâris et al., 2017) based on power spectra, taking advantage of the large overlap between these surveys and the large redshift range of the BOSS samples, and developing an independent pipeline that is able to incorporate observables with different masks (thus maximizing signal-to-noise). To this end, we build a Gaussian joint likelihood of auto- and cross- pseudo power spectra of large-scale structure probes –here, CMB lensing and the contrast densities of tracers– and implement it as part of the public Numerical Cosmology333http://numcosmo.github.io library (NumCosmo, Dias Pinto Vitenti & Penna-Lima, 2014). In particular, we use the Planck 2015 lensing map (Planck Collaboration et al., 2016f) and the three spectroscopic samples of BOSS –that is, LOWZ, CMASS and the uniform quasar sample (QSO). We then run several Markov Chain Monte-Carlo (MCMC) analyses to extract constraints on cosmological parameters. We complete our analysis by adding the likelihood of the auto power spectrum of CMB temperature anisotropies from Planck (thus neglecting the small CMB-LSS direct correlations sourced by the ISW effect, Planck Collaboration et al., 2016d, e), which allows us to perform a joint statistical analysis of CMB and LSS probes. We demonstrate the performance of such an analysis to constrain the $\Lambda$CDM model at first, fitting the six base parameters ($H_{0}$, $\omega_{\mathrm{b}}$, $\omega_{\mathrm{c}}$, $A_{\mathrm{s}}$, $n_{\mathrm{s}}$ and $z_{\mathrm{re}}$), and then extensions of this base model, fitting the total mass of neutrinos, ${\Sigma m_{\nu}}$, and the DE equation of state, $w=p/\rho$, separately and then jointly. Finally, we perform the joint analysis using two different cut-offs at small scales in order to quantify their impact on cosmological constraints.

Note that our work follows an approach similar to that presented in Nicola et al. (2016), which is also based on power spectra. However, our goals differ as we aim at constraining extensions of the $\Lambda$CDM model, which we do using, in part, different data sets (we use spectroscopic data from BOSS) and we propose a different formalism based on pseudo power spectra and a semi-analytical covariance matrix evaluated at each set of cosmological parameters during the MCMC analysis, rather than a fixed covariance matrix (either simulation-based or computed analytically for a fiducial model). We also provide a new, public implementation of multi-probe analysis as part of the public NumCosmo library.

The paper outline is as follows. In section 2, we develop the theoretical formalism, considering the Limber formula, to compute the pseudo angular power spectra of CMB lensing and galaxy overdensities. In section 3, we describe the data used in this work and prepare the galaxy and quasar density maps. In section 4, we develop the estimators and construct the likelihood used in this work, and then we perform the validation and null tests. In section 5, we detail the results of the statistical analyses and present constraints on cosmological parameters as well as their potential limitations related to theoretical uncertainties (e.g. related to the non-linear power spectrum), contamination of cross-spectra or data cuts. Finally, we draw conclusions in section 6.

2 Theoretical background

2.1 Angular power spectra

The matter distribution of the Universe is traced by cosmological probes such as galaxies, quasars (QSOs) and CMB lensing, among others. Their projected random fields on the observed direction $\hat{\mathbf{n}}$ can be written as

[TABLE]

where $W^{A}(z)$ is the kernel function of an observable $A$ (kernels of the probes used in this work are plotted in fig. 1), $\delta=\delta\rho/\rho$ is the matter density contrast, $\rho$ is the matter density, and $\chi(z)$ is the comoving distance at redshift $z$. The fields $A(\hat{\mathbf{n}})$ and $B(\hat{\mathbf{n}})$, associated with galaxy or QSO density contrast and CMB lensing, are assumed to be statistically homogeneous and isotropic. Therefore the correlation function $\expectationvalue{A(\hat{\mathbf{n}})B^{*}(\hat{\mathbf{n}}^{\prime})}$ only depends on ${\hat{\mathbf{n}}\cdot\hat{\mathbf{n}}^{\prime}}$ and can be expanded as

[TABLE]

where $P_{\ell}$ are the Legendre polynomials, which defines the angular power spectrum $C_{\ell}^{AB}$.444Expanding the fields in spherical harmonics, $A{(\hat{\mathbf{n}})=\sum_{{\ell}m}A_{{\ell}m}Y_{{\ell}m}(\hat{\mathbf{n}})}$, this implies that ${\expectationvalue{A_{{\ell}m}B_{{\ell}^{\prime}m^{\prime}}}=C_{\ell}^{AB}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}}$.

Using the inverse Fourier transform of the matter density field in eq. (1) and substituting it into eq. (2), we obtain

[TABLE]

where $P(k,z,z^{\prime})$ is the matter power spectrum and $j_{\ell}$ are the spherical Bessel functions. In this work, we adopt the Limber approximation (Limber, 1953; LoVerde & Afshordi, 2008) and assume that spatial sections of the Universe are flat555It is worth noting that Limber approximation can be also applied for curved space universe, see Lesgourgues & Tram (2014), for instance., thus eq. (3) becomes

[TABLE]

where $c$ is the speed of light and $H(z)$ is the Hubble parameter. This approximation is valid when $P(k,z,z^{\prime})$ varies slowly in comparison with the Bessel functions. In particular, the CMB lensing spectra are accurate for $\ell>10$ (Lesgourgues & Tram, 2014). The selection functions of galaxies and quasars are wider than the largest scales probed666For a given sample, the largest scale probed is ${\chi_{\mathrm{max}}\sim\pi/k_{\mathrm{min}}}$ with ${k_{\mathrm{min}}=(\ell_{\mathrm{min}}+1/2)/\chi(z_{\mathrm{eff}})}$, where $\chi(z_{\mathrm{eff}})$ is the comoving distance at the mean redshift of the sample and $\ell_{\mathrm{min}}=20$, see section 3. For LOWZ, CMASS and QSO, these scales are of order $110220630\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$, while the selection functions have widths of order $1\,080860970\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$. for the spectroscopic tracers used here. Therefore, we can safely make use of this approximation to compute theoretical power spectra, which are integrals of the matter power spectrum weighted by the kernel functions corresponding to each observable. We detail this in the following sections.

2.2 Cosmic microwave background gravitational lensing

The trajectories of the CMB photons are disturbed by the matter distribution such that, among other effects, the observed anisotropies of the temperature field in a direction $\hat{\mathbf{n}}$ correspond to the unlensed field deflected by $\boldsymbol{\alpha}$, i.e., $\tilde{T}(\hat{\mathbf{n}})=T(\hat{\mathbf{n}}+\boldsymbol{\alpha})$ (Lewis & Challinor, 2006). Assuming the small-angle Born approximation, $\boldsymbol{\alpha}$ comprises the variations in the gravitational potential $\Psi$ along the line of sight from today’s observer to the last scattering surface at redshift $z_{*}$, i.e.,

[TABLE]

The remapping of the CMB temperature anisotropies due to the weak lensing effect is described at lowest order by the convergence ${\kappa_{\mathrm{CMB}}}=-\frac{1}{2}\nabla_{\hat{\mathbf{n}}}\boldsymbol{\alpha}$. Given that on small angular scales $\nabla^{2}_{\hat{\mathbf{n}}}\simeq\nabla^{2}$ (Jain et al., 2000) and using the Poisson equation and eq. (1), we obtain the CMB convergence kernel

[TABLE]

where $\Omega_{\mathrm{m}}$ and $H_{0}$ are the present-day matter density and Hubble constant, respectively.

2.3 Large-scale structure tracers

Similarly to the CMB lensing, the galaxy or quasar overdensity in the direction $\hat{\mathbf{n}}$ is a function of $\delta(\chi\hat{\mathbf{n}},z)$, namely

[TABLE]

The kernel $W^{g}(z)$ is given by (Peiris & Spergel, 2000; Bonvin & Durrer, 2011)

[TABLE]

where

[TABLE]

The function $b(z)$ is the linear bias relating the galaxy overdensity to the matter overdensity at large scales as $\delta_{g}(\chi\hat{\mathbf{n}},z)=b(z)\delta(\chi\hat{\mathbf{n}},z)$, and $dn/dz$ is the normalised redshift distribution of the tracers, which also contains the survey selection function. The second term in eq. (8) is due to the effects of gravitational lensing, with two opposing contributions – the dilation of the apparent surveyed volume and the magnification bias effect for flux-limited samples777Lensed galaxies may appear brighter than they are and pass the luminosity threshold., where

[TABLE]

$N(<m)$ denotes the cumulative count of objects with a magnitude smaller than $m$ and the derivative is estimated at the faint end of the catalog (Scranton et al. (2005); Hui et al. (2007)). This term can be neglected when using LOWZ and CMASS samples, but it is relevant for quasars (Chisari & Dvorkin, 2013). Following Scranton et al. (2005), we use $s_{\textsf{QSO}}=0.2$ throughout this analysis.

2.4 Pseudo spectra

Many galaxy and CMB surveys cover only a fraction of the sky due to, for example, the limited field of view or galactic contamination, among others. In order to properly account for the partial sky coverage in the calculation of the angular power spectra, we define the mask function associated with the field $A(\hat{\mathbf{n}})$ as being $\mathcal{W}^{A}(\hat{\mathbf{n}})$ with value 1 if the direction $\hat{\mathbf{n}}$ lies in the observed region, and 0 otherwise.

The cross-pseudo spectrum of observables $A$ and $B$ is thus defined as the cross spectrum of the cut-sky fields ${\tilde{A}(\hat{\mathbf{n}})=\mathcal{W}^{A}(\hat{\mathbf{n}})A(\hat{\mathbf{n}})}$ and ${\tilde{B}(\hat{\mathbf{n}})=\mathcal{W}^{B}(\hat{\mathbf{n}})B(\hat{\mathbf{n}})}$, and its expected value can be related to the (true) full-sky cross spectrum in eq. (3) by (Brown et al., 2005)

[TABLE]

where $M_{\ell\ell^{\prime}}^{AB}$ is the mixing matrix which is given in terms of the Wigner-$3j$ symbols

[TABLE]

The cross spectra of the masks are

[TABLE]

where

[TABLE]

The mixing matrix introduces a scaling factor equal to $\left(f_{\rm sky}^{A}f_{\rm sky}^{B}\right)^{1/2}$, i.e., the geometric mean of the observed sky fractions for the observables $A$ and $B$, respectively, since the form of the masked function is constant (Hivon et al., 2002). It also couples the multipoles $\ell$ and $\ell^{\prime}$, that would be otherwise uncorrelated, especially at large scales. It is computed analytically for each pair of observables (see fig. 2).

3 Data

3.1 Planck data

Planck888http://sci.esa.int/planck/ (Planck Collaboration et al., 2016a), the fourth satellite to survey the CMB over the full sky, was launched on May 14th, 2009. Its scientific payload comprised two instruments: the Low Frequency Instrument (Mandolesi et al., 2010), which observed for four years in bands at 30, 44 and 70 GHz, and the High Frequency Instrument (Lamarre et al., 2010), which observed for almost two-and-a-half years in bands at 100, 143, 217, 353, 545 and 857 GHz.

In this work, we use Planck data for both the primordial CMB temperature anisotropies (Planck Collaboration et al., 2016b) and the CMB lensing (Planck Collaboration et al., 2016f). For temperature anisotropies, we use the two Planck likelihood codes: Plik for the high multipoles, $\ell\geq 30$, and Commander for low multipoles, $\ell<30$ (see Planck Collaboration et al., 2016d, for details on components separation).

We also use the CMB lensing convergence map from the Planck 2015 data release (Planck Collaboration et al., 2016f). The Planck Collaboration provides the convergence map999The convergence map and mask files are publicly available at http://pla.esac.esa.int/pla/. in the Healpix101010http://healpix.jpl.nasa.gov/ (Górski et al., 2005) format, with resolution parameter $N_{\rm side}=2048$, and the corresponding binary mask, with a sky fraction $f_{\rm sky}=0.67$. Lensing potential maps were reconstructed from foreground-cleaned temperature and polarisation maps, obtained from the SMICA code. These were used to form five quadratic estimators $\hat{\phi}^{TT}$, $\hat{\phi}^{TE}$, $\hat{\phi}^{EE}$, $\hat{\phi}^{EB}$ and $\hat{\phi}^{TB}$, combined into a minimum-variance estimator (see Planck Collaboration et al., 2016f, for specificities about the Planck reconstruction). The Wiener-filtered convergence map is shown in fig. 3a with its mask.

The full-sky lensing power spectrum was evaluated by the Planck collaboration following the methods detailed in Appendix A of Planck Collaboration et al. (2016f), which we briefly describe now: first, the pseudo power spectrum of the masked, reconstructed lensing map (itself using optimized temperature map masks) is upweighted by $1/f_{\rm sky}$. Then, different contributions must be subtracted: the Gaussian noise from the disconnected part of the four-point function (the N0 term), the non-Gaussian noise from the connected contribution (the N1 term) and the subdominant contribution of shot-noise from unresolved point sources (the PS term). Finally, a term (the MC term) derived from the mismatch between input and output power spectra derived from simulations corrects for errors in the normalization, mask-related mode mixing and computation of the N1 term. Finally, the likelihood of the binned, full-sky power spectrum includes small, linear corrections to take into account the dependence of the reconstructed lensing map and the N1 term on the fiducial cosmological model.

Given our goal to jointly analyze CMB lensing with tracers of the LSS, we cannot simply use the lensing likelihood code provided by the Planck Legacy Archive. However, an estimate of the unbinned (pseudo) lensing auto-power spectrum has not been released yet. Therefore, considering that redoing the Planck analysis goes well beyond the scope of this paper, we measured the lensing pseudo-power spectrum from the released map and estimated the dominant N0 and N1 noise terms from the simulated reconstructed maps (see section 4.1.1) and neglected the small MC term, the largely subdominant PS term and fiducial model correction terms. Indeed, the MC term partly corrects for differences between the full-sky and pseudo-power spectra while we only rely on the pseudo spectrum, and the corrections applied at the likelihood level are very small in the multipole ranges used for cosmology. Therefore, while a simplification, our method should provide reasonable results. Indeed, in appendix A, we compare cosmological constraints obtained from CMB temperature and lensing using our pipeline and that provided by the Planck collaboration and demonstrate that we retrieve unchanged cosmological constraints, except for a $0.5\sigma$ shift in the $A_{\mathrm{s}}-z_{\mathrm{re}}$ degeneracy, still below the statistical error. While the reader should bear in mind this small difference, we consider this comparison to justify our simplifications given the purposes of this work.

3.2 SDSS-III/BOSS data

The spectroscopic samples from the Baryon Oscillations Spectroscopic Survey (BOSS, Dawson et al., 2013) consists of two galaxy catalogs named LOWZ and CMASS and one quasar catalog, a subset of which has a uniform selection function. They are extensively described in Reid et al. (2016) and we only summarize relevant information in this section.

3.2.1 Luminous Red Galaxies: LOWZ & CMASS

LOWZ contains Luminous Red Galaxies (LRG) at low redshift ($z\lesssim 0.4$), and it aims at a constant number density of about $\bar{n}\sim$3\text{\times}{10}^{-4}\text{,}{\mathrm{\mathit{h}}}^{3}\text{,}{\mathrm{Mpc}}^{-3}$$ over the redshift range $[0.1,0.4]$. This is done using a redshift dependent magnitude cut. In this work, we use the twelfth data release (DR12) and select galaxies in the range 0-0.4, which contains $383\,876$ galaxies. The CMASS sample contains galaxies at higher redshifts $0.4\lesssim z\lesssim 0.8$ with a constant stellar mass in this redshift range. The twelfth data release contains $849\,637$ galaxies in the redshift range 0.4-0.8 used in this work. The normalised redshift distributions of the two samples are shown on fig. 1 (multiplied by their respective biases).

BOSS’s spectroscopic fibres are plugged into tiles of diameter $3\text{\,}\mathrm{\SIUnitSymbolDegree}$ to observe predetermined targets. The combined footprints of all tiles can be decomposed into non-overlapping sky sectors. Because of the finite size of fibres, galaxies closer than 62’ may not be observed even after multiple observations of the same field. The pipeline may also fail in determining the redshift of some galaxies (especially the faintest ones). Therefore, for each sector $i$, the completeness is defined as the ratio of observed galaxies with a measured redshift to the number of targets lying in that same sector

[TABLE]

The completeness maps are defined in the Mangle software111111See http://space.mit.edu/~molly/mangle/. format and are converted to Healpix maps with resolution parameter $N_{\rm side}=2048$. The mask functions of the galaxy samples are obtained by assigning 1 to pixels where the completeness is above 75% and then removing small areas that were vetoed for bad photometry, bright objects and stars and instrumental constraints, such as fibre centerposts and fibre collisions.

In order to correct for completeness, each galaxy is thus assigned a weight

[TABLE]

where $w_{\mathrm{star}}$ and $w_{\mathrm{seeing}}$ correct for non-cosmological fluctuation in the target selection due to stellar density (only for the CMASS sample) and atmospheric seeing, $w_{\mathrm{cp}}$ corrects for fibre collisions and $w_{\mathrm{noz}}$ corrects for redshift failures.

3.2.2 Quasars

The selection function (over the sky) of the full quasar sample of BOSS is not uniform due to the observing strategy, hence we shall use the so-called CORE sample which contains QSOs with redshift $z\geq 2.15$ that were uniformly selected by the XDQSO algorithm (Bovy et al., 2011). There are 94971 quasars in the CORE sample of DR12 within this redshift range. The completeness is computed using the BOSSQSSOMASK software121212See http://faraday.uwyo.edu/~admyers/bossqsomask/. from Eftekharzadeh et al. (2015) and is then combined with the veto mask to build the mask of the quasar density map.

3.2.3 Building the maps

For the LOWZ and CMASS samples, we build Healpix maps with resolution parameter $N_{\rm side}=2048$, where for each pixel $p$,

[TABLE]

$N_{w(p)}=\sum_{i\in p}w_{i}$ is the number of galaxies in pixel $p$ counted with their weights and $\overline{N}=\frac{1}{N_{\rm pix}}\sum_{p=1}^{N_{\rm pix}}N_{w(p)}$ is the mean pixel count (where the sum runs only on pixels in the observed area, i.e. where the mask function is equal to 1).

For the quasars, there is no weight provided in the BOSS DR12 catalog and the density map is computed as

[TABLE]

where $N_{p}$ denotes the number of QSOs lying in pixel $p$, $C_{s(p)}$ is the completeness of the sector $s(p)$ where the pixel $p$ lies and $\overline{N}$ denotes the mean pixel count (up-weighted by completeness) in the observed area.

The angular densities of the samples are $\bar{n}_{\textsf{LOWZ}}=$150\text{\times}{10}^{3}\text{,}{\mathrm{sr}}^{-1}, $\bar{n}_{\textsf{CMASS}}=$300\text{\times}{10}^{3}\text{\,}{\mathrm{sr}}^{-1} and $\bar{n}_{\textsf{QSO}}=$36\text{\times}{10}^{3}\text{,}{\mathrm{sr}}^{-1}$$. The maps of the estimated overdensity for the three samples are shown in figs. 3b, 3c and 3d.

3.3 CMB lensing–large-scale structure correlations data

In this work, we will use, in addition to CMB temperature data, auto and cross spectra of CMB lensing from Planck and spectroscopic tracers from BOSS. More precisely, we will use the auto-pseudo spectra of the CMB lensing map $\kappa_{\mathrm{CMB}}$ and of the density contrast maps of the LOWZ and CMASS samples. We also use the pseudo-cross-spectra of the CMB lensing map with the three LSS tracers. The collection of these six spectra (shown in fig. 4) will henceforth be referred to as “CMB lensing-LSS correlations”, and denoted “CMB lensing $\otimes$ BOSS tracers” in the figures.

We do not use the auto spectrum of the QSO map because it is completely shot-noise-dominated in multipole space. We do not use the galaxy cross spectra because their redshift ranges do not overlap and the cross spectra should therefore be zero in the Limber approximation (which we check in the next section).

We use different multipole ranges for the different spectra and describe here the cuts that were applied. At very large scales, the Limber approximation breaks (LoVerde & Afshordi, 2008), RSD becomes non-negligible (Padmanabhan et al., 2007) and observational systematics become more difficult to handle (see section 4.3). For these reasons, we chose a minimum multipole of ${\ell_{\rm min}=20}$ common to all spectra involving LSS tracers. At small scales, uncertainties in the non-linear power spectrum severely constrain the use of angular power spectra. We therefore considered two cut-offs in Fourier space, using a fixed smallest scale $k_{\rm max}$ that translates into a maximum multipole $\ell_{\rm max}\approx\ifrac{k_{\rm max}}{\chi\left(z_{\rm eff}\right)}$ where $\chi\left(z_{\rm eff}\right)$ is the comoving distance to the mean redshift of the samples. We will use a conservative cut-off at $k_{\rm max}=$0.1\text{,}{\mathrm{Mpc}}^{-1} and a more optimistic one at $k_{\rm max}=$0.15\text{\,}{\mathrm{Mpc}}^{-1}. For LOWZ, CMASS and QSO, we find respectively $\ell_{\rm max}=$120200480 and $\ell_{\rm max}=$180320720. At these scales, the updated halofit model reaches close to percent-level precision, even for models including massive neutrinos (Smith et al., 2003; Bird et al., 2012; Takahashi et al., 2012). Moreover, the scale-independent bias approximation has been found to work reasonably well for scales down to $20\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$, with variations of less than 5% (Cresswell & Percival, 2009; Rodríguez-Torres et al., 2016). We shall therefore, for both cuts, be in the regime where these assumptions are safe.

The CMB lensing power spectrum was measured in the range $\ell=8-2048$ (Planck Collaboration et al., 2014), but due to potential errors in modelling large-scale survey anisotropies and Gaussian noise at small scales (the N0 disconnected component), we only consider multipoles $\ell=40-400$ and $\ell=20-500$ for the conservative and optimistic cuts. The high-$\ell$ cut is justified by the facts that Planck, as of the 2015 results, has the best sensitivity among CMB experiments at $\ell<500$ and that it includes additional signal-to-noise while discarding the range $500<\ell<700$ where the estimated power spectrum is significantly lower that expected for the $\Lambda$CDM model. Uncertainties in the N1 term becomes more important at $\ell\gtrsim 600$ and we shall therefore also be in a safe regime for both cuts considered here.

The six observed pseudo spectra are shown in fig. 4 together with theoretical curves for our best-fit biases ($b_{\textsf{LOWZ}}=1.831\pm 0.048$, $b_{\textsf{CMASS}}=2.077\pm 0.029$ and $b_{\textsf{QSO}}=2.21\pm 0.44$) and with fixed cosmology.131313Parameters’ values are fixed at the best-fit cosmology for Planck “TT,TE,EE+lowP+lensing+ext” (Planck Collaboration et al., 2016e) for the flat $\Lambda$CDM model with a total mass of neutrinos ${{\Sigma m_{\nu}}}=$0.06\text{,}\mathrm{eV}$$. We report detections for the CMB lensing-galaxy density cross-correlations of 4.7 $\sigma$, 13.9 $\sigma$ and 10.6 $\sigma$ for LOWZ, CMASS and QSO, respectively.

4 Methodology

This section describes the method to analyse the data. We choose to use the pseudo-power spectrum formalism, which takes into account partial sky coverage by appropriately scaling the theoretical full-sky power spectra and comparing them to the observed pseudo-spectra. This forward-modelling method has the advantage that it provides simple ways to obtain unbiased prediction, without the need to reverse the effects of the mask on the observation, an operation that can be difficult and unstable given the complexity of the masks. Moreover, it naturally deals with observables using different masks, thus maximizing the signal-to-noise ratio of the Fourier coefficients of cross-spectra. The drawback is that the covariance matrix is impossible to compute analytically, even in the Gaussian field approximation: instead, we must either estimate it with Monte-Carlo simulations, or, as we do here, use a semi-analytical approximation. We perform several tests to validate this method and the statistical pipeline used in the next section. Finally, we also perform null tests searching for possible contamination of the power spectra by effects due to the masks or the SDSS photometry.

4.1 Likelihood

4.1.1 Pseudo spectra estimator

Pseudo spherical harmonic coefficients $\tilde{A}_{\ell m}$ and $\tilde{B}_{\ell m}$ (for $A,B\in\{\textsf{LOWZ},\textsf{CMASS},\textsf{QSO},{\kappa_{\mathrm{CMB}}}\}$) of the four maps are estimated with the map2alm function of Healpix, corrected for the Healpix pixel window function and summed to give an estimator of the pseudo spectra

[TABLE]

These pseudo spectra have a noise contribution and an expectation value

[TABLE]

where $\tilde{N}_{\ell}^{A}$ is the noise pseudo-spectrum of the measured field $A(\hat{\mathbf{n}})$, which needs to be subtracted. It is assumed here that different observables have uncorrelated noise, i.e. that noise cross spectra are null (for both full-sky and pseudo).

In principle, noise pseudo spectra can be computed using the mixing matrix and eq. (11).

[TABLE]

However, the sum over $\ell^{\prime}$ runs from 0 to infinity, so, in practice, it has to be cut at some maximum multipole $\ell_{\max}$. But convergence is not guaranteed, since noise spectra are increasing functions of $\ell$ (they are quasi-constant for galaxies and grow like $\sim\ell^{2}$ for CMB lensing). Therefore, we used instead simulated noise maps for CMB lensing and a shuffling technique for spectroscopic tracers positions (similar to that used in Nicola et al. (2016)).

For CMB lensing, we used the 100 simulated lensing reconstruction maps provided by the Planck Legacy Archive.141414See http://pla.esac.esa.int/pla/. Given a known, full-sky, input convergence map (to be masked) and a masked, reconstructed convergence map, one can compute the difference of the pseudo spectra in order to obtain an estimate of the noise pseudo power spectrum $\tilde{N}_{\ell}^{\kappa}$, which is then averaged over realisations of the simulation.

For the clustering of spectroscopic tracers, the full-sky shot noise spectrum is constant, equal to $N_{\ell}=1/\bar{n}$ where $\bar{n}$ is the angular density of objects (weighted and expressed in steradian-1). However, weights associated with each object in the galaxy spectroscopic samples to compensate incompleteness slightly increase the noise level (by up to 8%)151515This can be intuited by noting that observing two galaxies with weight 1 as more information that observing one galaxy with weight two. In Ho et al. (2012), this was dealt with by adding an additive free term.. Therefore, we randomly reposition objects within the masks, keeping their weights. This operation breaks the spatial, cosmological correlation and therefore the cosmological contribution to the spectrum, leaving only Poisson noise with appropriate weighting. Density maps are then built according to the procedure described in section 3.2.3 and their pseudo-spectra $\tilde{N}_{\ell}^{g}$ are evaluated with Healpix. This process is repeated one thousand times and the noise spectrum is averaged over realizations. The high resolution of the density maps allows us to measure the angular spectra at very large multipole values (up to $\ell=3N_{\rm side}-1=6143$) where the cosmological signal becomes negligible with respect to shot-noise. We find excellent agreement of our noise estimator and the one measured on the real density map. Moreover, we validate this process by fitting the measured spectra with a free additive constant as in Ho et al. (2012), found to be consistent with zero, within $1\sigma$ error bars.

Our estimator thus reads

[TABLE]

where $\hat{\tilde{N}}_{\ell}^{A,i}$ is the estimated pseudo-noise spectrum of simulation number $i$. The pseudo spectra used in this work are shown in fig. 4.

4.1.2 Covariance matrix and likelihood

The covariance matrix of the pseudo spectra used in this work assumes that the density field is Gaussian for the scales exploited in the analysis. It is computed using an extension of Efstathiou’s symmetrisation approximation (Efstathiou, 2004) following Brown et al. (2005) and is given by

[TABLE]

with

[TABLE]

where $C_{\ell}^{AB}$ and $N_{\ell}^{A}$ are the full-sky theoretical and noise spectra. $\textbf{{X}}_{\ell\ell^{\prime}}^{abcd}$ and $\textbf{{Y}}_{\ell\ell^{\prime}}^{abcd}$ are two matrices depending only on the masks of observables ${A,B,C,D}$, determined to arbitrary precision by a Monte-Carlo (MC) simulation (see appendix B for more details). The covariance matrix (estimated for a fiducial cosmology) is shown in fig. 5.

A Gaussian likelihood is used for the stacked pseudo spectra vector

[TABLE]

where

[TABLE]

$\tilde{\mn@boldsymbol{C}}^{\mathrm{obs}}$ is the stacked vector of observed pseudo spectra (see fig. 4) and $\tilde{\mn@boldsymbol{C}}^{\mathrm{th}}$ is the stacked vector of theoretical pseudo spectra computed from the Limber approximation (see eq. (4)) and multiplied by the mixing matrices. The covariance matrix Cov is that of the stacked vector as defined in eq. (23) and is shown in fig. 5.

4.2 Validation

In this section, we perform validation tests for the pseudo spectrum estimator, the covariance matrix and the statistical pipeline.

In order to validate the pseudo spectrum estimator and the semi-analytical expression of the covariance matrix given in eq. (23), we generate 1000 sets of four correlated full-sky maps with appropriate auto and cross spectra (for $\{\textsf{LOWZ},\textsf{CMASS},\textsf{QSO},{\kappa_{\mathrm{CMB}}}\}$) computed using eq. (4), using the synfast function of Healpy161616synfast generates independent identically distributed random normal variables and makes linear combinations of these variables to generate Gaussian distributed spherical Fourier coefficients $a_{\ell m}$ with appropriate covariances. following a procedure similar to Bianchini et al. (2015); Nicola et al. (2016). These maps are then masked and their pseudo spectra are compared to the analytical expected value from eq. 11, as shown in fig. 6. In order to validate the covariance matrix, a similar set of full-sky maps are then added realistic noise: for each lensing convergence map, we add an uncorrelated Gaussian noise with spectrum $N^{\kappa}_{\ell}$ given by the approximate spectrum delivered by the Planck Legacy Archive, which is precise enough for the covariance validation. For each galaxy density map, we also need to simulate Poisson sampling. To do so, we generate a map where the value in pixel $p$ is a Poisson random variable of mean $\lambda_{p}$, i.e.,

[TABLE]

where $\delta_{p}$ is the simulated overdensity at pixel $p$ and $\overline{N}$ is the mean number of galaxies per pixel (different for the three samples). A reconstructed density map is then built using eq. (17), which now incorporates Poisson shot-noise. These full-sky maps are then masked and their pseudo spectra are evaluated. The empirical covariance of the sets of pseudo spectra is finally computed and compared to the semi-analytical covariance we use throughout this analysis. The result in fig. 7 shows good agreement and validates the estimator and the simulation of the matrices X and Y (that were computed using generic spectra, see appendix B).

The statistical pipeline is validated by performing a Monte Carlo analysis similar to the one performed in Penna-Lima et al. (2014). Specifically, we want to check if the estimated parameters are unbiased. For that purpose, given the adopted fiducial model, we use the likelihood as the probability distribution of the pseudo spectra $\tilde{C}_{\ell}^{AB}$ to generate random samples (i.e. sets of stacked vectors $\tilde{\mn@boldsymbol{C}}^{\mathrm{obs}}$). For each sample, we fit all parameters to be tested, thus building a collection $\{\mn@boldsymbol{\theta}_{i}\}$ of best-fit values for these parameters. At step $n$, the means $\overline{\mn@boldsymbol{\theta}}_{n}=\sum_{i=1}^{n}\mn@boldsymbol{\theta}_{i}/n$ and standard deviations of the collection of best-fit values are computed. The largest relative error (LRE) over parameter means is computed and we repeat the process, adding more samples, until the LRE has reached a level of 0.1% and check that the fiducial values are within the error bars. For this test, we only use one sample of galaxies with the redshift distribution of the CMASS sample and generate samples of $\tilde{\mn@boldsymbol{C}}_{\ell}^{\kappa_{\mathrm{CMB}}\times\kappa_{\mathrm{CMB}}}$, $\tilde{\mn@boldsymbol{C}}_{\ell}^{\kappa_{\mathrm{CMB}}\times\delta_{\textsf{CMASS}}}$ and $\tilde{\mn@boldsymbol{C}}_{\ell}^{\delta_{\textsf{CMASS}}\times\delta_{\textsf{CMASS}}}$. Results of this test are plotted in figs. 8 and 9, showing respectively the evolution of the mean values of the best fit parameters $\overline{\mn@boldsymbol{\theta}}_{n}$ as a function of the number of realisations $n$ and the distribution of the best-fit parameters for those same realisations. They confirm that the parameters’ estimators are unbiased at least at the 0.1% level.

4.3 Null tests

We present in this section null tests that were performed to exclude potential systematic errors related to the masks and selection of the spectroscopic tracers.

In order to assess potential leakage of power in the cross spectra due to the masking, we cross-correlate the 100 simulated reconstructed lensing maps of the Planck Legacy Archive with the observed density maps of the three galaxy samples, and then correlate the observed lensing map with 100 simulated galaxy maps. This procedure removes cosmological angular correlation, and what correlation remains will be linked to the masks themselves. We find that all results are consistent with no correlation, excluding strong contamination from masking. We also measure the cross spectra between the galaxy and quasars sample and find marginal correlations, well below the auto-correlation signals.

Variable observational conditions during the SDSS photometric survey could potentially result in non-uniform selection functions of the galaxy and quasar samples, and introduce artificial power in the auto spectra at large scales. In order to exclude dramatic power leakage, we constructed maps of resolution $N_{\mathrm{side}}=128$ of the seeing, sky flux, extinction (for the $g$, $r$ and $i$ bands) and air mass of the photometric observations that were used to select galaxies and quasars in the spectroscopic catalogs,171717To do so, we made use of the CasJob service of the SDSS SkyServer, at http://skyserver.sdss.org/casjobs/. as well as a map of the stellar density using the same cuts as in Ross et al. (2017). We then verify that the cross power spectra with the density maps (built in section 3.2.3) are consistent with a null value. To do so, we use the covariance matrix of the pseudo cross spectra given in eq. 23, taking advantage of the fact that the masks are the same, such that the covariance matrix of reduced pseudo cross spectra, defined as

[TABLE]

is approximated by

[TABLE]

where $M_{\ell\ell^{\prime}}^{gg}$ is the mixing matrix associated with the masks of the LOWZ, CMASS and QSO samples, and $\tilde{D}_{\ell}$ denotes the pseudo auto power spectrum including noise. Here we have approximated the ratio of full-sky to partial sky spectra as $\tilde{D}_{\ell}^{AA}\approx f_{\rm sky}^{A}D_{\ell}^{AA}$ and have neglected the variance181818Note that we consider here the covariance over cosmological realizations at fixed realization of the observational systematics. of $\sqrt{\tilde{D}_{\ell}^{\delta_{g}\times\delta_{g}}}$ in the denominator in eq. 28 since it is measured at high signal-to-noise ratio (which, because of cosmic variance, is a ${\mathcal{O}(1/\sqrt{\ell})}$ approximation further tamed by the square root).

Figures 11, 12 and 13 show the measured cross power spectra $\tilde{C}_{\ell}^{\mathrm{syst}\times\delta_{g}}$ for multipoles ${0<\ell\leq 192}$. $\chi^{2}$ statistics are measured (values are given in the legends of the figures) and consistent with no correlation191919The stellar density shows significant anti-correlation with the quasar density, but only for the smallest multipoles which we discard. Above $\ell>20$, it is consistent with zero., with no discernable trend across multipoles. From these tests, we conclude that photometric systematics do not strongly correlate with the overdensity maps.

5 Analyses and results

5.1 Cosmological model

Our base model is the standard $\Lambda$CDM model with flat spatial sections (hence $\Omega_{k}=0$) and a DE component with equation of state $w=-1$. The base parameters are the present-day baryon and cold dark matter densities, $\omega_{\mathrm{b}}\equiv\Omega_{\mathrm{b}}h^{2}$ and $\omega_{\mathrm{c}}\equiv\Omega_{\mathrm{c}}h^{2}$, respectively – where $\Omega_{i}=\rho_{i}/\rho_{c}$ is the ratio of the component’s energy density to the critical energy density $\rho_{c}$–, the Hubble constant today $H_{0}=$100\text{,}\mathrm{\mathit{h}}\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}$$, the redshift of reionization $z_{\mathrm{re}}$, the logarithm of the primordial curvature $\zeta$ dimensionless power spectrum $\ln 10^{10}A_{\mathrm{s}}$ and its tilt $n_{\mathrm{s}}$ such that

[TABLE]

with the pivot scale $k_{0}=$0.05\text{,}{\mathrm{Mpc}}^{-1}. We include massive neutrinos, parametrised by the effective number of neutrinos in the relativistic limit $N_{\rm eff}=$3.046 (taking into account non-instantaneous decoupling), an effective temperature $T_{\nu}/T_{\gamma}=$0.716,11$$, where $T_{\gamma}$ is the photon temperature (slightly departing from $\left(4/11\right)^{1/3}$ to take into account neutrino heating from electron/positron annihilation, see Lesgourgues et al., 2009), and using one massive neutrino of mass {{m_{\nu}}=0.06\text{,}\mathrm{eV}} and two massless neutrinos, consistent with the Planck base $\Lambda$CDM model. The linear matter power spectrum $P_{\mathrm{m}}(k,z)$ and the CMB temperature power spectrum $C_{\ell}^{TT}$ are computed using the Cosmic Linear Anisotropy Solving System (CLASS) as a backend to NumCosmo. The non-linear matter power spectrum is computed using a halofit prescription (Smith et al., 2003) implemented in NumCosmo, with parameters from Takahashi et al. (2012), modified to take into account neutrinos as in CLASS. Reionization is modelled in a CAMB-like fashion (Lewis et al., 2000) and parametrised by the mid-point $z_{\mathrm{re}}$, fixed width $\Delta z_{\mathrm{re}}=0.5$, and includes Helium reionization at a fixed redshift $z_{\mathrm{re}}^{\mathrm{He}}=3.5$. Recombination is computed within CLASS and Big Bang nucleosynthesis is computed with PArthENoPE202020http://parthenope.na.infn.it/(Pisanti et al., 2008).

From the constraints that we will obtain in our analyses, we will also estimate the total matter density parameter212121The DE density parameter today is $\Omega_{\Lambda}\approx 1-\Omega_{\mathrm{m}}$ since we consider only flat space sections (neglecting radiation). ${\Omega_{\mathrm{m}}=\Omega_{\mathrm{b}}+\Omega_{\mathrm{c}}+\Omega_{\nu}}$, the optical depth to the last scattering surface $\tau$ and the variance of the linear matter density fluctuations $\sigma_{8}^{2}$ in spheres of radius $R=$8\text{,}{\mathrm{\mathit{h}}}^{-1}\text{,}\mathrm{Mpc}$$ extrapolated to $z=0$,

[TABLE]

where the top-hat window function is $W(k,R)=3j_{1}(kR)/kR$ and the matter power spectrum is computed from linear theory.

5.2 Statistical analysis

In this section, we describe our Bayesian statistical analysis and present constraints on cosmological parameters and BOSS spectroscopic tracers’ biases.

We first apply the Markov Chain Monte Carlo (MCMC) approach using only CMB lensing-LSS correlations data – the set of the six auto and cross spectra of CMB lensing and BOSS galaxy and quasar overdensities, as shown in fig. 4 – and varying only a subset of cosmological parameters in order to assess the constraining power of these. In particular, we also consider different combinations of the auto and cross spectra to measure the effects on the parameter constraints provided by these probes. Then, we add CMB temperature information and obtain constraints on the $\Lambda$CDM model and extensions including the total mass of neutrinos ${\Sigma m_{\nu}}$ – that impacts small-scale structure formation – and the DE equation of state $w$ – that impacts the expansion in the low redshift Universe.

We perform MCMC analyses using an ensemble sampler222222Some authors refers to ensemble samplers as population Monte Carlo. with many walkers (32 to 1000), moving their positions in the parameter space as an ensemble via a stretch move scheme (Goodman & Weare, 2010) implemented in NumCosmo. We monitored the convergence of the chains using three numerical tools, namely the Multivariate Potential Scale Reduction Factor (MPSRF, Gelman & Rubin, 1992; Brooks & Gelman, 2012), the Heidelberger-Welch test (Heidelberger & Welch, 1981, 1983) and the Effective Sample Size (ESS); see appendix E for more details.

These diagnostics can fail in different situations. For this reason, we also performed three different visual inspections for each parameter:

the parameter trace plot, i.e. the value of the parameter for a given walker vs iteration time. 2. 2.

the ensemble distribution trace plot, that is, the empirical ensemble distribution given by the walkers’ positions vs iteration time. This allows us to monitor the evolution of the ensemble mean and variance. 3. 3.

the total mean vs the cumulative sum of the ensemble means: if the chain has reached convergence, the difference (scaled by the spectral density at null frequency) is distributed as a brownian bridge, the $L_{2}$ norm of which is used in the Schruben test.

For all MCMCs, we ran them until all the relative errors of the means were smaller than $10^{-2}$; at this point, we applied all the tests above, and if the chains failed some of them, we continued the run until all tests were satisfied.

5.2.1 Constraints on $\sigma_{8}$ and $\Omega_{\mathrm{m}}$ from CMB lensing–LSS correlations only

Data from CMB lensing and spectroscopic tracers of matter alone cannot efficiently constrain all cosmological parameters. However, we want to highlight the cosmological information carried by these probes. To do so, we perform several MCMC analyses considering a subset of free cosmological parameters. These are only illustrative in the sense that the posterior distribution of cosmological parameters will be shrunk from fixing some others.

The theoretical spectra have different dependences on the cosmological parameters and galaxy/quasar biases. The most explicit dependencies of the angular spectra $\tilde{C}_{\ell}$ are on the power spectrum amplitude ($A_{\mathrm{s}}$ or $\sigma_{8}$), the matter density parameter $\Omega_{\mathrm{m}}$ and the galaxy/quasar biases (see the kernels in section 2.1):

[TABLE]

This system is closed, i.e. in principle, comparing the various auto and cross spectra should allow for an non-degenerate estimation of the parameters.

Therefore, we run MCMCs freeing $\omega_{\mathrm{c}}$, $\ln 10^{10}A_{\mathrm{s}}$ and the galaxy/quasar biases, and fixing all other cosmological parameters. Their fiducial values are from Planck 2015 “TT,TE,EE+lowP+lensing+ext” best fits (Planck Collaboration et al., 2016e). We assume flat prior distributions over wide ranges (larger than the sampled ranges). In order to distinguish and quantify the information contained in the various measured auto and cross spectra, we try different combinations. That is, we run an MCMC with the full dataset (“$\kappa\kappa$+\kappa$$\delta_{g}+\delta_{g}$$\delta_{g}”), and then repeat without the cross spectra (“$\kappa\kappa$+\delta_{g}$$\delta_{g}”), without the CMB lensing auto spectrum (“\kappa$$\delta_{g}+\delta_{g}$$\delta_{g}”) and without the galaxy auto spectra (“$\kappa\kappa$+\kappa$$\delta_{g}”). We run these chains with 100 walkers to ensure a good mixing. Their MPSRFs are below 1.02 and the correlation lengths are of order 20-40, varying amongst parameters.

The sampled posterior distributions of these parameters are shown in fig. 14 for the full dataset and the three subsets aforementioned. We note that the “\kappa$$\delta_{g}+\delta_{g}$$\delta_{g}” and “$\kappa\kappa$+\kappa$$\delta_{g}” subsets, dominated by respectively galaxy clustering and CMB lensing information, provide complementary information, since the correlations between the parameters, except those in the ($b_{\textsf{LOWZ}}$, $b_{\textsf{CMASS}}$) plane, present different alignments (see the blue and yellow confidence regions in fig. 14). Therefore, the constraints on the parameters are greatly improved when combining both auto spectra, “$\kappa\kappa$+\delta_{g}$$\delta_{g}”. Apart from the constraints on $b_{\textsf{LOWZ}}$, $b_{\textsf{CMASS}}$, which are already strongly determined by galaxy density auto spectra, the additional information contained in the cross spectra narrows the distribution, as can be observed in the (${\ln 10^{10}A_{\mathrm{s}}}$,$\omega_{\mathrm{c}}$) plane by comparing the confidence regions for “$\kappa\kappa$+\delta_{g}$$\delta_{g}” (in green) with “$\kappa\kappa$+\kappa$$\delta_{g}+\delta_{g}$$\delta_{g}” (in red). The addition of the cross spectra decreases the statistical error by 10% for ${\ln 10^{10}A_{\mathrm{s}}}$ and 20% for $\omega_{\mathrm{c}}$, and slightly shifts the best fits (by less than 1 $\sigma$). This plane is translated into the ($\sigma_{8}$, $\Omega_{\mathrm{m}}$) plane in fig. 15 where the degeneracy breaking expected from the joint analysis is highlighted.

5.2.2 Cosmological constraints from the joint analysis of Planck and BOSS data

In this section, we carry out the analysis combining CMB temperature and the joint likelihood of CMB lensing and galaxy/quasar densities used in the previous section (that is all six power spectra of fig. 4) to derive cosmological constraints. First, we constrain the 6-parameters $\Lambda$CDM model (with parameters $\omega_{\mathrm{b}}$, $\omega_{\mathrm{c}}$, $H_{0}$, $z_{\mathrm{re}}$, $\ln 10^{10}A_{\mathrm{s}}$ and $n_{\mathrm{s}}$) and compare constraints from the full joint analysis (“Planck $TT$ + lensing $\otimes$ BOSS tracers”) with that derived from CMB temperature anisotropies only (“Planck $TT$”) or CMB temperature and lensing (“Planck $TT$ + lensing”). For the joint analysis, we will also consider two cuts as mentioned in section 3.3. Note that we neglect the correlation between CMB temperature and the matter density at later times (either baryonic matter in galaxies and quasars or dark matter weighted by CMB lensing), i.e., we neglect the late ISW effect, as it is not yet detected with a strong significance (Planck Collaboration et al., 2016g; Nicola et al., 2016), and we discuss possible consequences in section 5.3. In practice, this means that we approximate the total likelihood by the product of the CMB temperature and CMB lensing-LSS correlations likelihood functions.

We use the Planck likelihood codes Plik and Commander (Planck Collaboration et al., 2016d) respectively for high and low multipoles of the temperature-only power spectrum $C_{\ell}^{TT}$. The likelihood codes introduce 15 additional nuisance parameters related to foreground and instrument models ($A^{\mathrm{CIB}}_{217}$, $\xi^{\mathrm{tSZ}\times\mathrm{CIB}}$, $A^{\mathrm{tSZ}}$, $A^{\mathrm{PS}}_{100}$, $A^{\mathrm{PS}}_{143}$, $A^{\mathrm{PS}}_{143\times 217}$, $A^{\mathrm{PS}}_{217}$, $A^{\mathrm{kSZ}}$, $A^{\mathrm{dust}TT}_{100}$, $A^{\mathrm{dust}TT}_{143}$, $A^{\mathrm{dust}TT}_{143\times 217}$, $A^{\mathrm{dust}TT}_{217}$, $c_{100}$, $c_{217}$ and $y_{\mathrm{cal}}$; see Planck Collaboration et al., 2016d). We use the profile likelihood to speed up our MCMC analyses, subfitting the nuisance parameters for each set of cosmological parameters. We describe this methodology in appendix C, and also show that it does not affect the results on cosmological parameters.

We use flat, large priors on all cosmological parameters and tracer biases, common to all MCMC runs. The lower and upper limits are given in table 1, together with constraints derived from the (optimistic) joint analysis.

Contraints on $\Lambda$CDM.

Figure 16 shows the constraints on the base $\Lambda$CDM model’s parameters for the three aforementioned data combinations. When using CMB $TT$ only, we find parameter constraints that are in perfect agreement with the Planck analysis232323The detailed results of the Planck MCMC analyses are available here: https://wiki.cosmos.esa.int/planckpla2015/images/f/f7/Baseline_params_table_2015_limit68.pdf.. The strong degeneracy observed between the power spectrum amplitude $A_{\mathrm{s}}$ and the reionization redshift $z_{\mathrm{re}}$ corresponds to the amplitude of the power spectrum of CMB temperature anisotropies, which is proportional to $A_{\mathrm{s}}e^{-2\tau}$ where $\tau$ is the optical depth to the last scattering surface, strongly dependent on reionization history. Adding CMB lensing drives $A_{\mathrm{s}}$ and $z_{\mathrm{re}}$ towards lower values along this degeneracy with a shift of about $1\sigma$ for each parameter as evinced by the one- and two-dimensional projection of the posterior distributions. Finally, adding information of LSS tracers (both the auto-correlations and cross-correlations with CMB lensing) provides only slightly smaller contours for these parameters and, therefore, they do not significantly help in breaking this degeneracy. We find no significant improvement for $\tau$ or $\sigma_{8}$, although it is consistent with the constraints from CMB lensing. In the ($\sigma_{8},\;\Omega_{\mathrm{m}}$) plane (see fig. 17), we observe that early Universe data favors bigger values of $\sigma_{8}$ than the late one, as repeatedly reported in the literature (Hildebrandt et al., 2017; Planck Collaboration et al., 2016i). Whether this is indication of new physics or a systematics artefact is beyond the scope of this work, but it might be an important issue in the future.

Nonetheless, there is an improvement of order 20% on the measurements of $H_{0}$ and $\omega_{\mathrm{c}}$ for both the optimistic and conservative cuts, and the total volume in parameter space is notably reduced. This can be quantified in the approximation of Gaussian posteriors by computing a figure of merit ${\mathrm{FoM}\equiv\absolutevalue{\mathrm{Cov}(\Theta_{\rm cosmo})}^{-1/2}}$, where $\mathrm{Cov}(\Theta_{\rm cosmo})$ is the empirical covariance matrix of cosmological parameters. Relative to the CMB $TT$ only posterior, adding CMB lensing increases the FoM by a factor 2.0, and the full joint analysis increases the FoM by a factor $2.9$. We find a best-fit of {H_{0}=68.8\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}}, slightly higher than CMB temperature alone – albeit still lower than distance measurements from supernovæ (Riess et al., 2016) or time delays in strong lensing (Bonvin et al., 2017) – and ${\omega_{\mathrm{c}}=0.117\pm 0.002}$ as the degeneracy between $A_{\mathrm{s}}$ and $\omega_{\mathrm{c}}$ is broken by the lensing–LSS correlations (see fig. 14). This results in a constraint on the matter density parameter ${\Omega_{\mathrm{m}}=0.296\pm 0.011}$. Additionally, we obtain strong constraints on the biases of the galaxy samples, respectively

[TABLE]

These 4% constraints are in general agreement with previous measurements using angular power spectra (Ho et al., 2012)), which have the advantage of being model-independent in the sense that estimating angular power spectra $\tilde{C}_{\ell}$ does not require any assumption on cosmology since we don’t measure distances. Moreover, all the cosmological parameters of the $\Lambda$CDM model are fitted alongside. Note, however, that our modelling assumes a constant bias, that can be interpreted as a redshift- and scale-averaged bias, when other analyses exploring the non-linear regime used a scale-dependent bias (in the form of a Taylor expansion, e.g. in Gil-Marín et al. (2017)) or simply more redshift bins. Interestingly, the analysis also shows significant correlations between the biases and cosmological parameters, in particular with $\omega_{\mathrm{c}}$, $H_{0}$ and $A_{\mathrm{s}}$. If one considers biases as effective parameters encoding structure formation and clustering of galaxies, these correlations can shed light on the astrophysical and cosmological processes governing the formation of such structures. Finally, we also obtain a broad constraint on the bias of the uniform sample of quasars from the cross-correlation with CMB lensing

[TABLE]

This value is in tension with other measurements (White et al., 2012; DiPompeo et al., 2015; Laurent et al., 2016, 2017) that found a bias of order 3 to 4 (although they did assume a cosmology), but in agreement with Alonso et al. (2017) that also uses the cross angular power spectrum. We found no difference when fitting for this bias when using data from only the northern or southern galactic caps, excluding a possible strong asymmetry and contamination of higher multipoles. We note a surprising trough in the $\tilde{C}_{\ell}^{\kappa_{\mathrm{CMB}}\times\mathrm{QSO}}$ cross spectrum (see fig. 4) around $\ell\sim 400$ that we could not explain. However, this bias directly depends on the amplitude of lensing that may be underestimated (Planck Collaboration et al., 2014).

Finally, we observe that both cuts yield very similar results for cosmological parameters, with little improvement from the optimistic cut for the $\Lambda$CDM model (which can arise from small tensions between data sets). The conservative cut yields lower clustering biases for LOWZ and CMASS galaxies, while remaining within $1\sigma$ of the optimistic cut (that includes more non-linear scales).

Constraints on the total mass of neutrinos ${\Sigma m_{\nu}}$ and the dark energy equation of state $w$.

In the next set of MCMCs, we additionally sample the total mass of neutrinos ${\Sigma m_{\nu}}$ (with one massive and two massless neutrinos) or the DE equation of state $w$ (where $w$ is constant over time) separately and compare the performance of the joint analysis in these extended models.

In the first case, we find that the joint analysis yields an upper bound on the total mass of neutrinos of {{\Sigma m_{\nu}}<0.28\text{,}\mathrm{eV}\ [68\%]} ($<$0.39\text{,}\mathrm{eV}$\ [95\%]$) with the optimistic cut, dividing the higher bound by a factor of two with respect to the constraint from CMB $TT$ alone (see fig. 18). We do not detect a total neutrino mass significantly different from zero, but the best fit we obtain around $m\sim$0.15\text{,}\mathrm{eV}$$ is in agreement with lower bounds around $0.05\text{\,}\mathrm{eV}$ derived from neutrino oscillations (Olive, 2014), and in agreement with cosmological upper bounds around $0.12\text{\,}\mathrm{eV}$, e.g. that derived from the combination of CMB, either with the Lyman-$\alpha$ forest power spectrum (Palanque-Delabrouille et al., 2015) or with BAO measurement (Vagnozzi et al., 2017). However, the joint analysis with the conservative cut does not improve the upper bound when compared to the “Planck $TT$ + lensing” case. As shown in fig. 19, the galaxy auto power spectra are sensitive to the total mass of neutrinos as they can probe relatively small scales at low redshift, where massive neutrinos tend to smooth out density fluctuations. Therefore, the improved constraints with the optimistic cut is likely to be due to small scale contributions. Moreover, this also means that the total mass of neutrinos should be positively correlated with galaxy biases, which is indeed observed in the lower panels of fig. 18. We also observe that adding LSS information significantly improves the constraints on the other cosmological parameters in this extended model, for both the conservative and optimistic cuts. Because of the anti-correlation between $H_{0}$ and the total mass of neutrinos, the joint analysis favours a higher expansion rate {H_{0}=67.6_{-1.4}^{+1.3}\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}} (optimistic case) than CMB data alone. It also noticeably shifts the posterior distributions for $z_{\mathrm{re}}$ and $A_{\mathrm{s}}$ towards lower values, resulting in a lower value of the reionization optical depth ${\tau=0.088\pm 0.020}$ (though still higher than CMB polarization).

In the second case, we release $w$, the sum of the neutrino masses being fixed to {{\Sigma m_{\nu}}=0.06\text{,}\mathrm{eV}}. CMB temperature anisotropies are only very weakly sensitive to DE and CMB lensing probes the Universe at redshift $z\sim 2$ where matter is still dominating. Therefore these probes do not contain much information on $w$. We observe a strong anti-correlation between $w$ and the Hubble parameter $H_{0}$, meaning that observations can be matched by a more slowly expanding Universe with a more negative DE pressure. As a consequence, CMB only posteriors hit the upper bound of the prior set at $H_{0}=$100\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}$$, thus artificially shrinking cosmological constraints. Adding LSS information becomes necessary and rewarding as it breaks the degeneracies of the constraints on $H_{0}$, $\ln 10^{10}A_{s}$, $z_{\mathrm{re}}$ and $w$. Constraints from the joint analysis (${w=-1.04^{+0.21}_{-0.32}}$) are consistent with a cosmological constant (${w=-1}$), while constraints from CMB favour a lower value of $w$. We also note a strong correlation between the biases and $w$ of 79% and 92% for CMASS and LOWZ respectively.

In summary, in both cases, constraints from the joint analysis are substantially better for almost all parameters because of its ability to break degeneracies related to the chosen new parameters. This result constitutes a forceful encouragement to perform this type of analysis when data from the next generation of surveys becomes available.

Constraints on $w\mathrm{CDM}+{\Sigma m_{\nu}}$.

Finally, in the last set of MCMC analyses, we release both the total mass of neutrinos ${\Sigma m_{\nu}}$ and the DE equation of state $w$ and demonstrate that a joint analysis of currently available data can constrain the 8-parameter cosmological model $w$CDM + ${\Sigma m_{\nu}}$, with free parameters $H_{0}$, $\omega_{\mathrm{b}}$, $\omega_{\mathrm{c}}$, $\ln 10^{10}A_{\mathrm{s}}$, $n_{\mathrm{s}}$, $z_{\mathrm{re}}$, ${\Sigma m_{\nu}}$ and $w$. Similarly to the previous cases, the results are presented on fig. 20 for the full joint analysis and for CMB data, allowing for comparison. As expected, CMB temperature data alone is found unable to constrain this model, and posterior distributions hit the prior upper bound on $H_{0}$. However, additional information extracted by the joint analysis (partially) breaks the $w-H_{0}$ degeneracy, enabling for control of all eight cosmological parameters, plus the biases, i.e. eleven parameters in total. For these reasons and for readability of the figures, we only show the posterior distribution from “Planck $TT$ + lensing” data and the joint analysis.

Both the conservative and optimistic cuts in the joint analysis yield constraints in agreement with the current picture of the $\Lambda$CDM model as well as those obtained in the previous sections, with a value of ${w=-0.93^{+0.16}_{-0.22}}$ consistent with a cosmological constant. The conservative cut does not improve the upper bound on the total neutrino mass, but it substantially improves constraints on $H_{0}$ and $w$ in comparison to those derived from CMB data alone. With the optimistic cut, we obtain a higher bound on the mass of the neutrinos of {{\Sigma m_{\nu}}<0.28\text{,}\mathrm{eV}\ [68\%]} ({<0.39\text{,}\mathrm{eV}\ [95\%]}) and a low value of the Hubble constant of {H_{0}=65.6^{+6.7}_{-4.8}\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}}, albeit with larger error bars with respect to the previous sections. The correlation coefficient matrix reveals a strong correlation of galaxy biases with $H_{0}$, $w$ and ${\Sigma m_{\nu}}$ (see fig. 21) indicating that upcoming surveys will require exquisite control of these biases to get tight constraints on $w$ and its possible time evolution.

Because of the degeneracy between $H_{0}$, $w$ and ${\Sigma m_{\nu}}$, precision is lost on $\Omega_{\mathrm{m}}$, even though the physical density ${\omega_{\mathrm{m}}\equiv\Omega_{\mathrm{m}}h^{2}}$ is well constrained by CMB $TT$ and CMB lensing even in this model (we find $\omega_{\mathrm{m}}=0.1411\pm 0.0023$). In the ($\sigma_{8},\;\Omega_{\mathrm{m}}$) plane (see fig. 22), we obtain constraints that are consistent from the joint analysis over the models tested here, with increasing degeneracy. We measure ${\sigma_{8}^{2.7}\Omega_{\mathrm{m}}=0.1709\pm 0.0068}$ from the joint analysis on the 8-parameter $w$CDM + ${\Sigma m_{\nu}}$ model.

5.3 Limits and perspectives

In this section, we discuss assumptions that were made and technical difficulties that we were able to pinpoint.

We neglected the correlation, generated by the ISW effect, between the CMB temperature map and the large-scale structure as traced by CMB lensing or spectroscopic tracers. This correlation originates in the net energy gain (loss) of photons crossing gravitational potentials wells (hills) evolving thanks to dark energy. In principle, this would lead to underestimation of error bars on cosmological parameters. However, this correlation is weak and affects only very large scales ${\ell\lesssim 40}$, and it has not been detected with a strong statistical significance on SDSS galaxies: the signal-to-noise ratios for the correlation with the LOWZ and CMASS samples reported by the Planck collaboration is of order 2.4, and that with the lensing map (corresponding to a temperature bispectrum) is of order 3.2 (Planck Collaboration et al., 2016g). Therefore, taking this cross-correlation term into account would not dramatically change our constraints.

One possible source of systematics in the galaxy-lensing cross-correlations comes from the tSZ component separation that is required to produce the lensing map (Planck Collaboration et al., 2016b, c, h). Free electrons in hot galaxy clusters imprint a specific local spectral distortion on the CMB temperature map. These clusters must be identified and removed before measuring the spatial distortion due to gravitational lensing. If these clusters hold some of the galaxies in the samples we use, this might lead to a systematic underestimation of the lensing signal in the direction of these galaxies. However, the SZ decrement from SDSS LRGs is small, as can be seen in Table 2 of Hand et al. (2011). Moreover, the residual SZ signal primarily increases the noise in the lensing map and is unlikely to produce appreciable bias (see the systematics checks in Madhavacheril et al., 2015).

We now discuss theoretical uncertainties. The first one comes from the Limber approximation (LoVerde & Afshordi, 2008): it fails at very large scales, the transition scale depending on the width of the redshift bin considered (Campagne et al., 2017). In this work, we used very broad redshift bins and discarded low multipoles $\ell<20$, so as to be in the safe regime of the approximation. We also did not consider redshift space distortion (RSD) for the same reasons, as they were shown to be negligible for $\ell\gtrsim 20$ (Alonso et al., 2015; Padmanabhan et al., 2007) and rapidly decreasing with the width of redshift bins (Saito, 2016). To test that RSD or other large-scale effects were not driving parameter constraints inconsistently, we performed an MCMC analysis of CMB lensing-galaxy correlations on a smaller multipole range (${100\leq\ell<400}$, see fig. 23) and found no significant deviation. However, future surveys aiming at measuring extremely large-scales (Alonso & Ferreira, 2015) will require better modelling, especially for tomographic studies with thin redshift bins.

In this paper, we used a Gaussian likelihood and a Gaussian covariance, i.e. we did not incorporate higher order statistics of the matter density field nor the so-called super-sample variance due to the finite size of the surveyed volume and inaccessible modes therein (Schaan et al., 2014; Krause & Eifler, 2017; Schaan et al., 2017). At the current level of signal-to-noise ratio, these simplifications are safe242424Super-sample variance, for instance, can introduce correlations between multipoles of order 10% at $z=0.1$ and decreasing with redshift and multipole number, see Lacasa & Rosenfeld (2016). but they should be lifted in future data analysis. One limitation of our method regarding the covariance matrix is that the computation of the X and Y matrices, even if it needs to be done only once, is numerically expensive since they grow linearly with the multipole range but as $n^{4}/4$ with the number $n$ of different masks, and the estimate remains noisy far from the diagonal. However, this method has the advantage of naturally taking care of partial sky coverage, without the need of inverting the mixing matrix to recover full-sky spectra, necessarily introducing numerical noise in the data. Devising a method that takes care of partial sky coverage whilst incorporating all relevant non-Gaussian terms will be an important task for future surveys (Lacasa et al., 2016).

Finally, the non-linear power spectrum of the matter density field suffers from theoretical uncertainties (Baldauf et al., 2016). Throughout this analysis, we have used a version of the halofit model that includes the effect of massive neutrinos and predicts the matter power spectrum past the non-linear transition scale (around $k_{\mathrm{nl}}\sim$0.1\text{,}\mathrm{\mathit{h}}\text{,}{\mathrm{Mpc}}^{-1}$$). So as to be agnostic, we considered two cuts at small scales at $0.1\text{\,}{\mathrm{Mpc}}^{-1}$ and $0.15\text{\,}{\mathrm{Mpc}}^{-1}$. The halofit model reaches percent agreement with simulations at these scales and we have demonstrated that we recover consistent constraints for both cuts, with small shifts of best-fit values for all parameters, though expectedly tighter in the optimistic case. The scientific gain from near future surveys in terms of cosmological constraints will depend on our ability to model these non-linearities. In particular, the suppression of power due to massive neutrinos (Lesgourgues & Pastor, 2006) and the contribution of baryonic and feedback processes (Leauthaud et al., 2017) at these scales will certainly be an important theoretical issue for future surveys.

6 Conclusion

Cosmological experiments carried out in the last few decades have enabled the construction of the $\Lambda$CDM model. In this picture, cold dark matter drives the formation of the large-scale structure of the Universe and dark energy fuels the recent accelerated expansion. The combination of independent observations, such as the map of the anisotropies of the cosmic microwave background, distances of type IA supernovæ and the measurement of the scale of the baryon acoustic oscillations, have set constraints on the content of the Universe. However, the analysis of currently available data cannot distinguish between various models of dark matter and dark energy. Going further and deciphering the nature of these components requires better constraints, and thus, more information. To this end, deep galaxy surveys – such as LSST, Euclid and WFIRST – and CMB imagers – such as CMB-S4 and the Simons Observatory – with wide sky coverage and high resolution are currently under development. In the coming decade, they will probe the matter density field with ground-breaking precision and significantly increase the amount of cosmological information. Independent cosmological analyses have a strong potential to reveal new science, but model comparison will rely on exhausting the cosmological information held in the measurements of different cosmic probes and all their cross-correlations. In other words, joint analyses of these probes are required, thus explaining the recent intense activity around this subject (see the introduction).

In this paper, we have presented a joint analysis of currently available data combining CMB measurements –both temperature anisotropies and gravitational lensing– from the Planck satellite and LSS tracers from the SDSS-III/BOSS spectroscopic survey, taking advantage of the large areas covered by these surveys and their large overlap (a requirement for measuring cross-correlations with a high signal-to-noise ratio). To this end, we developed a general framework in NumCosmo to compute and analyse the auto and cross-correlations between an arbitrary number of cosmological probes, which is publicly available. In particular, we applied our framework to analyse CMB lensing and galaxy clustering at once by measuring all relevant auto and cross angular power spectra (shown in fig. 4). We validated the likelihood and pseudo-spectra estimators in sections 4.2 and 4.3. Note that our approach required few simplifications on the CMB lensing auto power spectrum (as explained in section 3.1), but we demonstrate in appendix A that our pipeline yields unchanged cosmological constraints, except for a small, yet non-negligible, $0.5\sigma$ shift in the $A_{\mathrm{s}}-z_{\mathrm{re}}$ degeneracy, below the level of statistical errors nonetheless. In section 5.2.1, we showed how including cross-correlation information – already present in the data – improves constraints on cosmological parameters and decreases the statistical errors, for example, by 10% for ${\ln 10^{10}A_{\mathrm{s}}}$ and 20% for $\omega_{\mathrm{c}}$ (when other parameters are fixed). This highlights the fact that ignoring part of the cosmological information (in this case, the cross-correlations) could lead to inaccurate posterior distributions of the parameters.

Next, we included CMB temperature anisotropies information by adding the likelihood of the $C_{\ell}^{TT}$ power spectrum (thus neglecting the small ISW-induced CMB-LSS correlation as discussed in section 5.3) and carried out different MCMC analyses to constrain the base, 6-parameter, flat $\Lambda$CDM model. Finally, we explored constraints on the total mass of neutrinos and the DE equation of state, constraining four different cosmological models (see figs. 16, 18 and 20). We compared the performance of the joint analysis (using two different cuts) with analyses using only CMB data. As expected, constraints from the joint analysis are stronger than those obtained from CMB data only, in all cases. Because of the sensitivity of galaxy clustering and the CMB lensing-galaxy cross-correlations to ${\Sigma m_{\nu}}$ and $w$, we were able to study extended models and constrain up to eight cosmological parameters at once (that is, $H_{0}$, $\omega_{\mathrm{b}}$, $\omega_{\mathrm{c}}$, $A_{\mathrm{s}}$, $n_{\mathrm{s}}$, $z_{\mathrm{re}}$, ${\Sigma m_{\nu}}$ and $w$), which is impossible with either of the data sets considered separately.

As a result, we observe the (partial) breaking of several degeneracies and significantly better constraints for various parameters, although this depends upon exactly which parameters are constrained and which are assumed to be fixed. We thus obtained upper limits on the total mass of neutrinos of $0.28\text{\,}\mathrm{eV}$ [68%] as a result of its impact on galaxy clustering at small scales, which is similar to limits obtained with other comparable analyses. It is interesting to note here that combining CMB and BAO distance measurements currently yields tighter constraints on the total neutrino mass because BAO reconstruction includes more $k$-modes than power spectrum measurements that exclude non-linear scales (Vagnozzi et al., 2017; Cuesta et al., 2016). Better modelling of the power spectrum and clustering biases at these scales is thus a key to improve cosmological constraints. In addition, galaxy clustering information within the joint analysis enabled us to obtain constraints on the dark energy equation of state $w$, found to be consistent with $w=-1$. We have also identified strong correlations between clustering biases and some cosmological parameters, in particular $H_{0}$ and $w$. A downside is that future surveys will have to measure and marginalise over these biases with great precision in order to pin down the values of these parameters and to constrain a possible time dependence of dark energy. Interestingly, if we used a value of {H_{0}=72\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{Mpc}}^{-1}} consistent with distance measurements from type Ia supernovæ, then our constraints would favour a value of the DE equation of state of ${w\sim-1.1}$, i.e. a phantom dark energy, which is disfavoured by theoretical considerations.

The strength of the cosmological constraints derived here is naturally limited by theoretical uncertainties in the non-linear regime on the matter power spectrum and linear clustering biases, as discussed in section 5.3. To address this issue, we have tested two different cut-off scales and shown that it results in consistent cosmological constraints across all models tested, although small shifts are expectedly observed. For instance, the higher bound on the total mass of neutrinos decreases when including smaller scales that are more strongly impacted by neutrinos. Additionally, we have discussed in section 5.3 several other limitations, such as the impact of the ISW effect, potential contamination of cross power spectra with CMB lensing by the SZ effect and the impact of the Limber approximation and redshift-space distorsions.

In this work, we followed an approach similar to Nicola et al. (2016, 2017) that, however, used photometric data from the SDSS and DES, combined with geometric probes. Our approach is based on pseudo spectra with different masks for each observables and therefore exploits the full observed area for each probe, thus maximizing the signal-to-noise ratio of power spectra. We used a semi-analytic, cosmology-dependent covariance matrix that is less noisy than Monte Carlo estimates, but required the assumption that the fields we measure are Gaussian distributed. A unified framework that incorporates non-Gaussian terms in the covariance and handles partial sky coverage remains to be derived.

We focused on spectroscopic observations of galaxies and quasars, insulating us from uncertainties inherent to photometric redshifts. However, we did not use galaxy weak lensing measurements (as was done in Nicola et al. (2016)), a powerful probe of dark energy that will be measured by future deep surveys like LSST, Euclid and WFIRST. The trade-off between the precision of photometric redshifts and the much larger number of galaxies combined to this additional probe will certainly lead to even better results. Combining CMB lensing, galaxy lensing and galaxy clustering (both photometric and spectroscopic) in a fully joint analysis is a promising avenue for cosmological parameters estimation.

Finally, in this near-future scenario of large amounts of data and joint analyses, we will be able to study different cosmological models emerging from different theories of gravity, such as effective field theories of dark energy (Gleyzes et al., 2016) or non-local gravity (Dirian et al., 2016), and hopefully start to distinguish and rule out some models with strong statistical significance.

Acknowledgements

CD would like to thank David Spergel and Emmanuel Schaan for useful comments and the opportunity to present a preliminary version of this work, Dhiraj Hazra for useful comments on the manuscript, Antony Lewis for useful comments about the lensing power spectrum, Alexie Leauthaud, Martin Kunz, Fabien Lacasa, Michele Maggiore, Mathew Madhavacheril and Sunny Vagnozzi for interesting discussions and comments.

The authors thank the anonymous referee for their thorough review, which resulted in major improvements, making our conclusions and cosmological constraints more robust.

The authors acknowledge use of Healpy (Górski et al., 2005), CLASS (Lesgourgues & Tram, 2014) and the Planck likelihood codes Plik, Commander and LensLike available on the Planck Legacy Archive (http://pla.esac.esa.int/pla/).

MPL received support from the National Council for Scientific and Technological Development - Brazil (CNPq grant 202131/2014-9 and PCI/MCTIC/CBPF program), and from the Labex ENIGMASS.

SDPV acknowledges the financial support from CNPq (PCI/MCTIC/CBPF program) and BELSPO non-EU postdoctoral fellowship.

Funding for SDSS-III 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 SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III 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.

Based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.

Appendix A Lensing auto-power spectrum

Throughout this work, we have used the CMB lensing auto-power spectrum from Planck based on the released reconstructed convergence map, neglecting corrections that were applied in the Planck cosmological analysis but that are not yet available in a user-friendly format. Our method is based on pseudo spectra while the Planck collaboration used full-sky spectra, making it difficult to quantify the impact of each correction. Therefore, we compare cosmological constraints derived from CMB temperature and lensing, using either the Planck lensing likelihood code or our code based on NumCosmo using the conservative cut, $40\leq\ell<400$, and the optimistic cut, $20\leq\ell<500$, used throughout the rest of the analysis. We use the ensemble MCMC sampler of NumCosmo and show constraints on parameters of the base $\Lambda$CDM model in fig. 24. We find compatible constraints for the two codes, with very close maximum likelihood values, statistical uncertainties and shifts with respect to CMB $TT$ alone, except for a $0.5\,\sigma$ shift along the $A_{\mathrm{s}}e^{-2\tau}$ degeneracy with an impact on $A_{\mathrm{s}}$ and $z_{\mathrm{re}}$ (and thus $\sigma_{8}$ and $\tau$). This demonstrates that our likelihood code performs very well, but that reconstruction corrections applied to the estimated lensing power spectrum have a non-negligible effect, that is however smaller than statistical errors. We conclude that our constraints on $A_{\mathrm{s}}$ and $z_{\mathrm{re}}$ might be slightly biased, though at a reasonable level that does not impact our main conclusions.

Appendix B X/Y matrices in the covariance

The X and Y matrices appearing in eq. (23) have the following analytical expressions (Brown et al., 2005),

[TABLE]

where the $W^{\ A}_{\ell\ell^{\prime}mm^{\prime}}$ describe the convolution of the mask ($\overline{W}^{\ A}_{\ell\ell^{\prime}mm^{\prime}}$ is its complex conjugate), i.e. if the field $A(\hat{\mathbf{n}})$ has full-sky spherical harmonics coefficients $A_{\ell m}$ and pseudo-coefficients $\tilde{A}_{\ell m}$ then

[TABLE]

These cannot be analytically computed and MC simulations are therefore required. We take advantage of the fact that eq. (23) is exact if initial full-sky spectra do not depend on $\ell$ and that they need not have physically relevant values. The algorithm then proceeds as follows. First, we generate sets of four correlated maps with generic constant input auto and cross spectra, which we mask by the four masks used in our analysis. We then compute the spectra of the masked maps, thus building a collection of estimated pseudo spectra $\{\tilde{C}_{\ell}^{AB,i}\}_{i}$ where $i$ represents the simulation index. The empirical covariance of the set of pseudo spectra is finally computed. Knowing the input spectra, an estimate of $\textbf{{X}}_{\ell\ell^{\prime}}^{ABCD}$ and $\textbf{{Y}}_{\ell\ell^{\prime}}^{ABCD}$ can be obtained using eq. (23). In the case where ${A=B}$ or ${C=D}$ and only in this case, the terms in the square roots in eq. (23) are equal and $\textbf{{X}}_{\ell\ell^{\prime}}^{ABCD}$ and $\textbf{{Y}}_{\ell\ell^{\prime}}^{ABCD}$ cannot be distinguished, but for all the other cases, it requires two sets of simulations to disentangle them.

We estimate the error on the empirical covariance matrices by bootstrapping the pseudo spectra $\{\tilde{C}_{\ell}^{AB,i}\}_{i}$ and require that the ratio of the norms of the error matrix to that of the empirical covariance matrix is smaller than 1%, which in our analysis necessitated more than 200 000 simulations.

In the case where ${A=B=C=D}$, these matrices reduce to symmetrized mixing matrices

[TABLE]

which allows for comparison and validation of the MC simulations (see fig. 25). We find percent-level agreement on the diagonal, with a decreasing precision when moving further away from the diagonal as numerical noise (at least four orders of magnitude smaller than the diagonal elements) starts dominating.

Appendix C Profile likelihood

In order to accelerate our MCMC analyses, we choose to use the profile likelihood instead of the marginal likelihood for the nuisance parameters. The reason is that this procedure decreases the dimension of the parameter space and requires less calls to the Boltzmann code, resulting in an overall faster convergence of the posterior distribution of the cosmological parameters. In practice, it amounts to compute the maximum likelihood estimator value of the nuisance parameters $\hat{\mn@boldsymbol{A}}(\mn@boldsymbol{\theta})$ for each set of cosmological parameters $\mn@boldsymbol{\theta}$ given the data (which is fast), and use this value in the likelihood. The posterior distribution is then given by $\mathcal{L}_{\mathrm{profile}}\left(\mn@boldsymbol{\theta}\middle|C_{\ell}^{TT}\right)\propto\mathcal{L}\left(C_{\ell}^{TT}\middle|\mn@boldsymbol{\theta},\hat{\mn@boldsymbol{A}}(\mn@boldsymbol{\theta})\right)$ while the marginal likelihood is $\mathcal{L}_{\mathrm{marginal}}\left(\mn@boldsymbol{\theta}\middle|C_{\ell}^{TT}\right)\propto\int\mathcal{L}\left(C_{\ell}^{TT}\middle|\mn@boldsymbol{\theta},\mn@boldsymbol{A}\right)\mathrm{d}\mn@boldsymbol{A}$. We demonstrate that it doesn’t affect the results on the cosmological parameters by running two MCMC using only the CMB temperature power spectrum $C_{\ell}^{TT}$, one performing the nuisance parameters subfitting procedure and the other using the standard marginalisation procedure. fig. 26 shows the posterior distribution in these two cases. The mean value of each parameter in both runs is shown in the one-dimensional plots on the diagonal. In all cases, the variation of the mean is much smaller than the statistical variance, and the standard deviation is at worst decreased by 15% in the profile likelihood (for $n_{\mathrm{s}}$ and $\omega_{\mathrm{b}}$, two parameters which are poorly constrained by the other observations), with almost no difference for the other parameters. This indicates that we can use either likelihoods indifferently. Since the profile likelihood method is faster overall, and that we don’t have other nuisance parameters, we used it for all simulations in section 5.2.2.

Appendix D The NumCosmo library

In this section, we provide a short description of the Numerical Cosmology library (NumCosmo, available on GitHub252525https://numcosmo.github.io/, Dias Pinto Vitenti & Penna-Lima (2014)). Apart from the crude data (observational maps), for which we used some Healpix functions to generate the observed pseudo-$C_{\ell}^{AB}$ values, as mentioned in section 3, all other pieces of the pipeline made use of NumCosmo. For a complete description we refer the reader to Vitenti et al. (2018).

NumCosmo contains a comprehensive set of tools to compute cosmological observables and to perform statistical analysis. The library is written in C, but since it uses the GObject framework262626https://developer.gnome.org/gobject/stable/, it is developed in a object-oriented fashion. Additionally, it has automatic bindings for every language that supports GObject introspection (e.g. Python, Ruby or Perl).

Physical models are implemented via the abstract class NcmModel. In particular, the $\Lambda$CDM and $w$CDM models, and all respective relevant functions are implemented in NcHICosmoDE and child classes (such as NcHICosmoDEXcdm), the primordial power spectrum is implemented in NcHIPrim, the reionization model in NcHIReion. Data objects deriving from the abstract class NcmData encapsulate the observations and implement likelihood functions. A general object for statistical analysis NcmFit is then built from the data and the model.

We first address the computation of the theoretical angular power spectrum, $C_{\ell}^{AB}$ (see eqs. (2) and (3)) and the likelihood function eq. (25):

NcXCor: abstract class that comprises, among others, the methods to compute the auto and cross power spectra $C_{\ell}^{AA}$ and $C_{\ell}^{AB}$. 2. 2.

NcXCorLimberKernel: abstract class of the type NcmModel272727Being a NcmModel, each implementation of NcXCorLimberKernel can define a respective set of parameters. For instance, the linear bias, $b(z)$, in eq. (8). which defines the methods and general properties that a kernel $W^{A}(z)$ must implement, for any observable $A$. For instance, the computation of $W^{A}(z)$ at a given $z$ and for a set of cosmological parameters, and the number of multipoles to be calculated.

• •

NcXCorLimberKernelCMBLensing: implements the CMB lensing kernel, $W^{\kappa_{\mathrm{CMB}}}(z)$ [eq. (6)].

• •

NcXCorLimberKernelGal: implements the galaxy kernel, $W^{g}(z)$ [eq. (8)]. 3. 3.

NcPowSpecMNL: abstract class for the nonlinear matter power spectrum. Here we use the halofit approach (Smith et al., 2003; Takahashi et al., 2012), which we implemented in the NcPowspecMNLHaloFit class. The linear matter power spectrum is calculated using NcPowspecMLCBE, i.e., the NumCosmo backend for the Cosmic Linear Anisotropy Solving System (CLASS) (Lesgourgues & Tram, 2011). 4. 4.

NcDataXCor: this object builds the likelihood given by eq. (25), and it derives from NcmDataGaussCov, i.e., the object that describes Gaussian-distributed data with non-diagonal covariance matrix.

Regarding the statistical analyses performed in this work, we made use of the following NumCosmo tools:

NcmFit: implements various interfaces with best-fit finders. In this work, we used the interface with the NLOpt library282828http://ab-initio.mit.edu/nlopt and the Nelder-Mead algorithm. 2. 2.

NcmFitMC: implements the Monte Carlo analysis described in section 4.2, using the same best-fit finder. 3. 3.

NcmFitESMCMC: implements the Ensemble Sampler Markov Chain Monte Carlo (Goodman & Weare, 2010) analysis used throughout section 5 of this paper. It requires an initial sampler NcmMSetTransKern and another sampler to move the walkers NcmFitESMCMCWalker.

Finally, we used CMB temperature data from Planck:

NcPlanckFICorTT: implements Planck foreground and instrumental models for TT measurements. 2. 2.

NcDataPlanckLKL: implements the interface with Planck’s likelihood codes Plik and Commander.

Appendix E MCMC convergence tests

In this work, we checked the convergence of the MCMC chains using three different methods, which we implemented closely following the R package CODA (Plummer et al., 2006). The first is the Multivariate Potential Scale Reduction Factor (MPSRF, Gelman & Rubin, 1992; Brooks & Gelman, 2012). This method requires multiple chains, whose initial values must be over-dispersed in comparison with the posterior, and quantifies the mixing of the walkers by comparing the ensemble variance to the per-walker variance. Nevertheless, we do not know a priori the posterior and, for this reason, we may only guess what an over-dispersed distribution would be.

The second method is the Heidelberger-Welch diagnostic test (Heidelberger & Welch, 1981, 1983), which consists in applying the Schruben stationarity test (Schruben, 1982) to subsets of a chain to obtain one that satisfies the test for a given $p$-value. Since we are using an ensemble sampler, we can apply this test to each individual chain, or, more efficiently, to the ensemble mean of each parameter. We applied the individual approach only when the Markov chain presents convergence problems.

In the third approach, we calculated the autocorrelation time as proposed by Goodman & Weare (2010). However, in NumCosmo, instead of estimating the autocorrelation time directly from the autocorrelations, we fit an Auto Regressive (AR) model as in CODA. In the AR model fitting, we use the bias corrected Akaike Information Criterion (AICc) (HURVICH & TSAI, 1989) to choose the best AR order to use for a given parameter in a chosen chain. This provides a less noisy estimate of the autocorrelation time than the direct inference from the autocorrelations (see Goodman & Weare, 2010). The Effective Sample Size (ESS) is computed using that estimated autocorrelation time and provides an equivalent measure of the effective number of independent points in each chain. Finally, the variance of the sample mean of the parameters is given by the empirical variance of the sampled values divided by the ESS.