# Comparing Compressed and Full-modeling Analyses with FOLPS: Implications for
DESI 2024 and beyond

H. E. Noriega  0000-0002-3397-3998 A. Aviles  0000-0001-5998-3986 H. Gil-Marín
0000-0003-0265-6217 S. Ramirez-Solano  0009-0006-2096-7974 S. Fromenteau  M.
Vargas-Magaña  0000-0003-3841-1836 J. Aguilar  S. Ahlen  0000-0001-6098-7247
O. Alves  S. Brieden  0000-0003-3896-9215 D. Brooks  J. L. Cervantes-Cota
0000-0002-3057-6786 S. Chen  0000-0002-5762-6405 T. Claybaugh  S. Cole
0000-0002-5954-7903 K. Dawson  A. de la Macorra  0000-0002-1769-1640 A. de
Mattia  P. Doel  N. Findlay  0009-0007-0716-3477 J. E. Forero-Romero
0000-0002-2890-3725 E. Gaztañaga  S. Gontcho A Gontcho  0000-0003-3142-233X K.
Honscheid  J. Hou  C. Howlett  0000-0002-1081-9410 M. Ishak
0000-0002-6024-466X S. Juneau  Y. Lai  0000-0001-9054-4324 M. Landriau
0000-0003-1838-8528 M. Manera  0000-0003-4962-8934 M. Maus
0000-0002-9020-911X R. Miquel  G. Morales-Navarrete  E. Mueller  A. Muñoz-
Gutiérrez  A. D. Myers  S. Nadathur  0000-0001-9070-3102 G. Niz
0000-0002-1544-8946 N. Palanque-Delabrouille  0000-0003-3188-784X W. J.
Percival  0000-0002-0644-5727 C. Poppett  M. Rezaie  0000-0001-5589-7116 A.
Rocher  0000-0003-4349-6424 G. Rossi  E. Sanchez  0000-0002-9646-8198 D.
Schlegel  M. Schubnell  D. Sprayberry  G. Tarlé  0000-0003-1704-0781 L. Verde
0000-0003-2601-8770 S. Yuan  0000-0002-5992-7586 P. Zarrouk
0000-0002-7305-9578 H. Zou  0000-0002-6684-3997

###### Abstract

The Dark Energy Spectroscopic Instrument (DESI) will provide unprecedented
information about the large-scale structure of our Universe. In this work, we
study the robustness of the theoretical modelling of the power spectrum of
Folps, a novel effective field theory-based package for evaluating the
redshift space power spectrum in the presence of massive neutrinos. We perform
this validation by fitting the AbacusSummit high-accuracy $N$-body simulations
for Luminous Red Galaxies, Emission Line Galaxies and Quasar tracers,
calibrated to describe DESI observations. We quantify the potential systematic
error budget of Folps finding that the modelling errors are fully sub-dominant
for the DESI statistical precision within the studied range of scales.
Additionally, we study two complementary approaches to fit and analyse the
power spectrum data, one based on direct Full-Modelling fits and the other on
the ShapeFit compression variables, both resulting in very good agreement in
precision and accuracy. In each of these approaches, we study a set of
potential systematic errors induced by several assumptions, such as the choice
of template cosmology, the effect of prior choice in the nuisance parameters
of the model, or the range of scales used in the analysis. Furthermore, we
show how opening up the parameter space beyond the vanilla $\Lambda$CDM model
affects the DESI observables. These studies include the addition of massive
neutrinos, spatial curvature, and dark energy equation of state. We also
examine how relaxing the usual Cosmic Microwave Background and Big Bang
Nucleosynthesis priors on the primordial spectral index and the baryonic
matter abundance, respectively, impacts the inference on the rest of the
parameters of interest. This paper pathways towards performing a robust and
reliable analysis of the shape of the power spectrum of DESI galaxy and quasar
clustering using Folps.

##  1 Introduction

In the coming decade, large-scale structure (LSS) cosmological probes are
poised to significantly advance our understanding of the Universe. The wealth
of data expected from spectroscopic surveys, such as the Dark Energy
Spectroscopic Instrument111desi.lbl.gov (DESI) [1, 2, 3, 4], combined with
data from cosmic microwave background (CMB) experiments [5, 6, 7] and galaxy
photometric and spectroscopic surveys [8, 9, 10], holds promise not only for
the precise measurement of cosmological parameters but also for addressing
fundamental open questions in physics. These questions include elucidating the
nature of dark energy, determining the absolute mass scale of neutrinos,
testing gravitational theories at cosmological scales, and potentially
detecting primordial non-Gaussianities, among others. However, as data
accuracy and precision improve, it is essential to refine theoretical models
of structure formation to meet the standards demanded by new observations.

Spectroscopic galaxy surveys, which measure the redshifts and angular
positions of galaxies in the sky, stand out as one of the most powerful tools
for studying the properties of our Universe. These surveys probe several key
physical phenomena. These include the baryon acoustic oscillations (BAO) [11,
12], which result from baryonic fluctuations in the primordial plasma that
froze out after photons decoupled, and the redshift-space distortions (RSD)
[13, 14], introducing anisotropies in the observed galaxy distribution due to
the Doppler effect induced by their peculiar velocities. A key complication in
extracting cosmological information from the LSS arises from our inability to
directly observe dark matter. Instead, we study its effects indirectly using
biased tracers such as galaxies and the Lyman-$\alpha$ forest. Additionally,
non-linearities of the gravitational collapse become relevant at late times,
so we must consider models beyond linear order, such as perturbation theory
(PT) [15] and effective field theory (EFT) [16, 17, 18, 19, 20]. Recently,
these techniques have yielded some of the tightest constraints to date on
cosmological parameters, drawing exclusively from data collected by galaxy
surveys, as presented in [21, 22], among many other studies.

To extract cosmological information from data, there are two main kind of
methods: Compressed and Direct approaches. In the compressed approach, a
template power spectrum is constructed based on a fixed cosmology, which is
then adjusted by a set of compressed parameters (such as dilation parameters,
$f\sigma_{8}$, $\cdots$) to fit the data. These parameters effectively
condense the information of the power spectrum data vector, allowing, in
principle, for the estimation of cosmological parameters within a specific
cosmological model. This approach has the benefit that the step of compressing
the clustering information is cosmology-independent. The assumption of the
cosmological model only enters when converting from compressed to cosmological
parameters, avoiding the need for re-fitting when testing different
cosmological models. The most popular compressed approaches are the Standard
(or Classic) and ShapeFit [23]. The Standard compressed analysis was
previously used in the Baryon Oscillation Spectroscopic Survey (BOSS, [24,
25]) and extended BOSS (eBOSS, [26, 27]) collaborations and focused on
extracting cosmological information from BAO and RSD by introducing the
scaling parameters $\alpha_{\parallel,\perp}$ and the product of the
logarithmic growth of structure, $f$, and the amplitude of the dark matter
field fluctuations smoothed by a scale of $8\,h^{-1}\,{\rm Mpc}$,
$\sigma_{8}$. The ShapeFit approach aims to bridge the gap between the
constraining power of the standard compressed and direct-fitting analyses.
ShapeFit partially covers this gap by adding two shape parameters that extract
additional information from the broadband slope of the power spectrum. The
direct-fit methodology, also known as Full-Modelling, directly compares the
model to the data, varying the underlying cosmological parameters used to
generate a new linear power spectrum in each iteration of the posterior
exploration. The Full-Modelling analysis simultaneously captures the BAO and
RSD signals, as well as the rest of large-scale cosmological information
contained in the shape of the power spectrum, for an assumed cosmological
model, reaching a similar constraining power as the Planck experiment for some
cosmological parameters [21, 22, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38].

The present work examines the error budget associated with Full-Modelling and
compressed (ShapeFit) methodologies, resulting from the theoretical Modelling
of the full-shape power spectrum using the
Folps222https://github.com/henoriega/FOLPS-nu. JAX implementation:
https://github.com/cosmodesi/folpsax. code [39]. To do this, we use the
AbacusSummit high-accuracy $N$-body simulations targeted to produce DESI
galaxy mocks in three redshift bins: Luminous Red Galaxy (LRG) at $z=0.8$,
Emission Line Galaxy (ELG) at $z=1.1$, and Quasar (QSO) at $z=1.4$. The
simulations consist of 25 independent realizations per tracer, each with a
cubic box volume of $8\,\,h^{-3}\text{Gpc}^{3}$, for a total physical volume
of $200\,\,h^{-3}\text{Gpc}^{3}$ per tracer.

In this study, we explore various configurations encompassing a broad spectrum
of setups. These configurations include varying the maximum wave-number
$k_{\text{max}}$ utilized in fitting procedures, the use of the hexadecapole
in addition to the monopole and quadrupole of the power spectrum multipoles,
the effect of imposing priors on nuisance parameters, the assumption of
coevolution for the bias paramters, and the utilization of approximations to
accelerate computations. Additionally, we investigate the effect of opening
the parameter space on models beyond $\Lambda$CDM, such as those involving
massive neutrinos, dark energy (via the $w$CDM parametrization), curvature,
variation of the spectral index, $n_{s}$, and variation of the baryon
abundance $\omega_{b}$.

Our findings reveal that both the ShapeFit and Full-Modelling methodologies
yield very similar results for the vanilla $\Lambda$CDM model with fixed
$n_{s}$ and baryon abundance with a Big Bang Nucleosynthesis (BBN) prior
inspired by [40, 41], albeit exhibiting some variations for the extended
models. Moreover, we observe that all fitted parameters fall within 1 or
2-$\sigma$ of the true simulation values, even for the largest volume
considered, $200\,h^{-3}\text{Gpc}^{3}$. Importantly, our results show
significant agreement with those obtained by other DESI groups employing
PyBird333https://github.com/pierrexyz/pybird and
velocileptors444https://github.com/sfschen/velocileptors, as presented in
refs. [42] and [43], respectively. Detailed comparisons between the three
utilized codes/pipelines for Fourier space Full-Shape analysis within DESI are
presented in the companion paper [44], as well as the configuration space
analysis performed using the code gsm-
eft555https://github.com/alejandroaviles/gsm [45].

This work provides support to the Full-Shape analysis from galaxies and
quasars [46] of the DESI Data Release 1 [47] within the broader context of the
DESI Year 1 main studies [48, 49, 50, 51, 52, 53].

The rest of the paper is organized as follows. We begin by describing the
simulations used throughout the work in §2. We describe the theoretical
Modelling in §3 and the Folps pipeline in §4. In §5, we describe the different
methodologies for extracting information from galaxy clustering, namely Full-
Modelling, Standard, and ShapeFit. In §6 we present the results for different
settings on the maximal scale of the wave-number, priors on nuisance
parameters. Based on these results, we define our baseline settings. In §7 we
explore extension to the baseline analysis. Finally, in §8 we present our
conclusions.

The reader interested in the main results of this paper can skip directly to
§5, §6 and §7.

##  2 Datasets from mock catalogues

In this section, we introduce the mock catalogues used to generate the data
power spectrum multipoles and covariance matrices, which will serve us to
study the systematic effects of the pipeline and methodologies. To generate
the data multipoles, we use $N$-body simulations that provide a robust and
highly accurate representation of the clustering, but reproducing them is
computationally expensive. In contrast, to estimate the covariance matrix, we
turn to fast mocks that are computationally efficient and have the capacity to
provide an accurate representation of gravitational evolution up to
intermediate scales.

###  2.1 High-precision mocks

We use synthetic data obtained from the
AbacusSummit666https://abacussummit.readthedocs.io high-accuracy $N$-body
simulations [54], produced with the Abacus $N$-body code [55]. These
simulations were designed to meet (and exceed) the requirements set by the
DESI cosmological simulations [54]. Further, as mentioned in the Introduction,
the companion paper [44] shows a remarkable agreement between Folps,
velocileptors and PyBird codes, as well as the simulations, which persists
even when considering the full available volume.

Most of these AbacusSummit simulations consist of cubic boxes with a volume of
$8\,h^{-3}\text{Gpc}^{3}$ and $6912^{3}$ total particles, resulting in an
individual particle mass of approximately $2\times
10^{9}\,h^{-1}\text{M}_{\odot}$. The suite of simulations includes 97
cosmological models, including the Planck 2018 model and variations, as
outlined in [54]. Our analysis focuses on three types of tracers: LRGs at a
redshift of $z=0.8$, ELGs at $z=1.1$, and QSOs at $z=1.4$, which consist of a
collection of 25 realizations for each tracer, with each realization adopting
the fiducial cosmology provided by Planck 2018:
$\\{h=0.6736,\,\omega_{cdm}=0.1200,\,\omega_{b}=0.02237,\,\ln(10^{10}A_{s})=3.0364,\,n_{s}=0.9649,\,M_{\nu}=0.06\,\text{eV},\,w_{0}=-1\\}$.
Each individual realization has a volume of $8\,h^{-3}\text{Gpc}^{3}$, and the
$k$-bins have a spacing of $\Delta k=0.005\,h{\rm Mpc}^{-1}$. Consequently,
the total volume reaches $200\,h^{-3}\text{Gpc}^{3}$, ensuring remarkably
small statistical errors and fulfilling the desired statistical precision of
DESI mocks. During the analyses, to reduce sample variance ‘noise’, we fit the
average of the 25 realizations.777In contrast to the DESI-galaxy BAO analysis
in [49, 56], in this work we refrain from employing the control variate
technique designed to mitigate the noise stemming from sample variance [57,
58].

###  2.2 Covariance matrix

We employ the Effective Zel’dovich mocks (EZmocks) to estimate the covariance
matrices. These mocks consist of 1000 independent synthetic realizations for
LRG, ELG galaxies, and QSOs. As the name indicates, these mocks are generated
using fast techniques based on the Zel’dovich approximation [59], as presented
in [60, 61]. The advantage of these mocks lies in their computational
efficiency compared to $N$-body simulations while still providing a good
approximation to the two- and three-point statistics up to mildly non-linear
scales. Thus, they serve as a cheap and reliable method for generating
extensive mock catalogues, enabling the estimation of reliable covariance
matrices via

$\textit{Cov}_{\ell,\ell^{\prime}}(k_{i},k_{j})=\frac{1}{N_{m}-1}\sum_{n=1}^{N_{m}}\big{[}P_{\ell,n}(k_{i})-\mu_{\ell}(k_{i})\big{]}\big{[}P_{\ell^{\prime},n}(k_{j})-\mu_{\ell^{\prime}}(k_{j})\big{]},$
(2.1)

where $\mu_{\ell}$ indicates the mean power spectrum,

$\mu_{\ell}(k)\equiv\frac{1}{N_{m}}\sum_{n=1}^{N_{m}}P_{\ell,n}(k),$ (2.2)

while $N_{m}$ is the number of mock realizations and $P_{\ell,n}(k)$ denotes
the power spectrum multipole of the $n$-th mock. We opt to follow the
procedure of correcting the covariance by recipe provided by [62]. Although
this approach does not guarantee obtaining an unbiased posterior on the
parameters of the model (only an unbiased estimate of the inverse covariance),
for the precision of the settings here (i.e. 1000 realizations, a data vector
of 72 elements, which only corrects the amplitude of the covariance by 8%)
other more precise approaches [63, 64, 65] only represent minor corrections on
the posteriors of the parameters. Thus, we scale the inverted covariance
matrix by,

$\textit{Cov}^{-1}_{\ell,\ell^{\prime}\,\text{Hartlap}}=\frac{N_{m}-n_{b}-2}{N_{m}-1}\,\textit{Cov}^{-1}_{\ell,\ell^{\prime}},$
(2.3)

where $n_{b}$ is the number of data bins.

The volume of each simulation is $V_{1}=8\,\,h^{-3}\text{Gpc}^{3}$. In this
work we re-scale the covariance by factors 1/5 and 1/25, the latter being the
maximum available by the simulations, to have a total volume

$\displaystyle V_{5}$ $\displaystyle=40\,\,h^{-3}\text{Gpc}^{3},$ (2.4)
$\displaystyle V_{25}$ $\displaystyle=200\,\,h^{-3}\text{Gpc}^{3}.$ (2.5)

Throughout this paper, we primarily work with volumes of $V_{5}$ and $V_{25}$,
as well as the minimal volume of $V_{1}=8\,\,h^{-3}\text{Gpc}^{3}$ for a few
selected tests. One primary motivation for utilizing a volume larger than the
DESI physical volume is to detect potential systematic errors in the model
above the usual statistical noise ($\sim 2\sigma$) of the DESI Year-1 or
Year-5 datasets, DESI-Y1 or DESI-Y5, respectively. Under the assumption of a
perfectly unbiased model, we establish that sample variance fluctuations of a
simulation of volume $V_{\rm sim}$ can typically produce noise-shifts which
falls within a range of $2\sigma_{\rm sim}$. This difference should not exceed
the comfort limit of DESI, denoted as $\frac{1}{n}\sigma_{\text{DESI}}$, with
$1/n$ being the fractional threshold which defines when a systematic becomes
relevant given your statistical precision. We estimate that $n$ should be at
least 3, but preferably $\geq 5$. Thus, we establish that our simulations
should have a volume such that, $2\sigma_{\rm
sim}\leq\tfrac{1}{n}\sigma_{\text{DESI}}$. Hence, the optimal effective volume
of simulations is given by $V_{\rm sim}=(2n)^{2}\,V_{\text{DESI}}$; for the
DESI Y5 volume, estimated around $5\,\,h^{-3}\text{Gpc}^{3}$ for a given
tracer, we get that $V_{\rm sim}\simeq 180\,\,h^{-3}\text{Gpc}^{3}$, for the
conservative choice of $n=3$; while for the DESI Y1 volume, estimated to be
1/5 of the DESI Y5, we obtain $V_{\rm sim}\simeq 40\,\,h^{-3}\text{Gpc}^{3}$.
These numbers motivate the choices of volumes for re-scaling the covariance
provided above.

