Model-independent Constraints on Clustering and Growth of Cosmic Structures from BOSS DR12 Galaxies in Harmonic Space

Up to linear order in cosmological perturbation theory, we can write the observed galaxy number count fluctuations in configuration space and in the longitudinal gauge as (Yoo 2010; Bonvin & Durrer 2011; Challinor & Lewis 2011)

$$\begin{eqnarray}&&{\rm{\Delta }}=b\,\delta -\displaystyle \frac{{\partial }_{\parallel }^{2}V}{{ \mathcal H }}+{{\rm{\Delta }}}_{\mathrm{loc}}+{{\rm{\Delta }}}_{\mathrm{int}}.\end{eqnarray} \tag{ 1 }$$

The first term represents matter density fluctuations, with δ being the density contrast in the longitudinal gauge, modulated by the linear galaxy bias b. The second term is linear RSD, with ∂_∥ being the spatial derivative along the line-of-sight direction, $\hat{{\boldsymbol{r}}}$, V the peculiar velocity potential, and ${ \mathcal H }$ the conformal Hubble factor. The remaining two terms respectively collect all local and integrated contributions to the signal. Such terms are mostly subdominant and affect only the largest cosmic scales, ergo we neglect them hereafter. ⁴

Since measurements from a galaxy survey corresponds to angular positions of galaxies on the sky, we can naturally decompose the observed ${\rm{\Delta }}(\hat{{\boldsymbol{r}}})$ on the celestial sphere into its spherical harmonic coefficients, Δ_{ℓ m} . If additional information on the galaxies' redshift is available—as in the case of a spectroscopic galaxy survey—we can subdivide the data ${\rm{\Delta }}(\hat{{\boldsymbol{r}}},z)$, which is now also a function of redshift, into concentrical redshift shells. This effectively corresponds to dealing with a set of Δ_{i,ℓ m} = ∫dz n_i (z)Δ_{ℓ m} (z), where the index i runs over the number of redshift bins and n_i (z)dz is the number of sources per steradian in the ith redshift bin between z and z + dz. Note that n(z)dz = n(r)dr holds true, with r(z) being the radial comoving distance to redshift z.

Then, the harmonic-space tomographic power spectrum of galaxy number density fluctuations is

$$\begin{eqnarray}&&{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}\,:=\,\left\langle {{\rm{\Delta }}}_{i,{\ell }m}\,{{\rm{\Delta }}}_{j,{\ell }m}^{* }\right\rangle \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{2}{\pi }\,\int {\text{}}{dk}\,{k}^{2}\,{P}_{\mathrm{lin}}(k){{\rm{\Delta }}}_{i,{\ell }}(k){{\rm{\Delta }}}_{j,{\ell }}(k),\end{eqnarray} \tag{ 3 }$$

where k = ∣ k ∣ is the Fourier mode corresponding to the physical separation between a pair of galaxies, P_lin(k) := P_lin(k, z = 0) is the present-day linear matter power spectrum, and

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i,{\ell }}(k)=\int {\text{}}{dr}\,{n}_{i}(r)D(r)\left[b(r){j}_{{\ell }}({kr})-f(r)j{{\prime\prime} }_{{\ell }}({kr})\right].\end{eqnarray} \tag{ 4 }$$

In the equation above, D (normalized such that D = 1 at z = 0) is the growth factor of density perturbations, $f\,:=\,-{\text{}}d\mathrm{ln}D/{\text{}}d\mathrm{ln}(1+z)$ is the growth rate, j_ℓ is the ℓth-order spherical Bessel function, and primes denote derivation with respect to the argument of the function.

Here, we are interested in developing a new analysis framework for the harmonic-space power spectrum of galaxy clustering, which can allow for a direct comparison to the results obtained with standard analyses of the three-dimensional Fourier-space power spectrum and real-space two-point correlation function. To do so, we adopt the common approach followed in these studies, in which constraints are obtained primarily on the bias, the growth, and the so-called distortion parameters related to the Alcock–Paczynski test (see, e.g., Euclid Collaboration et al. 2020; Section 3.2.1). However, it should be noted that, in a harmonic-space analysis, we do not have the need to translate angles and redshifts into physical distances; we therefore are insensitive to the Alcock–Paczynski effect. Let us hence focus on expressing Equation (3) in terms of bias and growth.

First of all, it is useful to introduce the following two derived quantities:

$$\begin{eqnarray}&&b{\sigma }_{8}(z)\,:=\,b(z)D(z){\sigma }_{8},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&f{\sigma }_{8}(z)\,:=\,f(z)D(z){\sigma }_{8},\end{eqnarray} \tag{ 6 }$$

