© T. Dumerchat and J. E. Bautista 2022

1. Introduction

The study of the accelerated nature of the expansion of the Universe has seen significant progress in the past decade thanks to measurements of baryon acoustic oscillations (BAOs) in the three-dimensional distribution of galaxies. Similarly to type-Ia supernovae (see Brout et al. 2022 and references therein for the latest results), BAO measurements have shown that the expansion is accelerating, imposing the need of a dark energy component in the energy budget of the Universe. The next decade will see a large flow of data coming from galaxy surveys whose goal is the precise measurement of BAO over a large span of the cosmic history.

The latest BAO measurements from spectroscopic surveys were produced by the cosmological component of the fourth generation of the Sloan Digital Sky Survey (SDSS-IV, Blanton et al. 2017), named the extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016). The eBOSS project produced measurements of BAO using luminous red galaxies at an effective redshift z_eff = 0.7 (Bautista et al. 2021; Gil-Marín et al. 2020), using emission line galaxies at z_eff = 0.85 (Raichoor et al. 2021; Tamone et al. 2020; de Mattia et al. 2021), using quasars as tracers of the matter field at z_eff = 1.48 (Hou et al. 2021; Neveux et al. 2020), and quasars with visible Lyman-α forests at z_eff = 2.33 (du Mas des Bourboux et al. 2020). Two additional lower redshift measurements were made at z_eff = 0.38 and 0.51 from the third generation of SDSS, BOSS (Eisenstein et al. 2011; Dawson et al. 2013; Alam et al. 2017). At z_eff = 0.15, we have measurements from the SDSS-II Main Galaxy Sample (Ross et al. 2015a; Howlett et al. 2015). Measurements from other surveys include the 6dFGS (Beutler et al. 2011) at z_eff = 0.10 and the WiggleZ Dark Energy Survey (Kazin et al. 2014) with three redshift bins spanning 0.2 <  z <  1.0. The latest BAO measurements from photometric surveys were produced by the Dark Energy Survey (DES, DES Collaboration 2016). DES produced angular BAO measurements in five redshift intervals between 0.6 and 1.1 using three years of its data (DES Collaboration 2021).

Baryon acoustic oscillation measurements have traditionally been performed in configuration space and/or in Fourier space. Configuration space analysis is based on estimates of correlation function ξ as a function of separation r, commonly done with pair-counting techniques. Fourier space analyses assign galaxies into a regular grid so as to Fourier transform it in order to compute the power spectrum multipoles P_ℓ as a function of wavevector k. Both types of analyses perform a statistical measurement starting from the exact same dataset, defined as a list of angular positions, redshifts, and weights, as well as some definition of the window function of the survey which is often a purely Poisson set of points following the same angular and redshift distribution as the real data. In principle, both types of analyses should yield the same cosmological constraints since they have a common starting point. In practice, choosing a limited range of scales used when fitting models slightly breaks this perfect degeneracy. Noise properties or systematic effects on the estimated (and binned) statistics also differ between ξ(r) and P(k). The models of clustering used to fit two-point functions also commonly differ, that is to say when performing a full shape analysis the clustering models in configuration space can differ from just a Fourier transform of those in Fourier space. All these differences slightly reduce the correlations to levels below 100%. As an example, Bautista et al. (2021) show using the best-fit values for the dilation parameters on an ensemble of mock catalogs that BAO results in configuration and Fourier spaces are roughly 90% correlated, for their particular choices of scales and models. Recent analysis of BAO and redshift-space distortions (RSD) usually combine results from configuration and Fourier space analyses into a single consensus result.

Given that the correlations between cosmological results are not exactly 100% between configuration and Fourier space analyses, there is a slight statistical gain in producing a combined result. Systematic errors might as well be reduced by this combination, since they may manifest themselves differently in both spaces. In collaboration working groups, often two or more teams produce their analyses in a single space, making it difficult to decide which results to quote as final. Therefore, there are several advantages to produce joint results between Fourier and configuration spaces.

When calculating a consensus result, it is important to correctly take into account the strong correlations between measurements in configuration and Fourier space. Sánchez et al. (2017) described a method to compute a consensus result, which was used for the first time in SDSS-III BOSS (Alam et al. 2017) and extended to the latest measurements from SDSS-IV eBOSS survey. This method was used to produce consensus results on BAO and RSD analyses of a given survey, yielding constraints on D_H/r_d, D_M/r_d and fσ₈, where D_H is the Hubble distance c/H, D_M is the comoving angular diameter distance, r_d is the comoving sound horizon at drag epoch (the BAO scale), f is the growth-rate of structures and σ₈ is the normalization of the smoothed linear matter power spectrum. However, the method from Sánchez et al. (2017) assumes that the individual likelihoods on BAO parameters are Gaussian, which is not necessarily true, particularly in a regime of low signal-to-noise ratio.

In this work, we developed a framework to perform simultaneous BAO analyses in both configuration and Fourier spaces. The advantages of our framework over past work are (1) individual likelihoods on BAO parameters from each space are not assumed to be Gaussian, (2) the resulting posterior distribution is not necessarily Gaussian, resulting in more reliable uncertainties, (3) it is simpler, and (4) yields smaller systematic biases. Our method only relies on a sufficiently good estimate of the full covariance matrix, particularly the cross-covariance between two-point functions in Fourier and configuration spaces. We validated our method on realistic mock catalogs and performed comparison with previous work.

This paper is organized as follows. In Sect. 2, we describe the dataset used to validate our methods. In Sect. 3, we introduce the BAO modeling, which is mostly the same used in Bautista et al. (2021), Gil-Marín et al. (2020), as well as the methods to produce consensus results. Section 4 presents several statistical and systematical tests of both methods and finally the application on real data on Sect. 5.

2. Dataset

To validate our methodology, we used 1000 mock catalogs reproducing the Luminous Red Galaxy (LRG) sample from the eBOSS survey, though our methods are not survey specific. This dataset is the same employed in the cosmological analyses of Bautista et al. (2021), Gil-Marín et al. (2020), and we refer the reader to these references for further details.

2.1. The eBOSS survey

The extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016) was a 5-year observing program using multiobject fiber-fed spectrographs (Smee et al. 2013) mounted on the focal plane of the 2.5 m Sloan Foundation Telescope (Gunn et al. 2006) at the Apache Point Observatory. During eBOSS, 174 816 LRG redshifts were obtained over 4242 deg² of both northern and southern skies, in the redshift interval 0.6 <  z <  1.0. These were combined with BOSS galaxies in the same redshift range covering 9493 deg² of sky, resulting in a total of 377 458 LRG redshifts (Ross et al. 2020). The survey geometry was defined using a set of randomly distributed points, taking into account masked areas. The final data sample has correction weights accounting for the photometric and spectroscopic incompleteness, as well as some spurious correlations (e.g., with stellar density, Galactic extinction). The final catalog is the basis for mock catalog production.

2.2. Mock catalogs

In this work we use a sample of 1000 mock realizations of the eBOSS LRG sample to validate our methodology and to perform statistical tests. A set of 1000 realisations of the eBOSS LRG survey were produced using the EZMOCK method (Chuang et al. 2015), which employs the Zel’dovich approximation to compute the density field at a given redshift. This method is faster than n-body simulations and has been calibrated to reproduce the two- and three-point statistics of the given galaxy sample, to a good approximation and up to mildly nonlinear scales. The angular and redshift distributions of the eBOSS LRG sample were reproduced, including both photo and spectro-incompleteness effects as well as systematic effects introducing spurious correlations. These mocks aim to include all known features in real data. A detailed description of the EZMOCKS can be found in Zhao et al. (2021). These mocks were used in this work to estimate covariance matrices and to perform statistical studies of our method. As we will describe in Sect. 3.4, mocks are also required by the Gaussian approximation method for combining configuration and Fourier space results into a single consensus result.

