Unified galaxy power spectrum measurements from 6dFGS, BOSS, and eBOSS

We make use of recent developments in the analysis of galaxy redshift surveys to present an easy to use matrix-based analysis framework for the galaxy power spectrum multipoles, including wide-angle effects and the survey window function. We employ this framework to derive the deconvolved power spectrum multipoles of 6dFGS DR3, BOSS DR12 and the eBOSS DR16 quasar sample. As an alternative to the standard analysis, the deconvolved power spectrum multipoles can be used to perform a data analysis agnostic of survey specific aspects, like the window function. We show that in the case of the BOSS dataset, the Baryon Acoustic Oscillation (BAO) analysis using the deconvolved power spectra results in the same likelihood as the standard analysis. To facilitate the analysis based on both the convolved and deconvolved power spectrum measurements, we provide the window function matrices, wide-angle matrices, covariance matrices and the power spectrum multipole measurements for the datasets mentioned above. Together with this paper we publish a \code{Python}-based toolbox to calculate the different analysis components. The appendix contains a detailed user guide with examples for how a cosmological analysis of these datasets could be implemented. We hope that our work makes the analysis of galaxy survey datasets more accessible to the wider cosmology community.


Introduction
In the last two decades the analysis of galaxy redshift survey datasets has become one of the most powerful tools to constrain cosmological models [3,4]. The next generation of galaxy redshift surveys, such as the Dark Energy Spectroscopic Instrument (DESI [5]) and the Euclid space mission [6], aim to observe 20 to 40 million galaxy redshifts, increasing the largest current dataset (BOSS [7], 1 million galaxies) by more than an order of magnitude. While observations of the Cosmic Microwave Background (CMB) have reached the sample variance limit for some observables, galaxy surveys are far away from this limit. Extracting cosmological information from galaxy surveys does pose significant challenges. Most modes which can be accessed by galaxy surveys are contaminated by non-linear contributions along with redshift-space distortions and a complicated relation between galaxy density and matter density. Many different avenues are currently investigated to tackle these challenges (e.g. [8][9][10][11]). Here we will ignore these issues, but focus on the analysis formalism itself.
Extracting cosmological information from the galaxy power spectrum comes with technical complications, such as wide-angle effects (e.g., [2,12,13]) and the survey window function, which convolves the measured power spectrum. While the survey window is usually known, this convolution can add significant modeling challenges. The current standard analysis derives the survey window function from the data and convolves any power spectrum model with the survey window before comparing it with the power spectrum measurement [1,14,15]. Instead of convolving the model we can also deconvolve the measured power spectrum. Past attempts of deconvolution were done on a mode by mode basis [16] or assumed the global plane parallel approximation (often called distant observer approximation [17][18][19]). Here we present a deconvolution formalism, which can be applied to a wide-angle multipole representation of the power spectrum. Recently, [20] advocated deconvolution in the context of going back to the quadratic estimator formalism originally derived in [21,22]. For these estimators one can naturally quote results convolved with a window, or deconvolved, or something in between [23,24]. Note that there is a technical difference between that formalism, which assumes the k bands in which the quadratic averaging is performed are the same as the k bands used to model the theory power, so that the window matrix is automatically square and invertible, and the formalism in this and other recent papers, in which the theory bands can have arbitrarily fine k resolution, independent of the observational averaging band width.
First we express the power spectrum analysis as two matrix multiplications, one matrix accounting for wide-angle effects and a second matrix accounting for the window function. We derive these matrices for some of the largest galaxy redshift surveys currently available (6dFGS DR3, BOSS DR12 and the eBOSS DR16 QSOs sample) and make these matrices available together with the power spectrum measurements 1 . The appendix of this paper provides a step-by-step guide for a galaxy survey analysis, including Python-based examples. The aim is to make galaxy redshift survey datasets more easily accessible for the wider cosmology community.
As an example application we perform a BAO analysis on the deconvolved power spectrum of BOSS DR12. We show that with the products derived in this paper, the likelihood derived from the convolved and deconvolved power spectrum multipoles is identical. This paper is organized as follows. We start with a review of the current formulation of the survey window function, which we turn into a matrix multiplication in section 2. In section 3 we derive a similar matrix accounting for wide-angle effects. In section 4 we deriving our deconvolution procedure. In section 5 we introduce the 6dFGS, BOSS and eBOSS datasets followed by an application of our matrix based deconvolution procedure in section 6. We conclude in section 7. Appendix A provides a detailed user guide with Pythonbased examples. In appendix B we discuss the correlation matrices needed for the window function calculations. We also show consistency between the equations used in this paper with the equations in [1] and [2] in appendix C. In appendix D we derive analytic equations for wide-angle effects. Finally, appendix E contains a summary of the convolved and deconvolved power spectrum measurements for the different datasets.
Throughout the paper we use a (fiducial) flat ΛCDM cosmology with Ω m = 0.31 when transferring observables (RA, DEC, z) into cartesian coordinates (x, y, z). When analysing mock datasets we use the cosmological parameters of the underlying simulations, which are different for each of the galaxy samples studied here (see section 5). The Python code used to calculate the wide-angle and window function matrices discussed in this paper is available at https://github.com/fbeutler/pk_tools.