Despite potential uncertainties about the accuracy of simulations for
$V=200\,h^{-3}\text{Gpc}^{3}$, ref. [44] compares different Modelling
approaches with these simulations, finding remarkable agreements between the
models and with the simulations, even for the largest volume. This suggests
that both the models and AbacusSummit simulations are reliable, even in these
extreme scenarios, despite the expected systematic errors.

##  3 Redshift space power spectrum modelling

In redshift space, the position of a galaxy located at an Eulerian real space
coordinate $\bm{x}$ is distorted because of its peculiar velocities $\bm{v}$
relative to the Hubble flow, such that it appears to be located at a redshift
space coordinate $\bm{s}$. The map between these coordinates is given by

$\bm{s}(\bm{x})=\bm{x}+\frac{\hat{\bm{n}}\cdot\bm{v}(\bm{x})}{aH}\hat{\bm{n}}.$
(3.1)

We have adopted the distant observer approximation on which $\hat{\bm{n}}$ is
a unit vector in the direction of the sample of the observed survey, instead
of the direction of each individual galaxy $\hat{\bm{x}}$. Under the above
coordinate transformation, the galaxy number density fluctuation in redshift
space, $\delta_{s}$, is given by

$(2\pi)^{3}\delta_{\text{D}}(\bm{k})+\delta_{s}(\bm{k})=\int
d^{3}x\left(1+\delta_{g}(\bm{x})\right)e^{i\bm{k}\cdot\bm{s}(\bm{x})},$ (3.2)

where $\delta_{g}$ is the real space galaxy number density fluctuation, and
$\delta_{\text{D}}$ is the 3-dimensional Dirac delta function. Using the
momentum expansion approach presented in [66, 67], the power spectrum can be
written as

$(2\pi)^{3}\delta_{\text{D}}(\bm{k})+P_{s}(\bm{k})=\sum_{m=0}^{\infty}\frac{(-i)^{m}}{m!}(k\mu)^{m}\tilde{\Xi}^{(m)}(\bm{k})$
(3.3)

with $\mu\equiv\hat{\bm{k}}\cdot\hat{\bm{n}}$ the angle cosine of the wave-
number and the line of sight direction, and $\tilde{\Xi}^{(m)}$ is the
velocity weighted density momentum of degree $m$,

$\tilde{\Xi}^{(m)}(\bm{k})=\frac{1}{(aH)^{m}}\int
d^{3}x\,e^{-i\bm{k}\cdot\bm{x}}\,\langle\big{(}1+\delta_{g}(\bm{x}_{2})\big{)}\big{(}1+\delta_{g}(\bm{x}_{1})\big{)}(\hat{\bm{n}}\cdot\Delta\bm{v})^{m}\rangle,$
(3.4)

with $\bm{x}=\bm{x}_{2}-\bm{x}_{1}$ and
$\Delta\bm{v}=\bm{v}(\bm{x}_{2})-\bm{v}(\bm{x}_{1})$. Each one of these
moments can be brought into the integral form [68]

$\tilde{\Xi}^{(m)}(\bm{k})=\sum_{n}\int
d^{3}p\,(\hat{\bm{p}}\cdot\hat{\bm{n}})^{n}S_{n}^{(m)}(\bm{k},\bm{p}),$ (3.5)

for some rotational invariant functions $S_{n}^{(m)}(\bm{k},\bm{p})\equiv
S_{n}^{(m)}(k,p,x)$ with $x\equiv\hat{\bm{k}}\cdot\hat{\bm{p}}$. Then, one can
expand [68, 69]

$\int\frac{d^{3}p}{(2\pi)^{3}}(\hat{\bm{p}}\cdot\hat{\bm{n}})^{n}S_{n}^{(m)}(\bm{k},\bm{p})=\sum_{m=0}^{n}\mu^{m}\int\frac{d^{3}p}{(2\pi)^{3}}G_{nm}(x)S_{n}^{(m)}(k,p,x)$
(3.6)

with

$\displaystyle G_{nm}(x)$
$\displaystyle=\sum_{\ell=0}^{n}\frac{(1+(-1)^{\ell+n})(2\ell+1)}{2(1+\ell+n)}\binom{\ell}{m}\binom{2\ell}{\ell}\binom{\frac{\ell+m-1}{2}}{\ell}$
$\displaystyle\quad\times\,{}_{3}F_{2}\left(\frac{1-\ell}{2},-\frac{\ell}{2},\frac{1}{2}(-1-\ell-n);\frac{1}{2}-\ell,\frac{1}{2}(1-\ell-n);1\right)\mathcal{L}_{\ell}(x),$
(3.7)

where ${}_{3}F_{2}({\text{\bf a}};{\text{\bf b}};z)$ is the generalized
hypergeometric function of the kind ($p=3$, $q=2$) [70], and
$\mathcal{L}_{\ell}(x)$ is the Legendre polynomial of degree $\ell$. Equation
(3.6) allows us to write the power spectrum in eq. (3.3) as

$P_{s}(k,\mu)=\sum_{m=0}^{\infty}\sum_{n=0}^{m}\mu^{2n}f^{m}I_{mn}(k)$ (3.8)

with $f$ the linear growth rate and $I_{mn}(k)$ a set of functions that only
depend on the wave-number $k$. To obtain the one-loop power spectrum, one must
truncate the sum at $m=4$ ($m>4$ gives higher than one-loop terms, as can be
read from eq. (3.4)). The leading order, tree-level expansion reduces to the
Kaiser power spectrum [13], while the one-loop pieces of functions $I_{mn}(k)$
can be written either as [68]

$\displaystyle I_{mn}(k)$ $\displaystyle=\int
d^{3}p\,\mathcal{K}_{mn}(\bm{k},\bm{p})P_{L}(|\bm{k}-\bm{p}|)P_{L}(p),$ (3.9)
or $\displaystyle I_{mn}(k)$ $\displaystyle=P_{L}(k)\int
d^{3}p\,\mathcal{K}_{mn}(\bm{k},\bm{p})P_{L}(p),$ (3.10)

that we call $P_{22}$-type and $P_{13}$-type, respectively, since they have
the form of the contributions $P_{22}$ and $P_{13}$ to the one-loop real space
power spectrum [15]. In the case of Einstein-de Sitter evolution, or more
precisely when $f(a)=\Omega^{1/2}_{m}(a)$, the kernels $\mathcal{K}_{mn}$ are
known analytically, and the integrals can be solved via standard methods in
perturbation theory. In particular, the Folps code uses Fast Fourier Transform
in logarithmic spaced intervals (FFTLog) [71, 72, 73, 74, 75, 76, 77] to speed
up these computations (see section 4 and Appendix C).

Keeping only up to one-loop corrections one obtains

$P_{s}(k,\mu)=P_{\delta\delta}(k)+2f\mu^{2}P_{\delta\theta}(k)+f^{2}\mu^{4}P_{\theta\theta}(k)+A^{\text{TNS}}(k,\mu)+D(k,\mu),$
(3.11)

where $A^{\text{TNS}}$ is the standard function “$A$” in the Taruya-
Nishimichi-Saito (TNS) model (see eqs. (19) and (A3) of [78]) constructed out
of three-point correlators, and $D(k,\mu)$ is given by eq. (3.42) of [39],
which is constructed from four-point correlators and hence by only linear
fields at one-loop.

We adopt the definition

$\theta(\bm{x})=-\frac{\bm{\nabla}\cdot\bm{v}(\bm{x})}{aHf}$ (3.12)

for the divergence of the velocity, and hence, $P_{\delta\delta}$,
$P_{\delta\theta}$ and $P_{\theta\theta}$ are the density-density, density-
velocity and velocity-velocity one-loop power spectra, respectively. With the
definition of equation (3.12), the velocity and density fields become equal at
the linear order, i.e. $\delta^{(1)}=\theta^{(1)}$. This holds true for the
$\Lambda$CDM cosmology, where neutrinos are approximated as massless, and
hence, the growth rate becomes only a function of time. This is a reasonable
approximation when the total sum of their masses is small. Indeed, for a total
neutrino mass of $0.06\,\text{eV}$, using Einstein-de Sitter kernels produces
results that are almost indistinguishable from those obtained using the full
theoretical framework [79, 80].

###  3.1 Bias, Effective Field Theory and Shot noise

At large scales, galaxies are biased tracers of the combined cold dark matter
plus baryons field, and not the total matter field which also includes the
massive neutrinos [81]. To relate galaxy and matter (cold+baryonic, $cb$)
overdensities we use the biasing prescription of [17], which after
renormalization yields the equations [82]

$\displaystyle P_{\delta\delta}(k)$
$\displaystyle=b_{1}^{2}P^{\text{1-loop}}_{cb,\delta\delta}(k)+2b_{1}b_{2}P_{b_{1}b_{2}}(k)+2b_{1}b_{s^{2}}P_{b_{1}b_{s^{2}}}(k)+b_{2}^{2}P_{b_{2}^{2}}(k)$
$\displaystyle\quad+2b_{2}b_{s^{2}}P_{b_{2}b_{s^{2}}}(k)+b_{s^{2}}^{2}P_{b_{s^{2}}^{2}}(k)+2b_{1}b_{3\rm
nl}\sigma^{2}_{3}(k)P^{L}_{cb,\delta\delta}(k),$ (3.13) $\displaystyle
P_{\delta\theta}(k)$
$\displaystyle=b_{1}P^{\text{1-loop}}_{cb,\delta\theta}(k)+b_{2}P_{b_{2},\theta}(k)+b_{s^{2}}P_{b_{s^{2}},\theta}(k)+b_{3\rm
nl}\sigma^{2}_{3}(k)P^{L}_{cb,\delta\theta}(k),$ (3.14) $\displaystyle
P_{\theta\theta}(k)$ $\displaystyle=P^{\text{1-loop}}_{cb,\theta\theta}(k),$
(3.15)

which enter into eq. (3.11). Expressions for all the functions $P_{XY}$ can be
found in several articles (e.g. [83, 84]), although these can differ a little
bit in the definitions. Here, we use eqs. (3.43)–(3.49) of [68]. Further, the
function $\sigma^{2}_{3}$, which accompanies the third-order non-local bias
parameter $b_{3\rm nl}$, is given by

$\sigma^{2}_{3}(k)=\frac{105}{16}\int\frac{d^{3}p}{(2\pi)^{3}}P_{L}(p)\left[S_{2}(\bm{p},\bm{k}-\bm{p})\left(\frac{2}{7}S_{2}(-\bm{p},\bm{k})-\frac{4}{21}\right)+\frac{8}{63}\right],$
(3.16)

with $S_{2}(\bm{k}_{1},\bm{k}_{2})=\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}-1/3$.

On the other hand, the functions $A^{\text{TNS}}$ and $D$ are affected by the
linear bias in the following way,

$\displaystyle A^{\text{TNS}}(k,\mu;f)\quad\rightarrow$ $\displaystyle\quad
b_{1}^{3}A^{\text{TNS}}(k,\mu;f/b_{1}),$ (3.17) $\displaystyle
D(k,\mu;f)\quad\rightarrow$ $\displaystyle\quad b_{1}^{4}D(k,\mu;f/b_{1}).$
(3.18)

While the identification of the biased function $D$ is exact, $A^{\text{TNS}}$
is also modified in the presence of second-order local bias and tidal bias.
However, these effects are subdominant in comparison to those in the $P_{XY}$
spectra, so we keep them out for simplicity. They can be added using the
complete expression for the biased $A^{\text{TNS}}$ that can be found in eq.
(A.9) of Appendix A.1 of [68].

The power spectrum presented in the preceding subsection, along with various
summary statistics, originates from the theory of fluid equations for cold
dark matter. Nevertheless, the underlying theory that governs a collection of
particles subject to gravitational attraction is described by the Boltzmann
equation. However, the fluid approximation loses its validity beyond a
specific scale determined by the velocity dispersion and the breakdown of the
single-stream approximation. Consequently, we expect our theory to remain
applicable only well beyond this scale. Furthermore, while small scales evade
a perturbative description, understanding their impact on larger scales
through backreaction effects is essential for uncovering a wealth of
cosmological information. To accommodate these effects, the standard procedure
involves the regularization of loop integrals with a cutoff scale.
Additionally, a series of counterterms are introduced to eliminate the cutoff
scale dependence from the final expressions.888In practice we do not
regularize the integrals since we use FFTLog. This is allowed because the loop
integrals do converge, and UV modes yield only small contribution. This
approach, which is the most common in the literature, becomes equivalent to a
re-scaling of the EFT counterterms. This framework, known as Effective Field
Theory(ies) for Large Scale Structure (EFT), was initially proposed in [16,
17, 18] and has since been further developed in subsequent papers; see [85]
for a recent review.

Finally, one has to account for the stochastic noise associated with the
discreteness of observations and the finiteness of the number of galaxies in a
sample. Hence, our final expression for the EFT power spectrum is

$\displaystyle P^{\text{EFT}}_{s}(k,\mu)$
$\displaystyle=P_{s}(k,\mu)+(\alpha_{0}+\alpha_{2}\mu^{2}+\alpha_{4}\mu^{4})k^{2}P_{L}(k)+P_{\text{\rm
shot}}\Big{(}\alpha^{\text{\rm shot}}_{0}+(k\mu)^{2}\alpha^{\text{\rm
shot}}_{2}\Big{)},$ (3.19)

with $P_{s}$ the redshift space perturbation theory power spectrum given in
eq. (3.11), while $\alpha_{0,2,4}$ are the EFT counterterms, and $P_{\text{\rm
shot}}$ is a constant that its value can be fixed at choice. For example, it
can be instructive to set it to the inverse of the number density of galaxies,
$1/\bar{n}$, and since the constant part of the shot noise is given by
$P_{\text{\rm shot}}\times\alpha^{\text{\rm shot}}_{0}$, the parameter
$\alpha^{\text{\rm shot}}_{0}\neq 1$ quantifies the departure from a
Poissonian noise. This can be helpful to impose a Gaussian prior around the
Poisson noise since we expect $\alpha^{\text{\rm shot}}_{0}$ to be of order
unity.

###  3.2 IR-resummations and multipoles

The smearing and degradation of the BAO oscillations [86, 87, 88, 89] is not
described by SPT, though it is well captured by LPT even at its lower order,
the Zeldovich approximation [90, 91, 20, 92, 34, 56]. In order to model them,
we utilize IR-resummation methods [93, 20]. Specifically, we adopt the
approach outlined in [94, 95, 29], which involves splitting the linear power
spectrum into two components as $P_{L}=P_{nw}+P_{w}$: the non-wiggle power
spectrum ($P_{nw}$), devoid of the Baryon Acoustic Oscillations (BAO), and the
wiggle component ($P_{w}$). We construct the one-loop IR-resummed EFT
redshift-space power spectrum as [95]

$\displaystyle P_{s}^{\text{IR}}(k,\mu)$
$\displaystyle=e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}P_{s}^{\text{EFT}}(k,\mu)+\big{(}1-e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}\big{)}P_{s,nw}^{\text{EFT}}(k,\mu)$
$\displaystyle\quad+e^{-k^{2}\Sigma^{2}_{\text{tot}}(k,\mu)}P_{w}(k)k^{2}\Sigma^{2}_{\text{tot}}(k,\mu).$
(3.20)

The full component, denoted as $P_{s}^{\text{EFT}}(k,\mu)$, corresponds to the
one-loop power spectrum calculated using equation (3.19). On the other hand,
the non-wiggle part, represented as $P_{s,nw}^{\text{EFT}}(k,\mu)$, is also
obtained using equation (3.19), but with the non-wiggle linear power spectrum
$P_{nw}$ as the input. The function $\Sigma^{2}_{\text{tot}}$ is defined by

$\Sigma^{2}_{\text{tot}}(k,\mu)=\big{[}1+f\mu^{2}\big{(}2+f\big{)}\big{]}\Sigma^{2}+f^{2}\mu^{2}(\mu^{2}-1)\delta\Sigma^{2},$
(3.21)

with

