© The Authors 2026

1. Introduction

Over the past decades, the study of the statistical properties of the large-scale structure (LSS) of our Universe has led to major advances in our understanding of the cosmic amounts of dark energy and dark matter as well as the initial conditions that were the origin of such structures. Larger and denser maps of the LSS have been built thanks to the increasing power of multiplex fiber-fed spectroscopic surveys such as the cosmological programs of the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2011; Blanton et al. 2017) and the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016a,b, 2022). Such surveys measure millions of high precision galaxy redshifts and span large fractions of our cosmic history. These maps enable estimates of the galaxy two-point functions, such as the power spectrum (in Fourier space) or the correlation function (in configuration space). The shape of such statistical measurements heavily depends on the expansion history of the Universe and on the growth history of initial cosmological perturbations.

The observed two-point functions of a galaxy redshift survey is anisotropic with respect to the line-of-sight (LOS) due to redshift-space distortions (RSD). Such distortions appear when converting redshifts to distances: peculiar velocities modify the observed redshifts in a coherent manner. On large scales, the velocity field is simply a response to the growth of density perturbations. RSD analyses aim to measure such anisotropies in the two-point functions and provide constraints on the logarithmic growth rate of structure f ≡ dlnD(a)/dlna (where D(a) is the growth factor). Measurements of the growth rate of structures are a powerful cosmological probe of gravity and are complementary to measurements of the expansion rate from type-Ia supernovae or baryon acoustic oscillations (BAOs).

Modeling anisotropies on large linear scales (above ∼80 Mpc) is straighforward (Kaiser 1987), but on intermediate to small scales (below ∼50 Mpc) it can be a challenge due to nonlinear clustering and a more complex galaxy-dark-matter connection. Recently, effective field theories of large-scale structures (EFTofLSS; Ivanov 2022) have been able to model intermediate scale two-point clustering in redshift space. Such theories incorporate effective coefficients that depend on small-scale physics, which are then marginalised over.

Traditionally, RSD measurements have been conducted either in configuration space (Bautista et al. 2021; Hou et al. 2021; Tamone et al. 2020) or in Fourier space (Gil-Marín et al. 2020; de Mattia et al. 2020; Neveux et al. 2020). In configuration space, the two-point correlation function, ξ(r), is estimated with pair-counting algorithms. In Fourier space, the power spectrum, P(k), is estimated by placing galaxies on a regular mesh and performing fast Fourier transforms. Both approaches use the same initial dataset and, in principle, should yield the same cosmological results. In practice, however, they differ in their noise properties, on how they are impacted by systematic effects, and on the choice of scales used for inference¹. Therefore, Fourier and configuration space results are highly correlated but not identical, so there is a small gain in combining both results into a single consensus result.

Sánchez et al. (2017) proposed a method to compute consensus results between Fourier and configuration space analyses, which works at the level of the parameter posteriors. The method assumes that individual parameter posteriors are Gaussian and also yield Gaussian consensus posteriors. This assumption can be broken with noisier data where posteriors in each individual analysis space are non-Gaussian.

Dumerchat & Bautista (2022) developed an alternative framework to obtain consensus results, which is based on a joint fit of correlation function and power spectrum multipoles and accounts for their covariance at the inference level. Such a methodology can yield non-Gaussian consensus posteriors on the fitted parameters, so it relies on less assumptions than the method by Sánchez et al. (2017). This methodology was tested on BAO measurements of eBOSS LRG data and yielded consistent results with the official analysis.

In this work, we extend the joint configuration and Fourier space analysis by Dumerchat & Bautista (2022) to an RSD analysis. This paper is organized as follows. Section 2 describes the real and synthetic datasets employed in our analysis, followed by a presentation of different consensus methodologies in Sect. 3. The theoretical model of the two-point functions used in our inferences is described in Sect. 4. In Sect. 5 we present the results of the joint RSD analysis on mock data. Results on eBOSS LRG data are discussed in Sect. 6.

2. Dataset

In this section we summarize the datasets used in our RSD analyses. We also describe the synthetic samples employed to validate our methodology.

2.1. The BOSS and eBOSS sample of LRGs

The third and fourth generations of the SDSS contained two major cosmological surveys: the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013) and its extension, eBOSS (Dawson et al. 2016), each program observed for about five years each. The SDSS used the 2.5-meter Sloan Foundation Telescope at Apache Point Observatory (APO; Gunn et al. 2006), which is equipped with a focal plane able to host 1000 optical fibers and two spectrographs covering the 3600–10 400 Å range in wavelengths (Smee et al. 2013). The main goal of BOSS and eBOSS was the measurement of BAO and RSD to constrain cosmological parameters (Alam et al. 2017, 2021). For that, these surveys measured millions of galaxies and quasars over a large range in redshift (0 < z < 4). These samples are publicly available as part of SDSS Data Release 16 (Ahumada et al. 2020).

To apply our framework, we consider the sample of luminous of red galaxies (LRG) from the combination of BOSS and eBOSS redshits over 0.6 < z < 1.0. The full description of the construction of this sample for the purpose of clustering measurements can be found in Ross et al. (2020). This sample is composed of 174 816 eBOSS LRG redshifts over 4242 deg² and BOSS LRGs over 9493 deg² over the same redshift range, yielding a total of 377 458 galaxies. Several weights are computed and assigned to each galaxy to correct for spurious density fluctuations originating from target selection, spectroscopic completeness, fiber collisions, among others. The window function of the survey is characterized by a set of unclustered points, the random catalog, containing about 20 times the number of galaxies. This sample of LRG was used to measure BAO and RSD, both in configuration and Fourier space (Bautista et al. 2021; Gil-Marín et al. 2020).

2.2. Synthetic datasets

We employed two types of synthetic datasets in this work: N-body simulations and approximate methods. N-body simulations are fully nonlinear and are a perfect test bench for our theoretical clustering models (Sect. 4). Approximate mocks have less realistic clustering on small-scales but easily reproduce the survey geometry, the density of tracers and make possible the production of thousands of realisations, which are useful to test the statistical properties of our inferences and estimate sample covariance matrices.

2.2.1. N-body simulations:

The ABACUSSUMMITN-body simulation suite (Garrison et al. 2019) contains 25 independent realisations in a periodic cubic box, each covering a volume of 8 h⁻³ Gpc³. These mocks use the ABACUS code (Garrison et al. 2021) to evolve 6912³ matter particles with mass of 2 × 10⁹ h⁻¹ M_⊙ from initial conditions set at z = 99 from second-order Lagrangian perturbation theory. Halos are identified from groups of particles using the specialized spherical-overdensity-based halo finder CompaSO (Hadzhiyska et al. 2022). We used the z = 0.8 snapshot to match the redshift range of the LRG sample. The halos were populated with LRGs using a halo occupation distribution model (HOD; Yuan et al. 2023) matching the clustering properties of DESI Early Data Release (DESI Collaboration 2024). To obtain the covariance matrix of the clustering measurements of our N-body simulations, we used a set of 1000 realizations of (2 h⁻¹ Gpc)³ boxes run with the EZMOCK method (Chuang et al. 2015), tuned to match the clustering of the n-body simulations. Note that these EZmocks are different than those described next.

2.2.2. Approximate mocks:

We used 1000 realizations of the LRG eBOSS+BOSS survey geometry generated with the EZMOCK method (Chuang et al. 2015). This approach employs the Zel’dovich approximation to model the density field at a given redshift and populate it with galaxies. The method is computationally efficient and has been specifically calibrated to reproduce the two- and three-point statistics of the galaxy sample, including mildly nonlinear scales. The angular and redshift distributions of the eBOSS LRG sample, combined with the z > 0.6 CMASS sample, were accurately replicated in these mock catalogs. A detailed description of the EZMOCKS LRG samples can be found in Zhao et al. (2021). These mock catalogs were crucial to the analysis, enabling the estimation of error covariance matrices for clustering measurements, the evaluation of systematic observational effects on cosmology, and the assessment of correlations between different methodologies for deriving consensus results. In our analysis, we consider two fiducial cosmologies. Table 1 presents the fiducial cosmology used for DESI and eBOSS.

2.3. Clustering measurements

We used the Yamamoto estimator to estimate the power spectrum multipoles (Yamamoto et al. 2006). These are computed using PYPOWER². The survey geometry is described by a window function that defines the observed regions and accounts for selection effects. We apply it by directly convolving the model power spectrum with the window function in Fourier space, following Beutler et al. (2019), Beutler & McDonald (2021). The window function is computed with PYPOWER. The correlation function multipoles are estimated using the Landy–Szalay estimator (Landy & Szalay 1993), and are computed using PYCORR³.

3. Methodology for consensus results