3. Methodology

3.1. Measuring the clustering

We used correlation function and power spectrum multipoles estimated for each of the 1000 EZMOCK realisations. The correlation function is estimated using the Landy & Szalay (1993) estimator, while power spectra were calculated using the Yamamoto et al. (2006) estimator, as implemented by Bianchi et al. (2015), Scoccimarro (2015).

We attempt to remove the impact of non-linear evolution of galaxies and a fraction of redshift-space distortions by applying the reconstruction technique of Burden et al. (2015). Using the density field itself, we estimate the displacements using the Zeldovich approximation to move the galaxies “back in time”. This technique increases the precision of the measurement by sharpening the BAO feature. According to Carter et al. (2020) the BAO results are not sensitive to small variations in the cosmology used to perform the reconstruction. In Appendix B, we performed the analysis for the prereconstructed sample, where we derive the same conclusions as for the postreconstruction cases (though at lower signal-to-noise ratios).

3.2. BAO modelling

The models used in this work are the same employed in Bautista et al. (2021), Gil-Marín et al. (2020), which are themselves based on previous work (Alam et al. 2017; Gil-Marín et al. 2016, 2018; Ross et al. 2017; Bautista et al. 2018). We briefly summarize these models here. The code that produces the model and perform the fitting to the data is publicly available¹.

The aim is to model both the power spectrum P_ℓ(k) and correlation function multipoles ξ_ℓ(r) as a function of wave-number k and separations r, respectively, where ℓ is the order of the multipole expansion. We focus on scales relevant for the measurement of the baryon acoustic oscillation (BAO) feature, typically 0 <  k <  0.3 h⁻¹ Mpc for power-spectrum and 30 <  r <  180 h⁻¹ Mpc for correlation function.

The starting point is a linear-theory-based model for the redshift-space anisotropic galaxy power-spectrum P(k, μ),

$$\begin{array}{cc}\hfill P(k,{\mu }_k)=& \frac{b^2{[1+\beta (1-S(k)){\mu }_k^2]}^2}{(1+k^2{\mu }_k^2{\mathrm{\Sigma }}_{\mathrm{s}}^2/2)}\hfill \\ \hfill & \times [P_{\mathrm{no}\mathrm{peak}}(k)+P_{\mathrm{peak}}(k)e^{-k^2{\mathrm{\Sigma }}_{\mathrm{nl}}^2({\mu }_k)/2}],\hfill \end{array}$$(1)

where b is the linear density bias of the galaxy population, β is the redshift-space distortions parameter defined as the ratio between the growth-rate f and b, k is the modulus of the wave-vector k and μ_k is the cosine of the angle between the wave-vector and the line of sight. The broadening of the BAO peak caused by non-linear clustering is reproduced by applying an anisotropic Gaussian smoothing, that is by multiplying the “peak-only” power spectrum P_peak (see below) by a Gaussian function with dispersion given by ${\mathrm{\Sigma }}_{nl}^2({\mu }_k)={\mathrm{\Sigma }}_{\parallel }^2{\mu }_k^2+{\mathrm{\Sigma }}_{\perp }^2(1-{\mu }_k^2)$. The non-linear random motions on small scales are modeled by a Lorentzian distribution parametrized by Σ_s.

Following Seo et al. (2016), the term $S(k)=e^{-k^2{\mathrm{\Sigma }}_r^2/2}$ for the postreconstruction model and S(k)=0 for the prereconstruction, where Σ_r = 15 h⁻¹ Mpc is the smoothing parameter used when reconstructing the galaxy catalogs.

We follow the procedure from Kirkby et al. (2013) to decompose the BAO peak component P_peak from the full linear power-spectrum P_lin. We start by computing the correlation function by Fourier transforming P_lin, then we replace the correlations over the peak region by a polynomial function fitted using information outside the peak region (50 <  r <  80 and 160 <  r <  190 h⁻¹ Mpc). The resulting correlation function is then Fourier transformed back to get P_no peak. The linear power spectrum P_lin is obtained from CAMB² (Lewis et al. 2000) with cosmological parameters of our fiducial cosmology defined in Table 1.

The multipoles of the power-spectrum P_ℓ(k) are obtained by integrating over μ_k weighted by the Lengendre polynomials L_ℓ(μ_k):

$$\begin{array}{c}\hfill P_{\ell }(k)=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}P(k,\mu )L_{\ell }(\mu )\mathrm{d}\mu .\end{array}$$(2)

The correlation function multiples ξ_ℓ(r) are obtained by Hankel transforming the P_ℓ(k):

$$\begin{array}{c}\hfill {\xi }_{\ell }(r)=\frac{i^{\ell }}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{k}^2{j}_{\ell }(kr)P_{\ell }(k)\mathrm{d}k,\end{array}$$(3)

where j_ℓ are the spherical Bessel functions. These transforms are computed using HANKL³ that implements the FFTLog algorithm by Hamilton (2000).

The BAO peak position is parametrised via two dilation parameters, one scaling separations across the line of sight, α_⊥, and one scaling separations along the line of sight, α_∥. The observed k is related to the true k′ by $k_{\parallel }=k_{\parallel }^{\mathrm{\prime }}/{\alpha }_{\parallel }$ and k_⊥ = k′/α_⊥. Therefore, the observed power spectrum relates to the true power spectrum as

$$\begin{array}{c}\hfill P(k_{\Vert },k_{\perp })=\frac{1}{{{\alpha }_{\perp }}^2{\alpha }_{\Vert }}{P}^{\prime }(\frac{k_{\Vert }}{{\alpha }_{\Vert }},\frac{k_{\perp }}{{\alpha }_{\perp }})·\end{array}$$(4)

We note that this volume scaling does not apply for the correlation function:

$$\begin{array}{c}\hfill \xi (r_{\Vert },r_{\perp })={\xi }^{\prime }({\alpha }_{\Vert }{r}_{\Vert },{\alpha }_{\perp }{r}_{\perp }).\end{array}$$(5)

These BAO dilation parameters are related, respectively, to the comoving angular diameter distance, D_M = (1 + z)D_A(z), and to the Hubble distance, D_H = c/H(z), by

$$\begin{array}{cc}& {\alpha }_{\perp }=\frac{D_{\mathrm{M}}(z_{\mathrm{eff}})/r_{\mathrm{d}}}{D_{\mathrm{M}}^{\mathrm{fid}}(z_{\mathrm{eff}})/r_{\mathrm{d}}^{\mathrm{fid}}},\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\alpha }_{\Vert }=\frac{D_{\mathrm{H}}(z_{\mathrm{eff}})/r_{\mathrm{d}}}{D_{\mathrm{H}}^{\mathrm{fid}}(z_{\mathrm{eff}})/r_{\mathrm{d}}^{\mathrm{fid}}},\hfill \end{array}$$(7)

where r_d is the comoving sound horizon at drag epoch, and z_eff is the effective redshift of the survey. The r_d rescaling is here to account for the choice of template cosmology used to compute the fixed linear power spectrum. For simplicity, the template cosmology is chosen to match the fiducial cosmology used to estimate distances from redshifts.

We apply the scaling factors exclusively to the peak component of the power spectrum P_peak(k), effectively removing any dependency of these parameters on the smooth part P_{no peak}(k) (Kirkby et al. 2013).

In order to properly marginalize over any mis-modelling of the smooth part of P_ℓ(k) and ξ_ℓ(r), we add to our model a linear combination of smooth functions of scale, with free amplitudes to be marginalised over. These smooth functions can also account for potential unknown systematic correlations that contaminating our measurements. Furthermore, since there are no accurate analytical models for correlations postreconstruction (the S(k) term in Eq. (1) is generally not sufficient), these smooth functions can also account for this mis-modelling. As these smooth functions are highly correlated with the finger of god Lorentzian, we fix the parameter Σ_s to zero. Our final template can be written as:

$$\begin{array}{cc}& P_{\ell }^t(k)=P_{\ell }({\alpha }_{\perp },{\alpha }_{\Vert },k)+{\displaystyle \sum _{i=i_{\min }}^{i_{\max }}}{a}_{\ell ,i}^P{k}^i,\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\xi }_{\ell }^t(r)={\xi }_{\ell }({\alpha }_{\perp },{\alpha }_{\Vert },r)+{\displaystyle \sum _{i=i_{\min }}^{i_{\max }}}{a}_{\ell ,i}^{\xi }{r}^i,\hfill \end{array}$$(9)

where $a_{\ell ,i}^P$ and $a_{\ell ,i}^{\xi }$ are the amplitudes for each power i of scale, k or r, and multipole order ℓ. In Gil-Marín et al. (2020), the BAO analysis in Fourier space used (i_min, i_max)=(−1, 1), while in Bautista et al. (2021) the configuration space analysis used (i_min, i_max)=(−2, 0). For both, this corresponds to three free parameters per multipole. In this work, our baseline choice is (i_min, i_max)=(−2, 1) when performing joint fits unless stated otherwise.

Our baseline BAO analysis uses the monopoles P₀, ξ₀ and quadrupoles P₂, ξ₂ of the power spectrum and correlation functions (see Appendix A for results using the hexadecapole). We fix β = 0.35 and fitting b with a flat prior between b = 1 and 4. For all fits, the broadband parameters are free, while both dilation parameters are allowed to vary between 0.5 and 1.5.

3.3. Parameter inference

The cosmological parameter inference is performed by means of the likelihood analysis of the data. The likelihood ℒ is defined such that

$$\begin{array}{c}\hfill -2ln\mathcal{L}(\theta )={\displaystyle \sum _{i,j}^{N_p}}{\mathrm{\Delta }}_i(\theta ){\widehat{\mathrm{\Psi }}}_{\mathit{ij}}{\mathrm{\Delta }}_j(\theta )+C={\chi }^2+C,\end{array}$$(10)

where θ is the vector of parameters, Δ is a vector containing residuals between observed multipoles and their model, N_p is the total number of elements in Δ. An estimate of the precision matrix $\widehat{\mathrm{\Psi }}=(1-D){\widehat{C}}^{-1}$ is obtained from the unbiased estimate of the covariance from 1000 realisation of EZmocks: ${\widehat{C}}_{\mathit{ij}}={(N_{\mathrm{mock}}-1)}^{-1}{\displaystyle \sum _k^{N_{\mathrm{mock}}}}(X_i^k-⟨X_i⟩)(X_j^k-⟨X_j⟩)$ where $X_i^k$ is the measured P_ℓ(k_i) or ξ_ℓ(r_i) for the kth realisation. The factor D = (N_p + 1)/(N_mocks − 1) accounts for the skewed nature of the Wishart distribution (Hartlap et al. 2007).

The best-fit BAO parameters (α_⊥, α_∥) are determined by minimizing χ² of Eq. (10) using a quasi-Newton minimum finder algorithm IMINUIT⁴ which marginalizes over nuisance parameters while sampling for the parameters of interest. The uncertainties in α_∥ and α_⊥ are estimated with the MINOS function provided by IMINUIT. This functions computes the intervals where χ² increases by unity, which corresponds to a 68% confidence interval. Gaussianity is not assumed in this calculation and uncertainties can be asymmetric with respect to the best-fit value. The two-dimensional confidence contours in (α_⊥, α_∥), are estimated using the CONTOUR function from IMINUIT. Similarly to MINOS, this function scans χ² values in two dimensions, looking for the contours yielding a Δχ² = 2.3 or 5.9 for 68 and 95% confidence levels, respectively.

3.4. Consensus via Gaussian approximation

In previous work, configuration and Fourier space results were combined into a single consensus result using the method presented in Sánchez et al. (2017).

The idea of the method is to generate a consensus result from M different measurement vectors x_m, each containing p elements, and their covariance matrices C_mm, each of size p × p. For example, we want to combine M = 2 measurements of the vector x_m = [α_⊥, α_∥] containing p = 2 parameters from Fourier (m = 1) and configuration space (m = 2), each with its own 2 × 2 error matrix C_mm derived from their posteriors. The consensus is a single vector x_c with covariance C_c, for which the expressions assume that the χ² between individual measurements is the same as the one from the consensus result. The expression for the combined covariance matrix is

$$\begin{array}{c}\hfill {\mathbf{C}}_c\equiv {\left({\displaystyle \sum _{m=1}^M\sum _{n=1}^M}{\mathbf{C}}_{\mathit{mn}}^{-1}\right)}^{-1},\end{array}$$(11)

and the combined data vector is

$$\begin{array}{c}\hfill {\mathit{x}}_c={\mathbf{C}}_c{\displaystyle \sum _{m=1}^M}\left({\displaystyle \sum _{n=1}^M}{\mathbf{C}}_{\mathit{nm}}^{-1}\right){\mathit{x}}_m,\end{array}$$(12)

where C_mn is a p × p block from the full covariance matrix between all parameters and methods 𝒞, containing pM × pM elements, defined as

$$\begin{array}{c}\hfill \mathcal{C}=\left(\begin{array}{cccc}{\mathbf{C}}_{11}& {\mathbf{C}}_{12}& \cdots & {\mathbf{C}}_{1M}\\ {\mathbf{C}}_{21}& {\mathbf{C}}_{22}& \cdots & {\mathbf{C}}_{2M}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{C}}_{M1}& {\mathbf{C}}_{M2}& \cdots & {\mathbf{C}}_{\mathit{MM}}\end{array}\right).\end{array}$$(13)

The diagonal blocks C_mm come from each measurement method M, therefore, assuming Gaussian likelihoods for the p parameters. The off-diagonal blocks C_mn with m ≠ n cannot be in principle estimated from the data itself. These off-diagonal blocks are commonly built from mock catalogs. Using a set of many realisations, one can build 𝒞 from all the realisations of x_m for each method. We obtain the correlation coefficients ${\rho }_{p_1,p_2,m,n}^{\mathrm{mocks}}$, that is the covariance 𝒞 normalized by its diagonal elements, between parameters p₁, p₂ and methods m, n. We scale these coefficients by the diagonal errors from the data, to obtain the final matrix 𝒞 for the data.

It is worth emphasizing that the matrix 𝒞 for our data, and more specifically its off-diagonal blocks, depend on that particular realisation of the data in principle. However, the ones derived from mock measurements are ensemble averages. We account for this fact by scaling the correlations coefficients from the mocks in order to match the maximum correlation coefficient that would be possible with the data (Ross et al. 2015b). For the same parameter p₁ measured by two different methods m and n, we assume that the maximum correlation between them is given by ρ_max = σ_p1, m/σ_p1, n, where σ_p is the error of parameter p. This number is computed for the data realisation ${\rho }_{\max }^{\mathrm{data}}$ and for the ensemble of mocks ${\rho }_{\max }^{\mathrm{mocks}}$. We can write the adjusted correlation coefficients for one particular realisation as

$$\begin{array}{c}\hfill {\rho }_{p_1,p_1,m,n}^{\mathrm{data}}={\rho }_{p_1,p_1,m,n}^{\mathrm{mocks}}\frac{{\rho }_{\max }^{\mathrm{data}}}{{\rho }_{\max }^{\mathrm{mocks}}}·\end{array}$$(14)

The equation above accounts for the diagonal terms of the off-diagonal block C_mn. For the off-diagonal terms, we use

