# Cosmological measurements from angular power spectra analysis of BOSS DR12 tomography

Arthur Loureiro, Bruno Moraes, Filipe B. Abdalla, Andrei Cuceu, Michael Mcleod, Lorne Whiteway, Sreekumar T. Balan, Aurelién Benoit-Lévy, Ofer Lahav, Marc Manera, et al.

► To cite this version:

Arthur Loureiro, Bruno Moraes, Filipe B. Abdalla, Andrei Cuceu, Michael Mcleod, et al.. Cosmological measurements from angular power spectra analysis of BOSS DR12 tomography. Monthly Notices of the Royal Astronomical Society, 2019, 485 (1), pp.326-355. 10.1093/mnras/stz191. hal-01891232

HAL Id: hal-01891232

https://hal.science/hal-01891232

Submitted on 1 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmological measurements from angular power spectra analysis of BOSS DR12 tomography

Arthur Loureiro<sup>1,1*</sup> Bruno Moraes,<sup>1,1*</sup> Filipe B. Abdalla,<sup>1,2*</sup> Andrei Cuceu,<sup>1</sup>

Michael McLeod,<sup>1</sup> Lorne Whiteway,<sup>1</sup> Sreekumar T. Balan,<sup>1</sup> Aurélien Benoit-Lévy,<sup>3</sup>

Ofer Lahav,<sup>1</sup> Marc Manera,<sup>4,5</sup> Richard P. Rollins<sup>6</sup> and Henrique S. Xavier<sup>7</sup>

Accepted 2018 December 25. Received 2018 November 16; in original form 2018 September 5

## ABSTRACT

We constrain cosmological parameters by analysing the angular power spectra of the Baryon Oscillation Spectroscopic Survey DR12 galaxies, a spectroscopic follow-up of around 1.3 million SDSS galaxies over 9376 deg<fnref n="2" /> with an effective volume of ~6.5 (Gpc h⁻<fnref n="1" />)<fnref n="3" /> in the redshift range 0.15 ≤ z < 0.80. We split this sample into 13 tomographic bins (Δz = 0.05); angular power spectra were calculated using a Pseudo-Cℓ estimator, and covariance matrices were estimated using lognormal simulated maps. Cosmological constraints obtained from these data were combined with constraints from Planck cosmic microwave background experiment as well as the JLA supernovae compilation. Considering a wCDM cosmological model measured on scales up to kmax = 0.07h Mpc⁻<fnref n="1" />, we constrain a constant dark energy equation-of-state with an ~4 per cent error at the 1σ level: w₀ = -0.993<sup>+0.046</sup><sub>-0.043</sub>, together with Ωm = 0.330 ± 0.012, Ωb = 0.0505 ± 0.002, S8 ≡ σ8√Ωm/0.3 = 0.863 ± 0.016, and h = 0.661 ± 0.012. For the same combination of data sets, but now considering a ΛCDM (lambda cold dark matter) model with massive neutrinos and the same scale cut, we find Ωm = 0.328 ± 0.009, Ωb = 0.05017<sup>+0.0009</sup><sub>-0.0008</sub>, S8 = 0.862 ± 0.017, and h = 0.663<sup>+0.006</sup><sub>-0.007</sub>, and a 95 per cent credible interval (CI) upper limit of Σmν < 0.14 eV for a normal hierarchy. These results are competitive if not better than standard analyses with the same data set, and demonstrate that this should be a method of choice for future surveys, opening the door for their full exploitation in cross-correlation probes.

**Key words:** neutrinos – dark energy – large-scale structure of Universe – cosmology: observations.

## 1 INTRODUCTION

Since the discovery that the expansion of the Universe is accelerating, made by two independent Type Ia Supernovae analyses in the late 1990s (Riess et al. 1998; Perlmutter et al. 1999), the main questions in cosmology have concerned the nature of dark

energy and dark matter. Uncertainties in some of the properties of these components have been reduced to the few per cent level, and recent results from an array of different probes support a standard cosmological model with no spatial curvature and a cosmological constant driving acceleration (DES Collaboration 2017; Hildebrandt et al. 2017; Planck Collaboration 2018). Despite some suggested external and internal tensions in cosmological data sets (Efstathiou & Lemos 2018; Planck Collaboration 2018; Riess et al. 2018), this so-called flat ΛCDM (lambda cold dark matter) model is supported by observations from galaxy photometric surveys such

<footnote n="1"><sup>1</sup>Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown, 6140, South Africa</footnote>

<footnote n="3"><sup>3</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="4"><sup>4</sup>Institut de Física d'Altes Energies, The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="5"><sup>5</sup>Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="6"><sup>6</sup>Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="7"><sup>7</sup>Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão, São Paulo 05508-090, Brazil</footnote>

<footnote>* E-mail: arthur.loureiro.14@ucl.ac.uk (AL); b.moraes@ucl.ac.uk (BM); fba@star.ucl.ac.uk (FBA)</footnote>
as KiDS (Hildebrandt et al. 2017) and DES (DES Collaboration 2017), cosmic microwave background (CMB) fluctuations (Planck Collaboration 2018), type Ia supernovae (Scolnic et al. 2018), and galaxy cluster counts (de Haan et al. 2016), among others. At present, there is no overwhelming statistical evidence from cosmological data requiring any model extensions.

Cosmology has proven a fertile area for progress in measuring known matter components of our Universe. Measuring neutrino masses and the neutrino hierarchy is one of the great modern challenges in physics (Kajita 2016; McDonald 2016). The best current constraints from particle physics experiments place a lower limit of 0.06 eV to the sum of neutrino masses in the normal hierarchy (e.g. Esteban et al. 2017). Due to their effect of smoothing matter perturbations on the primordial Universe, the neutrino masses can potentially be measured using probes of the large-scale structure of the Universe (Lesgourgues et al. 2013). Upper bounds on the sum of their masses are currently reaching a level where strong constraints on their hierarchy are possible (Planck Collaboration 2018), and a measurement of their masses is tantalisingly close. Innovative strategies for the analysis of cosmological data could help to overcome the remaining hurdles.

Recent years have seen increased interest in measuring cross-correlations of distinct cosmological probes. Simultaneously modelling and fitting auto- and cross-correlations of observable cosmological fields can improve the dark energy figure-of-merit of surveys (Wang 2008), provide better control of systematic errors, and potentially unveil new physics (e.g. Kirk et al. 2015). Examples of this approach include combinations of CMB primary and secondary anisotropies with galaxy clustering and cosmic shear signals that help to constrain galaxy bias and intrinsic alignments (Hand et al. 2015; Giannantonio et al. 2016), ‘3x2pt’ correlations between galaxy clustering and lensing signals that provide the strongest low-redshift constraints on cosmological models (DES Collaboration 2017; Hildebrandt et al. 2017), and also between galaxy clustering and CMB (Nicola, Refregier & Amara 2016a, 2017; Doux et al. 2017).

A consistent treatment of all probes requires a common theoretical framework for the analysis of the data and covariance matrices across the different correlations. A natural candidate for this is the angular power spectrum. It has been in widespread use by the CMB community for decades (Fixsen et al. 1996; Szapudi et al. 2001; Wandelt, Hivon & Górski 2001; Hivon et al. 2002; Górski et al. 2005), providing several advantages over other statistical estimators. Spherical harmonic decompositions are particularly suited to the analysis of data on the sphere, as they are easily connected to the underlying linear cosmological perturbations in a statistically isotropic and homogeneous Universe, and possess a simple covariance structure for most practical cases despite mode mixing from partial sky observations. Construction of the estimator from galaxy survey data does not require any de-projection using cosmological information, and covariance estimation from lognormal simulations can be estimated in a cosmology-independent way. This allows for a consistent end-to-end analysis. Last, but not least, self-calibration of photometric redshift distributions using cross-correlations with spectroscopic surveys is more readily implemented, and more robust to potential systematic errors (McQuinn & White 2013; McLeod, Balan & Abdalla 2017) when compared to other methods such as $P(k)$, $\xi(r)$, and $w(\theta)$ (Ross et al. 2017; Salazar-Albornoz et al. 2017). We argue that this should be the case because methods that live in angular space such as angular power spectra and $w(\theta)$ can be naturally binned finely and hence more information about the redshift evolution can be extracted without further modelling and

further assumptions. We further argue that non-linearities are better separated in our angular power spectra method than they would be if using the data in configuration space.

Spectroscopic surveys give precise information about the radial distances to galaxies, since the redshifts can be precisely measured from the spectra. In light of the precision in redshift for such galaxy surveys, the usual cosmological approach is the use the 3D power spectrum, $P(k)$, or the 3D correlation function in real space, $\xi(r)$ (Percival et al. 2001; Beutler et al. 2017a; Ross et al. 2017; Wang et al. 2017). Although these approaches have some advantages related to exploring the full radial information from spectroscopic surveys, a fiducial cosmology always needs to be assumed in order to translate from redshift space to real space. This choice of fiducial cosmology may potentially bias cosmological measurements, justifying once more the choice of a tomographic angular power spectra analysis.

However, there are difficulties involved in using angular power spectrum estimators on a spectroscopic galaxy survey. First, it is not simple to ensure that all of this radial information is contained in the angular power spectra of projected redshift bins – even if a fine redshift binning strategy is employed. A second and more relevant issue is that spectroscopic surveys have a much lower galaxy density due to necessarily long integration times and targeting of specific galaxies with fibre spectrographs. This leads to a low signal-to-noise ratio of galaxies to Poisson noise once the data are projected in several tomographic redshift bins. A judicious choice of redshift bin width and Fourier scales can ensure that all relevant linear cosmological information is retrieved (Asorey et al. 2012; Gaztañaga et al. 2012; Eriksen & Gaztañaga 2015; Kirk et al. 2015), but no consistent application of 2D angular power spectra tomography with multiple narrow bins has been attempted on real spectroscopic survey data.<fnref n="1" />

In this work, we apply the angular power spectrum formalism to the Baryon Oscillation Spectroscopic Survey (BOSS) 12th and final public data release. The BOSS is one of the components of the third phase of the Sloan Digital Sky Survey (SDSS-III). Its main aim is to measure the preferred scale of baryonic acoustic oscillations in the primordial baryon-photon plasma, as imprinted in the late-time galaxy distribution. The DR12 data release contains the largest spectroscopic catalogue to date (Alam et al. 2015). It is based on observations of around 2.5 million objects of which around 1.5 million were classified as galaxies, which are further selected to form a large-scale structure galaxy sample ready for cosmological analysis (Reid et al. 2016). We choose to work with this data set because of its constraining power, its public availability, and because of the possibility of comparing our results to those previously obtained by the BOSS collaboration with this same data set (see Alam et al. 2016 and the BOSS publications website<fnref n="2" /> for a list of cosmological publications from the collaboration).

Using the BOSS large-scale structure sample, we show that it is possible not only to measure the full shape of the angular power spectra in very thin tomographic redshift bins, but also to obtain reliable cosmological constraints for $\Lambda$CDM, $w$CDM, and $\Lambda$CDM with $\sum m_\nu$ cosmologies using such a survey alone. The method presented here uses the full shape of the angular power spectra – not just the BAO scale. This is achieved by separating the galaxy samples into tomographic redshift bins with $\Delta z = 0.05$, and using

<footnote n="1">¹Salazar-Albornoz et al. (2017) perform a similar analysis in real space with the BOSS DR12 galaxies.</footnote>

<footnote n="2">²http://www.sdss3.org/science/boss publications.php</footnote>
both the auto power spectra and the cross-power spectra of adjacent bins to extract information from the radial correlation of galaxies. Further combination with external CMB (Planck Collaboration XIII 2016c, as the Planck 2018 likelihood was not available by the time we submitted this work) and SNIa (Betoule et al. 2014) data sets achieves competitive constraints on the models mentioned above.

This paper is organized as follows. Section 2 describes the BOSS LSS sample selection criteria, the mask creation, and the construction of the galaxy overdensity maps. Section 3 describes the Pseudo-C$_{l}$ estimator used for the angular power spectrum analysis. Section 4 describes the theoretical modelling of the angular power spectrum and the use of lognormal mocks for covariance matrix estimation. Section 5 describes our analysis of potential systematic errors using the cross-power spectra between the data and different sources of systematic effects. Section 6 explains the Bayesian modelling for cosmological parameter estimation, describes a series of consistency checks performed on the data, the covariance matrix, and the pipelines, and finally presents cosmological parameter constraints for flat ΛCDM, wCDM, and ΛCDM + $\sum m_v$ models using the BOSS C$_{ls}$ alone and in combination with external CMB results from the Planck collaboration XI 2016b) and type Ia supernovae results from JLA (Betoule et al. 2014). This is a publication by the ZXCorr Collaboration. All data products, data vectors, covariances, and cosmological chains will be made available in the ZXCorr website: http://zxcorr.org.

## 2 BOSS DR12 DATA

The BOSS Data Release 12 (BOSS DR12) is the result of one of the experiments in the third phase of the Sloan Digital Sky Survey (SDSS-III); it is the largest spectroscopic redshift galaxy survey to date. See Alam et al. (2015) for a full description of BOSS DR12 (and in particular for descriptions of the target selection criteria and of the object weighting scheme for offsetting various systematic effects).

The BOSS DR12 is subdivided into two main samples: LOWZ and CMASS. The BOSS collaboration created these samples by applying colour-magnitude and colour-colour cuts to the SDSS photometric catalogue in order to generate lists of targets for spectroscopic observation. The LOWZ sub-sample is designed as a simple extension of the original SDSS Luminous Red Galaxy (LRG) sample (Eisenstein et al. 2001) at low redshifts, while the CMASS sample is defined to select a stellar mass-limited sample of galaxies of all colours – hence its name, for ‘constant stellar mass’ – complemented by a colour cut whose goal is to select higher-redshift objects. The targets were then observed spectroscopically and objects that revealed themselves not to be galaxies (e.g. stars or quasars) were discarded. For a comprehensive discussion of the photometric cuts, selection criteria, and the terminology used, see Dawson et al. (2013).

## 2.1 Galaxy catalogues

To facilitate comparison of our results with the official BOSS collaboration results (Alam et al. 2016), the construction of the catalogues used in this work followed a procedure similar to that outlined in Reid et al. (2016). The data set was downloaded from the BOSS DR12 website.<fnref n="3" /> We have further restricted these samples by applying redshift cuts of $0.15 \le z < 0.45$ for LOWZ and 0.45

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The redshift distribution of the BOSS samples. The darker histograms correspond to the total samples for LOWZ ($0.15 \le z < 0.45$) and CMASS ($0.45 \le z < 0.80$). The overlap between samples was excluded using the redshift information, leaving a total number of 1157755 galaxies.</figcaption></figure>

$z < 0.80$ for CMASS. These cuts ensure that our LOWZ and CMASS samples do not overlap in redshift, which simplifies our tomographic analysis. We use $z = 0.45$ (and not a lower $z$) as the dividing point between the two samples because the LOWZ sample has around 12 per cent more galaxies in $0.4 < z < 0.45$ than does CMASS. See Fig. 1 for the resulting redshift distributions. Note also that our upper limit of $z < 0.8$ for CMASS is greater than the $z < 0.75$ limit used in Reid et al. (2016). As a result of these factors, our redshift ranges differ from those quoted in Reid et al. (2016) and Alam et al. (2016). The subsections that follow outline the main characteristics of the samples.

### 2.1.1 LOWZ sample

The LOWZ sample contains LRGs with redshifts up to around 0.45 as a extension of the SDSS-I/II LRG Cut I sample (Eisenstein et al. 2001). The targets are selected at low redshifts by a cut around the predicted colour locus (equation 1), and a selection of bright red objects is done at each redshift by a variable colour-magnitude cut in the r band (equation 2). This is the main cut in the LOWZ sample as it produces a constant number density on the redshift range of this sample. According to Reid et al. (2016), the number of galaxies in the sample is extremely sensitive to this cut (see also Ross et al. 2013). Star-galaxy separation is done, for LRGs, with a cut on the r-band magnitudes as shown in equation (3). Finally, to guarantee a high spectroscopic redshift success rate, a cut is performed on the r band to impose a brightness limit, as shown in equation (4).

In summary, the photometric target selection criteria for the LOWZ sample are

$$|c_{\perp}| < 0.2 \qquad (1)$$

$$r_{\text{cmod}} < 13.5 + c_{\parallel}/0.3 \qquad (2)$$

<footnote n="3">³http://data.sdss3.org/sas/dr12/boss/</footnote>
### 2.1.3 Total galaxy weights

Various observational effects, such as fibre collisions, will introduce bias into clustering statistics calculated from raw DR12 data. To offset this, the BOSS collaboration provides a weighting scheme for each object; clustering statistics calculated using object counts weighted by this scheme are then expected to be unbiased by such effects. The scheme is described in Reid et al. (2016), which in turn was based on Anderson et al. (2014). We use the same weighting scheme.

For each targeted galaxy, the BOSS collaboration provides three components to the weighting scheme, corresponding to different observational effects (Anderson et al. 2014; Reid et al. 2016; Ross et al. 2017):

(i) $w_{\text{systot}}$, a combination of stellar density with airmass, sky flux, reddening, and other seeing conditions;

(ii) $w_{\text{cp}}$, which is due to close-pair objects, i.e. pairs of objects that cannot have both their spectra measured due to fibre collisions;

(iii) $w_{\text{noz}}$, which takes into account nearest neighbours following a redshift failure by up-weighting such galaxies.

We follow equation (49) in Reid et al. (2016) to combine these into a single weight for each galaxy (Ross et al. 2017):

$$w_{\text{tot}} = w_{\text{systot}}(w_{\text{cp}} + w_{\text{noz}} - 1). \quad (13)$$

The default values of $w_{\text{cp}}$ and $w_{\text{noz}}$ are unity. By construction the term inside the parentheses in equation (13) conserves the total number of targeted galaxies. A more detailed study of the impact of observational systematics is presented in Section 5.

## 2.2 Masks and map making

We now describe the construction of the maps and masks that are our final data products. In this construction we rely on the HEALPix<fnref n="4" /> software for pixelizing the celestial sphere (Górski et al. 2005). The procedures described here were used for both CMASS and LOWZ.

### 2.2.1 Masks and angular selection function

The BOSS collaboration provides<fnref n="5" /> an acceptance mask and several veto masks; these are in MANGLE format (Swanson et al. 2008).

The acceptance mask is continuous (i.e. takes values between 0 and 1), the value for a given region reflecting the completeness of observations there i.e. the extent to which spectra were obtained for all targets. The precise value is

$$C_{\text{BOSS}} = \frac{N_{\text{obs}} + N_{\text{cp}}}{N_{\text{obs}} + N_{\text{cp}} + N_{\text{missed}}} \quad (14)$$

where:

(i) $N_{\text{obs}}$ is the number of spectroscopically observed objects including galaxies, stars, and unclassified objects;

(ii) $N_{\text{cp}}$ is the number of close-pair objects; and

(iii) $N_{\text{missed}}$ is the number of targeted objects with no spectra.

The veto masks are binary maps (i.e. regions are marked as either good or bad); these maps mask out regions affected by observational factors such as centerpost collisions, collision priorities, bright

$$r_{\text{psf}} - r_{\text{cmod}} > 0.3 \quad (3)$$

$$16 < r_{\text{cmod}} < 19.6 \quad (4)$$

In the first months of observation, the BOSS collaboration used a different star-galaxy separation criterion compared to that used later (see Appendix A from Reid et al. 2016). As a result, some sky regions from the LOWZ sample (specifically LOWZE2 and LOWZE3) have a redshift distribution that differs from that in other regions. Our method relies on having a consistent redshift distribution across the sky, and therefore we excluded these regions from our LOWZ sample (see Fig. 2a).