$\displaystyle\Sigma^{2}$
$\displaystyle=\frac{1}{6\pi^{2}}\int_{0}^{k_{\text{IR}}}dp\,P_{nw}(p)\left[1-j_{0}\left(p\,\ell_{\text{BAO}}\right)+2j_{2}\left(p\,\ell_{\text{BAO}}\right)\right],$
(3.22) $\displaystyle\delta\Sigma^{2}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{k_{\text{IR}}}dp\,P_{nw}(p)j_{2}\left(p\,\ell_{\text{BAO}}\right),$
(3.23)

where $j_{n}$ is the spherical Bessel functions of degree $n$, and
$\ell_{\text{BAO}}\approx 105\,h^{-1}\text{Mpc}$ is the BAO scale. The wave-
number $k_{\text{IR}}$ represents a transition scale between long and short
modes, and it is somewhat arbitrary; however, the final results exhibit weak
dependence on it as long as $k_{\text{IR}}\gtrsim 0.1\,h\,\text{Mpc}^{-1}$, as
shown in figure 22 of Appendix C.

Equation (3.2) is our ultimate model for the power spectrum. To fit the data,
we take its monopole, quadrupole and hexadecapole from

$P_{\ell}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}d\mu\;P_{s}^{\text{IR}}(k,\mu)\mathcal{L}_{\ell}(\mu),$
(3.24)

where $\mathcal{L}_{\ell}$ are the Legendre polynomial of degree $\ell$.

##  4 FOLPS code

To compute the multipoles of the redshift space power spectrum we use
Folps$\nu$ (Fast One Loop Power Spectrum in the presence of massive neutrinos)
[39].999https://github.com/henoriega/FOLPS-nu. JAX implementation:
https://github.com/cosmodesi/folpsax. Folps is a code written entirely in
Python with the aid of built-in numpy and scipy functions [96, 97]. It
computes the multipoles of the IR-resummed power spectrum of eq. (3.24) in a
fraction of second.101010In a standard personal computer, this time is about
$\sim 0.2\,\text{sec}$, using $N_{\text{FFT}}=128$. We refer the reader to
Appendix C or fig. 9 of [39]. The code can run in three different setups:
Full-Modelling, Standard, and ShapeFit fits. Broadly speaking, the former
involves a direct-fitting of cosmological parameters to the data. On the other
hand, the Standard and ShapeFit approaches employ two-step methodologies. In
the first step, information from clustering is compressed into a set of
parameters, which are subsequently used to infer the cosmological parameters.
Section §5 explores in detail each of these methodologies.

The general pipeline of the Folps code consists of the following steps:

  1. 1.

It receives as input the linear matter power spectrum of the $cb$ field, and
the set of cosmological parameters: $\Omega_{m}$, $h$, and $M_{\nu}$, all of
which are required for beyond Einstein-de Sitter (EdS) kernels, as explained
below in §4.1. Note that the growth rate $f$ can be treated as a derived
parameter of the model.111111As it will be explicitly described in Section §5,
$f$ is a derived parameter in the case of the direct-fit approach but is
treated as a free parameter for the case of the compressed fit approach.
Additionally, one needs to feed the code with biases, counterterms and
stochastic parameters,

$\\{b_{1},b_{2},b_{s^{2}},b_{3\rm
nl},\alpha_{0},\alpha_{2},\alpha_{4},\alpha_{0}^{\text{\rm
shot}},\alpha_{2}^{\text{\rm shot}}\\}.$ (4.1)

  2. 2.

The code extrapolates the input linear matter power spectrum and computes the
non-wiggle linear power spectrum using the fast sine transform method
presented in [98, 29].

  3. 3.

It computes the EFT power spectrum, given by eq. (3.19), where the loop
integrals of the form (3.9) and (3.10) are calculated using FFTLog methods
[72, 73, 74, 75, 76, 77, 29]. This is done separately for the regular linear
power spectrum and for the linear power spectrum with the wiggles removed, and
thereafter mixed together in the IR-resummed power spectrum of eq. (3.2) [95].

  4. 4.

Finally, it calculates the Legendre multipoles using the integral in eq.
(3.24).

In practice, eq. (3.24) includes scaling parameters to account for the
distance dilation of scales of the measured power spectrum along and across
the line of sight. This dilation is caused by the – a priori – unknown mapping
between redshift-to-distances in a catalogue provided in terms of angles and
redshifts. This is the case for analyses of real data (or synthetic data
mimicking observational effects, e.g., a simulated lightcone), where the power
spectrum is computed by assuming a fiducial cosmology to transform redshifts
and angles into Cartesian distances, which does not necessarily match the true
cosmology. Conversely, for analyses of simulated data where the positions are
given in physical Cartesian coordinates (for e.g., a cubic mock), there is no
such dilation. However, in the direct-fit analysis, it is also common to
include these scaling parameters to fully mimic the effect on the parameter
space expected from actual catalogues. On the other hand, the Standard
compressed fits and ShapeFit compressed fits setups naturally include these
dilation parameters among their compressed set of parameters, as we will
discuss in more detail in §5.

Folps incorporates the option to reduce the dimensionality of the parameter
space within a Markov Chain Monte Carlo (MCMC) pipeline by analytically
marginalizing over some nuisance parameters. Then, the above flowchart is
slightly modified, as explained in Appendix A.

Finally, the code fixes the parameters related to IR-resummations to
$\ell_{\text{BAO}}=104\,h^{-1}\text{Mpc}$ and
$k_{\text{IR}}=0.4\,h\,\text{Mpc}^{-1}$. In Appendix C we explore the
robustness of the code by changing these setups.

###  4.1 Beyond Einstein-de Sitter kernels

In the presence of additional physical scales, such as the neutrino mass, the
linear growth function $D_{+}$ and growth rate $f$ are no longer scale-
independent. This implies that EdS, or more generally $\Lambda$CDM, kernels
are no longer strictly correct. Folps uses the fkPT method introduced in [80,
39, 99] for dealing with the new scale introduced by the neutrino mass.

The main difference from the standard EdS treatment comes from the fact that
linear velocities and overdensities relate non-locally in the presence of
massive neutrinos. Indeed, in Fourier space we have [79]

$\theta^{(1)}(\bm{k},t)=\frac{f(k,t)}{f_{0}(t)}\delta^{(1)}(\bm{k},t)$ (4.2)

where $f(k,t)$ is the time- and scale-dependent linear growth rate, and
$f_{0}$ is a very large scale limit, where neutrinos behave as cold dark
matter. This linear relation is inherited to PT velocity kernels $G_{n}$. The
most easy to grasp phenomenon is the advection of large-scale velocity fields:
e.g., at second order, the velocity field receives corrections from the large-
scale Lagrangian displacements $\bm{\Psi}$, such that
$\theta^{(2)}(\bm{x})\ni\theta^{(1)}(\bm{x}+\bm{\Psi})-\theta^{(1)}(\bm{x})=\Psi^{(1)}_{i}(\bm{x})\partial_{i}\theta^{(1)}(\bm{x})=\partial_{i}\big{[}\nabla^{-2}\delta^{(1)}(\bm{x})\big{]}\partial_{i}\theta^{(1)}(\bm{x})$.
Transforming to Fourier space, symmetrizing over outward wavevectors and using
eq. (4.2), we obtain the contribution to the second order $G_{2}$ kernel

$\frac{1}{2}\left(\frac{f(k_{1})}{f_{0}}\frac{k_{1}}{k_{2}}+\frac{f(k_{2})}{f_{0}}\frac{k_{2}}{k_{1}}\right)\in
G_{2}(\bm{k}_{1},\bm{k}_{2}).$ (4.3)

The full expressions for $G_{2}$ and $G_{3}$ “fk-kernels” and corresponding
modifications to the standard FFTLog method can be found in [80, 39].

We compute the factor $f(k)/f_{0}$ using the approximation of Hu-Eisenstein
presented in [100], while to isolate $f(k)$ we obtain $f_{0}$ by solving the
growth function differential equation in the limit $k=0$. It turns out that
the shape of $f(k)$ changes very little with redshift, and only the overall
normalization given by $f_{0}$ is relevant. When performing a joint analysis
of tracers, we exploit this fact and compute the linear power spectrum and
loop corrections at a single redshift, and then re-scale them with the linear
growth function at large scales $D_{+}(k\rightarrow 0,z)$.

The above formalism seems very cumbersome and the differences in the results
using EdS are indeed very small for $m_{\nu}\sim 0.06\,\text{eV}$. However,
when neutrino masses are sampled over wider ranges, the differences may be
important, as we explore in §7.4. Furthermore, the implementation of fk-
kernels is straightforward and computationally efficient. This is because with
the FFTLog formalism, we transform the evaluation of the loop integrals
$I_{mn}$ from eq. (3.8), into matrix multiplications of the form $P_{22}$-type
or $P_{13}$-type, depending on whether they arise from $P_{22}$-type or
$P_{13}$-type loop integrals [77, 39]. This formalism allows us to
significantly speed up the loop computations compared to the usual direct
integration, where calculations are performed with brute force. This is made
possible by employing highly optimized algorithms for matrix multiplication,
which are commonly found in standard numerical libraries. To incorporate
effects beyond EdS, the code only requires computing $2\times 26$
$P_{22}$-type matrix multiplications in the FFTLog formalism, in contrast to
the total of $2\times 24$ matrix multiplications of this type needed for EdS
kernels.121212The factor of 2 appears because we have to perform matrix
multiplications for power spectra with and without wiggles. Similarly, the
number of operations for the computationally less expensive $P_{13}$-type
contributions increases from $2\times 7$ to $2\times 11$ when including
effects beyond EdS. Therefore, the time saved when using EdS kernels is
negligible [39].

##  5 Extracting cosmological information from galaxy clustering