where ${\sigma }_{8}^{2}$ is the rms variance of matter density fluctuations on spheres of 8 h⁻¹ Mpc radius; it effectively acts as a normalization of P_lin(k) at redshift zero. The definitions of Equations (5) and 6 come from the consideration that, in a given redshift bin centered on $\bar{z}$, the three quantities b, f, and σ₈ are indistinguishable from each other in a measurement of the observed Fourier-space galaxy power spectrum

$$\begin{eqnarray}&&{P}_{{\rm{\Delta }}{\rm{\Delta }}}(k,\mu ;\bar{z})\simeq {\left[b(\bar{z})+f(\bar{z}){\mu }^{2}\right]}^{2}\,{D}^{2}(\bar{z}){P}_{\mathrm{lin}}(k),\end{eqnarray} \tag{ 7 }$$

where $\mu =\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{k}}/k$.

To recast Equation (3) in terms of b σ₈ and f σ₈, we require that these two quantities vary slowly across the redshift bin. Such a condition can be easily satisfied when dealing with spectroscopic redshift accuracy, which allows (in theory) for almost arbitrarily narrow redshift bins. Then, we factorize b σ₈ and f σ₈ out of the integral in Equation (4) and thus write

$$\begin{eqnarray}&&{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}\simeq b{\sigma }_{8}^{i}\,b{\sigma }_{8}^{j}\,{T}_{{ij},{\ell }}^{\delta \delta }+f{\sigma }_{8}^{i}\,f{\sigma }_{8}^{j}\,{T}_{{ij},{\ell }}^{{VV}}-2\,b{\sigma }_{8}^{(i}\,f{\sigma }_{8}^{j)}\,{T}_{{ij},{\ell }}^{\delta V},\end{eqnarray} \tag{ 8 }$$

where, for a generic quantity X(z), we call Xⁱ the value at the central redshift of the ith bin, round brackets around indices denote symmetrization, and we define the templates

$$\begin{eqnarray}&&{T}_{{ij},{\ell }}^{\delta \delta }:=\displaystyle \frac{2}{\pi }\,\int {\text{}}{dk}\,{\text{}}{dr}\,{\text{}}d\tilde{r}\,{k}^{2}\,\displaystyle \frac{{P}_{\mathrm{lin}}(k)}{{\sigma }_{8}^{2}}\,{n}_{i}(r){n}_{j}(\tilde{r}){j}_{{\ell }}({kr}){j}_{{\ell }}(k\tilde{r}),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{T}_{{ij},{\ell }}^{{VV}}:=\displaystyle \frac{2}{\pi }\,\int {\text{}}{dk}\,{\text{}}{dr}\,{\text{}}d\tilde{r}\,{k}^{2}\,\displaystyle \frac{{P}_{\mathrm{lin}}(k)}{{\sigma }_{8}^{2}}\,{n}_{i}(r){n}_{j}(\tilde{r})j{{\prime\prime} }_{{\ell }}({kr})j{{\prime\prime} }_{{\ell }}(k\tilde{r}),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{T}_{{ij},{\ell }}^{\delta V}:=\displaystyle \frac{2}{\pi }\,\int {\text{}}{dk}\,{\text{}}{dr}\,{\text{}}d\tilde{r}\,{k}^{2}\,\displaystyle \frac{{P}_{\mathrm{lin}}(k)}{{\sigma }_{8}^{2}}\,{n}_{(i}(r){n}_{j)}(\tilde{r}){j}_{{\ell }}({kr})j{{\prime\prime} }_{{\ell }}(k\tilde{r}).\end{eqnarray} \tag{ 11 }$$

Operatively, to compute the spectra defined above, we assume a Planck cosmology with parameters (Ade et al. 2016): Ω_c = 0.2603, Ω_b = 0.0484, H₀ = 67.74 km s⁻¹ Mpc⁻¹, σ₈ = 0.8301, τ = 0.667, and n_s = 0.9667, which are, respectively, the cold dark matter abundance today, the baryon abundance today, the Hubble constant, the clustering amplitude, the optical depth to reionization, and the slope of the primordial power spectrum of curvature perturbations.

Moreover, we also need to specify the input redshift distribution of sources. For this, we consider both the LOWZ and CMASS spectroscopic galaxy samples of BOSS DR12 (for details, see Section 2.2). Top panel of Figure 1 shows the two observed values of n(z) and the binning we adopt for LOWZ (blue hues) and CMASS (pink hues). We construct the top-hat bins as