The survey window function
Most studies involving the galaxy power spectrum include the survey window function by convolving the model for the power spectrum before comparing it to the measurement. Any asymmetry in the survey window distributes power between the multipoles. This can bias measurements of anisotropic observables like redshift-space distortions, if not taken into account correctly. The convolution of the power spectrum multipoles with the window function multipoles has been laid out in [1] and is given by where W (k, k ) describes the contribution of multipole to multipole due to the survey window 2 . The second term on the right hand side describes the integral constraint correction [25] and includes the 1D survey window function multipoles in Fourier space Q (k). While the treatment of the integral constraint in eq. (2.1) is a good approximation, [26] showed that there are additional contributions depending on how the reference (random) catalog has been generated. While [1] calculated W (k, k ) through pair counting (see eq. 33 of that reference), recently it was pointed out that this quantity can be obtained using a double Bessel integral [27] W (k, k ) = 2 π (−i) i ds s 2 j (ks)j (k s) where the matrices C L are given in section 3.1 of [27]. The 1D window function multipoles in configuration-space are defined as Q L (s) = (2L + 1) A dΩ s 4π ds 1nw (s 1 )n w (s + s 1 )L L (ŝ ·d) , (2.3) withn w (x) =n(x)w(x), wheren(x) is the local mean density of galaxies (i.e., expected number after selection effects), and w(x) represents a weight applied to the density measured at x (completeness and signal to noise weight). The normalization factor A here is the same number that appears in the power spectrum estimate, e.g., eq. (3) of [15]. However, we do not use the standard value of A given in eq. (13) of [15]. Instead we set A to the value necessary to enforce Q 0 (s → 0) ≡ 1. The same A value is then also used to normalize the power spectrum estimate, as required for consistency (as recently emphasized by [26]). We use this definition because we find it gives much closer to the desired unit normalized window, i.e., window where integration over theory k gives 1. This is desirable because it makes the convolved power close to the true power, i.e., just smeared slightly by the window instead of adding an offset in the overall normalization. Note that, at the continuum limit of eq. (2.3), our normalization convention can be written as However, the traditional definition as represented by eq. (13) of [15] is also intended to represent this equation! In the traditional calculation the integral over position and then(x) factors are represented by a sum over randoms timesn(z), wheren(z) is supposed to be a redshift dependent average ofn(x). This is an approximation that in practice only achieves the goal of Q 0 (s → 0) = 1 to ∼ 10%. We guarantee Q 0 (s → 0) = 1 by computing Q L (s) first and using it directly to fix A. For comparison with past results we give the relative normalization between the two methods in Table 1. To be clear, if the normalisation between the power spectrum and window function are consistent, it will not impact any likelihood analysis, so it is essentially cosmetic.
In the following section we will demonstrate how our window function formalism can be extended to include wide-angle effects before moving on to develop a matrix based convolution and deconvolution procedure.