So far we have described the details of how Folps models the non-linear galaxy
power spectrum given a set of free nuisance parameters (see eq. (4.1)), and
given the linear cdm+baryon power spectrum, $P_{L}^{cb}(k)$, which the codes
take as input for evaluating the one-loop integrals as described in §3. The
information in terms of cosmological parameters can be essentially extracted
in two approaches: performing a direct-fit, or performing a fit in terms of
compressed parameters. On one hand, the direct-fit approach starts assuming a
cosmological model (for eg., $\Lambda$CDM, $\omega\Lambda$CDM, $k\Lambda$CDM,
…) with a set of free cosmological parameters, $\bf\Omega$ (typically
$\Omega_{m}$, $h$, $A_{s}$,…), and a set of informative priors – flat or
Gaussian – that may be applied on those parameters that are not very well
constrained by the data (typically $\Omega_{b}$, $n_{s}$, $\tau$, …). From
this given set of parameters of the model, the full
$P_{L}^{cb}(k|{\bf\Omega})$ is inferred through a Boltzmann code (in this case
CLASS131313https://lesgourg.github.io/class_public/class.html [101]), and the
full non-linear power spectrum is computed (given a set of nuisance
parameters), compared to the data, and eventually the posteriors of the
parameters of the models extracted. On the other hand, the compressed
parameter fit approach relies on defining the linear power spectrum function,
from a fixed-template linear power spectrum - evaluated at some reference
cosmology - plus a set of free parameters $\bf p$ which modify that template
and the redshift-space and non-linear terms. Similar to the direct-fit case,
the posteriors of the compressed parameters ${\bf p}$, as well as the rest of
the nuisance parameters, are extracted by comparing the non-linear power
spectrum to the data.

It is important to note that in both approaches the form of the linear power
spectrum is changed, either by ${\bf\Omega}$ or by ${\bf p}$, and the
difference between the two is the type of freedom that each approach provides
to the linear power spectrum. Once a model is assumed, the compressed set of
${\bf p}$ parameters can be mapped to the parameters of the model $\bf\Omega$
and vice-versa. However, depending on the specific choice of the
parametrization ${\bf p}$ and the details of the model which $\bf\Omega$
describe, this mapping may not be invertible. The ${\bf p}$ parametrization
may not be lossless (hence the $\bf\Omega$ parameters may contain more
information than $\bf p$), but internal model priors can also add information
not contained by ${\bf p}$. This was explored for the first time in [102] for
a vanilla $\Lambda$CDM model.

As we will see in this section, in the direct-fit approach the power spectrum
shape changes in tandem with the Alcock-Paczynski effect and the growth rate
in a highly constrained manner. In the compressed parameter approach instead,
the power spectrum shape is in general less tight to these physical effects,
and therefore both approaches will be making fairly different assumptions at
the fitting step. However, at the interpretation step, one can end up imposing
equivalent physical conditions, which should bring both approaches to retrieve
similar results.

Finally, we highlight that each of these approaches has its own advantages and
disadvantages. Importantly the two approaches address slightly different
questions. What can the data tell us about the parameters of a model (assuming
this is the actual model for the data being analysed), for the direct-fits?
What are the model-agnostic features that the data constrain, and how do these
features relate to the parameters of a given model, for the compressed type of
fits?

Below we describe in more detail how these different types of fits are
performed. In Table 1, we classify and summarize the methodological aspects of
each of these types of fits.

###  5.1 Full-Modelling analysis

The Full-Modelling or direct-fit approach aims to extract all the cosmological
information contained by the full shape of the power spectrum. This includes
not only the BAO and RSD but also the rest of the features contained by the
shape such as the large-scale amplitude, its slope or the turnaround peak.
Each of these features encodes different physical processes at different
epochs: recombination epoch for the BAO, late-time physics for the RSD,
matter-radiation equality for the turnaround, etc. Thus, in this approach it
is not easy to isolate the constraints coming from the different pieces
contained within $\Lambda$CDM-like models.

In the literature, we can find implementations of the Effective Field Theory
[18, 19, 20, 103] on extracting cosmological information from spectroscopic
surveys (such as BOSS and eBOSS) employing the Full-Modelling approach [22,
21]. In this approach, the linear shape of the power spectrum changes as the
code explores the posterior of the parameters of interest. This requires
calling in each step of the analysis, an Einstein-Boltzmann code to generate
the linear power spectrum, and a non-linear code, Folps in our case, to
evaluate the high-order loop corrections to eventually produce the non-linear
power spectrum, which makes the whole analysis very time-consuming. This issue
is notably mitigated by using efficient methods such as FFTLog [72], which
boosts the performance of the one-loop order calculations, and an analytic
marginalization over some of the nuisance parameters, which significantly
reduces the number of necessary steps for reaching convergence. These points
are described in more detail in section §4 and Appendix A. Furthermore, one
can always expedite computational performance even more by emulating the
linear and non-linear power spectra, as used e.g. in [42, 43] and [45], using
Taylor series emulators and neural network accelerators, respectively.

In the iteration process described above, the parameters sensitive to the
background cosmology (such as $\Omega_{i}$, with $i=m,\,\Lambda,\,k,\ldots$)
change. These parameters determine the redshift-to-distance conversion of the
data catalogue, where the positions are given in terms of redshift and two
angles, into the physical catalogue, where the positions are given in
Cartesian comoving coordinates, and from which the power spectrum is computed.
To avoid recomputing the power spectrum in each step, it is common to fix the
cosmology at which the data catalogue is converted into the physical catalogue
and the power spectrum measured. This means that the power spectrum is
computed in terms of the ‘observed’ wave vectors $k$’s, instead of the true
wave vectors $k^{\prime}$’s. If we decompose the 3D wave-vector on modes along
and across the line-of-sight, $k_{\parallel}$ and $k_{\perp}$, respectively,
we can write,

$k^{\prime}_{\parallel}=k_{\parallel}/q_{\parallel}\,,\quad\text{and}\quad
k_{\perp}^{\prime}=k_{\perp}/q_{\perp},$ (5.1)

where true and observed wave-vectors are related by the scaling parameters
$q_{\parallel,\perp}$. These parameters can be inferred by knowing the
cosmology at which the physical catalogue has been generated (referred to as
‘fiducial’ cosmology) and the explored cosmology at each step of the posterior
exploration, ${\bf\Omega}$. These parameters are then the ratio of the Hubble
distance, $D_{H}(z)=c/H(z)$ and the angular diameter distance, $D_{A}(z)$ at
the redshift $z$ of the sample,

$\displaystyle
q_{\parallel}(z,{\bf\Omega})=\frac{D_{H}^{\bf\Omega}(z)}{D_{H}^{\text{fid}}(z)},\qquad
q_{\perp}(z,{\bf\Omega})=\frac{D_{A}^{\bf\Omega}(z)}{D_{A}^{\text{fid}}(z)}.$
(5.2)

These parameters can then be used to express the theory power spectrum in
terms of the observed wave-vectors and compare them to the data. Before doing
so, it is convenient to write the parallel and perpendicular components of the
wave-vectors as the $k$-modulus and the cosine of the line-of-sight variables,
$(k_{\parallel},\,k_{\perp})\rightarrow(k,\mu)$,

$\displaystyle
k^{\prime}=\frac{k}{q_{\perp}}\left[1+\mu^{2}\left(F^{-2}-1\right)\right]^{1/2},\qquad\mu^{\prime}=\frac{\mu}{F}\left[1+\mu^{2}\left(F^{-2}-1\right)\right]^{-1/2},$
(5.3)

where $F\equiv q_{\parallel}/q_{\perp}$ is the Alcock-Paczyński (AP)
parameter. With all this, we can write the theory power spectrum as a function
of the observed wave-vectors for each value of the explored $\bf\Omega$
parameters as,

$P_{\ell}(k)=\frac{2\ell+1}{2q_{\perp}^{2}q_{\parallel}}\int_{-1}^{1}d\mu\,P^{\text{th.}}_{s}\big{(}k^{\prime}(k,\mu),\mu^{\prime}(\mu)\big{)}\mathcal{L}_{\ell}(\mu),$
(5.4)

where $P^{\text{th.}}_{s}$ stands for any model used for the theoretical power
spectrum in redshift space, in this case, the one given by eq. (3.2). The pre-
factor $1/({q_{\perp}^{2}q_{\parallel}})$ accounts for the isotropic change in
volume when normalizing the true volume, given by the $k^{\prime}$’s, in terms
of the observed volume, given by the $k$’s.

###  5.2 Standard compression analysis

The Standard Compression technique focuses on extracting cosmological
information from two features imprinted in the $P(k)$ data vector: the BAO and
the RSD. The BAO information is contained in the oscillatory feature of the
power spectrum, $P_{w}$. Ideally, to extract that information one just needs
to measure the position of the BAO peak in the data power spectrum relative to
some known BAO pattern used as a reference, $P_{w}^{\rm ref}$,

$P_{w}(k)=P_{w}^{\rm ref}(k/s),$ (5.5)

where $P_{w}^{\rm ref}$ is the fixed-template at the reference cosmology for
the oscillatory component of the power spectrum. In this case, $s$ is just the
ratio of sound horizon scales at the drag epoch, of the cosmology of the data
and the reference cosmology, $s=r_{d}/r_{d}^{\rm ref}$.141414Note that if
$P^{\rm ref}(k)$ is expressed in units of $\,h^{-1}\text{Mpc}$ then $r_{d}$
must also be expressed in the same units in order to account for the shift in
$k$ due to the ratio between $h/h^{\rm ref}$. Note that the arbitrary choice
of this reference cosmology should not impact the determination of $r_{d}$:
choosing a template whose cosmology has a low/high $r_{d}^{\rm ref}$ is
compensated by obtaining a high/low best-fit value of $s$. Only a mismatch in
the shape and amplitude of the wiggles (not their position) between data and
template could in principle produce a bad fit to the data and a systematic
shift on $s$, although this could be mitigated by adding extra parameters that
enable the template to be modified for those features. We will come back to
this question later in §5.6.

As for the Full-Modelling case, the data power spectrum is expressed as a
function of the observed wave-vector and we must account for this in the
theory model. Unlike for the Full-Modelling approach, the parameters that
describe the scaling between true and observed wave-vectors, $q_{\parallel}$
and $q_{\perp}$ are free parameters of the model in the compressed analysis
approach, and therefore their best-fit values represent the ratio between the
distances expressed in the true underlying cosmology and the fiducial
cosmology of the catalogue,151515Note that since the comoving distances
physical catalogues are given in units of $\,h^{-1}\text{Mpc}$ one does not
need to assume any fiducial value for $h^{\rm fid}$ to transform from data to
physical catalogue. Consequently, the inferred $q$’s are independent of
$h^{\rm fid}$ as the $h$ value in both $D_{H,A}$ true and fiducial distances
is the same and cancel out. Thus, $H^{\rm fid}$ stands for $100h^{\rm
data}E({\bf\Omega}^{\rm fid})\,\text{km/s/Mpc}$.

$\displaystyle q_{\parallel}(z)=\frac{D_{H}(z)}{D_{H}^{\text{fid}}(z)},\qquad
q_{\perp}(z)=\frac{D_{A}(z)}{D_{A}^{\text{fid}}(z)}.$ (5.6)

At this point, we note that the $q$’s re-scale the full theory power spectrum,
and $s$ re-scales the isotropic wiggle power spectrum such that,

$P(k,\mu)=P^{\rm ref}_{\rm nw}(k^{\prime}(k,\mu),\mu^{\prime}(\mu))+P^{\rm
ref}_{\rm w}(k^{\prime}(k,\mu)/s,\mu^{\prime}(\mu))$ (5.7)

In this case the amplitude of the $P^{\rm ref}_{\rm nw}$ is parametrized
through another two parameters $f(z)$ and $\sigma_{8}(z)$. In particular,
$f(z)$ accounts for the effect of RSD, whereas $\sigma_{8}^{2}(z)$ is related
to the amplitude of the linear power spectrum at a given redshift. The linear
Kaiser terms will only depend on $f(z)\sigma_{8}(z)$, whereas higher-order
loop terms (see eqs. (3.9)–(3.15)) will display a combined dependence of
$f(z)$ and $\sigma_{8}(z)$ at different powers. Certainly, the $P(k)$ function
of eq. (3.8) displays a very strong correlation between $f(z)$ and
$\sigma_{8}(z)$ but a weak constraining power on those parameters
individually. In this scenario is convenient to consider the re-
parametrization
$\\{f(z),\,\sigma_{8}(z)\\}\rightarrow\\{f(z)\sigma_{8}(z),\,\sigma_{8}(z)\\}$.
In this parametrization, $\sigma_{8}(z)$ shows an extremely weak constraining
power, and fixing it to a fiducial value does not impact in a significant way
the results. We highlight that since $f(z)\sigma_{8}(z)$ does vary within the
analysis, fixing $\sigma_{8}(z)$ only fixes effectively the part of
$\sigma_{8}(z)$ beyond a power of the Kaiser term:
$f^{m}\sigma_{8}^{n}\rightarrow(f\sigma_{8})^{m}(\sigma_{8})_{\rm
fixed}^{n-m}$.

In practice since $P^{\rm ref}_{\rm nw}$ is very close to a power-law, $P_{\rm
nw}(k)\sim k^{-n}$, a shift in the amplitude of $P_{\rm nw}$ is perfectly
compensated by a shift in the argument. This implies that the isotropic shift
from $k^{\prime}\rightarrow k$ given by eq. (5.3) is indistinguishable from a
shift of $k\rightarrow k/s$. Although in reality, $P_{\rm nw}(k)$ is not a
power law, it has been shown that almost all the constraining power on the
$q$’s comes from the wiggle part of the power spectrum. Thus the following
approximation holds,

$P(k,\mu)\simeq P^{\rm ref}(k^{\prime}(k,\mu)/s,\mu^{\prime}(\mu)).$ (5.8)

For simplicity, one could choose the reference cosmology of the template to be
equal to the fiducial cosmology used to generate the physical catalogue. In
that case, the above expression consists of a re-scaling of $k_{\parallel}$
and $k_{\perp}$ (or $k$ and $\mu$), by,

$k^{\prime}_{\parallel}=k_{\parallel}/\alpha_{\parallel},\quad\text{and}\quad
k^{\prime}_{\perp}=k_{\perp}/\alpha_{\perp},$ (5.9)

with the scaling parameters given by,161616The eq. (5.10) is valid when the
reference cosmology of the template coincides with the fiducial cosmology used
to convert the redshift-to-distance when creating the physical catalogue. When
the reference and fiducial cosmologies differ, the sound horizon scale at the
drag epoch should be expressed in terms of the reference cosmology,
$r^{\text{ref}}_{d}$ .

$\displaystyle\alpha_{\parallel}=\frac{D_{H}(z)/r_{d}}{D_{H}^{\text{fid}}(z)/r^{\text{fid}}_{d}},\qquad\alpha_{\perp}=\frac{D_{A}(z)/r_{d}}{D_{A}^{\text{fid}}(z)/r^{\text{fid}}_{d}}.$
(5.10)

From these expressions, one can analogously write the expressions of eq. (5.3)
just changing $q_{\parallel,\,\perp}$ by $\alpha_{\parallel,\,\perp}$. With
all this, we write the theory power spectrum in terms of the observed wave-
vectors as,

$P_{\ell}(k)=\frac{2\ell+1}{2\alpha_{\perp}^{2}\alpha_{\parallel}}\int_{-1}^{1}d\mu\hskip
2.84544ptP^{\text{th.}}_{s}\big{(}k^{\prime}(k,\mu),\mu^{\prime}(\mu)\big{)}\mathcal{L}_{\ell}(\mu),$
(5.11)

where $P^{\text{th.}}_{s}$ is computed from the fixed linear power spectrum
template at the cosmology of reference (chosen to be equal to the fiducial
one) which includes $f(z)\sigma_{8}(z)$ as a free parameter. Since this
template does not change in each iteration (only the parameters
$\\{\alpha_{\parallel},\alpha_{\perp},f\sigma_{8}\\}$) the loop integrals of
eqs. (3.9)–(3.10) can be pre-computed which boosts the efficiency of this
approach. Finally, it is worth noting that the above expression re-scales the
volume of the power spectrum as $1/(\alpha_{\parallel}\alpha_{\perp}^{2})$.
This term is not fully correct, as the BAO shift (represented by the $s$
parameter) does not change the scale of the observed wave-modes. Therefore
this pre-factor should in fact read as, $(r_{d}^{\rm
fid}/r_{d})/(\alpha_{\parallel}\alpha_{\perp}^{2})$. However, as far as the
fiducial choice of the template has a $r_{d}^{\rm fid}$ value reasonably close
to the actual $r_{d}$ this mismatch in volume will be re-absorbed by other
bias or nuisance parameters.

At this point, it is worth mentioning that the inferred $\sigma_{8}$ parameter
(as part of $f\sigma_{8}$) corresponds to the fluctuations of the linear power
spectrum smoothed by a top-hat function at a scale of $8\,{\rm Mpc}/h^{\rm
ref}$. It is then the template which provides a ‘standard amplitude’ to
measure the data, in a similar fashion to the standard ruler $r_{d}$ for the
BAO peak measurements. In order to make the notation clearer we follow [23]
(see eq. 2.9) and refer to $\sigma_{s8}$ to the measured amplitude at the
fiducial cosmology of the template and reserve the $\sigma_{8}$ to denote the
usual amplitude fluctuations smoothed by $8\,{\rm Mpc}h^{-1}$ in the cosmology
of the data. Later in §5.4 we will discuss the actual relation between
$\sigma_{8}$ and $\sigma_{s8}$ in the light of ShapeFit.

In summary, the standard compression analysis employs a fixed template to
compress the BAO and RSD information into the parameters
$\\{\alpha_{\parallel},\alpha_{\perp},f\sigma_{8}\\}$, allowing us to
drastically reduce the computational time during the parameter inference
process because the linear power spectrum and the non-linear corrections
(which are the most expensive pieces) are computed only once for the reference
cosmology. This is because, in this analysis, the cosmological parameters are
fixed to the reference cosmology, while the compression parameters are
constrained.

The constraints on the compressed parameters could be interpreted in terms of
the cosmological parameters. This procedure is presented in section §5.4 and
essentially consists of using eqs. (5.10) in combination with
$\Omega_{m}(z)=f(z)^{\gamma}$ to obtain constraints on $\Lambda$CDM-like
parameters.

###  5.3 ShapeFit compression analysis

Standard compression analysis has proven to be a powerful tool for extracting
and analyzing information from BOSS and eBOSS surveys [104, 105, 106, 107,
108]. It perfectly compresses the most relevant information contained by LSS
galaxy catalogues, which is the BAO as the standard ruler and the RSD as a
probe of gravity. However, it completely ignores the information contained by
the transfer function. This limitation is historically mitigated by combining
the LSS galaxy catalogue with CMB anisotropic experiments, such as WMAP or
Planck. The latter are not sensitive to late-time effects such as the galaxy
BAO or RSD but are extremely good at constraining the shape of the linear
matter power spectrum. However, when the LSS analysis is considered alone (or
with just a few CMB priors) the Standard Compression is not lossless, and the
constraints significantly differ from those retrieved by the Full-Modelling
analysis. This limitation has been recently addressed by [23], in which the
authors propose a two-parameter extension to the standard compression that
improves the constraining power of the standard compression approach. In [23,
109] the authors show that for the $\Lambda$CDM model and for a volume of few
$\,h^{-3}\text{Gpc}^{3}$ this compression is close to being lossless. In
addition, under a tight prior on $\omega_{b}$ the compression has shown also
to be lossless for $\Lambda$CDM for volumes of few hundred
$\,h^{-3}\text{Gpc}^{3}$ [102]. In §7, this methodology is tested for several
cases, such as relaxing the prior on $\omega_{b}$ and $n_{s}$, as well as for
models beyond $\Lambda$CDM.

ShapeFit relies on the same main idea as the standard compression, consisting
of taking a linear power spectrum template evaluated at a reference cosmology
and introducing a set of variables that modify the form of that power spectrum
template. In addition to the BAO-wiggle modification (parametrized by
$\\{\alpha_{\parallel},\,\alpha_{\perp}\\}$), and the anisotropic amplitude
modification, ShapeFit introduces a novel modification of the broadband shape
of the linear power pivoting at a certain scale. Such modification aims to
capture information both from matter-radiation equally epoch, as well as the
spectral index of the primordial power spectrum [23],

$P^{\prime}_{L}(k)=P^{\text{ref}}_{L}(k)\exp\left\\{\frac{m}{a}\tanh\left[a\ln\left(\frac{k}{k_{p}}\right)\right]+n\ln\left(\frac{k}{k_{p}}\right)\right\\},$
(5.12)

where $a$ and $k_{p}$ are fixed to the values $a=0.6$ and
$k_{p}=0.03\,h^{\text{ref}}\,\text{Mpc}^{-1}\approx\pi/r^{\text{ref}}_{d}$
(although other values could also be chosen). The hyperbolic tangent is chosen
as a generic sigmoid function between two regimes where $k\rightarrow 0$ and
$k\rightarrow\infty$, as it was shown in figure 4 of [23]. The parameters $m$
and $n$ are freely varied. The $m$ parameter (or Shape) is the maxim slope at
the pivot scale $k_{p}$ and has a physical relation with the matter-radiation
equality epoch. Under certain assumptions for a $\Lambda$CDM model, this
parameter is related to $\Omega_{m}h^{2}$. Conversely, the $n$ parameter
controls the power-law index of the power spectrum at large $k$ values and is
related to the primordial spectral index parameter $n_{s}$. As we will discuss
later, in §7.3, both $m$ and $n$ present an extremely high anticorrelation,
such that one of them can be fixed (for e.g., $n=0$) and proceed just with one
ShapeFit parameter $m$. In the cosmology interpretation step (when one goes
from $\bf p$ variables to $\bf\Omega$) one takes $n_{s}$ to be fixed (at the
value of the chosen reference template, $P^{\rm ref}$, or alternative,
interpret $m$ as if it were $m-n$ being varied and $m+n$ being fixed to zero).

Hence, ShapeFit allows an extension of the Standard Compressed set of
parameters which allows access to the BAO and RSD information, but also to
early-time physics features through the $m$ and $n$ parameters. In its full
extension, the ShapeFit compression set of parameters consists of
$\\{\alpha_{\parallel},\alpha_{\perp},f\sigma_{s8},m,n\\}$ per redshift bin,
although in many cases $n$ will be kept fixed to zero. In Appendix B we
discuss some details on how the ShapeFit template variation of eq. (5.12) is
included in Folps.

| Compressed Fits | Direct-Fit
---|---|---
Analysis aspects |  Standard |  ShapeFit |  Full-Modelling
Source of information |  BAO and RSD |  BAO, RSD and slope of $P(k)$ |  Full-shape of $P(k)$
Parameters being varied |  $\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{8}\\}$ |  $\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{s8},m,n\\}$ |  Model parameters e.g., $\\{h,\Omega_{m},A_{s}\\}$
Linear $\vec{P_{L}(k)}$ |  Fixed (only amplitude variation through $\sigma_{8}$). |  Changes via eq. (5.12) |  Varies according to the model parameters
Non-linear corrections |  Computed only once at $P_{L}^{\rm ref}$ |  Computed via (B.1) |  Computed using the $P_{L}$ of the model
Scaling parameters |  Free parameters ($\alpha_{\parallel},\,\alpha_{\perp}$) |  Free parameters ($\alpha_{\parallel},\,\alpha_{\perp}$) |  Derived parameters ($q_{\parallel},q_{\perp}$) calculated via eq. (5.2)

Table 1: Summary of methodological aspects of the Standard, ShapeFit, and
Full-Modelling fits.

###  5.4 Interpreting compressed parameters in terms of cosmological
parameters

We have already presented the compressed step for both the Standard and
ShapeFit analyses, in which clustering information is condensed into the
parameters
$\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{s8},m,\,n\\}$.171717The
parameters $m$ and $n$ enter in the ShapeFit analysis via the eq. (5.12). They
are not considered in the Standard analysis. To interpret these results in
terms of cosmological parameters, we need to perform an additional inference
step. In this step, the parameters to be fitted will be the cosmological
parameters of a given model, typically $\\{h,\omega_{b},\omega_{cdm},A_{s}\\}$
for the $\Lambda$CDM model. Note, however, that the assumption of the
cosmological model is done a posteriori, i.e. after the multipole data has
been fitted during the compression step. Therefore, this inference step does
not depend on the perturbative model. Furthermore, the same set of compressed
data-vector parameters can be mapped into several models within the same
family of models (e.g., $\Lambda$CDM, $w\Lambda$CDM, $k\Lambda$CDM, etc.)
without needing to re-fit the data $P_{\ell}(k)$. In this sense, the
compressed set of parameters,
$\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{s8},m,\,n\\}$, are model-
agnostic, as they are valid for a wide family of cosmological models.