In this section we describe our framework to derive consensus results from a joint RSD analysis in Fourier and configuration space. We first introduce the inference algorithms. Second, we remind the basics of the Gaussian approximation method to obtain consensus results. Lastly, we describe our joint analysis in detail.

3.1. Parameter inference

In this work, we are ultimately interested in the inference of the growth rate of structure f, or its combination with the amplitude of linear perturbations σ₈, as well as the dilation parameters q_⊥, q_∥ describing deviations from the assumed fiducial cosmology. These parameters are defined in Sect. 4.

The likelihood function ℒ of observing a data vector given a model obtained from a set of parameters θ can be written as

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

where N_bins is the dimension of the data vector, Δ_i(θ)≡d_i − m_i(θ) is the data-model difference in bin i. The $\widehat{\mathrm{\Psi }}$ is an estimate of the precision matrix, computed as $\widehat{\mathrm{\Psi }}=(1-D_{\mathrm{H}}){\widehat{C}}^{-1}$, where $\widehat{C}$ is an estimate of the data vector covariance matrix and (1 − D_H) is a correction factor, described next.

In this work, our data vector is either the correlation function or power spectrum multipoles. In our joint analysis, the data vector is a concatenation of the multipoles from both spaces. Our covariance matrices C are estimated from a set of synthetic datasets (also referred to as mocks hereafter) reproducing our survey geometry and clustering, so accounts for both cosmic variance and shot-noise contributions. From N_mock mock measurements of two-point functions, the sample covariance is simply defined as

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{ij}}=\frac{1}{N_{\mathrm{mock}}-1}{\displaystyle \sum _k^{N_{\mathrm{mock}}}}(d_i^k-\overline{d_i})(d_j^k-\overline{d_j}),\end{array}$$(2)

where $\overline{d_i}$ is the sample average of the data vector over all mocks.

Since we use a finite number of mocks, the correction factor takes the form D_H ≡ (N_bins + 1)/(N_mock − 1) to correct for the biased estimate of the precision matrix (Hartlap et al. 2007). Also, each element of $\widehat{C}$ is affected by its own uncertainty which shall be appropriately propagated to the parameter uncertainties. Following Percival et al. (2014), this correction is done by scaling factors given by

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

$$\begin{array}{cc}\hfill m_2& =\frac{m_1}{1-D_{\mathrm{H}}},\hfill \end{array}$$(4)

where D_H is the Hartlap factor defined above, and

$$\begin{array}{cc}\hfill A& \equiv \frac{2}{(N_{\mathrm{mock}}-N_{\mathrm{bins}}-1)(N_{\mathrm{mock}}-N_{\mathrm{bins}}-4)}\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill B& \equiv \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}$$(6)

The factor m₁ is applied directly to the estimated covariance matrix of the parameters for a specific measurement. Conversely, the factor m₂ is used to scale the standard deviation of a parameter across a series of mocks. The values of the parameters m₁ and m₂ of our baseline analysis are given in Table 2 for configuration space (CS), Fourier space (FS) and joint space (JS).

We assumed the covariance matrix to be independent of model parameters and it remains fixed during the inference. Therefore, the c term in Eq. (1) is simply an additive constant. We use two approaches for parameter estimation: a frequentist and a Bayesian one.

3.1.1. Frequentist Approach:

Best-fit parameters are determined by minimizing χ² using a quasi-Newton minimizer implemented in IMINUIT⁴. Parameter uncertainties are estimated by identifying intervals where χ² increases by one unit, corresponding to a 68% confidence level. This approach does not assume Gaussianity of the parameter posterior distribution, so uncertainties can be asymmetric.

3.1.2. Bayesian Approach:

The posterior probability distribution of the parameters θ given the observed data d is given by the Bayes’ theorem, P(θ|d) = ℒ(d|θ)P(θ)/P(d), where ℒ is the likelihood from Eq. (1), P(θ) is the prior and P(d) is the evidence. The likelihood is explored using a Markov chain Monte Carlo (MCMC) algorithm, implemented in EMCEE⁵ (Foreman-Mackey et al. 2013). We run chains in parallel and assess their convergence using the Gelman-Rubin criterion (Gelman & Rubin 1992).

For most test cases, we checked that both approaches give consistent results. The first approach is faster and used to quickly obtain constraints of a small subset of parameters, which can be run on large samples of mocks, while the second one is slower but accurately describes the correlated posterior distributions of all parameters.

We used the DESILIKE⁶ framework, a public tool developed by the DESI collaboration that integrates theory modeling, observable definitions, likelihood computation, theory emulation, and posterior sampling and profiling. These tools were employed in the main cosmological results from DESI (DESI Collaboration 2025a,b,c; DESI Collaboration 2025d,f).

3.2. Consensus results from the Gaussian approximation method

Sánchez et al. (2017) proposed a method to obtain consensus results from configuration-space (CS) and Fourier-space (FS) analyses. This method combines posteriors distributions of the parameters of interest, assuming they are Gaussian and also yielding Gaussian consensus porsteriors. We refer to this method as the Gaussian approximation (GA) hereafter.

In the case of RSD analysis, the parameters of interest are the normalized growth rate of structures fσ₈ and the dilation parameters q_⊥, q_∥. Each analysis produces its own posterior which is approximated by a three-dimensional Gaussian centered in θ = [fσ₈, q_∥, q_⊥] with a 3 × 3 parameter covariance C_θ. The consensus result is obtained by constructing a 6 × 6 parameter covariance, where the off-diagonal blocks (i.e., the cross-covariance between parameters in CS and FS) are estimated from mock catalogs. The consensus result is given by θ^GA and its 3 × 3 covariance C_θ^GA.

Furthermore, Bautista et al. (2021) argued that the cross-covariance block obtained from a set of mocks is representative of the ensemble, but not of a particular realisation. An adjustement of these off-diagonal blocks was proposed, based on the uncertainties of the data realisation itself. This technique was used to obtain consensus between configuration and Fourier space analyses as well as between BAO and RSD constraints (each yielding dilation parameters). A full description of the method we use in this work can be found in Bautista et al. (2021).

3.3. Consensus results from a joint space analysis

An alternative to the Gaussian approximation for obtaining consensus results is to perform a joint fit of both configuration-space and Fourier-space two-point functions, using a single set of parameters. The cross-covariance between multipoles in CS and FS is estimated from mock catalogs. We refer to this approach as the Joint space (JS), hereafter. Dumerchat & Bautista (2022) investigated the performance of JS fits in the context of BAO analyses. In this work, we extend that methodology to RSD analyses.

In addition to the cosmological parameters of interest (fσ₈, q_⊥, q_∥), our models contain several physical parameters which are simultaneously adjusted and marginalised over. Ideally, these nuisance parameters should be consistent across CS and FS analyses. In practice, we see discrepancies on the constraints of such parameters when analysing mock catalogs in either CS or FS. These discrepancies motivated us to create versions of the JS analysis where some of the nuisance parameters are not shared, i.e., each space contains separate versions of the same nuisance parameter. We refer to this model as JS_sep or simply JS, as it is our baseline choice.

Table 3 summarizes the variations considered in our JS analysis. The parameters are described in Sect. 4.

4. Clustering model

We use a Lagrangian perturbation theory (LPT) model to describe the observed clustering in both configuration and Fourier space. This model is based on the Lagrangian formalism developed in Chen et al. (2020, 2021). A more detailed presentation is provided in Appendix B.

4.1. EFTofLSS and LPT in redshift space

In the Lagrangian picture, cosmological structure formation is modeled by tracking the trajectories of fluid elements, rather than evolving the density and velocity fields directly. The key quantity in LPT is the displacement field Ψ(q, η), which maps the initial Lagrangian coordinates q (at some early conformal time η₀) to Eulerian positions x(q) at time η:

$$\begin{array}{c}\hfill \mathbf{x}(\mathit{q},\eta )=\mathit{q}+\mathbf{\Psi }(\mathit{q},\eta ),\end{array}$$(7)

with the initial condition Ψ(q, η₀) = 0. The field Ψ(q, η) fully captures the evolution of the fluid element labeled by q. At first order, the displacement is given by the Zel’dovich approximation (Zel’dovich 1970):

$$\begin{array}{c}\hfill {\mathbf{\nabla }}_{\mathit{q}}{\mathbf{\Psi }}^{(1)}=-D_{+}(\eta ){\delta }_0(\mathit{q}),\end{array}$$(8)

or equivalently,

$$\begin{array}{c}\hfill \mathbf{\Psi }(\mathit{q},\eta )=i{D}_{+}(\eta ){\displaystyle \int \frac{d^{3}k}{{(2\pi )}^3}}\frac{\mathit{k}}{k^2}{\delta }_0(\mathit{k})e^{i\mathit{k}·\mathit{q}},\end{array}$$(9)

where D₊(η) is the linear growth factor and δ₀(q) the initial density contrast.