Including wide-angle effects
Any measurement of the anisotropic power spectrum has to make a choice regarding the line-of-sight (LOS) direction for a galaxy pair (see, e.g., [2,12,13]). Based on this choice the triangle configuration between the observer and the galaxy pair is reduced to a separation amplitude s and a multipole dependent weighting given by the Legendre polynomials L . This geometric simplification introduces wide-angle effects in the power spectrum multipoles > 0 as well as the associated window function multipoles. Crucially, wide-angle effects can break the symmetry between the galaxy pair and hence can introduce odd multipoles, like the dipole and octopole. These wide-angle effects can be absorbed into the window function formalism described above. Based on [2] we can extend eq. (2.2) to get where the index n describes the order in the wide-angle expansion and the C (n) L matrices are given in appendix B. Eq. (2.5) includes new window function multipoles Q (n) L at each order of the wide-angle expansion n given by A derivation of eq. (2.5) is included in Appendix C.2. The constraints imposed by FFTs usually enforce the end-point LOS definitiond =ŝ 1 where the line of sight is chosen along the distance vector of one of the galaxies (see [2] for details). Eq. (2.5) can be calculated using one 1D Fourier transforms (FT) for each k bin with the complexity O(N k × N k log N k ), where N k and N k are the number of bins in k and k , respectively. Additionally we need to calculate the Q (k) using a 3D FT, which has the complexity O(N log N ), where N is the number of grid cells in which the random galaxies are binned (see appendix E.1 of [2] Mpc 2 ] Figure 1: The window function multipoles of the low redshift bin 0.2 < z < 0.5 of BOSS DR12, NGC, calculated following eq. (2.5) at zero order in the wide-angle expansion (n = 0). The corresponding window function multipoles at first order in the wide-angle expansion (n = 1) are shown in figure 2. The three columns correspond to the contributions of the three even multipoles to the five non-zero multipoles, including the dipole and octopole. Each window function is plotted multiple times as a function of k with fixed values of k (vertical dashed lines). This window function does not account for any bin averaging, but is evaluated at specific values k with 16 384 values of k (here we focus on 0 < k < 0.12 h Mpc −1 , while in practice we calculate the full range of 0 < k < 0.4 h Mpc −1 ). For visual purposes, each of the sub-panels is re-normalized by the factor written in the right hand corner. The final window function matrix used for the actual analysis requires bin averaging in both k o and k th (see eq. 2.16 for more details).
paper we will use eq. (2.5) to calculate the window function using N k = 16 384 and N k = 400 (from 0 < k < 0.4 h Mpc −1 in bins of ∆k = 0.001 h Mpc −1 ), which typically takes < 1 minute on a single state of the art laptop. Figure 1 shows the window function multipoles W (n) (k, k ) for the low redshift bin of BOSS DR12 in the North Galactic Cap (NGC). This figure only shows the window function multipoles at zero order in the wide-angle expansion (n = 0 in eq. 2.5). At this order, only even multipoles are generated, meaning we have a contribution from the monopole (left column), quadrupole (middle column) and hexadecapole (right column). Each window function is plotted several times with fixed k (dashed lines). These plots agree with the plots shown in Figure 7 of [1], where these multipoles have first been investigated, but they now include contributions to the odd multipoles. Figure 2 shows the corresponding window function multipoles at first order in the wideangle expansion (n = 1 in eq. 2.5), in which case we only have odd contributions from the dipole (left column) and octopole (right column). Also note that any window function multipole W (n) in which the combination + is a odd number, results in a complex quantity, Mpc 2 ] Figure 2: The window function multipoles of the low redshift bin 0.2 < z < 0.5 of BOSS DR12, NGC, calculated following eq. (2.5) at first order in the wide-angle expansion (n = 1).
The corresponding window function multipoles at zero order in the wide-angle expansion (n = 0) are shown in figure 1. The two columns correspond to the contributions of the two odd multipoles to the five non-zero multipoles. Each window function is plotted multiple times as a function of k with fixed values of k (vertical dashed lines). For visual purposes, each of the sub-panels is re-normalized by the factor written in the right hand corner.
to account for the fact that the estimated odd power spectrum multipoles are complex [2,28]. Based on the discussion above we can extend eq. (2.1) in section 2, which describes the convolution of the power spectrum multipoles with the window function multipoles by including wide-angle effects: (2.7) We now define a new multipole expansion including wide-angle effects as where d is the LOS vector and d = |d|. Based on this definition the power spectrum multipoles, P (n) (k), at different order in the wide-angle expansions n are given by [2,13] where ξ (n) are the configuration space multipoles at wide-angle order n.
For the remainder of this paper we will use the end-point LOS definition, meaning the LOS for each galaxy pair is oriented along the distance vector of one of the galaxies (d = s 1 ). In this case the zero order term (n = 0) is equivalent to the commonly used power spectrum multipoles and at this order only the even multipoles are present. At n = 1 only the odd multipoles are present and can be obtained as linear combinations of the even multipoles in configuration space The second order terms (n = 2) are given in eq. (2.16) -(2.18) of [2]. The analysis in this paper will be limited to n ≤ 1, even though including higher order terms (n > 1) is straightforward. Note that the dipole and octopole power spectra can be calculated directly in Fourier space without any need for a Hankel transform (see also eq. 3.56 of [29]): We included a derivation of these equations in appendix D. Using linear theory we obtain

Window function convolution as matrix multiplication
If we assume that we calculate our power spectrum model in bins of ∆k th and intend to compare this model with power spectrum measurements in bins of 'observed' bandpowers ∆k o , the window function needs to be integrated over k and k as where the step function Θ is defined as We will test the bin averaging employed in our analysis in the next subsection. Before we get to that, let us define our matrix nomenclature, starting with the vectors  where P true is a vector of model power spectrum multipoles and the convolution matrix is given by where each sub-matrix W (n) is of size N th × N o . Based on the definitions above, we can write the convolution of the power spectrum multipoles with the window function multipoles as a matrix multiplication where Q Q 0 (0) W 0 (0, k ) is a matrix representing the integral constraint correction. For the remainder of this paper we will absorb the integral constraint into W meaning we re-define This works exactly because the integral constraint effect is linear in the power spectrum just like the standard window function. In this paper the window function matrix W is generally defined in the k-range 0.0 < k < 0.4 h Mpc −1 , while the recommended k-range for any fit is 0.01 < k < 0.3 h Mpc −1 . The increased k-range accounts for the redistribution of power due to the window function, which connects modes inside the fitting range with modes beyond that range. The equations above include wide-angle effects at order n = 0 and n = 1, even though this formalism could easily be extended to n > 1. Figure 3 shows the matrix W after performing the integral in eq. (2.16) (and including the integral constraint correction).

Bin averaging the window function
In eq. (2.16) we wrote down how to take the window function (as derived in eq. 2.5 and plotted in figure 1 and 2) an05d account for the bin averaging to match the layout of the data vector. The default choices for our window functions are ∆k th = 0.001 h Mpc −1 and ∆k o = 0.01 h Mpc −1 . Here we will demonstrate how we implemented the bin averaging of the window function and we will test certain assumptions made in our formalism.

The observational binning ∆k o
All power spectrum multipole measurements provided with this paper use k-bins of ∆k o = 0.001 h Mpc −1 . However, for practical purposes 3 we use 10 times larger bins of ∆k o = 0.01 h Mpc −1 . The measured power spectra can easily be re-binned into larger bins using (2.23) One crucial point to highlight here is that the fundamental mode for all galaxy datasets included in this paper is smaller than our choice for ∆k o . This means we need to bin average our theoretical model in accordance with the average of eq. (2.22). For that reason we measure the window function in bins of ∆k o = 0.001 h Mpc −1 and average as