To accomplish the interpretation in terms of cosmological parameters, we need
to extract the mean and covariance matrix from the compressed parameters
obtained earlier. These extracted values will be used as the “data-vector” and
the corresponding covariance matrix, which we denote as $\bf{p}_{\text{mean}}$
and $\textit{Cov}_{\textbf{p}}$, respectively. After selecting our
cosmological model, during the inference process, we calculate, at each step,
the expected compressed parameters
$\mathbf{p}_{\bm{\Omega}}=\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{s8},m,\,n\\}_{\text{exp}}$
as a function of the cosmological parameters
$\bm{\Omega}=\\{h,\omega_{cdm},A_{s},...\\}$ of the chosen model.
Subsequently, we extract the posterior distributions of the cosmological
parameters. At this stage, we assume that the compressed parameters follow a
Gaussian distribution, and the likelihood function $L$ can be expressed as,

$L=\exp{\left[\left(\bf{p}_{\bm{\Omega}}-\bf{p}_{\text{mean}}\right)\textit{Cov}^{-1}_{\textbf{p}}\left(\bf{p}_{\bm{\Omega}}-\bf{p}_{\text{mean}}\right)^{T}\right]}.$
(5.13)

The interpretation of the scaling parameters $\alpha_{\perp,\parallel}$ is the
same for both Standard and ShapeFit analyses, and is given by eq. (5.10).
However, the interpretation of $f\sigma_{8}$ is slightly different in ShapeFit
since the parameter $\sigma_{8}$ becomes a function of the compressed
parameter $m$ according to eq. (5.12). We refer to this parameter as
$f\sigma_{s8}$, and the interpretation is provided by181818Note that equations
(5.14)-(5.17) include the ‘ref’ label, which indicates the reference template
cosmology. This cosmology may differ from the fiducial cosmology used in
constructing the physical catalogue. The implications of using different
template cosmologies will be further explored in §5.6. However, the results
presented in sections §6 and §7 are based on fixing the template cosmology to
the fiducial cosmology of simulations. [23, 110]

$f\sigma_{s8}=f\sigma^{\rm
ref}_{s8}A^{1/2}\times\exp{\left\\{\frac{m}{2a}\tanh{\left[a\,\ln{\left(\frac{r^{\text{ref}}_{d}\,[\,h^{-1}\text{Mpc}]}{8\,\,h^{-1}\text{Mpc}}\right)}\right]}\right\\}},$
(5.14)

with $A=A_{sp}/A_{sp}^{\rm ref}$, where $A_{sp}$ represents the fluctuation
amplitude given by

$A_{sp}=\frac{1}{s^{3}}P_{L,\text{nw}}\left(\frac{k_{p}}{s},\vec{\Omega}\right),\,\,\,s=\frac{r_{\text{d}}}{r^{\text{ref}}_{\text{d}}},$
(5.15)

where $a=0.6$ and $k_{p}=0.03\,h^{\text{ref}}\,\text{Mpc}^{-1}$. Meanwhile,
the effective ShapeFit parameter $m$ that controls the slope of the power
spectrum can be interpreted in terms of cosmological parameters as [23]

$m=\left.\frac{d}{d\ln
k}\left(\ln\left[\frac{P_{L,\text{nw}}\left(\frac{k_{p}}{s},\vec{\Omega}\right)/\mathcal{P}_{\mathcal{R}}(k_{p},\vec{\Omega})}{P_{L,\text{nw}}\left(\frac{k_{p}}{s},\vec{\Omega}^{\text{ref}}\right)/\mathcal{P}_{\mathcal{R}}(k_{p},\vec{\Omega}^{\text{ref}})}\right]\right)\right|_{k=k_{p}},$
(5.16)

with $\mathcal{P}_{\mathcal{R}}(k)=A_{s}\left(k/k_{p}\right)^{n_{s}-1}$ the
primordial power spectrum, which is given in terms of the global amplitude
$A_{s}$ and the spectral index $n_{s}$. For most of the tests, we set the
ShapeFit parameter $n$ to zero during the compression step due to its strong
degeneracy with $m$, as presented in §7.3. In this case, when interpreting the
results in terms of the cosmological parameters, we should keep the spectral
index fixed to its reference value $n^{\text{ref}}_{s}$. When varying the
power spectrum spectral index in the analysis, the shaft variable $n$ is
related to $n_{s}$ by

$n=n_{s}-n^{\text{ref}}_{s}.$ (5.17)

###  5.5 Comparison of Standard, ShapeFit and Full-Modelling approaches

Here we compare the performance of the three fitting methods discussed
earlier. To do this, in figure 1 we show the fittings to synthetic
AbacusSummit LRG mocks with a total volume of $200\,h^{-3}\text{Gpc}^{3}$. In
the left panel, we vary the compressed parameters for both Standard and
ShapeFit compression analyses, observing quite similar contours and
marginalized posterior distributions for the common parameters. Both
approaches successfully recover the expected parameters within 68% limits,
corresponding to $\\{\alpha_{\perp,\parallel}=1.0$, $f\sigma_{s8}=0.4501$,
$m=0.0\\}$. For this test, we set the reference template cosmology to the true
value of the simulation. The systematic effects introduced by the choice of
the reference cosmology are addressed below in §5.6.

The right panel of figure 1 displays the result when interpreting the
compressed parameters in terms of the $\Lambda$CDM cosmology. We also include
the results obtained with the direct Full-Modelling analysis. For all three
methods, we vary the cosmological parameters $\omega_{cdm}$, $h$ and
$\ln(10^{10}A_{s})$ using flat priors, and $\omega_{b}$ with a very
restrictive Gaussian prior around the value on the simulations, while the
spectral index $n_{s}$ is kept fixed. We notice the performance for this
setting is quite similar for the cases of Full-Modelling and ShapeFit,
outperforming the Standard approach. The Standard analysis has limited power
compared to other methods because it only relies on information from the BAO
peak position and RSD. On the other hand, ShapeFit provides additional
information from the broadband slope power spectrum, while Full-Modelling
extracts information from the full-shape of the spectra. For this reason,
throughout the rest of the paper, we will only focus on the ShapeFit and Full-
Modelling methods.

Figure 1:  Comparison of the Standard (gray), ShapeFit (orange), and Full-
Modelling (blue) posteriors when fitting the monopole and quadrupole of the
LRG AbacusSummit mocks. The left panel shows 2D contours of compressed
parameters for the Standard and ShapeFit approaches. Meanwhile, the right
panel compares the cosmological parameters for compressed and direct-fits. We
use the range of $k=[0.02,0.18]\,\,h^{-1}\text{Mpc}$ and the maximum available
volume by the simulations, $V_{25}=200\,h^{-3}\text{Gpc}^{3}$ . Black-dashed
lines represent the expected values.

###  5.6 Effect of the reference template in ShapeFit

ShapeFit builds upon the same fundamental concept as standard compressed
analysis, which involves defining a linear power spectrum fixed template that
is evaluated at a reference cosmology. This reference cosmology is usually
assumed to be close to the values reported by previous observations. However,
in realistic analyses, the assumed reference cosmology may differ from the
actual cosmology of the Universe. This can lead to differences between the
observed spectra from the data and those obtained from the fixed template. To
address this issue, ShapeFit introduces a set of parameters that modify the
power spectrum template, enabling it to capture cosmological information of
the clustering, as outlined in §5. At this point, we ask ourselves: What
happens if the reference template cosmology is very different from the true
cosmology? What impact does it have on the cosmological parameters?

Fits on LRGs with $V_{1}=8\,h^{-3}\text{Gpc}^{3}$

Figure 2:  Influence on the compressed parameters (top panel) and on the
cosmological parameters (bottom panel) caused by the choice of the reference
cosmology used to create the fixed template for the ShapeFit analysis. We
compare the results when the reference cosmology corresponds to the true
cosmology of the simulations (solid orange line) and when the cosmological
parameters $h$, $\omega_{\text{cdm}}$, and $A_{s}$ are biased with respect to
the true simulation values by a factor of +3% (dash-dotted purple line) and
+10 % (solid blue line), simultaneously in $h$, $\omega_{cdm}$ and $A_{s}$. We
use the range of $k=[0.02,0.18]\,\,h^{-1}\text{Mpc}$ and fit it to the LRG
cubic box with $V_{1}=8\,h^{-3}\text{Gpc}^{3}$. The black dashed lines
represent the expected values for each parameter.

Here we explore the impact of varying the reference template cosmology on
ShapeFit by changing it relative to the known true cosmology of the
simulations. To accomplish this, we utilize the LRG cubic mocks with a volume
of $V=8\,h^{-3}\text{Gpc}^{3}$ and conduct fits by incrementing the parameters
$h$, $w_{cdm}$, and $A_{s}$ by +3% and +10% relative to the true cosmology.
The differences in terms of the compressed (cosmological) parameters are
presented in the top (bottom) panel of figure 2, where we define
$\Delta\alpha_{\perp,\parallel}=\alpha_{\perp,\parallel}^{\text{obs}}-\alpha_{\perp,\parallel}^{\text{exp}}$,
with the labels ‘obs’ and ‘exp’ referring to the observed and expected values
of the parameters, respectively. Similarly, these definitions apply to the
other parameters.

From the figure, we can infer that as the reference template deviates further
from the true cosmology, the bias in the results becomes more pronounced.
These differences are primarily driven by the ShapeFit parameter $m$, which
significantly affects $w_{cdm}$ since both parameters control the slope of the
power spectrum. In contrast, the differences for $h$ and $A_{s}$ are minimal
(even for the case of the template with +10% shifts), as expected from the
compressed $\Delta\alpha^{\prime}s$ and $\Delta
f\sigma^{\text{ref}}_{s8}A^{1/2}$. This suggests that the BAO peak position
and RSD are minimally impacted by the reference cosmology.

The main motivation for using the LRG tracer and volume
$V=8\,\,h^{-3}\text{Gpc}^{3}$ is that, as will be shown in section 6, the LRG
tracer produces the smallest systematic errors, allowing for a clearer
visualization of the effect introduced purely by the choice of the reference
template cosmology. Additionally, $V=8\,\,h^{-3}\text{Gpc}^{3}$ represents the
minimum volume allowed by the simulations, so it is expected that for larger
volumes and when combining different tracers and redshift bins, the systematic
error introduced by the choice of the reference cosmology will be more
pronounced. This allows us to determine how relevant this systematic error in
ShapeFit will be for future observations.

It is important to note that the systematic effect shown in figure 2 is solely
associated with the choice of the reference template cosmology in ShapeFit.
This systematic effect is not present in Full-Modelling because it iteratively
varies the template as a function of the cosmological parameters. However,
both ShapeFit and Full-Modelling are affected by an additional systematic
introduced by the fiducial cosmology used to convert from redshift to
distances when generating the physical catalogue. However, the latter effect
is beyond the scope of this work and will be examined in a separate study
[111, 112].

##  6 Baseline analysis

In this section, we establish our baseline configurations, including the range
of wave-numbers utilized, chosen parameters, and priors. We perform MCMC runs
and present the findings of these baseline analyses employing a simplified
$\Lambda$CDM model. Subsequently, in §7, we explore models that extend beyond
these established settings.

Fits with: $V_{25}=200\,h^{-3}\text{Gpc}^{3};\,\ell=0,2;\,\text{Min. F.}$

Figure 3:  Dependence on the maximal wave-number $k_{\text{max}}$ for ShapeFit (left panel) and the direct Full-Modelling (right panel) analyses. We fix $k_{\text{min}}=0.02\,h^{-1}\text{Mpc}$ and consider the monopole and quadrupole with the full covariance of $V=200\,h^{-3}\text{Gpc}^{3}$ for LRGs ($z=0.8$), ELGs ($z=1.1$), and QSOs ($z=1.4$). We set the non-local biases $b_{s^{2}}$ and $b_{3\rm nl}$ to their coevolution prediction given by eqs. (6.1); however, this condition will be later relaxed in §6.2. The dotted lines represent the expected values for the simulations, while the grey shadow indicates the optimal maximum wave-number choice for our analysis.  Parameter type | Range | Parameter type | Range
---|---|---|---
Compressed |  | Cosmological |
$\alpha_{\perp}$ | $\mathcal{U}[0.8,1.4]$ | $h$ | $\mathcal{U}[0.5,0.9]$
$\alpha_{\parallel}$ | $\mathcal{U}[0.8,1.4]$ | $\omega_{cdm}$ | $\mathcal{U}[0.05,0.2]$
$f\sigma_{s8}$ | $\mathcal{U}[0.0,1.0]$ | $\omega_{b}$ | $\mathcal{N}(0.02237,0.00037)$
$m$ | $\mathcal{U}[-3.0,3.0]$ | $\ln(10^{10}A_{s})$ | $\mathcal{U}[2,4]$
Nuisances Min. F. |  | Nuisances Max. F. |
$b_{1}$ | $\mathcal{U}[10^{-5},10]$ | $b_{1}$ | $\mathcal{U}[10^{-5},10]$
$b_{2}$ | $\mathcal{U}[-50,50]$ | $b_{2}$ | $\mathcal{U}[-50,50]$
$b_{s^{2}}$ | Fixed to coev. | $b_{s^{2}}$ | Uninformative*
$b_{3\rm nl}$ | Fixed to coev. | $b_{3\rm nl}$ | Uninformative*
Table 2: Priors on compressed (used for ShapeFit), cosmological (used for
Full-Modelling), and bias parameters for minimal and maximal freedom settings
(Min.F. and Max.F, respectively). We adopt also uninformative priors on EFT
counterterms $\alpha_{0}$, $\alpha_{2}$ and $\alpha_{4}$ (when the
hexadecapole is considered) and in stochastics $\alpha_{0}^{\rm shot}$ and
$\alpha_{2}^{\rm shot}$. Still, we marginalize analytically over these
parameters as explained in Appendix A. The neutrino mass is fixed to
$M_{\nu}=0.06\,\text{eV}$. The asterisks accompanying Uninformative indicate
that these priors are in practice Gaussian, but sufficiently wide so they do
not alter the results from truly uninformative-flat priors, as illustrated in
figure 6.

###  6.1 k-range

First, we study the dependence on the maximal wave-number for ShapeFit and
Full-Modelling analyses. In these tests, we keep the minimal wave-number fixed
at $k_{\text{min}}=0.02\,h\,\text{Mpc}^{-1}$, while varying the upper limit
across the values of $0.12\leq k\,[\,h\,\text{Mpc}^{-1}]\leq 0.22$. These
analyses involve fitting the monopole and quadrupole of our model to the LRG,
ELG, and QSO datasets, using the full re-scaled covariance, corresponding to a
volume of $V_{25}=200\,h^{-3}\text{Gpc}^{3}$. The results are presented in
figure 3, illustrating the compressed ShapeFit parameters on the left and the
cosmological parameters obtained through the direct Full-Modelling approach on
the right, where the color of the symbol stands for the type of tracer, and
indicated. The horizontal dotted lines represent the expected values. We
observe that QSO datasets tend to have more significant systematics and larger
error bars compared to the LRG and ELG samples. The larger error bars are due
to the low number density of objects in the QSO catalog, compared to the ELG
and LRGs, resulting in a shot-noise-dominated sample, with a smaller effective
volume, but also limiting the range at high-$k$.

From figure 3, we observe that at higher wave-numbers, we obtain tighter
constraints but also larger systematic biases in some parameters. This is
because small scales are primarily influenced by non-linear effects,
complicating the model description, and these effects become more pronounced
in a larger volume. Our tests show that for the range of scales
$k_{\text{max}}=0.16-0.20\,h\,\text{Mpc}^{-1}$, we are able to recover the
parameters around the 68% confidence limits, but overall,
$k_{\text{max}}=0.18\,h\,\text{Mpc}^{-1}$ seems to have a better performance
with smaller standard deviations and systematic errors. Therefore, this will
be employed as the optimal maximum wave-number for our analysis, and hereafter
our default choice unless otherwise stated.

Figure 4:  The impact on the ShapeFit compressed parameters when altering the
non-local bias setups on the LRGs AbacusSummit sample. The top panel and
bottom panel display the impact of relaxing the coevolution condition (navy)
on $b_{s^{2}}$ and $b_{3\rm nl}$, respectively. The options explored are to
leave the remaining bias non-local bias term with an uninformative prior or
with Gaussian priors with mean values derived from coevolution and standard
deviations $\sigma$ set as multiples of the coevolution values, as indicated.
We re-scale the covariance by the factor $1/25$, hence the effective volume of
the sample is $V=200\,h^{-3}\text{Gpc}^{3}$ and fit the data over the interval
$k=[0.02,0.18]\,h\,\text{Mpc}^{-1}$.