### 2.1.2 CMASS sample

The CMASS sample was designed to be closer to a mass limited sample, extending the Cut-II LRGs from SDSS-I/II (Eisenstein et al. 2001) to bluer and fainter objects using a sliding colour-magnitude cut as shown in equation (5). The cut in the quantity $d_{\perp}$ (equation 6) results in an increase in the number density of objects for the redshift range of $0.45 < z < 0.80$ (see Fig. 1). Model and magnitude limit cuts (equations 7 and 8) ensure high redshift success rates while preventing low redshift objects from being erroneously targeted. Outliers and problematic blended objects are excluded using cuts in i- and r-band magnitudes (equations 9 and 10). Finally, star-galaxy separation was done by performing a varying cut in $i_{\text{psf}} - i_{\text{mod}}$ and $z_{\text{psf}} - z_{\text{mod}}$ based on Cannon et al. (2006) (equations 11 and 12).

In summary, the CMASS sample photometric target selection criteria for most of the survey are

$$i_{\text{mod}} < \min(19.86 + 1.6(d_{\perp} - 0.8), 19.9) \quad (5)$$

$$d_{\perp} > 0.55 \quad (6)$$

$$17.5 < i_{\text{cmod}} < 19.9 \quad (7)$$

$$i_{\text{fib2}} < 21.5 \quad (8)$$

$$r_{\text{mod}} - i_{\text{mod}} < 2 \quad (9)$$

$$r_{\text{dev,i}} < 20.0 \text{ pix} \quad (10)$$

$$i_{\text{psf}} - i_{\text{mod}} > 0.2(21 - i_{\text{mod}}) \quad (11)$$

$$z_{\text{psf}} - z_{\text{mod}} > 0.46(19.8 - z_{\text{mod}}) \quad (12)$$

Although around 25 000 galaxies were targeted with slightly different selection criteria (see section 3.3.1 from Reid et al. 2016 for further details), this does not affect significantly the sample’s redshift distribution (in the way that it did for LOWZE2 and LOWZE3 samples), and therefore we retain these galaxies.

<footnote n="4">⁴http://healpix.sourceforge.net</footnote>

<footnote n="5">⁵http://data.sdss3.org/sas/dr12/boss/lss/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) LOWZ final pixel completeness angular mask with $N_{side} = 512$. We excluded the LOWZE2 and LOWZE3 regions (the holes in the NGC) due to the non-standard $N(z)$ in these regions (the result of an initially different observing strategy that affected these regions). After performing a pixel completeness cut of 0.8, the total used area of the mask is around 8529.58 deg$^2$ ($f_{sky} = 0.2067$). (b) CMASS final pixel completeness angular mask with $N_{side} = 512$. After performing a pixel completeness cut of 0.8, the total used area of the mask is around 9444.63 deg$^2$ ($f_{sky} = 0.2286$).</figcaption></figure>

stars, bright objects, seeing cuts, extinction cuts, and others (see section 5.1 in Reid et al. 2016).

We transform the BOSS acceptance and veto masks into a high-resolution HEALPix pixelization with $N_{side} = 16384$. Using this pixelization scheme, we combine the BOSS masks to yield a high-resolution binary mask. This is done by accepting pixels in which the acceptance mask value $C_{BOSS}$ exceeds 0.7 and which are not marked as bad in any of the veto masks; other pixels are rejected. This choice of completeness cut is based on the BOSS LSS catalogue algorithm from Reid et al. (2016). This high-resolution binary mask is then degraded to a lower resolution ($N_{side} = 512$) continuous mask with values $C_{pix}$ (the pixel completeness factor), defined for a given pixel to be the fraction of high-resolution sub-pixels that are marked as good in the high-resolution binary mask. This is our final mask product and can be seen in Figs 2(a) for LOWZ and (b) for CMASS. The masks used for the pseudo angular power spectrum estimator (PCL) measurements in Section 3 contains a hard cut in $C_{pix} \ge 0.8$: values <0.8 are set to 0 and values $\ge 0.8$ are set to 1.

### 2.2.2 HEALPix galaxy overdensity maps

From the galaxy catalogues, we generate the final data products to be used in our analysis: the galaxy overdensity HEALPix maps. First, we bin both data catalogues into tomographic redshift bins of $\Delta z = 0.05$. This gives six tomographic bins for LOWZ ($0.15 \le z < 0.45$) and seven for CMASS ($0.45 \le z < 0.80$). Details about each redshift bin can be found in Table 1. According to Asorey et al. (2012), $\Delta z = 0.05$ is the thickest possible redshift bin size a spectroscopic redshift survey with $z < 1$ can have in order to keep sufficient radial information without suppressing the radial BAO information due to averaging originating from mode projection. Smaller bin sizes could improve the quality of radial information; however, the trade-off between bin size and shot-noise per bin for the case considered in this work is such that the shot-noise would then be too high for the considered scales. The use of the cross-power spectra between adjacent bins allows for redshift-space distortion (RSD) information to be properly probed as explained in Section 4.1.2.

Next, we create a weighted number counts map that contains the number of objects in each HEALPix pixel, $n_p$, weighted by the total galaxy weight ($w_{tot}$) given by equation (13). To create the final galaxy overdensity maps, we up-weight the maps by the inverse of the pixel completeness factor from the masks, $C_{pix}$. Here, objects in

<table><caption>**Table 1.** Details of each tomographic redshift bin: redshift limits, the number of objects, and shot noise. Note that shot noise is calculated after applying the galaxy weights (Section 2.1.3, equation 13).</caption><thead><tr><th>Sample bin</th><th>z<sub>min</sub></th><th>z<sub>max</sub></th><th>Num galaxies</th><th>Shot noise (gal/strd)<sup>-1</sup></th></tr></thead><tbody><tr><td>LOWZ-0</td><td>0.15</td><td>0.20</td><td>43 265</td><td>6.143 &times; 10<sup>-5</sup></td></tr><tr><td>LOWZ-1</td><td>0.20</td><td>0.25</td><td>51 271</td><td>5.156 &times; 10<sup>-5</sup></td></tr><tr><td>LOWZ-2</td><td>0.25</td><td>0.30</td><td>59 713</td><td>4.416 &times; 10<sup>-5</sup></td></tr><tr><td>LOWZ-3</td><td>0.30</td><td>0.35</td><td>85 394</td><td>3.064 &times; 10<sup>-5</sup></td></tr><tr><td>LOWZ-4</td><td>0.35</td><td>0.40</td><td>83 537</td><td>3.136 &times; 10<sup>-5</sup></td></tr><tr><td>LOWZ-5</td><td>0.40</td><td>0.45</td><td>57 165</td><td>4.605 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-6</td><td>0.45</td><td>0.50</td><td>177 383</td><td>1.577 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-7</td><td>0.50</td><td>0.55</td><td>217 636</td><td>1.275 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-8</td><td>0.55</td><td>0.60</td><td>179 571</td><td>1.545 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-9</td><td>0.60</td><td>0.65</td><td>114 398</td><td>2.435 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-10</td><td>0.65</td><td>0.70</td><td>57 537</td><td>4.850 &times; 10<sup>-5</sup></td></tr><tr><td>CMASS-11</td><td>0.70</td><td>0.75</td><td>23 631</td><td>1.182 &times; 10<sup>-4</sup></td></tr><tr><td>CMASS-12</td><td>0.75</td><td>0.80</td><td>7253</td><td>3.839 &times; 10<sup>-4</sup></td></tr></tbody></table>

pixels with $C_{pix} < 0.8$ are now considered outside the footprint, i.e. the pixel value is set to zero. Thus, the expression for the overdensity maps is

$$ \delta_{i,p}^{g} = \begin{cases} \left( \frac{1}{C_{pix,p}} \frac{n_{i,p}^{g}}{\bar{n}_{i}} \right) - 1 & \text{if } C_{pix,p} > 0.8 \\ 0 & \text{otherwise,} \end{cases} \qquad (15) $$

where $\bar{n}_i$ is the mean number of weighted galaxies per observed pixel, in each tomographic redshift bin. Note that the weights to which we are referring here are the ones mentioned in equation 13; the $\bar{n}_i$ are not weighted by the pixel completeness weight. After these procedures are applied to all 13 tomographic redshift bins, we are ready to measure the power spectra of these maps using the Pseudo-$C_\ell$ estimator described in the next sections.

## 3 ANGULAR POWER SPECTRA ESTIMATORS AND MEASUREMENTS

The first proposed method for estimating the angular power spectrum $C_\ell$ (Peebles 1973) consists of projecting the density field on to the celestial sphere, decomposing this projected field into spherical harmonics, and then analysing statistically the coefficients of this
decomposition. We refer to this method of estimating the power spectrum as a pseudo power spectrum estimator (PCL).

In Appendix A, we describe the PCL estimator for overdensity in detail, following recent approaches as presented in e.g. Scharf et al. (1992), Fisher, Scharf & Lahav (1994), Wright et al. (1994), Huterer, Knox & Nichol (2001), Hivon et al. (2002), Blake, Ferreira & Borrill (2004), Blake et al. (2007), Thomas, Abdalla & Lahav (2010), Thomas, Abdalla & Lahav (2011a), Thomas, Abdalla & Lahav (2011b), Ho et al. (2012) and Balaguera-Antolinez et al. (2018). In this section, we will briefly summarize how the estimator works and its main properties.

## 3.1 Pseudo-C<sub>ℓ</sub> estimator

Our observations contain both signal and (Poisson) noise; spatial variations in the latter contribute to the measured auto power spectrum and this effect must be removed when estimating the power spectrum of the underlying signal.

As the signal and noise are uncorrelated, the angular power spectra of the signal ($S_ℓ$), data ($D_ℓ$), and noise ($N_ℓ$) are related by

$$S_{\ell} = D_{\ell} - N_{\ell}. \qquad (16)$$

Here $D_ℓ$ is measured from performing a decomposition in spherical harmonics of the pixelized galaxy overdensity field (equations A9 and A7). For most tomographic bins, we can (using the variance of a Poisson distribution) approximate the angular power spectrum of the noise as

$$N_{\ell} \approx \frac{\Delta\Omega_{\text{tot}}}{n_{\text{tot}}^{g}} = \frac{1}{\bar{n}}, \qquad (17)$$

where $\bar{n}$ is the mean number of galaxies per steradian.

Amending equation (A6) to account for pixelization and shot noise yields an estimator $\hat{S}_{\ell}^{ij}$ for the partial sky signal power spectrum between two redshift bins i and j:

$$\hat{S}_{\ell}^{ij} = \frac{1}{w_{\ell}^2} \left[ \left( \frac{1}{(2\ell + 1)} \sum_{m=-\ell}^{\ell} D_{\ell m}^{ij} \right) - N_{\ell} \delta_{ij} \right], \qquad (18)$$

where $w_ℓ$ is the pixel window function given by equation (A11). Note that the estimator in equation (18) is symmetric in i and j. Note also that there is no shot noise contribution for the cross-power spectra ($i \neq j$). The PCL estimator described here uses galaxy overdensity maps instead of the galaxy counts maps used in Blake et al. (2007), Thomas et al. (2011b), and others; Appendix B describes the correspondence between the two approaches. This estimator is unbiased (Peebles 1973) but does not have minimum variance: maximum likelihood estimators such as QML (e.g. Efstathiou 2004) have smaller variance. However, these maximum likelihood estimators are computationally expensive to use; this is why, we use PCL.

## 3.2 Bandwidth binning and measurements

We bin the $ℓ$ values into bins $\Deltaℓ$ of width 8 (so e.g. the first bin is 2 $\le$ $ℓ$ $\le$ 9). For each bin we calculate a weighted average of the $\hat{S}_{\ell}^{ij}$ (weighted by the number of spherical harmonic coefficients):

$$\hat{S}_{\Delta\ell}^{ij} = \frac{\sum_{\ell \in \Delta\ell} (2\ell + 1) \hat{S}_{\ell}^{ij}}{\sum_{\ell \in \Delta\ell} (2\ell + 1)}. \qquad (19)$$

This binning acts on the measurement in a way that decorrelates mixed modes (which arise from the convolution of the true measurement and the survey's angular window function).

We measure the PCL estimator up to $ℓ_{max} = 510$; Fig. 3 shows the results for the auto- and cross-power spectra for LOWZ and CMASS. We do not consider in this work any cross-correlations between the two samples. The figure also shows error bars given by the diagonal of the covariance (estimated in Section 4.3), as well as the splines used to generate the lognormal mocks (Section 4.3). The figure shows that the last two CMASS bins are dominated by shot noise (due to their small density of galaxies). Uncertainty in the characterization of this noise will be included into the theoretical forward modelling presented in Section 4 and marginalized over during the cosmological parameter estimation (Section 6).

# 4 MODELLING OF THEORY AND ESTIMATION OF COVARIANCE MATRICES

## 4.1 Theoretical angular power spectra

Our goal is to use observations to constrain cosmological parameters; as part of this we describe the theory that connects the statistics of the underlying matter field with the measured angular power spectra. Our treatment is similar to that found in the literature (Scharf et al. 1992; Huterer et al. 2001; Padmanabhan et al. 2007; Thomas et al. 2011b; Asorey et al. 2012).

Let $\delta_g(\mathbf{k}, z)$ denote the galaxy density function. Let $\delta_g(\mathbf{k}, z)$ be its Fourier transform; we can write this in terms of the growth function $D(z)$, the bias $b(z)$ (assumed here to be scale-independent), and the Fourier components $\delta(\mathbf{k}, 0)$ of the underlying matter distribution at the current time:

$$\delta_g(\mathbf{k}, z) = D(z)\delta_g(\mathbf{k}) = D(z)b(z)\delta(\mathbf{k}, 0). \qquad (20)$$

The correlation structure of the Fourier transform is

$$\langle \delta_g(\mathbf{k}, z) \delta_g^*(\mathbf{k}', z) \rangle = (2\pi)^3 \delta^{(D)}(\mathbf{k} - \mathbf{k}') P_g(\mathbf{k}, z) \qquad (21)$$

where $P_g(\mathbf{k}, z) = b(z)^2 P(k, z)$ is the power spectrum of the galaxy density field and $P(k, z)$ is the power spectrum of the underlying matter density field.

Integrating the galaxy density along the line of sight, $\hat{\mathbf{n}}$, yields

$$\delta_g(\hat{\mathbf{n}}) = \int_0^\infty \delta_g(\chi(z)\hat{\mathbf{n}}, z) n(z) dz \qquad (22)$$

where $n(z)$ is the normalized redshift-dependent selection function and $\chi(z)$ is the comoving distance. The spherical harmonic components $a_{\ell m}$ of this projected galaxy distribution are

$$a_{\ell m} = \int Y_{\ell m}(\hat{\mathbf{n}}) \delta_g(\hat{\mathbf{n}}) d\Omega \qquad (23)$$

$$= \int Y_{\ell m}(\hat{\mathbf{n}}) \int \delta_g(\chi(z)\hat{\mathbf{n}}, z) n(z) dz d\Omega \qquad (24)$$

$$= \frac{4\pi}{(2\pi)^3} \int b(z) n(z) D(z) \int \delta(\mathbf{k}, 0)^i j_\ell(k\chi(z)) Y_{\ell m}(\hat{\mathbf{k}}) d^3k dz \qquad (25)$$

where $j_\ell(k\chi(z))$ are the spherical Bessel functions (Thomas et al. 2010, 2011b). The final step uses the plane wave expansion and the spherical harmonic addition theorem. We may collect the z dependencies into a window function:

$$W_{g,\ell}(k) = \int b(z) n(z) D(z) j_\ell(k\chi(z)) dz. \qquad (26)$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Measured signal auto- and cross-power spectra for the LOWZ (blue) and CMASS (red) samples (equation 19). The black dashed lines show the estimated Poissonian shot noise (equation 17). The solid grey line shows the deconvolved spline used in Section 4.3 to generate the lognormal mocks for covariance matrix estimation from which the error bars in this figure were estimated. Even though the measured $\hat{S}_{\Delta t}$ had the shot noise removed, note that the last two CMASS bins have a significant part of their signals below the level of Poissonian shot noise.</figcaption></figure>

Using the window function in equation (25) yields a simple expression for the angular power spectrum:

$$C_{\ell}^{ij} \equiv \langle a_{\ell m}^{i} a_{\ell m}^{j*} \rangle \qquad (27)$$

$$= \frac{2}{\pi} \int W_{g, \ell}^{i}(k) W_{g, \ell}^{j}(k) k^2 P(k, 0) dk. \qquad (28)$$

Here we have introduced superscripts i and j to denote different redshift shells and the equation above therefore defines both auto-$C_{\ell}$ (for $i = j$) and cross-$C_{\ell}$ (for $i \neq j$). The same

formalism can be used to obtain the $C_{\ell}$ between two different tracers, between photometric and spectroscopic redshift shells, etc.

This work uses the *Unified Cosmological Library for C_{ell}s*, also referred to as the UCLCL code (Cuceu et al., in preparation). This code obtains the primordial power spectra and transfer functions from the *CLASS Boltzmann code* (Blas, Lesgourgues & Tram 2011), and then applies equation (28) to obtain the angular power spectrum. UCLCL deals with the redshift distribution in more flexible ways than do *CLASS* and *CAMB* (Lewis, Challinor &
Lasenby 2000): it allows for the input $n(z)$ distribution to be a spline and also allows it to be convolved with a Gaussian error function to take into account redshift systematic effects (equation 29 in Section 4.1.1). A comparison between these codes is presented in Appendix C.

### 4.1.1 Spectroscopic redshift distribution and shot-noise modelling

The spectroscopic selection provides a full (un-normalized) $n(z)$ function for both the LOWZ and CMASS samples (see Fig. 1). Binning is achieved by hard cuts on each of these samples in intervals of $\Delta z = 0.05$, with no overlap or gaps between bins. Despite the impressive precision of spectroscopy, to suggest that these bins have no overlap (i.e. that there is no error in the spectroscopic measurement) is unrealistic; there is overlap, and it has a significant impact on the cross-correlations between bins. Spectroscopic errors are modelled within the distribution functions by a convolution with a narrow Gaussian function representing the uncertainty on a given measurement. Such a convolution is given by

$$n^i(z) = \int n_*^i(z - z') \exp \left( -\frac{z'^2}{2\sigma_s^2} \right) dz', \quad (29)$$

where $n_*^i(z)$ is the raw redshift distribution, $\sigma_s$ (the standard deviation of the Gaussian) is a proxy for the spectroscopic measurement error, and $n^i(z)$ is the final redshift distribution to be used in calculations. In practice, the convolution is achieved by means of a fast Fourier transform algorithm, multiplication of the functions in Fourier space, and reverse transform.

This convolution is also used to account for the Fingers-of-God (FoG) effect (Kang et al. 2002; Percival et al. 2011), which could dominate the measurements on $\sigma_s(z)$. This effect, which is similar to that of RSD, arises from the peculiar motions of galaxies within virialized structures. These motions elongate structures in redshift space, smearing out the redshift distribution by the addition of Doppler shift to cosmological redshift (Kaiser 1987). The convolution width $\sigma_s$ models the combined impact of spectroscopic redshift errors and of the FoG effect; $\sigma_s$ is then varied and marginalized over during the cosmological analysis. Due to the sensitivity of the cross- angular power spectra to these effects, a separate $\sigma_s$ is used for each redshift bin (for more details, see Section 6.1).