The theoretical binning ∆k th
The multiplication of the power spectrum model with the window function matrix effectively calculates a convolution. The question we want to address here is, what would be the required resolution for this convolution, meaning how sensitive is the convolved power spectrum to the choice of ∆k th ? In figure 5 we show convolved power spectrum multipoles using our default binning of ∆k o = 0.01 h Mpc −1 but different ∆k th . The reference model is based on the current standard method, where the power spectrum model is (1) Hankel transformed into configuration-space, (2) multiplied by the configuration-space window and (3) Hankel transformed back into Fourier space [14,30]. The result is than bin-averaged to ∆k o = 0.01 h Mpc −1 . The differences between our default choice of ∆k th = 0.001 h Mpc −1 and the reference model are < 1% of the measurement uncertainties for the monopole on most scales and even smaller for all other multipoles. Finally we note that the above binning tests for the power spectrum and window function are based on a ΛCDM model (with the default cosmology given at the end of the introduction). Our default binning choices might not be optimal for non-standard analysis that searches for variations in the power spectrum smaller than the chosen bin width. Examples for such cases are primordial features [31] or primordial non-Gaussianity [32]. However, our binning choice should be perfectly suitable for the standard RSD and BAO analysis.

Wide-angle effects as matrix multiplication
The odd power spectrum multipoles sourced by wide-angle effects are linearly related to the even multipoles through eq. (2.14) and (2.15). This implies that we can define a linear transformation between 'flat-sky' statistics (no wide-angle effects) and 'curved-sky' statistics (including wide-angle effects) as where P true,flat-sky = (P 0 , P 2 , P 4 ) is a vector of the predicted even multipoles and M is a 5N k th × 3N k th matrix, which transforms this vector into a new vector P true , which contains 5 multipoles, including the dipole and octopole. We can define where I is the identity matrix of size N th × N th and with Θ defined in eq. (2.17). As shown in these equations, we can implement derivatives within a matrix multiplication by including off-diagonal terms. Here we use two-sided derivatives except for the first and last bin in the data vector, where we use forward and backwards derivatives instead. The transformation matrix M does depend on the order of wide-angle effects n, since at n = 2 there are new even multipoles sourced by wide-angle effects. Following the rest of this paper, we only include terms at n ≤ 1. The plot on the right in figure 3 shows the matrix M for the low redshift bin of BOSS DR12 NGC. The only survey specific parameter in this matrix is the amplitude of the LOS vector d. The code to calculate the matrix M is publicly available 4 .

Deconvolution as matrix multiplication
Usually the aim of a power spectrum analysis is to obtain the likelihood where C −1 conv is the inverse covariance matrix of the (convolved) power spectrum multipoles, P conv o represents the measured power spectrum and P true is the unconvolved power spectrum model.
Using eq. (4.1) together with the maximum likelihood condition ∂χ 2 /∂P true = 0 implies that the power spectrum before convolution with the window function is related to the convolved power spectrum by   if W is a square matrix and non-singular. Note that the deconvolution process implied in eq. (4.2) cannot recover the true underlying power spectrum (which can be calculated in theory), since the window function did erase all information below the fundamental mode. The covariance matrix of the deconvolved power spectrum is given by The correlation matrix, R = C ij / C ii C jj , for BOSS DR12 NGC in the low redshift bin is plotted in figure 7 before (left) and after (right) deconvolution. Deconvolution often leads to anti-correlated bins [33], which can be intuitively explained by imagining a power spectrum model with very fine k binning (∆k fundamental mode k f ). If one moves the power spectrum in one bin up, but compensates by moving the power spectrum in the adjacent bin down, the convolved version of this power spectrum will look very similar to the original power spectrum, since the window function averages neighbouring bins. Hence such a mode of variation is not well constrained in a deconvolved estimate, leading to anti-correlation.
Generally however, there is a decrease in the correlation between bandpowers as clearly visible in figure 8, which implies an increase in the variance of modes within that bin, given by (4.5) Figure 6 shows the behaviour of W (n) (k, k) as a function of k for different ∆k in the low redshift bin of BOSS DR12. This clearly shows that deconvolution can significantly increase the uncertainties in the power spectrum multipoles if ∆k is small. However, this increase in the uncertainties does not reflect any loss of information, but only shows that for smaller ∆k the window function has a larger impact. The larger uncertainty is compensated by having less correlation (or anti-correlation) between the bandpowers. The information content of It should be noted, however, that our deconvolution procedure is based on a Gaussian likelihood (as given in eq. 4.1), which might not be true on the largest scales of the survey (see e.g. [34]). This could bias large-scale signals like primordial non-Gaussianity. Alternative to our procedure one could derive the deconvolved bandpowers using an MCMC approach without the need to assume a Gaussian likelihood or include a Gaussianisation step as e.g. suggested in [35].

Wide-angle compression
In eq. (4.2) we can replace W withW ≡ WM to compress our vector of 5 multipoles to a vector of 3 even multipoles, i.e., while the corresponding covariance matrix is This is equivalent to a minimum χ 2 fit of the band-parameterized model for P true,flat-sky to all multipoles. This is an over-constrained fit with the degrees of freedom equal to the number of bins in the odd multipoles 5 . Because the prediction of the odd multipoles from P true,flat-sky is model-independent, the reduced χ 2 (χ 2 /ν) for this fit is a model-independent test of the various quantities we report (C, M, etc.). Unfortunately, we find that when we try to fit all multipoles on all scales this way, χ 2 is unacceptably bad, i.e., the measured points are not consistent with any model, given our M, W, and C. As shown in figure 9 (left), this problem is driven by the high k bins in the odd multipoles. Figure 9 shows that χ 2 /ν is near unity if we drop the odd multipoles with k 0.1 h Mpc −1 , but quickly grows when including higher k. This is not a straightforward thing to understand, or even see, by looking at the measured points, covariance matrix, etc. It appears to be driven by a few special linear combinations of even and odd multipoles which have very low variance in the mocks used to compute the covariance matrix. These are then not well-predicted by 5 Assuming the theory is band-parameterized using the same set of k bins as the measurement, i.e., the number of free parameters is equal to the number of bins in the even power spectrum multipoles so that they cancel in the calculation of the degrees of freedom.