$$\begin{array}{c}\hfill {\rho }_{p_1,p_2,m,n}^{\mathrm{data}}=\frac{1}{4}({\rho }_{p_1,p_1,m,n}^{\mathrm{data}}+{\rho }_{p_2,p_2,m,n}^{\mathrm{data}})({\rho }_{p_1,p_2,m,m}^{\mathrm{data}}+{\rho }_{p_1,p_2,n,n}^{\mathrm{data}}).\end{array}$$(15)

This adjustment of correlation coefficients is an addition to the method proposed by Sánchez et al. (2017) and was first implemented in Bautista et al. (2021) to obtain consensus results.

In summary, the method proposed by Sánchez et al. (2017) can provide consensus results from strongly correlated measurements of the same quantities. However, it assumes Gaussian likelihoods and yields Gaussian posteriors on the final parameters, which is not a good approximation when measurements are noisy. Also, it relies on mock measurements twice: once to derive the covariance matrix of two-point functions, from which we derive each parameter vector x_m; and once again to derive the full parameter-method covariance 𝒞, using Eqs. (14) and (15). We call this method Gaussian approximation (GA). In the following section, we study joint fits in Fourier and configuration space, such that mock catalogs are used only once to derive a consensus result.

3.5. Joint analysis

We studied a new alternative method to get a unique set of constraints on the parameters α_∥ and α_⊥, taking into account the full correlations between the configuration space (CS) and Fourier space (FS) measurements. We refer to this method as Joint Space (JS) analysis. We concatenate the measured ξ_ℓ multipoles with the P_ℓ multipoles such that we obtain a data vector [ξ₀,ξ₂,P₀,P₂] with (2×r_bins+2×k_bins) data points. We then fit the data using a unique set of parameters, except for the nuisance parameters defining the polynomial smooth functions (Eqs. (8) and (9)).

In order to perform the joint inference, one needs to estimate the full covariance matrix. We use the measured P_ℓ and ξ_ℓ on the 1000 EZmocks to compute the covariance matrix $\widehat{\mathit{C}}$. As in Bautista et al. (2021) and Gil-Marín et al. (2020), we apply the following scale cuts: r ∈ [50, 150] h⁻¹ Mpc and k ∈ [0.02, 0.3] h⁻¹ Mpc. Figure 1 shows the resulting correlation matrix, defined as $\mathit{R}={\mathit{\sigma }}^{-1}\widehat{\mathit{C}}{[{\mathit{\sigma }}^{-1}]}^T$, where σ is the vector composed by the square-root of the diagonal elements of $\widehat{\mathit{C}}$. The diagonal blocks show the well-known correlations between the monopole and quadrupole in the same space, while the off diagonal blocks reveal the correlation patterns between the two spaces. The correlations between the two spaces monopole-monopole and quadrupole-quadrupole present the same non-linear features. A given physical scale r is correlated (and anti-correlated) with several k modes, meaning that the statistical information contained in a given scale is spread out over several spectral modes. A physical modeling of these features is proposed in Appendix C.

Fig. 1. Correlation matrix of the correlation function ξℓ and power spectrum multipoles Pℓ estimated using the 1000 independent measurements of EZmocks.

We then need an estimate for the precision matrix $\widehat{\mathbf{\Psi }}$ in Eq. (10). Inverting the covariance mixes the different modes and scales across CS and FS. The top panel of Fig. 2 presents the resulting correlations in the normalized precision matrix (for convenience we call the coefficients of the precision matrix ‘correlations’). Green arrows indicate the regions detailed in the bottom panels, where we focus on the correlations between the CS ξ₂ and the FS P₀ and P₂, respectively by ploting the amplitude of the precision matrix at two fixed scaled r_min = 52.5 h⁻¹ Mpc and r_max = 147.5 h⁻¹ Mpc. The error bars of the precision matrix are computed according to Eq. (29) of Taylor et al. (2013). We see that at both scales, ξ₂ is weakly correlated with P₀ while being strongly correlated with P₂. This result is in agreement with the correlations between ξ₀ and ξ₂. Furthermore, for large r, while the number of correlated k modes increases, the amplitude of the correlations decreases distinctly for ξ₂ − P₂, and increases for ξ₂ − P₀. Note that the correlations between the ξ₀ and P_ℓ are not shown here but behave in a similar fashion.

Fig. 2. Precision matrix of the correlation function ξℓ and power spectrum multipoles Pℓ. Top panel: normalized estimated precision matrix $R_{\mathit{ij}}^{\xi P}={\mathrm{\Psi }}_{\mathit{ij}}/\sqrt{{\mathrm{\Psi }}_{\mathit{ii}}{\mathrm{\Psi }}_{\mathit{jj}}}$, obtained from 1000 EZmock realisations. Bottom panel: slices of $R_{\mathit{ij}}^{\xi P}$ showing the cross-correlation between ξℓ and Pℓ for particular values of scales as shown by the green arrows in the left panel. Errors are estimated from Eq. (29) of Taylor et al. (2013).

While capturing additional information about the scale dependence of the joint space correlations, our new methodology has some drawbacks. Indeed, because the covariance matrix $\widehat{C}$ we use to infer our set of parameters is an estimate of the real precision matrix drawn from a Wishart distribution (finite sample of mocks), each element is affected by its own uncertainty that should be correctly propagated to the uncertainty of the estimated parameters. It has been shown in Dodelson & Schneider (2013) that this additional “noise” is directly proportional to the parameter covariance. We therefore need to apply correction factors to the obtained C_θ. These factors are given in Percival et al. (2014) as:

$$\begin{array}{cc}& m_1=\frac{1+B(N_{\mathrm{bins}}-N_{\mathrm{par}})}{1+A+B(N_{\mathrm{par}}+1)},\hfill \\ & m_2=\frac{m_1}{1-D},\hfill \end{array}$$(16)

with D the Hartlap factor defined in Sect. 3.3 and

$$\begin{array}{cc}& A=\frac{2}{(N_{\mathrm{mock}}-N_{\mathrm{bins}}-1)(N_{\mathrm{mock}}-N_{\mathrm{bins}}-4)},\hfill \\ & B=\frac{(N_{\mathrm{mock}}-N_{\mathrm{bins}}-2)}{(N_{\mathrm{mock}}-N_{\mathrm{bins}}-1)(N_{\mathrm{mock}}-N_{\mathrm{bins}}-4)}·\hfill \end{array}$$(17)

Note that this is also true for a classic FS or CS analysis but since the correction factors only scale with the number of parameters, data bins and mocks used to estimate the covariance, the enlargement on the constraints is smaller. The factor m₁ is to be directly applied to the estimated parameter covariance matrix for a given measurement, while the factor m₂ scales the standard deviation of a given parameter over a set of mocks. The values of the parameters m₁ and m₂ are given in Table 2. Since the enlargement of the parameter constraints ($\sqrt{m_1}$) expected in CS is about 1%, 1.7% in FS and 3% in JS, we expect slightly looser constraints for the joint analysis.

4. Results on mock catalogs

In this section, we use postreconstruction EZMOCKS to validate our parameter inference and error estimation (see Appendix B for a discussion on the analysis of the prereconstruction EZMOCKS). The aim is to compare results from configuration space (CS), Fourier space (FS), Joint space (JS) with the consensus results from the Gaussian approximation (GA).

4.1. Fits on average correlations

Using the inference methodology described in the Sect. 3, we fit the average P_ℓ and ξ_ℓ of the 1000 EZmocks in order to study potential biases and compare the GA and JS methods.