### 4.1.2 Redshift space distortions

The full large-scale structure window function needs to take into account Redshift Space Distortions (RSD) (Padmanabhan et al. 2007; Thomas et al. 2011b). This effect tends to increase the power for large scales, $\Delta \ell < 60$, due to the mix of redshift and peculiar velocities of galaxies. This creates the illusion that local peculiar motion of galaxies moving towards us appear closer (i.e. appear to be at lower redshifts), while galaxies with peculiar motion moving away from us appear to be ever further away (i.e. they appear to be at higher redshifts). This effect can be easily taken into account by adding the RSD window function (Scharf et al. 1992; Fisher et al. 1994; Padmanabhan et al. 2007; Kirk et al. 2015; McLeod et al. 2017) to equation (28):

$$W_{\ell}^{\text{Tot},i}(k) = W_{g,\ell}^{i}(k) + W_{\text{RSD},\ell}^{i}(k). \quad (30)$$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** A series of effects that impact the angular power spectrum ($0.45 < z \le 0.50$) in a variety of ways. The two solid lines are the linear and non-linear power spectra (Section 4.1.3); they diverge as a function of scale, $\ell$. Dashed lines include the redshift space distortions effect (Section 4.1.2), which increases the power for larger scales (see sub-panel). Dot-dashed lines show the effect of the mixing matrix convolution (Section 4.1.4), which tends to suppress power on all scales. Finally, the black dotted line is a combination of all effects: RSDs, non-linearities, and mixing matrix convolution. Theoretical power spectra were calculated using the UCLCL pipeline (Cuceu et al., in preparation), using a flat cosmology: $b = 1$, $h = 0.6725$, $\Omega_b = 0.0492$, $\Omega_{cdm} = 0.265$, $w_0 = -1.0$, $\tau_r = 0.079$, $\log A_s = 3.093 \times 10^{-10}$, $n_s = 0.965$.</figcaption></figure>

Here the RSD window function is given by

$$ 
\begin{aligned}
W_{\text{DRS},\ell}^{i}(k) &= \frac{\beta^{i}}{k} \int d\chi \frac{dn^{i}}{d\chi} j_{\ell}^{i}(k\chi(x)) \\
&= \beta^{i} \int n^{i}(\chi(z)) \left[ \frac{(2\ell^{2} + 2\ell - 1)}{(2\ell + 3)(2\ell - 1)} j_{\ell}^{i}(k\chi(z)) \right. \\
&\qquad \left. + \frac{\ell(\ell - 1)}{(2\ell - 1)(2\ell + 1)} j_{\ell-2}^{i}(k\chi(z)) \right. \\
&\qquad \left. - \frac{(\ell + 1)(\ell + 2)}{(2\ell + 1)(2\ell + 3)} j_{\ell+2}^{i}(k\chi(z)) \right] d\chi \quad (31)
\end{aligned}
$$

where we defined the redshift distortion parameter, $\beta^i(z) = (d \ln D(z)/d \ln a)/b^i(z) \approx \Omega_m^i(z)/b^i(z)$, to be dependent on the bias of the given redshift shell or tracer. The RSD window function does not account for the FoG effect, which affects small scales due to the virial motion of galaxies inside clusters (Kang et al. 2002); instead, as discussed in Section 4.1.1, the FoG effect is subsumed into the spread of the spectroscopic redshift distribution.

Fig. 4 shows the impact on the angular power spectrum of different effects considered in this section: redshift space distortions, non-linearities (Section 4.1.3), and partial sky convolution with the mixing matrix (Section 4.1.4). Note that for some of these effects, the scale at which they have impact varies with redshift.
### 4.1.3 Non-linear angular power spectra: Halofit

In the UCLCL pipeline, the $C_\ell$ estimation may be extended some way into the non-linear regime by introducing the scale-dependent non-linear overdensity $\delta_{NL}(k, \chi)$, and therefore the corresponding non-linear growth function

$$D_{NL}(k, \chi) = \frac{\delta_{NL}(k, \chi)}{\delta_{NL}(k, 0)}. \qquad (32)$$

The calculation of this non-linear density is extracted from the CLASS code (see Blas et al. 2011; Di Dio et al. 2013), which expresses a ratio

$$R_{NL}(k, \chi) = \frac{\delta_{NL}(k, \chi)}{\delta_{L}(k, \chi)} = \left( \frac{P_{NL}(k, \chi)}{P_{L}(k, \chi)} \right)^{\frac{1}{2}} \qquad (33)$$

of the non-linear perturbations to the linear ($\delta_L(k, \chi)$); the second equality follows from $P = \langle\delta\delta^*\rangle$. This ratio is calculated using the modified HALOFIT of Takahashi et al. (2012) (also employed by CAMB SOURCES; Challinor & Lewis 2011), with additional corrections from Bird, Viel & Haehnelt (2012) for neutrino effects.

The window function in equation (30) contains both the selection function and the growth function, which tracks the ratio of the power spectrum at different redshifts. The non-linear power spectrum is related to the linear, present-day power spectrum by

$$\begin{aligned} P_{NL}(k, \chi) &= R_{NL}^2(k, \chi) P_L(k, \chi) \\ &= R_{NL}^2(k, \chi) D_L^2(\chi) P_L(k, 0). \end{aligned} \qquad (34)$$

This means that the window functions in equation (30) should have an additional factor of $R_{NL}(k, \chi)$ inside the integral over $\chi$. In the case of these very narrow spectroscopic redshift bins, we take

$$R_{NL}(k, z) = R_{NL}(k, \bar{z}) \qquad (35)$$

where $\bar{z}$ is the mean of the redshift bin i.e. we assume that the non-linear ratios vary negligibly over the width of a single bin (but may vary between different bins). This simplifies the calculation of the window function considerably, and is a good approximation when the width of the bin is small. In this case the window functions for the redshift bins are related in a straightforward way to their linear counterparts:

$$W_{NL,\ell}^i(k) = R_{NL}(k, \bar{z}^i) W_{g,\ell}^i(k). \qquad (36)$$

The rest of the calculation may proceed as usual.

### 4.1.4 Partial sky: mixing matrix convolution

When dealing the PCL estimator measurements, partial sky effects mean that we must calculate the convolution of the theory and the survey's angular selection function. It is computationally expensive and unstable to deconvolve this effect from the measurements. This leads to forward modelling, where the experimental systematics are modelled and introduced into the theoretical predictions (Scharf et al. 1992; Fisher et al. 1994; Thomas et al. 2011b). This effect is taken into account through a convolution with the mixing matrix, $R_{\ell\ell'}$ (Hauser & Peebles 1973; Hivon et al. 2002; Brown, Castro & Taylor 2005; Blake et al. 2007):

