Baryon acoustic oscillations from a joint analysis of the large-scale clustering in Fourier and configuration space

Baryon acoustic oscillations (BAOs) are a powerful probe of the expansion history of our Universe and are typically measured in the two-point statistics of a galaxy survey, either in Fourier space or in configuration space. In this work, we report a first measurement of BAOs from a joint fit of power spectrum and correlation function multipoles. We tested our new framework with a set of 1000 mock catalogs and showed that our method yields smaller biases on BAO parameters than individually fitting power spectra or correlation functions, or when combining them with the Gaussian approximation method. Our estimated uncertainties are slightly larger than those from the Gaussian approximation, likely due to noise in our sample covariance matrix, the larger number of nuisance parameters, or the fact that our new framework does not rely on the assumption of Gaussian likelihoods for the BAO parameters. However, we argue that our uncertainties are more reliable since they rely on fewer assumptions, and because our method takes correlations between Fourier and configuration space at the level of the two-point statistics. We performed a joint analysis of the luminous red galaxy sample of the extended baryon oscillation spectroscopic survey (eBOSS) data release 16, obtaining $D_H/r_d = 19.27 \pm 0.48$ and $D_M/r_d = 17.77 \pm 0.37$, in excellent agreement with the official eBOSS consensus BAO-only results $D_H/r_d = 19.33 \pm 0.53$ and $D_M/r_d =17.86 \pm 0.33$.


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 on D H /r d , D M /r d and f σ 8 , 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 σ 8 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.

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.

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 2 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 2 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.

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 spectroincompleteness 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.

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

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 A80, page 2 of 15 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 available1 .
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 −1 Mpc for power-spectrum and 30 < r < 180 h −1 Mpc for correlation function.
The starting point is a linear-theory-based model for the redshift-space anisotropic galaxy power-spectrum P(k, µ), 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 Σ 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 −k2 Σ 2 r /2 for the postreconstruction model and S (k) = 0 for the prereconstruction, where Σ r = 15 h −1 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 powerspectrum 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 −1 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 2 (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 ): The correlation function multiples ξ (r) are obtained by Hankel transforming the P (k): where j are the spherical Bessel functions.These transforms are computed using Hankl3 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 = k /α and k ⊥ = k /α ⊥ .Therefore, the observed power spectrum relates to the true power spectrum as We note that this volume scaling does not apply for the correlation function: 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 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 A80, page 3 of 15 A&A 667, A80 (2022) 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: where a P ,i and a ξ ,i 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 0 , ξ 0 and quadrupoles P 2 , ξ 2 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.