In this work, we use the Effective Field Theory of Large-Scale Structure (EFTofLSS) in its Lagrangian formulation. This involves smoothing the displacement field and systematically incorporating small-scale physics through an effective expansion (see Eq. (B.8)).

This leads to the Zel’dovich power spectrum:

$$\begin{array}{c}\hfill P_{\mathrm{Zel}}(\mathit{k})={\displaystyle \int d^3}\mathit{q}{e}^{i\mathit{k}·\mathit{q}}[e^{-k_i{k}_j[{\mathrm{\Psi }}_{\mathit{ij}}(\mathbf{0})-{\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})]}-1],\end{array}$$(10)

where Ψ_ij(q) is related to the linear power spectrum via

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})=\frac{1}{{(2\pi )}^3}{\displaystyle \int d^3}k\frac{k_i{k}_j}{k^4}{P}_{\mathrm{lin}}(k)e^{-i\mathit{k}·\mathit{q}}.\end{array}$$(11)

Equation (10) describes the nonlinear mapping of the initial linear power spectrum into the evolved power spectrum via LPT.

Redshift-space distortions (RSD) are incorporated in the Lagrangian formalism by modifying displacements along the line-of-sight (LOS) direction $\widehat{n}$, accounting for galaxy peculiar velocities:

$$\begin{array}{c}\hfill \mathbf{\Psi }\to {\mathbf{\Psi }}^s=\mathbf{\Psi }+\frac{\widehat{n}(\mathit{v}·\widehat{n})}{\mathcal{H}},\end{array}$$(12)

where v is the peculiar velocity and Ψ^s denotes the redshift-space displacement vector.

Finally, the density of biased tracers is modeled using a Lagrangian bias prescription, where the tracer density is given by a bias functional of the initial fields. This treatment is detailed in Appendix B.3.

The final expression of the theoretical model for the redshift-space galaxy power spectrum, accounting for bias in LPT and EFTofLSS, can be found in Chen et al. (2021), Vlah et al. (2016), Carlson et al. (2013). In a compact form it is written

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_s(\mathit{k},\mu )& =P_{\mathrm{s},\mathrm{zel}}(\mathit{k})+P_{\mathrm{s},1-\mathrm{loop}}(\mathit{k})\hfill \\ \hfill & +k^2({\alpha }_0+{\alpha }_2{\mu }^2+{\alpha }_4{\mu }^4)P_{\mathrm{s},\mathrm{zel}}(\mathit{k})\hfill \\ \hfill & +R_h^3(1+{\sigma }_2{k}^2{\mu }^2+{\sigma }_4{k}^4{\mu }^4)\hfill \end{array},\end{array}$$(13)

where P_s, zel(k) is the linear in LPT (Zeldovich) power spectrum, P_{s, 1 − loop}(k) is the 1-loop LPT power spectrum, α_i are the counterterms, and σ_i represent stochastic contributions, which scale with the typical size of halo/galaxy formation R_h. The parameter $\mu =\widehat{k}·\widehat{n}$ represents the cosine of the angle between the wavevector k and the line of sight. In what follows, the stochastic terms will be simplified using sn_i notation, where sn₀ = R_h³, sn₂ = R_h³σ₂, and sn₄ = R_h³σ₄. Counterterms and stochastic contributions include components such as the shot noise and finger-of-God (FoG) effects.

In our work, the model is computed using the VELOCILEPTORS⁷ package, integrated within DESILIKE⁸. The Fourier-based VELOCILEPTORS code implements the redshift-space power spectrum at one-loop order in Lagrangian Perturbation Theory (LPT) with resummation up to a scale k_IR (Chen et al. 2021). In addition to the three stochastic parameters and three counterterm parameters, the bias model includes three parameters: b₁, b₂, and b_s. We decided to fix b₃ to zero as it is expected to be small in our analysis and is quite degenerate with the counterterms (DESI Collaboration 2025c; Maus et al. 2024).

When performing cosmological inference each likelihood evaluation requires calculating loop-correction terms, which involve computationally expensive 2D integrals (∼1s per integral). Since these calculations are repeated tens of thousands of times, a fourth-order Taylor expansion emulator is employed to accelerate the process. This emulator interpolates pre-evaluated grid models, avoiding repeated loop integral evaluations and significantly reducing computation time. More information about the VELOCILEPTORS model can be found in DESI Collaboration (2025c). We also speed up our MCMC fits by analytically marginalizing over the linear nuisance parameters in our model, i.e., the parameters of the stochastic and counterterm.

The choice of priors, particularly for nuisance parameters, can influence cosmological posteriors through two main effects: the prior weight effect (PWE), where the prior pulls the posterior away from the data, and the prior volume effect (PVE), where marginalization over many nuisance parameters shifts the posterior due to large-volume regions in parameter space. To mitigate these effects, we follow the conservative strategy of DESI Collaboration (2025e), Maus et al. (2024), applying Gaussian priors on nuisance parameters. We reparameterize galaxy bias in terms of (1 + b₁)σ₈(z), b₂σ₈²,  b_sσ₈², reducing projection effects and improving robustness at low σ₈. For stochastic terms (sn_{i = 0, 2, 4}) and counterterms (α_{i = 0, 2, 4}), we use Gaussian priors motivated by physical arguments. Shot-noise terms are scaled from Poisson expectations and effective velocity dispersion, while counterterms are centered at zero with width 12.5 h⁻² Mpc², corresponding to 50% of the correction at k_max = 0.20 h Mpc⁻¹. The priors are given in Table 4.

4.2. Shapefit

In this work, we use the Shape-Fit approach, which builds on the standard fit method by taking a linear power spectrum template evaluated at a reference cosmology with a set of parameters. ShapeFit captures early-time physics through the slope of the power spectrum, as highlighted by Brieden et al. (2021). They proposed an extension to the standard method by modeling the slope of the linear power spectrum using a logarithmic rescaling, modifying the reference template as

$$\begin{array}{c}\hfill ln\left(\frac{P{\prime }_{\mathrm{ref}}(k)}{P_{\mathrm{ref}}(k)}\right)=\frac{m}{a}tanh[aln\left(\frac{k}{k_p}\right)]+nln\left(\frac{k}{k_p}\right),\end{array}$$(14)

where a = 0.6 and k_p = 0.03 h Mpc⁻¹ are fixed values following Brieden et al. (2021). The parameters m and n correspond to the scale-dependent and scale-independent slopes, respectively. They encode the slope of the power spectrum (“Shape”). This ansatz was designed to replicate the effects of varying ω_b, ω_m, and n_s on the power spectrum’s shape. In this work, we vary only m, keeping n = 0 fixed, as both show a strong anti-correlation (DESI Collaboration 2025e).

In addition to the free parameters of the LPT EFT model presented in Eq. (13), we include dm = m − m_fid as an additional free parameter in our analysis. This parameter accounts for deviations from the fiducial value of the parameter m.

4.3. The distance–redshift relationship

In the context of full shape analysis, we want to model the anisotropies of the two-point measurements. To do so, we define the multipoles of the power-spectrum P_ℓ(k) by integrating the model power spectrum (Eq. (13)) over μ_k weighted by the Legendre 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 )d\mu .\end{array}$$(15)

The correlation function multipoles ξ_ℓ(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)dk,\end{array}$$(16)

where j_ℓ are the spherical Bessel functions.

The fiducial cosmology used to convert angular positions and redshifts into distances might not match the true underlying cosmology of the observational or simulated data. Such discrepancies result in additional anisotropies often referred to as Alcock-Paczynski (AP) distortions (Alcock & Paczynski 1979). We introduce the scaling parameters (q_⊥, q_∥), which adjust the observed radial and transverse wavenumbers (k_∥, k_⊥) to the true wavenumbers (k′_∥, k′_⊥) = (k_∥/q_∥, k_⊥/q_⊥). Ballinger et al. (1996) gave the transformation that maps the space $(k,\mu )=(\sqrt{k_{\Vert }^2+k_{\perp }^2},k_{\Vert }/k)$ to the space (k′,μ′) :

$$\begin{array}{cc}& k\prime =\frac{k}{q_{\perp }}{[1+{\mu }^2(\frac{q_{\perp }^2}{q_{\Vert }^2}-1)]}^{1/2},\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mu \prime =\frac{\mu q_{\perp }}{q_{\Vert }}{[1+{\mu }^2(\frac{q_{\perp }^2}{q_{\Vert }^2}-1)]}^{-1/2}.\hfill \end{array}$$(18)

Then we include the AP effect in the power spectrum multipoles as

$$\begin{array}{c}\hfill P_{\ell }(k)=\frac{2\ell +1}{2q_{\Vert }{q}_{\perp }^2}{\displaystyle {\int }_{-1}^1}d\mu P_g(k\prime (\mu ),\mu \prime (\mu )){\mathcal{L}}_{\ell }(\mu ).\end{array}$$(19)