Survey/Sample
where w i includes any completeness weight as well as the FKP weight [36]. To calculate the FKP weight as well as the effective volume, we assume P 0 = 10 000h −3 Mpc 3 for BOSS and 6dFGS and P 0 = 6 000h −3 Mpc 3 for eBOSS. The effective volume in this table might differ from other references due to (1) a different fiducial cosmology and (2) only the raw (unweighted) number density of galaxies is used in our calculation. k f represents the fundamental mode, which is given by the wavenumber at which the monopole window Q 0 (k → 0) (column 6) or 95% (column 7). N m represents the number of mock realizations available for each sample, and the last column gives the ratio between our A definition, enforcing Q 0 (s → 0) = 1, and the normalization used in many previous papers (see section 2). For BOSS we also have catalogs after applying density field reconstruction. However, only 1000 of the 2048 BOSS mock catalogs have been processed in this way 6 . More information about the different samples can be found in the relevant data release papers of 6dFGS [37], BOSS [38] and eBOSS [39]. the theory as multiplied by M and W. E.g., for BOSS DR12 NGC, z3, k < 0.2 h Mpc −1 , we find χ 2 = 1400, of which 1127 is contributed by the single worst eigenvector of the covariance matrix. We were unable to find any flaw in the calculations that would fix this. Since we expect the wide-angle effects to be primarily important at low-k, our solution is to simply drop the odd multipoles at k > 0.1 h Mpc −1 from fits. The right panel of figure 9 shows an example of the impact of the lower k odd multipoles on the inferred deconvolved even multipole variance. As expected the impact of the odd multipoles is primarily at very low k.

Datasets
Here we will introduce the three galaxy redshift survey datasets we analyze in this paper, namely the 6dFGS DR3 sample, the BOSS DR12 sample and the eBOSS DR16 quasar sample.

6dFGS DR3
The 6-degree Field Galaxy Survey (6dFGS [37]) is a K-band selected, magnitude limited (K ≤ 12.9) galaxy survey, based on the 2MASS Extended Source Catalog (2MASS XSC; [40]). 6dFGS covers nearly the entire southern sky and is the lowest redshift sample included in this paper. The survey made use of the Six-Degree Field (6dF) multi-fibre instrument on the UK The curves in the plot on the right hand side are normalized so that the maximum of the curve is at unity (which allows comparison between the different surveys).
Schmidt telescope at the Siding Spring Observatory. The three data release papers [37,41,42] describe 6dFGS in full detail, including comparisons between 6dFGS, 2dFGRS and SDSS.
Mock catalogs: The 6dFGS mock catalogs are based on COLA simulations [49] with 1728 3 particles in boxes with 1.2h −1 Gpc on each side. The simulations use 20 time steps down to z = 0 and a mass resolution of 2.8 × 10 10 h −1 M . A friends-of-friends (FoF) finder is used to locate halos with a minimum of 32 dark matter particles per halo. The cosmology used in these simulations is Ω m = 0.3, Ω b = 0.0478, h = 0.68, σ 8 = 0.82 and n s = 0.96. The derived halo catalogs are populated with galaxies using an HOD measured on the dataset itself [43] (for more details see section 3 of [48]).
The 6dFGS sample has by far the smallest volume of all the datasets included in this analysis (see table 1 and figure 10). However, it does occupy a unique redshift range, which is not covered by BOSS or eBOSS.

BOSS DR12
The Baryon Oscillation Spectroscopic Survey (BOSS) was part of SDSS-III [7,50] and measured spectroscopic redshifts of 1 198 006 million galaxies [38]. The survey covers 10 252 deg 2 divided in two patches on the sky, the North Galactic Cap (NGC) and the South Galactic Cap (SGC), over a redshift range of 0.2 -0.75. Here we split this redshift range into two redshift bins defined by 0.2 < z < 0.5 and 0.5 < z < 0.75 with the effective redshifts z eff = 0.38 and 0.61, respectively. We also include the different incompleteness weights as which account for redshift failures (w rf ), fibre collisions (w fc ) and photometric systematics related to the observational seeing conditions and correlations with stellar density (w sys ) [51,52].
Mock catalogs: The BOSS collaboration provided mock catalogs for the final BOSS DR12 dataset (MD-Patchy mock catalogs [53]). These catalogs have been produced using approximate gravity solvers and analytical-statistical biasing models calibrated to a reference sample from the BigMultiDark simulations [54]. The BigMultiDark simulation is based on gadget-2 [55] with 3840 3 particles in a volume of (2.5h −1 Gpc) 3 assuming a ΛCDM cosmology with Ω m = 0.307115, Ω b = 0.048206, σ 8 = 0.8288, n s = 0.9611 and h = 0.6777. The mock catalogs use halo abundance matching to reproduce the observed BOSS two-and three-point clustering measurements [56]. This technique is applied as a function of redshift to reproduce the BOSS DR12 redshift evolution.
We already saw the BOSS DR12 NGC window function W (n) of the low redshift bin (z1) in figure 1 and 2. We also include the monopole window function, Q