$$\begin{eqnarray}&&{n}_{i}(z)=n(z)\times \displaystyle \frac{1}{2}\left[1-\tanh \left(\displaystyle \frac{2\,| z-{\bar{z}}_{i}| -{\rm{\Delta }}{z}_{i}}{\sigma \,{\rm{\Delta }}{z}_{i}}\right)\right],\end{eqnarray} \tag{ 12 }$$

where Δz_i is the bin width, ${\bar{z}}_{i}$ the ith bin center, and σ a smearing for the bin edges necessary to ensure numerical stability of the integration. For this, we choose σ = 0.002, and we have checked that the specific choice does not impact the results.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 1, we present the three factorized spectra of Equation (9)–(11), modulated by their corresponding parameters as in the right-hand side of Equation (8), thus accounting for the density, RSD, and density–RSD contributions to the total signal (left-hand side of Equation (8)). As a benchmark, we choose the CMASS auto-bin correlation with i = j = 5, and we plot the ratio between the factorized spectra and the total of ${S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}$. As expected, the main contribution comes from matter density fluctuations, with RSD being important only on very large scales. However, the cross-correlation between density and RSD is a significantly non-negligible term, especially at ℓ ≤ 50, where it contributes to the total signal between 20% and 40%.

To test the validity of our factorization, in Figure 2 we demonstrate the performance of the approximated expansion described in Equation (8) against the correct output of Equation (3), for LOWZ in the top panel and CMASS in the bottom panel. (Redshift-bin ranges can be found in the second column of Table 1). The factorized spectra are calculated by modifying the publicly available code CLASS (Lesgourgues 2011). It is evident that our approximation is in excellent agreement with the exact output at the subpercent level.

Zoom In Zoom Out Reset image size

Table 1. Redshift Range, Number of Sources, Shot Noise, and Maximum Multipole for the Redshift Bins Considered in the Analysis

Bin ID	$({z}_{\min },{z}_{\max })$	# of Gals	Shot Noise (sr)	${{\ell }}_{\max }^{i}$
LOWZ-1	[0.150, 0.190]	33, 906	8.08 × 10−5	49
LOWZ-2	[0.190, 0.230]	38, 728	7.07 × 10−5	60
LOWZ-3	[0.230, 0.270]	44, 291	6.18 × 10−5	71
LOWZ-4	[0.270, 0.310]	51, 781	5.27 × 10−5	81
LOWZ-5	[0.310, 0.350]	70, 879	3.86 × 10−5	91
LOWZ-6	[0.350, 0.390]	68, 701	3.98 × 10−5	101
LOWZ-7	[0.390, 0.430]	53, 191	5.15 × 10−5	111
CMASS-1	[0.450, 0.477]	81, 836	3.71 × 10−5	123
CMASS-2	[0.477, 0.504]	112, 589	2.67 × 10−5	129
CMASS-3	[0.504, 0.531]	118, 915	2.55 × 10−5	135
CMASS-4	[0.531, 0.558]	113, 139	2.68 × 10−5	141
CMASS-5	[0.558, 0.585]	99, 672	3.04 × 10−5	147
CMASS-6	[0.585, 0.612]	80, 914	3.75 × 10−5	153
CMASS-7	[0.612, 0.639]	61, 551	4.93 × 10−5	159
CMASS-8	[0.639, 0.670]	44, 057	6.89 × 10−5	164

Download table as:  ASCIITypeset image

Finally, we emphasize that Equations (3) and (4) follow directly from Equation (1) and are therefore correct only at first order in cosmological perturbation theory. Therefore, we restrict our analysis to strictly linear scales and fix the maximum wavenumber valid for linear theory to ${k}_{\max }=0.1\,h\,{\mathrm{Mpc}}^{-1}$, independent of redshift. Then, we implement a conservative redshift-dependent maximum multipole ℓfor each bin given by ${{\ell }}_{\max }^{i}={k}_{\max }\,r({\bar{z}}_{i})$. As for the minimum multipole, we define it as ${{\ell }}_{\min }=\pi /(2{f}_{\mathrm{sky}})$, where f_sky is the survey sky coverage. In our specific case of f_sky = 0.2081 for LOWZ and f_sky = 0.2416 (see Section 2.2), this corresponds to a common ${{\ell }}_{\min }=7$.