The scaling parameters, which might indicate a shift of the BAO feature, are defined as

$$\begin{array}{cc}& q_{\Vert }=\frac{D_{\mathrm{H}}(z_{\mathrm{eff}})}{D_{\mathrm{H}}^{\mathrm{fid}}(z_{\mathrm{eff}})},\hfill \end{array}$$(20)

$$\begin{array}{cc}& q_{\perp }=\frac{D_M(z_{\mathrm{eff}})}{D_M^{\mathrm{fid}}(z_{\mathrm{eff}})},\hfill \end{array}$$(21)

where the “fid” superscript stands for the values that correspond to fiducial cosmology, $D_{\mathrm{H}}^{\mathrm{fid}}(z_{\mathrm{eff}})$ the Hubble distance, and D_M^fid(z_eff) the comoving angular diameter distance given at the effective redshift of the galaxy sample z_eff.

5. Results on mock catalogs

In this section we use the average of the 25 LRG cubicbox mocks from ABACUSSUMMITN-body simulations and the 1000 eBOSS EZMOCK LRG realizations (described in Sect. 2.2). We use these samples to test our methodology.

5.1. Fits to N-body cubic-box LRG mocks

The AbacusSummit suite of N-body simulations have realistic clustering so we use them to test the performance of our models, quantifying their biases and determining a valid range of scales. The halos are populated with luminous red galaxies using an HOD model that matches the clustering measured on early DESI data.

Our data vectors are the multipoles of the correlation function and power spectrum, averaged over the 25 realisations, reducing statistical noise and providing a more robust test of the theoretical model. The covariance matrix is obtained by performing the same measurements on a set of 1000 cubic-box EZMOCKS, which aim at reproducing the same LRG samples of the N-body simulations. We scale this covariance matrix to match the volume corresponding to 25 times 8 (Mpc/h)³. We used the Bayesian approach for our inferences (see Sect. 3). The best-fit parameter values correspond to the median of the posterior distribution and the uncertainties are derived from the central 68% confidence level.

We start by analysing the scales used in our fits. In configuration space (CS), we use separations between r_min and r_max, while in Fourier space (FS), we use wavenumbers between k_min and k_max. To test the performance of our nonlinear modeling, we increasingly fit for smaller scales: we vary r_min in CS and k_max in FS. The large-scale limit is fixed to r_max = 150 h⁻¹ Mpc in CS and k_min = 0.02 h Mpc⁻¹ in FS.

Figure 1 presents the results of RSD measurements as a function of r_min in CS. The first four rows show the best-fit values for q_∥, q_⊥, dm and df = f/f_fid, where f_fid = 0.838 is the growth rate of the simulation, computed from their corresponding fiducial parameters (Table 1). These best-fit values are compared against their expected values: q_∥^exp, q_⊥^exp = 1, df^exp = 1 while dm is expecting to be 0. Both q_∥ and q_⊥ present residuals below 1 percent for all r_min values, while df shows deviations below 3 percent, except for the smallest r_min where the deviation slightly exceeds 3 percent.

Fig. 1. Configuration space RSD results from fits to the average clustering of 25 LRG cubic-box N-body simulations as a function of the minimum separation scale, rmin. Best-fit values of q∥, q⊥, df, and dm are compared to their expected values in the top panels. The mid panel shows the reduced chi-squared rχ2, and the bottom three panels display estimated uncertainties σq∥, σq⊥, σdf, and σdm. Filled markers represent the baseline fit while empty markers indicate fits that include the additional ShapeFit parameter m.

The mid panel of Fig. 1 displays the reduced chi-squared, _rχ² ≡ χ²/(N_bins − N_par). As expected, the _rχ² values increase for decreasing r_min. Scales below 30 h⁻¹ Mpc exhibit _rχ² > 3, indicating a poor fit and suggesting that the model breaks down on such scales. As r_min increases to approximately 40 h⁻¹ Mpc, the _rχ² stabilizes around 1, reflecting improved agreement between the model and the data.

The last four rows of Fig. 1 display the uncertainties σ_q∥, σ_q⊥, σ_df and σ_dm of the parameters of interest. Uncertainties naturally become larger with increasing r_min. Consistent with past works, all values of σ_q∥ are approximately twice as large as those for σ_q⊥.

Our baseline choice for CS analysis aims at maximizing the information extracted from the correlation function while obtaining unbiased results on cosmological parameters and a good quality of fit. Based on the results just discussed, we set our baseline choice for our CS analysis to r_min = 30 h⁻¹ Mpc. It is worth noting that we also tested the inclusion of stochastic terms, similar to those used in Fourier space, within the configuration space analysis. This addition helps regularize small scales beyond 30 h⁻¹ Mpc, reducing the reduced χ² for r_min = 25 h⁻¹ Mpc to 1.6.

Figure 2 presents the similar analysis as Fig. 1 but for Fourier space (FS) and both joint space analyses (JS and JS_sep). Here, we are interested in the impact of varying the maximum wavenumber k_max.

Fig. 2. Same as Fig. 1 but for Fourier space (FS, red circles) and the joint space analyses (JS and JSsep, gray circles and green squares, respectively). For reference, CS results for rmin = 30 h−1 Mpc are shown in blue.

Focusing first on FS results in Fig. 2 (red points), both q_∥ and q_⊥ show deviations below 1 percent and df below 3 percent for all k_max. The reduced chi-squared _rχ² remains stable, slightly below 1.4, comparable to CS results at r_min = 30 h⁻¹ Mpc. The uncertainties on σ_q⊥ and σ_f are notably larger than in CS, expect for k_max ≥ 0.22 h Mpc⁻¹ where they are comparable. Note that FS includes three additional parameters (shot-noise terms), which reduce its constraining power relative to the CS analysis. In the following, we restrict our analysis to k_max = 0.20 h Mpc⁻¹, a choice that is justified later when discussing the joint results. Figure 2 also presents the joint-space results, which is discussed later.

Focusing now in the JS analysis in Fig. 2 (gray circles), we observe small shifts on the best-fit values of all three parameters compared to FS or CS, though not significant. Regarding estimated uncertainties, JS yields uncertainties comparable to, or smaller than, those in CS, indicating some small gain in the combination with FS. However, _rχ² values are significantly higher for JS, with values above 2.5 and increasing with k_max, indicating a poor fit at this statistical precision level. As we discuss next, the quality of the JS fit is degraded due to tensions in some of the nuisance parameters between CS and FS. For this reason, we perform a modified joint analysis where parameters b₂, b_s, α₀, α₂, and α₄ are allowed to vary independently in the power spectrum and correlation function models. This alternative approach is referred to as JS_sep. Figure 2 presents the JS_sep results (gray squares), as a function of k_max. The overall performance of the fit is comparable to the standard JS analysis, though we note significant a reduction in _rχ² to ∼1.5, comparable to CS and FS fits. For k_max > 0.20 h Mpc⁻¹, the JS_sep results for f and dm exhibit some discrepancies. Given those results, we set the baseline choice of k_max = 0.20 h Mpc⁻¹ for FS, JS and JS_sep.

Figure 3 presents the best-fit models with our baseline choice of separations (30 < r < 150 h⁻¹ Mpc and 0.02 < k < 0.20 h Mpc⁻¹) of our two-point multipoles (ℓ = 0, 2, 4), along with their respective residuals in the lower panels. The differences between the models and the measurements remain within a 2σ deviation, except for the hexadecapole in CS. We remind that neighboring separations in CS measurements are strongly correlated (Fig. A.1) so these deviations are not significant.

Fig. 3. Fourier space (FS) power spectrum multipoles (left panel) and correlation function (CS) multipoles (right panel) with their best-fit models. We consider the monopole (blue), quadrupole (orange), and hexadecapole (green) in both spaces. The three bottom panels show the normalized residuals and the shaded areas correspond to 2σ. The best-fit model of FS and CS are shown in solid lines while the joint-space JSsep model is represented with a dashed line.