At first, we model ξ_ℓ and P_ℓ in the separation ranges r ∈ [50, 150] h⁻¹ Mpc and k ∈ [0.02, 0.3] h⁻¹ Mpc, by setting the broad band smooth polynomial functions from Eqs. (8) and (9) to the shape (i^min, i^max)=(−2, 1) and letting all fitting parameters free. Then, we assess the robustness of the best fit parameters with respect to variations in Σ_⊥, Σ_∥, k ranges, r ranges, number of broadband terms, and template cosmology for the linear power spectrum. Each time we vary one of those settings we keep the other ones fixed in all spaces (CS, FS and JS), facilitating the comparisons of results. Since we fit the mean of the mocks, we normalize the covariance matrix by the total number of mocks N_mocks = 1000. As the covariance is estimated with the same 1000 realisations, we do not expect it to be highly accurate, but sufficient for our purposes.

Figure 3 presents the impact of different analysis settings on the best fit values of α_∥ and α_⊥ for the GA (in black) and JS (in green) analysis. In every plot, a gray shaded area represents the 1 percent deviation from the expected value, which was the tolerance used in previous analysis. First, we tested different set of values for the broadening of the BAO peak, (Σ_∥, Σ_⊥), corresponding to the best values found in CS (7.90, 5.58), FS (7.23, 4.72) and JS (7.43, 5.21), all in units of h⁻¹ Mpc. Both GA and JS analyses are robust to changes in those parameters. Both GA and JS analyses are also robust for different ranges of scales used in the fit, presenting small and coherent shifts. This indicates that the information is correctly extracted only from the BAO feature. When varying the broadband terms, the JS method shows stable results except for α_∥ in one setting (−1, 1), where the broadband is not flexible enough to fit the residuals. When increasing the number of polynomial terms, the systematic shifts are notably smaller in JS than GA. The last two panels of Fig. 3 show the influence of the template cosmology used to compute the linear power spectrum (not the one to convert redshifts into distances). We chose to vary separately the two parameters Ω_b and Ω_cdm. Each time, we compute the new expected values for the α parameters by rescaling $r_{\mathrm{d}}^{\mathrm{fid}}$ in Eqs. (6) and (7) (the fiducial distances are unchanged as the cosmology used to compute distances does not vary). We find that the systematic shift behave in the same way for the two methods, and do not exceed 1% when varying Ω_cdm or Ω_b by 10% from the baseline cosmology. Note however that we performed the fits while letting free the parameters Σ_∥ and Σ_⊥, which helps reducing the shift amplitude between the different cosmologies.

Fig. 3. Impact of the choice of fitting scales, fiducial cosmology, broadening Σ’s parameters and polynomial broadband order in the recovered values of the parameters α∥ and α⊥. Each point is the best-fit from the average of the 1000 EZmocks. The JS results are in green and GA in black. The grey shaded areas correspond to a 1% error and red shared areas indicate the fiducial choices of our analysis.

Overall, the inferred value for the parameter α_⊥ is the more robust than α_∥ for both methods. Globally, GA and JS behave in the same way when varying the fitting configuration, proving that these are not caused by the new methodology. Best values and uncertainties are consistent between GA and JS. The systematic shifts for α_∥ are smaller in JS, while no significant difference is seen for α_⊥. The chosen baseline for the rest of the analysis is highlighted in red. The systematic errors may partially result from the covariance matrix that would require more realisations, but as they are much smaller than statistical errors budget for the real data (σ_α/α ∼ 2%), we can safely neglect these.

4.2. Fits on individual mocks

We now focus on the statistical properties of best-fit values and their uncertainties, based on fits of individual realisations of EZmocks. When fitting individual mocks, we fix (Σ_∥,Σ_⊥) to the best-fit values of the Joint analysis in the stack of 1000 measurements (see previous section), while letting all other parameters free. Table 3 summarises the parameters and flat priors used in these fits. The constraints on (α_∥,α_⊥) are then obtained by marginalizing over b and all other nuisance parameters.

For each method (FS, CS and JS) we remove results from nonconverging likelihoods (unsuccessful contour estimation) and extreme best-fit values at 5σ level (σ is defined as half of the range covered by the central 68% values). The best-fit values for which the estimated uncertainties touch the prior boundaries are removed as well. The remaining number of realisations referred as N_good is given for each method in Table 4 and is consistent with the results of Gil-Marín et al. (2020) and Bautista et al. (2021).

To obtain GA results, we extract for each mock posterior profile the 1 sigma contour of the parameter space. We fit the likelihood contours corresponding to 68% C.L. with an ellipse, which we translate into a parameter covariance matrix C_mm, where m refers to either CS or FS. We scale the resulting covariance with the parameter m₁ (see Table 2). We then construct the total covariance matrix 𝒞 from Eq. (13) obtained from the 1000 best-fit (α_⊥, α_∥), adjusting each time the coefficients (according to Eqs. (14) and (15)) to account for the observed errors of a given realisation. The corresponding correlation matrix before the individual adjustments is shown in Fig. 4. Using the combination method, we compute for each mock the consensus data vector x_c and covariance C_c from Eqs. (11) and (12). Here we emphasize that the strong assumption of Gaussian elliptic contours for the parameters (α_∥,α_⊥) of every inference highly depends on the statistic of the considered tracer. The presence of non-Gaussian contours in our sample of mocks introduces biases in the total covariance matrix 𝒞 and each individual consensus data vector x_c and covariance C_c. This is one of the limitations of the GA method that we can avoid with a JS fit.

Fig. 4. Correlation coefficients between α∥ and α⊥ in CS and FS obtained from fits to the 1000 EZmock realisation.

Figure 5 compares the distributions of α_∥ (left) and α_⊥ (right) and and their 1σ uncertainties as estimated by GA (top four panels) and JS (bottom four panels) versus the same quantities obtained from CS and FS fits. We see that the distributions are nicely correlated and scattered around the identity line (dashed). We find that the scatter (especially for errors) is less important between FS and JS than CS and JS. This would indicate that the Fourier space information has slightly more weight than the configuration space in the minimization of the likelihood. Moreover, while the GA errors are almost systematically inferior to the CS and FS ones, the JS errors are scattered on both sides of the dashed lines. Due to statistical fluctuations, JS analysis may result in looser constraints than FS or CS. We believe these fluctuations might origin from the finite number of mock realisations used to build the covariance matrix, though it is hard to test this hypothesis without a larger number of mocks.

Fig. 5. Comparison between the distributions of the best-fit BAO parameters (α∥,α⊥) and their estimated errors obtained by fitting individual mock realisations. The four top (bottom) panels compare the CS (blue), FS (red) distributions with the GA (JS). The dashed lines represent a perfect correlation.

Table 4 summarizes the statistical properties of (α_∥,α_⊥) for the different analyses, CS, FS, GA and JS, performed on 1000 EZmock realisations. For each parameter, we show six quantities: the average bias Δ_α = ⟨α − α_exp⟩ with respect to the expected value (α_⊥,exp = α_∥,exp = 1), the standard deviation of best-fit values σ, the mean estimated error ⟨σ_i⟩, the mean asymmetry of the estimated error distribution $⟨A\left({\sigma }_i\right)⟩=⟨2({\sigma }_i^{\sup }-{\sigma }_i^{\inf })/({\sigma }_i^{\sup }+{\sigma }_i^{\inf })⟩$, the mean of the pull Z_i = (α_i−⟨α_i⟩)/σ_i and its standard deviation. If errors are correctly estimated and follow a Gaussian distribution, we expect that σ = ⟨σ_i⟩, ⟨Z⟩=0 and σ(Z_i)=1. Table 4 also shows the number N_good of valid realisations after removing undetermined likelihoods, extreme values and errors, along with the mean value of the reduced chi-square ${}_r{\chi }_{\min }^2$.

Firstly, one can notice that the number of valid realisations (as defined above) for the JS analysis is larger than the one of the combined analysis GA. The GA method requires both FS and CS fits to converge, so the number of valid realisations for GA is necessarily the intersection of valid FS and CS realisations. By joining the information of both spaces, the JS fit is able to constrain the acoustic scale even on noisy mocks, where FS and CS separately fail. For each fitting procedure, we find the minimum reduced chi-square $⟨{}_r{\chi }_{\min }^2⟩\sim 1$, showing that the majority of the mocks are accurately modeled by our templates.