As already mentioned, in this study we use LOWZ and CMASS spectroscopic galaxies from the DR12 of BOSS. The LOWZ sample consists of Luminous Red Galaxies (LRGs) in the redshift range 0.1 < z < 0.45, while CMASS galaxies are LRGs at higher redshift, 0.45 < z < 0.8. Then, we have applied further redshift cuts on both samples, namely 0.15 ≤ z ≤ 0.43 for LOWZ and 0.45 ≤ z ≤ 0.67 for CMASS. This choice is justified as follows. First, we ensure that the redshift ranges of the two samples do not overlap. Second, we make sure that the available mock catalogs from BOSS DR12 cover the same redshift range with our data selection. These mocks, as we will explain in detail in Section 3.2.3, are essential for the investigation of systematic effects and for checking the pipeline's internal consistency. Finally, we avoid the inclusion of sources at z > 0.67 for CMASS so that we can ignore the small number of sources at that redshift range, which would introduce a considerable Poisson noise far surpassing the measured signal.

Let us now describe the procedure followed to construct the masks and the galaxy overdensity maps for LOWZ and CMASS. All details are summarized in Table 1.

Besides the galaxy catalogs, the BOSS collaboration has made publicly available random catalogs that are 50 times denser than the galaxy ones. Those were constructed after taking into account the completeness and the veto masks. The completeness masks indicate to what extent the observations under consideration are complete, while the veto masks account for observational effects such as bright objects and stars, extinction and seeing cuts, fiber collisions, fiber centerposts, and others. We built the final binary mask based on the random catalogs, assigning zero where there are no observations and one otherwise, using the HEALPix pixelization scheme (Gorski et al. 2005) with N_side = 1024. We deem this resolution more than sufficient for the scales of interest in this work (last column of Table 1), given that the largest scale that can be safely resolved by a HEALPix map is ${{\ell }}_{\max }=2{N}_{\mathrm{side}}$ (Ando et al. 2017).

Then, we proceed with the the construction of the galaxy overdensity maps (where the mask has a value of one) for the seven redshift bins of LOWZ and the eight bins of CMASS as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{p}=\displaystyle \frac{{n}_{{\rm{g}},p}-\langle {n}_{{\rm{g}}}{\rangle }_{p}}{\langle {n}_{{\rm{g}}}{\rangle }_{p}},\end{eqnarray} \tag{ 13 }$$

where n_g,p is the number of weighted galaxies in a given pixel p and 〈n_g〉_p is the mean number of weighted galaxies per pixel. We should note that the BOSS systematic weights (Reid et al. 2015) are not included in our analysis because we checked that they do not impact our final results.

The resulting maps are shown in Figure 3, where for simplicity we present only the cumulative galaxy angular distribution for each sample. As previously stated, the total sky areas correspond to f_sky = 0.2081 for LOWZ and f_sky = 0.2416 for CMASS, excluding the masked regions (shown in gray).

Zoom In Zoom Out Reset image size

The harmonic-space tomographic power spectrum can be estimated from the galaxy overdensity maps as ⁵

$$\begin{eqnarray}&&\widehat{{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}}=\displaystyle \frac{1}{2{\ell }+1}\,\displaystyle \sum _{m=-{\ell }}^{{\ell }}{{\rm{\Delta }}}_{i,{\ell }m}\,{{\rm{\Delta }}}_{j,{\ell }m}-\displaystyle \frac{{\delta }_{{ij}}^{({\rm{K}})}}{{\bar{n}}_{{\rm{g}}}^{i}},\end{eqnarray} \tag{ 15 }$$

where Δ_{i,ℓ m} are the harmonic coefficients of the pixelized overdensity map of Equation (13), and we have subtracted the shot-noise term (fourth column of Table 1), which is diagonal in i − j, with δ^(K) being the Kronecker symbol.

To relate the estimated power spectrum to the underlying one, we need to account for the mask, which introduces a coupling between different multipoles. Because the unmasked field is statistically isotropic, it is related to the measured one via

$$\begin{eqnarray}&&\left\langle \widehat{{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}}\right\rangle =\displaystyle \sum _{{{\ell }}^{{\prime} }}{{\mathsf{R}}}_{{\ell }{{\ell }}^{{\prime} }}\,\widehat{{S}_{{ij},{{\ell }}^{{\prime} }}^{{\rm{\Delta }}{\rm{\Delta }}}},\end{eqnarray} \tag{ 16 }$$

where ${{\mathsf{R}}}_{{\ell }{{\ell }}^{{\prime} }}$ is the coupling matrix. It is defined as