eBOSS DR16 QSO
The extended Baryon Oscillation Spectroscopic Survey ( [60], eBOSS) is part of SDSS-IV [61] and relies on the same optical spectrographs [62] as the SDSS-III BOSS survey. In addition to observing luminous red galaxies (LRGs) and emission line galaxies (ELGs), eBOSS collected redshifts for ∼ 500 000 quasars. While the Quasar density is comparatively low, this sample has the distinction of covering the largest cosmic volume, leading to the smallest fundamental mode of all samples discussed in this paper (see figure 10). Here we focus on the eBOSS quasar sample to avoid any overlap with the other samples.
The eBOSS targets [63,64] are selected from the DR7 [65] and DR8 [66] photometric catalogs as well as the Wide Field Infrared Survey Explorer (WISE, [67]), as described in [68]. Just like BOSS, the eBOSS quasars are split in two angular regions, the North Galactic Cap (NGC) and South Galactic Cap (SGC). The effective areas of these regions are 2924 deg 2 and 1884 deg 2 , respectively (see Table 1 for more details). Each eBOSS object has a completeness weight including corrections for fibre collisions, redshift failures [69] and photometric systematic effects [70]: We refer to section 5.6 of [71] for details about these weights. The cosmology results for the DR16 quasar sample were presented in [4,72] 7 , which reported a BAO distance measurement in the range 0.8 < z < 2.2.
Mock datasets: We make use of a set of 1000 mock catalogs to estimate the covariance matrix of the eBOSS DR16 quasar power spectrum. The mocks are based on the Extended Zel'dovich (EZ) approximate N-body simulation scheme [73]. These mocks rely on the Zel'dovich approximation to generate a density field, while including nonlinear and halo biasing effects through the use of free parameters. These free parameters are tuned to produce two-point and three-point clustering of a desired data set. The EZ mock catalogs account for the redshift evolution of the eBOSS quasars by constructing a light-cone out of 7 redshift shells, generated from periodic boxes of side length L = 5000h −1 Mpc at different redshifts. These mock catalogs also mimic the fibre collisions and redshift failures, so that each object has an associated w foc and w cp . The cosmology of the EZmocks is ΛCDM with Ω m = 0.307115, Ω b = 0.048206, h = 0.6777, σ 8 = 0.8255, and n s = 0.9611 (same cosmology as the BOSS DR12 MD-Patchy mocks discussed above). Figure 10 shows the window function monopole for the different samples where the NGC sample of eBOSS DR16 indicates the smallest fundamental mode of all samples discussed in this paper. The fact that the eBOSS quasar sample probes the largest scales currently accessible with galaxy redshift surveys, makes it a valuable tool to test primordial non-Gaussianity [32]. However, the high shot noise level does mean that the effective volume of the eBOSS DR16 QSO sample at most wavenumbers is below BOSS DR12.

Data analysis
Here we will apply the deconvolution formalism developed in section 4 to the three datasets introduced in the previous section.

Setup of power spectrum measurements
All power spectra discussed in this section are measured using the estimator of [74] and [75] including the odd multipoles as discussed in appendix E of [2]. For all power spectrum measurements we use a cubic grid with L = 2000, 3500 h −1 Mpc and N = 600, 700 for 6dFGS and BOSS, respectively. For the eBOSS QSO sample we use L = 5400, 6600 h −1 Mpc with N = 1040, 1270 for the SGC and NGC, respectively. This setup ensures that for all three surveys the Nyquist frequency is > 0.6 h Mpc −1 . The galaxies are assigned to the grid using the triangular shape cloud procedure and we correct for the associated pixel window function [76]. To further reduce aliasing effects we included the interlacing procedure of [77]. As defined in Section 2, we use a different normalization for the power spectra than other recent papers, which brings the convolved power closer to the true/deconvolved power. For comparison with past work, Table 1 gives the ratio, of our definition of the normalisation to the one that we would compute using eq. (13) of [15]. Since deconvolution requires a square matrix in eq. (4.2), we also use ∆k th = 0.01 h Mpc −1 . In figure 11 we can see the increase in the variance caused by deconvolution (comparison of the red and black shaded regions). At the same time deconvolution reduces the correlation between bandpowers. The odd multipoles in figure 12 also include a best fitting model (green dashed lines) as well as the intrinsic dipole and octopole (blue dashed lines) as given by eq. (2.14) and eq. (2.15). These models are based on a fit to the even multipoles as discussed in section 3.6 of [2]. One can clearly see that deconvolving the dipole removes the window function contributions, which dominate the dipole on most scales and recovers the intrinsic dipole expected due to wide-angle effects. We included the corresponding results for 6dFGS and eBOSS in appendix E. Figure 14 in appendix E shows the convolved and deconvolved power spectrum multipoles of 6dFGS. From these plots we can see that the monopole power spectrum of the mock catalogs does not perfectly match the data power spectrum amplitude, while higher order multipoles agree well [78]. Even though 6dFGS is the lowest redshift sample, the wide-angle effects seem far less important compared to BOSS DR12, and the dipole is (by eye) consistent with zero. This agrees with [45] where wide-angle effects are discussed in appendix C. Window function effects also seem to be much smaller than the statistical noise, which might be caused by the very compact geometry of the 6dFGS survey (see [37] for details about the angular and redshift distribution of 6dFGS).