Secondly, we see good agreement between ⟨σ_i⟩ and σ for both parameters in all analyses. While the JS errors lies in between the CS and the FS errors for α_∥, the combined results present in average smaller errors. Note that the standard deviation σ of the best-fit values are scaled with the appropriate correction factor $\sqrt{m_2}$ (see Table 2). A crucial issue with the GA method is m₂ is ill-defined.

Overall, systematic shifts Δ_α for all methods are below half a per cent. The GA and JS methods yield reduced systematic shifts Δ_α than each individual space alone. While the shifts between GA and JS are the same for α_⊥, the JS method yields slightly smaller shift for α_∥ as also observed in the previous section.

If we assume that σ is a good estimate of the standard deviation of α, then the uncertainty of Δ_α should be ${\sigma }_{\mathrm{\Delta }}=\sigma /\sqrt{N_{\mathrm{good}}}$. In this case, all systematic shifts are smaller than 3σ_Δ except for α_∥ in FS, which reaches a 3.5σ_Δ discrepancy. Shifts for GA are (−1.7, −1.8)σ_Δ for (α_⊥, α_∥) and (−1.5, −1.0)σ_Δ for JS fits. Note that our results are slightly different from those reported in Bautista et al. (2021) and Gil-Marín et al. (2020) due to a few differences in the analyses, such as how we decompose peak and smooth parts of the template, the number of polynomial terms, the values of non-linear damping terms.

Figure 6 displays the distributions of the pull Z_i for α_∥ and α_⊥. For all four methods, the mean of the pull ⟨Z_i⟩ is centered in zero to a few percent level, suggesting again that there are no systematic bias in the estimated α. Moreover the standard deviation of the pull σ(Z_i) is smaller than one for both parameters in all analyses indicating a slight overestimation of the errors. This overestimation is more important for the parameter α_∥, and in JS in general. As the standard deviation of the pulls are only a few percent away from unity, we do not attempt to correct these effects. If any overestimation of uncertainties is real, not correcting for it can be considered as an conservative approach. Yet, we see that for the parameter α_∥ the overestimation of the errors reaches 10% in JS. While this seems to be significant, this result can be a consequence of the Gaussian assumption for the distribution of the parameters.

Fig. 6. Distributions of the pull Zi = (xi−⟨xi⟩)/σi for the parameters (α∥,α⊥) obtained obtained from the CS, FS, GA and JS analyses of the 1000 EZmocks sample. Orange curves describing a Gaussian distribution with zero mean and the standard deviation of the corresponding pull are shown for comparison. The dotted lines are centered on zero.

To investigate the Gaussian nature of the distribution of the α’s parameters, we perform a D’Agostino-Pearson’s test (D’Agostino 1971; D’Agostino & Pearson 1973) on the sample of best-fit parameters using the SCIPY.STATS library⁵. The test combines the high-order statistical moments skewness and kurtosis to test the hypothesis that the two parameters are normally distributed. The statistical quantity K₂ is constructed as a combination of the renormalized skewness and kurtosis such that K₂ asymptotically follows a χ² law with two degrees of freedom. When the skewness and kurtosis simultaneously deviate from 0 (Fisher kurtosis) K₂ gets larger. The p-value for the K₂ statistic is then to be compared with the risk a usually set to 5% for a Pearson test. As the a is the risk to reject the null hypothesis H₀ while it is true, the p-value should be less than a to reject H₀. Here the null hypothesis is that the α’s sample comes from a normal distribution. Table 5 shows the results of the Pearson tests on the α_⊥, and α_∥ distributions. Independently of the method, the Gaussian distribution hypothesis can be rejected for the parameter α_∥ only, for which p/a <  1. This result indicates that for α_∥, the smaller values of σ(Z_i) observed in Table 4 might be partially induced by the non-Gaussianity of the α distributions and not simply by an overestimation of the uncertainties.

Figure 7 shows two dimensional confidence levels for a few realisations of mocks. Red contours are for FS results, blue for CS, green for JS and black for GA. The left panels are cases where the likelihoods can be approximated by a 2D Gaussian distribution while the right panels the opposite. In these examples, we can see how the JS fits are a better description for the non-Gaussian cases, where the final consensus are not necessarily Gaussian.

Fig. 7. Comparison between the different methods, using some EZmock realisations (mock tag above each plot), when the contours found in configuration or Fourier space are Gaussian (left) and non Gaussian (right). Red contours are for FS results, blue for CS, green for JS and black for GA. The expected value is the intersection of the dotted black lines, and the best fit values are described by a star for the JS and GA methods. We can see how the JS method yields better combined results that are not necessarily Gaussian.

5. Application to eBOSS DR16 data

In this section, we apply our methodology to the eBOSS LRG sample, described in Sect. 2.1, letting free only α_∥, α_⊥, the linear density bias b and the broadband parameters.

Figure 8 shows the best-fit models for the power spectrum and correlation function multipoles, for three analyses: FS, CS and JS (note that we do not show a GA best-fit model as GA results are derived using the CS and FS best-fit). The residuals between in the bottom panels show excellent agreement between the JS model and the individual models for FS and CS, particularly over the BAO features. Furthermore, as the global amplitudes driven by the parameter b are sensibly the same, the small differences between the models arise mainly from the different broadbands. This is specially true for the correlation function as it is composed of fewer data bins.

Fig. 8. Best monopole and quadrupole models for the CS, FS and JS analysis. The residues with respect to the measurements are standardized for each point by the errorbar. The grey shaded area correspond to a 1σ difference.

Figure 9 shows the two dimensional 68 and 95% confidence levels for CS, FS, GA and JS. Once again we see an excellent agreement between the four analysis, with the JS giving slightly looser constraints. We also notice small differences in the inclination of the contours. The correlation coefficients for the CS, FS, GA and JS are respectively (−0.395, −0.404, −0.396, −0.445). Those values are consistent with the correlations found with 1000 EZmock realisations: (−0.380, −0.395, −0.383, −0.420). Note that this coefficient for the GA is highly dependent on the quality of the estimated correlation matrix of the parameters estimated from the EZmock sample and shown in Fig. 4. We find a higher correlation coefficient using the JS analysis, resulting in a thinner and steeper contours. This have for effect to increase the errors on the individual parameters without varying much the Figure-of-Merit (area of the contour, see Appendix D for a discussion on the FoM).

Fig. 9. Confidence contours of the parameters posterior distributions for the four analysis. For the CS and JS, the 1σ contours are filled in blue and red respectively. For the JS and GA we give the 1σ and 2σ contours for comparison in green and black. The best values are designated with a star.

Table 6 summarises the results of the four different analysis. The reduced χ² values are 1.28, 1.27 and 1.21 respectively for FS, CS and JS indicating relatively good adjustment to the data without over fitting from the broadband. In each case, the p_value indicates a valid fit (if compared with a risk of 5%). The best fit values for the parameters (α_⊥, α_∥) are all in agreement according to their respective 1σ error bar. For both parameters, the JS analysis gives the larger uncertainty. Rescaling with the fiducial cosmology, we also derive the corresponding physical quantities D_M/r_d and D_H/r_d. Note that our best fit results do not exactly match the BAO constraints obtained in Bautista et al. (2021) and Gil-Marín et al. (2020), neither in configuration space nor in Fourier space, where the GA consensus gives D_H/r_d = 19.33 ± 0.53 and D_M/r_d = 17.86 ± 0.33 (third entry of Table 14 of Bautista et al. 2021). Those discrepancies arise from the slightly different non linear broadening parameters (Σ_∥, Σ_⊥) and a more flexible broadband polynomial expansion (i_min, i_max)=(−2, 1) to be compared to ( − 2, 0) for Bautista et al. (2021) and ( − 1, 1) for Gil-Marín et al. (2020).