Figure 5:  The impact on the Full-Modelling parameters when altering the non-
local bias setups on the LRGs AbacusSummit. Same notation as in figure 4.

###  6.2 Priors on nuisance parameters

We now turn our attention to choosing the parameters to be varied and their
priors. We will use two different settings. The first one is the simplest
implementation which consists of assuming coevolution for the bias parameters
and letting flat, uninformative prior in the rest of the nuisances. We call
this choice Minimal Freedom (from now on Min.F.), and we already used it to
generate the results shown in figures 1 and 3. That is, we fix $b_{s^{2}}$ and
$b_{3\rm nl}$ to their values obtained from coevolution theory [113, 114, 82]

$\text{Min.F.:}\qquad b_{s^{2}}=-\frac{4}{7}(b_{1}-1)\quad\text{and}\quad
b_{3\rm nl}=\frac{32}{315}(b_{1}-1).$ (6.1)

We employ coevolution theory for describing the clustering of galaxies,
although strictly is only valid for halos fulfilling a local Lagrangian bias.
We consider as a reasonable approximation that both haloes and galaxies should
have similar properties in terms of the locality of their bias. We already
showed for the LRG and ELG AbacusSummit mocks in figure 3, that that was
indeed a very good approximation. The bottom left panel of table 2 summarizes
all the nuisance parameters in Min.F. and their priors.

We want to investigate how relaxing the coevolution assumption impacts the
constraints on the compressed and cosmological parameters. To this end, we
adopt different Gaussian priors
$\mathcal{N}\left(\text{coev},N\times\text{coev}\right)$ with the mean given
by the coevolution value in eq. (6.1) and with standard deviation equals to
$N$ times that same coevolution value. These prior choices can be interpreted
as giving a slight preference for the coevolution setting, but departing out
of it as $N$ grows. At some point, the standard deviation of the priors is
sufficiently large so they become uninformative.

Figure 6:  Posteriors corresponding to the LRG AbacusSummit mocks, for the
ShapeFit (left panel) and Full-Modelling (right panel) analysis for the
Maximal (red) and Minimal Freedom (blue) setups, as defined in table 2. We
additionally display the contours for the choice of completely uninformative
priors in gray demonstrating that it matches quite accurately with the Max. F.
setting.

In figures 4 and 5 we show 1-dimensional marginalized posteriors when applying
the above priors to $b_{s^{2}}$ (top panel) and to $b_{3\rm nl}$ (bottom
panel), for ShapeFit and Full-Modelling, respectively. We consider the case of
the LRG AbacusSummit mocks. This analysis reveals that the posteriors exhibit
high sensitivity to restrictive priors applied to the third-order bias
$b_{3\rm nl}$, whereas the discrepancies arising from $b_{s^{2}}$ are
comparatively less pronounced. Further, we show that for scaling factors
$N=12$ for $b_{s^{2}}$, and $N=6$ for $b_{3\rm nl}$, we recover the
uninformative priors with high precision. We then take this choice as our
Max.F. case for LRG (in table 2, we indicate this with an asterisk mark). We
can, instead, simply let vary non-local biases over flat priors with large
enough intervals. However our option of Max.F. have better convergence
properties when sampling the parameters with MCMC, while the difference of the
posteriors between this option and the true Max.F. (which has uninformative
priors) are negligible. This behaviour is illustrated in figure 6, where we
compare between Min.F., Max.F. and uninformative priors. This analysis is
repeated for QSO and ELG datasets, such that we can guarantee that the all our
chains attain convergence and our priors are indeed uninformative, and that
the choice of Gaussian priors is not affecting the posterior distribution. In
table 2, we summarize the parameters and their priors for Max.F.

Gaussian priors in counterterms and stochastic parameters

Figure 7:  The impact on ShapeFit compressed parameters (left panel) and
cosmological parameters obtained with the direct Full-Modelling approach
(right panel) when using Gaussian priors on the counterterm
$\alpha_{0},\alpha_{2}$ and stochastic parameters $\alpha_{0}^{\text{\rm
shot}},\alpha_{2}^{\text{\rm shot}}$. The Gaussian priors are centered at zero
with a variance that goes from $0.2\times$ up to $30\times$ the mean value of
that parameter given in table 9, represented by the orange, purple, blue and
green symbols. Alternatively, the black crosses represent the effect of
shifting the center of the Gaussian prior from zero to the mean values given
by table 9 and keeping the same variance as in the overlapping color symbol.
The navy star represents the effect of setting a flat and uninformative prior
in each of these parameters. The vertical dotted lines display the expected
value for each parameter. For all the cases $b_{s^{2}}$ and $b_{3\rm nl}$ are
kept to their coevolution prediction.

The effects of priors have been warned in [36, 115, 116]. The choice of priors
on nuisance parameters may affect the inference of cosmological parameters,
which would indicate that the posteriors are dominated by the priors and not
by the data. Prior weights, for example, can explain some of the discrepancies
found for different analyses of BOSS data. In this work our main results adopt
uninformative priors. However, for illustration, in figure 7 (left panel on
ShapeFit parameter; right panel on Full-Modelling parameters) we show the
impact of adopting more tight restrictions on our baseline settings, a
practice commonly adopted in numerous studies within this field. First, for
each counterterm $\alpha_{i}$, we choose a family of Gaussian priors all
centred with mean $\bar{\alpha}_{i}=0$ and standard deviations
$\sigma_{\alpha_{i}}=N\times\text{mean}$; that is, a factor $N$ times the mean
values found in the case of uninformative priors obtained from a previous
analysis. We repeat these experiments, but using Gaussian priors centered at
the mean values, $\bar{\alpha}_{i}=\text{means}$, not surprisingly, showing
even better constraints. This analysis shows that one can obtain tighter
constraints on cosmological parameters by choosing ad hoc priors. One can
argue that EFT counterterms are expected to operate above a non-linear scale
$k_{\text{NL}}$, such their contributions $c(k/k_{\text{NL}})^{2}P_{L}(k)$
should have a parameter $c$ of order unity. In this sense, it is reasonable to
impose a prior $c\sim\mathcal{N}(0,1)$. However, although the above argument
may sound compelling, the choice of the non-linear scale $k_{\text{NL}}$ is
arbitrary and the posterior distributions can be very sensitive to this value.
That is, the results may become highly prior dependent, which has now
translated to the specific chosen value for $k_{\text{NL}}$. For this reason,
in this work, we rely only on uninformative priors.

However, we foresee situations where imposing informative priors may be
necessary to obtain some results. This happens, for example, for the
$w_{0}w_{a}$CDM model where the equation of state, in addition to being
degenerate with $A_{s}$, becomes highly degenerate with the scale-independent
piece of the shot noise in the absence of external data. In this case, an
uninformative prior in $\alpha^{\text{\rm shot}}_{0}$ simply wipes out any
constriction on $w_{0}$ and $w_{a}$. To gain some information in this extreme
situation one may want to assume a physical (and informative) Gaussian prior
on the shot noise parameter centered around the Poissonian prediction with a
moderate variance. However, we must bear in mind that the information on
cosmology will be informed by the arbitrary choice of this variance.

###  6.3 Baseline analysis results

Having discussed the nuisance parameters and their priors, we are now in a
position to fix a baseline analysis. This serves as an anchor from which we
think analysis beyond it should be compared. The baseline analysis considers
the compressed and cosmological parameters,

SF:
$\displaystyle\quad\\{\alpha_{\parallel},\,\alpha_{\perp},\,f\sigma_{s8},\,m\\}.$
FM:
$\displaystyle\quad\\{\omega_{cdm},\,\omega_{b},\,h,\,\ln(10^{10}A_{s})\\}.$

Meanwhile, for the nuisance parameters, we use two different setups: 1)
Min.F., which contains all the parameters consistent with local Lagrangian
biases and uses the coevolution relations given by eqs.(6.1) for tidal and
non-local third-order biases. 2) Max.F., which contains all the parameters
included in the full one-loop expansion.

Min.F.:
$\displaystyle\quad\\{b_{1},\,b_{2},\,\alpha_{0},\,\alpha_{2},\,\alpha_{0}^{\text{\rm
shot}},\,\alpha_{2}^{\text{\rm shot}}\\}.$ Max.F.:
$\displaystyle\quad\\{b_{1},\,b_{2},\,b_{s^{2}},\,b_{3\rm
nl},\,\alpha_{0},\,\alpha_{2},\,\alpha_{0}^{\text{\rm
shot}},\,\alpha_{2}^{\text{\rm shot}}\\}.$

Furthermore, in the baseline analysis, all the parameters have uninformative
priors, except for $\omega_{b}$ for which a Gaussian prior is chosen, as
summarized in table 2.

In figure 8 we show a corner plot of the posteriors of our fittings to the
AbacusSummit mocks. The top panels correspond to the volume
$V_{5}=40\,\,h^{-3}\text{Gpc}^{3}$, while the bottom panels to the maximum
volume allowed by the simulations $V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$. The
left panels display the ShapeFit analysis once we have converted the
compressed to cosmological parameters, whereas the right panels display the
Full-Modelling fits. In all these cases we have adopted the Max.F. setup.
Table 3 display the corresponding numerical results on which we report the
mean and 68% confidence interval (c.i.) of the 1-dimensional marginalized
posterior. Additionally, this table also reports the results under the Min.F.
setup.

The way the different tracers are combined is the following. In Full-
Modelling, we sample cosmological parameters as usual and calculate power
spectrum multipoles using different sets of nuisance parameters for each
redshift bin. In ShapeFit, we compress parameters separately for each tracer
and then combine them during interpretation in terms of cosmological
parameters. In both cases, the likelihood is computed from the combined
datasets, following $L=L_{\rm LRG}+L_{\rm ELG}+...$, which does not assume
correlation across $z$-bins. In fact, we expect that the $z$-bins that arise
from the same initial conditions (same realization index) present a certain
correlation across them. However, since the redshift snapshots are non-
contiguous, and the bias of the tracers is very distinct, for simplicity we
approximate them to be uncorrelated. The only potential effect of doing so is
to slightly underestimate the true posteriors, which would result in enlarging
the systematic errors (hence our systematic budget would only be more
conservative). By comparing the offsets on the parameters on the individual
tracers, and in the combined sample (see figure 8 and Table 3) we do not
particularly see significant differences, suggesting that our $z$-independent
approximation is reasonable.

Figure 8: Baseline Analyses. Contour plots for the posterior distributions
computed with the ShapeFit (left) and Full-Modelling (right) analyses for the
volumes of $V_{5}=40\,\,h^{-3}\text{Gpc}^{3}$ (top) and
$V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$ (bottom), comparing the constraints on
the cosmological parameters of the LRG $(z=0.8)$ in orange, ELG $(z=1.1)$ in
green, QSO $(z=1.4)$ tracers in gray. We also plot the constraints of the
joint LRG + ELG + QSO analysis.

We have performed separate fittings to LRGs, ELGs, QSOs, and the joint
analysis includes the three tracers. We also performed the joint analysis
including only LRGs and ELGs and obtain results very similar the full, three
tracers, joint analyisis. To avoid cluttering, we do not plot these results.

Figure 9 displays the Full-Modelling and ShapeFit-inferred posteriors on the
$\Lambda$CDM parameters for better visualization. Both panels display the
results on the combined $\text{LRG}+\text{ELG}+\text{QSO}$ analysis for the
Min.F. and Max.F. setups, in the left and right panels, respectively.

From Table 3 and figure 9 we notice that in some configurations, the ShapeFit
analysis yields tighter constraints on some cosmological parameters. This is
counter-intuitive since by construction one expects the ShapeFit analysis to
be, at most, as constraining as the Full-Modelling one. Notice, however, that
the posterior for the shape parameter $m$ is in general not Gaussian, and not
even symmetric, as can be seen in the figure 4. In general, this is true for
all compressed parameters, but more evident for $m$. However, these
distributions are treated, for simplicity, as if they were Gaussian in the
conversion step to cosmological parameters. Therefore, this approximation in
the posterior of $m$ can yield some situations in which the ShapeFit analysis
seems more constraining than the Full-Modelling one. To avoid this, instead of
performing a Gaussian approximation, we could take the full posterior surface
of the compressed parameter distribution in the interpretation.

Figure 9:  Baseline Analyses. Contour plots for the posterior distributions of the cosmological parameters obtained with ShapeFit (SF) and Full-Modelling (FM) approaches for the joint LRG $(z=0.8)$ \+ ELG $(z=1.1)$ \+ QSO $(z=1.4)$. The left panel corresponds to the Min. F. setting, while the right panel corresponds to the Max. F. setting, both with a total volume of $V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$.  |  | LRG |  | ELG |  | QSO |  | LRG+ELG+QSO
---|---|---|---|---|---|---|---|---
|  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$ |  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$ |  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$ |  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$
FM,Min F.$V_{25}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6741\pm 0.0036$ | $0.15$ |  | $0.6744\pm 0.0035$ | $0.24$ |  | $0.6780\pm 0.0048$ | $0.92$ |  | $0.6738^{+0.0031}_{-0.0030}$ | $0.05$
$\omega_{cdm}$ |  | $0.1197^{+0.0018}_{-0.0015}$ | $0.19$ |  | $0.1198^{+0.0014}_{-0.0012}$ | $0.16$ |  | $0.1182^{+0.0022}_{-0.0017}$ | $0.88$ |  | $0.1196\pm 0.0010$ | $0.44$
$\ln(10^{10}A_{s})$ |  | $3.032^{+0.017}_{-0.022}$ | $0.20$ |  | $3.021^{+0.010}_{-0.015}$ | $1.07$ |  | $3.022^{+0.022}_{-0.032}$ | $0.48$ |  | $3.028^{+0.008}_{-0.010}$ | $0.87$
SF,Min F.$V_{25}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6747^{+0.0034}_{-0.0035}$ | $0.31$ |  | $0.6752\pm 0.0034$ | $0.46$ |  | $0.6779^{+0.0042}_{-0.0043}$ | $0.98$ |  | $0.6751\pm 0.0031$ | $0.49$
$\omega_{cdm}$ |  | $0.1197^{+0.0015}_{-0.0014}$ | $0.21$ |  | $0.1197\pm 0.0012$ | $0.24$ |  | $0.1185^{+0.0016}_{-0.0018}$ | $0.85$ |  | $0.1195\pm 0.0010$ | $0.45$
$\ln(10^{10}A_{s})$ |  | $3.038^{+0.020}_{-0.018}$ | $0.08$ |  | $3.030\pm 0.014$ | $0.45$ |  | $3.029^{+0.029}_{-0.028}$ | $0.26$ |  | $3.032\pm 0.011$ | $0.37$
FM,Max F.$V_{25}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6748^{+0.0040}_{-0.0043}$ | $0.28$ |  | $0.6747^{+0.0037}_{-0.0036}$ | $0.31$ |  | $0.6775^{+0.0050}_{-0.0049}$ | $0.78$ |  | $0.6744\pm 0.0033$ | $0.23$
$\omega_{cdm}$ |  | $0.1209^{+0.0024}_{-0.0028}$ | $0.34$ |  | $0.1213^{+0.0019}_{-0.0021}$ | $0.64$ |  | $0.1188^{+0.0026}_{-0.0032}$ | $0.40$ |  | $0.1205\pm 0.0015$ | $0.34$
$\ln(10^{10}A_{s})$ |  | $3.025^{+0.038}_{-0.035}$ | $0.31$ |  | $3.020\pm 0.025$ | $0.65$ |  | $3.032^{+0.044}_{-0.043}$ | $0.10$ |  | $3.030\pm 0.018$ | $0.34$
SF,Max F.$V_{25}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6745^{+0.0036}_{-0.0040}$ | $0.22$ |  | $0.6752^{+0.0035}_{-0.0036}$ | $0.44$ |  | $0.6775^{+0.0050}_{-0.0049}$ | $0.77$ |  | $0.6749\pm 0.0032$ | $0.40$
$\omega_{cdm}$ |  | $0.1203^{+0.0021}_{-0.0022}$ | $0.14$ |  | $0.1210^{+0.0018}_{-0.0019}$ | $0.54$ |  | $0.1192^{+0.0032}_{-0.0031}$ | $0.24$ |  | $0.1205\pm 0.0014$ | $0.34$
$\ln(10^{10}A_{s})$ |  | $3.041^{+0.028}_{-0.026}$ | $0.16$ |  | $3.037\pm 0.024$ | $0.04$ |  | $3.039^{+0.041}_{-0.040}$ | $0.06$ |  | $3.039\pm 0.017$ | $0.18$
FM,Min F.$V_{5}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6744^{+0.0050}_{-0.0049}$ | $0.17$ |  | $0.6753\pm 0.0048$ | $0.34$ |  | $0.6793^{+0.0087}_{-0.0085}$ | $0.66$ |  | $0.6750\pm 0.0039$ | $0.37$
$\omega_{cdm}$ |  | $0.1191^{+0.0032}_{-0.0028}$ | $0.31$ |  | $0.1193^{+0.0028}_{-0.0023}$ | $0.25$ |  | $0.1172^{+0.0039}_{-0.0034}$ | $0.74$ |  | $0.1193\pm 0.0016$ | $0.41$
$\ln(10^{10}A_{s})$ |  | $3.040^{+0.037}_{-0.048}$ | $0.08$ |  | $3.028^{+0.021}_{-0.034}$ | $0.26$ |  | $3.029^{+0.051}_{-0.057}$ | $0.13$ |  | $3.032^{+0.017}_{-0.021}$ | $0.20$
SF,Min F.$V_{5}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6739^{+0.0047}_{-0.0051}$ | $0.06$ |  | $0.6753\pm 0.0049$ | $0.34$ |  | $0.6779\pm 0.0085$ | $0.50$ |  | $0.6749^{+0.0040}_{-0.0034}$ | $0.35$
$\omega_{cdm}$ |  | $0.1191^{+0.0023}_{-0.0027}$ | $0.38$ |  | $0.1194\pm 0.0021$ | $0.30$ |  | $0.1180\pm 0.0031$ | $0.63$ |  | $0.1190^{+0.0016}_{-0.0015}$ | $0.66$
$\ln(10^{10}A_{s})$ |  | $3.048\pm 0.040$ | $0.30$ |  | $3.035\pm 0.026$ | $0.05$ |  | $3.042^{+0.054}_{-0.055}$ | $0.11$ |  | $3.042\pm 0.021$ | $0.25$
FM,Max F.$V_{5}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6753\pm 0.0056$ | $0.30$ |  | $0.6754^{+0.0048}_{-0.0050}$ | $0.36$ |  | $0.6787^{+0.0090}_{-0.0089}$ | $0.57$ |  | $0.6748\pm 0.0042$ | $0.29$
$\omega_{cdm}$ |  | $0.1212^{+0.0045}_{-0.0052}$ | $0.24$ |  | $0.1225^{+0.0038}_{-0.0044}$ | $0.59$ |  | $0.1201^{+0.0053}_{-0.0070}$ | $0.01$ |  | $0.1211^{+0.0028}_{-0.0027}$ | $0.40$
$\ln(10^{10}A_{s})$ |  | $3.023^{+0.073}_{-0.066}$ | $0.19$ |  | $3.011^{+0.055}_{-0.054}$ | $0.46$ |  | $3.008^{+0.094}_{-0.086}$ | $0.31$ |  | $3.024\pm 0.035$ | $0.34$
SF,Max F.$V_{5}$: |  |  |  |  |  |  |  |  |  |  |  |
$h$ |  | $0.6740\pm 0.0053$ | $0.07$ |  | $0.6749^{+0.0051}_{-0.0050}$ | $0.25$ |  | $0.6769\pm 0.0089$ | $0.37$ |  | $0.6748^{+0.0039}_{-0.0040}$ | $0.30$
$\omega_{cdm}$ |  | $0.1200^{+0.0033}_{-0.0037}$ | $0.01$ |  | $0.1222^{+0.0035}_{-0.0041}$ | $0.58$ |  | $0.1206^{+0.0056}_{-0.0062}$ | $0.10$ |  | $0.1208^{+0.0024}_{-0.0023}$ | $0.32$
$\ln(10^{10}A_{s})$ |  | $3.049^{+0.053}_{-0.049}$ | $0.25$ |  | $3.037^{+0.054}_{-0.045}$ | $0.02$ |  | $3.041^{+0.091}_{-0.081}$ | $0.05$ |  | $3.044^{+0.033}_{-0.034}$ | $0.22$
Table 4: Baseline Analyses. 1-dimensional posterior means and 68% confidence
interval (c.i.) for Full-Modelling (FM) and ShapeFit (SF), and volumes $V_{5}$
and $V_{25}$. We also show the deviation from the value of the simulations
$\Delta\bm{\Omega}/\sigma$.