Figure 4 shows the 68 and 95 percent contours of all parameter posteriors for FS, CS and JS. Posteriors in both FS and CS are in good agreement for all parameters, though CS constraints tend to be tighter (e.g., b₂, b_s, α₀). For JS, contours for q_∥ and q_⊥ exhibit slightly larger biases compared to individual analyses, although they remain within 2σ. Surprisingly, the JS constraints on b₁, b₂, α₂ and α₄ are inconsistent with the FS and CS results. This might be due to the breaking of degeneracies in a multi-dimensional parameter space. An important caveat is that our JS fits use a sample covariance from 1000 mocks, though the data vector is larger. This covariance is thus noisier for JS compared to FS or CS (see Table 2), which might shift contours. Given the large cross-covariance between FS and CS (as illustrated in Fig. A.1), posteriors might be more sensitive to sample noise. To test this hypothesis, we computed JS posterior while neglecting the cross-covariance between FS and CS multipoles. These posteriors are shown in dashed lines in Fig. 4, now in good agreement with the FS-CS contours, as expected when combining two independent likelihoods. We also tested fits on the average of 1000 EZMOCKS (both DESI and eBOSS versions), using covariance matrices derived from the same set of mocks. This guarantees a correct covariance treatment in the fit, as it is directly derived from the same mocks used for the average. For both EZMOCKS sets, no significant bias is observed between CS, FS, and JS analyses. All three posteriors are fully consistent, in contrast to the bias seen in Fig. 4 for JS when using cubic box mocks. We then compared the results of Fig. 4 with the best-fit values obtained using IMINUIT. For all analyses, the IMINUIT best-fit points lie at the center of the contours, except for b_s and b₂ which can be due to the presence of local maximum in the likelihood. This supports the bias observed for b₁, α₂, and α₄. We also test the impact of the finite number of mocks used to compute the covariance by using 500 instead of 1000 mocks. No significant difference was seen, which discards the impact of the limited number of mocks. These tests could indicate that the cross terms in a covariance matrix derived from 1000 EZMOCKS are not accurate enough when fitting N-body simulations.

Fig. 4. Likelihood posterior (1σ − 2σ) for FS (red), CS (blue), and JS (gray) obtained from fits to the average of the 25 LRG CubicBox mocks. JS contours without the cross covariance (labelled “w/o crosscov”) are displayed in dashed gray. The dotted black lines represent the expected values. Results are shown for kmax = 0.20 h Mpc−1 and rmin = 30 h−1 Mpc.

Figure 5 presents posteriors for FS, CS (same as in Fig. 4) as well as those from the JS_sep analysis and the Gaussian approximation method (GA; see Sect. 3.2). All four contours are consistent with each other within 1σ. The GA contours tend to fall between the FS and CS contours, while the JS_sep do not show significant bias in the four parameters of interest to the separate analyses.

Fig. 5. Contours of 68 and 95 percent of the posterior distributions for FS (red), CS (blue), JSsep (green), and GA (orange) derived from fits to the average clustering of 25 LRG cubic-box mocks from N-body simulations.

Additionally, JS_sep does not show any notable improvement over GA when tested on the average of the 25 mocks, as both yield similar results.

These tests validate our baseline choices for the analyses and we turn to the study of a large set of mock catalogs.

5.2. Fits to 1000 eBOSS LRG EZmocks

We evaluated the statistical properties of our best-fit values and their uncertainties with fits to 1000 realizations of eBOSS LRG EZMOCKS, which reproduce realistic angular and redshift distributions (see Sect. 2). These mocks allow us to assess how the results on individual realizations compare to their corresponding distributions, and check that our uncertainties are correctly estimated. We computed two-point multipoles for all realizations and derived their corresponding sample covariance matrix (Eq. (2)).

We used the frequentist approach to fit individual mocks (see Sect. 3) as it is faster than the Bayesian approach.

Figure 6 compares the distributions of best-fit parameters in FS (y-axis) and CS (x-axis) for the 1000 individual eBOSS EZMOCKS. The top row displays the best-fit values for q_∥, q_⊥, and f, while the bottom row shows their estimated errors σ_q∥, σ_q⊥, and σ_f. The density of points is represented by grayscale shading, and the diagonal dashed line is the identity line. The distribution of the best-fit values is well correlated and mostly scattered around the identity line, except for q_∥, which tends to have slightly higher values in FS than in CS.

Fig. 6. Comparison of the RSD measurements in FS (y-axis) and CS (x-axis) for the individual 1000 EZMOCKS of the eBOSS LRG sample. The top row displays the best-fit values for q∥, q⊥, and f, while the bottom row shows their estimated errors σq∥, σq⊥, and σf. The density of points is represented by grayscale shading, with dashed lines indicating the line of perfect correlation where the compared parameters are equal. We display the Pearson correlation coefficients between CS and FS. The plus symbol represents the results from the real eBOSS data, as listed in Table 6.

The discrepancy seen for q_∥ in FS may be due to non-Gaussian effects in the posteriors of parameters, which could be caused by the volume difference between the large effective volume of the covariance and the smaller volume of individual mocks. To mitigate this effect, we combine the 1000 EZMOCKS in groups of 5, making a total of 200 independent EZMOCKS with a volume 5 times larger. The discrepancy in FS is still present, which discards this hypothesis. We simply think that this bias is due to the approximate nature of EZMOCKS and that it is not guaranteed that the recovered cosmological parameters are aligned with the expected value, unlike ABACUS cubic boxes. We also performed a BAO fit on the average of the 1000 EZMOCKS, and both CS and FS recovered the expected value, while the full-shape fit on this average still shows the bias. This may indicate that when the fit is limited to higher scales as in the BAO case, the results are not biased, which confirms that the issue is related to the approximations in EZMOCKS. The uncertainties show good agreement overall for q_∥ and q_⊥, but for f, CS uncertainties are larger than FS on average.

We computed consensus results using the Gaussian approximation method (GA) for each mock realization. First, we construct the 6 × 6 parameter covariance matrix from the 1000 eBOSS EZMOCKS best-fit values of q_∥, q_⊥ and f. Figure A.2 displays the resulting correlation matrix, where we can notice in particular that the correlations for q_∥, q_⊥ and fσ₈ are 61, 51 and 81 percent, respectively, between CS and FS. Second, we computed consensus parameters and their corresponding covariances for each mock realization, adjusting each time the off-diagonal blocks of the parameter covariance (see Sect. 3.2).

Figure 7 shows the distributions of q_∥, q_⊥, f, and their associated errors σ_q∥, σ_q⊥, σ_f, as estimated with the GA method, compared to the same quantities obtained from CS and FS. The distribution of the best-fit values is well correlated and generally scattered around the identity line, except for q_∥, where the scatter between GA and CS is slightly smaller than that between GA and FS. This suggests a higher weight of CS in the GA analysis compared to FS. Additionally, the errors estimated by GA are systematically smaller than those obtained in either CS or FS analyses. It is worth noting that the distribution of σ_f splits into two distinct clouds, as shown in Fig. 6.

Fig. 7. Comparison of the RSD measurements in FS (red) and CS (blue) with GA on the individual 1000 eBOSS EZmock. FS and CS results are shown on the y-axis, while GA results are on the x-axis. The top row displays the best-fit values for q∥, q⊥, and fσ8, while the bottom row shows their estimated errors σq∥, σq⊥, and σfσ8. The density of points is represented by the color-scale shading, and the dashed lines indicate the line of perfect correlation where the compared parameters are equal. We display the Pearson correlation coefficients between GA and FS, as well as GA and CS. The plus symbols represent the best-fit results from the real eBOSS data as listed in Table 6.

Figure 8 is analogous to Fig. 7 but for the results from the JS_sep method. Similar to the GA results, the best-fit values are well-correlated and closely cluster around the identity line. However, the scatter is noticeably more pronounced in JS_sep compared to GA.

Fig. 8. Same as Fig. 7 but for the JSsep method.

While JS_sep uncertainties are consistently smaller than or similar to those observed in CS and FS, they exhibit greater scatter relative to the GA results shown in Fig. 7. This scatter is due to the different nature of the fit, which now incorporates the cross-covariance between FS and CS at the level of the two-point functions. These noisy terms (due to a limited number of mocks) are propagated to the final infered parameters. The similar or smaller uncertainties show that there is gain in combining FS and CS results at the two-point level, even though they are highly correlated.

Table 5 summarizes the statistical properties of (q_∥, q_⊥, f) for CS, FS, GA, and JS_sep. For each parameter p, we show the average bias with respect to its expected value, Δp = ⟨p_i − p_exp⟩, the standard deviation of best-fit values σ_p, the mean estimated error ⟨σ_pi⟩, the mean of the pull $Z_{p_i}\equiv \frac{p_i-⟨p_i⟩}{{\sigma }_{p_i}}$, and its standard deviation. If the parameter uncertainties are correctly estimated and follow a Gaussian distribution, we expect σ_p = ⟨σ_pi⟩, ⟨Z_pi⟩ = 0 and σ_Zpi = 1. The table also includes the mean value of the reduced chi-square, _rχ_min². From this table, we observe:

• the number of valid realizations, N_good, is nearly maximal for all analyses, indicating that most fits have successfully converged. The discarded cases in FS and GA are due to the removal of outliers beyond the 5σ threshold. For each method, the average minimum reduced chi-square ⟨_rχ_min²⟩ is close to 1, suggesting that the majority of the mocks are well-modeled, with JS_sep exhibiting the closest value from 1.