$$\begin{eqnarray}{{\mathsf{R}}}_{{\ell }{{\ell }}^{{\prime} }}=\displaystyle \frac{2{{\ell }}^{{\prime} }+1}{4\pi }\,\displaystyle \sum _{{\ell }^{\prime\prime} }(2{\ell }^{\prime\prime} +1){W}_{{\ell }^{\prime\prime} }{\left(\begin{array}{ccc}{\ell } & {{\ell }}^{{\prime} } & {\ell }^{\prime\prime} \\ 0 & 0 & 0\end{array}\right)}^{2},\end{eqnarray} \tag{ 17 }$$

with the matrix in square parentheses being the Wigner 3-j symbol and W_ℓ'' the harmonic-space power spectrum of the mask. The latter reads

$$\begin{eqnarray}&&{W}_{{\ell }}=\displaystyle \frac{1}{2{\ell }+1}\,\displaystyle \sum _{m=-{\ell }}^{{\ell }}{\left|{v}_{{\ell }m}^{}\right|}^{2},\end{eqnarray} \tag{ 18 }$$

where v_{ℓ m} are the coefficients of the harmonic decomposition of the binary mask $v(\hat{{\boldsymbol{r}}})$. The coupling matrices for the LOWZ and CMASS masks are shown in the two panels of Figure 4.

Zoom In Zoom Out Reset image size

Generally, the direct inversion of the coupling matrix is not possible, because the loss of information due to the mask makes the coupling matrix singular and therefore the inversion is ill-conditioned. One way to overcome this problem is to introduce bandpowers, which we do by utilizing the public code pymaster (Alonso et al. 2019). Our bandpowers, s, are a set of eight multipoles with weights ${w}_{s}^{{\ell }}$ normalized such that ${\sum }_{{\ell }\in s}{w}_{s}^{{\ell }}=1$. Now, for the coupled pseudo-C_ℓ in the sth bandpower, we define