Other differences can be seen in the profile of the posteriors, particularly
in those for $\ln(10^{10}A_{s})$. In the Full-Modelling fits, the posterior is
clearly asymmetric with pronounced right tails; see, e.g. the left panel of
figure 9. A similar behaviour is observed in [43], which uses Lagrangian PT
with velocileptors code, hinting that this is a property of the data itself.
However, this feature is completely lost in the ShapeFit analysis, where the
posteriors are highly symmetric.

Figure 10:  Baseline Analyses. Constraints on cosmological parameters obtained
from the compressed ShapeFit (in orange) and direct Full-Modelling (in blue)
analyses using both Min. F. (solid lines) and Max. F. (dashed lines) settings,
as presented in Table 3. The left panel displays the results based on a
covariance matrix scaled to match a volume of
$V_{5}=40\,\,h^{-3}\text{Gpc}^{3}$, while the right panel is scaled to
$V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$. We compare the constraints of the LRG
$(z=0.8)$, ELG $(z=1.1)$, QSO $(z=1.4)$ tracers, as well as the combined LRG +
ELG and ALL = LRG + ELG + QSO sets.

In figure 10 we present 1-dimensional plots comparing the ShapeFit and Full-
Modelling approaches for all considered tracers using both Max.F and Min.F.,
and volumes $V_{5}$ and $V_{25}$. The majority of configurations/parameters
show very similar constraints, except for a few of them. Differences in this
case are more evident in $\omega_{cdm}$ and $\ln(10^{10}A_{s})$ for LRGs,
where the constrictions given by ShapeFit are considerably tighter in the
Max.F setup.

As stated above, in Table 3, we present the posterior means along with their
corresponding 0.68 c.i. for the FM and SF analyses. Additionally, we assess
the accuracy of recovering the cosmological parameters as quantified by
$\Delta\Omega/\sigma_{\Omega}$. Here,
$\Delta\Omega=|\Omega_{i}^{\text{mean}}-\Omega_{i}^{\text{simulations}}|$, the
difference between the mean of the 1-dimensional marginalized posterior
distribution of parameter $\Omega_{i}$ and its value on the mocks. Meanwhile,
$\sigma_{\Omega}$ corresponds to the standard deviation of such a parameter.
In all our cases, all the parameters fall within $1\sigma$ of their true
values. Specifically, for the Max.F. parameter in the volume $V_{5}$, which is
the least statistically restrictive configuration we consider, all the
parameters are within $0.6\sigma$. Our least accurate estimations are observed
in QSO mocks, which constitute the data that exhibit the highest level of
scattering, in particular, we quote here the Min.F. case on which we find a
maximum deviation of $0.9\sigma$.

##  7 Beyond baseline analyisis

In this section, we explore some extensions to the baseline analysis: such as
including the hexadecapole (§7.1), relaxing the prior on the baryon density
(§7.2), letting the spectral index parameter-free (§7.3), allowing the total
neutrino mass to vary (§7.4), relaxing the flatness assumption on the spatial
curvature parameter (§7.5), and treating the equation of state of dark energy
as a free parameter of the model (§7.6). We maintain the same settings and
priors as in the baseline analysis, as defined in Table 2, and perform the
analysis on the LRG sample with the $V_{5}=40\,h^{-3}\text{Gpc}^{3}$ and
$V_{25}=200\,h^{-3}\text{Gpc}^{3}$ scaling-covariance choices.

The results are summarized in tables 5 and 7, for the $V_{25}$ and $V_{5}$
choices, respectively.

###  7.1 Including the hexadecapole

As an extension to section 6, we examine the effect of including the next non-
zero multipole in the expansion, the hexadecapole, which is also fitted up to
$k_{\text{max}}=0.18\,h{\rm Mpc}^{-1}$.

Figure 11: Including the Hexadecapole. Comparison of the posterior
distributions when fitting the monopole and quadrupole ($\ell=0,2$) and when
including the hexadecapole ($\ell=0,2,4$) in the analysis. The left panel
displays the constraints on the compressed ShapeFit parameters, while the
right panel compares the constraints obtained on the $\Lambda$CDM parameters
for the ShapeFit and Full-Modelling analyses. Both cases fit the LRG sample,
assume the Min.F. setup and scale the covariance matrix with the $V_{25}$
volume, as indicated.

In the left panel of figure 11, we compare the constraints on the compressed
ShapeFit parameters when fitting the monopole and quadrupole (orange lines)
and when including the hexadecapole (green lines) in the fit. With the
addition of the hexadecapole, the $\alpha$ parameters, and $f\sigma_{s8}$ tend
to shift slightly closer to their true values. Furthermore, the inclusion of
the hexadecapole results in a reduction of about $20\%$ in the
$\alpha_{\parallel}$ error bar. A similar behavior is presented in figure 7 of
[102]. In this analysis, which improves upon [102] by accounting for the
inhomogenous mode sampling in the hexadecapole (eqs. 40, 41) of [107]), we
find that $\alpha_{\parallel}$ is correctly recovered, confirming that the
shift seen in [102] was indeed due to this mode- sampling affect.

The right panel of figure 11 compares the constraints on the cosmological
parameters for ShapeFit (orange lines) and Full-Modelling (blue lines),
considering both the inclusion and exclusion of the hexadecapole. For
$\Lambda$CDM the effect of adding the hexadecapole is very mild ( $\lesssim
3\%$ reduction on the error bars, see table 5): the internal priors of the
3-parameters $\Lambda$CDM model make the anisotropic information carried by
the hexadecapole fully redundant with that of the quadrupole. However, as we
will see in §7.5 and §7.6, this might be not the case for cosmologies beyond
$\Lambda$CDM, where including the hexadecapole could improve the constraining
power because the internal model priors change.

###  7.2 Dependence on the $\omega_{b}$ prior

In the baseline analysis we have applied a Gaussian prior on the baryon
density centered on the true underlying value of the AbacusSummit simulations
and with a very narrow width given by BBN observations. In this section, we
study the impact of relaxing this assumption. We focus on two different priors
scenarios: one with a less restrictive Gaussian prior represented by
$\omega_{b}=\mathcal{N}(0.02237,0.001)$, and another with a uniform prior
given by $\omega_{b}\in\mathcal{U}(0.01,0.04)$. The results for the volume
$V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$ are presented in figure 12 and table 5,
while results with $V_{5}=40\,\,h^{-3}\text{Gpc}^{3}$ are summarized in table
7.

Figure 12 shows that for both ShapeFit (top panel) and Full-Modelling (bottom
panel) the constraining power on $h$ is reduced by approximately a factor of
$2$ when transitioning from the BBN to the less restrictive Gaussian prior
with width $\sigma=0.001$ (purple and red contours). In fact, the ability to
determine $h$ is directly related to the ability to infer the sound horizon
scale, $r_{d}$ (and thus calibrate the distances in our analysis). This
happens because the baryon density is a key ingredient to determining $r_{d}$
(given the early-time physics of the photon-baryon plasma, see eq (C.1)).
Thus, increasing the uncertainty on $\omega_{b}$, translates into increasing
the uncertainty on $r_{d}$ and consequently on $h$. In the case of the Full-
Modelling analysis, some information on $\omega_{b}$ is extracted from the
amplitude of the BAO wiggles, and therefore, the system is not completely
blind to $r_{d}$ in the absence of $\omega_{b}$ priors (as for the ShapeFit
analysis). This helps to diminish the uncalibration of the system and to
partially mitigate the degeneration between $h$ and $\omega_{b}$. This is
clearer when adopting a uniform prior (green and grey contours). For the Full-
Modelling analysis the constraining power on $h$ decreases by a factor of
$\sim 3.5$ with respect to the baseline analysis, while for ShapeFit, without
an informative prior on $\omega_{b}$ the analysis can’t inform at all on
$r_{d}$ values by construction, which causes almost flat posteriors on the
variables that require physical distances ($h$ and $\omega_{cdm}$). However,
$\Omega_{m}$ is still informed through the AP effect on the BAO peak position
(radial BAO relative to angular BAO), as the uncalibrated expansion history is
sensitive to the background variables [23].

This limitation of ShapeFit with respect to the Full-Modelling analysis has no
practical impact. In the case of the absence of a $\omega_{b}$ prior, the
Full-Modelling is not able to provide competitive results on $h$, and once a
BBN or CMB prior on $\omega_{b}$ is added, the extra information on $r_{d}$
from the amplitude of the BAO wiggles become irrelevant.

Figure 12: Dependence on the $\omega_{b}$ prior. Impact on the posterior
distributions of the cosmological parameters when relaxing the prior on the
baryon abundance $\omega_{b}$ with respect to the baseline analysis, where
$\omega_{b}=\mathcal{N}(0.02237,0.00037)$. We explore a weaker Gaussian prior
of $\omega_{b}=\mathcal{N}(0.02237,0.001)$, and a uniform prior
$\omega_{b}\in\mathcal{U}(0.01,0.04)$. On the top panel, are for the
compressed ShapeFit, and the bottom panel, for the direct Full-Modelling.

###  7.3 Impact of varying $n_{s}$

Figure 13: Impact of varying $n_{s}$. Comparison of the posterior
distributions between the baseline analysis and an extended analysis where the
spectral index $n_{s}$ is relaxed to be a free parameter with a flat prior
$n_{s}\in\mathcal{U}(-4,4)$. We present this comparison for the ShapeFit and
the Full-Modelling analyses. For ShapeFit, it was necessary to include the two
shape parameters $m$ and $n$ during the compression stage, following equation
(5.12). Subsequently, they were transformed into cosmological parameters using
equations (5.16) and (5.17). The top right box shows the strong
anticorrelation between $m$ and $n$.

We explore now the impact of varying the spectral index $n_{s}$ while keeping
the BBN prior on $\omega_{b}$ in both the ShapeFit and Full-Modelling
analyses. We compare this extension to our baseline results, where the
spectral index was set to the value from simulations. To do so, it was
necessary to recompute the pipeline from the beginning in both the compressed
ShapeFit and direct Full-Modelling analyses. For the ShapeFit, we include the
additional compressed parameter $n$, presented in eq. (5.12). Then, we convert
the set of compressed parameters
$\\{\alpha_{\perp},\alpha_{\parallel},f\sigma_{s8},m,n\\}$ into the
cosmological parameters {$h,\omega_{cdm},\omega_{b},A_{s},n_{s}$} by using the
prescription outlined in §5.4. On the other hand, for the Full-Modelling
analysis, we simply allow the spectral index to be a free parameter of the
model. In both analyses, we use uninformative prior on the spectral index,
given by $n_{s}\in\mathcal{U}(-4,4)$. The results are displayed in figure 13
and tables 5 and 7.

The parameters $n_{s}$ and $w_{cdm}$ are somewhat degenerate because they both
affect the slope of the power spectrum after the turnaround [117, 118].
Therefore, when $n_{s}$ is treated as a free parameter, a degeneracy with
$w_{cdm}$ emerges, as evidenced in figure 13. This degeneracy worsens the
constraint on $w_{cdm}$, which in turn, due to its degeneracy with $A_{s}$,
negatively impacts the measurement of $A_{s}$. We also note that $n_{s}$ has a
relatively small impact on $h$, which primarily depends on the $\omega_{b}$
prior, as shown in §7.2.

The upper right panel in figure 13 shows the strong anti-correlation between
the ShapeFit parameters $m$ and $n$. This can be seen by expanding eq. (5.12)
with $n=-m$ around $k=k_{p}$,

$\exp\left\\{\frac{m}{a}\tanh\left[a\ln\left(\frac{k}{k_{p}}\right)\right]-m\ln\left(\frac{k}{k_{p}}\right)\right\\}=1+\frac{1}{3}a^{2}m\left(1-\frac{k}{k_{p}}\right)^{3}+\cdots,$
(7.1)

showing a second-order degeneracy between $n$ and $m$. This degeneracy could
potentially be the reason behind the poor constraining power of $n_{s}$ in
ShapeFit. The degeneracy between $n$ and $m$ directly affects the measurements
of $n_{s}$ through eq. (5.17), and consequently, $\omega_{cdm}$ and $A_{s}$
are also affected. In Full-Modelling, this effect is less pronounced because
Full-Modelling extracts additional information that helps break the degeneracy
between $n_{s}$ and $w_{cdm}$.

###  7.4 Impact of varying $M_{\nu}$

Neutrinos are the most abundant massive elementary particles in our Universe
and play a fundamental role in nature, ranging from the smallest to the
largest scales in the Universe. We know from neutrino oscillation experiments
that they have non-zero masses [119, 120, 121]. These observations have led to
the identification of two different neutrino mass splitting patterns known as
the normal hierarchy, where the smaller mass split occurs between the first
and second lightest masses, and the inverted hierarchy, where the smaller mass
differences occur between the heaviest masses. While these experiments provide
insights into mass-squared differences, they cannot directly measure
individual masses or, equivalently, the neutrino mass scale. They however,
establish lower limits on the sum of the neutrino masses:
$M_{\nu}=\sum_{i}m_{\nu_{i}}>0.06\,\text{eV}$ for the normal hierarchy and
$M_{\nu}>0.1\,\text{eV}$ for the inverted hierarchy. Nevertheless, alternative
methods are needed to determine the total neutrino mass and distinguish
between their mass hierarchies.