6. Conclusions

This work introduced a new BAO analysis of the DR16 eBOSS LRG sample, using Fourier and configuration-space information simultaneously. We compared the joint space (JS) analysis with the commonly used Gaussian approximation (GA) method, which combines results from the two spaces (Fourier and configuration) at the parameter level. The main advantage of the JS method is that it does not require any Gaussian assumption for the likelihood profiles. While the GA method is accurate only if the individual likelihoods to be combined are both Gaussian, yielding only Gaussian posterior distributions.

We assessed the systematic biases and errors of both methods by applying JS and GA to a set of 1000 EZMocks, which reproduce the eBOSS DR16 LRG sample properties. Compared to GA, JS provides a more accurate estimation of the acoustic scale by lowering the systematic shift of α_∥ with respect to its expected value while GA has by construction a better precision (smaller error bars). Moreover, the JS offers a better control over variations of the analysis. Indeed the same mocks as the ones used to estimate the covariance are commonly used to perform systematic studies. Because of this, the standard deviation of any parameter should be rescaled with the correction factor m₂ (see Table 2), which is not properly defined for the GA. One should note that because of larger size of the data vector, the JS methods requires a sufficiently large mock sample to estimate the covariance. However we found that the constraints are quite stable, and the number of mocks needed is not much larger than for a regular configuration or Fourier analysis.

We applied the different analysis to the eBOSS LRG sample and found consistent results (see Table 6). As expected from the statistical study the JS gives slightly looser constraints and a larger correlation coefficient. Despite providing looser constraints on cosmological parameters than the standard GA, we believe that JS is a more robust and reliable method for modelling the clustering signal by correctly accounting for the correlation between configuration and Fourier space and not relying on the Gaussianity of the parameter likelihoods.

One specific feature of the JS in a BAO analysis is that the broadband terms are independent between spaces. Hence, the JS introduces more nuisance parameters to be fitted simultaneously.

We are currently working on extending joint-space fits for measurements of the growth-rate of structures using information from the full shape of the power spectrum and correlation function. This work can also be extended to joint fits between pre and postreconstruction catalogs, as already performed in Fourier space by Gil-Marín (2022).

Acknowledgments

We would like to thank Hector Gil-Marín, Ashley Ross, Elena Sarpa, for useful discussions. The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University – A*MIDEX, a French “Investissements d’Avenir” program (AMX-20-CE-02 – DARKUNI).

References

Appendix A: Fit with hexadecapole

Table A.1 shows the results for the fit of the averaged 1000 EZmocks, with and without using the hexadecapole. For every analysis the use of the additional information results in a slight increase of the systematic bias and of the errors bars. The looser constraints mainly arise from the larger size of the data vector when using hexadecapole. Indeed a larger amount of data bins results in an increase of the Whishart bound through the rescaling parameter m₁ (especially for the JS analysis). Thus we decide not to use the hexadecapole information in this work.

Appendix B: Prereconstruction results

We perform fits on average correlations and individual mocks for the prereconstruction sample. Figure B.1 presents the impact of different parameters choice on the best fit values of α_∥ and α_⊥ for the GA (in black) and JS (in green) analysis. In every plot, a grey shaded area represents the 1 percent deviation from the expected value. Both analysis appears to be more dependant on the range of scales used (particularly in fourier space) as the small scales non linearities are not well modeled. However, without the reconstruction procedure that requires a fiducial cosmology, the dependence on Ω_b and Ω_cdm is reduced. The two analysis gives different results for the parameter α_∥, the JS analysis tends to give a positive bias while the GA gives a negative.

Fig. B.1. Impact of the choice of fitting scales, fiducial cosmology, broadening Σ’s parameters and polynomial broadband order in the recovered values of the parameters α∥ and α⊥. Each point is the best-fit from the average of the 1000 prereconstruction EZmocks. The JS method is in green and GA method is in black. The grey shaded areas correspond to a 1 percent error and red shared areas indicate the fiducial choices of our analysis.

Figure B.2 shows the correlation matrix between the prereconstructed 2PCF and PS. The main differences between prereconstruction and postreconstruction covariance matrices are in the off-diagonal blocks within each space (either FS or CS), which are significantly reduced in postreconstruction. We construct the α’s parameters prereconstruction total covariance matrix obtained from the 1000 best-fit (α_⊥, α_∥). The corresponding correlation matrix before the individual adjustments is shown in Figure B.3. The reconstruction procedure enhances the correlations for a given parameter in two spaces, but reduces the correlations between the two parameters themselves (both in a given space and across spaces). Overall, the measurements depend on the signal-to-noise ratio of the BAO features, which are lower for prereconstruction data.

Table B.1 summarises the statistical properties of (α_∥,α_⊥) for the different analyses, CS, FS, GA and JS, performed on 1000 EZmock realisations. For each parameter, we show six quantities: the average bias Δ_α = ⟨α − α_exp⟩ with respect to the expected value (α_⊥,exp = α_∥,exp = 1), the standard deviation of best-fit values σ, the mean estimated error ⟨σ_i⟩, the mean asymmetry of the estimated error distribution

Fig. B.2. Correlation matrix of the correlation function ξℓ and power spectrum multipoles Pℓ estimated using the 1000 independent measurements of prereconstructed EZmocks.

Fig. B.3. Correlation coefficients between α∥ and α⊥ in CS and FS obtained from fits to the prereconstructed 1000 EZmock realisation.

$⟨A\left({\sigma }_i\right)⟩=⟨2({\sigma }_{\mathrm{i}}^{\sup }-{\sigma }_{\mathrm{i}}^{\inf })/({\sigma }_{\mathrm{i}}^{\sup }+{\sigma }_{\mathrm{i}}^{\inf })⟩$, the mean of the pull Z_i = (α_i−⟨α_i⟩)/σ_i and its standard deviation. If errors are correctly estimated and follow a Gaussian distribution, we expect that σ = ⟨σ_i⟩, ⟨Z⟩=0 and σ(Z_i)=1. For both parameters, the JS analysis presents larger systematic shifts Δ_α and larger errors than the GA. For the four analysis, σ(Z_i) is larger than one, indicating an underestimation of the errors. This effect is more important for the GA. As our modelling is not well suited to describe the non linearities of the prereconstruction two point statistics, for every analysis the precision and accuracy are degraded. Here the JS method takes into account the correlations between two observables that are not well modeled in the first place. However, using an appropriate full shape modeling and removing the broadbands nuisance parameters, the CS and FS should be described by the exact same set of parameters. Hence the JS could improve the analysis.

Appendix C: Modelisation of the cross correlation pattern

We propose an analytical description of the cross correlations between the correlation function and power spectrum multipoles. The 2-point statistic between the fluctuation field in configuration and fourier space is

$$\begin{array}{cc}\hfill ⟨{\delta }_g(\mathit{k}){\delta }_g(\mathit{r})⟩& =\frac{1}{{(2\pi )}^3}{\displaystyle \int d^3}{k}^{\prime }{e}^{i{\mathit{k}}^{\prime }·\mathit{r}}⟨{\delta }_g(\mathit{k}){\delta }_g({\mathit{k}}^{\prime })⟩\hfill \\ \hfill & =e^{i\mathit{k}·\mathit{r}}P(\mathit{k}),\hfill \end{array}$$(C.1)

while assuming gaussian fields, the cross covariance is given by