$$\begin{eqnarray}&&\widehat{{{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}}=\displaystyle \sum _{{\ell }\in s}{w}_{s}^{{\ell }}\,\langle \widehat{{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}}\rangle ,\end{eqnarray} \tag{ 19 }$$

whose expectation value is then

$$\begin{eqnarray}&&\left\langle \widehat{{{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}}\right\rangle =\displaystyle \sum _{{\ell }\in s}\,{w}_{s}^{{\ell }}\,\displaystyle \sum _{{{\ell }}^{{\prime} }}{{\mathsf{R}}}_{{\ell }{{\ell }}^{{\prime} }}\,\widehat{{S}_{{ij},{{\ell }}^{{\prime} }}^{{\rm{\Delta }}{\rm{\Delta }}}}.\end{eqnarray} \tag{ 20 }$$

By doing so, we implicitly assume that the true power spectrum is also a stepwise function, which is related to the bandpowers via

$$\begin{eqnarray}&&\widehat{{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}}=\displaystyle \sum _{s}\widehat{{{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}}\,{\rm{\Theta }}({\ell }\in s),\end{eqnarray} \tag{ 21 }$$

with Θ the Heaviside step function.

Finally, after inserting the above relation in Equation (19), we can derive the unbiased estimator

$$\begin{eqnarray}&&{{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}=\displaystyle \sum _{{s}^{{\prime} }}{{\mathsf{M}}}_{{{ss}}^{{\prime} }}^{-1}\,\widehat{{{/}\!\!\!\!{S}}_{{ij},{s}^{{\prime} }}^{{\rm{\Delta }}{\rm{\Delta }}}},\end{eqnarray} \tag{ 22 }$$

where M is the binned coupling matrix, i.e.,

$$\begin{eqnarray}&&{{\mathsf{M}}}_{{{ss}}^{{\prime} }}=\displaystyle \sum _{{\ell }\in s}\displaystyle \sum _{{{\ell }}^{{\prime} }\in {s}^{{\prime} }}{w}_{s}^{{\ell }}\,{{\mathsf{R}}}_{{\ell }{{\ell }}^{{\prime} }}.\end{eqnarray} \tag{ 23 }$$

Finally, the theory spectra ${S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}$ should also be binned according to the bandpower choice in order to be compared with the data (for in general the theory curve is not a stepwise function). This translates into constructing

$$\begin{eqnarray}&&{{/}\!\!\!\!{S}}_{{ij},s}^{\mathrm{gg},\mathrm{th}}=\displaystyle \sum _{{s}^{{\prime} }}{{\mathsf{M}}}_{{{ss}}^{{\prime} }}^{-1}\,\displaystyle \sum _{{{\ell }}^{{\prime} }\in {s}^{{\prime} }}{w}_{{s}^{{\prime} }}^{{{\ell }}^{{\prime} }}\,{{\mathsf{R}}}_{{{\ell }}^{{\prime} }{\ell }}\,{S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}},\end{eqnarray} \tag{ 24 }$$

where we remind the reader that our theoretical model for ${S}_{{ij},{\ell }}^{{\rm{\Delta }}{\rm{\Delta }}}$ is defined in Equation (8).

Here, we shall describe the various approaches we have followed to estimate the data covariance matrix.

3.2.1. Theory Covariance

The simplest scenario is to assume that the data covariance matrix is Gaussian, which then reads

$$\begin{eqnarray}&&{{\mathsf{C}}}_{{{ss}}^{{\prime} }}^{{ij},{kl}}=\displaystyle \frac{{\delta }_{{{ss}}^{{\prime} }}^{({\rm{K}})}}{(2s+1){\rm{\Delta }}s\,{f}_{\mathrm{sky}}}\left({{/}\!\!\!\!{C}}_{{ik},s}^{{\rm{\Delta }}{\rm{\Delta }}}\,{{/}\!\!\!\!{C}}_{{jl},s}^{{\rm{\Delta }}{\rm{\Delta }}}+{{/}\!\!\!\!{C}}_{{il},s}^{{\rm{\Delta }}{\rm{\Delta }}}\,{{/}\!\!\!\!{C}}_{{jk},s}^{{\rm{\Delta }}{\rm{\Delta }}}\right),\end{eqnarray} \tag{ 25 }$$

with Δs = 1 being the multipole range in our bandwidth binning and

$$\begin{eqnarray}&&{{/}\!\!\!\!{C}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}={{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}+\displaystyle \frac{{\delta }_{{ij}}^{({\rm{K}})}}{{\bar{n}}_{{\rm{g}}}^{i}}.\end{eqnarray} \tag{ 26 }$$

3.2.2.  PolSpice Covariance

The Gaussian covariance is, by construction, diagonal in multipoles, which is then reflected by the structure of Equation (25)—i.e., by the presence of the Kronecker symbol. However, as we have described above when introducing pseudo-C_ℓ s, convolution with the survey mask induces a coupling between modes. This, in turn, reflects onto the covariance matrix, and a way to estimate it is provided by the PolSpice package (Chon et al. 2004). This is described in detail in Efstathiou (2004), and the corresponding covariance matrix has following form:

$$\begin{eqnarray}&&{{\mathsf{C}}}_{{{ss}}^{{\prime} }}={{\mathsf{M}}}_{{{ss}}^{{\prime} }}^{-1}\,{{\mathsf{V}}}_{{{ss}}^{{\prime} }}{\left[{\left({{\mathsf{M}}}^{-1}\right)}_{{{ss}}^{{\prime} }}\right]}^{{\mathsf{T}}},\end{eqnarray} \tag{ 27 }$$

with

$$\begin{eqnarray}&&{{\mathsf{V}}}_{{{ss}}^{{\prime} }}=\displaystyle \frac{2\,{{\boldsymbol{d}}}_{s}\,{{\boldsymbol{d}}}_{{s}^{{\prime} }}\,{{\mathsf{M}}}_{{{ss}}^{{\prime} }}}{2\,{s}^{{\prime} }+1}.\end{eqnarray} \tag{ 28 }$$

We remind the reader that ${{\mathsf{M}}}_{{{ss}}^{{\prime} }}$ is the binned coupling matrix, given in Equation (23), and d _s is the data vector corresponding to bandpower s.

3.2.3. Mocks Covariance

The most realistic way to obtain a realistic data covariance matrix is through simulating the data themselves and extracting their covariances from the sample. The BOSS collaboration provides us with simulated galaxy catalogs, namely the PATCHY (Kitaura et al. 2016) and the QPM (White et al. 2013; Wang et al. 2017) mocks. They are constructed assuming a fixed cosmological model and can be used to estimate accurate covariance matrices because they include various observational and systematic effects. In this work, we have selected the QPM mocks to calculate the covariance matrix for the LOWZ sample and the PATCHY mocks for the CMASS sample. This decision is made due to the different redshift ranges covered by the simulated catalogs. Indeed, the PATCHY mocks contain galaxies at redshifts 0.2 < z < 0.75, while our selected LOWZ sample starts from z > 0.15. The selection of two different sets of mocks is a further test to validate our analysis on top of the consistency checks that we make in Appendix 5.

The resulting covariance matrix is then

$$\begin{eqnarray}&&{{\mathsf{C}}}_{{{ss}}^{{\prime} }}=\displaystyle \frac{1}{{N}_{{\rm{m}}}-1}\,\displaystyle \sum _{m=1}^{{N}_{{\rm{m}}}}\left({{\boldsymbol{d}}}_{s}^{m}-{\bar{{\boldsymbol{d}}}}_{s}\right){\left({{\boldsymbol{d}}}_{{s}^{{\prime} }}^{m}-{\bar{{\boldsymbol{d}}}}_{{s}^{{\prime} }}\right)}^{{\mathsf{T}}},\end{eqnarray} \tag{ 29 }$$

with m = 1...N_m running through the N_m = 1000(2048) available QPM(PATCHY) mocks, with ${{\boldsymbol{d}}}_{s}^{m}=\{{{/}\!\!\!\!{{\boldsymbol{S}}}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }},m}\}$ being the mock data vector, and

$$\begin{eqnarray}&&{\bar{{\boldsymbol{d}}}}_{s}^{m}=\displaystyle \frac{1}{{N}_{{\rm{m}}}}\,\displaystyle \sum _{m=1}^{{N}_{{\rm{m}}}}{{\boldsymbol{d}}}_{s}^{m}.\end{eqnarray} \tag{ 30 }$$

Some extra care has to be taken when the covariance matrix is estimated using simulations, and the final likelihood should be corrected accordingly. This is because the inverse of the covariance matrix derived from simulations can be a biased estimator of the inverse of the true covariance matrix (Hartlap et al. 2007). Here, we present two methods that account for this correction:

• 1.  
  The former was first proposed by Hartlap et al. (2007), and it is based on the following reasoning. While the covariance matrix inferred from simulations can be an unbiased estimator of the true covariance C, its inverse (entering the χ² in Equation (14)) is not, and therefore should be rescaled according to Parvin (2004)

  $$\begin{eqnarray}&&{{\mathsf{C}}}^{-1}\to \displaystyle \frac{{N}_{{\rm{m}}}-{N}_{{\rm{d}}}-2}{{N}_{{\rm{m}}}-2}\,{{\mathsf{C}}}^{-1}.\end{eqnarray} \tag{ 31 }$$

  By doing so, the assumption of a Gaussian likelihood can be maintained.
• 2.  
  The latter has been proposed by Sellentin & Heavens (2015), where the authors state that the Gaussian likelihood should be replaced by the Student's t-distribution,

  $$\begin{eqnarray}&&{ \mathcal L }\propto {\left[1+\displaystyle \frac{{\chi }^{2}}{{N}_{{\rm{m}}}-1}\right]}^{-{N}_{{\rm{m}}}/2},\end{eqnarray} \tag{ 32 }$$

  replacing as well the Gaussian covariance in the χ² of Equation (14) with Equation (29). The proportionality in Equation (33) is set by

  $$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Gamma }}\left({N}_{{\rm{m}}}/2\right){[\det ({\mathsf{C}})]}^{-1/2}}{{[\pi ({N}_{{\rm{m}}}-1)]}^{{N}_{{\rm{d}}}/2}\,{\rm{\Gamma }}\left[({N}_{{\rm{m}}}-{N}_{{\rm{d}}})/2\right]},\end{eqnarray} \tag{ 33 }$$

  where Γ is the Euler Gamma function, and assuming N_m > N_d.