The LSS offers an alternative approach to placing constraints on the sum of
neutrino masses. This is achieved through indirect measurements of the
features that neutrinos imprint in cosmic structures, as first noticed in
[122]. One of the most relevant features is the introduction of an additional
scale, known as the free-streaming scale of neutrinos
$k_{\text{fs}}$.191919For realistic neutrino masses, this scale is currently
around $0.1\,h\,\text{Mpc}^{-1}$, at the onset of non-linear scales. This
scale plays a significant role in clustering because, on distance scales
smaller than $k_{\text{fs}}$, the large velocity dispersion of neutrinos
prevents them from clustering together with the CDM-baryon fluid, resulting in
the washout of small-scale perturbations. This effect is manifested as a
suppression in the power spectrum in regions where $k>k_{\text{fs}}$, allowing
to extract information about the total neutrino mass. Neutrinos also impact
other observables, such as the CMB, and hence the more stringent constraints
to their mass will come from combined analyses; see [123, 124] for a review of
neutrinos in cosmology.

To date, the Planck 2018 legacy data have placed rigorous limits on the
neutrino mass through CMB temperature and polarization anisotropies, yielding
$\sum_{i}m_{\nu,i}\lesssim 0.26$ eV (95% c.i.). This limit tightens to $0.12$
eV when combined with additional probes such as CMB-lensing, and BAO, among
others [125].

Here we evaluate the ability of the Full-Modelling analysis to recover the
total neutrino mass in a controlled scenario, where we already know the true
value used in the simulations; see also [80, 39]. In this case, we exclusively
present this test utilizing the Full-Modelling approach, while the results for
the ShapeFit will be detailed in future research. For this test, we allow the
total neutrino mass to vary within the wide flat range
$M_{\nu}\in\mathcal{U}(0,5)$ eV. We utilize one massive and two massless
neutrino states. Other than that, the settings remain the same as for the
baseline analysis, where the total neutrino mass and the ultra-relativistic
species are fixed to $M_{\nu}=0.06\,\text{eV}$ and $N_{\text{ur}}=2.0328$,
respectively [54, 55]. The results are shown in Figure 14 and Table 5.

Our first observation is that the total neutrino mass has a minimal impact on
the measurements of $h$, as also noticed in [21, 126]: While for the CMB data
$M_{\nu}$ and $h$ are highly degenerated, for clustering alone this seems not
to be the case. Additionally, we observe that the total neutrino mass exhibits
an expected degeneracy with the primordial power spectrum amplitude $A_{s}$,
worsening considerably the constraints of this parameter and $\omega_{cdm}$.
However, we see below that this degeneracy can be broken by adding more
tracers. Another interesting finding is that the total neutrino mass exhibits
a weak dependence on the tidal bias parameter $b_{s^{2}}$ and the non-local
third order bias parameter $b_{3\rm nl}$, demonstrating robustness in the
constraints on the total neutrino mass under both the minimal freedom and
maximal freedom cases.

The synthetic data are generated using a different code than Folps,
specifically we use velocileptors.

One of the most important features of Folps is its accurate treatment of
massive neutrinos via the fkPT method [80, 39, 99]. Here, we compare the
neutrino mass constraints by switching between fk\- and EdS-kernels in Folps,
both of which are included in our code [39].202020To do this, the code simply
sets $f(k,z)=f_{0}\equiv f(k=0,z)$. In particular, we estimate how large are
the differences compared to DESI Y1 and DESI Y5 error bars. For this purpose,
we employ noiseless synthetic data generated with an EFT model at the DESI
fiducial cosmology and a set of given nuisance parameters. The synthetic data
presented here are produced with a distinct code from Folps, namely
velocileptors. Nonetheless, based on the findings of [44, 45], we don’t expect
difference when employing any EFT code.

Figure 14: Impact of varying $M_{\nu}$. Comparison of the posterior
distributions between the $\Lambda$CDM baseline analysis (configured as
described in table 2) and an extended analysis where the total neutrino mass
$M_{\nu}$ is relaxed to be a free parameter with a flat prior of
$M_{\nu}\in\mathcal{U}(0,5)$ eV. We present this comparison for the Full-
Modelling analysis, showing that the constraint on $M_{\nu}$ remains largely
unaffected by the tidal bias $b_{s^{2}}$ and the non-local bias $b_{3\rm nl}$,
included as free parameters of the model in the maximal freedom case.

Figure 15: Impact of varying $M_{\nu}$. Comparison of the 1-dimensional
marginalized posterior distribution for the total neutrino mass $M_{\nu}$ (top
panel) and the other cosmological parameters (bottom panel) resulting from the
fitting of the synthetic data for DESI Y1 and DESI Y5 redshift bins with
corresponding analytical covariance matrices.

Figure 15 shows the estimated differences in cosmological parameters between
EdS- and fk-kernels for data mimicking DESI Y1 (solid lines) and DESI Y5 (dot-
dashed lines). In general, all the parameters are recovered within $1\sigma$,
with essentially no differences in $h,\omega_{cdm}$, and $A_{s}$, hence we
only report the constraints on $M_{\nu}$. The $95\%$ limits estimated for the
neutrino mass using the noiseless data are, for DESI Y1,

DESI Y1: $\displaystyle\qquad M_{\nu}<0.368\,\text{eV}\quad\text{(EdS-
kernels)},$ (7.2) DESI Y1: $\displaystyle\qquad
M_{\nu}<0.323\,\text{eV}\quad\text{({fk}-kernels)},$ (7.3)

while for DESI Y5,

DESI Y5: $\displaystyle\qquad M_{\nu}<0.170\,\text{eV}\quad\text{(EdS-
kernels)},$ (7.4) DESI Y5: $\displaystyle\qquad
M_{\nu}<0.149\,\text{eV}\quad\text{({fk}-kernels)}.$ (7.5)

This represents an enhanced constraining power of around $14\%$ when switching
from the usual EdS- to our fk-kernels. A similar effect in the mean and error
bars of the posterior distribution was previously obtained with simulations in
[39], showing the reliability of the fkPT method. The observed differences in
neutrino mass constraints are potentially relevant for extracting neutrino
mass information through full-shape analyses.

###  7.5 Curvature: $k\Lambda$CDM

Figure 16: Curvature $k\Lambda$CDM. Comparison of the posterior distributions
between Full-Modelling and ShapeFit for the extended $k\Lambda$CDM model,
where the curvature parameter $\Omega_{k}$ is a free parameter of the model
with a flat prior $\Omega_{k}\in\mathcal{U}(-0.2,0.2)$. Solid lines represent
fits to the monopole and quadrupole, while dash-dotted lines include the
hexadecapole. All other settings remain consistent with the baseline
$\Lambda$CDM analysis as outlined in Table 2.

Recent studies have drawn attention to a potential tension in early Universe
curvature constraints [127, 128, 129, 130, 131]. So, it’s useful to compare
these results using new approaches like the PT/EFT of LSS. In this section, we
relax the assumption that the Universe is globally flat and let the curvature
parameter $\Omega_{k}$ to vary freely. We use a uniform prior,
$\Omega_{k}\in\mathcal{U}(-0.2,0.2)$, while keeping the other parameter
settings as in the baseline analyses. Using both Full-Modelling and ShapeFit
methods, we obtained constraints showing that the curvature is well centered
at $\Omega_{k}=0$, recovering the simulations’ settings. We present these
results in Figure 16 and Tables 5 and 7, for volumes
$V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$ and $V_{5}=40\,\,h^{-3}\text{Gpc}^{3}$,
respectively. It’s interesting that for the $k\Lambda$CDM cosmological model,
Full-Modelling and ShapeFit give similar results, though ShapeFit contours are
a bit smaller. We find that dropping the flatness assumption introduces
correlations among $\Omega_{k}$, $h$, and $A_{s}$, weakening the constraints
on these parameters but does not affect $w_{cdm}$ much. This is expected
because the curvature does not affect the shape of the power spectrum at the
scales of interest: it only enters through the background evolution.
Additionally, we included the hexadecapole in the analysis, improving the
constraints on $\Omega_{k}$, $h$, and $A_{s}$. For these parameters, the error
bars are reduced by about 12–20% in both the Full-Modelling and ShapeFit
fitting methods.

###  7.6 Dark Energy: $w$CDM

The $\Lambda$CDM model is the most successful cosmological model to date,
supported by numerous observations. However, there are still some fundamental
questions that remain unanswered, such as the nature of dark matter or dark
energy. In order to address these questions, there has been a considerable
improvement in observational instruments in recent times, leading to a larger
number of very accurate observations.

On the theoretical side, there have also been numerous extensions to the
standard $\Lambda$CDM model that aim to address more fundamental questions. In
this section, we will focus on one of these simplest extensions, known as the
$w$CDM model. Here, the dark energy component is described as a barotropic
fluid with an equation of state given by $P=w\rho$, with $P$ and $\rho$ the
pressure and density of the dark energy and $w$ a constant. Hence, under this
parametrization, the energy density of the dark energy dilutes with the scale
factor of the Universe as $\rho\propto a^{-3(1+w)}$. We will assume, following
standard practice, that fluctuations in dark energy density are negligible.
Although $w$ can take arbitrary values, it only affects the background
evolution of the Universe. We allow the equation of state to vary as a free
parameter of the model within the flat range of $w\in\mathcal{U}(-1.8,-0.3)$,
while keeping the same settings as in the $\Lambda$CDM baseline analysis,
which is presented in table 2. We consider the case of the maximum volume
available $V_{25}=200\,\,h^{-3}\text{Gpc}^{3}$.

Figure 17: Dark Energy $w$CDM. Comparison of the posterior distributions
between the baseline $\Lambda$CDM and the extended $w$CDM model, where $w$ is
a free parameter of the model with a flat prior $w\in\mathcal{U}(-1.8,-0.3)$.
All other settings remain consistent with the baseline $\Lambda$CDM analysis
as outlined in Table 2.  Figure 18: Dark Energy $w$CDM. Comparison of the
1-dimensional marginalized posterior distributions between the compressed
ShapeFit and direct Full-Modelling fits within the extended $w$CDM model. This
figure serves as an extension of Figure 17, illustrating the contrast in
constraints when incorporating (dot-dashed lines) or omitting (solid lines)
the hexadecapole in the fits.

In figure 17 and Tables 5 and 7, we present the results for the $w$CDM
cosmology for volumes $V_{25}$ and $V_{5}$, and compare them with the baseline
$\Lambda$CDM model. The results are presented for both compressed ShapeFit and
direct Full-Modelling analyses. Introducing $w$ affects the measurements of
$h$ and $A_{s}$, due to expected degeneracies, but has a small impact on the
constraints for $\omega_{cdm}$. Introducing $w$ shifts the mean to slightly
higher (lower) values for $h$ ($A_{s}$). Additionally, the constraining power
on $h$ and $A_{s}$ worsens by a factor of around 2, compared to the
$\Lambda$CDM case.

Comparing the ShapeFit and Full-Modelling results, we notice that both
approaches provide similar contours and 1-dimensional posterior distributions.
However, ShapeFit tends to have slightly tighter constraints with more
symmetrical distributions.

We further study the impact of adding the hexadecapole power spectrum to the
data vector. In Figure 18 we highlight the contrast in constraining power
between the inclusion and exclusion of $\ell=4$ in the fits. We find that the
constraints on the cosmological parameters are in general slightly better both
in precision and accuracy when including the hexadecapole, which seems
particularly evident for $A_{s}$. The error bars for $w$, $h$, and $A_{s}$ are
reduced by approximately 9–20% in both the Full-Modelling and ShapeFit cases.

|  | FM |  | SF |
---|---|---|---|---|---
|  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$ |  | 68% c.i. | $\,\Delta\Omega/\sigma_{\Omega}$ |
Baseline, $V_{25}$: |  |  |  |  |  |  |
$h$ |  | $0.6741\pm 0.0036$ | $0.15$ |  | $0.6747^{+0.0034}_{-0.0035}$ | $0.31$ |
$\omega_{cdm}$ |  | $0.1197^{+0.0018}_{-0.0015}$ | $0.19$ |  | $0.1197^{+0.0015}_{-0.0014}$ | $0.21$ |
$\ln(10^{10}A_{s})$ |  | $3.032^{+0.017}_{-0.022}$ | $0.20$ |  | $3.038^{+0.020}_{-0.018}$ | $0.08$ |
Hexadecapole, $V_{25}$: |  |  |  |  |  |  |
$h$ |  | $0.6736\pm 0.0037$ | $0.01$ |  | $0.6740\pm 0.0035$ | $0.13$ |
$\omega_{cdm}$ |  | $0.1195^{+0.0018}_{-0.0014}$ | $0.33$ |  | $0.1194\pm 0.0014$ | $0.45$ |
$\ln(10^{10}A_{s})$ |  | $3.040^{+0.014}_{-0.024}$ | $0.17$ |  | $3.048\pm 0.020$ | $0.56$ |
$\omega_{b}$ (flat prior), $V_{25}$: |  |  |  |  |  |  |
$\omega_{b}$ |  | $0.0211^{+0.0012}_{-0.0016}$ | $0.92$ |  | — | — |
$h$ |  | $0.664^{+0.011}_{-0.014}$ | $0.77$ |  | — | — |
$\omega_{cdm}$ |  | $0.1171\pm 0.0033$ | $0.87$ |  | — | — |
$\ln(10^{10}A_{s})$ |  | $3.036^{+0.015}_{-0.027}$ | $0.01$ |  | $3.011^{+0.070}_{-0.12}$ | $0.26$ |
Varying $n_{s}$, $V_{25}$: |  |  |  |  |  |  |
$n_{s}$ |  | $0.954\pm 0.013$ | $0.87$ |  | $0.973\pm 0.024$ | $0.36$ |
$h$ |  | $0.6754\pm 0.0039$ | $0.47$ |  | $0.6740\pm 0.0045$ | $0.09$ |
$\omega_{cdm}$ |  | $0.1216\pm 0.0029$ | $0.56$ |  | $0.1179\pm 0.0063$ | $0.33$ |
$\ln(10^{10}A_{s})$ |  | $3.011\pm 0.034$ | $0.76$ |  | $3.051\pm 0.048$ | $0.31$ |
$\nu\Lambda$CDM, $V_{25}$: |  |  |  |  |  |  |
$M_{\nu}$ |  | $0.106^{+0.049}_{-0.089}$ | $0.67$ |  |  |  |
$h$ |  | $0.6742\pm 0.0035$ | $0.16$ |  |  |  |
$\omega_{cdm}$ |  | $0.1207^{+0.0019}_{-0.0023}$ | $0.35$ |  |  |  |
$\ln(10^{10}A_{s})$ |  | $3.042^{+0.024}_{-0.028}$ | $0.21$ |  |  |  |
$k\Lambda$CDM ($\ell=0,2$), $V_{25}$: |  |  |  |  |  |  |
$\Omega_{k}$ |  | $-0.008\pm 0.015$ | $0.52$ |  | $-0.006^{+0.012}_{-0.015}$ | $0.42$ |
$h$ |  | $0.6765\pm 0.0063$ | $0.46$ |  | $0.6765\pm 0.0053$ | $0.54$ |
$\omega_{cdm}$ |  | $0.1197^{+0.0017}_{-0.0013}$ | $0.20$ |  | $0.1198\pm 0.0014$ | $0.17$ |
$\ln(10^{10}A_{s})$ |  | $3.009\pm 0.049$ | $0.56$ |  | $3.019\pm 0.046$ | $0.37$ |
$k\Lambda$CDM ($\ell=0,2,4$), $V_{25}$: |  |  |  |  |  |  |
$\Omega_{k}$ |  | $-0.006\pm 0.012$ | $0.49$ |  | $-0.003\pm 0.011$ | $0.25$ |
$h$ |  | $0.6755\pm 0.0054$ | $0.36$ |  | $0.6750\pm 0.0047$ | $0.29$ |
$\omega_{cdm}$ |  | $0.1196^{+0.0017}_{-0.0013}$ | $0.28$ |  | $0.1195\pm 0.0014$ | $0.33$ |
$\ln(10^{10}A_{s})$ |  | $3.022^{+0.038}_{-0.043}$ | $0.35$ |  | $3.037\pm 0.040$ | $0.36$ |
$w$CDM ($\ell=0,2$), $V_{25}$: |  |  |  |  |  |  |
$w$ |  | $-1.019\pm 0.034$ | $0.56$ |  | $-1.016\pm 0.033$ | $0.50$ |
$h$ |  | $0.6780\pm 0.0080$ | $0.55$ |  | $0.6783\pm 0.0079$ | $0.60$ |
$\omega_{cdm}$ |  | $0.1197^{+0.0016}_{-0.0013}$ | $0.19$ |  | $0.1198\pm 0.0013$ | $0.17$ |
$\ln(10^{10}A_{s})$ |  | $3.015^{+0.031}_{-0.037}$ | $0.61$ |  | $3.024\pm 0.034$ | $0.35$ |
$w$CDM ($\ell=0,2,4$), $V_{25}$: |  |  |  |  |  |  |
$w$ |  | $-1.014^{+0.029}_{-0.026}$ | $0.49$ |  | $-1.008\pm 0.028$ | $0.30$ |