• the mean estimated error, ⟨σ_pi⟩, is generally consistent with the standard deviation of the best-fit values, σ_p, across all parameters and analyses, indicating no significant issues with the uncertainty estimations. Notably, the mean estimated errors for JS_sep are smaller than those for CS, FS, and GA (except for q_⊥) which means that there is small gain in performing the joint analysis.

• the standard deviations of the best-fit values, σ_p, are consistently smaller in GA compared to JS_sep for all three parameters. It is important to note that these standard deviations are scaled by the correction factor $\sqrt{m_2}$ (Sect. 3.1) and these corrections are different between FS, CS and JS.

• the average bias for q_∥ is larger in JS_sep and FS (more than 3 percent) when compared to GA and CS. The bias in JS_sep might be pulled by the large bias seen in FS. The GA results lie between FS and CS as expected, for all parameters. For the growth rate f, JS_sep has the smallest bias among all methods.

• the mean of the pull ⟨Z_xi⟩ is consistent with zero for all methods. For f, the smallest value of the mean pull is achieved by JS_sep. The standard deviations of the q_⊥ and q_∥ pulls are slightly larger than 1 in all cases, suggesting that errors are slightly underestimated, with the most significant underestimation for q_⊥ in JS_sep, of about 30 percent. For f, the standard deviation of JS_sep is slightly smaller than 1, indicating a slight overestimation of the error. However, JS_sep provides results closest to 1 compared to the other methods.

• The high value of σ_Zq⊥ observed in JS_sep may be due to the limited number of mocks used for the covariance. We recomputed the same statistics presented in Table 5 using a covariance matrix derived from twice fewer mocks, i.e., 500 mocks. The standard deviation of the pull for the three parameters in JS_sep increased to ∼4, which may indicate the necessity of using a larger number of mocks for the covariance when performing JS analyses to obtain more reliable uncertainty estimates.

6. Application to the BOSS+eBOSS LRG sample

In this section we apply our methodology to the BOSS+eBOSS LRG sample presented in Sect. 2.1. We present results for Fourier space (FS), configuration space (CS), and our two consensus methods, the Gaussian approximation (GA) and joint space (JS_sep).

Figure 9 displays the best-fit models for the power spectrum and correlation function multipoles for the separate analyses (FS and CS, shown as solid lines) and the joint analysis (JS_sep, shown as dashed lines). The residuals are shown in the bottom panels, highlighting good agreement between JS_sep and CS/FS, which we attribute to differences in the inferred bias and counterterms.

Fig. 9. Best-fit models for the monopole, quadrupole, and hexadecapole of BOSS+eBOSS LRG sample in Fourier-space (left panel, solid line) and configuration space (right panel, solid line). The bottom three subpanels display the normalized residuals, with shaded areas indicating a 2σ range. The best-fit model for joint space fit, JSsep, is shown as a dashed line.

Figure 10 shows the contours of the posterior distribution (1σ and 2σ) for the FS, CS, GA, and JS_sep analyses, all of which are in good agreement. The posteriors are consistent with one another (within 1σ), displaying comparable widths except for f in CS, which shows a higher value with a significantly larger width, though consistent with the behavior seen in eBOSS mocks (Fig. 6), and the posterior of m for GA which is completly biased by more than 2σ) The JS_sep posterior is consistent with the other posteriors and comparable to GA, though we can see how these contours are not exactly Gaussian, which is one of the advantages of the method.

Fig. 10. Likelihood posterior for FS (red), CS (blue), JSsep (green), and GA (orange) derived from fits on the eBOSS LRG data.

Figure 10 shows the contours of the posterior distribution (1σ and 2σ) for the FS, CS, GA, and JS_sep analyses, all of which are in good agreement. The posteriors are consistent with one another (within 1σ), displaying comparable widths except for f in CS, which shows a higher value with a significantly larger uncertainty. This behavior is consistent with what we observe in the eBOSS mocks (Fig. 6). The posterior of m in GA, however, is significantly biased, deviating by more than 2σ, illustrating the limitations of the GA method when the posteriors are non-Gaussian, as is the case for this parameter in CS. The JS_sep posterior remains consistent with the other estimates and closely tracks GA, except for m, although its contours exhibit non-Gaussian features–highlighting one of the key strengths of the joint-space approach.

Table 6 summarizes the results of the four different analyses. The reduced χ² values are 1.82, 1.10, and 1.34 for FS, CS, and JS_sep, respectively. The high χ² value of FS represents a deviation greater than 3σ compared to the distribution of χ² from EZMOCK eBOSS. The best-fit values for the parameters q_∥, q_⊥, and fσ₈ are all consistent within their respective 1σ errors. For q_∥ and q_⊥, JS_sep provides the smallest uncertainties while yielding the highest and lowest best-fit values, respectively, compared to the other methods. Additionally, JS_sep yields a slightly smaller value for fσ₈, similar to the one from CS. Still in Table 6, we compare our main JS_sep results with the official RSD results from Gil-Marín et al. (2020), FS eBOSS), Bautista et al. (2021), CS eBOSS), and the consensus between the two (GA eBOSS). The main difference between our analysis and official ones lies in the modeling approach. Previous analyses employed the TNS model (Taruya et al. 2010), which is based on an Eulerian description and incorporates a damping term for fingers-of-God rather than counterterms. The scale ranges are the same for both CS analyses while for FS, we increased from k_max = 0.15 to 0.20 h Mpc⁻¹.

Our final result on this dataset using our JS_sep method yields

$$\begin{array}{c}\hfill \{\begin{array}{c}{q}_{\Vert }=1.031\pm 0.044\hfill \\ q_{\perp }=0.995\pm 0.021\hfill \\ f{\sigma }_8=0.463\pm 0.052\hfill \end{array}.\end{array}$$(22)

Our results are compatible with all three analyses but show a slightly lower value for q_∥, a higher value for q_⊥, and for fσ₈. These discrepancies can largely be attributed to differences in the clustering models, the scale ranges, the inclusion of cross-covariance terms in the JS analysis, and the inclusion of systematic errors (in the case of the official ones). Uncertainties on q_∥ are notably larger than in the official eBOSS results, likely due to the inclusion of additional counterterms and shot noise terms in our model, which introduce μ² and μ⁴ dependencies (see Eq. (13)).

7. Conclusion

This work presents a new RSD analysis of the DR16 LRG sample (BOSS+eBOSS), where we extract cosmological parameters simultaneously from the configuration-space and Fourier-space two-point multipoles. We compared the joint-space (JS) approach to obtain consensus results with the commonly used Gaussian approximation (GA), which combines CS and FS results at the parameter level, assuming Gaussian posteriors.

To validate our methodology, we first applied it to the average of 25 N-body LRG cubic-box mocks from the AbacusSummit simulation. We assessed parameter consistency between CS and FS and compared JS results with each individual space. Unexpected behavior was observed for parameters such as b₂, α₂, and α₄, with 1σ–2σ deviations from CS and FS, possibly due to noisy cross-terms in the covariance matrix. Allowing parameters b₂, b_s, and the counterterms to vary independently in JS (JS_sep) led to unbiased results with a reduced _rχ² = 1.5. However, JS did not show significant improvements over CS, FS, or GA when tested on the mock average. Further tests were conducted on 1000 LRG EZmocks, which reproduce the DR16 LRG sample geometry. We found JS results to be consistently smaller than those of CS and FS but with greater scatter compared to GA.

After validation of our framework, we applied it to real data of the DR16 LRG multipoles. The JS method provided the tightest constraints on q_∥, q_⊥, and fσ₈, remaining consistent with the official 2020 results from Bautista et al. (2021) and Gil-Marín et al. (2020). However, JS_sep yielded the highest value for fσ₈ = 0.463 ± 0.052, which we attribute to differences in modeling and the inclusion of cross-covariance terms between FS and CS.

We believe our framework is more robust than FS, CS, or GA approaches and should be used in future surveys. Systematic effects in either space are averaged out at the multipole level, providing unbiased constraints on the growth rate and dilation parameters. Future work employing our methods will pay attention to the precision of covariance matrices estimates, as the large noise might affect the final constraints.

Acknowledgments

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). Thanks to NERSC.

References

Appendix A: Covariance matrix

Figure A.1 shows the correlation matrix R, defined as

$$\begin{array}{c}\hfill \mathit{R}={\mathit{\sigma }}^{-1}\widehat{\mathit{C}}{[{\mathit{\sigma }}^{-1}]}^T,\end{array}$$(A.1)