We note that the results after implementing both methods are equivalent, and therefore we present results of the latter correction alone.

As a first test, we want to replace the data vector d _s with a mock data vector ${{\boldsymbol{d}}}_{s}^{m}$ randomly chosen from the QPM or PATCHY sets of available mocks for LOWZ or CMASS. If our pipeline does not suffer from systematic effects, the analysis using the mock vector is expected to produce results with constraining power comparable to that of the same samples using the real data vector, but the points will be scattered differently around the theory prediction. As shown in Figure 6, this is indeed the case, with the mock data vector performing similarly to the real data. In particular, values with the real data results for the b σ₈ (blue for mock data in Figure 6 and blue for real data in the right panel of Figure 7) and for f σ₈ (red for mock data in Figure 6 and red for real data in the right panel of Figure 7) all randomly scattered around the ΛCDM prediction (solid gray line).

As discussed in Section 2, our method allows for almost model-independent measurements of the bias and growth of galaxies directly in harmonic space. However, we still have to assume a cosmology to compute the various ingredients of Equation (9)–(11).

In order to validate our working hypothesis, we here change the underlying fiducial cosmology that is assumed for the theory spectra ${T}_{{ij},{\ell }}^{\delta \delta }$, ${T}_{{ij},{\ell }}^{{VV}}$, and ${T}_{{ij},{\ell }}^{\delta V}$ calculated in CLASS. In particular, we first vary the cosmological parameters at the edge of their 95% C.L. for Planck ΛCDM best-fit values of Ade et al. (2016), namely Ω_b = 0.0621, Ω_c = 0.2718, and H₀ = 68.67 km s⁻¹ Mpc⁻¹. Furthermore, we perform the analysis assuming the best-fit values of the Planck wCDM cosmology of Ade et al. (2016) with Ω_b = 0.04802, Ω_c = 0.2568, H₀ = 68.1 km s⁻¹ Mpc⁻¹, and w = −1.019. Finally, we consider a more significantly different cosmology: the ΛCDM best-fit from WMAP7 (Komatsu et al. 2011), with Ω_b = 0.045, Ωc = 0.227, and H₀ = 70.4 km s⁻¹ Mpc⁻¹.