Deconvolution and BAO
The convolution of the measured power spectrum with the survey window function could smear out a signal which exists at the same scale as the fundamental mode (or the k o band width, whichever is larger). The BAO signal has a wavelength of ∼ 0.06 h Mpc −1 , which is much larger than the fundamental mode of the surveys discussed in this paper (see table 1). Nevertheless, deconvolving the power spectrum does increase the BAO signature as clearly visible in figure 13. This figure shows the BAO signal in the low and high redshift bins of the BOSS DR12 MD-Patchy mock catalogs together with the best fitting models. The data points represent the mean and variance of 1000 post-reconstruction MD-Patchy mock catalogs.
It is important to note that even though the BAO signal appears to be enhanced postdeconvolution, our deconvolution procedure is based on eq. (4.1) and hence the convolved and  deconvolved analysis should lead to the same likelihood and the same model parameters, if the same cuts and theory binning are used. To demonstrate this point we perform an isotropic BAO analysis based on the mean of 1000 BOSS DR12 MD-Patchy mock power spectrum monopoles (post-reconstruction). For the pre-deconvolution fit we build the monopole model following the current standard analysis pipeline described in [30]. We also provide higher order multipoles for the quadrupole and hexadecapole based on the simple Kaiser model. We use the window function and wide-angle matrices (W and M) and limit the likelihood evaluation to the monopole only (see the second example in appendix A). Our fitting range is 0.01 < k < 0.3 h Mpc −1 and we are jointly fitting the NGC and SGC. For the postdeconvolution fit we only fit the monopole without any window function or wide-angle matrix, again using the model of [30]. Our fitting procedure is using the Python-based MCMC sampler zeus [79,80] 9 . In the low redshift bin of BOSS DR12 we find α = 0.999 ± 0.013 before deconvolution and an identical value after deconvolution. Figure 13 (left) compares the best fitting model and the mean of the MD-Patchy mocks (the plot shows the weighted average of the NGC and SGC). The equivalent values for the high redshift bin shown in figure 13 (right) are α = 1.003 ± 0.012 before deconvolution and α = 1.002 ± 0.012 after deconvolution. Any observed differences are consistent with noise in the MCMC chain, the slightly different cuts implied by using the monopole within 0.01 < k < 0.3 h Mpc −1 in convolved vs. deconvolved space, and coarsening of the theory side of the window function used in the deconvolution. The window matrix used for deconvolution has to be a square matrix with relatively large k-bins (∆k = 0.01 h Mpc −1 in our case). Such large k-bins cannot capture small-scale features in the theory power spectrum (see figure 4). While these effects should not introduce any issues within a smooth ΛCDM power spectrum (as demonstrated here for the BAO case), it could become relevant for non-ΛCDM models especially if small scale features are present.

Conclusion
When analysing galaxy redshift surveys in Fourier space, one needs to account for the survey window function as well as wide-angle effects. In this paper we leverage recent new developments dealing with the survey window function and wide-angle effects, to lay out a simple power spectrum analysis framework based on matrix multiplications. The main results of this paper are: (1) We derive a matrix to account for wide-angle effects in the power spectrum multipoles.
We use a new analytic approach rather than the commonly used Hankel transforms.
(2) We expand the window function matrix approach presented in [27] by including wideangle effects.
(3) We use this matrix-based analysis framework for the power spectrum multipoles to demonstrate two possible analysis pipelines, one using the standard path of convolving the model vector and one based on the deconvolution of the data vector.
(4) We apply the deconvolution procedure to a set of existing galaxy redshift surveys, namely 6dFGS DR3, BOSS DR12 and eBOSS DR16. Using a BAO analysis we demonstrate that our deconvolution analysis framework leads to the same likelihood as the standard analysis.
(5) We provide the power spectrum multipoles as well as the window function matrices, wide-angle matrices and covariance matrices for 6dFGS DR3, BOSS DR12 and eBOSS DR16. In the appendix we also provide Python-based examples and a general user guide for a clustering analysis. These easy to use components hopefully simplify the analysis of these datasets and make them more accessible for the wider cosmology community.
The deconvolution framework outlined in this paper does not suffer from limitations inherent to other methods presented in the literature, such as the assumption of a global plane parallel approximation. Nevertheless, the inversion of the window function does require a square window matrix, which enforces large ∆k th bins, a limitation which is not present when convolving the model vector.
Our analysis focuses on the key science targets of galaxy redshift surveys such as RSD and BAO. Other observables such as primordial non-Gaussianity through the scale-dependent bias in the power spectrum, naturally requires to focus on the largest scales of the survey. The products provided with this paper are sub-optimal for such an observable. However, the formalism presented in this paper can easily be adapted to suit such an observable.
A similar matrix-based analysis approach could also be developed for higher order statistics, like the bispectrum. The extension of our analysis framework to higher order statistics will be addressed in future work.
regarding the window function normalisation. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 853291). FB is a Royal Society University Research Fellow. PM was supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Contract no. DE-AC02-05CH11231.

