ZXCorr: Cosmological Measurements from Angular Power Spectra Analysis of BOSS DR12 Tomography

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 9,376 deg$^2$ with an effective volume of $\sim 6.5$ (Gpc $h^{-1}$)$^3$ in the redshift range $0.15 \leq z<0.80$. We split this sample into 13 tomographic bins ($\Delta z = 0.05$); angular power spectra were calculated using a Pseudo-$C_{\ell}$ estimator, and covariance matrices were estimated using log-normal simulated maps. Cosmological constraints obtained from these data were combined with constraints from Planck CMB experiment as well as the JLA supernovae compilation. Considering a $w$CDM cosmological model measured on scales up to $k_{max} = 0.07h$ Mpc$^{-1}$, we constrain a constant dark energy equation-of-state with a $\sim 4\%$ error at the 1-$\sigma$ level: $w_0 = -0.993^{+0.046}_{-0.043}$, together with $\Omega_m = 0.330\pm 0.012$, $\Omega_b = 0.0505 \pm 0.002$, $S_8 \equiv \sigma_8 \sqrt{\Omega_m/0.3} = 0.863 \pm 0.016$, and $h = 0.661 \pm 0.012$. For the same combination of datasets, but now considering a $\Lambda$CDM model with massive neutrinos and the same scale cut, we find: $\Omega_m = 0.328 \pm 0.009$, $\Omega_b = 0.05017^{+0.0009}_{-0.0008}$, $S_8 = 0.862 \pm 0.017$, and $h = 0.663^{+0.006}_{-0.007}$ and a 95\% credible interval (CI) upper limit of $\sum m_{\nu}<0.14$ eV for a normal hierarchy. These results are competitive if not better than standard analyses with the same dataset, and demonstrate this should be a method of choice for future surveys, opening the door for their full exploitation in cross-correlations probes.


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-percent 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 (Hildebrandt et al. 2017;DES Collaboration et al. 2017;Planck Collaboration et al. 2018). Despite some suggested external and internal tensions in cosmological data sets (Riess et al. 2018;Planck Collaboration et al. 2018;Efstathiou & Lemos 2018), this so-called flat ΛCDM model is supported by observations from galaxy photometric surveys such as KiDS (Hildebrandt et al. 2017) and DES (DES Collaboration et al. 2017), cos-mic microwave background (CMB) fluctuations (Planck Collaboration et al. 2018), type Ia supernovae , 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 . Upper bounds on the sum of their masses are currently reaching a level where strong constraints on their hierarchy are possible (Planck Collaboration et al. 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 crosscorrelations 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 (Giannantonio et al. 2016;Hand et al. 2015), '3x2pt' correlations between galaxy clustering and lensing signals which provide the strongest low-redshift constraints on cosmological models (Hildebrandt et al. 2017;DES Collaboration et al. 2017), and also between galaxy clustering and CMB (Nicola et al. 2016a(Nicola et al. , 2017Doux 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;Górski et al. 2005;Wandelt et al. 2001;Szapudi et al. 2001;Hivon et al. 2002), 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 log-normal 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 read-ily implemented, and more robust to potential systematic errors (McQuinn & White 2013;McLeod et al. 2017) when compared to other methods such as P(k), ξ(r) and w(θ) (Ross et al. 2017;Salazar-Albornoz et al. 2017). We argue this should be the case because methods which live in angular space such as our method and w(θ) 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 method that 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, ξ(r) (Percival et al. 2001;Ross et al. 2017;Beutler et al. 2017a;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. Firstly, 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 signalto-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. 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 Baryon Oscillation Spectroscopic Survey (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 baryonphoton 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 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 ΛCDM, wCDM and ΛCDM with m ν 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 ∆z = 0.05, and using both the auto power spectra and the cross power spectra of adjacent bins to extract information from the radial correlation of galaxies. We refer to this method as 'ZXCORR'. Further combination with external CMB (Planck Collaboration et al. 2016c, as the Planck 2018 likelihood was still 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 organised 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 estimator used for the angular power spectrum analysis. Section 4 describes the theoretical modelling of the angular power spectrum and the use of log-normal 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 + m ν models using the BOSS C s alone and in combination with external CMB results from the Planck collaboration (Planck Collaboration et al. 2016b) and type Ia supernovae results from JLA (Betoule et al. 2014).

BOSS DR1DATA
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 colourcolour 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 2 http://www.sdss3.org/science/boss publications.php 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).

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. 3 We have further restricted these samples by applying redshift cuts of 0.15 ≤ z < 0.45 for LOWZ and 0.45 ≤ 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% more galaxies in 0.4 < z < 0.45 than does CMASS. See figure 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.

LOWZ sample
The LOWZ sample contains luminous red galaxies (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: 16 < r cmod < 19.6 (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 to 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 Figure 2(a)).

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 ⊥ (Equation 6) results in an increase in the number density of objects for the redshift range of 0.45 < z < 0.80 (see Figure 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 ps f − i mod and z ps f − z 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: 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.

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): • w systot , a combination of stellar density with airmass, sky flux, reddening, and other seeing conditions; • w cp , which is due to close-pair objects, i. e., objects that can not have both their spectra measured due to fibre collisions; • w noz , which takes into account nearest neighbours following a redshift failure by up-weighting such galaxies.
We follow Equation 50 in Reid et al. (2016) to combine these into a single weight for each galaxy (Ross et al. 2017): The default values of w cp and w 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.

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 4 software for pixelising the celestial sphere (Górski et al. 2005). The procedures described here were used for both CMASS and LOWZ. 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 s k y = 0.2067). (b) CMASS final pixel completeness angular mask with N si d e = 512. After performing a pixel completeness cut of 0.8, the total used area of the mask is around 9444.63 deg 2 ( f s k y = 0.2286).

Masks and angular selection function
The BOSS collaboration provides 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: where: • N obs is the number of spectroscopically observed objects including galaxies, stars, and unclassified objects; • N cp is the number of close-pair objects; • N 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 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 pixelisation with N side = 16384. Using this pixelisation 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 5 http://data.sdss3.org/sas/dr12/boss/lss/ binary mask. This is our final mask product and can be seen in Figures 2(a) for LOWZ and 2(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 ≥ 0.8: values < 0.8 are set to 0 and values ≥ 0.8 are set to 1.

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 ∆z = 0.05. This gives six tomographic bins for LOWZ (0.15 ≤ z < 0.45) and seven for CMASS (0.45 ≤ z < 0.80). Details about each redshift bin can be found in table 1. Next, we create a weighted number counts map which 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 upweight the maps by the inverse of the pixel completeness factor from the masks, C pix . Here, objects in 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: wheren 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; then i are not weighted by the pixel completeness weight. After these procedures are applied to all 13 tomographic redshift bins, we are ready measure the power spectra of these maps using the Pseudo-C estimator described in the next sections.

ANGULAR POWER SPECTRA ESTIMATORS AND MEASUREMENTS
The first proposed method for estimating the angular power spectrum C (Peebles 1973) consists of projecting the density field onto 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).

Pseudo-C Estimator
We describe the PCL estimator, following recent approaches as presented in e.g. Scharf et al. (1992), Fisher et al. (1994), Wright et al. (1994), Huterer et al. (2001), Hivon et al. (2002), Blake et al. (2004), Blake et al. (2007), Thomas et al. (2010), Thomas et al. (2011a), Thomas et al. (2011b), and Ho et al. (2012). Our aim is to measure the angular power spectrum of the galaxy overdensity field, δ g . Letρ g be the average of ρ g over the sky and define the galaxy overdensity field to be This field may be represented using spherical harmonic expansion: where the spherical harmonic coefficients d m are defined by Here and in what follows we have fixed a coordinate system (θ, φ) 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 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 Ω tot be the survey region and define This is a normalization factor due to the average of modes in the partial sky coverage; note that J m = 1 for a fullsky 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 A for details. We can repeat this analysis for galaxy density fields ρ g,i and ρ g, j defined in tomographic bins i and j. Combining the partial sky effect and tomographic binning results in an estimatorD i j for the cross-(i j) or auto-(i = j) power spectrum of the datâ Here we take the real part Re() of a quantity whose expectation value will have no imaginary part. In reality we work with a pixelised celestial sphere and we measure not ρ g but rather a galaxy count n g p per pixel p. From this we derive the per-pixel galaxy overdensity where n g tot is the total galaxy count, ∆Ω p the solid angle subtended by pixel p, and ∆Ω tot the total solid angle of the survey region Ω tot .
On the pixelated sphere, the spherical harmonic coefficients are estimated by where (θ p , φ p ) are the coordinates at the centre of pixel p, ∆Ω p is the area of p, and the sum is over all pixels in the survey region. Pixelisation is a smoothing operator, and hence suppresses power at small scales. We summarise here the standard treatment of this effect; see Górski et al. (2005); Leistedt et al. (2013) and the HEALPix documentation for details. 6 Consider the contribution of a given pixel p to d m both for the (measured) pixelised field and for the (desired) ideal continuous field; the ratio of these quantities is This quantity depends sensitively on : for small , Y m is  ). The solid grey line shows the deconvolved spline used in Section 4.3 to generate the log-normal mocks for covariance matrix estimation from which the error bars in this figure were estimated. Even though the measuredŜ ∆ 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.
slowly varying and hence w p m will be close to unity while for large the rapidly varying Y 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: The ratio of the power spectra of the (measured) pixelised overdensity field to that of the (desired) continuous field will then be w 2 . This study uses a HEALPix resolution of N side = 512; which means that at max = 510 this ratio of powers (C pix /C unpix ) is then 0.911 due to the pixel window function.
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: For most tomographic bins we can (using the variance of a Poisson distribution) approximate the angular power spectrum of the noise as: wheren is the mean number of galaxies per steradians. Amending (21) to account for pixelisation and shot noise yields an estimatorŜ i j for the partial sky signal power spectrum between two redshift bins i and j : The estimator is symmetric in i and j. Note also that there is no shot noise contribution for the cross-power spectra (i 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 A 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.

Bandwidth Binning and Measurements
We bin the values into bins ∆ of width 8 (so e.g. the first bin is 2 ≤ ≤ 9). For each bin we calculate a weighted average of theŜ i j (weighted by the number of spherical harmonic This binning acts on the measurement in a way that decorrelates mixed modes (that arise from the convolution of the true measurement and the survey's angular window function). We measure the PCL estimator up to max = 510; Figure 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 log-normal 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 characterisation of this noise will be included into the theoretical forward modelling presented in Section 4 and marginalised over during the cosmological parameter estimation (Section 6).

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 δ g (x, z) denote the galaxy density function. Let δ g (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 δ(k, 0) of the underlying matter distribution at the current time: The correlation structure of the Fourier transform is where P g (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,n, yields: where n(z) is the normalised redshift-dependent selection function and χ(z) is the comoving distance. The spherical harmonic components a m of this projected galaxy distribution are: where j (k χ(z)) are the Spherical Bessel functions (Thomas et al. 2010(Thomas et al. , 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: Using the window function in Equation (36) yields a simple expression for the angular power spectrum: Here we have introduced superscripts i and j to denote different redshift shells and the equation above therefore defines both auto-C (for i = j) and cross-C (for i j). The same formalism can be used to obtain the C between two different tracers, between photometric and spectroscopic redshift shells, etc. This work uses the UCLCL code (Cuceu et al., in prep.). This code obtains the primordial power spectra and transfer functions from the CLASS Boltzmann code (Blas et al. 2011), and then applies Equation (39) to obtain the angular power spectrum. UCLCL deals with the redshift distribution in more flexible ways than does CLASS and CAMB (Lewis et al. 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 40 in Section 4.1.1). A comparison between these codes is presented in Appendix.B.

Spectroscopic redshift distribution and shot-noise modelling
The spectroscopic selection provides a full (unnormalised) n(z) function for both the LOWZ and CMASS samples (see Figure 1). Binning is achieved by hard cuts on each of these samples in intervals of ∆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 where n i * (z) is the raw redshift distribution, σ 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 (FFT) 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 σ s (z). This effect, which is similar to that of redshift-space distortion (RSD), arises from the peculiar motions of galaxies within virialised 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 σ s models the combined impact of spectroscopic redshift errors and of the FoG effect; σ s is then varied and marginalised over during the cosmological analysis. Due to the sensitivity of the crossangular power spectra to these effects, a separate σ s is used for each redshift bin (for more details see Section 6.1).

Redshift space distortions
The full large scale structure window function needs to take into account Redshift Space Distortions (RSD) (Blake et al. 2007;Padmanabhan et al. 2007;Thomas et al. 2011b). This effect tends to increase the power for large scales, ∆ < 60, due to the mix of redshift and peculiar velocities of galaxies. This local peculiar motion of galaxies creates the illusion that the ones moving towards us appear closer (i. e., they 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 lower 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 (39): Here the RSD window function is given by: where we defined the redshift distortion parameter, , 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 4.1.1, the FoG effect is subsumed into the spread of the spectroscopic redshift distribution. Figure 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. . Slices through the mixing matrices (figs. 5(a) and 5(b)) for LOWZ (resp. CMASS) using = 200(250) and = 300(350). Amplitudes were normalized by R . As expected, the maximum amplitude peaks in the fixed and approaches zero in a given ∆ . The shape of both matrices remains the same as a function of . This shape indicates the correlation introduced due to -mixing by the convolution between the mask and the true signal, present in the PCL measurements ( Figure 3).

Non-linear angular power spectra: Halofit
In the UCLCL pipeline, the C estimation may be extended some way into the non-linear regime by introducing the scale-dependent non-linear overdensity δ N L (k, χ), and therefore the corresponding non-linear growth function 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 of the non-linear perturbations to the linear (δ L (k, χ)), the second equality follows from P = δδ * . This ratio is calculated using the modified halofit of Taka The window function in Equation (41) contains both the selection function and the growth function, which tracks the ratio of the power spectrum at different redshifts. The nonlinear power spectrum is related to the linear, present day power spectrum by: This means that the window functions in Equation (41) should have an additional factor of R N L (k, χ) inside the integral over χ. In the case of these very narrow spectroscopic redshift bins, we take wherez 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: The rest of the calculation may proceed as usual.

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 (Hauser & Peebles 1973;Hivon et al. 2002;Brown et al. 2005;Blake et al. 2007): The mixing matrix (see Figures 5 and 6) itself depends only on the survey's geometry through the mask's angular power spectrum where (see Appendix A) the mixing matrix is then The 2 × 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 Figure 5 and in more detailed slices in Figure 6 which gives an intuition about the size of ∆ -bands used to bin the measuredŜ s as it shows the range of multipoles that are mixed due to the survey's mask. This small correlations between the multipoles can be "washed away" by binning our measurements. In addition, as can be seen in Figure 4, the mixing matrix convolution tends to suppress power in all scales. Finally, after being convolved with the mixing matrix (Equation 48), the theoretical S is binned in the same way as the data in Equation (30):

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:  where P is the projector in pixel space, defined as: The projector satisfies the following identity in the full sky case: Using this identity, the Pseudo-C estimator from Equation (29) can be written in terms of the projector and data covariance from (53): where ∆Ω p is the area the pixels assumed to have an equal area. Assuming a Gaussian signal (Blake et al. 2007), the covariance matrix for the angular power spectra estimator between different multipoles and can be expressed as: The symmetry of the P and D matrices, together with the definition Cov(X, X ) = XX − X X , allows us to rewrite the covariance as: This expression works for both full and partial sky cases. The difference between the two cases appears on the projector identity from Equation (55). Using the definition of the pixel space projector (Equation 54), the I m expression from Equation (50), 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 (57) for a partial sky observation with area ∆Ω tot as: where the last equality uses the definition of the mixing matrix from Equation (51) and f sky is the observed fraction of the sky. This expression is similar to ones used in Blake et al. (2007); Padmanabhan et al. (2007); 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 generalise the data angular power spectra by changing S + N to and we include the effect of cross-correlation in the covariance by changing in Equation (60) (Rassat et al. 2007). The final expression for the angular power spectra theoretical covariance matrix is therefore: By performing these modifications, Equation (61) 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. Figure 7 shows a comparison between the variance (the diagonal) of equation (61) with the variance from the estimated covariance matrix from Section 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 realisations (Blake et al. 2007;Thomas et al. 2011b;Nicola et al. 2016b) or the mocks provided by the BOSS Collaboration (Kitaura et al. 2016;Manera et al. 2013). However this work instead uses log-normal 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 7 7 http://www.astro.iag.usp.br/~flask/ (Xavier et al. 2016), a publicly available code that produces log-normal simulations of correlated fields on the sphere. We used the dataŜ 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 s measured from the data. For each sample we produced 6,000 log-normal 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,S( ), using theŜ ∆ measurements ( Figure 3) and a Gaussian filter to smooth the measurements.
2. Deconvolve the mixing matrix R from the splines to 3. Monotonically extrapolate the splines to 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 = max = 2048. 8 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 2 (as described in Equations (29) and (30)) 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: Here N S is the number of simulations. To validate the estimated covariance matrix, we compared the diagonal of the covariance matrix in Equation (63) with the expression for the theoretical variance for the measured angular power spectra in Equation (61); Figure 7 shows a typical result.

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 datasets (Blake et al. 2007;Ross et al. 2011;Thomas et al. 2011b;Leistedt et al. 2013;Ho et al. 2012;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 deviates 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. Plots showing the complete analysis of systematics contamination for all 13 tomographic bins can be found in Appendix C.

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 pixelisation of the observed sky in MANGLE format and another that records the observational information for each MANGLE polygon. 9 The first step is to reconstruct the MANGLE maps for each observational systematic from these files. We use the MANGLE python wrapper 10 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 9 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. 10 https://github.com/mollyswanson/manglepy SDSS SkyServer Schema Browser. 11 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: where the extinction-corrected magnitude cut ensures robust star selection (Padmanabhan et al. 2007), the type selection is the standard SDSS star-galaxy classifier, 12 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. 13 We use the original SFD scaling (Schlegel et al. 1998).
The MANGLE masks created from the SDSS FITS files are not appropriately snapped, pixelised and balkanised, 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 pixelisation anew, which corrects whatever imperfections remained in the first pass. From these masks, we create full-sky HEALPix maps at resolution N side = 16384, which defines an angular scale much smaller than the average resolution of the mask features. For each observational systematic, we populate the subresolution 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 side = 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 pix = 0.8 in the maps. This up-weighting is due to the pixelisation of the degraded mask, as explained in Section 2.2.1. From these post-processed maps, we create the systematics overdensity maps as:   Figure 8. 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 δ S y s 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 Figure  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 similar to the results shown in this figure and can be found in Appendix C. Note also that the first −band in the stellar overdensity cross-C is completely out of the acceptable range, which lead us to exclude this data point on all bins for both samples.
where n Sys p is the pixel value for a given systematic andn 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 are shown in Figure C1 using the CMASS mask as an example.

Cross-power spectra between data and systematic maps
For the systematics analysis using cross-power spectra, we will follow a data analysis in a similar fashion as the one performed for the galaxy overdensity maps in Section 3. Using Equation (24) we decompose the systematics overdensity maps, δ Sys , 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 (29): 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Ŝ gs (Equation 66) 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 crosspower spectra. Figure 8 shows an example for LOWZ-3 and CMASS-9 bins and cross-power spectra for all systematics. The full results are shown in Appendix C. From all of these, the majority of measurements are consistent with the variance of the data measured from the log-normal simulation (Section 4.3), which lead us to be confident in using the full shape of our measured C s. Note, however, that a few of the large scale measurements (low-) 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 -band, we decided to exclude this first point from our cosmological analysis (see Section 6 for details on the range used). As for the second -band ( = 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.

COSMOLOGICAL ANALYSIS
In this section, we present the cosmological implications from the measured angular power spectra of BOSS galaxies for flat ΛCDM, wCDM, and a ΛCDM with m ν 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 prep) nested sampler. Table 2. Maximum multipole considered in the cosmological analysis for each tomographic redshift bin. All the samples start at mi n = 13.5 and have a bandwidth of ∆ = 8. When considering the cross-power spectra between bins, the lower ma x is used. The 5% ma x column corresponds to a k ma x 0.07h Mpc −1 , and the 10% ma x column corresponds to a k ma x 0.

Likelihoods, Priors & 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) The posterior distribution of the cosmological parameters, Θ Θ Θ, given the measured angular power spectra,Ŝ ∆ , and a model M can be written as a marginalisation of the full posterior over the nuisance parameters, ν ν ν: The full posterior distribution can be written as: where L(Ŝ ∆ |Θ Θ Θ, ν ν ν, M) is the likelihood, π(Θ Θ Θ, ν ν ν) is the prior on the sampled parameters, and Z(Ŝ ∆ |M) is the evidence, which is calculated using PLINY , a nested sampler used in our analysis (Rollins 2015 and Rollins et al., in prep;Feroz & Hobson (2008)). If we had access to the true covariance matrix Σ, the likelihood, assumed here to be Gaussian, would be: where S th ∆ (Θ Θ Θ, ν ν ν) is the theoretical angular power spectra after being convolved with the mixing matrix (Equation 48) and binned into the same bandwidths as the data (Equation 52).
However, this is not the case when estimating the covariance matrix C from simulations (Equation 63). Even though C can be an unbiased estimator of the true covariance Σ, its inverse C −1 is not necessarily an unbiased estimator of the inverse covariance Σ −1 , needed to estimate the likelihood in Equation (69). Hartlap et al. (2007) proposed to keep using the Gaussian likelihood, and to apply a simple rescaling to BOSS C(l) SH16 BOSS C(l) Hartlap BOSS C(l) Gaussian Figure 9. Comparison between the three likelihood methods mentioned in Section 6.1 using the BOSS C data only for Ω 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 log-normal 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.
the inverse of the estimated covariance matrix in order to de-bias it (Anderson 2003).
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 marginalise 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): and Γ is the gamma function. Even though the non-linear model described in Section  Figure 10. 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 (61) (blue contours). Note that the same parameters were sampled in both cases as the theoretical covariance also depends on the same nuisance parameters. These marginalised credible intervals (CI) 1and 2-σ plots indicate both the estimated covariance matrix and the UCLCL pipeline robustness.  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 s and performed a preliminary cut in max where the percent deviation between the linear and non-linear models were smaller than 5%. We performed robustness checks on the 5% deviation cut choice by extending the cuts to max which had a deviation up to 20% and concluded that our cosmological results could be trusted up to a 15% deviation between linear and non-linear theories for this fiducial test. In this paper we present results where this percentage cut is 5% and 10%. For avoidance of doubt, for the fiducial cosmology of choice, applying a 5% implies rejecting the majority of modes k 0.07h Mpc −1 , whereas 10% implies rejecting the majority of modes k 0.1h 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.
We used a total of 28 nuisance parameters (ν ν ν) in the BOSS C likelihood analysis in most of our results for a 5% cut on max . These parameters are: a scale independent bias, b(z), for each redshift bin; a redshift error dispersion, σ s (z) (Equation 40), 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:Ŝ where N is a constant that takes into account extra shot noise due to the very low number of galaxy in these two redshift bins. In the case of a 10% we used two further shot noise nuisance parameters for bins 9 and 10 respectively as we go further into the non linear regime where the shot noise in those bins dominates over the signal. The prior ranges can be found in table 3 for all parameters considered in the cosmological analysis in this section: the baryonic matter density (Ω b ), the cold dark matter density (Ω cdm ), the amplitude of the primordial power spectrum (A s ), the spectral index (n s ), the Hubble constant (h), the equation-of-state of dark energy (w 0 ), the sum of neutrino mass species ( m ν ), the optical depth at reionisation epoch (τ reio ), the bias of BOSS galaxies as a function of redshift (b(z)), the redshift dispersion parameter for BOSS galaxies (σ s (z)), the extra shot-noise for BOSS galaxies (N i ), the Planck absolute calibration parameter (y Planck cal ), and the Absolute magnitude of SNe Ia at peak light in blue band (M J L A B ). Finally, to perform consistency checks between BOSS DR12 and the external datasets described in Section 6.3, we used the Bayes factor. The Bayes factor for the consistency of two datasets (A and B) is given by: or, for three datasets (A, B, and C): where M is the model, P( ì A, ì B|M) is the evidence when the model is fitted to both datasets simultaneously and P( ì A|M)P( ì B|M) is the product of the evidences when the model is fitted to each dataset individually. Since PLINY is a nested sampler, all our cosmological estimations lead to values for the evidences of each model, dataset, and combination of datasets. We are aware that this method has received some criticism as the value of R A, B depends crucially on the prior volume chosen for the sampling. Although this is true, choosing priors which are physically motivated should lead to values of R A,B which can be reflective of the real consistency between datasets A and B. However, a host of new methods have been proposed in the literature which we have not fully investigated in this paper, but which we intend to implement in further analysis with this dataset. 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 datasets. We refer the reader to these methods in Charnock et al. (2017a); Raveri & Hu (2018); 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. Figure 9 shows a comparison between the three likelihood for a wCDM model, using the BOSS C s only, for w 0 and Ω m . In all of the following results in this Section, we use the SH16 likelihood.

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.

Parameter-dependant theoretical covariance matrix
We implemented an alternative likelihood based on a theoretical expression for the covariance matrix (Section 4.2, see Equation (61)) where the signal angular power spectra, S , also depends on the sampled parameters. Most standard cosmological analysis in the literature (Alam et al. 2016;DES Collaboration et al. 2017;Hildebrandt et al. 2017) assume a covariance matrix which is independent of cosmology and which 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. Figure 10 shows the results for this test for a ΛCDM cosmological model.

Controlled cosmology pipeline test
For this test, we generated theory autos-and cross-C s to mimic the BOSS dataset using a Planck-like cosmology: h = 0.6725, Ω b = 0.0492, Ω 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 Figure 1. We chose the nuisance parameters to match the best fit values obtained in Section 6.4 from the combination of the entire cosmological dataset available to us. Using BOSS masks as input, we created FLASK mocks like as 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 Figure 11, where the recovered parameters are within the errors with no indications of biases in the entirety of the pipeline.

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 Figure 12 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).

Distribution of residuals
For a dataset with uncorrelated errors (diagonal covariance matrix), the vector of normalised residuals is given by: 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 et al. 2010). On the other hand, if the errors are estimated from the data, the residuals are given by a Student's t-distribution. This distribution approaches a Gaussian with 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, we either found the true model, or the current data is 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 (77) 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 : 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Λ −1 Q −1 , which transforms the χ 2 into: If we treat Λ as the new (diagonal) covariance matrix, it follows that the normalised residuals are now given by: where Ξ now contains the square roots of the eigenvalues. We use Equation (80) to calculate the residuals at our best-fit point in a flat ΛCDM Cosmology and plot the results in a histogram ( Figure 13). There are no significant deviations from a Gaussian with µ = 0 and σ = 1, which means the model seems to be a very good representation of the data.

External data
We compared and combined our results with results obtained from the Planck satellite CMB experiment (Planck Collaboration et al. 2016a), and Type Ia Supernovae from the Joint Light curve Analysis (JLA) collaboration (Betoule et al. 2014). The relevant likelihood codes for these probes were implemented and tested in the UCLCL 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 et al. 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-   The SN data from JLA was added to the UCLCL 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 4 nuisance parameters (α, β, M B , ∆ 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 the sampled cosmological parameters). We sample the absolute magnitude at peak brightness (M B ) as part of our analysis, and keep the other 3 nuisance parameters fixed to their best-fit values found by Betoule et al. (2014) since these have a small impact in the cosmological parameters when combined with the BOSS dataset.
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 measurements use 3 redshift bins with z e f f = 0.38, 0.51, and 0.61, where we used the full shape post-reconstruction measurements from the correlation function and 3D power spectra, which contains additional information from measurements of f σ 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 sessions to give the reader an impression of how our results compare with the BOSS alone results from Alam et al. (2016).

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. Ω k = 0, and varied five cosmological parameters: the baryonic density, Ω b ; the dark matter density, Ω 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 Λ, we fixed the w 0 parameter to a cosmological constant (w 0 = −1), therefore: Ω Λ = 1 − (Ω b + Ω cdm ). Here, we also fixed the sum of neutrino masses to the minimum found from neutrino oscillation experiments, m ν = 0.06 eV -the lower bound as suggested from neutrino oscillation experiments (Lesgourgues & Pastor 2006. From these, we also obtained derived parameters as the total matter density, Ω m ≡ Ω b + Ω cdm ; and the fluctuation of amplitude at 8 h −1 Mpc, σ 8 or S 8 = σ 8 Ω m /0.3. Finally, as described in the previous section, we also varied the BOSS, Planck and JLA  Figure 13. Comparison between the histogram of the distribution of residuals (given by Equation (80)) calculated at the best-fit 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.  2D Ω m × S 8 marginalised credible intervals for a ΛCDM Cosmology. In this figure we show in detail the cosmological results from Section 6.4 for BOSS C s only (blue); BOSS C s plus JLA (green); BOSS C s plus JLA and Planck (purple); together with results using the post-reconstruction full shape (incl. f σ 8 (z)) from Alam et al. (2016) consensus results (pink), and Planck alone (red). For details on the external datasets, see Section 6.3. nuisance parameters. For this analysis, we used the 5% max cuts (see table 2).
We checked the consistency of the datasets 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 (76), for combinations of the considered datasets indicates consistency between all three probes: these indicate that the datasets are compatible for the considered model, given the chosen priors. Finally, when considering the combination of BOSS C s, Planck and JLA, we find results consistent with Alam et al.  Figure 14 shows the Ω m × S 8 2D plane for this analysis and comparisons with Planck and the BOSS fullshape post-reconstruction from Alam et al. (2016). Despite the Bayes factors showing no significant reason to be concerned about the compatibility of these datasets we see an interesting trend in this figure insofar as a tension and the BOSS dataset in this paper preferring a smaller S 8 than the Planck analysis. We argue here that this method would prove potentially very useful in resolving any S 8 tensions which exist currently between CMB and weak lensing data (MacCrann et al. 2015;Charnock et al. 2017b).
The results for the 1D marginalised cosmological constrains for BOSS C and combinations, together with the 68% credible intervals can be found in table 4. The 1-σ and 2-σ contour levels can be found in Figure 15 where the nuisance parameters have been marginalised over. Figure 16 shows the best fit theory C using the parameters estimated from this analysis with a χ 2 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 ΛCDM with a fixed zero neutrino mass, m ν = 0 eV. We compared it with the model used in this section, ΛCDMwith m ν fixed to 0.06 eV, using the Bayes factor for model selection. Consider ì D representing the combination of data vectors; the Bayes factor is given by (85) This indicates that, for a ΛCDM model, the data prefers massive neutrinos over no neutrino mass at all.

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, Λ. The procedure for this analysis followed in similar fashion as the one outlined in Section 6.4, varying six parameters instead of five: Ω b , Ω cdm , n s , ln 10 10 A s , h, and the extra w 0 . Note that, for  Figure 15. Marginalised 1 & 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: Ω m , Ω b , S 8 , h, and n s (marginalising over τ r ei o for the Planck combinations). The blue contours where estimated from the BOSS C 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 , JLA and Planck (also high-TT, TE, EE and low-P). Note that none of the results here use Planck Lensing data.
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 m ν = 0.06 eV (Lesgourgues & Pastor 2006. Here, we used the same 5% max cuts as in the last section (see table 2). Figures 17 and 18 show in detail the contours for S 8 × Ω m and w 0 × Ω m , respectively, and comparisons with previous measurements in the literature. From the Figure 18 and from the complete set of results in 19 we show that a ∼ 4% error (1-σ CI) on the equation of state of dark energy is obtained from this cosmological analysis: This results is consistent with the ΛCDM scenario of standard cosmology, i. e., it is consistent with dark energy being a cosmological constant, Λ. Note from Figure 19 that we find a small value of h (compared to Planck Collaboration et al. 2016c) when combining BOSS C s, Planck, and JLA: This value is lower than the quoted Planck value alone; putting further tension in of this result if compared to the Hubble constant result from Cepheid Variables (Riess et al. 2016(Riess et al. , 2018). As the model in this section is different from the previous section, we performed an evidence analysis using the Bayes factor, Equation (76), in order to be sure that our measurements can be combined with the the external data described in Section 6.3. The following measurements indicate that such combinations are consistent for a wCDM model: Finally, we used the ratio of the evidences, the Bayes factor, to perform a model selection between wCDM and Shot-Noise BestFit Used data data from sample 1 data from sample 2 Figure 16. 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 estimates, described in Section 3, for LOWZ and CMASS. The solid blue lines, generated with UCLCL, reflect the best fit auto-and cross-power spectra for the Λ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 χ 2 for this fit is χ 2 r e d ≈ 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. 2D Ω m × S 8 marginalised 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).
ΛCDM using the final dataset combination. Assuming ì D to be the combination of data vectors for all the datasets, the Bayes factor between the two models is

Flat ΛCDM + m ν 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, m ν . In the previous sections, we fixed the sum of neutrino masses to m ν = 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 m ν for one remaining massive species. We do not investigate details of how the prior on the hierarchy or on the absolute mass change this result and we leave this to a future analysis. We fix N e f f = 3.046 by changing the values of massive neutrinos and ultra-relativist particles for the case considered, i. e., N ν = 1 and N ur = 2.0328.
We perform an analysis using the same -range as in the previous sections, 5% max from table 2. A summary for the marginalised 1D credible intervals from the cosmological 0 . estimation made with this cut can be found in the third part of table 4 showing the one sigma intervals for the ΛCDM parameters plus the 95% upper limit for m ν . The 1D and 2D marginalised credible intervals for this analysis can be found in Figure 20. When considering an approximation for the normal hierarchy, for a combination of BOSS C s, Planck CMB data, and supernovae data from JLA, the 95% upper limit for sum of neutrino masses is: m ν < 0.14 eV (BOSS + Planck + JLA -5% ma x cut) (93) From Figure 21 and even more so from Figure 22, one can notice that we are not so far from excluding zero total neutrino mass using cosmological data alone. As the power of such datasets 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 perform a consistency of datasets by using the evidence of each cosmological parameter estimation for these model to calculate the Bayes factor (Equation The blue contours where estimated from the BOSS C s data alone using the SH16 likelihood; green contours are a combination the BOSS likelihood and JLA data; red contours are the Planck high-TT, TE, EE and low-P likelihood results; finally, the purple contours are a combination of the three probes: BOSS C , JLA and Planck (also high-TT, TE, EE and low-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.
scales considered for the 10% max cuts (see table 2 for details). This allow us to access smaller scales that are still in the beginning of the so-called the weak nonlinear 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;Hildebrandt et al. 2017;DES Collaboration et al. 2017) -DES Collaboration et al. (2017, for example, uses scales up to 0.78h Mpc −1 . In other words, one can be confident that the 10% max 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 ΛCDM model with one massive species of neutrino, the approximation to the normal hierarchy. The 1D and 2D marginalised credible intervals for these final analysis can be find in Figure 23 and the marginalised 68% credible intervals for the ΛCDM parameters and the 95% credible interval upper limit for m ν using this cut can be found in table 4. The Bayes factors for this choice are shown below (note that we have 2 further nuisance parameters in the 10% max cut and we checked that failure to add these reduces the Bayes factor significantly). . 1D and 2D marginalised credible intervals for a ΛCDM Cosmology with m ν when using scales up to k ma x ≈ 0.07h Mpc −1 ( 5% ma x cut). Here we show the Ω m , Ω b , S 8 , h, n s , and m ν contours for BOSS C s alone (blue); BOSS C s plus JLA (green); Planck high-TT, TE, EE and low-P (red); and BOSS C s plus JLA and Planck high-TT, TE, EE and low-P (purple). As most of the scales that contain clean information on the neutrino masses are cut off, the 95% CI upper bound found is m ν < 0.14 eV.

CONCLUSIONS
In this work, we have taken a different approach 14 to obtain galaxy clustering information from the BOSS DR12 large scale structure catalogue (Reid et al. 2016). This approach consisted in using a pseudo angular power spectra estimator (PCL) applied to 13 tomographic redshift bins ranging from 0.15 ≤ 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 UCLCL (Cuceu et al., in prep). 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, in no moment in this analysis a fiducial cosmology was assumed. This is, by itself, a great advantage over methods that use P(k) or ξ(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 systematic 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 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 3 different models: ΛCDM, wCDM, and ΛCDM with m ν . 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 big collaboration results such as Dark Energy Survey and KiDS with errors as small as these collaborations (DES Collaboration et al. 2017;Hildebrandt et al. 2017).
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 a ∼4% error when combining BOSS C s, Planck CMB, and JLA Supernovae. This was never achieved before using C with a spectroscopic survey and has constraints 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 et al. 2017).
3. For the models and datasets considered, we find very high values for the Bayes factor, R, when combining BOSS C s, Planck, and JLA. We would like to highlight: R BOSS+PLANCK+JLA 4 × 10 4 for ΛCDM and R BOSS-10%+PLANCK+JLA 3 × 10 5 for ΛCDM varying neutrino masses.
4. The Bayes factor can also be used for model selection. Considering the combination of datasets, the Bayes factor between ΛCDM and wCDM is where ì D here just represents the overall combination of data vectors. This indicates that this combination prefers slightly more ΛCDM than wCDM, although no strong conclusion can be made.
5. We find a small tension between BOSS C s and Planck for S 8 in all models considered with BOSS preferring smaller values. For example, for ΛCDM: although the combination of these datasets prefers higher values such as Planck (see table 4) and the Bayes factor suggest the datasets are compatible. We conclude that such tension can be investigated further with this method as LSS data increases in size and depth. 6. Even though we do show these results, we performed a cosmological analysis using a ΛCDM model but fixing m ν = 0eV and compared with the ΛCDM results from Section 6.4, which has a m ν fixed to 0.06 eV. Using the Bayes factor for model selection, it is clear that the data prefers massive neutrinos against no neutrino mass at all: 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 m ν < 0.474 eV (95% CI). The reason we obtain much tighter constraints ( m ν < 0.14 eV (95% 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 "shellcrossing" (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. In Salazar-Albornoz et al. (2017), it is considered 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 et al. 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 m ν . In a future paper, we will study the impact of the model in the sum of neutrino masses and their hierarchy.
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. Finally, we would like to point out that the biggest 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 datasets from across the entire electromagnetic spectrum. The angular power spectrum approach offers a unified framework for coherently combining different datasets 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 . Cosmological constraints for the ΛCDM + m ν model, using the 10% ma x cut, varying now 6 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: Ω m , Ω b , S 8 , h, n s , and m ν . Note that the Planck chains also varied τ r ei o . The blue contours where estimated from the BOSS C s data alone using the SH16 likelihood; the green contours are a combination the BOSS likelihood and JLA data; the red contours are the Planck high-TT, TE, EE and low-P likelihood results; finally, the purple contours are a combination of the three probes: BOSS C , JLA and Planck (also high-TT, TE, EE and low-P). For this scale cut, the combination of datasets yields an upper bound on m ν < 0.16eV. 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.
Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No 6655919. We would also like to thank Niall Jeffrey and Dr. Andreu Font-Ribera for their great comments on the final draft of this work. All cosmological contour plots were generated using ChainConsumer 15 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 B: CODE COMPARISON
We produce our results using CLASS (Blas et al. 2011) (background evolution and perturbations) and the C estimation code UCLCL (projected statistics). Here, we show a comparison for C 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 withz = {0.5, 0.6} and σ 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.

B1 Auto-and cross-correlation precision
The auto-power spectrum for a bin withz = 0.5 and σ z = 0.05 is shown in Figure B1, calculated in each of the three codes for a flat ΛCDM cosmology with Ω b = 0.05, Ω cdm = 0.25, h = 0.67, log(A s × 10 −10 ) = 3.2, n s = 0.95. Codes are in sub-percent level agreement up to ≈ 200 (although the CLASS low RSDs disagree to a slightly larger extent), after which there is a small discrepancy between CLASS and CAMBSources non-linear density perturbations. As one might expect, the differences between uclcl and CAMBSources trace the differences between CLASS and CAMB-Sources as the former two share the same perturbations, i.e. P(k).
The same trend is observed in the cross-correlations in Fig B2, with a notable wobble in the CLASS crosscorrelation presumably when transitioning between approximation schemes (and thus possibly remedied by adjusting accuracy parameters away from the default).

B2 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 be overwhelmed by the (approximately percent level) uncertainty in the C 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 the C s change consistently as one moves away from the fiducial cosmology.
In Figure B3 one can see that the C s for h = 0.64, 0.67, 0.70 are clearly delineated and their differences significantly larger than the differences between the C s from different codes. With w over the range −1.1 to −0.9, shown in Figure B5, one can see that at low the C s are well distinguished from each other, but at high w 0 has little effect, and thus is unlikely to be distinguished from the uncertainties inherent in the non-linear regime. In Figure B6 one can  Figure B2. Cross-correlation C (z i = 0.5,z j = 0.6, σ z = 0.05) comparison the three codes UCLCL , CLASS , and CAMBSources. The upper panel shows the three C s over-plotted, whilst the lower panel shows the percentage difference between UCLCL / CLASS compared to CAMBSources: ×100. Again UCLCL follows CLASS closely, except in the RSDs and in a distinctive wobble around l ≈ 50 where CLASS is presumably transitioning between some approximation schemes. also see that the variation at high is significantly different for UCLCL and CAMBSources, likely originating from the difference in the perturbations between CLASS and CAMB-Sources. Nevertheless, the shape of the derivatives w.r.t. to w 0 up to ≈ 200, and w.r.t. h throughout the range, look consistent with CAMBSources. This shows that the C s are changing in the correct way around this fiducial cosmology, and will yield the correct shape of posterior contours.

APPENDIX C: FULL SYSTEMATICS ANALYSIS
We show in this appendix the full cross correlation analysis we have performed in order to confirm that there is no strong evidence for systematic effects which would bias the power spectrum analysis. We do that by cross-correlating selected redshift bins with the systematic maps we have produced. In Figure C1 we show all the systematic overdenstity maps, using the CMASS mask as an example, created as mentioned in Section 5.1. We have already shown these results for one of the bins in the body of the text and, for clarity, we show here, in Figures C2, C3, C4, C5, C6, and C7 the remaining cross-power spectra.  Figure B6. Comparison of C derivatives dC d w between UCLCL and CAMBSources. Here we see a more significant difference at high , which can also been seen in Figure B5. This characteristic bump appears to come from a difference in the CLASS and CAMBSources non-linear perturbations.
This paper has been typeset from a T E X/L A T E X file prepared by the author.  Figure C2. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ-0(CMASS-6) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the δ S y s 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.  Figure C3. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ-1(CMASS-7) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the δ S y s 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.  Figure C4. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ-2(CMASS-8) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the δ S y s 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.  Figure C5. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ-4(CMASS-10) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the δ S y s 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.  Figure C6. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1, and LOWZ-5(CMASS-11) tomographic bins in blue dots (red squares). The error-bars were obtained by cross-correlating the δ S y s 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.  Figure C7. Cross-power spectra between the 18 systematics overdensity maps produced in 5.1 CMASS-12 tomographic bin in red squares. The error-bars were obtained by cross-correlating the δ S y s 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.