Constraints for these three cosmologies are respectively shown in the left panel of Figure 7 in cyan, green, and brown for LOWZ, and in gold, purple, and pink for CMASS. The reconstructed values of $b{\sigma }_{8}^{i}$ and $f{\sigma }_{8}^{i}$ are well within the 68% C.L. intervals among all cosmologies (shown with blue and red), implying that the majority of the cosmological information is contained in b σ₈(z) and f σ₈(z) themselves.

A further test for our pipeline is to check the effect of neglecting cross-bin terms in the data vector, i.e., ${{/}\!\!\!\!{S}}_{{ij},s}^{{\rm{\Delta }}{\rm{\Delta }}}$ with i ≠ j, and the corresponding entries in the covariance matrix—namely, terms like {i − j, i − j}, {i − j, j − j}, {j − i, i − j}, etc. Basically, we use the full covariance matrix as defined in Equation (25). Marginalized constraints presented in the right panel of Figure 7, illustrated in brown for LOWZ and pink for CMASS, respectively, show a satisfactory consistency with the default case where we do not consider the cross-bin-no-pair correlations at all.

Finally, we perform the following test to validate our derived cosmological measurements (see also Loureiro et al. 2019). This test is described as follows. In the case that we have a diagonal covariance matrix, we are able to construct a vector

$$\begin{eqnarray}&&{\boldsymbol{R}}={{\mathsf{O}}}^{-1}\,[{\boldsymbol{d}}-{\boldsymbol{t}}],\end{eqnarray} \tag{ 34 }$$

which is the normalized residuals, where O is the diagonal matrix constituted of the square root of the variances, d is the data vector, and t is the theory vector. Following, e.g., Andrae et al. (2010), we note that the residuals are described by a standard normal distribution if the model given by t represents the data and the errors are known exactly. Nonetheless, the errors are usually estimated from the limited number of samples provided by the experiments, and the residual distribution is therefore given by the Student's t-distribution, which approaches the standardized Gaussian as the number of samples increases. Hence, if the residuals are described by a standard distribution, we can either infer that either we found the correct model or that the data are not good enough to show any model preference. On the other hand, is the residual distribution deviates from Gaussianity, then the model is ruled out.

In the realistic scenario, we know that the errors are correlated and the covariance matrix is no longer diagonal. However, we can diagonalize it as

$$\begin{eqnarray}&&{\mathsf{C}}={\mathsf{Q}}\,{\mathsf{Y}}\,{{\mathsf{Q}}}^{-1},\end{eqnarray} \tag{ 35 }$$

where Q is the eigenvector matrix and Y is a diagonal matrix constructed from the eigenvalues of the original covariance matrix, C. By treating the new diagonal covariance matrix as Y, the residual distribution now reads

$$\begin{eqnarray}&&{\boldsymbol{R}}={{\mathsf{O}}}^{-1}\,{{\mathsf{Q}}}^{-1}\,[{\boldsymbol{d}}-{\boldsymbol{t}}],\end{eqnarray} \tag{ 36 }$$

where O contains the square roots of the eigenvalues.

In Figure 8, we show the normalized residuals for LOWZ (left) and CMASS (right) after we input the mean values of the constrained parameter set $\{{\theta }_{\alpha }\}=\{b{\sigma }_{8}^{i},f{\sigma }_{8}^{i}\}$. By considering a Kolmogorov–Smirnov test, we accept the null hypothesis that both histograms are consistent with being derived from a standard normal distribution, since the p-values are 0.18 and 0.13, respectively, at the 95% significance level, and we conclude that our model represents the data well.