Parameter inference
The cosmological parameter inference is performed by means of the likelihood analysis of the data.The likelihood L is defined such that 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 Ψ = (1 − D) Ĉ−1 is obtained from the unbiased estimate of the covariance from 1000 realisation of EZmocks: where X k i 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 χ 2 of Eq. ( 10) using a quasi-Newton minimum finder algorithm iMinuit 4 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 χ 2 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 χ 2 values in two dimensions, looking for the contours yielding a ∆χ 2 = 2.3 or 5.9 for 68 and 95% confidence levels, respectively.

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 χ 2 between individual measurements is the same as the one from the consensus result.The expression for the combined covariance matrix is and the combined data vector is where C mn is a p × p block from the full covariance matrix between all parameters and methods C, containing pM × pM elements, defined as 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 C from all the realisations of x m for each method.We obtain the correlation coefficients ρ mocks p 1 ,p 2 ,m,n , that is the covariance C normalized by its diagonal elements, between parameters p 1 , p 2 and methods m, n.We scale these coefficients by the diagonal errors from the data, to obtain the final matrix C for the data.
It is worth emphasizing that the matrix C 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 1 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 ρ data max and for the ensemble of mocks ρ mocks max .We can write the adjusted correlation coefficients for one particular realisation as The equation above accounts for the diagonal terms of the offdiagonal block C mn .For the off-diagonal terms, we use 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 C, 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.

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 ξ 0 , ξ 2 , P 0 , P 2 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 Ĉ.As in Bautista et al. (2021) and Gil-Marín et al. (2020), we apply the following scale cuts: r ∈ [50, 150] h −1 Mpc and k ∈ [0.02, 0.3] h −1 Mpc. Figure 1 shows the resulting correlation matrix, defined as R = σ −1 Ĉ[σ −1 ] T , where σ is the vector composed by the square-root of the diagonal elements of Ĉ.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 monopolemonopole and quadrupole-quadrupole present the same nonlinear features.A given physical scale r is correlated (and anticorrelated) 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.
We then need an estimate for the precision matrix Ψ in Eq. ( 10).Inverting the covariance mixes the different modes and scales across CS and FS.The top panel of Fig. 2 presents A80, page 5 of 15 A&A 667, A80 (2022) 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 ξ 2 and the FS P 0 and P 2 , respectively by ploting the amplitude of the precision matrix at two fixed scaled r min = 52.5 h −1 Mpc and r max = 147.5 h −1 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, ξ 2 is weakly correlated with P 0 while being strongly correlated with P 2 .This result is in agreement with the correlations between ξ 0 and ξ 2 .Furthermore, for large r, while the number of correlated k modes increases, the amplitude of the correlations decreases distinctly for ξ 2 − P 2 , and increases for ξ 2 − P 0 .Note that the correlations between the ξ 0 and P are not shown here but behave in a similar fashion.
While capturing additional information about the scale dependence of the joint space correlations, our new methodology has some drawbacks.Indeed, because the covariance matrix Ĉ 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: with D the Hartlap factor defined in Sect.3.3 and 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 1 is to be directly applied to the estimated parameter covariance matrix for a given measurement, while the factor m 2 scales the standard deviation of a given parameter over a set of mocks.The values of the parameters m 1 and m 2 are given in Table 2. Since the enlargement of the parameter constraints ( √ 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.

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

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.Notes.m 1 is the factor to be applied to the estimated covariance matrix of the parameters and m 2 is the factor that scales the scatter of best-fit parameters of a set of mocks (if these were used in the calculation of the covariance matrix).N mock is the number of mocks used in the estimation of the covariance matrix, N par is the total number of parameters fitted and N bins is the total size of the data vector.The derivation of m 1 and m 2 can be found in Percival et al. (2014).
At first, we model ξ and P in the separation ranges r ∈ [50, 150] h −1 Mpc and k ∈ [0.02, 0.3] h −1 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 −1 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 fid d 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.
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

-No
Notes.See Eq. ( 1) and text for detailed description of the parameters.
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.

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 1 (see Table 2).We then construct the total covariance matrix C 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 C 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.
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 A80, page 7 of 15 A&A 667, A80 ( 2022 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.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 , 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 χ 2 min .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 χ 2 min ∼ 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 √ m 2 (see Table 2).A crucial issue with the GA method is m 2 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 σ ∆ = σ/ N good .In this case, all systematic shifts are smaller than 3σ ∆ except A80, page 8 of 15 T. Dumerchat and J. E. Bautista: BAO from a joint Fourier and configuration-space analysis Table 4. Statistics on the fit of the 1000 EZmocks realisations.Notes.N good is the number of valid realisations after removing undefined contours and extreme values and errors.We show the mean value of the best-fit reduced r χ 2 min .For each parameter, we show the average bias ∆ α ≡ α i − α exp , the standard deviation of best-fit values σ ≡ α 2 i − α i 2 , the average of the per-mock estimated uncertainties σ i , the asymmetry of the estimated error distribution where σ sup i and σ inf i are the superior and inferior one sigma errors estimated from the likelihood profile, the average of the pull Z i ≡ (α i − α i )/σ i and its standard deviation σ(Z i ). 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) andGil-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 s is a combination of the skewness and kurtosis normalized coefficients.The p value of the test is a 2-sided chi squared probability for the Gaussianity hypothesis test.The parameter a that we set at 5% is the risk (usually named α in the literature) that we reject the hypothesis H 0 of normality while it is true.
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.
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.statslibrary 5 .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 2 is constructed as a combination of the renormalized skewness and kurtosis such that K 2 asymptotically follows a χ 2 law with two degrees of freedom.When the skewness and kurtosis simultaneously deviate from 0 (Fisher kurtosis) K 2 gets larger.The p-value for the K 2 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 0 while it is true, the p-value should be less than a to reject H 0 .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  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.
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.

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.
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).
Table 6 summarises the results of the four different analysis.The reduced χ 2 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) andGil-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).

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 2 (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).

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

Fig. 2 .
Fig. 2. Precision matrix of the correlation function ξ and power spectrum multipoles P .Top panel: normalized estimated precision matrixR ξP i j = Ψ i j / Ψ ii Ψ j j, obtained from 1000 EZmock realisations.Bottom panel: slices of R ξP i j 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) ofTaylor et al. (2013).

Fig. 3 .
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.

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

Fig. 5 .
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.

Fig. 6 .
Fig.6.Distributions of the pull Z i = (x i − x i ) /σ 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.
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.

Table 1 .
Fiducial cosmologies used in this work.All models are parameterized by their fraction of the total energy density in form of total matter Ω m , cold dark matter Ω cdm , baryons Ω b , and neutrinos Ω ν , the Hubble constant h = H 0 /(100 km s −1 Mpc −1 ), the primordial spectral index n s and primordial amplitude of power spectrum A s of scalar perturbations.With these parameters we compute the normalisation of the linear power spectrum σ 8 at z = 0 and the comoving sound horizon scale at drag epoch r d .

Table 2 .
Correction factors for the three analysis performed in this work.

Table 3 .
Parameters used in the fits, with their flat priors.

Table 5 .
Results of D'Agostino and Pearson's test of normality over the α's distributions.
A80, page 10 of 15 T. Dumerchat and J. E. Bautista: BAO from a joint Fourier and configuration-space analysis 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.