$$S_{\ell} = \sum_{\ell'} R_{\ell\ell'} C_{\ell'}. \qquad (37)$$

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Slices through the mixing matrices for LOWZ (resp. CMASS) using $\ell' = 200$ (250) and $\ell' = 300$ (350). Amplitudes were normalized by $R_{\ell\ell'}. As expected, the maximum amplitude peaks in the fixed $\ell'$ and approaches zero in a given $\Delta\ell$. The shape of both matrices remains the same as a function of $\ell$. This shape indicates the correlation introduced due to $\ell$-mixing by the convolution between the mask and the true signal, present in the PCL measurements (Fig. 3).</figcaption></figure>

The mixing matrix itself depends only on the survey's geometry through the mask's angular power spectrum

$$W_{\ell} = \sum_{m=-\ell}^{\ell} \frac{|I_{\ell m}|^2}{(2\ell + 1)} \qquad (38)$$

where (see Appendix B)

$$I_{\ell m} = \sum_{p}^{N_{pix}} Y_{\ell m}^{*}(\theta_p, \phi_p) \Delta\Omega_p; \qquad (39)$$

the mixing matrix is then

$$R_{\ell\ell'} = \frac{2\ell' + 1}{4\pi} \sum_{\ell''} (2\ell'' + 1) W_{\ell''} \begin{pmatrix} \ell & \ell' & \ell'' \\ 0 & 0 & 0 \end{pmatrix}^2. \qquad (40)$$

The $2 \times 3$ matrix above is the Wigner $3j$ function; these coefficients were calculated using the WIGXJPF library (Johansson & Forssén 2015). The mixing matrices are shown in detailed slices in Fig. 5, which gives an intuition about the size of $\Delta\ell$ bands used to bin the measured $\hat{S}_\ell$ as it shows the range of multipoles that are mixed due to the survey's mask. These small correlations between the multipoles can be 'washed away' by binning our measurements. In addition, as can be seen in Fig. 4, the mixing matrix convolution tends to suppress power in all scales.

Finally, after being convolved with the mixing matrix (equation 37), the theoretical $S_\ell$ is binned in the same way as the data in equation (19):

$$S_{\Delta\ell}^{ij} = \frac{1}{\sum_{\ell'}^{\ell' + \Delta\ell} (2\ell' + 1)} \sum_{\ell'} (2\ell' + 1) S_{\ell}^{ij}. \qquad (41)$$

## 4.2 Data theoretical covariance

We follow here the formalism developed in Dahlen & Simons (2008) for the covariance of spectral estimation on a sphere. For clarity we first re-derive some of the results from Section 3 from a different perspective, that of projectors in pixel space.

Consider a data vector $d$ that is a sum of signal and noise ($d(r) = s(r) + n(r)$) and that has a covariance, $D$, that is a combination of
signal covariance S and a noise covariance N. In pixel space, the
data covariance can be expressed as

$$
\mathcal{D} = \langle \mathbf{s}\mathbf{s}^T \rangle + \langle \mathbf{n}\mathbf{n}^T \rangle = \sum_{\ell} (S_{\ell} + N_{\ell})\mathcal{P}_{\ell} \quad (42)
$$

where $\mathcal{P}_\ell$ is the projector in pixel space, defined as

$$
\mathcal{P}_{\ell} = \sum_{m} Y_{\ell m}(\mathbf{r}) Y_{\ell m}^{*}(\mathbf{r}). \quad (43)
$$

The projector satisfies the following identity in the full sky case:

$$
\mathrm{tr}(\mathcal{P}_{\ell} \mathcal{P}_{\ell'}) = (\Delta \Omega_{\mathrm{p}})^{-2} (2\ell + 1) \delta_{\ell\ell'}. \quad (44)
$$

Using this identity, the Pseudo-C<sub>ℓ</sub> estimator from equation (18) can
be written in terms of the projector and data covariance from (42):

$$
\hat{S}_{\ell} = \frac{\Delta \Omega_{\mathrm{p}}^{2}}{(2\ell + 1)} [\mathbf{d}^{\mathrm{T}} \mathcal{P}_{\ell} \mathbf{d} - \mathrm{tr}(\mathcal{N} \mathcal{P}_{\ell})]. \quad (45)
$$

where $\Delta\Omega_p$ is the area of the pixels. Assuming a Gaussian signal (Blake et al. 2007), the covariance matrix for the angular power spectra estimator between different multipoles $\ell$ and $\ell'$ can be expressed as

$$
\Sigma_{\ell\ell'} = \mathrm{Cov}(\hat{S}_{\ell}, \hat{S}_{\ell'}). \quad (46)
$$

The symmetry of the $\mathcal{P}_\ell$ and $\mathcal{D}$ matrices, together with the definition
$\mathrm{Cov}(\mathbf{X}, \mathbf{X}') = (\mathbf{X}\mathbf{X}^\top) - (\mathbf{X})\langle\mathbf{X}\rangle\langle\mathbf{X}'\rangle$, allows us to rewrite the covariance
as

$$
\Sigma_{\ell\ell'} = \frac{2(\Delta\Omega_p)^4}{(2\ell+1)(2\ell'+1)} \mathrm{tr}(\mathcal{D}\mathcal{P}_\ell\mathcal{D}\mathcal{P}_{\ell'}) \quad (47)
$$

This expression works for both full and partial sky cases. The
difference between the two cases appears on the projector identity
from equation (44). Using the definition of the pixel space projector
(equation 43), the $I_{\ell m}$ expression from equation (39), and the fact
that the spectra we consider are moderately coloured, which means
that the spectra do not vary drastically within the range considered
(Dahlen & Simons 2008), one can rewrite equation (46) for a partial
sky observation with area $\Delta\Omega_{\text{tot}}$ as

$$
\begin{align}
\Sigma_{\ell\ell'} &= \frac{1}{2\pi} \left( \frac{4\pi}{\Delta\Omega_{\text{tot}}} \right)^2 (S_{\ell} + N_{\ell})(S_{\ell'} + N_{\ell'}) \nonumber \\
&\quad \times \sum_{\ell''} (2\ell'' + 1) W_{\ell''} \begin{pmatrix} \ell & \ell' & \ell'' \\ 0 & 0 & 0 \end{pmatrix}^2 \nonumber \\
&= \frac{2}{f_{\text{sky}}(2\ell' + 1)} (S_{\ell} + N_{\ell})(S_{\ell'} + N_{\ell'}) R_{\ell\ell'} \tag{48}
\end{align}
$$

where the last equality uses the definition of the mixing matrix
from equation (40) and $f_{\text{sky}}$ is the observed fraction of the sky.
This expression is similar to ones used in Blake et al. (2007),
Padmanabhan et al. (2007), and Thomas et al. (2011b), but has
been extended to account for the mixing of modes due to the mask.

However, this expression (derived by Dahlen & Simons 2008) accounts neither for the pixel window function effect nor for cross-correlations between tomographic redshift bins. To include these effects, we generalize the data angular power spectra by changing $S_{\ell} + N_{\ell}$ to $w_{\ell}^2 S_{\ell}^{ij} + N_{\ell} \delta_{ij} = D_{\ell}^{ij}$ and we include the effect of cross-correlation in the covariance by changing $(S_{\ell} + N_{\ell})(S_{\ell'} + N_{\ell'})$ to $\frac{1}{2}[D_{\ell}^{ij} D_{\ell'}^{ij} + D_{\ell}^{ii} D_{\ell'}^{jj}]$ in equation (48) (Rassat et al. 2007).

The final expression for the angular power spectra theoretical
covariance matrix is therefore

$$
\Sigma_{\ell\ell'}^{ij} = \frac{1}{f_{\text{sky}}(2\ell' + 1)} [D_{\ell}^{ij} D_{\ell'}^{ij} + D_{\ell}^{ii} D_{\ell'}^{jj}] R_{\ell\ell'} \quad (49)
$$

<figure><img src="image_7.png" /><figcaption>**Figure 6.** A typical example of the validation of the covariance matrix obtained from the FLASK lognormal simulations; here we show the result for the first auto-power spectrum for CMASS. We show the analytical expression for the angular power spectrum variance (equation 50) and the variance from the simulations. The bottom panel shows the relative error for this example. All 49 measured $C_\ell$s from Fig. 3 were validated in this way and no trends were apparent.</figcaption></figure>

$$
= \frac{R_{\ell\ell'}}{f_{\text{sky}}(2\ell' + 1)} \left[ (w_{\ell}^2 S_{\ell}^{ij} + N_{\ell}\delta_{ij})(w_{\ell'}^2 S_{\ell'}^{ij} + N_{\ell'}\delta_{ij}) \right. \\
\left. \qquad + (w_{\ell}^2 S_{\ell}^{ii} + N_{\ell}\delta_{ii})(w_{\ell'}^2 S_{\ell'}^{jj} + N_{\ell'}\delta_{jj}) \right] \quad (50)
$$

By performing these modifications, equation (50) recovers the
variance expression for cross-power spectra from Rassat et al.
(2007) when considering just the diagonal; and recovers the original
expression by Dahlen & Simons (2008) when considering just
the auto-power spectrum. Fig. 6 shows a comparison between the
variance (the diagonal) of equation (50) with the variance from
the estimated covariance matrix from Section 4.3. Note that the
covariance between different angular power spectra is considered
to be zero, i.e. $\Sigma_{\ell\ell'}^{ij,i'j'} = \Sigma_{\ell\ell'}^{ij}\delta_{ii'}\delta_{jj'}$.

4.3 Covariance matrices using log-normal mocks

We seek to constrain cosmological parameters using observations; one of the requirements of this process is accurate covariance matrices. Covariances can be estimated using galaxy clustering simulations that reflect not only the cosmology but also systematic effects and observational artefacts. Previous works have used either Gaussian realizations (Blake et al. 2007; Thomas et al. 2011b; Nicola, Refregier & Amara 2016b) or the mocks provided by the BOSS Collaboration (Manera et al. 2013; Kitaura et al. 2016). However this work instead uses lognormal simulations. The decision not to use the official BOSS PATCHY mocks from Kitaura et al. (2016) was made due to the different choice of redshift ranges for our samples: the CMASS PATCHY mocks do not contain galaxies beyond redshift $z = 0.75$, whereas the samples we selected extend to $z = 0.80$ (as described in Section 2).

We generated our mocks using FLASK<fnref n="6" /> (Xavier, Abdalla & Joachimi 2016), a publicly available code that produces lognormal simulations of correlated fields on the sphere. We used the data

<footnote n="6">⁶http://www.astro.iag.usp.br/~flask/</footnote>
$\hat{S}_\ell$ measurements (Section 3) as inputs for the simulations; this allows us to reproduce systematic effects, RSD, non-linear power spectra, and other known and unknown effects that may be present in the data (with no need to model the effects nor to assume any fiducial cosmology). This is a main benefit of this approach to covariance estimation: any effects present in the measured angular power spectra will be reproduced in FLASK's simulations via the $S_\ell$ s measured from the data.

For each sample we produced 6000 lognormal mocks to estimate the data covariance matrix. These mocks were also Poisson sampled to reproduce noise properties and radial and angular selection effects.

The data covariance matrix was produced as follows:

(1) Produce a spline, $\tilde{S}(\ell)$, using the $\hat{S}_{\Delta\ell}$ measurements (Fig. 3) and a Gaussian filter to smooth the measurements.

(2) Deconvolve the mixing matrix $R_{\ell\ell'}$ from the splines to obtain

$$ \tilde{C}^{ij}(\ell) = \sum_{\ell'} R_{\ell\ell'}^{-1} \tilde{S}^{ij}(\ell). \qquad (51) $$

(3) Linearly extrapolate the splines to $\ell_{max} = 8192$ (necessary to allow FLASK to create high-resolution HEALPix maps).

(4) For each tomographic redshift bin, produce FLASK partial sky galaxy number count mocks with $N_{side} = \ell_{max} = 2048$.<fnref n="7" />

(5) Degrade the mocks to $N_{side} = 512$ to match the $N_{side}$ used when analysing the data.

(6) Produce up-weighted galaxy overdensity maps using the pixel completeness factor (as described in Section 2.2.2).

(7) Run the partial sky PCL estimator; include here the pixel window function correction $w_\ell^2$ (as described in equations (18) and (19)) that arises from the degrading of the maps at step 5.

(8) Measure the covariance of the ensemble of angular power spectra obtained from the simulated data:

$$ C_{\Delta\ell\Delta\ell'}^{ij} = \frac{1}{N_S - 1} \sum_{s=1}^{N_S} \left( S_{\Delta\ell}^{ij,s} - \langle S_{\Delta\ell}^{ij} \rangle \right) \left( S_{\Delta\ell'}^{ij,s} - \langle S_{\Delta\ell'}^{ij} \rangle \right)^T . \quad (52) $$

Here $N_S$ is the number of simulations. To validate the estimated covariance matrix, we compared the diagonal of the covariance matrix in equation (52) with the expression for the theoretical variance for the measured angular power spectra in equation (50); Fig. 6 shows a typical result.

# 5 SYSTEMATICS TESTS

Large-scale survey observations, spread over thousands of observation hours, are taken under a variety of conditions. Turbulence in the atmosphere, sky background brightness, and telescope inclination angle are amongst the factors that can influence image quality and object detection. Other than those atmospheric effects, galactic properties are also at play: extinction from dust within the Milky Way and variations of stellar density, as well as the presence of bright stars, are position-dependant and also have an impact on our ability to detect galaxies. Jointly, those observational factors can create small density fluctuations in the galaxy distribution that can

imprint a statistical signal easily confused with the cosmological large-scale structure fluctuations that we are attempting to measure. This effect has been detected and corrected for in several previous analyses with a range of data sets (Blake et al. 2007; Ross et al. 2011; Thomas et al. 2011b; Ho et al. 2012; Leistedt et al. 2013; Leistedt & Peiris 2014; Elsner, Leistedt & Peiris 2016, 2017; Doux et al. 2017)

In this section, we present the analysis performed on the data to ensure that the measured power spectra are not significantly dominated by any known observational systematic effects. We consider a systematic to have a significant effect on the observed power spectra if the cross-power spectra between them deviate from zero, with a deviation that is bigger than both the data variance and the cross-power spectra variance. We describe the systematic effects considered in our analyses, describe our methods for map creation and cross-spectrum measurement, and give some representative results.

## 5.1 Systematic maps

The Sloan Digital Sky Survey monitors and records observational conditions for every tile of the survey. This information is available as a combined set of two files, one that defines a pixelization of the observed sky in MANGLE format and another that records the observational information for each MANGLE polygon.<fnref n="8" /> The first step is to reconstruct the MANGLE maps for each observational systematic from these files. We use the MANGLE python wrapper<fnref n="9" /> to perform this transformation. Since there is potentially more than one observation in a given region of the sky, there can be multiple values for a given polygon. The SDSS files indicate which amongst multiple options is to be taken as the primary value for the field. We select the IDs from those primary fields, match them to their observational properties in the fields list, and create a new MANGLE mask for each of those properties, which are recorded in the weight of the masks.

We create MANGLE masks of sky background flux, sky variance, and average PSF FWHM in all five photometric bands. We also create a mask of the score of each field, defined by the SDSS collaboration to express ‘observational quality’ as an empirical combination of observational values with processing status flags. Additional observational properties can be found in the Field Table, available from the SDSS SkyServer Schema Browser.<fnref n="10" /> The choice of which systematics to take into account is somewhat arbitrary, as there are correlations between observational properties that make information redundant (Leistedt et al. 2013). We choose to add stellar density and galactic extinction to the systematics listed above, as those have been shown to correlate with galaxy density in several previous analyses (e.g. Thomas et al. 2011b; Elvin-Poole et al. 2017). We construct a bright star catalogue from the SDSS object catalogue with the following cuts:

$$ \begin{align} 18 < r_{\text{psf}} &< 19.5, \notag \\ \text{type} &= 6, \notag \\ r_{\text{psf}} - r_{\text{model}} &< 0.25, \tag{53} \end{align} $$

<footnote n="7">⁷The signal realization maps were sampled using a lognormal transformation. Due to the transformation’s non-linearity, we had to generate mocks with a higher $N_{side}$ & $\ell_{max}$ than the data as the lognormal realizations introduce a damping after a certain $\ell$ (see fig. 18 from Xavier et al. 2016). The simulated data maps also used a $N_{side} = 2048$ version of the masks presented in Section 2.2.1.</footnote>

<footnote n="8">⁸The files, window_unified.fits and window_flist.fits, can be found in http://www.sdss.org/dr12/algorithms/resolve/, together with a detailed description of the construction of the survey geometry and of the score quantity described further in the text.</footnote>

<footnote n="9">⁹https://github.com/mollyswanson/manglepy</footnote>

<footnote n="10">¹⁰http://skyserver.sdss.org/dr12/en/help browser browser.aspx</footnote>
where the extinction-corrected magnitude cut ensures robust star selection (Padmanabhan et al. 2007), the type selection is the standard SDSS star–galaxy classifier,<fnref n="11" /> and the magnitude-difference cut is an additional point-source selection performed by the GAMA survey (e.g. Christodoulou et al. 2012). For galactic extinction, we create a map directly in HEALPix format. For simplicity, we take advantage of a PYTHON implementation of extinction E(B − V) value retrieval and map creation.<fnref n="12" /> We use the original SFD scaling (Schlegel, Finkbeiner & Davis 1998).

The MANGLE masks created from the SDSS FITS files are not appropriately snapped, pixelized, and balkanized, which breaks the local character of the MANGLE procedure (Swanson et al. 2008). As a consequence, further operations suffer from impractically large processing times. We therefore run all the steps of the MANGLE pixelization anew, which corrects whatever imperfections remained in the first pass. From these masks, we create full-sky HEALPix maps at resolution N<sub>side</sub> = 16384, which defines an angular scale much smaller than the average resolution of the mask features. For each observational systematic, we populate the sub-resolution HEALPix pixels with values from the associated MANGLE mask. The resulting HEALPix maps encapsulate all the information contained in the original footprint description.

Once the HEALPix systematics maps are created, the next step is to transform them into overdensity maps using the same procedure outlined in Section 2.2.2 for the data (Leistedt et al. 2013). The idea is to treat the systematic maps in the same way as the data in order to apply the statistical estimators consistently. Therefore, we degrade the high-resolution maps to the data resolution (N<sub>side</sub> = 512) and up-weight the maps according to the pixel completeness mask that takes into account the holes in the footprint (see Section 2.2.1); we then perform a cut in pixel completeness C<sub>pix</sub> = 0.8 in the maps. This up-weighting is due to the pixelization of the degraded mask, as explained in Section 2.2.1. From these post-processed maps, we create the systematics overdensity maps as

$$
\delta_p^{\text{Sys}} = \begin{cases} \left( \frac{1}{C_{\text{pix},p}} \frac{n_p^{\text{Sys}}}{\bar{n}^{\text{Sys}}} \right) - 1, & \text{if } C_{\text{pix},p} \ge 0.8 \\ 0, & \text{otherwise} \end{cases} \tag{54}
$$

where $n_p^{\text{Sys}}$ is the pixel value for a given systematic and $\bar{n}^{\text{Sys}}$ is the mean value of the map in the observed fraction of the sky. The systematics overdensity maps were created using both the CMASS and LOWZ masks presented in Section 2.2.1. The resulting systematics overdensity maps will be available at the ZXCorr Collaboration website.

5.2 Cross-power spectra between data and systematic maps

For the systematics analysis using cross-power spectra, we will
follow a data analysis similar to that performed for the galaxy
overdensity maps in Section 3. Using equation (9), we decompose
the systematics overdensity maps, δ<sup>sys</sup>, into spherical harmonics.

The estimator for the cross-power spectra between the data
overdensity maps, δ^g, and the systematics can be written as a
modified version of equation (18):

$$
\hat{S}_{\ell}^{\text{gs}} = \frac{1}{(2\ell + 1)w_{\ell}^2} \sum_{m=-l}^{l} \frac{\frac{1}{2} |d_{\ell m}^{\text{s*}} + d_{\ell m}^{\text{s}*} d_{\ell m}^{\text{s}}|}{J_{\ell m}} \quad (55)
$$

where the index *g* stands for a data map and *s* for a systematics map. We then obtained the estimates for the variance of the systematics cross-power spectra by measuring the *Ŝ*<sub>*g*</sub><sup>*gs*</sup> (equation 55) between the systematics maps and the data mocks described in Section 4.3.

We cross-correlated all 13 tomographic redshift bins with all 18 systematic maps, resulting in a total of 234 cross-power spectra. Fig. 7 shows an example for LOWZ-3 and CMASS-9 bins and cross-power spectra for all systematics. From all of these measurements, the majority are consistent with the variance of the data measured from the lognormal simulation (Section 4.3), which lead us to be confident in using the full shape of our measured C<sub>ℓ</sub>s. Note, however, that a few of the large-scale measurements (low-<i>ℓ</i>) in CMASS sample have a small excess in cross-power spectra with stellar overdensity. The first point on the cross-power spectra between some of the systematic maps and most BOSS bins is clearly more than one sigma away from the data’s variance. Due to this excess in correlation with stellar overdensity and the level of cosmic variance on the first <i>ℓ</i>-band, we decided to exclude this first point from our cosmological analysis (see Section 6 for details on the <i>ℓ</i> range used). As for the second <i>ℓ</i>-band (<i>ℓ</i> = 13.5) presenting an excess of correlation between a few bins and stellar overdensity, we found it to be sub-dominant, with no significant impact from this measurement in our cosmological analysis; therefore, we decided to keep it.

6 COSMOLOGICAL ANALYSIS

In this section, we present the cosmological implications from the
measured angular power spectra of BOSS galaxies for flat ΛCDM,
wCDM, and ΛCDM with Σm<sub>ν</sub> models. Using the theoretical
framework and having estimated covariance matrices as described
in Section 4, we performed a Bayesian analysis using the PLINY
(Rollins 2015 and Rollins et al., in preparation) nested sampler and
the *Unified Cosmological Library for Parameter Inference* code, or
UCLPI (Cuceu et al., in preparation). All analyses considered in
this section use the auto-power spectra and the cross-power spectra
from adjacent tomographic bins using the measurements presented
in Section 3. Cross-power spectra between distant bins are not a
part of our final BOSS-C<sub>ℓ</sub> data vector.

6.1 Likelihoods, priors, and Bayes factor

The cosmological analysis performed in this work follows a
standard Bayesian analysis framework as commonly performed in
the literature (e.g. Blake et al. 2007; Thomas et al. 2011b; DES
Collaboration et al. 2017; Hildebrandt et al. 2017).

The posterior distribution of the cosmological parameters, Θ,
given the measured angular power spectra, Ŝ<sub>Δℓ</sub>, and a model M
can be written as a marginalization of the full posterior over the
nuisance parameters, v:

$$
\mathrm{Pr}(\mathbf{\Theta} | \hat{S}_{\Delta \ell}, \mathcal{M}) = \int \mathrm{Pr}(\mathbf{\Theta}, \mathbf{v} | \hat{S}_{\Delta \ell}, \mathcal{M}) \mathrm{d}\mathbf{v} \quad (56)
$$

The full posterior distribution can be written as

$$
\mathrm{Pr}(\mathbf{\Theta}, \mathbf{v} | \hat{S}_{\Delta \ell}, \mathcal{M}) = \frac{\mathcal{L}(\hat{S}_{\Delta \ell} | \mathbf{\Theta}, \mathbf{v}, \mathcal{M}) \pi(\mathbf{\Theta}, \mathbf{v})}{Z(\hat{S}_{\Delta \ell} | \mathcal{M})} \quad (57)
$$

where L(̂S<sub>Δℓ</sub> | Θ, v, M) is the likelihood, π(Θ, v) is the prior on
the sampled parameters, and Z(̂S<sub>Δℓ</sub> | M) is the evidence, which
is calculated using PLINY, a nested sampler used in our analysis
(Feroz & Hobson 2008; Rollins 2015; Rollins et al., in preparation).

<footnote n="11">¹¹http://www.sdss.org/dr12/algorithms/classify/</footnote>

<footnote n="12">¹²https://github.com/kbarbary/sfdmap</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** An example of the systematics analysis described in Section 5.2. Here, we show the cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ–3(CMASS–9) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the $\delta^{\text{Sys}}$ maps with the FLASK mocks produced in Section 4.3; the shaded region shows the variance of the data, which was also obtained from the same mocks. This figure indicates that the shape of the measured power spectra in Fig. 3 is not dominated by any of the systematics considered, as the variance of the cross-power spectra between data and systematics is consistent with the variance of the data’s auto-power spectra. The results for the other bins are similar to the results shown in this figure. Note also that the first $\ell$ band in the stellar overdensity cross-$C_\ell$ is completely out of the acceptable range, which leads us to exclude this data point on all bins for both samples.</figcaption></figure>
If we had access to the true covariance matrix Σ, the likelihood, assumed here to be Gaussian, would be

$$
\mathcal{L}_G(\hat{S}_{\Delta\ell} | \mathbf{\Theta}, \mathbf{v}) = \frac{1}{\sqrt{|2\pi \Sigma|}} \exp \left\{ -\frac{1}{2} [\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{\text{th}}(\mathbf{\Theta}, \mathbf{v})]^T \times \Sigma^{-1} [\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{\text{th}}(\mathbf{\Theta}, \mathbf{v})] \right\} \quad (58)
$$

where $S_{\Delta\ell}^{\text{th}}(\mathbf{\Theta}, \mathbf{v})$ is the theoretical angular power spectra after being convolved with the mixing matrix (equation 37) and binned into the same bandwidths as the data (equation 41).

However, this is not the case when estimating the covariance matrix C from simulations (equation 52). Even though C can be an unbiased estimator of the true covariance Σ, its inverse C⁻¹ is not necessarily an unbiased estimator of the inverse covariance Σ⁻¹, needed to estimate the likelihood in equation (58). Hartlap, Simon & Schneider (2007) proposed to keep using the Gaussian likelihood, and to apply a simple rescaling to the inverse of the estimated covariance matrix in order to de-bias it (Anderson 2003).

$$
\Sigma^{-1} \rightarrow \alpha C^{-1} \qquad (59)
$$

where

$$
\alpha = \frac{N_s - p - 2}{N_s - 1} \qquad (60)
$$

and $N_s$ is the number of simulations and $p$ is the size of the data vector.

More recently, Sellentin & Heavens (2016, hereafter SH16) showed that when replacing the true covariance Σ with an estimated C, one should marginalize over the true covariance conditioned on the estimated one from simulations. The resulting likelihood is no longer Gaussian; instead, the likelihood is now given by a multivariate t-distribution (SH16):

$$
\mathcal{L}_{SH}(\hat{S}_{\Delta\ell} | \mathbf{\Theta}, \mathbf{v}) = \frac{\bar{c}_p}{|C|^{1/2}} \left[ 1 + \frac{(\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{\text{th}})^T C^{-1} (\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{\text{th}})}{N_s + 1} \right]^{-\frac{N_s}{2}} \quad (61)
$$

where

$$
\bar{c}_p = \frac{\Gamma\left(\frac{N_s}{2}\right)}{\left[\pi(N_s - 1)\right]^{p/2} \Gamma\left(\frac{N_s - p}{2}\right)} \quad (62)
$$

and Γ is the gamma function.

Even though the non-linear model described in Section 4.1.3 is sufficiently reliable, we performed cuts in $ℓ_{max}$ for each of the tomographic redshift bins in order to exclude non-linear scales. In order to make this choice, we used a fiducial cosmology (the same used in Section 6.2.2) to generate theory $C_ℓ$ and performed a preliminary cut in $ℓ_{max}$ where the percent deviation between the linear and non-linear models was smaller than 5 per cent. We performed robustness checks on the 5 per cent deviation cut choice by extending the cuts to $ℓ_{max}$, which had a deviation up to 20 per cent, and concluded that our cosmological results could be trusted up to a 15 per cent deviation between linear and non-linear theories for this fiducial test. In this paper we present results where this percentage cut is 5 per cent and 10 per cent. For avoidance of doubt, for the fiducial cosmology of choice, applying a 5 per cent implies rejecting the majority of modes $k \gtrsim 0.07h \text{ Mpc}^{-1}$, whereas 10 per cent implies rejecting modes where $k \gtrsim 0.1h \text{ Mpc}^{-1}$. The resulting cuts can be found in Table 2. As for the $ℓ$ cuts for cross-power spectra, we chose the lowest $ℓ_{max}$ between the two relevant bins in order to keep a consistent and conservative cut for each bin.

<table><caption>**Table 2.** Maximum multipole considered in the cosmological analysis for each tomographic redshift bin. All the samples start at $ℓ_{min} = 13.5$ and have a bandwidth of $\Delta ℓ = 8$. When considering the cross-power spectra between bins, the lower $ℓ_{max}$ is used. The $ℓ_{max}^{5 \text{ per cent}}$ column corresponds to a $k_{max} \lesssim 0.07h \text{ Mpc}^{-1}$, and the $ℓ_{max}^{10 \text{ per cent}}$ column corresponds to a $k_{max} \lesssim 0.10h \text{ Mpc}^{-1}$.</caption>
  <thead>
    <tr>
      <th>Sample bin</th>
      <th>z<sub>min</sub></th>
      <th>z<sub>max</sub></th>
      <th>ℓ<sub>max</sub><sup>5</sup> per cent</th>
      <th>ℓ<sub>max</sub><sup>10</sup> per cent</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>LOWZ-0</td>
      <td>0.15</td>
      <td>0.20</td>
      <td>53</td>
      <td>69</td>
    </tr>
    <tr>
      <td>LOWZ-1</td>
      <td>0.20</td>
      <td>0.25</td>
      <td>77</td>
      <td>93</td>
    </tr>
    <tr>
      <td>LOWZ-2</td>
      <td>0.25</td>
      <td>0.30</td>
      <td>93</td>
      <td>109</td>
    </tr>
    <tr>
      <td>LOWZ-3</td>
      <td>0.30</td>
      <td>0.35</td>
      <td>109</td>
      <td>133</td>
    </tr>
    <tr>
      <td>LOWZ-4</td>
      <td>0.35</td>
      <td>0.40</td>
      <td>125</td>
      <td>157</td>
    </tr>
    <tr>
      <td>LOWZ-5</td>
      <td>0.40</td>
      <td>0.45</td>
      <td>141</td>
      <td>173</td>
    </tr>
    <tr>
      <td>CMASS-6</td>
      <td>0.45</td>
      <td>0.50</td>
      <td>157</td>
      <td>221</td>
    </tr>
    <tr>
      <td>CMASS-7</td>
      <td>0.50</td>
      <td>0.55</td>
      <td>165</td>
      <td>237</td>
    </tr>
    <tr>
      <td>CMASS-8</td>
      <td>0.55</td>
      <td>0.60</td>
      <td>189</td>
      <td>261</td>
    </tr>
    <tr>
      <td>CMASS-9</td>
      <td>0.60</td>
      <td>0.65</td>
      <td>197</td>
      <td>277</td>
    </tr>
    <tr>
      <td>CMASS-10</td>
      <td>0.65</td>
      <td>0.70</td>
      <td>213</td>
      <td>317</td>
    </tr>
    <tr>
      <td>CMASS-11</td>
      <td>0.70</td>
      <td>0.75</td>
      <td>245</td>
      <td>333</td>
    </tr>
    <tr>
      <td>CMASS-12</td>
      <td>0.75</td>
      <td>0.80</td>
      <td>261</td>
      <td>381</td>
    </tr>
  </tbody>
</table>

We used a total of 28 nuisance parameters (v) in the BOSS C<sub>ℓ</sub> likelihood analysis in most of our results for a 5 per cent cut on ℓ<sub>max</sub>. These parameters are a scale independent bias, b(z), for each redshift bin; a redshift error dispersion, σ<sub>s</sub>(z) (equation 29), for each redshift bin that takes into account spectroscopic redshift error and Finger-of-God effects due to shell-crossing (see Section 4.1.1 for details); and an extra shot noise term for bins 11 and 12 that is forward modelled into the theoretical angular power spectrum inside the likelihood as

$$
\hat{S}_{\Delta\ell}^{th} \rightarrow \hat{S}_{\Delta\ell}^{th} + \mathcal{N} \qquad (63)
$$

where N is a constant that takes into account extra shot noise due to the very low number of galaxies in these two redshift bins. In the case of a 10 per cent cut we used two further shot noise nuisance parameters for bins 9 and 10 as we go further into the non-linear regime where the shot noise in those bins dominates over the signal.

We chose flat priors in all Bayesian analysis. These were based on priors used in Betoule et al. (2014), Alam et al. (2016), Planck Collaboration XIII (2016c), and DES Collaboration et al. (2017) and were set equally for all analyses. The prior ranges can be found in Table 3 for all parameters considered in the cosmological analysis in this section: the baryonic matter density (Ω<sub>b</sub>), the cold dark matter density (Ω<sub>cdm</sub>), the amplitude of the primordial power spectrum (A<sub>s</sub>), the spectral index (n<sub>s</sub>), the Hubble constant (h), the equation of state of dark energy (w<sub>0</sub>), the sum of neutrino mass species (Σm<sub>v</sub>), the optical depth at reionization epoch (τ<sub>reio</sub>), the bias of BOSS galaxies as a function of redshift (b(z)), the redshift dispersion parameter for BOSS galaxies (σ<sub>s</sub>(z)), the extra shot-noise for BOSS galaxies (N<sub>f</sub>), the Planck absolute calibration parameter (γ<sub>cal</sub><sup>Planck</sup>), and the absolute magnitude of SNe Ia at peak light in blue band (M<sub>B</sub><sup>JLA</sup>).

Finally, to perform consistency checks between BOSS DR12 and the external data sets described in Section 6.3, we used the Bayes factor. The Bayes factor for the consistency of two data sets (A and B) is given by

$$
R_{A,B} = \frac{P(\vec{A}, \vec{B} | M)}{P(\vec{A} | M) P(\vec{B} | M)} \qquad (64)
$$

or, for three data sets (A, B, and C):

$$
R_{A,B,C} = \frac{P(\vec{A}, \vec{B}, \vec{C} | M)}{P(\vec{A} | M) P(\vec{B} | M) P(\vec{C} | M)} \quad (65)
$$
<table><caption>**Table 3.** Ranges on the flat priors used in all Bayesian analysis. Parameters are divided into two groups: cosmological and nuisance.</caption><thead><tr><th>Parameter</th><th>Prior Range</th></tr></thead><tbody><tr><td>Ω<sub>b</sub></td><td>1 × 10<sup>-3</sup>, 0.3</td></tr><tr><td>Ω<sub>cdm</sub></td><td>0.0, 0.8</td></tr><tr><td>ln 10<sup>10</sup>A<sub>s</sub></td><td>2.0, 4.0</td></tr><tr><td>n<sub>s</sub></td><td>0.87, 1.07</td></tr><tr><td>h</td><td>0.55, 0.91</td></tr><tr><td>w<sub>0</sub></td><td>-3, -0.3</td></tr><tr><td>Σm<sub>v</sub></td><td>0.0, 1.0 eV</td></tr><tr><td>τ<sub>reio</sub><sup>Planck</sup></td><td>0.0, 0.8</td></tr><tr><td>b(z)</td><td>1.1, 3.3</td></tr><tr><td>σ<sub>s</sub>(z)</td><td>1 × 10<sup>-6</sup>, 9 × 10<sup>-3</sup></td></tr><tr><td>N<sub>9</sub>, N<sub>10</sub></td><td>0.0, 1 × 10<sup>-4</sup></td></tr><tr><td>N<sub>11</sub></td><td>0.0, 8 × 10<sup>-5</sup></td></tr><tr><td>N<sub>12</sub></td><td>0.0, 4 × 10<sup>-4</sup></td></tr><tr><td>y<sub>cal</sub><sup>Planck</sup></td><td>0.99, 1.01</td></tr><tr><td>M<sub>B</sub><sup>LA</sup></td><td>-20.0, -18.5.</td></tr></tbody></table>

where *M* is the model, *P*(*A*, *B*|*M*) is the evidence when the model is fitted to both data sets simultaneously, and *P*(*A*|*M*)*P*(*B*|*M*) is the product of the evidences when the model is fitted to each data set individually. Since PLINY is a nested sampler, all our cosmological estimations lead to values for the evidences of each model, data set, and combination of data sets. We are aware that this method has received some criticism as the value of *R*<sub>A,B</sub> depends crucially on the prior volume chosen for the sampling. Although this is true, choosing priors that are physically motivated should lead to values of *R*<sub>A,B</sub> which can be reflective of the real consistency between data sets *A* and *B*. However, a host of new methods have been proposed in the literature that we have not fully investigated in this paper, but that we intend to implement in further analysis with this data set. These methods involve a range of overlap integrals between the separate posteriors and can be better suited to determine the consistency between two or more data sets. We refer the reader to these methods in Charnock, Battye & Moss (2017a), Raveri & Hu (2018), and Feeney et al. (2018).

In order to perform a robust analysis, we implemented the three likelihoods mentioned above: Gaussian, Gaussian plus the Hartlap correction, and the SH16 *t*-distribution likelihood. We observed that the Gaussian + Hartlap correction, the SH16, and the Gaussian likelihoods led to very similar cosmological contours. This is most likely as we have sufficient mocks to reliably estimate our covariance matrix. Fig. 8 shows a comparison between the three likelihoods for a wCDM model, using the BOSS $C_\ell$ only, for $w_0$ and $\Omega_m$. In all of the following results in this section, we use the SH16 likelihood.

## 6.2 Consistency checks

In this section, we perform a series of consistency checks in order to assess the validity of our cosmological parameter estimation pipelines and data samples.

### 6.2.1 Parameter-dependent theoretical covariance matrix

We implemented an alternative likelihood based on a theoretical expression for the covariance matrix (Section 4.2, see equation 50) where the signal angular power spectra, $S_\ell$, also depends on the sampled parameters. Most standard cosmological analysis in the literature (Alam et al. 2016; DES Collaboration 2017; Hildebrandt

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison between the three likelihood methods mentioned in Section 6.1 using the BOSS $C_\ell$ data only for $\Omega_m$ and $w_0$ in a wCDM model: Gaussian (red), Gaussian with Hartlap correction (black), and the SH16 (blue) likelihoods. Note how, given the high number of lognormal simulation used to estimate the inverse of the covariance, the Hartlap correction likelihood, SH16, and Gaussian have equivalent contours even though the sampled parameters and likelihood values are different. It is clear from this analysis that our estimated covariance matrix from Section 4.3 is robust and was estimated with a sufficient number of simulations.</figcaption></figure>

et al. 2017) assume a covariance matrix that is independent of cosmology and thta is estimated for a fiducial simulation. We do not expect to obtain the same cosmological contours from this method as those presented in the sections that follow; however, we do not expect the contours from this method to disagree significantly with the ones obtained with our estimated covariance matrix. Fig. 9 shows the results for this test for a ΛCDM cosmological model.

### 6.2.2 Controlled cosmology pipeline test

For this test, we generated theory auto- and cross-$C_\ell$ s to mimic the BOSS data set using a Planck-like cosmology: $h = 0.6725$, $\Omega_b = 0.0492$, $\Omega_m = 0.314$, $w_0 = -1.0$, $S_8 = 0.830$, $n_s = 0.96575$. We simulated these fiducial power spectra using the BOSS redshift distribution $n(z)$ from Fig. 1. We chose the nuisance parameters to match the best fit values obtained in Section 6.4 from the combination of the entire cosmological data set available to us. Using BOSS masks as input, we created FLASK mocks like described in Section 4.3: generating the mocks at higher resolution, degrading them, and creating galaxy overdensity maps. We applied the PCL estimator on the 13 overdensity maps, calculating the auto- and cross- power spectra as described for the data in Section 3.

Finally, we ran a cosmological parameter estimation for a ΛCDM cosmology, varying also the 28 BOSS nuisance parameters and using the theoretical covariance matrix as in the previous section. The results are shown in Fig. 10, where the recovered parameters are within the errors with no indications of biases in the entirety of the pipeline.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison between ΛCDM cosmologies recovered using the covariance matrix estimated in Section 4.3 (green contours, same as the ones from Section 6.4), and the cosmology dependent theoretical covariance matrix from equation (50) (blue contours). Note that the same parameters were sampled in both cases as the theoretical covariance also depends on the same nuisance parameters. These marginalized credible intervals (CI) 1σ and 2σ plots indicate both the estimated covariance matrix and the UCLPI pipeline robustness.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Cosmological constraints recovered from a controlled cosmology pipeline test. The dashed lines denote the Planck-like cosmology used as input in the simulations analysed in the blue contours. All parameters agree within the estimated errorbars.</figcaption></figure>

### 6.2.3 Internal checks: single redshift bin consistency

To test the data's internal consistency, we performed a full cosmological analysis in each individual redshift bin from our BOSS samples. The test was performed using a ΛCDM model, varying the same nuisance parameters as described in Section 6.1 for each bin: redshift dispersion, bias, and extra shot noise for the last two CMASS bins. If each individual bin is consistent with all others, this indicates that one can obtain cosmological constraints from the combination of the individual bins. This is shown in Fig. 11 for the posterior projections of Ω_m and Ω_b. In these figures, all contours overlap and, even though some tomographic redshift bins prefer a secondary peak, they are consistent across the redshift bins. This secondary peak is due to a known cosmological parameter degeneracy (Percival et al. 2001).

### 6.2.4 Distribution of residuals

For a data set with uncorrelated errors (diagonal covariance matrix), the vector of normalized residuals is given by

$$R = \Xi^{-1}(D - T(\theta)) \qquad (66)$$

where Ξ is a diagonal matrix containing the square roots of the variances, D is the data vector, and T(θ) is the theory vector for a given parameter vector θ. If T(θ) represents the true model, and the true errors are known, the residuals are by definition given by a Gaussian with μ = 0 and σ = 1 (Andrae, Schulze-Hartung & Melchior 2010). On the other hand, if the errors are estimated from the data, the residuals are given by Student's t-distribution. This distribution approaches a Gaussian with an increasing number of data points. If this distribution shows a significant deviation from a Gaussian, the model is ruled out. If it follows a Gaussian distribution, either we found the true model, or the current data are not enough to distinguish between the model we found and the true model (Andrae et al. 2010).

When the covariance matrix is not diagonal (the errors are correlated), equation (66) is no longer true and we have to deal with the full covariance matrix. In order to get back to a diagonal matrix, we write the covariance matrix in terms of its eigen-decomposition:

$$C = Q \Lambda Q^{-1} \qquad (67)$$

where Q is the matrix of eigenvectors and Λ is the diagonal matrix containing the eigenvalues of C. The inverse is then given by $C^{-1} = Q \Lambda^{-1} Q^{-1}$, which transforms the $\chi^2$ into:

$$\chi^2 = [\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{th}(\Theta, \nu)]^T Q \Lambda^{-1} Q^{-1} [\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{th}(\Theta, \nu)]. \qquad (68)$$

If we treat Λ as the new (diagonal) covariance matrix, it follows that the normalized residuals are now given by

$$R = \Xi^{-1} Q^{-1} [\hat{S}_{\Delta\ell} - S_{\Delta\ell}^{th}(\Theta, \nu)] \qquad (69)$$

where Ξ now contains the square roots of the eigenvalues. We use equation (69) to calculate the residuals at our best-fitting point in a flat ΛCDM cosmology and plot the results in a histogram (Fig. 12). There are no significant deviations from a Gaussian with μ = 0 and σ = 1, which means that the model seems to be a very good representation of the data.

## 6.3 External data

We compared and combined our results with results obtained from the Planck satellite CMB experiment (Planck Collaboration I 2016a) and Type Ia Supernovae from the Joint Light curve Analysis (JLA)
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Consistency checks for single tomographic redshift bins for all LOWZ and CMASS bins. Here, we show the $\Omega_b \times \Omega_m$ contours taken from a ΛCDM cosmological inference, i.e. varying the same cosmological parameters as the ones from Section 6.4. Panel a shows LOWZ-0, LOWZ-1, LOWZ-2, and LOWZ-3 tomographic bins; panel b shows LOWZ-3, LOWZ-4, LOWZ-5, and CMASS-6 tomographic bins; panel c shows CMASS-6, CMASS-7, CMASS-8, and CMASS-9 tomographic bins; and panel d shows CMASS-9, CMASS-10, CMASS-11, and CMASS-12 tomographic bins.</figcaption></figure>

collaboration (Betoule et al. 2014). The relevant likelihood codes for these probes were implemented and tested in the UCLPI pipeline. We checked that the official cosmological results from the relevant collaborations were recovered in order to use them.

The CMB data from Planck was added through the Planck likelihood codes Commander and Plik (Planck Collaboration XI 2016b). For low multipoles, in the range $l = 2–29$, Commander is used with temperature (TT) and polarization auto- and cross-power spectra (BB, TB, EB). For high multipoles, in the range $l = 30–2508$, Plik is used with temperature (TT) and polarization auto- and cross-power spectra (TE, EE). This configuration is commonly referred to as Planck TT,TE,EE + lowP. Plik also introduces 94 additional nuisance parameters. In order to reduce this large parameter space, we use the lite version of the data offered by the Planck Collaboration. This lite version allows us to compute

a nuisance marginalized CMB likelihood. The only CMB nuisance parameter left is the Planck absolute calibration parameter ($y_{cal}$). We sample this parameter in all the runs that include Planck data. The Planck likelihood codes were added to UCLPI and all the Planck results presented have been obtained using this pipeline. We show in the cosmological contours and in Table 4 that this pipeline reproduces the cosmological results found by Planck Collaboration XIII (2016c).

The SN data from JLA was added to the UCLPI pipeline through the likelihood code provided by the JLA Collaboration (Betoule et al. 2014). This likelihood code needs the luminosity distances to the 740 Supernovae in the sample and four nuisance parameters ($\alpha$, $\beta$, $M_B$, and $\Delta_M$) described in Betoule et al. (2014). The luminosity distances are calculated by CLASS (Blas et al. 2011) using the redshifts of the supernovae within a given cosmology (set by
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison between the histogram of the distribution of residuals (given by equation 69) calculated at the best-fitting point in a flat ΛCDM cosmology (see Section 6.4) and a Gaussian with $μ = 0$ and $σ = 1$. As this distribution does not show any significant deviation from the Gaussian, the model is either the truth or the current data cannot make any further distinction between the model and the truth.</figcaption></figure>

the sampled cosmological parameters). We sample the absolute magnitude at peak brightness ($M_B$) as part of our analysis, and keep the other three nuisance parameters fixed to their best-fitting values found by Betoule et al. (2014) since these have a small impact

in the cosmological parameters when combined with the BOSS data set.

We also implemented a BAO likelihood in order to compare our results with the official *BOSS Consensus BAO* results from Alam et al. (2016). This measurement uses three redshift bins with $z_{eff} = 0.38, 0.51, \text{and } 0.61$, where we used the full shape post-reconstruction measurements from the correlation function and 3D power spectra, which contain additional information from measurements of $f\sigma_8$ (see table 7 from Alam et al. 2016). We have not combined these results with our BOSS results, but we plot them alongside our results alone in the next section to show how our results compare with the BOSS alone results from Alam et al. (2016).

## 6.4 Flat ΛCDM constraints

We obtain constraints for the standard cosmological model, a flat ΛCDM cosmology. We fixed the curvature of the universe to zero, e.g. $\Omega_k = 0$, and varied five cosmological parameters: the baryonic density, $\Omega_b$; the dark matter density, $\Omega_{cdm}$; the amplitude of the primordial power spectra, $A_s$; the spectral index, $n_s$; and the Hubble constant, $h$. As this model considers dark energy as the cosmological constant $\Lambda$, we fixed the $w_0$ parameter to a cosmological constant ($w_0 = -1$); therefore: $\Omega_\Lambda = 1 - (\Omega_b + \Omega_{cdm})$. Here, we also fixed the sum of neutrino masses to the minimum found from neutrino oscillation experiments, $\sum m_\nu = 0.06$ eV (Lesgourgues & Pastor

<table><caption>**Table 4.** Marginalized cosmological constraints and 68% credible intervals for the models considered in this work using a variety of data sets and combinations. The contours for these results are shown in Figs 14 for ΛCDM, 18 for wCDM, 19 for the ΛCDM + $\sum m_\nu$ with $\ell_{\text{max}}^{5 \text{ per cent}}$ cut, and 22 for ΛCDM + $\sum m_\nu$ with $\ell_{\text{max}}^{10 \text{ per cent}}$ cut.</caption><thead><tr><th>Model</th><th>Parameter</th><th>BOSS</th><th>BOSS + JLA</th><th>BOSS + JLA + Planck</th><th>Planck</th></tr></thead><tbody><tr><td rowspan="3">ΛCDM</td><td>Ω<sub>m</sub></td><td>0.315<sup>+0.034</sup><sub>-0.033</sub></td><td>0.317<sup>+0.022</sup><sub>-0.021</sub></td><td>0.327 ± 0.008</td><td>0.315 ± 0.011</td></tr><tr><td>Ω<sub>b</sub></td><td>0.0404<sup>+0.010</sup><sub>-0.009</sub></td><td>0.0381<sup>+0.007</sup><sub>-0.008</sub></td><td>0.0502 ± 0.0006</td><td>0.0492 ± 0.0009</td></tr><tr><td>S<sub>8</sub></td><td>0.715<sup>+0.072</sup><sub>-0.064</sub></td><td>0.745<sup>+0.059</sup><sub>-0.052</sub></td><td>0.862<sup>+0.015</sup><sub>-0.016</sub></td><td>0.850<sup>+0.023</sup><sub>-0.021</sub></td></tr><tr><td rowspan="5">wCDM</td><td>h</td><td>0.716<sup>+0.088</sup><sub>-0.069</sub></td><td>0.699 ± 0.039</td><td>0.663 ± 0.005</td><td>0.672 ± 0.008</td></tr><tr><td>n<sub>s</sub></td><td>0.929<sup>+0.064</sup><sub>-0.045</sub></td><td>0.955<sup>+0.052</sup><sub>-0.048</sub></td><td>0.960 ± 0.004</td><td>0.964 ± 0.006</td></tr><tr><td>Ω<sub>m</sub></td><td>0.277<sup>+0.050</sup><sub>-0.042</sub></td><td>0.308<sup>+0.021</sup><sub>-0.018</sub></td><td>0.330 ± 0.012</td><td>0.213<sup>+0.062</sup><sub>-0.039</sub></td></tr><tr><td>Ω<sub>b</sub></td><td>0.0318<sup>+0.0117</sup><sub>-0.0098</sub></td><td>0.0429 ± 0.007</td><td>0.0505 ± 0.002</td><td>0.0334<sup>+0.009</sup><sub>-0.006</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.726<sup>+0.072</sup><sub>-0.061</sub></td><td>0.743<sup>+0.079</sup><sub>-0.068</sub></td><td>0.863 ± 0.016</td><td>0.811<sup>+0.037</sup><sub>-0.034</sub></td></tr><tr><td rowspan="4">wCDM</td><td>h</td><td>0.767<sup>+0.069</sup><sub>-0.091</sub></td><td>0.745<sup>+0.049</sup><sub>-0.052</sub></td><td>0.661 ± 0.012</td><td>0.816<sup>+0.073</sup><sub>-0.101</sub></td></tr><tr><td>n<sub>s</sub></td><td>0.939<sup>+0.057</sup><sub>-0.049</sub></td><td>0.957<sup>+0.049</sup><sub>-0.050</sub></td><td>0.960 ± 0.004</td><td>0.964 ± 0.006</td></tr><tr><td>w<sub>0</sub></td><td>-1.36<sup>+0.36</sup><sub>-0.38</sub></td><td>-1.030<sup>+0.073</sup><sub>-0.076</sub></td><td>-0.993<sup>+0.046</sup><sub>-0.043</sub></td><td>-1.45<sup>+0.32</sup><sub>-0.23</sub></td></tr><tr><td>Ω<sub>m</sub></td><td>0.326<sup>+0.038</sup><sub>-0.035</sub></td><td>0.304<sup>+0.022</sup><sub>-0.021</sub></td><td>0.328 ± 0.009</td><td>0.326<sup>+0.028</sup><sub>-0.021</sub></td></tr><tr><td rowspan="5">ΛCDM + Σm<sub>v</sub> [ℓ<sup>5</sup><sub>max</sub> cut]</td><td>Ω<sub>b</sub></td><td>0.040<sup>+0.009</sup><sub>-0.010</sub></td><td>0.0432 ± 0.008</td><td>0.05017<sup>+0.0009</sup><sub>-0.0008</sub></td><td>0.0506<sup>+0.0039</sup><sub>-0.0026</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.723<sup>+0.069</sup><sub>-0.063</sub></td><td>0.700<sup>+0.065</sup><sub>-0.056</sub></td><td>0.862 ± 0.017</td><td>0.836<sup>+0.031</sup><sub>-0.035</sub></td></tr><tr><td>h</td><td>0.730<sup>+0.075</sup><sub>-0.078</sub></td><td>0.814<sup>+0.054</sup><sub>-0.064</sub></td><td>0.663<sup>+0.006</sup><sub>-0.007</sub></td><td>0.662<sup>+0.018</sup><sub>-0.026</sub></td></tr><tr><td>n<sub>s</sub></td><td>0.933<sup>+0.066</sup><sub>-0.046</sub></td><td>0.941<sup>+0.055</sup><sub>-0.049</sub></td><td>0.960 ± 0.042</td><td>0.962<sup>+0.006</sup><sub>-0.007</sub></td></tr><tr><td>Σm<sub>v</sub> (95% CI)[eV]</td><td>&lt;0.75</td><td>&lt;0.71</td><td>&lt;0.14</td><td>&lt;0.76</td></tr><tr><td rowspan="5">ΛCDM + Σm<sub>v</sub> [ℓ<sup>10</sup><sub>max</sub> cut]</td><td>Ω<sub>m</sub></td><td>0.345<sup>+0.033</sup><sub>-0.030</sub></td><td>0.324<sup>+0.034</sup><sub>-0.029</sub></td><td>0.333<sup>+0.014</sup><sub>-0.012</sub></td><td>0.326<sup>+0.050</sup><sub>-0.029</sub></td></tr><tr><td>Ω<sub>b</sub></td><td>0.045 ± 0.009</td><td>0.040 ± 0.013</td><td>0.0510<sup>+0.0016</sup><sub>-0.0014</sub></td><td>0.0506<sup>+0.0069</sup><sub>-0.0033</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.751<sup>+0.062</sup><sub>-0.057</sub></td><td>0.768<sup>+0.097</sup><sub>-0.092</sub></td><td>0.864<sup>+0.030</sup><sub>-0.029</sub></td><td>0.839<sup>+0.058</sup><sub>-0.067</sub></td></tr><tr><td>h</td><td>0.689<sup>+0.076</sup><sub>-0.066</sub></td><td>0.661<sup>+0.067</sup><sub>-0.063</sub></td><td>0.658<sup>+0.010</sup><sub>-0.011</sub></td><td>0.662<sup>+0.024</sup><sub>-0.044</sub></td></tr><tr><td>n<sub>s</sub></td><td>0.930<sup>+0.062</sup><sub>-0.044</sub></td><td>1.011<sup>+0.056</sup><sub>-0.086</sub></td><td>0.958 ± 0.006</td><td>0.962 ± 0.013</td></tr><tr><td></td><td>Σm<sub>v</sub> (95% CI)[eV]</td><td>&lt;0.72</td><td>&lt;0.66</td><td>&lt;0.16</td><td>&lt;0.76</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** 2D $\Omega_m \times S_8$ marginalized credible intervals for a $\Lambda$CDM cosmology. In this figure we show in detail the cosmological results from Section 6.4 for BOSS $C_\ell$s only (blue); BOSS $C_\ell$s plus JLA (green); BOSS $C_\ell$s plus JLA and Planck (purple); together with results using the post-reconstruction full shape (incl. $\sigma_8(z)$) from Alam et al. (2016) consensus results (pink, dot-dashed), and Planck alone (red). In order to compare our results to a weak-lensing probe, we also show the results from Hildebrandt et al. (2017) (grey, dashed). For details on the external data sets, see Section 6.3.</figcaption></figure>

2006, 2014). From these, we also obtained derived parameters: the total matter density, $\Omega_m \equiv \Omega_b + \Omega_{cdm}$, and the fluctuation of amplitude at 8 $h^{-1}$Mpc, $\sigma_8$ or $S_8 = \sigma_8\sqrt{\Omega_m/0.3}$. Finally, as described in the previous section, we also varied the BOSS, Planck, and JLA nuisance parameters. For this analysis, we used the $\ell_{\text{max}}^5$ per cent cuts (see table 2).

We checked the consistency of the data sets by running the same analysis for these probes alone and combined (Planck, JLA, and Planck plus JLA) and calculating respective Bayes factors for these runs. The Bayes factor, equation (65), for combinations of the considered data sets indicates consistency between all three probes:

$$R_{\text{BOSS+JLA}}^{\Lambda\text{CDM}} \approx 18 \qquad (70)$$

$$R_{\text{BOSS+PLANCK}}^{\Lambda\text{CDM}} \approx 74 \qquad (71)$$

$$R_{\text{PLANCK+JLA}}^{\Lambda\text{CDM}} \approx 11 \qquad (72)$$

$$R_{\text{BOSS+PLANCK+JLA}}^{\Lambda\text{CDM}} \approx 4 \times 10^4; \qquad (73)$$

these indicate that the data sets are compatible for the considered model, given the chosen priors.

Finally, when considering the combination of BOSS $C_\ell$s, Planck, and JLA, we find results consistent with Alam et al. (2016), Hildebrandt et al. (2017), and DES Collaboration (2017). Fig. 13 shows the $\Omega_m \times S_8$ 2D plane for this analysis and comparisons with Planck and the BOSS full-shape post-reconstruction from Alam et al. (2016). Despite the Bayes factors showing no significant reason to be concerned about the compatibility of these data sets,

we see an interesting trend in this figure in so far as a tension and the BOSS data set in this paper preferring a smaller $S_8$ than the Planck analysis. We argue that this method would prove potentially very useful in resolving any $S_8$ tensions that exist currently between CMB and weak-lensing data (MacCrann et al. 2015; Charnock, Battye & Moss 2017b).

The results for the 1D marginalized cosmological constraints for BOSS $C_\ell$ and combinations, together with the 68 per cent credible intervals, can be found in Table 4. The $1\sigma$ and $2\sigma$ contour levels can be found in Fig. 14, where the nuisance parameters have been marginalized over. Fig. 15 shows the best-fitting theory $C_\ell$ using the parameters estimated from this analysis with a $\chi^2_{\text{red}} = 1.08$, which also indicates reliability and robustness of the analysis performed.

Even though we do not show the results in this work, we performed a cosmological analysis using a $\Lambda$CDM with a fixed zero neutrino mass, $\sum m_\nu = 0$ eV. We compared it with the model used in this section, $\Lambda$CDM with $\sum m_\nu$ fixed to 0.06 eV, using the Bayes factor for model selection. Consider $\tilde{D}$ representing the combination of data vectors; the Bayes factor is given by

$$R_{A,B} = \frac{P(\tilde{D}_{\text{BOSS+PLANCK+JLA}} | \sum m_\nu = 0.06 \text{ eV})}{P(\tilde{D}_{\text{BOSS+PLANCK+JLA}} | \sum m_\nu = 0 \text{ eV})} = 8 \times 10^5 \qquad (74)$$

This indicates that, for a $\Lambda$CDM model, the data prefers massive neutrinos over no neutrino mass at all.

## 6.5 Flat wCDM constraints

In this section, we allowed the equation of state of dark energy, $w_0$, to vary. This is a trivial extension of the standard model of cosmology with just one extra parameter. If $w_0 = -1$, the solution indicates that the nature of dark energy is actually the cosmological constant, $\Lambda$. The procedure for this analysis followed in a similar fashion as the one outlined in Section 6.4, varying six parameters instead of five: $\Omega_b, \Omega_{cdm}, n_s, \ln 10^8 A_s, h$, and the extra $w_0$. Note that, for this case, we are not varying $w_a$, i.e. we do not consider a redshift evolution in the equation of state of dark energy. Again, we fixed the neutrino parameter to $\sum m_\nu = 0.06$ eV (Lesgourgues & Pastor 2006, 2014). Here, we used the same $\ell_{\text{max}}^5$ per cent cuts as in the last section (see Table 2).

Figs 16 and 17 show in detail the contours for $S_8 \times \Omega_m$ and $w_0 \times \Omega_m$, respectively, and comparisons with previous measurements in the literature. From Fig. 17 and from the complete set of results in 18, we show that an ~4 per cent error ($1\sigma$ CI) on the equation of state of dark energy is obtained from this cosmological analysis:

$$w_0 = -0.993^{+0.046}_{-0.043}. \qquad (75)$$

This result is consistent with the $\Lambda$CDM scenario of standard cosmology, i.e. it is consistent with dark energy being a cosmological constant, $\Lambda$. Note from Fig. 18 that we find a small value of $h$ (compared to Planck Collaboration XIII 2016c) when combining BOSS $C_\ell$s, Planck, and JLA:

$$h^{\Lambda\text{CDM}} = 0.661 \pm 0.012. \qquad (76)$$

This value is lower than the quoted Planck value alone; this puts further tension in these measurements if compared to the Hubble constant result from Cepheid Variables (Riess et al. 2016, 2018).

As the model in this section is different from the previous section, we performed an evidence analysis using the Bayes factor, equation (65), in order to be sure that our measurements can be combined with the external data described in Section 6.3. The following measurements indicate that such combinations are
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Marginalized 1D and 2D cosmological constraints for the ΛCDM model varying five cosmological parameters with 1σ (darker) and 2σ (lighter) contour levels. We show here a combination of sampled and relevant derived parameters: Ω<sub>m</sub>, Ω<sub>b</sub>, S<sub>8</sub>, h, and n<sub>s</sub> (marginalizing over τ<sub>reio</sub> for the Planck combinations). The blue contours were estimated from the BOSS C<sub>ℓ</sub>s data alone using the SH16 likelihood; the green contours are a combination the BOSS likelihood and JLA data (see Section 6.3); the red contours are the Planck high-ℓ TT, TE, EE, and low-ℓ P likelihood results (see Section 6.3); finally, the purple contours are a combination of the three probes: BOSS C<sub>ℓ</sub>, JLA, and Planck (also high-ℓ TT, TE, EE, and low-ℓ P). Note that none of the results here use Planck Lensing data.</figcaption></figure>

consistent for a wCDM model:

$$
R_{\text{BOSS+JLA}}^{w\text{CDM}} \approx 2 \times 10^2 \qquad (77)
$$

$$
R_{\text{BOSS+PLANCK}}^{w\text{CDM}} \approx 4 \times 10^3 \qquad (78)
$$

$$
R_{\text{PLANCK+JLA}}^{w\text{CDM}} \approx 2 \qquad (79)
$$

$$
R_{\text{BOSS+PLANCK+JLA}}^{w\text{CDM}} \approx 3 \times 10^5 . \qquad (80)
$$

Finally, we used the ratio of the evidences, the Bayes factor, to perform a model selection between wCDM and ΛCDM using the final data set combination. Assuming $\tilde{D}$ to be the combination of data vectors for all the data sets, the Bayes factor between the two models is

$$
R_{w, \Lambda} = \frac{P(\tilde{D}_{\text{BOSS+Planck+JLA}} | w\text{CDM})}{P(\tilde{D}_{\text{BOSS+Planck+JLA}} | \Lambda\text{CDM})} = 0.67 \qquad (81)
$$

6.6 Flat ΛCDM + ∑m<sub>v</sub> constraints

For the last model considered in this work, we assume a flat ΛCDM with variable neutrino masses, varying the sum of the species' masses, $\sum m_v$. In the previous sections, we fixed the sum of neutrino masses to $\sum m_v = 0.06$ eV due to results from neutrino oscillation experiments for the lower bound of the normal neutrino mass ordering (Hannestad 2003; Lesgourgues & Pastor 2006; Hannestad & Schwetz 2016).

In this section, we considered one of the two different scenarios regarding different neutrino hierarchies, the normal hierarchy. To approximate the normal hierarchy, one can approximate the two lower masses to be zero and vary $\sum m_v$ for one remaining massive species. We do not investigate details of how the prior on the hierarchy or on the absolute mass changes this result and we leave this to a future analysis. We fix $N_{\text{eff}} = 3.046$ by changing the values of massive neutrinos and ultra-relativist particles for the case considered, i.e. $N_v = 1$ and $N_{\text{ur}} = 2.0328$.

We perform an analysis using the same ℓ-range as in the previous sections, $\ell_{\text{max}}^{\text{per cent}}$ from Table 2. A summary of the marginalized 1D
Total $\chi^2_{red} = 1.08319$

N. Data Points: 441, Numb Params: 33

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Auto- and cross-angular power spectra for the 13 tomographic redshift bins considered for the BOSS DR12 samples: LOWZ (sample 1) and CMASS (sample 2). The shaded blue regions show the scales considered in the cosmological parameter estimation in Section 6. The data points are the pseudo-$C_\ell$ estimates, described in Section 3, for LOWZ and CMASS. The solid blue lines, generated with UCLCL, reflect the best-fitting auto- and cross-power spectra for the $\Lambda$CDM model estimated in Section 6.4. Finally, the black dashed lines show both shot noise and sampled shot noise (for bins 11 and 12). The overall reduced $\chi^2$ for this fit is $\chi^2_{red} \approx 1.08$, where the number of data points is 441 and the total number of sampled parameters is 33 – 5 cosmological parameters and 28 nuisance parameters. The title on each individual plot reflects the bins i & j for each $C_\ell^{ij}$, the $\chi^2$ per data point ($\chi^2/N_d$), and the number of data points for that individual angular power spectrum, $N_d$. The $\ell$-ranges used in this figure correspond to $\ell_{max}^5$ per cent in Table 2. Most of the constraining power comes from the auto-power spectra. The cross-power spectra serve to constrain parameters related to the RSD by helping to break the degeneracy between the bias and $A_s$ while also probing the redshift dispersion due to the peculiar motion of galaxies (FoG).</figcaption></figure>

credible intervals from the cosmological estimation made with this cut can be found in the third part of Table 4 showing the $1\sigma$ intervals for the $\Lambda$CDM parameters plus the 95 per cent upper limit for $\sum m_\nu$. The 1D and 2D marginalized credible intervals for this analysis can be found in Fig. 19. When considering an approximation for the

normal hierarchy, for a combination of BOSS $C_\ell$, Planck CMB data, and supernovae data from JLA, the 95 per cent upper limit for sum of neutrino masses is

$$ \sum m_\nu < 0.14 \text{ eV (BOSS + Planck + JLA} - \ell_{\text{max}}^5 \text{ per cent cut).} \quad (82) $$
<figure><img src="image_17.png" /><figcaption>**Figure 16.** 2D $\Omega_m \times S_8$ marginalized credible intervals for a wCDM Cosmology. This shows in detail the cosmological results from Section 6.5 for BOSS Cℓs only (blue); BOSS Cℓs plus JLA (green); BOSS Cℓs plus JLA and Planck (purple); together with results using just the full shape (pre-reconstruction) from Alam et al. (2016) consensus results (pink), and Planck alone (red).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** 2D $w_0 \times \Omega_m$ marginalized credible intervals for a wCDM Cosmology. This shows in detail the cosmological results from Section 6.5 for BOSS Cℓs only (blue); BOSS Cℓs plus JLA (green); BOSS Cℓs plus JLA and Planck (purple); together with results using just the full shape (pre-reconstruction) from Alam et al. (2016) consensus results (pink), JLA (yellow), and Planck alone (red).</figcaption></figure>

From Fig. 20 and even more so from Fig. 21, one can notice that we are not so far from excluding zero total neutrino mass using cosmological data alone. As the power of such data sets increase we should be able, using the correct analysis and tools, to measure and detect neutrino masses independently from atmospheric experiments.

We then proceed to check the consistency of data sets by using the evidence of each cosmological parameter estimation for these models to calculate the Bayes factor (equation 65):

$$R_{\text{BOSS+JLA}}^{\Delta \text{CDM} + \sum m_v \text{ 5 per cent}} \approx 1 \times 10^2 \quad (83)$$

$$R_{\text{BOSS+PLANCK}}^{\Delta \text{CDM} + \sum m_v \text{ 5 per cent}} \approx 4 \times 10^2 \quad (84)$$

$$R_{\text{PLANCK+JLA}}^{\Delta \text{CDM} + \sum m_v} \approx 40 \quad (85)$$

$$R_{\text{BOSS+PLANCK+JLA}}^{\Delta \text{CDM} + \sum m_v \text{ 5 per cent}} \approx 3 \times 10^2. \quad (86)$$

These values, again, indicate the consistency of data sets for the considered model.

For the final analysis in this work, we extended the scales considered for the $\ell_{\text{max}}^{10 \text{ per cent}}$ cuts (see Table 2 for details). This allows us to access smaller scales that are still in the beginning of the so-called the weak non-linear regime (Thomas et al. 2010; Bird et al. 2012). Note that these scales are still larger than the scales that most collaborations use for power spectra or correlation function cosmological analysis (Ho et al. 2012; Alam et al. 2016; DES Collaboration 2017; Idebrandt et al. 2017) – DES Collaboration (2017), for example, uses scales up to $0.78 h \text{ Mpc}^{-1}$. In other words, one can be confident that the $\ell_{\text{max}}^{10 \text{ per cent}}$ cuts are safe to be used, not using scales outside the weak non-linear regime.

We then proceed to perform a similar cosmological analysis for a $\Delta \text{CDM}$ model with one massive species of neutrino, approximating the normal hierarchy. The 1D and 2D marginalized credible intervals for these final analyses can be found in Fig. 22 and the marginalized 68 per cent credible intervals for the $\Delta \text{CDM}$ parameters and the 95 per cent credible interval upper limit for $\sum m_v$ using this cut can be found in Table 4. The Bayes factors for this choice are shown below (note that we have two further nuisance parameters in the $\ell_{\text{max}}^{10 \text{ per cent}}$ cut and we checked that failure to add these reduces the Bayes factor significantly).

$$R_{\text{BOSS+JLA}}^{\Delta \text{CDM} + \sum m_v \text{ 10 per cent}} \approx 70 \quad (87)$$

$$R_{\text{BOSS+PLANCK}}^{\Delta \text{CDM} + \sum m_v \text{ 10 per cent}} \approx 6 \times 10^2 \quad (88)$$

$$R_{\text{BOSS+PLANCK+JLA}}^{\Delta \text{CDM} + \sum m_v \text{ 10 per cent}} \approx 3 \times 10^5. \quad (89)$$

This extended scale analysis demonstrates the robustness of the results presented in this section as the 95 per cent CI upper bound for $\sum m_v$ remains robust to these cuts (see Figs 20 and 21):

$$\sum m_v < 0.16 \text{ eV (BOSS + Planck + JLA} - \ell_{\text{max}}^{10 \text{ per cent}} \text{ cut)}. \quad (90)$$
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Cosmological constraints for the wCDM model varying now six cosmological parameters. This figure contains a combination of sampled and relevant derived parameters from the chains: $\Omega_m, \Omega_b, S_8, h, n_s$, and $w_0$. Note that the Planck chains also varied $\tau_{\text{reio}}$. The blue contours were estimated from the BOSS C$_l$ data alone using the SH16 likelihood; green contours are a combination of the BOSS likelihood and JLA data; red contours are the Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P likelihood results; finally, the purple contours are a combination of the three probes: BOSS C$_l$, JLA, and Planck (also high-$\ell$ TT, TE, EE, and low-$\ell$ P). The apparent cuts in the Planck alone contours are due to the prior in $h$. Note, again, that none of the results here use Planck Lensing data.</figcaption></figure>

## 7 CONCLUSIONS

In this work, we have taken a different approach<fnref n="13" /> obtaining galaxy clustering information from the BOSS DR12 large-scale structure catalogue (Reid et al. 2016). This approach consisted of using a pseudo angular power spectra estimator (PCL) applied to 13 tomographic redshift bins ranging from $0.15 \le z < 0.8$ with a redshift-dependent bias, a redshift dispersion, and extra shot-noise as nuisance parameters to be sampled with the cosmological parameters using UCLPI (Cuceu et al., in preparation). In this approach, we also used splines of the data as input for the simulation used for covariance matrix estimation.

The tomographic approach in redshift space and the covariance matrix estimation method used in this work allowed us to perform a cosmology-free inference from the data. In other words, nowhere in this analysis was a fiducial cosmology assumed. This is, by itself, a great advantage over methods that use $P(k)$ or $\xi(r)$ as

these need to assume a fiducial cosmology in order to transform from redshift space to radial distances. The impact of such strong assumption in the cosmological inference is still unknown. We performed systematic and consistency checks with the data and the method itself with satisfactory results. From the 18 different sources of systematics considered in Section 5, none demonstrated worrying excess of power in the scales considered in Table 2. Consistency checks demonstrated the robustness of our estimated covariance matrix, since the recovered cosmology was the same under different estimation methods (through simulation and theory); the likelihood, given that it returns the same contours under three different approaches; and of our whole method, since we recovered the right cosmology from a controlled simulation.

Cosmological parameters were obtained for three different models: $\Lambda$CDM, wCDM, and $\Lambda$CDM with $\sum m_\nu$. We highlight the following main points regarding the results obtained in Section 6:

(1) The constraints obtained for all three models considered, using a tomographic analysis in harmonic space, are extremely competitive in comparison to the ones obtained by the BOSS Collaboration (Alam et al. 2016) and other recent large collabo-

<footnote n="13">¹³Compared to the approaches from the official BOSS Collaboration papers: Ross et al. (2017); Beutler et al. (2017a,b); Satpathy et al. (2017); Sánchez et al. (2017); Grieb et al. (2017); Salazar-Albornoz et al. (2017); Wang et al. (2017); and Zhao et al. (2017).</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 19.** 1D and 2D marginalized credible intervals for a ΛCDM Cosmology with $\sum m_\nu$ when using scales up to $k_{\text{max}} \approx 0.07h \text{ Mpc}^{-1}$ ($\ell_{\text{max}}^{5 \text{ per cent}}$ cut). Here we show the $\Omega_m, \Omega_b, S_8, h, n_s$, and $\sum m_\nu$ contours for BOSS $C_\ell$s alone (blue); BOSS $C_\ell$s plus JLA (green); Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P (red); and BOSS $C_\ell$s plus JLA and Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P (purple). As most of the scales that contain clean information on the neutrino masses are cut off, the 95 per cent CI upper bound found is $\sum m_\nu < 0.14 \text{ eV}$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** 1D marginalized with 95 per cent credible intervals for $\sum m_\nu$ in three different cases: (red solid) Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P (yellow dashed) BOSS $C_\ell$s with the $\ell_{\text{max}}^{10 \text{ per cent}}$ cut plus Planck and JLA, and (purple dotted) BOSS $C_\ell$s with the $\ell_{\text{max}}^{5 \text{ per cent}}$ cut plus Planck and JLA. The 95 per cent upper limit for each case is, respectively: (red) 0.76 eV, (yellow) 0.16 eV, and (purple) 0.14 eV.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21.** 2D marginalized $1\sigma$ and $2\sigma$ credible intervals for the $\sum m_\nu - \Omega_m$ plane for three different cases: (blue solid) Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P, (red dashed) BOSS $C_\ell$s with the $\ell_{\text{max}}^{10 \text{ per cent}}$ cut plus Planck and JLA, and (orange dotted) BOSS $C_\ell$s with the $\ell_{\text{max}}^{5 \text{ per cent}}$ cut plus Planck and JLA.</figcaption></figure>

ration results such as DES (DES Collaboration 2017) and KiDS (Hildebrandt et al. 2017) with errors as small as those obtained by these collaborations.

(2) Even though information along the line-of-sight is 'washed away' due to projecting the data into tomographic bins, we obtain one of the tightest constraints for the equation-of-state of dark energy with an ~4 per cent error when combining BOSS $C_\ell$s, Planck CMB, and JLA Supernovae. This has not been achieved before using $C_\ell$ with a spectroscopic survey and the constraint is as tight as the one obtained from the state-of-the-art Dark Energy Collaboration analysis, using a combination of DES galaxy clustering & weak lensing, Planck, JLA, and BAO (DES Collaboration 2017).

(3) For the models and data sets considered, we find very high values for the Bayes factor, $R$, when combining BOSS $C_\ell$s, Planck, and JLA. We highlight: $R_{\text{BOSS+PLANCK+JLA}} \simeq 4 \times 10^4$ for ΛCDM and $R_{\text{BOSS-10 per cent+PLANCK+JLA}} \simeq 3 \times 10^5$ for ΛCDM varying neutrino masses.

(4) The Bayes factor can also be used for model selection. Considering the combination of data sets, the Bayes factor between ΛCDM and wCDM is

$$R_{w,\Lambda} = \frac{P(\bar{D}_{\text{BOSS+Planck+JLA}} | w\text{CDM})}{P(\bar{D}_{\text{BOSS+Planck+JLA}} | \Lambda\text{CDM})} = 0.67,$$

where $\bar{D}$ here represents the overall combination of data vectors. This indicates that this combination prefers slightly ΛCDM to wCDM, although no strong conclusion can be made.

(5) We find a small tension between BOSS $C_\ell$s and Planck for $S_8$ in all models considered, with BOSS preferring smaller values. For example, for ΛCDM:

$$S_8 = 0.715^{+0.072}_{-0.064} \quad (\text{BOSS})$$

$$S_8 = 0.850^{+0.023}_{-0.021} \quad (\text{Planck})$$
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Cosmological constraints for the $\Lambda$CDM + $\sum m_\nu$ model, using the $\ell_{\text{max}}^{10 \text{ per cent}}$ cut, varying now six cosmological parameters, including the sum of neutrino masses considering only one massive species. This figure contains a combination of sampled and relevant derived parameters from the chains: $\Omega_m$, $\Omega_b$, $S_8$, $h$, $n_s$, and $\sum m_\nu$. Note that the Planck chains also varied $\tau_{\text{reio}}$. The blue contours were estimated from the BOSS $C_\ell$ data alone using the SH16 likelihood; the green contours are a combination of the BOSS likelihood and JLA data; the red contours are the Planck high-$\ell$ TT, TE, EE, and low-$\ell$ P likelihood results; finally, the purple contours are a combination of the three probes: BOSS $C_\ell$, JLA, and Planck (also high-$\ell$ TT, TE, EE and low-$\ell$ P). For this scale cut, the combination of data sets yields an upper bound on $\sum m_\nu < 0.16$ eV.</figcaption></figure>

although the combination of these data sets prefers higher values such as Planck (see Table 4) and the Bayes factor suggest the data sets are compatible. We conclude that such tension can be investigated further with this method as LSS data increase in size and depth.

(6) Even though we do not show these results, we performed a cosmological analysis using a $\Lambda$CDM model but fixing $\sum m_\nu = 0$eV and compared with the $\Lambda$CDM results from Section 6.4, which has a $\sum m_\nu$ fixed to 0.06 eV. Using the Bayes factor for model selection, it is clear that the data prefer massive neutrinos against no neutrino mass at all:

$$R_{0.06\text{eV},0} = \frac{P(\vec{D}_{\text{BOSS+Planck+JLA}} | \Lambda\text{CDM} + \sum m_\nu = 0.06)}{P(\vec{D}_{\text{BOSS+Planck+JLA}} | \Lambda\text{CDM} + \sum m_\nu = 0)} = 8 \times 10^5.$$

(7) The neutrino mass constraints we obtain here can be compared to the tomographic analysis in real space done by Salazar-

Albornoz et al. (2017), which obtains an upper bound of $\sum m_\nu < 0.474$ eV (95 per cent CI). The reason we obtain much tighter constraints ($\sum m_\nu < 0.14$ eV (95 per cent CI)), even though we are also performing a tomographic analysis, is due to a series of decisions, including the approach we take to model the redshift dispersion and galaxy ‘shell-crossing’ (see Section 4.1.1), bias, and extra-shot noise. It is possible that the main difference between the results is due to different approach in modelling the neutrino mass hierarchy. Salazar-Albornoz et al. (2017) consider a model where the three neutrino species have degenerate mass hierarchy, i.e. the three masses are equal. This is already ruled out by particle physics experiments that measure the mass splitting from neutrino oscillation experiments (see Gonzalez-Garcia, Maltoni & Schwetz 2014 for an update on the neutrino mass splitting fits). The approach we took in this work (see Section 6.6) naturally yields smaller upper bounds in $\sum m_\nu$. In Loureiro et al. (2018), we studied the impact of model choice in the sum of neutrino masses and their hierarchy, showing that this upper bound was driven by our model choices.
We have shown here that it is possible to recover very powerful and competitive cosmological constraints from spectroscopy alone with a 2 + 1D analysis – as good as the constraints from doing a 3D analysis. There are several reasons why a tomographic analysis in harmonic space performs as well as the standard 3D analysis. The precision of spectroscopic redshifts allows for very fine radial binning, probing the evolution of structure in more detail. The ‘effective redshift’ approximation impacts 3D analyses, complicating the proper treatment of ‘light-cone’ effects. By contrast, this is not a problem with $C_ℓ$s as we are dealing with data in redshift space using a fine binning. As we are using data in redshift space, no transformation to comoving coordinates is necessary, and hence no fiducial cosmology is required. We argue that it is much simpler in Fourier space than in real space to correct certain types of systematic errors and to perform scale cuts to deal with non-linearities; the arguments here are similar to the ones used to compare $P(k)$ and $\xi(r)$.

Finally, we point out that the most promising advantage of performing this 2+1D tomographic analysis in harmonic space is the practicality of combining spectroscopic and photometric probes, which include also cosmic shear, as the latter also ‘lives’ in a 2 + 1D space. Using the method proposed by McLeod et al. (2017), the spectroscopic sample has the potential of ‘fixing’ the photometric redshift limitations when probing the photometric clustering redshift distribution together with the cosmological parameters.

With increasingly large future photometric and spectroscopic surveys such as LSST (LSST Dark Energy Science Collaboration 2012), DESI (Levi et al. 2013), J-PAS (Benitez et al. 2014), and Euclid (Laureijs et al. 2011), the future of precision cosmology lies in our ability to combine data sets from across the entire electromagnetic spectrum. The angular power spectrum approach offers a unified framework for coherently combining different data sets in order to obtain maximal information from each (Joachimi & Bridle 2010; Kirk et al. 2015; McLeod et al. 2017). We believe that the approach used in this paper, in which cosmological information is extracted from the projected distribution of the galaxies in a spectroscopic survey, is a useful step towards achieving this unified framework. We further claim that this approach leads to a better understanding of the evolution of structure in the Universe as it provides more information on the redshift evolution of galaxy bias.

## ACKNOWLEDGEMENTS

We would like to acknowledge the BOSS Collaboration for their great work and for making their data public. AL and BM thank the support of the Brazilian people through the Science without Borders CNPq Brazil fellowships. BM and FBA acknowledge support from the European Community through the DEDALE grant (contract no. 665044) within the H2020 Framework Program of the European Commission. FBA acknowledges the Royal Society for support via a Royal Society URF. AC acknowledges the Royal Astronomical Society for support via a summer bursary. OL acknowledges support from an European Research Council Advanced Grant FP7/291329 and support from the UK Science and Technology Research Council (STFC) Grant No. ST/M001334/1. MM acknowledged the support from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No 6655919. We also thank Niall Jeffrey and Andreu Font-Ribera for their very helpful comments on the final draft of this work. All cosmological contour plots were generated using ChainConsumer (Hinton 2016).

## REFERENCES

* Alam S. et al., 2015, ApJS, 219, 12

* Alam S. et al., 2017, MNRAS, 470, 2617

* Anderson L. et al., 2014, MNRAS, 441, 24

* Anderson T. W., 2003, An Introduction to Multivariate Statistical Analysis, 3rd edn. Wiley, New York

* Andrae R., Schulze-Hartung T., Melchior P., 2010, preprint (arXiv:1012.3754)

* Asorey J., Crocce M., Gaztañaga E., Lewis A., 2012, MNRAS, 427, 1891

* Balaguera-Antolinez A., Bilicki M., Branchini E., Postiglione A., 2018, MNRAS, 476, 1050

* Benitez N. et al., 2014, preprint (arXiv:1403.5237)

* Betoule M. et al., 2014, A&A, 568, A22

* Beutler F. et al., 2017a, MNRAS, 464, 3409

* Beutler F. et al., 2017b, MNRAS, 466, 2242

* Bird S., Viel M., Haehnelt M. G., 2012, MNRAS, 420, 2551

* Blake C., Collister A., Bridle S., Lahav O., 2007, MNRAS, 374, 1527

* Blake C., Ferreira P. G., Borrill J., 2004, MNRAS, 351, 923

* Blas D., Lesgourgues J., Tram T., 2011, J. Cosmol Astropart. Phys., 7, 034

* Brown M. L., Castro P. G., Taylor A. N., 2005, MNRAS, 360, 1262

* Cannon R. et al., 2006, MNRAS, 372, 425

* Challinor A., Lewis A., 2011, Phys. Rev. D, 84, 043516

* Charnock T., Battye R. A., Moss A., 2017a, Phys. Rev. D, 95, 123535

* Charnock T., Battye R. A., Moss A., 2017b, Phys. Rev. D, 95, 123535

* Christodoulou L. et al., 2012, MNRAS, 425, 1527

* Dahlen F. A., Simons F. J., 2008, Geophys. J. Int., 174, 774

* Dawson K. S. et al., 2013, AJ, 145, 10

* de Haan T. et al., 2016, ApJ, 832, 95

* DES Collaboration et al., 2018, Phys. Rev. D, 98, 043526

* Di Dio E., Montanari F., Lesgourgues J., Durrer R., 2013, J. Cosmol. Astropart. Phys., 11, 044

* Doux C., Penna-Lima M., Vitenti S. D. P., Tréguer J., Aubourg E., Ganga K., 2018, MNRAS, 480, 5386

* Efstathiou G., 2004, MNRAS, 349, 603

* Efstathiou G., Lemos P., 2018, MNRAS, 476, 151

* Eisenstein D. J. et al., 2001, AJ, 122, 2267

* Elsner F., Leistedt B., Peiris H. V., 2016, MNRAS, 456, 2095

* Elsner F., Leistedt B., Peiris H. V., 2017, MNRAS, 465, 1847

* Elvin-Poole J. et al., 2018, Phys. Rev. D, 98, 042006

* Eriksen M., Gaztañaga E., 2015, MNRAS, 452, 2168

* Esteban I., Gonzalez-Garcia M. C., Maltoni M., Martinez-Soler I., Schwetz T., 2017, J. High Energy Phys., 1, 87

* Feeney S. M., Peiris H. V., Williamson A. R., Nissanke S. M., Mortlock D. J., Alsing J., Scolnic D., 2018, preprint (arXiv:1802.03404)

* Feroz F., Hobson M. P., 2008, MNRAS, 384, 449

* Fisher K. B., Scharf C. A., Lahav O., 1994, MNRAS, 266, 219

* Fixsen D. J., Cheng E. S., Gales J. M., Mather J. C., Shafer R. A., Wright E. L., 1996, ApJ, 473, 576

* Gaztañaga E., Eriksen M., Crocce M., Castander F. J., Fosalba P., Marti P., Miquel R., Cabré A., 2012, MNRAS, 422, 2904

* Giannantonio T. et al., 2016, MNRAS, 456, 3213

* Gonzalez-Garcia M. C., Maltoni M., Schwetz T., 2014, J. High Energy Phys., 11, 52

* Grieb J. N. et al., 2017, MNRAS, 467, 2085

* Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

* Hand N. et al., 2015, Phys. Rev. D, 91, 062001

* Hannestad S., 2003, Phys. Rev. D, 67, 085017

* Hannestad S., Schwetz T., 2016, J. Cosmol. Astropart. Phys., 11, 035

* Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

* Hauser M. G., Peebles P. J. E., 1973, ApJ, 185, 757

* Hildebrandt H. et al., 2017, MNRAS, 465, 1454

* Hinton S. R., 2016, J. Open Source Softw., 1, 00045
Hivon E., Górski K. M., Netterfield C. B., Crill B. P., Prunet S., Hansen F., 2002, ApJ, 567, 2

Ho S. et al., 2012, ApJ, 761, 14

Huterer D., Knox L., Nichol R. C., 2001, ApJ, 555, 547

Joachimi B., Bridle S. L., 2010, A&A, 523, A1

Johansson H. T., Forssén C., 2016, SIAM J. Sci. Comput., 38, 376

Kaiser N., 1987, MNRAS, 227, 1

Kajita T., 2016, Rev. Mod. Phys., 88, 030501

Kang X., Jing Y. P., Mo H. J., Börner G., 2002, MNRAS, 336, 892

Kirk D., Lahav O., Bridle S., Jouvel S., Abdalla F. B., Frieman J. A., 2015, MNRAS, 451, 4424

Kitaura F.-S. et al., 2016, MNRAS, 456, 4156

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Leistedt B., Peiris H. V., 2014, MNRAS, 444, 2

Leistedt B., Peiris H. V., Mortlock D. J., Benoit-Lévy A., Pontzen A., 2013, MNRAS, 435, 1857

Lesgourgues J., Mangano G., Miele G., Pastor S., 2013, Neutrino Cosmology. Cambridge University Press, Cambridge

Lesgourgues J., Pastor S., 2006, Phys. Rep., 429, 307

Lesgourgues J., Pastor S., 2014, New J. Phys., 16, 065002

Levi M. et al., 2013, preprint (arXiv:1308.0847)

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Loureiro A. et al., 2018, preprint (arXiv:1811.02578)

LSST Dark Energy Science Collaboration, 2012, preprint (arXiv:1211.0310)

MacCrann N., Zuntz J., Bridle S., Jain B., Becker M. R., 2015, MNRAS, 451, 2877

Manera M. et al., 2013, MNRAS, 428, 1036

McDonald A. B., 2016, Rev. Mod. Phys., 88, 030502

McLeod M., Balan S. T., Abdalla F. B., 2017, MNRAS, 466, 3558

McQuinn M., White M., 2013, MNRAS, 433, 2857

Nicola A., Refregier A., Amara A., 2016a, Phys. Rev. D, 94, 083517

Nicola A., Refregier A., Amara A., 2016b, Phys. Rev. D, 94, 083517

Nicola A., Refregier A., Amara A., 2017, Phys. Rev. D, 95, 083523

Padmanabhan N. et al., 2007, MNRAS, 378, 852

Peebles P. J. E., 1973, ApJ, 185, 413

Percival W. J., Samushia L., Ross A. J., Shapiro C., Raccanelli A., 2011, Phil. Trans. R. Soc. Lond. Ser. A, 369, 5058

Percival W. J. et al., 2001, MNRAS, 327, 1297

Perlmutter S. et al., 1999, ApJ, 517, 565

Planck Collaboration et al., 2018, preprint (arXiv:1807.06209)

Planck Collaboration I 2016a, A&A, 594, A1

Planck Collaboration et al., 2016b, A&A, 594, A11

Planck Collaboration et al., 2016c, A&A, 594, A13

Rassat A., Land K., Lahav O., Abdalla F. B., 2007, MNRAS, 377, 1085

Raveri M., Hu W., 2018, preprint (arXiv:1806.04649)

Reid B. et al., 2016, MNRAS, 455, 1553

Riess A. G. et al., 1998, AJ, 116, 1009

Riess A. G. et al., 2016, ApJ, 826, 56

Riess A. G. et al., 2018, ApJ, 861, 126

Rollins R. P., 2015, PhD thesis, University College London, London, UK

Ross A. J. et al., 2011, MNRAS, 417, 1350

Ross A. J. et al., 2013, MNRAS, 428, 1116

Ross A. J. et al., 2017, MNRAS, 464, 1168

Salazar-Albornoz S. et al., 2017, MNRAS, 468, 2938

Satpathy S. et al., 2017, MNRAS, 469, 1369

Scharf C., Hoffman Y., Lahav O., Lynden-Bell D., 1992, MNRAS, 256, 229

Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Scolnic D. M. et al., 2018, ApJ, 859, 101

Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132 (SH16)

Swanson M. E. C., Tegmark M., Hamilton A. J. S., Hill J. C., 2008, MNRAS, 387, 1391

Szapudi I., Prunet S., Pogosyan D., Szalay A. S., Bond J. R., 2001, ApJ, 548, L115

Sánchez A. G. et al., 2017, MNRAS, 464, 1640

Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152

Thomas S. A., Abdalla F. B., Lahav O., 2010, Phys. Rev. Lett., 105, 031301

Thomas S. A., Abdalla F. B., Lahav O., 2011a, Phys. Rev. Lett., 106, 241301

Thomas S. A., Abdalla F. B., Lahav O., 2011b, MNRAS, 412, 1669

Wandelt B. D., Hivon E., Górski K. M., 2001, Phys. Rev. D, 64, 083003

Wang Y., 2008, Phys. Rev. D, 77, 123525

Wang Y. et al., 2017, MNRAS, 469, 3762

Wright E. L., Smoot G. F., Bennett C. L., Lubin P. M., 1994, ApJ, 436, 443

Xavier H. S., Abdalla F. B., Joachimi B., 2016, MNRAS, 459, 3693

Zhao G.-B. et al., 2017, MNRAS, 466, 762

## APPENDIX A: THE OVERDENSITY PSEUDO-$C_\ell$ ESTIMATOR

Our aim is to measure the angular power spectrum of the galaxy overdensity field, $\delta^g$. Let $\bar{\rho}^g$ be the average of $\rho^g$ over the sky and define the galaxy overdensity field to be

$$ \delta^g = \frac{\rho^g - \bar{\rho}^g}{\bar{\rho}^g} = \frac{\rho^g}{\bar{\rho}^g} - 1. \qquad (A1) $$

This field may be represented using spherical harmonic expansion:

$$ \delta^g(\theta, \phi) = \sum_{\ell=0}^{\ell_{\max}} \sum_{m=-\ell}^{\ell} d_{\ell m} Y_{\ell m}(\theta, \phi), \qquad (A2) $$

where the spherical harmonic coefficients $d_{\ell m}$ are defined by

$$ d_{\ell m} = \int \delta^g(\theta, \phi) Y_{\ell m}^*(\theta, \phi) d\Omega. \qquad (A3) $$

Here and in what follows we have fixed a coordinate system $(\theta, \phi)$ for the celestial sphere; the spherical harmonic functions are defined with respect to this coordinate system. Our estimator of the angular power spectrum of the data is then

$$ \hat{D}_{\ell} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} d_{\ell m} d_{\ell m}^* . \qquad (A4) $$

The averaging over $m$ is motivated by the assumed isotropy of the probability distribution governing the location of galaxies.

To handle the partial-sky case, let $\Omega_{tot}$ be the survey region and define

$$ J_{\ell m} = \int_{\Omega_{tot}} |Y_{\ell m}|^2 d\Omega. \qquad (A5) $$

This is a normalization factor due to the average of modes in the partial sky coverage; note that $J_{\ell m} = 1$ for a full-sky survey. There will also be a term correcting for bias introduced by the partial sky measurement. However, this term is proportional to the average field value; in our case this average vanishes, so the bias correction need not be made. See Appendix B for details.

We can repeat this analysis for galaxy density fields $\rho^{g,i}$ and $\rho^{g,j}$ defined in tomographic bins $i$ and $j$. Combining the partial sky effect and tomographic binning results in an estimator $\hat{D}_{\ell}^{ij}$ for the cross- $(i \neq j)$ or auto- $(i = j)$ power spectrum of the data

$$ \hat{D}_{\ell}^{ij} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} D_{\ell m}^{ij} \qquad (A6) $$

where

$$ D_{\ell m}^{ij} = \frac{\mathrm{Re}(d_{\ell m}^{i} d_{\ell m}^{*})}{J_{\ell m}}. \qquad (A7) $$

Here we take the real part $\mathrm{Re}(\cdot)$ of a quantity whose expectation value will have no imaginary part.

In reality we work with a pixelized celestial sphere and we measure not $\rho^g$ but rather a galaxy count $n_p^g$ per pixel $p$. From
this we derive the per-pixel galaxy overdensity

$$ \delta_p^g = \frac{n_p^g}{\Delta\Omega_p} \frac{\Delta\Omega_{\text{tot}}}{n_{\text{tot}}^g} - 1, \qquad (A8) $$

where $n_{\text{tot}}^g$ is the total galaxy count, $\Delta\Omega_p$ the solid angle subtended by pixel $p$, and $\Delta\Omega_{\text{tot}}$ the total solid angle of the survey region $\Omega_{\text{tot}}$.

On the pixelated sphere, the spherical harmonic coefficients are estimated by

$$ d_{\ell m} = \sum_{p} \delta_{p}^{g} Y_{\ell m}^{*}(\theta_{p}, \phi_{p}) \Delta\Omega_{p}, \qquad (A9) $$

where $(\theta_p, \phi_p)$ are the coordinates at the centre of pixel $p$, $\Delta\Omega_p$ is the area of $p$, and the sum is over all pixels in the survey region.

Pixelization is a smoothing operator, and hence suppresses power at small scales. We summarize here the standard treatment of this effect; see Górski et al. (2005), Leistedt et al. (2013), and the HEALPix documentation for details.<fnref n="14" /> Consider the contribution of a given pixel $p$ to $d_{\ell m}$ both for the (measured) pixelized field and for the (desired) ideal continuous field; the ratio of these quantities is

$$ w_{\ell m}^{p} = \frac{\int_{p} Y_{\ell m}(\theta, \phi) d\Omega}{Y_{\ell m}(\theta_{p}, \phi_{p}) \Delta\Omega_{p}}. \qquad (A10) $$

This quantity depends sensitively on $\ell$: for small $\ell$, $Y_{\ell m}$ is slowly varying and hence $w_{\ell m}^p$ will be close to unity while for large $\ell$ the rapidly varying $Y_{\ell m}$ will have vanishing integral over $p$. However the dependence on $m$ and $p$ will be small and can be averaged out (in quadrature), yielding:

$$ w_{\ell}^{2} = \frac{1}{N_{\text{pix}}(2\ell + 1)} \sum_{p,m} |w_{\ell m}^{p}|^{2}. \qquad (A11) $$

The ratio of the power spectra of the (measured) pixelized overdensity field to that of the (desired) continuous field will then be $w_{\ell}^2$. This study uses a HEALPix resolution of $N_{\text{side}} = 512$; this means that at $\ell_{\text{max}} = 510$ this ratio of powers ($C_{\ell}^{\text{pix}}/C_{\ell}^{\text{unpix}}$) due to the pixel window function is then 0.911.

## APPENDIX B: CORRESPONDENCE BETWEEN OVERDENSITY AND NUMBER COUNTS PSEUDO-C<sub>ℓ</sub> ESTIMATORS

In Section 3, we showed the Pseudo-C<sub>ℓ</sub> estimator for galaxy overdensity maps. To link this with what is most commonly done in the literature, one can show that this galaxy overdensity measure is closely related to the more familiar galaxy number counts estimator as seen in Peebles (1973), Scharf et al. (1992), Fisher et al. (1994), Blake et al. (2007), and Thomas et al. (2011b). For the purpose of this section, we define the galaxy overdensity quantities with an upper $\delta$ index and the number counts quantities with a $n$ upper index. Example: the galaxy overdensity angular power spectra is represented by $C_{\ell}^{\delta}$.

We start the derivation by multiplying the overdensity spherical harmonics coefficients from equation (A9) by $\tilde{n}_i^g = n_{\text{tot},i}^g / \Delta\Omega_{\text{tot}}$.

Equation (A7) then becomes:

$$ \tilde{n}_i^g \tilde{n}_j^g D_{\ell m}^{\delta, ij} \approx \left[ \sum_{p}^{N_{\text{pix}}} \delta_{p,i}^g \Delta\Omega_p \tilde{n}_i^g Y_{\ell m}^*(\theta_p, \phi_p) \right] \times \left[ \sum_{p}^{N_{\text{pix}}} \delta_{p,j}^g \Delta\Omega_p \tilde{n}_j^g Y_{\ell m}^*(\theta_p, \phi_p) \right] \qquad (B1) $$

where we bear in mind that different subsamples $i$ and $j$ can have different total numbers of galaxies and different galaxies in each pixel, but use the same pixels.

Using equation (A8), we can write

$$ \delta_p^g \Delta\Omega_p \tilde{n}_i^g = n_p^g - \Delta\Omega_p \tilde{n}_i^g. \qquad (B2) $$

Now, we can use the above expression to rewrite equation (B1) as

$$ \tilde{n}_i^g \tilde{n}_j^g D_{\ell m}^{\delta, ij} \approx \left[ \sum_{p}^{N_{\text{pix}}} (n_{p,i}^g - \tilde{n}_i^g \Delta\Omega_p) Y_{\ell m}^*(\theta_p, \phi_p) \right] \times \left[ \sum_{p}^{N_{\text{pix}}} (n_{p,j}^g - \tilde{n}_j^g \Delta\Omega_p) Y_{\ell m}^*(\theta_p, \phi_p) \right]. \qquad (B3) $$

Let us take one of the square brackets and analyse the sub-terms individually:

$$ \sum_{p}^{N_{\text{pix}}} Y_{\ell m}(\theta_p, \phi_p) \Delta\Omega_p \approx \int Y_{\ell m}^*(\theta, \phi) d\Omega \equiv I_{\ell m}. \qquad (B4) $$

We can therefore see that this term is approximately equivalent to the shot noise correction term $I_{\ell m}$ from Blake et al. (2007) and Thomas et al. (2011b). The second term can also be re-expressed as

$$ \sum_{p}^{N_{\text{pix}}} Y_{\ell m}(\theta_p, \phi_p) n_p^g \approx \sum_{g'} Y_{\ell m}(\theta_{g'}, \phi_{g'}) = \sum_{g'} \int \delta_D(x_{g'} - x) Y_{\ell m}(\theta, \phi) d\Omega \qquad (B5) $$

where the index $g'$ runs over galaxies in the sample that have not been excluded by the mask and $\delta_D(x)$ is the Dirac delta function. We can reverse the order of summation and integration, and express the number count function as

$$ \sigma_1 = \sum_{g'} \delta(x_{g'} - x), \qquad (B6) $$

i.e. the galaxy distribution is a sum of delta functions at the locations of the galaxies, and hence the integral over this function is the total number of galaxies in that area. The function $\sigma_1$ is the filtered galaxy distribution, which has been masked. It is related to the full galaxy distribution $\sigma_0$ by

$$ \sigma_1(\theta, \phi) = \sigma_0(\theta, \phi) W(\theta, \phi), \qquad (B7) $$

where $W: S^2 \to \mathbb{B}$ is a binary filter, and

$$ \sigma_0 = \sum_{g} \delta(x_g - x) \qquad (B8) $$

runs over the full underlying set of galaxies. We can therefore write

$$ \sum_{p}^{N_{\text{pix}}} Y_{\ell m}(\theta_p, \phi_p) n_p^g \approx \int \sigma_0(\theta, \phi) W(\theta, \phi) Y_{\ell m}(\theta, \phi) d\Omega \qquad (B9) $$

<footnote n="14">¹⁴For details on the pixel window function: https://healpix.jpl.nasa.gov/htm/l/intronode14.htm</footnote>
$$= a_{\ell m} \qquad (B10)$$

where $a_{\ell m}$ are the spherical harmonic coefficients of the filtered galaxy number count field. Finally, we end up with

$$\bar{n}_i^g \bar{n}_j^g D_{\ell m}^{\delta,ij} \approx \frac{[a_{\ell m}^i - \bar{n}_i^g I_{\ell m}][a_{\ell m}^j - \bar{n}_j^g I_{\ell m}]}{J_{\ell m}} \qquad (B11)$$

$$= D_{\ell m}^{n,ij}; \qquad (B12)$$

in other words, the overdensity and number count power spectra differ only by a factor of the number density of galaxies in each tomographic bin involved.

## APPENDIX C: CODE COMPARISON

We produce our results using CLASS (Blas et al. 2011) (background evolution and perturbations) and the Unified Cosmological Library for $C_\ell$ s code, UCLCL (projected statistics). Here, we show a comparison for $C_\ell$ s calculated with both CLASS (integrated functionality from the former CLASSGAL code (Di Dio et al. 2013)) and CAMBSources (Challinor & Lewis 2011), matching cosmologies as closely as possible. We also show the derivatives calculated with respect to key cosmological parameters. In this comparison we used Gaussian redshift bins, since this is the functionality provided in CLASS and CAMBSources. Two redshift bins are chosen with $\bar{z} = \{0.5, 0.6\}$ and $\sigma_z = 0.05$ to be of comparable size to the redshift bins used in the body of the paper; auto and cross-correlations are calculated. Codes are run with their default accuracy parameters.

### C1 Auto- and cross-correlation precision

The auto-power spectrum for a bin with $\bar{z} = 0.5$ and $\sigma_z = 0.05$ is shown in Fig. C1, calculated in each of the three codes for a flat $\Lambda$CDM cosmology with $\Omega_b = 0.05$, $\Omega_{cdm} = 0.25$, $h = 0.67$, $\log(A_s \times 10^{-10}) = 3.2$, $n_s = 0.95$. Codes are in sub-percent level agreement up to $\ell \approx 200$ (although the CLASS low $\ell$ RSDs disagree to a slightly larger extent), after which there is a small discrepancy between CLASS and CAMBSources non-linear

<figure><img src="image_24.png" /><figcaption>Figure C1. Auto-correlation $C_\ell$ ($\bar{z} = 0.5$, $\sigma_z = 0.05$) comparison the three codes UCLCL, CLASS, and CAMBSources. The upper panel shows the three $C_\ell$ s over-plotted, whilst the lower panel shows the percentage difference between UCLCL/CLASS compared to CAMBSources: $\frac{C_{\text{CAMB}} - C_{\text{UCLCL/CLASS}}}{C_{\text{CAMB}}} \times 100$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure C2. Cross-correlation $C_\ell$ ($\bar{z}^i = 0.5$, $\bar{z}^j = 0.6$, $\sigma_z = 0.05$) comparison of the three codes UCLCL, CLASS, and CAMBSources. The upper panel shows the three $C_\ell$ s over-plotted, whilst the lower panel shows the percentage difference between UCLCL/CLASS compared to CAMBSources: $[C_{\text{CAMB}} - C_{\text{UCLCL/CLASS}}]/C_{\text{CAMB}} \times 100$. Again UCLCL follows CLASS closely, except in the RSDs and in a distinctive wobble around $\ell \approx 50$ where CLASS is transitioning away from the Limber approximation scheme.</figcaption></figure>

density perturbations. As one might expect, the differences between UCLCL and CAMBSources trace the differences between CLASS and CAMBSources as the former two share the same perturbations, i.e. $P(k)$.

The same trend is observed in the cross-correlations in Fig C2, with a notable wobble in the CLASS cross-correlation presumably when transitioning between approximation schemes (and thus possibly remedied by adjusting accuracy parameters away from the default).

### C2 Sensitivity to cosmological parameters

In order to check that the accuracy of the codes is not strongly cosmology dependent, the comparison are also made for variations on $h$ and $w_0$ over sensible ranges of the parameters. It is crucial that the sensitivity to the cosmological parameters not to be overwhelmed by the (approximately per cent level) uncertainty in the $C_\ell$ calculation itself. It is also important to check that the derivatives w.r.t. the cosmological parameters are consistent between the codes, as this will ensure $C_\ell$ s change consistently as one moves away from the fiducial cosmology.

In Fig. C3 one can see that $C_\ell$ s for $h = 0.64, 0.67, 0.70$ are clearly delineated and their differences significantly larger than the differences between $C_\ell$ s from different codes. With $w$ over the range $-1.1$ to $-0.9$, shown in Fig. C5, one can see that at low $\ell$ the $C_\ell$ s are well distinguished from each other, but at high $\ell$ $w_0$ has little effect, and thus is unlikely to be distinguished from the uncertainties inherent in the non-linear regime. In Fig. C6 one can also see that the variation at high $\ell$ is significantly different for UCLCL and CAMBSources, likely originating from the difference in the perturbations between CLASS and CAMBSources. Nevertheless, the shape of the derivatives w.r.t. to $w_0$ up to $\ell \approx 200$, and w.r.t. $h$ throughout the $\ell$ range, look consistent with CAMBSources. This shows that the $C_\ell$ s are changing in the correct way around this fiducial cosmology, and will yield the correct shape of posterior contours.
<figure><img src="image_26.png" /><figcaption>**Figure C3.** The top panel shows the auto-correlations ($\bar{z} = 0.5$) for three values of $h$ calculated in UCLCL. The lower panel shows the percentage difference of each of these $C_\ell$s with the corresponding $C_\ell$s from CAMB-Sources (matching values of $h$).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure C4.** Comparison of $C_\ell$ derivatives $\frac{dC_\ell}{dh}$ between UCLCL and CAMBSources.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure C5.** The top panel shows the auto-correlations ($\bar{w} = 0.5$) for three values of $w$ calculated in UCLCL. The lower panel shows the percentage difference of each of these $C_\ell$s with the corresponding $C_\ell$s from CAMB-Sources (matching values of $w$).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure C6.** Comparison of $C_\ell$ derivatives $\frac{dC_\ell}{dw}$ between UCLCL and CAMBSources. Here we see a more significant difference at high $\ell$, which can also been seen in Fig. C5. This characteristic bump appears to come from a difference in the CLASS and CAMBSources non-linear perturbations.</figcaption></figure>

This paper has been typeset from a T\textsc{eX}/I\textsc{a}\textsc{t}\textsc{eX} file prepared by the author.