where σ is derived from the square root of the diagonal elements of the covariance matrix $\widehat{\mathit{C}}$. Beyond the well-known correlations between multipoles in the same space, the structure of the off-diagonal blocks reveals the correlation dynamics between the two distinct spaces. The multipole-multipole correlations between these spaces exhibit similar patterns. This arises due to the Fourier transform, which links a specific physical scale r with multiple spectral modes k. At a given scale r, there are notable correlations (and anti-correlations) with several k modes. This indicates that statistical information from a physical scale is distributed across multiple spectral modes rather than being confined to a single mode. These oscillations seen in the covariance can be derived analytically following Grieb et al. (2016) and shown in Appendix C of Dumerchat & Bautista (2022).

Fig. A.1. Normalized covariance matrix of the power spectrum multipoles P0, 2, 4 and the correlation function multipoles ξ0, 2, 4, measured from the 1000 cubic-box EZMOCKS of LRGs.

Figure A.2 shows the correlation matrix of the parameters q_∥, q_⊥ and fσ₈, obtained from 1000 fits to the eBOSS LRG EZmocks. This matrix is used to obtain the Gaussian approximation consensus results in Sect. 5.

Fig. A.2. Correlation coefficients between q⊥, q∥, and fσ8 for CS and FS, derived from fits to the 1000 EZMOCKS of the eBOSS LRG sample.

Appendix B: Clustering model

We use a Lagrangian perturbation theory (LPT) model to describe the observed clustering in both configuration and Fourier space. This model is based on the Lagrangian formalism developed in Chen et al. (2020, 2021). In this section we provide a brief overview of this formalism.

B.1. Lagrangian perturbation theory

In the Lagrangian picture, cosmological structure formation is modeled by following the trajectories of fluid elements rather than studying the dynamics of density and velocity fields. The variable of interest in LPT is the displacement field Ψ(q, η), which corresponds to the difference between the Lagrangian coordinates q (at some initial time η₀) and the Eulerian coordinates x(q) of a fluid element at a given conformal time η:

$$\begin{array}{c}\hfill \mathbf{x}(\mathit{q},\eta )=\mathit{q}+\mathbf{\Psi }(\mathit{q},\eta ),\end{array}$$(B.1)

where Ψ(q, η₀) = 0. Each element of the fluid is labeled by q, and Ψ(q, η) fully describes its evolution. The Eulerian density field ρ(x, η) satisfies the continuity relation:

$$\begin{array}{c}\hfill \rho (\mathbf{x},\eta )d^{3}x=\overline{\rho }(\eta )d^{3}q,\end{array}$$(B.2)

where $\overline{\rho }(\eta )$ is the mean density in comoving coordinates.

We can define the Jacobian of the transformation between Eulerian and Lagrangian coordinates as $J=det\left[\frac{\partial \mathbf{x}}{\partial \mathit{q}}\right]$:

$$\begin{array}{c}\hfill d^{3}x=J(\mathit{q},\eta )d^{3}q.\end{array}$$(B.3)

Using the definition of the density contrast $\rho (\mathbf{x},\eta )=\overline{\rho }(\eta )[1+\delta (\mathbf{x},\eta )]$ and Eqs. B.2 and B.3, we can relate the density contrast to the displacement field:

$$\begin{array}{c}\hfill \delta (\mathbf{x},\eta )=\frac{1}{J(\mathit{q},\eta )}-1.\end{array}$$(B.4)

This relation is equivalent to the following equation:

$$\begin{array}{c}\hfill 1+\delta (\mathbf{x},\eta )={\displaystyle \int d^3}q{\delta }_D[\mathbf{x}-\mathit{q}-\mathrm{\Psi }(\mathit{q},\eta )].\end{array}$$(B.5)

Fourier transforming Eq. B.5 (in x-space) gives

$$\begin{array}{c}\hfill \delta (\mathit{k},\eta )={\displaystyle \int d^3}q{e}^{i\mathit{k}·\mathit{q}}(e^{i\mathit{k}·\mathrm{\Psi }(\mathit{q},\eta )}-1).\end{array}$$(B.6)

The equation of motion for particle trajectories x(η) is

$$\begin{array}{c}\hfill \ddot{\mathbf{x}}+\mathcal{H}\dot{\mathbf{x}}=-{\mathrm{\nabla }}_{\mathbf{x}}\mathrm{\Phi }(\mathbf{x}).\end{array}$$(B.7)

From which we can infer the evolution equation of Ψ using Eq. B.1:

$$\begin{array}{c}\hfill \ddot{\mathbf{\Psi }}+\mathcal{H}\dot{\mathbf{\Psi }}=-{\mathrm{\nabla }}_{\mathbf{x}}\mathrm{\Phi }(\mathbf{x}),\end{array}$$(B.8)

where overdots indicate derivatives w.r.t. conformal time η. This equation is solved order by order in a perturbative way such that:

$$\begin{array}{c}\hfill \mathbf{\Psi }={\mathbf{\Psi }}^{\mathbf{(}\mathbf{1}\mathbf{)}}+{\mathbf{\Psi }}^{\mathbf{(}\mathbf{2}\mathbf{)}}+{\mathbf{\Psi }}^{\mathbf{(}\mathbf{3}\mathbf{)}}+...\end{array}$$(B.9)

The first-order solution of LPT Ψ⁽¹⁾ is known as the Zel’dovich approximation (Zel’dovich 1970),

$$\begin{array}{c}\hfill {\mathbf{\nabla }}_{\mathit{q}}{\mathbf{\Psi }}^{(1)}=-D_{+}(\eta ){\delta }_0(\mathit{q}),\end{array}$$(B.10)

or equivalently,

$$\begin{array}{c}\hfill \mathbf{\Psi }(\mathit{q},\eta )=i{D}_{+}(\eta ){\displaystyle \int \frac{d^{3}k}{{(2\pi )}^3}}\frac{\mathit{k}}{k^2}{\delta }_0(\mathit{k})e^{i\mathit{k}·\mathit{q}},\end{array}$$(B.11)

where D₊(η) is the linear growth factor and δ₀(q) is the initial density contrast. The Zel’dovich approximation uses the linear displacement field to model cosmic structure formation. One can find more details on the Zel’dovich approximation in White (2014). Higher-order solutions are expressed as integrals of higher powers of the linear density field (Matsubara 2015).

Unlike in SPT, there is no simple recursive solution for order-by-order expression in LPT due to the fact that the displacement fields in LPT are not irrotational in general (Rampf 2012; Matsubara 2015).

B.2. Effective field theory of large-scale structures

In this work we use the Effective Field Theory of large-scales structures (EFTofLSS). The main insight is that the LSS has various hierarchical distance scales. In the large scale limit, only the large-scale degrees of freedom are required, while the effects of unknown short-scale physics can be systematically captured by effective operators. This decoupling framework forms the core of the EFTofLSS (Baumann et al. 2012; Carrasco et al. 2012; Cabass et al. 2022; Ivanov 2022). The Lagrangian form of EFT involves smoothing Eq. B.8 and incorporating small-scale physics through a series of counterterms.

To achieve this, we used a filter W_R(q, q′), which separates the system into long-wavelength (L) and short-wavelength (S) modes:

$$\begin{array}{cc}& {\mathbf{\Psi }}_L(\mathit{q})={\displaystyle \int d^3}q\prime W_R(\mathit{q},\mathit{q}\mathbf{\prime })\mathbf{\Psi }(\mathit{q}\mathbf{\prime })\hfill \end{array}$$(B.12)

$$\begin{array}{cc}& {\mathbf{\Psi }}_S(\mathit{q},\mathit{q}\mathbf{\prime })=\mathbf{\Psi }(\mathit{q}\mathbf{\prime })-{\mathbf{\Psi }}_L(\mathit{q}).\hfill \end{array}$$(B.13)

We can write the displacement as

$$\begin{array}{c}\hfill \mathbf{\Psi }(\mathit{q})={\mathbf{\Psi }}_L^{(1)}(\mathit{q})+{\mathbf{\Psi }}_L^{(2)}(\mathit{q})+{\mathbf{\Psi }}_L^{(3)}(\mathit{q})+\cdots +\frac{1}{2}{\alpha }_1\mathbf{\nabla }{\delta }_0+\mathit{S}+\cdots ,\end{array}$$(B.14)

where the first three terms come from the perturbative treatment of the long-wavelength evolution and the last two terms parameterize the impact of the short-wavelength modes on the evolution where α₁ is a counterterm coefficient and S is equivalent to a stochastic term

B.3. Redshift-space distortions

Before writing the power spectrum we need to include bias scheme. The density of biased tracers can be modeled, assuming Lagrangian bias, by multiplying the δ_D in the above equation by the bias functional $F[{\delta }_0(\mathit{q}),{\mathbf{\nabla }}^2{\delta }_0(\mathit{q}),...]$ depending on linear density and its derivative. The bias functionnal is described by Eq. 4.4 in Chen et al. (2020), more details about the modeling of biased tracers in the Lagrangian picture can be found in Matsubara (2008a), McDonald & Roy (2009) and Vlah et al. (2016). Equation B.6 simply becomes for biased tracers:

$$\begin{array}{c}\hfill \delta (\mathit{k},\eta )={\displaystyle \int d^3}q{e}^{i\mathit{k}·\mathit{q}}(e^{i\mathit{k}·\mathrm{\Psi }(\mathit{q},\eta )}F(\mathit{q})-1).\end{array}$$(B.15)

The power spectrum P(k) can be expressed by the displacement field using Eq.(B.15):

$$\begin{array}{c}\hfill P(\mathit{k})={\displaystyle \int d^3}\mathit{q}{e}^{i\mathit{k}·\mathit{q}}(\langle e^{i\mathit{k}·\mathbf{\Delta }}\rangle -1)\end{array}$$(B.16)

where we have defined the pairwise Lagrangian displacement ${\mathbf{\Delta }}_i={\mathbf{\Psi }}_i({\mathit{q}}_1)-{\mathbf{\Psi }}_i({\mathit{q}}_2)$. The cumulant theorem allows us to write the expectation value of the exponential in terms of the exponential of expectation values. The cumulant theorem is well defined in Bernardeau et al. (2002), Matsubara (2008a,b) and can be expressed as

$$\begin{array}{c}\hfill \langle e^{-iX}\rangle =exp\left[{\displaystyle \sum _{N=1}^{\infty }}\frac{{(-i)}^N}{N!}{\langle X^N\rangle }_c\right]\end{array}$$(B.17)

where 〈X^N〉_c is the cumulant of a random variable X. The bracketed average in Eq. (B.16) can be expressed as

$$\begin{array}{c}\hfill ln\langle e^{i\mathit{k}·\mathbf{\Delta }}\rangle =-\frac{1}{2}{k}_i{k}_j{A}_{\mathit{ij}}-\frac{i}{6}{k}_i{k}_j{k}_k{W}_{\mathit{ijk}}+\dots \end{array}$$(B.18)

where $A_{\mathit{ij}}={\langle {\mathrm{\Delta }}_i{\mathrm{\Delta }}_j\rangle }_c$ and $W_{\mathit{ijk}}={\langle {\mathrm{\Delta }}_i{\mathrm{\Delta }}_j{\mathrm{\Delta }}_k\rangle }_c$ are the cumulants of the pairwise displacements. More details of those terms can be found on Vlah et al. (2015), Chen et al. (2021).

The power spectrum in Eq. B.16 is valid at all times, not only on the linear regime. From Taylor & Hamilton (1996) we can express the power spectrum in the Zel’dovich approximation. In the Zel’dovich approximation, the displacement field is assumed to be Gaussian at all times. The expectation value of the exponential in Eq. B.16

$$\begin{array}{c}\hfill \langle e^{i\mathit{k}·\mathbf{\Delta }}\rangle =exp(-k_i{k}_j[{\mathrm{\Psi }}_{\mathit{ij}}(0)-{\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})]),\end{array}$$(B.19)

with q = q₁ − q₂ and

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})=\langle {\mathrm{\Psi }}_i({\mathbf{q}}_1){\mathrm{\Psi }}_j({\mathbf{q}}_2)\rangle .\end{array}$$(B.20)

Using the Zel’dovich solution Eq.B.11 we can relate the Ψ_ij(q) to the linear power spectrum (Taylor & Hamilton 1996):

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})=\frac{1}{{(2\pi )}^3}{\displaystyle \int d^3}k\frac{k_i{k}_j}{k^4}{P}_{\mathit{lin}}(k)e^{-i\mathit{k}·\mathit{q}}\end{array}$$(B.21)

This results to the Zel’dovich power spectrum:

$$\begin{array}{c}\hfill P(\mathit{k})={\displaystyle \int d^3}\mathit{q}{e}^{i\mathit{k}·\mathit{q}}\times [e^{(-k_i{k}_j[{\mathrm{\Psi }}_{\mathit{ij}}(\mathbf{0})-{\mathrm{\Psi }}_{\mathit{ij}}(\mathit{q})])}-1]\end{array}$$(B.22)

This equation describes the nonlinear transformation of the initial power spectrum into the evolved power spectrum.

Finally, in the Lagrangian formalism, the RSD can be accounted by boosting displacements along the line-of-sight (LOS) direction $\widehat{n}$ by the velocities of galaxies:

$$\begin{array}{c}\hfill \mathbf{\Psi }\to {\mathbf{\Psi }}^s=\mathbf{\Psi }+\frac{\widehat{n}(\mathit{v}·\widehat{n})}{\mathcal{H}}\end{array}$$(B.23)

Where the superscript s is referring to vectors boosted into redshift space and v is the peculiar velocity of galaxy.

The final expression of the theoretical model for the redshift-space galaxy power spectrum, accounting for bias in LPT and EFTofLSS is giving by Eq. 13.

All Tables

All Figures

Fig. 1. Configuration space RSD results from fits to the average clustering of 25 LRG cubic-box N-body simulations as a function of the minimum separation scale, rmin. Best-fit values of q∥, q⊥, df, and dm are compared to their expected values in the top panels. The mid panel shows the reduced chi-squared rχ2, and the bottom three panels display estimated uncertainties σq∥, σq⊥, σdf, and σdm. Filled markers represent the baseline fit while empty markers indicate fits that include the additional ShapeFit parameter m.

In the text

	Fig. 2. Same as Fig. 1 but for Fourier space (FS, red circles) and the joint space analyses (JS and JSsep, gray circles and green squares, respectively). For reference, CS results for rmin = 30 h−1 Mpc are shown in blue.
In the text

Fig. 3. Fourier space (FS) power spectrum multipoles (left panel) and correlation function (CS) multipoles (right panel) with their best-fit models. We consider the monopole (blue), quadrupole (orange), and hexadecapole (green) in both spaces. The three bottom panels show the normalized residuals and the shaded areas correspond to 2σ. The best-fit model of FS and CS are shown in solid lines while the joint-space JSsep model is represented with a dashed line.

In the text

	Fig. 4. Likelihood posterior (1σ − 2σ) for FS (red), CS (blue), and JS (gray) obtained from fits to the average of the 25 LRG CubicBox mocks. JS contours without the cross covariance (labelled “w/o crosscov”) are displayed in dashed gray. The dotted black lines represent the expected values. Results are shown for kmax = 0.20 h Mpc−1 and rmin = 30 h−1 Mpc.
In the text

	Fig. 5. Contours of 68 and 95 percent of the posterior distributions for FS (red), CS (blue), JSsep (green), and GA (orange) derived from fits to the average clustering of 25 LRG cubic-box mocks from N-body simulations.
In the text

Fig. 6. Comparison of the RSD measurements in FS (y-axis) and CS (x-axis) for the individual 1000 EZMOCKS of the eBOSS LRG sample. The top row displays the best-fit values for q∥, q⊥, and f, while the bottom row shows their estimated errors σq∥, σq⊥, and σf. The density of points is represented by grayscale shading, with dashed lines indicating the line of perfect correlation where the compared parameters are equal. We display the Pearson correlation coefficients between CS and FS. The plus symbol represents the results from the real eBOSS data, as listed in Table 6.

In the text

Fig. 7. Comparison of the RSD measurements in FS (red) and CS (blue) with GA on the individual 1000 eBOSS EZmock. FS and CS results are shown on the y-axis, while GA results are on the x-axis. The top row displays the best-fit values for q∥, q⊥, and fσ8, while the bottom row shows their estimated errors σq∥, σq⊥, and σfσ8. The density of points is represented by the color-scale shading, and the dashed lines indicate the line of perfect correlation where the compared parameters are equal. We display the Pearson correlation coefficients between GA and FS, as well as GA and CS. The plus symbols represent the best-fit results from the real eBOSS data as listed in Table 6.

In the text

	Fig. 8. Same as Fig. 7 but for the JSsep method.
In the text

	Fig. 9. Best-fit models for the monopole, quadrupole, and hexadecapole of BOSS+eBOSS LRG sample in Fourier-space (left panel, solid line) and configuration space (right panel, solid line). The bottom three subpanels display the normalized residuals, with shaded areas indicating a 2σ range. The best-fit model for joint space fit, JSsep, is shown as a dashed line.
In the text

	Fig. 10. Likelihood posterior for FS (red), CS (blue), JSsep (green), and GA (orange) derived from fits on the eBOSS LRG data.
In the text

	Fig. A.1. Normalized covariance matrix of the power spectrum multipoles P0, 2, 4 and the correlation function multipoles ξ0, 2, 4, measured from the 1000 cubic-box EZMOCKS of LRGs.
In the text

	Fig. A.2. Correlation coefficients between q⊥, q∥, and fσ8 for CS and FS, derived from fits to the 1000 EZMOCKS of the eBOSS LRG sample.
In the text