$$\begin{array}{cc}\hfill Cov[P(\mathit{k}),\xi (\mathit{r})]& =⟨{\delta }_g(\mathit{k}){\delta }_g(-\mathit{k}){\delta }_g(\mathit{r}){\delta }_g(-\mathit{r})⟩\hfill \\ \hfill & =⟨{\delta }_g(\mathit{k}){\delta }_g(\mathit{r})⟩⟨{\delta }_g(-\mathit{k}){\delta }_g(-\mathit{r})⟩+\hfill \\ \hfill & ⟨{\delta }_g(\mathit{k}){\delta }_g(-\mathit{r})⟩⟨{\delta }_g(-\mathit{k}){\delta }_g(\mathit{r})⟩\hfill \\ \hfill & ={2cos(2\mathit{k}·\mathit{r})|P(\mathit{k})|}^2.\hfill \end{array}$$(C.2)

Then, including the shot noise contribution P_N and the volume rescaling, the cross covariance between the multipoles P_l(k) and ξ_l′(r) can be written

$$\begin{array}{cc}\hfill \mathit{Cov}& [P_l(k),{\xi }_{l^{\prime }}(r)]=\frac{2l+1}{2}\frac{2l^{\prime }+1}{2}\frac{1}{\mathit{Vs}}{\displaystyle {\int }_{-1}^1}d{\mu }_k{L}_l({\mu }_k){|P(k,{\mu }_k)+P_N|}^2\hfill \\ \hfill & \times {\displaystyle {\int }_{-1}^1}d{\mu }_r{L}_l({\mu }_r)2cos(2rk{\mu }_k{\mu }_r+2kr\sqrt{1-{\mu }_k^2}\sqrt{1-{\mu }_r^2}).\hfill \end{array}$$(C.3)

Fig. C.1. Analytical modeling of the cross correlations between the correlation function and the power spectrum monopoles for a given range of scales. We do not display the colorbar as we set the shot noise to zero and the volume of the survey to one.

Figure C.1 shows the modeling of the cross covariance features using Eq C.2 for a given range of scales. We find this description satisfying as we only aim for a physical understanding of these patterns, and use the full mocks sample to estimate the covariance. A more reliable modeling would require bin averaged integral as performed in Grieb et al. (2016).

Appendix D: Mock dependant constraints

We tested the robustness of the different analysis against the number of available mock realisations. We fit the averaged stack of mocks for different values of N_mocks used to estimate the sample covariance matrix.

Figure D.1 shows the figure of merit defined as $Fom=1/\sqrt{|C_{\theta }|}$ with C_θ the covariance matrix of the inferred parameters, and the quadratic shift from the expected parameters $\sqrt{{\mathrm{\Delta }}_{{\alpha }_{\Vert }}^2+{\mathrm{\Delta }}_{{\alpha }_{\perp }}^2}$ for different N_mock. As expected we find that, due to the Whishart bound, the JS precision highly decreases when the number of mocks get small, typically the Fom drops when N_mock <  4N_bins. However the JS constraint is compatible with the GA at 1σ, and appears to be robust as the Fom does not decrease very much in the range 1000-600 N_mock. While the GA analysis introduces an additional dependence on N_mock through the covariance C_c defined in Eq 11 (which has to be calculated for each point using the set of best fit parameters found with the appropriated sample covariance), the Fom is only decreased by a few percent when using 100 realisations.

We find that for the JS the quadratic shift from the expected parameters monotonically decreases as the amount of mocks increases. For a large enough amount of realisations, the JS gives the smaller systematic bias.

Fig. D.1. Figure of merit (Fom) and quadratic shift from the expected parameters resulting from the fit of the stack of EZmocks for different Nmock used to estimated the sample covariance. Errors for the α’s are scaled with $\sqrt{N_{\mathrm{mock}}}$ and errors for the Fom are derived from Eq 68 of Taylor et al. (2013).

All Tables

All Figures

	Fig. 1. Correlation matrix of the correlation function ξℓ and power spectrum multipoles Pℓ estimated using the 1000 independent measurements of EZmocks.
In the text

Fig. 2. Precision matrix of the correlation function ξℓ and power spectrum multipoles Pℓ. Top panel: normalized estimated precision matrix $R_{\mathit{ij}}^{\xi P}={\mathrm{\Psi }}_{\mathit{ij}}/\sqrt{{\mathrm{\Psi }}_{\mathit{ii}}{\mathrm{\Psi }}_{\mathit{jj}}}$, obtained from 1000 EZmock realisations. Bottom panel: slices of $R_{\mathit{ij}}^{\xi P}$ showing the cross-correlation between ξℓ and Pℓ for particular values of scales as shown by the green arrows in the left panel. Errors are estimated from Eq. (29) of Taylor et al. (2013).

In the text

Fig. 3. Impact of the choice of fitting scales, fiducial cosmology, broadening Σ’s parameters and polynomial broadband order in the recovered values of the parameters α∥ and α⊥. Each point is the best-fit from the average of the 1000 EZmocks. The JS results are in green and GA in black. The grey shaded areas correspond to a 1% error and red shared areas indicate the fiducial choices of our analysis.

In the text

	Fig. 4. Correlation coefficients between α∥ and α⊥ in CS and FS obtained from fits to the 1000 EZmock realisation.
In the text

	Fig. 5. Comparison between the distributions of the best-fit BAO parameters (α∥,α⊥) and their estimated errors obtained by fitting individual mock realisations. The four top (bottom) panels compare the CS (blue), FS (red) distributions with the GA (JS). The dashed lines represent a perfect correlation.
In the text

	Fig. 6. Distributions of the pull Zi = (xi−⟨xi⟩)/σi for the parameters (α∥,α⊥) obtained obtained from the CS, FS, GA and JS analyses of the 1000 EZmocks sample. Orange curves describing a Gaussian distribution with zero mean and the standard deviation of the corresponding pull are shown for comparison. The dotted lines are centered on zero.
In the text

Fig. 7. Comparison between the different methods, using some EZmock realisations (mock tag above each plot), when the contours found in configuration or Fourier space are Gaussian (left) and non Gaussian (right). Red contours are for FS results, blue for CS, green for JS and black for GA. The expected value is the intersection of the dotted black lines, and the best fit values are described by a star for the JS and GA methods. We can see how the JS method yields better combined results that are not necessarily Gaussian.

In the text

	Fig. 8. Best monopole and quadrupole models for the CS, FS and JS analysis. The residues with respect to the measurements are standardized for each point by the errorbar. The grey shaded area correspond to a 1σ difference.
In the text

	Fig. 9. Confidence contours of the parameters posterior distributions for the four analysis. For the CS and JS, the 1σ contours are filled in blue and red respectively. For the JS and GA we give the 1σ and 2σ contours for comparison in green and black. The best values are designated with a star.
In the text

Fig. B.1. Impact of the choice of fitting scales, fiducial cosmology, broadening Σ’s parameters and polynomial broadband order in the recovered values of the parameters α∥ and α⊥. Each point is the best-fit from the average of the 1000 prereconstruction EZmocks. The JS method is in green and GA method is in black. The grey shaded areas correspond to a 1 percent error and red shared areas indicate the fiducial choices of our analysis.

In the text

	Fig. B.2. Correlation matrix of the correlation function ξℓ and power spectrum multipoles Pℓ estimated using the 1000 independent measurements of prereconstructed EZmocks.
In the text

	Fig. B.3. Correlation coefficients between α∥ and α⊥ in CS and FS obtained from fits to the prereconstructed 1000 EZmock realisation.
In the text

	Fig. C.1. Analytical modeling of the cross correlations between the correlation function and the power spectrum monopoles for a given range of scales. We do not display the colorbar as we set the shot noise to zero and the volume of the survey to one.
In the text

	Fig. D.1. Figure of merit (Fom) and quadratic shift from the expected parameters resulting from the fit of the stack of EZmocks for different Nmock used to estimated the sample covariance. Errors for the α’s are scaled with $\sqrt{N_{\mathrm{mock}}}$ and errors for the Fom are derived from Eq 68 of Taylor et al. (2013).
In the text