A User guide
Together with this publication we provide: All products listed above are available at https://fbeutler.github.io/hub/deconv_paper. html and Python-based modules to read these quantities are available at https://github. com/fbeutler/pk_tools. We note that due to the Nyquist frequency in the measured power spectra, any analysis should be limited to k max < 0.3 h Mpc −1 . The purpose for the addition of the k-range 0.3 < k < 0.4 h Mpc −1 is mainly to allow a sensible window function contribution from outside the fitting range. We generally advice to exclude the first k-bin (k min > 0.01 h Mpc −1 ), which could suffer from large scale systematics not studied in detail in this analysis.
Here is an example of a likelihood analysis using Python: Note that you can speed up these calculations by multiplying the matrices W and M beforehand, which reduces the required matrix multiplications in the likelihood evaluation to one.
How can I limit the k-range? The power spectrum model should always have a k-range of 0 < k th < 0.4 h Mpc −1 in bins of ∆k th = 0.001 h Mpc −1 . If your model cannot predict the power spectrum up to k max = 0.4 h Mpc −1 , you should provide some "sensible" estimate, so that the window function contributions from those scales can be included. The likelihood analysis should be limit to k max < 0.3 h Mpc −1 , since the power spectra outside this k-range suffer from aliasing given that k Ny ≈ 0.6 h Mpc −1 for all measured power spectra. One can specify the k-range as  This implementation is about a factor of two times faster than the brute-force implementation shown in the first code example. Of course that is only significant if the window function convolution is dominating the likelihood evaluation. In many cases most of the time will be spend in your_favourite_model().
How can I account for uncertainties in the covariance matrix? When deriving the covariance matrix from a finite set of mock realisations, the resulting likelihood is no longer Gaussian, but follows a t-distribution. If assuming a Gaussian likelihood the parameter inference will be biased and this bias depends on the ratio of bins in the data vector and the number of mock realisations [85]. We can account for this by scaling the likelihood as where N m is the number of mock realisations (given in table 1) and N d is the size of the data vector. Alternatively one can directly account for the non-Gaussian likelihood as proposed in [86]. For all cases discussed in this paper this approach agrees very well with eq. (A.6). Assuming your MCMC sampler expects log L as a return value you could implement this equation as If the final parameter uncertainty is derived from the likelihood itself, you need to account for a bias caused by the mock based covariance estimate. We can do that with a re-scaling of the parameter errors [87,88] by the square root of , (we find that the approach of [86] does not take care of this factor).
How can I deconvolve a power spectrum measurement? To perform a deconvolution, all you have to do is to follow eq. (4.2). The difficulty here is that this equation only holds for square matrices W. The window functions using ∆k o = ∆k th = 0.01 h Mpc −1 are available.  where the wide-angle correction terms at n = 0 and n = 2 have the same shape.

C Derivation of the 2D window function
In this section we derive the equation for the 2D window function W (n) (k, k ). We first derive this equation for n = 0 showing consistency with [1], where this equation first appeared. We than include the wide-angle correction terms following [2], which leads to our eq. (2.5).

C.1 Excluding wide-angle terms and consistency with [1]
Here we show the relation between eq. (33) of [1] and eq. (2.2). The convolution of the power spectrum can be written as (using a LOS ofd =ŝ 1 and following eq. B.1 of [1]) e ik·s e −ik ·s L (k ·ŝ 1 ) .

Now using
Nran i,j,i=j and the definition (see eq. 2.21 of [2]) as well as

C.2 Including wide-angle terms and consistency with [2]
The convolution of the power spectrum multipoles including the wide-angle correction terms has first been derived in [2] and is given by 11 Now using ξ (n) The factors C D Analytic calculation of the odd power spectrum multipoles Here we derive eq. (2.14) and eq. (2.15) used in section 2.1. The dipole power spectrum is given by Following eq. (F.1), (F.6) and (F.10) of [12] we can write where δ D is the Dirac delta function. Using to replace the integral with a derivative, results in where the second line assumes linear theory. This final equation is numerically much easier to evaluate compared to the equation we started with. Equivalently we can derive the octopole (D.10)   Figure 14: Comparison between the power spectrum multipoles of 6dFGS DR3 measured in the mock catalogs (gray and red shaded area) and in the data (data points). The results before deconvolution are shown as the red shaded area and solid red line (mocks) and the green data points. The deconvolved results are shown as the gray shaded area and solid black line (mocks) and the blue data points. The residuals in the lower panel show that wide-angle and window function effects are sub-dominant in 6dFGS.  Figure 15: Comparison between the power spectrum multipoles of BOSS DR12 measured in the mock catalogs (gray and red shaded area) and in the data (data points). The results before deconvolution are shown as the red shaded area and solid red line (mocks) and the green data points. The deconvolved results are shown as the gray shaded area and solid black line (mocks) and the blue data points. The equivalent BOSS DR12 NGC results for the low redshift bin are included in the main text (see figure 11 and 12). Pdeconv(k) Figure 16: Comparison between the power spectrum multipoles of the eBOSS DR16 QSO sample measured in the mock catalogs (gray and red shaded area) and in the data (data points). The results before deconvolution are shown as the red shaded area and solid red line (mocks) and the green data points. The deconvolved results are shown as the gray shaded area and solid black line (mocks) and the blue data points. The increasing in the dipole power spectrum at high k seems to be an aliasing effect, even though the Nyquist frequency is twice as high as the scale range plotted here (k Ny > 0.6 h Mpc −1 ).