# Full Modeling and Parameter Compression Methods in configuration space for
DESI 2024 and beyond

S. Ramirez-Solano  0009-0006-2096-7974 M. Icaza-Lizaola  0000-0002-2547-3184
H. E. Noriega  0000-0002-3397-3998 M. Vargas-Magaña  0000-0003-3841-1836 S.
Fromenteau  0000-0002-3724-4768 A. Aviles  0000-0001-5998-3986 F. Rodríguez-
Martínez  J. Aguilar  S. Ahlen  0000-0001-6098-7247 O. Alves  S. Brieden
0000-0003-3896-9215 D. Brooks  T. Claybaugh  S. Cole  0000-0002-5954-7903 A.
de la Macorra  0000-0002-1769-1640 Arjun Dey  0000-0002-4928-4003 B. Dey
0000-0002-5665-7912 P. Doel  K. Fanning  0000-0003-2371-3356 J. E. Forero-
Romero  0000-0002-2890-3725 E. Gaztañaga  H. Gil-Marín  0000-0003-0265-6217 S.
Gontcho A Gontcho  0000-0003-3142-233X K. Honscheid  C. Howlett
0000-0002-1081-9410 S. Juneau  Y. Lai  M. Landriau  0000-0003-1838-8528 M.
Manera  0000-0003-4962-8934 M. Maus  R. Miquel  E. Mueller  A. Muñoz-Gutiérrez
A. D. Myers  S. Nadathur  0000-0001-9070-3102 J. Nie  0000-0001-6590-8122 W.
J. Percival  0000-0002-0644-5727 C. Poppett  M. Rezaie  0000-0001-5589-7116 G.
Rossi  E. Sanchez  0000-0002-9646-8198 D. Schlegel  M. Schubnell  H. Seo
0000-0002-6588-3508 D. Sprayberry  G. Tarlé  0000-0003-1704-0781 L. Verde
0000-0003-2601-8770 B. A. Weaver  R. H. Wechsler  0000-0003-2229-011X S. Yuan
0000-0002-5992-7586 P. Zarrouk  0000-0002-7305-9578 H. Zou
0000-0002-6684-3997

###### Abstract

In the contemporary era of high-precision spectroscopic surveys, led by
projects like DESI, there is an increasing demand for optimizing the
extraction of cosmological information from clustering data. This work
conducts a thorough comparison of various methodologies for modeling the full
shape of the two-point statistics in configuration space. We investigate the
performance of both direct fits (Full-Modeling) and the parameter compression
approaches (ShapeFit and Standard). We utilize the ABACUS-SUMMIT simulations,
tailored to exceed DESI’s precision requirements. Particularly, we fit the
two-point statistics of three distinct tracers (LRG, ELG, and QSO), by
employing a Gaussian Streaming Model in tandem with Convolution Lagrangian
Perturbation Theory and Effective Field Theory. We explore methodological
setup variations, including the range of scales, the set of galaxy bias
parameters, the inclusion of the hexadecapole, as well as model extensions
encompassing varying $n_{s}$ and allowing for $w_{0}w_{a}$CDM dark energy
model. Throughout these varied explorations, while precision levels fluctuate
and certain configurations exhibit tighter parameter constraints, our pipeline
consistently recovers the parameter values of the mocks within $1\sigma$ in
all cases for a 1-year DESI volume. Additionally, we compare the performance
of configuration space analysis with its Fourier space counterpart using three
models: PyBird, FOLPS and velocileptors, presented in companion papers. We
find good agreement with the results from all these models.

##  1 Introduction

The study of the large-scale structure (LSS) of the Universe is a fundamental
pillar in modern cosmology, and it provides an important observable to
investigate its subsequent evolution over time. Observations of the LSS
through the distribution of galaxies enable the extraction of valuable
information that can be used to compare to theoretical predictions and impose
constraints on cosmological parameters.

Spectroscopic galaxy surveys collect data on the angular position and redshift
of millions of galaxies across a wide range of cosmic epochs. Two key
phenomena, which have left distinct features in the clustering of galaxies,
are extensively utilized for extracting cosmological insights.

The first phenomenon involves fluctuations in the matter density, arising from
acoustic density waves in the primordial plasma of the early universe. These
effects are discernible through the peak of the correlation function at
comoving scales of $\sim 100\,h^{-1}\,\textrm{Mpc}$, and are commonly referred
to as baryon acoustic oscillations (BAO).

The second key observable arises from apparent anisotropies in the clustering
of galaxies along the line-of-sight (LoS), attributed to their peculiar
velocities. This phenomenon is known as redshift space distortions (RSD).
While RSD introduce complexities in interpreting the observed positions of
galaxies as a 3-D map, they also provide a unique opportunity to access
information concerning the growth rate of the LSS.

In the last decade, the clustering of galaxies has been analyzed through data
collected by spectroscopic galaxy surveys such as the 2dF Galaxy Redshift
Survey (2dFGRS) [1, 2], the Baryon Oscillation Spectroscopic Survey (BOSS) [3,
4], the extended Baryon Oscillation Spectroscopic Survey (eBOSS) [5, 6], or
the Dark Energy Survey (DES) [7, 8]. These surveys have significantly
increased both precision and volume over the years. At the time of
publication, the next generation of surveys is being developed with the EUCLID
mission [9, 10] launching in the summer of 2023, and The Rubin Observatory
Large Synoptic Survey Telescope (LSST) [11] survey expected to have its first
light in the summer of 2024. Simultaneously, the Dark Energy Spectroscopic
Instrument (DESI) 111www.desi.lbl.gov [12, 13, 14, 15] is currently ongoing
and is mapping the universe over an unprecedented range of
$14,000\,\textrm{deg}^{2}$ with an effective volume of
$40\,(h^{-1}\,\textrm{Gpc})^{3}$, improving the number of observed galaxies by
an order of magnitude compared to previous surveys. Plans for a continuing
survey, DESI-II [16], that will focus on the $z>2$ universe are being made,
and observations are expected to start in 2026.

DESI targets different tracers of matter at various redshifts, including
emission line galaxies (ELG) [17] at $0.6<z<1.6$, luminous red galaxies (LRG)
[18] at $0.3<z<1.0$, and quasars at (QSO) $0.9<z<2.1$ [19] during the dark-
time program. The exceptional accuracy of these surveys will necessitate
models capable of predicting clustering statistics with great precision.

Perturbation theory (PT) when applied to data in either the Lagrangian
approach [20, 21, 22, 23, 24, 25, 26] or Eulerian approach [27, 28, 29, 30]
has demonstrated to accurately predict galaxy densities, velocities and
redshift-space distortions within the quasi-linear regime. Recently, PT has
been complemented with effective field theory (EFT) [31, 32, 33, 34, 35, 36]
parameters to account for back-reaction of small scale physics over quasi-
linear scales. This has proven effective in confronting theoretical
predictions of clustering statistics with observations or simulations. In this
work, we follow the Lagrangian PT (LPT)/EFT approach [37, 38, 35, 39, 40].

In anticipation of the DESI Year 1 data release, the DESI collaboration aims
at testing systematic effects of their methodologies using the new ABACUS-
SUMMIT [41] suite of simulations, which are designed to surpass the resolution
and volume requirements for cosmological simulations within DESI. Currently,
four perturbations approaches for modeling the redshift space clustering
statistics are being analyzed and compared within the DESI collaboration.
These models include the configuration space EFT-GSM [42] code, which is
utilized throughout this work, as well as three Fourier space codes, namely
PyBird [43, 44, 45, 46], FOLPS$\nu$ [47, 48] and velocileptors [49, 50, 51]. A
comparison for the three Fourier space codes can be found in the companion
paper [52] and, in this paper, we compare our configuration space model and
the other Fourier space models.

In this work, we assess potential systematics within the pipeline of our EFT-
GSM model for the DESI first-year data release. EFT-GSM model by construction
is designed to model the redshift space correlation function up to one-loop
corrections working on the LPT framework and using the Gaussian Streaming
Model [53, 54, 55, 56, 35], in addition to EFT parameters. This framework
allows us to model small-scale physics down to scales as low as
$20\,h^{-1}\,\textrm{Mpc}$, where non-linear processes appear. Our methodology
has been introduced and optimized with neural networks in [42] we test it
using the Nseries simulations [4], a suite of N-body simulations that were
designed to analyze the eBOSS data and have a lower resolution than the
ABACUS-SUMMIT simulations. The main difference between our EFT-GSM methodology
and the other three methodologies is that we fit the correlation function
directly, whereas the other methodologies work in Fourier space.222Reference
[46] also works with the correlation function, through Fourier transforming
the power spectrum obtained with PyBird.

Two different approaches have been used for extracting cosmological
information from the shape of two-point statistics, the compressed
methodologies (as the Standard and Shape Fit methodologies), and the direct
fits (also called Full-Modelling). The Standard methodology employed for
analyzing the clustering of redshift space tracers throughout SDSS-III BOSS
[57, 3, 58] and SDSS-IV eBOSS [59, 60, 61] requires a one-time calculation of
the linear power spectrum within a chosen fiducial cosmology. This power
spectrum is then coupled with a PT framework, and the relevant cosmological
information are compressed into three free parameters: the linear growth rate
of structures $f\sigma_{8}$ and two scaling parameters, parallel and
perpendicular to the line of sight $\alpha_{\perp}$ and $\alpha_{\parallel}$,
that parameterize the Alcock-Paczynski effect [62].

The ShapeFit methodology of [63] suggests adding a new free parameter to the
Standard compressed methodology of BOSS and eBOSS, this parameter models the
transition of the linear power spectrum from large to small scales, by
including a new effective parameter $m$, that modifies the slope of the linear
power spectrum. The agnostic nature of compressed methodologies, like the
Standard approach utilized for BOSS and eBOSS, or the more recent ShapeFit
methodology, exhibit a relatively high degree of model independence compared
to those derived from Full Modeling analysis. Moreover, their significantly
lower computational costs make them more appealing from a practical
standpoint. Although the ShapeFit methodology has been tested in Fourier space
[e.g. 63], its performance in configuration space is still an ongoing area of
investigation.

An alternative approach opposed to the compressed nature of the Standard and
Shape Fit methodologies is to directly explore the posterior distribution of
all cosmological parameters. This involves calculating the linear power
spectrum at each step of the Markov Chain Monte Carlo (MCMC) exploration.
Referred to as the Full Modeling or Direct Fit approach, it offers the
advantage of tighter constraints on cosmological parameters, which are
directly constrained. This method has been applied in the analysis of both
BOSS [64, 65] and eBOSS [66, 67]. While these Full Modeling approaches can be
resource-intensive to construct, advancements in computing capabilities and
the evolution of new machine-learning algorithms have accelerated analytical
calculations in recent years. Consequently, it has become feasible to model
the full shape of the correlation function or power spectrum and
simultaneously vary cosmological, bias, and EFT (Effective Field Theory)
parameters.

Throughout this work, we evaluate the performance of the three methodologies
for extracting cosmological information Standard and Shape Fit, and Full
Modeling), and we compare the analysis of these three methods across a range
of conditions and across different DESI tracers to have a better understanding
of the differences among them.

This paper gives support to the DESI 2024: Full-Shape analysis from galaxies
and quasars [68] which is part of the DESI Year 1 main results [69, 70, 71,
72, 73, 74].

The structure of this paper is as follows. We begin in section 2 by
introducing and describing the simulations used throughout this work. Then in
Section 3 we describe the theory behind the EFT-GSM model. Section 4 describes
the implementation of our three different methods: the Standard methodology,
the Shape-Fit methodology and the Full Modeling methodology. The material
presented up until this point is not new to this work. So, a busy reader can
directly skip sections 3 and 4. Section 5 shows the results of our pipeline
for these three distinct methodologies, and explores different configurations
of our model. First, we define our baseline analysis and investigate how
different variations to it affect our results, these variations include
differences in the fitting range explored and different settings for the
galaxy bias terms. We continue by exploring the effects on our models
constraints under different extensions to our baseline configuration, this
extensions include adding the hexadecapole of the correlation function,
extending the parameter space explored by adding $n_{s}$ as a free parameter,
and including the parameters of the equation state of dark energy as free
parameters. We finalize this section by presenting the results for our three
different tracers and presenting a joint analysis where all tracers are
analyzed simultaneously. The last subsection presents the comparison between
Shape-Fit and Full Modeling methodologies for our baseline analysis. Finally,
Section 6 has a comparison between our GSM-EFT model and other EFT models in
Fourier space and configuration space. The conclusions of this paper are in
section 7.

##  2 Simulations mocks and covariance

As we have stated, the main goal of our work is to test the performance and
analyze the potential systematics of our methodology on a mock representation
of the DESI ELG, LRG, and QSO surveys. We begin our work in section 2.1 by
introducing our mock datasets, generated with the ABACUS-SUMMIT simulation.
The various analyses presented in this work involve estimating the likelihood
that the correlation functions produced by our EFT-GSM model agree with those
measured from the ABACUS-SUMMIT simulation across a range of parameters and
across different model configurations. The computation of these likelihoods
requires estimates of the covariance of our correlation function. The
methodology and simulations we utilize to compute this covariance are
introduced in Section 2.2.

###  2.1 Mock Data

In this study, we test the correlation function model using a set of 25 cubic
boxes, with a volume of $\left(2\,h^{-1}\mathrm{Gpc}\right)^{3}$, which were
generated utilizing the ABACUS-SUMMIT high-accuracy $N$-body simulations [41].
Consequently, the combined volume extends to $200\,h^{-3}\mathrm{Gpc}^{3}$.
Each ABACUS-SUMMIT box contains a total of $6912^{3}$ particles, with each
particle carrying a mass of $2\times 10^{9}\,h^{-1}\mathrm{M}_{\odot}$. These
simulations adopt the cosmological parameters from Planck 2018 [75]:
$h=0.6736$, $\omega{\mathrm{cdm}}=0.1200$, $\omega_{\mathrm{b}}=0.02237$,
$\ln(10^{10}A_{s})=3.0364$, $n_{s}=0.9649$, $M_{\nu}=0.06\,\mathrm{eV}$,
$w_{0}=-1$.

Our mocks are based on HOD fitted on the Survey Validation 3 (SV3)
spectroscopic data from DESI and consists of three independent sets of
catalogs, each representing a different tracer: LRG at $z=0.8$, ELG at
$z=1.1$, and the QSO at $z=1.4$. In order to populate our simulations with
galaxies, we utilize the halo occupation distribution (HOD) methodology
described in [76] for ELG and in [77] for LRG and QSO. The HOD models are
calibrated using small-scale wedges (below $5\,h^{-1}\textrm{Mpc}$) in
combination with large-scale bias evolution information, whenever the former
is available. These calibrations are done independently for each snapshot of
the simulation.

Through this work we utilize the correlation function that we measure from
these Abacus mock simulations, these correlation functions are utilized to fit
the parameters of our models. Our estimates of the correlation function are
defined as the mean of these 25 mock realizations for each tracer, where the
$s$-bins in our analysis are spaced at intervals of $\Delta s=4$
$h^{-1}\mathrm{Mpc}$. In Figure 1 we show the mean of the ABACUS multipoles
for $\ell=0,2,4$ for the three tracers, LRG (left panel), ELG(intermediate
panel) and QSO(right panel).

###  2.2 Covariance Mocks

One common approach for computing covariance matrices is to utilize the
correlation function measured from thousands of different realizations. These
simulations are usually computed by methods that are optimized for efficiency
rather than accuracy. In this work, we utilize 1000 approximate mocks
generated with the Extended Zel’dovich methodology [78] with the same
cosmology of ABACUS-SUMMIT. We compute the covariance matrix $C$ of each
tracer between bins $i$ and $j$, using the 1000 EZ mocks, as follows:

$\displaystyle
C_{s}^{(ij)}=\frac{1}{N_{\mathrm{mocks}}-1}\sum_{m=1}^{N_{\mathrm{mocks}}}\left(\xi_{i}^{m}-\bar{\xi}_{i}\right)\left(\xi_{j}^{m}-\bar{\xi}_{j}\right),$
(2.1)

where $N_{\text{mocks}}=1000$ represents the total number of mocks, and
$\bar{\xi}_{i}$ denotes the average of the $i^{th}$ bin in the analysis. The
Hartlap factor induces a slight change in the inverse covariance matrix,
typically around 0.94 to 0.95, depending on the scales. Hence, its influence
in this study is thought to be negligible, as the values vary by less than
0.5%.

In this work, we use three re-scaling for the covariance: without re-scaling
that corresponds to a volume of one simulation $8\,h^{-3}\,\textrm{Gpc}^{3}$,
re-scaling by a factor 1/5 that corresponds to a volume of
$V_{5}=40\,h^{-3}\,\textrm{Gpc}^{3}$ and re-scaling 1/25 corresponding to
$V_{25}=200\,h^{-3}\,\textrm{Gpc}^{3}$. In Figure 1 we show as shaded region
the square root of the diagonal terms of the EZ mocks covariance without re-
scaling for the three tracers, LRG (left panel), ELG (intermediate panel), and
QSO (right panel). We also add error bars corresponding to the 1/25 re-scaled
covariance.

Figure 1: ABACUS-SUMMIT simulations for the multipoles of the correlation
function for three tracers: LRGs (left), ELGs (middle), and QSOs (right). The
shadow region represents the errors of the covariance matrix, while the error
bar represents the same error re-scaled by $1/\sqrt{25}$.

##  3 Modeling Redshift Space Distortions with GSM-EFT

In this section, we introduce the methodology we follow to compute the
theoretical templates of the correlation function using our EFT-GSM model. We
begin in section 3.1 by briefly introducing the theoretical background behind
GSM, which in addition to the convolution Lagrangian perturbation theory
(CLPT) and EFT, allows us to model the redshift space two-point correlation
function. Then, in section 3.2, we summarize the methodology followed to build
the machine learning emulators that enhance the efficiency and speed of our
analysis.

###  3.1 Gaussian streaming model

The GSM-EFT code333https://github.com/alejandroaviles/gsm [79] is based on the
LPT approach, where each cold dark matter fluid element position is
parameterized by its initial location $\bm{q}$ and the Lagrangian displacement
field $\mathbf{\Psi}(\bm{q},t)$, as
$\bm{x}(\bm{q},t)=\bm{q}+\mathbf{\Psi}(\bm{q},t)$. Here, the displacement is
governed by the Euler equation
$\partial_{t}^{2}{\bf\Psi}+2H\partial_{t}{\bf\Psi}=-\frac{1}{a^{2}}\nabla{\bf\Phi}({\bf
q}+{\bf\Psi})$ and, which we solve perturbatively by expanding the
displacement as,

$\bf{\Psi}(\bf{q})=\bf{\Psi}^{(1)}(\bf{q})+\bf{\Psi}^{(2)}(\bf{q})+\cdots.$
(3.1)

Since we work with the observed position $\bm{s}$, we must consider the LoS
anisotropies introduced by particles’ peculiar velocities. By doing this, the
map between Lagrangian and redshift space coordinates is given by,

$\bm{s}=\bm{q}+\mathbf{\Psi}+\hat{\bm{n}}\frac{\mathbf{\dot{\Psi}}\cdot\hat{\bm{n}}}{H},$
(3.2)

where $\hat{\bm{n}}$ is a unitary vector on the LoS direction. The correlation
function in redshift space is obtained through the conservation of the number
of tracers $X$, expressed as
$\big{[}1+\delta_{s}(\bm{s})\big{]}d^{3}s=\big{[}1+\delta_{X}(\bm{x})\big{]}d^{3}x$.
Following [54, 80], one obtains

$1+\xi_{s}(\bm{s})=\int\frac{d^{3}k}{(2\pi)^{3}}d^{3}x\,e^{i\bm{k}\cdot(\bm{s}-\bm{r})}\exp\left[\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}k_{i_{1}}\cdots
k_{i_{n}}\mathcal{C}^{(n)}_{i_{1}\cdots i_{n}}(\bm{r})\right],$ (3.3)

where the relation between cumulants $\mathcal{C}$ and moments $\Xi$, of the
pairwise density weighted velocity distribution function, up to second order
is given by

$\displaystyle\mathcal{C}^{(0)}(\bm{r})$ $\displaystyle=\log[1+\xi(r)],$ (3.4)
$\displaystyle\mathcal{C}^{(1)}_{i}(\bm{r})$
$\displaystyle=\frac{\Xi^{(1)}_{i}(\bm{r})}{1+\xi(r)}\equiv v_{12,i},$ (3.5)
$\displaystyle\mathcal{C}^{(2)}_{ij}(\bm{r})$
$\displaystyle=\frac{\Xi^{(2)}_{ij}(\bm{r})}{1+\xi(r)}-\mathcal{C}^{(1)}_{i}(\bm{r})\mathcal{C}^{(1)}_{j}(\bm{r})\equiv\hat{\sigma}^{2}_{12,ij}-v_{12,i}v_{12,j}=\sigma^{2}_{12,ij},$
(3.6)

where $\xi(r)$ is the tracers real space correlation function and the the
moments are given by

$\Xi^{(n)}_{i_{1}\cdots
i_{n}}(\bm{r})=\langle(1+\delta(\bm{r}))(1+\delta(0))\Delta
u_{i_{1}}\cdots\Delta u_{i_{n}}\rangle,$ (3.7)

where $\Delta u_{i}=u_{i}(\bm{r})-u_{i}(0)$, and
$u_{i}=\hat{n}_{i}\hat{n}_{j}\dot{\Psi}_{j}/H$ is the peculiar velocity field
along the LoS in Hubble function units.

We have further introduced the pairwise velocity along the LoS $v_{12,i}$, and
the pairwise velocity dispersion moment and cumulant,
$\hat{\sigma}^{2}_{12,ij}$ and $\sigma^{2}_{12,ij}$, respectively. Throughout
this work we employ CLPT to compute the moments [37, 81, 38, 82, 39, 83]. A
detailed description of our implementation of CLPT into the GSM can be found
in [42]. The precise equations we use are given below, in §4.2.1.

By truncating the sum in the exponential in eq. (3.3) up to the second moment,
we obtain the standard expression for the Gaussian Streaming Model correlation
function:

$1+\xi_{s}(s_{\parallel},s_{\perp})=\int_{-\infty}^{\infty}\frac{dr_{\parallel}}{\big{[}2\pi\sigma^{2}_{12}(r,\mu)\big{]}^{1/2}}\big{[}1+\xi(r)\big{]}\exp\left[-\frac{\left(s_{\parallel}-r_{\parallel}-\mu
v_{12}(r)\right)^{2}}{2\sigma^{2}_{12}(r,\mu)}\right],$ (3.8)

where the pairwise velocity and pairwise velocity dispersion are projected
along the separation vector: $v_{12}=v_{12,i}\hat{r}_{i}$ and
$\sigma_{12}=\sigma_{12,i}\hat{r}_{i}\hat{r}_{j}$. The labels ∥ and ⟂ are used
to indicate components parallel and perpendicular to the LoS direction and
$\mu$ is the cosine of the angle between the separation direction of the two
galaxies and the LoS.

As noticed in [39], the second moment of the density distribution is very
sensitive to small scale physics, hence the velocity dispersion should be
corrected by an EFT parameter, such that

$\sigma_{12}^{2}(r)\rightarrow\sigma_{12}^{2}(r)+\sigma^{2}_{\text{EFT}}\frac{1+\xi^{\text{ZA}}(r)}{1+\xi(r)},$
(3.9)

where $\xi^{\text{ZA}}(r)$ is the Zel’dovich Approximation matter correlation
function. One may notice that for sufficiently large $r$, one has
$1+\xi(r)\approx 1$, recovering the phenomenological description of the
dispersion due to the Fingers of God effect:
$\sigma^{2}_{\text{EFT}}(1+\xi^{\text{ZA}}(r))/(1+\xi(r))\rightarrow\sigma^{2}_{\text{FoG}}$.
In our final modeling we will use a second counterterm as shown in section
4.2.1.

To compute the real space correlation function $\xi$, the pairwise velocity
$v_{12}$ and velocity dispersion $\sigma_{12}$ we utilize CLPT/EFT [38, 35] as
implemented in [39]. Specifically, we use eqs. (4.2.1), (4.2.1) and (4.2.1)
displayed below. Those equations also specify how the biases and counter-terms
(apart from $\sigma^{2}_{\text{EFT}}$ introduced above) enter our correlation
function modeling.

###  3.2 Neural Network Accelerators

The methodology outlined in the previous section has been implemented into our
GSM-EFT code, enabling the prediction of the galaxy correlation function for a
given set of cosmological parameters in just a couple of seconds. Despite this
relative speed, the sheer volume of evaluations required for a MCMC analysis
typically falls on the order of $10^{5}$. Combined with the numerous
independent MCMC analyses presented throughout this work, this incentivizes us
to accelerate the evaluation time of our models using neural network
emulators.

These emulators work as surrogate models that learn to predict the correlation
function of the parameters in approximately 0.015 seconds, resulting in a
speedup of two orders of magnitude once trained. Throughout this work, we
construct a distinct neural network for each redshift and for each multipole
of the correlation function. Each of these emulators requires around 30
minutes to train. Our neural network approach is an adaptation of the code
presented in [84]444Their code is publicly available and can be found at
https://github.com/sfschen/EmulateLSS/tree/main, which we modified to work in
configuration space instead of Fourier space, as originally designed. The
neural network architecture utilizes a Multi-Layer Perceptron with four hidden
layers, each containing 128 neurons. In this work, we construct one individual
neural network per multipole.

Figure 2 shows a scheme of our methodology. We implemented a decreasing
learning rate strategy, starting from $10^{-2}$ and decreasing to $10^{-6}$ in
steps of one order of magnitude. We reduce the learning rate with a patience
of 1000 epochs, i.e., if there is no improvement for 1000 epochs. At each
step, the batch size is doubled. Additionally, we employ the following
activation function, as suggested by [85, 84]:

$a(X)=\left[\gamma+(1+e^{-\beta\odot X})^{-1}(1-\gamma)\right]\odot X.$ (3.10)

These emulators were first introduced and thoroughly tested in [42]. Here, we
provide a brief description of the emulator methodology, with a reference to
our original work for a more detailed discussion.

To train a neural network, we initially generate 60,000 points in parameter
space. This holds true for most results presented in this work, with the
exception of certain extensions to our baseline analysis where the models are
more complex, which requires a larger training set, as discussed in sections
§5.2.2 and §5.2.3 below. Of these points, 50,000 are utilized as the training
data, fed into the training algorithm. Additionally, we allocate 5,000 points
as the validation set, which the neural network employs to assess the accuracy
of the model at each training epoch. Based on the performance of the model,
the algorithm may adjust the learning rate or finalize the model. The
remaining 5,000 points constitute the test set, which we employ to evaluate
the model’s accuracy on unseen data. All three sets are independently
generated using a Korobov design [86] to construct uniform samples in a large
dimensional space.

$\omega_{\textrm{cdm}}$$\omega_{b}$⋮$b_{s^{2}}$$y_{0}^{(1)}$$y_{1}^{(1)}$⋮$y_{128}^{(1)}$$\ldots$$\ldots$$\ldots$$y_{0}^{(4)}$$y_{1}^{(4)}$⋮$y_{128}^{(4)}$$\xi_{\ell}(s_{1})$$\xi_{\ell}(s_{2})$⋮$\xi_{\ell}(s_{N})$Model
parameters$1^{\text{st}}$ hidden layer$4^{\text{th}}$ hidden layerMultipole
values Figure 2: Architecture representation of our three Neural Networks, one
identical for each multipole ($\ell=0,2,4$) of the correlation function.
Figure 3: Schemes for the Full Modeling (left panel) and the ShapeFit (right
panel) analysis, allowing to show the differences in the modeling and the
cosmological parameters extraction of the information.

##  4 Extracting Cosmological Information from Clustering

As we illustrate in Figure 3, we can distinguish two different approaches for
extracting cosmological information from clustering:

1) Compressed approaches [e.g. 57, 3, 58, 87, 59, 60, 61] shown in the right
panel, where the philosophy is to extract the cosmological information from
the measured clustering data $\\{\xi(s,\mu),p(k,\mu)\\}$ through the
constraints on compressed variables, $\Theta_{\texttt{obs}}$, that are related
with physical observables. These compressed variables are further translated
into constraints to cosmological parameters using different cosmological
models. In this approach, the constraints on the physical observables are
proven to be mildly affected by the cosmological model and, for this reason,
it is considered agnostic to the cosmological model. 2) Full Modeling or
Direct fits [e.g. 88, 89, 90, 91, 92, 93] shown in the left panel, where the
extraction of cosmological information from the measured clustering data
$\\{\xi(s,\mu),p(k,\mu)\\}$ is done by comparing them directly with
theoretical predictions computed for each cosmological model constraining
directly the cosmological (and nuisance) parameters $\theta_{\Omega}$ of the
model, without any intermediate step. Full Modeling has received increasing
attention post-Stage-III survey’s period, given the optimization of models and
codes and the use of sophisticated statistical tools in the last decade, that
allows us to compute efficiently the posterior in a highly dimensional space
while varying directly the cosmological and nuisance parameters of the model.

Within the compressed approaches, we can distinguish two methodologies: 1)
Classical or Standard methodology proposed and developed within spectroscopic
survey collaborations (as SDSS-II, BOSS and eBOSS), which is based on
compressing the cosmological information into measuring physical observables
denoted as BAO and RSD; and 2) ShapeFit [63] approach that is an extension of
the classical approach aiming to achieve more constraining power through the
inclusion of two additional parameters that bring information about the early
universe. In this section, we briefly describe the three different approaches
and their implementation in configuration space for the model used in this
work, the GSM-EFT.

###  4.1 Classical methodology

The Classical or Standard methodology

is based on a fixed template, which is generated using a particular
perturbative model, computed for a fixed reference cosmology. The fixed
template is modified by three compressed variables associated with late-time
dynamics:

$\\{\alpha_{\parallel},\alpha_{\perp},f\sigma_{8}\\},$ (4.1)

$\alpha_{\perp}$ and $\alpha_{\parallel}$ parameterize distortions parallel
and perpendicular to the line of sight and control the variations along the
$r$-axis. The $f\sigma_{8}$ parameter varies the relative amplitude of our
multipoles. The quantity $f\sigma_{8}$ is the product of the linear growth
rate of structures $f$ defined by,

$f=\frac{d\ln D(a)}{d\ln a},$ (4.2)

where $D$ is the growing mode solution. The amplitude of the matter power
spectrum fluctuations on scales up to $8\,h^{-1}\,\text{Mpc}$, $\sigma_{8}$,
defined as

$\sigma_{8}^{2}(z)=\int_{0}^{\infty}\mathrm{d}(\mathrm{ln}\,k)\,k^{3}P(k,z)W_{TH}^{2}(k\,8h^{-1}\mathrm{Mpc})$
(4.3)

where $W_{TH}$ is a spherical top hat filter. In the standard approach,
$\sigma_{8}$ is usually fixed by the reference cosmology and is completely
degenerate with $f$.

The dilation parameters, $\alpha$’s, are defined as follows:

$\alpha_{\perp}\left(z\right)=\frac{D_{A}\left(z\right)r_{\mathrm{d}}^{\mathrm{ref}}}{D_{A}^{\mathrm{ref}}\left(z\right)r_{\mathrm{d}}},\quad\alpha_{\|}\left(z\right)=\frac{H^{\mathrm{ref}}\left(z\right)r_{\mathrm{d}}^{\mathrm{ref}}}{H\left(z\right)r_{\mathrm{d}}},$
(4.4)

where $D_{A}$ denotes the angular diameter distance, $H$ represents the Hubble
parameter, and $r_{\mathrm{d}}$ stands for the sound horizon at the drag
epoch. The superscript ref means that the estimation is conducted within the
reference cosmology555Here we refer with reference cosmology as a synonym of
fiducial cosmology, we are not distinguishing between grid cosmology or
template cosmology.. In our model, we adopt the alternative parametrization of
the distortion parameters:

$\alpha=\alpha_{||}^{1/3}\alpha_{\perp}^{2/3},\quad\epsilon=\left(\frac{\alpha_{||}}{\alpha_{\perp}}\right)^{1/3}.$
(4.5)

We incorporate distortions in the clustering by substituting $s$ with
$s\prime(s_{\mathrm{ref}},\mu_{\mathrm{ref}})$ and $\mu$ with
$\mu\prime(\mu_{\mathrm{ref}})$. These substitutions are computed using the
re-scaling parameters $\alpha$ and $\epsilon$ in the following way:

$\displaystyle
s\prime(s_{\mathrm{ref}},\mu_{\mathrm{ref}})=s_{\mathrm{ref}}\;\alpha\sqrt{(1+\epsilon)^{4}\mu_{\mathrm{ref}}^{2}+(1-\mu_{\mathrm{ref}}^{2})(1+\epsilon)^{-2}},$
(4.6)
$\displaystyle\quad\mu\prime^{2}(\mu_{\mathrm{ref}})=\left[1+\left(\frac{1}{\mu_{\mathrm{ref}}^{2}}-1\right)(1+\epsilon)^{-6}\right]^{-1}.$
(4.7)

The multipoles, denoted as $\xi_{\ell}(s_{\mathrm{ref}})$, are estimated
within the reference cosmology using the transformed coordinates
$s\prime(s_{\mathrm{ref}},\mu_{\mathrm{ref}})$ and
$\mu\prime(\mu_{\mathrm{ref}})$.

To incorporate the dilation parameters into our implementation, we employ
interpolation for each modeled multipole, $\xi_{\ell}^{\mathrm{model}}$,
utilizing the transformed coordinates
$s\prime(s_{\mathrm{ref}},\mu_{\mathrm{ref}})$. Additionally, we calculate the
observed Legendre polynomials, $\mathcal{L}^{\mathrm{obs}}(\mu\prime)$,
employing the transformed coordinate $\mu\prime(\mu_{\mathrm{ref}})$. We
formulate
$\xi^{\mathrm{obs}}(s\prime(s_{\mathrm{ref}},\mu_{\mathrm{ref}}),\mu\prime(\mu_{\mathrm{ref}}))$
by combining the product of each multipole and its corresponding Legendre
polynomial. The observed multipoles, are expressed as
$\xi_{\ell}(s\prime)=\sum_{\ell\prime}a_{\ell\ell\prime}\xi_{\ell\prime}(s)$,
where the coefficients $a_{\ell\ell\prime}$ are involved in this
summation666As shown in previous work [42], an excellent approximation is
achieved by truncating the summation at $\ell\prime=4$..

In our Standard methodology, we fit the data multipoles to the observed
multipoles, and then, we get constraints on the compressed parameters
$\Theta_{\texttt{obs}}$ and their covariance $C_{\Theta}$. The conversion of
compressed parameter constraints into cosmological constraints is done
separately and is described in §4.3 below.

###  4.2 ShapeFit methodology

The ShapeFit methodology [94, 95, 96, 97] emerges as an extension to the
Standard methodology, aiming to preserve its agnostic nature while
incorporating additional information to enhance its constraining power, thus
making it competitive with Full Modeling.

This is done by incorporating two additional parameters into the Standard
methodology variables, compressing the data vector into five parameters
${\alpha_{||},\alpha_{\perp},f\sigma_{s8},m,n}$. The two new parameters affect
the slope of the power spectrum, in two different ways: first, an overall re-
scaling of the power spectrum slope, through the parameter $n$, which in
principle accounts for a wrong fiducial spectral index $n_{s}$; and second, a
scale-dependent slope, through the parameter $m$, that models the shape of the
power spectrum mainly at wave-numbers around and above the turnaround. This
re-scaling process is applied to the template power spectrum
$P_{\mathrm{ref}}$ to produce a trial spectrum $P_{\mathrm{ref}}^{\prime}$,
defined as follows

$\ln\left(\frac{P_{\mathrm{ref}}^{\prime}(k)}{P_{\mathrm{ref}}(k)}\right)=\frac{m}{a}\tanh\left[a\ln\left(\frac{k}{k_{p}}\right)\right]+n\ln\left(\frac{k}{k_{p}}\right),$
(4.8)

where the pivot scale $k_{p}=\pi/r_{d}$ corresponds to the location where the
baryon suppression reaches its maximum.

The two slope re-scaling parameters recover information from the early
universe that is imprinted in the broadband shape of the matter power
spectrum, this information is omitted in the Standard methodology. In order to
extract cosmological information, the slope $m$ is interpreted in terms of the
cosmological parameters through [63]

$\displaystyle m$ $\displaystyle=\left.\frac{d}{d\ln
k}\left(\ln\left[\frac{T^{2}_{\text{no-wiggle
}}\left(k_{p}/s,\Theta\right)}{T^{2}_{\text{no-wiggle
}}\left(k_{p},\Theta^{\mathrm{ref}}\right)}\right]\right)\right|_{k=k_{p}}$
$\displaystyle=\left.\frac{d}{d\ln k}\left(\ln\left[\frac{P_{\text{no-wiggle
}}^{\text{lin }}\left(k_{p}/s,\bm{\Theta}\right)}{P_{\text{no-wiggle
}}^{\text{lin
}}\left(k_{p},\bm{\Theta^{\mathrm{ref}}}\right)}\frac{\mathcal{P}_{\mathcal{R}}(k,\Theta^{\mathrm{ref}})}{\mathcal{P}_{\mathcal{R}}(k/s,\Theta)}\right]\right)\right|_{k=k_{p}}$
(4.9)

where $\mathcal{P}_{\mathcal{R}}$ is the primordial power spectrum and $T$ the
transfer function. To compute the de-wiggled linear power spectrum,
$P_{\text{no-wiggle }}^{\operatorname{lin}}$ we use the Eisenstein-Hu formulae
of [98].

In most situations, the parameter $n$ is fixed to zero. However, by letting it
free, one can obtain estimations for $n_{s}$ through:

$n=n_{s}-n^{\text{ref}}_{s},$ (4.10)

where $n^{\text{ref}}_{s}$ is the spectral index of the reference cosmology.

Another modification to the Standard methodology introduced by ShapeFit
involves redefining the amplitude of the clustering. Instead of using the
usual $\sigma_{8}$ parameterization, we introduce $\sigma_{s8}$ as a re-
scaling in units of the sound horizon:

$\sigma_{s8}^{2}(z,\Omega)=\int_{0}^{\inf}\mathrm{d}(\mathrm{ln}\,k)\,(ks)^{3}P(ks,z)W_{\text{TH}}^{2}(ks\,8h^{-1}\mathrm{Mpc})$
(4.11)

where $s=r_{d}/r_{d}^{\mathrm{ref}}$. By construction, this re-scaling does
not vary under the fitting process, and all the compressed variables are
expressed in units of the sound horizon.

The templates used in ShapeFit include non-linearities of the power spectrum,
which are commonly computed using Eulerian Perturbation Theory. However, to
speed up computations, one can use approximations as,

$\displaystyle P_{22}(k)$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}}P_{\mathrm{ref}}^{\prime}(p)P_{\mathrm{ref}}^{\prime}(|\bm{k}-\bm{p}|)F_{2}(\bm{k},\bm{k}-\bm{p})$
(4.12)
$\displaystyle\approx\left(\frac{P_{\mathrm{ref}}^{\prime}(k)}{P_{\mathrm{ref}}(k)}\right)^{2}\int\frac{d^{3}p}{(2\pi)^{3}}P_{\mathrm{ref}}(p)P_{\mathrm{ref}}(|\bm{k}-\bm{p}|)F_{2}(\bm{k},\bm{k}-\bm{p}),$
(4.13)

Here, $P_{\mathrm{ref}}(p)$ denotes the linear power spectrum of the reference
template, while $P_{\mathrm{ref}}^{\prime}(p)$ represents the linear power
spectrum modified by the $m$ ShapeFit parameter. Similarly for the rest of the
loop integrals. A validity test of such approximations can be found in
[48]777This approximation is important for ShapeFit approach because this
computation need to be done at each step of the MCMC chains..

####  4.2.1 ShapeFit in Configuration Space

As we illustrated in Figure 3, for the computation of the correlation
function, we use the GSM presented in §3.1 and the Convolution LPT framework
[38] to obtain the necessary ingredients entering eq. (3.8), namely the real
space correlation function $\xi(r)$, the pairwise velocity $v_{12}(r)$ and the
pairwise velocity dispersion $\sigma_{12}(r)$.

For example, let us consider the simplified version of the one-loop dark
matter (unbiased) real space correlation function [38] that is given by

$\displaystyle
1+\xi(r)=\int\frac{d^{3}q}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}\Bigg{\\{}1-\frac{1}{2}A_{ij}^{\text{loop}}G_{ij}+\frac{1}{6}\Gamma_{ijk}W_{ijk}\Bigg{\\}},$
(4.14)

with $q$-functions
$A_{ij}(\bm{q})=\langle\Delta_{i}(\bm{q})\Delta_{j}(\bm{q})\rangle_{c}$, and
$W_{ijk}=\langle\Delta_{i}(\bm{q})\Delta_{j}(\bm{q})\Delta_{k}(\bm{q})\rangle_{c}$,
where $\Delta_{i}=\Psi_{i}(\bm{q})-\Psi_{i}(0)$. In the case of matrix
$\mathbf{A}$, which is split into linear and loop pieces

$\displaystyle A_{ij}(\bm{q})$
$\displaystyle=2\int\frac{d^{3}k}{(2\pi)^{3}}\big{(}1-e^{i\bm{p}\cdot\bm{q}}\big{)}\frac{k_{i}k_{j}}{k^{4}}\left(P_{L}(k)+\frac{3}{7}Q_{1}(k)+\frac{10}{21}R_{1}(k)\right)$
$\displaystyle=A^{L}_{ij}(\bm{q})+A^{\text{loop}}_{ij}(\bm{q}),$ (4.15)

where the $k$-functions $Q_{1}(k)$ and $R_{1}(k)$ are functions constructed
out of four-point functions of linear density fields, defined in [37] as,

$\displaystyle Q_{1}(k)$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}}q_{1}(\bm{k},\bm{p})P_{L}(p)P_{L}(|\bm{k}-\bm{p}|),$
(4.16) $\displaystyle R_{1}(k)$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}}r_{1}(\bm{k},\bm{p})P_{L}(k)P_{L}(p),$
(4.17)

for kernels $q_{1}$ and $r_{1}$. In ShapeFit we produce the trial matrix
$A\prime_{ij}$, by substituting the trial power spectrum $P\prime_{ref}$ of
eq. (4.8) into eq. (4.2.1). To speed up the computations, as it is done in the
case of the non-linear power spectrum, we evaluate functions $Q_{1}$ and
$R_{1}$ as

$\displaystyle Q_{1}(k;m,n)$
$\displaystyle=\left(\frac{P_{\mathrm{ref}}^{\prime}(k;m,n)}{P_{\mathrm{ref}}(k)}\right)^{2}\int\frac{d^{3}p}{(2\pi)^{3}}q_{1}(\bm{k},\bm{p})P_{\mathrm{ref}}(p)P_{\mathrm{ref}}(|\bm{k}-\bm{p}|),$
(4.18) $\displaystyle R_{1}(k;m,n)$
$\displaystyle=\left(\frac{P_{\mathrm{ref}}^{\prime}(k;m,n)}{P_{\mathrm{ref}}(k)}\right)^{2}\int\frac{d^{3}p}{(2\pi)^{3}}r_{1}(\bm{k},\bm{p})P_{\mathrm{ref}}(k)P_{\mathrm{ref}}(p),$
(4.19)

where we have written explicitly the dependence of the trial functions on the
ShapeFit parameters.

This method is analogously applied for all $k$-functions, from which the
$q$-functions are constructed, and finally for the real space correlation
function $\xi(r)$, the pairwise velocity $v_{12}(r)$, and velocity dispersion
$\sigma_{12}(r)$, which are necessary to construct the redshift space GSM
correlation function.

These ingredients are computed for biased tracers using CLPT, obtaining, The
real space correlation function $\xi_{X}(r)$, which corresponds to the zeroth-
order moment for tracer $X$, is obtained within CLPT, [37, 38, 35, 82, 80,
40],

$\displaystyle
1+\xi_{X}(r)=\int\frac{d^{3}q}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}\Bigg{\\{}1-\frac{1}{2}A_{ij}^{loop}G_{ij}+\frac{1}{6}\Gamma_{ijk}W_{ijk}$
$\displaystyle\quad+b_{1}(-2U_{i}g_{i}-A^{10}_{ij}G_{ij})+b_{1}^{2}(\xi_{L}-U_{i}U_{j}G_{ij}-U_{i}^{11}g_{i})+b_{2}(\frac{1}{2}\xi_{L}^{2}-U_{i}^{20}g_{i}-U_{i}U_{j}G_{ij})$
$\displaystyle\quad-2b_{1}b_{2}\xi_{L}U_{i}g_{i}-2b_{s^{2}}g_{i}V_{i}^{10}+b_{s^{2}}^{2}\zeta-2b_{1}b_{s^{2}}g_{i}V_{i}^{12}+b_{2}b_{s^{2}}\chi^{12}-\frac{1}{2}c_{1,\text{EFT}}\operatorname{tr}G\Bigg{\\}},$
(4.20)

for bias parameters $b_{1}$, $b_{2}$ and $b_{s^{2}}$ and EFT counterterm
$c_{1,\text{EFT}}$. The linear correlation function is

$\xi_{L}(q)=\int\frac{d^{3}p}{(2\pi)^{3}}e^{i\bm{p}\cdot\bm{q}}P_{L}(p),$
(4.21)

and we used the functions

$\displaystyle
A_{ij}^{mn}=\langle\delta^{m}(\bm{q})\delta^{n}(0)\Delta_{i}\Delta_{i}\rangle_{c},\qquad
U^{mn}_{i}=\langle\delta^{m}(\bm{q})\delta^{n}(0)\Delta_{i}\rangle_{c}.$
(4.22)

The involved $r$ and $q$ dependent tensors are
$g_{i}=(\mathbf{A}_{L}^{-1})_{ij}(r_{j}-q_{j})$,
$G_{ij}=(\mathbf{A}_{L}^{-1})_{ij}-g_{i}g_{j}$, and
$\Gamma_{ijk}=(\mathbf{A}_{L}^{-1})_{\\{ij}g_{k\\}}-g_{i}g_{j}g_{k}$.

The pairwise velocity is given by

$\displaystyle
v_{12}(\bm{r})=\frac{f\,\hat{r}_{i}}{1+\xi_{X}(r)}\int\frac{d^{3}q\,e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}\Bigg{\\{}-g_{r}\dot{A}_{ri}-\frac{1}{2}G_{rs}\dot{W}_{rsi}$
$\displaystyle\quad+b_{1}\left(2\dot{U}_{i}-2g_{r}\dot{A}^{10}_{ri}-2G_{rs}U_{r}\dot{A}_{si}\right)+b_{1}^{2}\left(\dot{U}^{11}_{i}-2g_{r}U_{r}\dot{U}_{i}-g_{r}\dot{A}_{ri}\xi_{L}\right)$
$\displaystyle\quad+b_{2}\left(\dot{U}^{20}_{i}-2g_{r}U_{r}\dot{U}_{i}\right)+2b_{1}b_{2}\xi_{L}\dot{U}_{i}\Bigg{\\}},$
(4.23)

and the pairwise velocity dispersion

$\displaystyle\sigma^{2}_{12}(\bm{r})=\frac{f^{2}\,\hat{r}_{i}\hat{r}_{j}}{1+\xi_{X}(r)}\int\frac{d^{3}q\,e^{-\frac{1}{2}(\bm{r}-\bm{q})^{\mathbf{T}}\mathbf{A}_{L}^{-1}(\bm{r}-\bm{q})}}{(2\pi)^{3/2}|\mathbf{A}_{L}|^{1/2}}\Bigg{\\{}\ddot{A}_{ij}-g_{r}\ddot{W}_{rij}-G_{rs}\dot{A}_{ri}\dot{A}_{sj}$
$\displaystyle\quad+2b_{1}\left(\ddot{A}^{10}_{ij}-g_{r}\dot{A}_{r\\{i}\dot{U}_{j\\}}-g_{r}U_{r}\ddot{A}_{ij}\right)+b_{1}^{2}\left(\xi_{L}\ddot{A}_{ij}+2\dot{U}_{i}\dot{U}_{j}\right)+2b_{2}\dot{U}_{i}\dot{U}_{j}\,\Bigg{\\}},$
(4.24)

with $q$-coordinate dependent correlators

$\displaystyle\dot{A}_{ij}^{mn}(\bm{q})=\frac{1}{fH}\langle\delta^{m}_{1}\delta^{n}_{2}\Delta_{i}\dot{\Delta}_{j}\rangle,\qquad\ddot{A}_{ij}^{mn}(\bm{q})=\frac{1}{f^{2}H^{2}}\langle\delta^{m}_{1}\delta^{n}_{2}\dot{\Delta}_{i}\dot{\Delta}_{j}\rangle,$
$\displaystyle\dot{W}_{ijk}=\frac{1}{fH}\langle\Delta_{i}\Delta_{j}\dot{\Delta}_{k}\rangle,\qquad\ddot{W}_{ijk}=\frac{1}{f^{2}H^{2}}\langle\Delta_{i}\dot{\Delta}_{j}\dot{\Delta}_{k}\rangle,$
$\displaystyle\dot{U}^{mn}(\bm{q})=\frac{1}{fH}\langle\delta^{m}_{1}\delta^{n}_{2}\dot{\Delta}_{i}\rangle,$
(4.25)

where $\dot{A}_{ij}\equiv\dot{A}^{00}_{ij}$ and
$\ddot{A}_{ij}\equiv\ddot{A}^{00}_{ij}$.

Finally, we would like to highlight that our EFT model is not exhaustive, as
it does not encompass all PT corrections and counterterms up to 1-loop.
Nonetheless, we adopt this simplified approach for two main reasons: firstly,
it has been extensively validated in prior research [79]; and secondly, we
demonstrate in this work that it satisfactorily meets the requirements of DESI
Y1 data. Nevertheless, complete LPT/EFT frameworks are available in the
literature [99, 100].

###  4.3 Inferring cosmological parameters from compressed parameters

As depicted in the bottom row of the schematic in the right panel of Figure 3,
our compressed methodology incorporates an additional step to derive
cosmological parameters from our final compressed parameters. In this section,
we briefly summarize the methodology employed to achieve this.

Let $\Theta_{\texttt{obs}}$ represent our final set of compressed parameters.
We aim to use them to determine a set of cosmological parameters
$\Omega=\\{\omega_{m},A_{s},h,...\\}$. For a given value of these parameters
we can compute a set of compressed parameters $\Theta_{\Omega}$ from theory
following the equations presented in sections 4.1 and 4.2. We use these
equations to run an MCMC chain that finds the regions of the parameter space
of $\Omega$ for which $\Theta_{\Omega}$ is a good fit of
$\Theta_{\texttt{obs}}$. Let $C_{\Theta}$ be the covariance of a given
compressed parameter, and then if we assume a Gaussian likelihood888A Gaussian
likelihood is a good approximation when the distributions of $\Theta_{obs}$
are gaussian distributed, what can be verified in the posteriors of the
compressed parameters for our MCMC exploration, the $\chi^{2}$ is defined as:

$\chi^{2}(\Theta_{\Omega}-\Theta_{\texttt{obs}})=(\Theta_{\Omega}-\Theta_{\texttt{obs}})^{T}C_{\Theta}^{-1}(\Theta_{\Omega}-\Theta_{\texttt{obs}}),$
(4.26)

As mentioned earlier, we employ two compression methodologies, namely the
Standard and Shape-fit methods. In our Standard approach, $\Theta_{\Omega}$ is
as follows:

$\Theta_{\Omega}=\\{D_{A}(z,\Omega)/r_{d}(\Omega),D_{H}(z)/r_{d}(\Omega),f\sigma_{8}(z,\Omega)\\}$
(4.27)

where the equations used for computing the compressed variables are 4.5, 4.2
and 4.3, and we have explicitly included their dependence on the cosmological
models into our notation.

Within our Shape-Fit methodology, in addition to the equations of the Standard
approach, we also need to consider the equations for the additional parameters
n and m, thus $\Theta_{\Omega}$ are given by:

$\Theta_{\Omega}=\\{D_{A}(z,\Omega)/r_{d}(\Omega),D_{H}(z)/r_{d}(\Omega),f\sigma_{s8}(z,\Omega),m(\Omega),n(\Omega)\\},$
(4.28)

where m, n, and $\sigma_{s8}$ are respectively given by equations 4.2, 4.10,
and 4.11.

Also in the case of ShapeFit, we have to keep in mind that the $m$ parameter
changes the tilt and the amplitude of the power spectrum, which means that we
are varying $\sigma_{s8}$. This deformation of the power spectrum is captured
through [63]:

$\begin{split}f\sigma_{s8}&=f\frac{\left(f\sigma_{s8}\right)^{\text{ref
}}}{\left(f\left(P_{\operatorname{lin}}\left(k_{p}\right)\right)^{1/2}\right)^{\text{ref
}}}\left(\frac{1}{s^{3}}P_{\text{lin }}\left(k_{p}/s\right)\right)^{1/2}\\\
&=\frac{\left(f\sigma_{s8}\right)^{\mathrm{ref}}}{\left(fA_{sp}^{1/2}\right)^{\mathrm{ref}}}(fA_{sp}^{1/2})\times\exp\left(\frac{m}{2a}\tanh\left(a\ln\left(\frac{r_{d}^{\mathrm{ref}}}{8}\right)\right)\right)\\\
&\approx\frac{(fA_{sp}^{1/2})}{\left(fA_{sp}^{1/2}\right)^{\mathrm{ref}}}\left(f\sigma_{s8}\right)^{\mathrm{ref}}\end{split}$
(4.29)

where

$A_{sp}=s^{-3}P_{\text{no-wiggle }}^{\text{lin
}}\left(k_{p}/s,\bm{\Theta}\right),$ (4.30)

with $\bm{\Theta}$ the cosmological parameters of the cosmology of reference
or the new one. The last step of equation 4.29 can be done thanks to the fact
that we only use the amplitude at the pivot scale during the MCMC exploration.

###  4.4 Full Modeling in Configuration space

In the Full Modeling approach, which is illustrated in the left side of Figure
3, we use a non-fixed template, which consists of modifying the shape of the
linear power spectrum through the variation of the cosmological parameters
$\omega_{\mathrm{cdm}}$, $\omega_{b}$, $h$, $A_{s}$, $n_{s}$, $N_{\text{eff
}}$, $\Omega_{\mathrm{ncdm}}$, parameters that also depends on early-time
physics, in contrast to the Standard approach. This methodology involves using
a Boltzmann solver code like CLASS or CAMB, and, computing the 1-loop
corrections for the perturbative approach used, in each step of the MCMC
exploration. Despite its higher computational cost, this allows us to extract
cosmological information not only from late-time physics but also from early-
time physics (in contrast with Standard methodology). Hence, the Full Modeling
approach, at least in principle, can impose tighter constraints than the
Standard methodology as it captures more information. As mentioned previously,
Shapefit also looks for extract additional information from early-time physics
through the additional parameters.

However, the high dimensionality of the parameter space requires the use of
new tools to speed up the template generation. In [42], we presented a neural
network approach that enhanced the efficiency of convergence time during MCMC
exploration dramatically. As clustering measurements employ a cosmology of
reference to transform redshift into distance measurements, we need to add two
distortion parameters that compensate the Alcock-Paczynski effect, denoted
them as $q_{||}$ and $q_{\perp}$, and following a similar description of 4.1
with the difference that we do not introduce the BAO scale at the drag epoch,
$r_{\text{d}}$, as its effect is already included in the variation of the
linear power spectrum.

The geometric distortion parameters, perpendicular, and parallel to the line
of sight, in the Full Modeling approach are defined as

$q_{\perp}\left(z\right)=\frac{D_{A}\left(z\right)}{D_{A}^{\mathrm{ref}}\left(z\right)},\quad
q_{\|}\left(z\right)=\frac{H^{\mathrm{ref}}\left(z\right)}{H\left(z\right)},$
(4.31)

here $D_{A}$ is the angular diameter distance, $H$ is the Hubble parameter,
and the $\mathrm{ref}$ superscript indicates that the estimate is done in the
reference or fiducial cosmology of the data multipoles.

The Full Modeling methodology offers the advantage of capturing a broad range
of dynamic-times, as it does not need to fix the early-time dynamics of the
universe. However, it requires the adoption of a specific cosmological model,
and for exploring alternative models of interest, a complete rerun of the MCMC
is needed.

In Table 1 we show the varied parameters and their priors for our baseline
Full Modeling analyses. We keep fixed the slope $n_{s}=0.97$, the effective
number of relativistic degrees of freedom $N_{\mathrm{eff}}=3.046$, the
massive neutrino abundance $\omega_{\mathrm{ncdm}}=0.00064$ and we impose a
big-bang nucleosynthesis (BBN) prior $\mathcal{N}[\mu=0.02237,\sigma=0.00037]$
to $\omega_{b}$ [101].

###  4.5 GSM-EFT code for compressed and direct fits: Users guide

We present our GSM-EFT C-code that computes the multipoles of the one-loop GSM
two-point correlation function in redshift
space999https://github.com/alejandroaviles/gsm.. The GSM-EFT code receives as
an input, the linear power spectrum, obtained from a Boltzmann Solver like
CAMB or CLASS, as well as a set of nuisance parameters: including biases
$\\{b_{1}$, $b_{2}$, $b_{s^{2}}$ $\\}$, the EFT parameters
$\\{\sigma_{\text{EFT}}^{2}$, $c_{1}^{\text{EFT}}\\}$ and the cosmological
parameter $\Omega_{m}$101010$\Omega_{m}$ is necessary to calculate the growth
rate $f$ at the output redshift. The code produces as GSM-EFT outputs the
correlation function multipoles $\ell=0,2,4$ in redshift space.

  * •

GSM-CosmoExtractor. The repository includes a python wrapper and has two
distinct approaches implemented for extracting the cosmological information
explored in this work (developed in section 4). The wrapper also includes a
function to compute the CAMB/CLASS linear power spectrum required as an input
by our GSM-EFT code. The methodology for modeling the effects of the re-
scaling parameters $\\{\alpha_{||},\alpha_{\perp}\\}$ (Eqs. 4.4) required for
compressed parameters, as well as the re-scaling parameters in Full Modeling,
the $\\{q_{||},q_{\perp}\\}$ are also implemented in the python wrapper.

  * •

GSM-EFT modules The GSM-EFT code is composed of two separate modules: CLPT and
GSM. The CLPT module computes the real-space correlation function, pairwise
velocity, and velocity dispersion. Then, the GSM module receives the CLPT
output statistics as inputs along with values of the nuisance parameters (bias
and EFT terms), and outputs the redshift-space two-point correlation function.
In Section 4.2.1 below, we describe the implementation of the compressed
approaches in the wrapper as well as in the GSM-EFT code. The GSM-EFT code
takes only $\sim 1$ second in computing the redshift multipoles, using either
the Full Modeling or ShapeFit  approaches. While this is a relatively small
time, given the sheer number of evaluations required for our MCMC analysis and
the large number of analyses we require, there is an incentive to accelerate
the model evaluation time even further. As introduced in section 3.2, we
generated an emulator to reduce the processing time to around $0.015$ seconds
per evaluation.

  * •

ShapeFit implementation in GSM-EFT. As indicated in Figure 3, the
implementation of ShapeFit in GSM-EFT code implies a modification of the CLPT
module. This is done through two modifications to the original pipeline:
First, the k-funtions that depend on the $P_{L}(k)$ and loop corrections, are
pre-computed using the modified linear power spectrum $P_{\texttt{ref}}(k)$.
Then, the q-functions are computed from k-functions, this is done trough the
re-scaling factor
$\left(\frac{P_{\mathrm{ref}}^{\prime}(k)}{P_{\mathrm{ref}}(k)}\right)^{2}$.

###  4.6 Likelihood and Priors

Our exploration of the parameter space is done through the Bayes law, with the
Bayesian evidence neglected. For a given model, the posterior probability of a
point in the parameter space, denoted as $\mathcal{P}(\bm{\theta})$, is
defined as follows:

$\mathcal{P}(\bm{\theta}\mid\bm{D})\propto\mathcal{L}(\bm{D}\mid\bm{\theta})\times\mathcal{\pi}(\bm{\theta}),$
(4.32)

where, $\mathcal{L}(\bm{D}\mid\bm{\theta})$ is the likelihood and
$\mathcal{\pi}(\bm{\theta})$ is the prior distributions. If we assume that the
parameter space follows a normal distribution, we can express the likelihood
as,

$\mathcal{P}(\bm{D}\mid\bm{\theta})\propto(\chi^{2})^{\frac{\nu-2}{2}}\exp\left(-\frac{\chi^{2}}{2}\right),$
(4.33)

where $\nu$ is the number of degrees of freedom, and $\chi^{2}$ is defined as
$\chi^{2}=(\bm{m}-\bm{d})^{T}C^{-1}(\bm{m}-\bm{d}),$ and the vectors $\bm{m}$
and $\bm{d}$ represent the model and data respectively. The covariance matrix
of the data is denoted as $C$.

The exploration of parameter space is done using the publicly available MCMC
code emcee [102] which is an ensemble sampler with affine invariance defined
in [103]. This code has some advantages, like enabling parallelization,
allowing us to take advantage of the computational power of large clusters. In
our work, we employ the priors outlined in Table 1 for all of our three
independent methodologies. In all the cases we fix the spectral index
$n_{s}=0.97$ unless stated otherwise, the effective number of relativistic
degrees of freedom $N_{\mathrm{eff}}=3.046$, and the massive neutrino
abundance $\omega_{\mathrm{ncdm}}=0.00064$.

To guarantee convergence, we employ an integrated auto-correlation time, which
is monitored every 100 steps [103]. Convergence is considered achieved when
two conditions are simultaneously satisfied: the chain’s length surpasses 100
times the estimated auto-correlation time, and the alteration in this
estimation remains below $1\%$.

Free parameters and priors
---
Cosmological |  | Compressed |
$\,h$ | $\mathcal{U}[0.55,0.91]$ | $f\sigma_{8}$ | $\mathcal{U}[0,1]$
$\,\omega_{b}$ | $\quad\mathcal{N}[0.02237,0.00037]\quad$ | $\alpha_{\perp}$ | $\mathcal{U}[0.5,1.5]$
$\,\omega_{\mathrm{cdm}}$ | $\mathcal{U}[0.08,0.16]$ | $\alpha_{\parallel}$ | $\mathcal{U}[0.5,1.5]$
$\,\log(10^{10}A_{s})$ | $\mathcal{U}[2.0,4.0]$ | $m$ | $\mathcal{U}[-1,1]$
Nuisances
$b_{1}$ | $\mathcal{U}[0,2.0]$ | $b_{2}$ | $\mathcal{U}[-5,10]$
$\sigma^{2}_{\text{EFT}}$ | $\mathcal{U}[-20,100]$ | $c_{1,\text{EFT}}$ | $\mathcal{U}[-100,100]$
$b_{s^{2}}$ | $\mathcal{U}[-10,10]$ |  |
Table 1: Free parameters and priors we use for the baseline Full Modeling
analyses of ABACUS-SUMMIT simulations. The nuisance parameters with their
priors are the same for all three methods. We fix the spectral index of the
power spectrum to the Planck 2019 value $n_{s}=0.9649$.

##  5 Cosmological constraints from configuration-space clustering statistics

In this section, we discuss and compare the results obtained with our pipeline
following three different methodologies: Standard, ShapeFit, and Full
Modeling. However, our attention will primarily be on ShapeFit and Full
Modeling, given their stronger constraining potential. In Section 5.1, we
explore the settings for defining our baseline analysis, in particular the
range of the fits, and the bias settings. In the Section 5.2, we explore
extensions to the baseline analysis, as including the hexadecapole, the
extension of the parameter space including $n_{s}$ and the
$\omega_{0}-\omega_{a}$. Section 5.3 is devoted to exploring the performance
of the analysis for the different DESI tracers and their combination of them.
Finally, in the last section, we compare the constraining power of the two
methodologies, ShapeFit, and Full Modeling, in the cosmological parameter
space.

###  5.1 Baseline Analysis

We start by exploring the behavior of the GSM model for different volumes with
Full Modeling and Shape Fit. Figure 4 shows the posteriors for the LRG tracer
with a minimal scale of $s_{m}in=30\,h^{-1}\,\textrm{Gpc}$ for three different
volumes: $8\,h^{-3}\,\textrm{Gpc}^{3}$, which represents approximately the
volume for DESI first year, $40\,h^{-3}\,\textrm{Gpc}^{3}$ which represents
the volume of DESI by the end of the fifth year, and
$200\,h^{-3}\,\textrm{Gpc}^{3}$ with the total combined volume of the 25
simulations and which represents an unrealistic volume for a survey but allows
us to test models with very small error bars. In Table 2 we report the central
values and 1$\sigma$ errors for the cosmological parameters. As expected, when
we extend to larger volumes, our constraints become tighter, maintaining
approximately the same central values of the posteriors. In Full Modeling all
our constraints tend to fall within the range of $1\sigma$ for $h$
(0.42$\sigma$ for $V_{1}$, 0.61$\sigma$ for $V_{5}$ and 1.14$\sigma$ for
$V_{25}$) and within the range of $1-2\sigma$ for $\mathrm{ln}(10^{10}A_{s})$
(0.46$\sigma$ for $V_{1}$, 0.11$\sigma$ for $V_{5}$ and 1.4$\sigma$ for
$V_{25}$), while $\Omega_{m}$ consistently remains below $1\sigma$
(0.05$\sigma$ for $V_{1}$, 0.12$\sigma$ for $V_{5}$ and 0.28$\sigma$ for
$V_{25}$). On the other hand, in Shape Fit, we obtain constraints within the
range of $1-2\sigma$ $f\sigma_{8}$ (0.04$\sigma$ for $V_{1}$, 0.50$\sigma$ for
$V_{5}$ and 1.7$\sigma$ for $V_{25}$) and below $1\sigma$ for
$\alpha_{\parallel}$ (0.28$\sigma$ for $V_{1}$, 0.48$\sigma$ for $V_{5}$ and
0.94$\sigma$ for $V_{25}$), while $\alpha_{\perp}$ remains around $1-2\sigma$
(0.44$\sigma$ for $V_{1}$, 0.97$\sigma$ for $V_{5}$ and 2.42$\sigma$ for
$V_{25}$) and lastly, for $m$ we obtain constrains around $1\sigma$
(0.03$\sigma$ for $V_{1}$, 0.41$\sigma$ for $V_{5}$ and 1.14$\sigma$ for
$V_{25}$).

The increasing shift for the largest volume can be understood by considering
the limitations of our 1-loop corrections, however even for the largest volume
we do not detect a significant shift in any parameter. Considering the
limitations in accuracy of N-body simulations that we start being sensitive to
at these big volumes [e.g. 104], whose systematics potentially could be
indistinguishable from inaccuracies in the modeling, and given the interest of
comparing the errors with the current DESI year one analysis, we decided
throughout the rest of this work, to report the results for a volume of
$8\,h^{-3}\,\textrm{Gpc}^{3}$ (unless otherwise stated).

Volume (Max. Freedom).
---
Full Modeling
Volume | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$ |
$V_{1}=8h^{-3}\mathrm{Mpc^{3}}$ | $0.3147\pm 0.0099$ | $0.6771\pm 0.0083$ | $2.98\pm 0.12$ |
$V_{5}=40h^{-3}\mathrm{Mpc^{3}}$ | $0.3146\pm 0.0046$ | $0.6768\pm 0.0052$ | $3.003\pm 0.058$ |
$V_{25}=200h^{-3}\mathrm{Mpc^{3}}$ | $0.3146\pm 0.0021$ | $0.6776\pm 0.0035$ | $2.999\pm 0.026$ |
Shape Fit
Volume | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | $m$
$V_{1}=8h^{-3}\mathrm{Mpc^{3}}$ | $0.449\pm 0.025$ | $1.005\pm 0.018$ | $0.9965\pm 0.0079$ | $-0.001\pm 0.031$
$V_{5}=40h^{-3}\mathrm{Mpc^{3}}$ | $0.444\pm 0.012$ | $1.0037\pm 0.0077$ | $0.9960^{+0.0039}_{-0.0043}$ | $0.007\pm 0.017$
$V_{25}=200h^{-3}\mathrm{Mpc^{3}}$ | $0.4408\pm 0.0054$ | $1.0033\pm 0.0035$ | $0.9954\pm 0.0019$ | $0.0107^{+0.0088}_{-0.0099}$

Table 2: Summary of the constraints obtained when fitting the ABACUS LRG
sample with our Full Modeling and Shape Fit methodology for
$s_{\text{min}}=30h^{-1}\mathrm{Mpc}$ using different volumes.

.

Figure 4: The triangular plot shows the posterior distribution of our EFT-GSM
templates from section 3, obtained when we fit the mean of 25 ABACUS-SUMMIT
simulations using a different rescaling number for the covariance matrix. On
the left, we have the fit resulting from Full-Modeling, while on the right, we
have the Shape Fit. The non-rescaled covariance matrix is in maroon, the
covariance matrix rescaled by 5 in green, and the covariance matrix rescaled
by 25 in blue. All the plots are obtained for the LRG tracer with
$s_{min}=30h^{-1}\mathrm{Gpc}$.

####  5.1.1 Optimizing the fitting range

In what follows, we present the results of our fits to the mean of the 25
ABACUS-SUMMIT simulations, employing the methodologies detailed in section 4.
We define our baseline methodology as using the monopole and quadrupole of the
correlation function, using the covariance for a volume of
$8\,h^{-3}\,\textrm{Gpc}^{3}$, with a maximum scale of
$s_{max}=130h^{-1}\,\textrm{Mpc}$ for the three methodologies. For the tests
of the baseline settings, we focus on the LRG sample. We explore in this
section the performance of the GSM-EFT model for different minimum scales for
all methodologies Standard analysis, ShapeFit and Full Modeling.

In the Standard analysis, we constrain $f\sigma_{8}$, $\alpha_{\perp}$ and
$\alpha_{\parallel}$ parameters, introduced in the section 4.1, and, given
that the reference template is computed based on the ABACUS-SUMMIT cosmology,
we expect $\alpha_{\perp}=1$ and $\alpha_{\parallel}=1$. In the ShapeFit
methodology, we also include the $m$ parameter, which, we expect to be $0$. In
Figure 5 we show the constraints for both methodologies in our baseline
analysis. We compare the results of the Standard methodology (left) and of the
ShapeFit (right) for three different minimum scale $s_{min}$ equal to:
$22\,h^{-1}\,\textrm{Mpc}$, $30\,h^{-1}\,\textrm{Mpc}$, and
$38\,h^{-1}\,\textrm{Mpc}$. In Table 3 we present the quantitative results. We
observe that all the results are in agreement with the cosmology of the
ABACUS-SUMMIT within $1\sigma$. We observe a similar performance among these
two methodologies for the parameters in common, which is something we expect,
given the parameter $m$ is intended to provide independent information. This
is the extra information that $m$ captures about the slope of the linear power
spectrum, which makes ShapeFit as constraining as Full-Modeling when
interpreting the set of compressed parameters from ShapeFit into the
cosmological parameters of $\Lambda$CDM, as shown in [94, 95].

Figure 5: The triangular plot shows the posterior distribution obtained from
the mean of 25 ABACUS-SUMMIT simulations using a non-rescaled covariance
matrix for the LRG tracer. On the left side, we have constraints for the
Standard methodology, while on the right side, we have constraints for the
ShapeFit methodology. We plot three different ranges for the LRG tracers, a
minimum scale $s_{min}=22\,h^{-1}\mathrm{Mpc}$ in red,
$s_{min}=30\,h^{-1}\mathrm{Mpc}$ in green, and
$s_{min}=38\,h^{-1}\mathrm{Mpc}$ in blue. Baseline
---
ShapeFit for Maximal Freedom
$s_{\text{min}}$ | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | $m$
$22\,h^{-1}\mathrm{Mpc}$ | $0.457^{+0.021}_{-0.025}$ | $1.010\pm 0.017$ | $0.9968\pm 0.0074$ | $-0.006^{+0.025}_{-0.022}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.449\pm 0.025$ | $1.005\pm 0.018$ | $0.9965\pm 0.0079$ | $-0.001\pm 0.031$
$38\,h^{-1}\mathrm{Mpc}$ | $0.448\pm 0.028$ | $1.006\pm 0.019$ | $0.9951\pm 0.0083$ | $0.002\pm 0.034$
Standard for Maximal Freedom
$22\,h^{-1}\mathrm{Mpc}$ | $0.456\pm 0.021$ | $1.009\pm 0.016$ | $0.9960\pm 0.0074$ | -
$30\,h^{-1}\mathrm{Mpc}$ | $0.447\pm 0.022$ | $1.006\pm 0.017$ | $0.9961\pm 0.0076$ | -
$38\,h^{-1}\mathrm{Mpc}$ | $0.452\pm 0.026$ | $1.005\pm 0.018$ | $0.9960^{+0.0079}_{-0.0087}$ | -
Table 3: Constraints on our compressed parameters when fitting our ABACUS LRG
sample for three minimum scales when using Standard and ShapeFit methodology.

In Full Modeling we define as our baseline analysis, four free cosmological
parameters, $\\{\Omega_{m}$, $h$, $\mathrm{ln}(10^{10}A_{s})$, $\omega_{b}$ },
however as we explained in the methodology we impose a Gaussian BBN prior on
$\omega_{b}$, as well as flat priors in the nuisance parameters which are
described in Table 1. In Figure 7 we plot the posterior distribution contours
for our baseline analysis varying $s_{{min}}$ for three different values with
a bin separation of $8\,h^{-1}\mathrm{Mpc}$, on the right side we have a
triangular plot for the Full Modeling methodology while on the left side, we
have the ShapeFit methodology. We show the values in Table 4. We notice that
both methodologies are in good agreement with the true values of the
cosmological parameters for all scales, within $\leq 0.1\sigma$ for
$\Omega_{m}$, $\leq 0.6\sigma$ for $h$ and, $\leq 0.7\sigma$ for $A_{s}$ when
$s_{min}=22\,h^{-1}\mathrm{Mpc}$ for Full modeling and ShapeFit. For ShapeFit
we observe a similar trend, all deviations from the true value are $\leq
0.2\sigma$ for $\Omega_{m}$, $\leq 0.4\sigma$ for $h$ and $\leq 0.5\sigma$ for
$\mathrm{ln}(10^{10}A_{s})$. We conclude, we do not detect any significant
shift with either of the methodologies employed. In Figure 6, we illustrate
the shape of the correlation function. This depiction represents the mean of
the values derived from the fit of the baseline analysis, with
$s_{min}=30,h^{-1}\mathrm{Mpc}$, and a $\chi^{2}=60.57$ for 58 bins and 9 free
parameter.

Figure 6: This plot shows the correlation function obtained from the mean of
25 ABACUS-SUMMIT simulations and the mean of the posterior distribution in the
range $s_{min}=30\,h^{-1}\mathrm{Mpc}$ and $s_{max}=130\,h^{-1}\mathrm{Mpc}$.
The shadow regions show the error bars from the non-rescaled covariance matrix
for the LRG tracer.

Figure 7: The triangular plot shows the posterior distribution obtained from
the mean of 25 ABACUS-SUMMIT simulations using a non-rescaled covariance
matrix for the LRG tracer following the Full Modeling methodology on the left
side and the ShapeFit methodology on the right side. We plot three different
ranges for the LRG tracers, a minimum scale $s_{min}=22\,h^{-1}\mathrm{Mpc}$
in red, $s_{min}=30\,h^{-1}\mathrm{Mpc}$ in green, and
$s_{min}=38\,h^{-1}\mathrm{Mpc}$ in blue.

As we can notice from the results in this section, although we obtain results
below $1\sigma$ for all the choices of minimum scales, the results with
$s_{min}=22\,h^{-1}\mathrm{Mpc}$ provide the smallest error bars without any
detectable bias, which is in agreement with what other works have concluded
[42]. Therefore, our optimal fitting range is
$s_{\text{range}}=22-130\,h^{-1}\mathrm{Mpc}$. Several crucial results will be
exclusively presented with this range, unless otherwise specified. Our
baseline analysis is extended and compared in further sections.

BASELINE
---
Full Modeling (Max Freedom)
$V$ | $s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
$V_{1}$ | $22\,h^{-1}\mathrm{Mpc}$ | $0.3138\pm 0.0081$ | $0.6760\pm 0.0074$ | $3.056\pm 0.093$
$V_{1}$ | $30\,h^{-1}\mathrm{Mpc}$ | $0.3147\pm 0.0099$ | $0.6771\pm 0.0083$ | $2.98\pm 0.12$
$V_{1}$ | $38\,h^{-1}\mathrm{Mpc}$ | $0.316\pm 0.011$ | $0.6786\pm 0.0086$ | $2.95\pm 0.13$
$V_{5}$ | $22\,h^{-1}\mathrm{Mpc}$ | $0.3146\pm 0.0037$ | $0.6775\pm 0.0042$ | $3.038\pm 0.044$
$V_{5}$ | $30\,h^{-1}\mathrm{Mpc}$ | $0.3146\pm 0.0046$ | $0.6768\pm 0.0052$ | $3.003\pm 0.058$
$V_{5}$ | $38\,h^{-1}\mathrm{Mpc}$ | $0.3144\pm 0.0050$ | $0.6771\pm 0.0052$ | $2.984\pm 0.067$
ShapeFit  (Max Freedom)
$V$ | $s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
$V_{1}$ | $22\,h^{-1}\mathrm{Mpc}$ | $0.3140\pm 0.0055$ | $0.6747\pm 0.0082$ | $3.09\pm 0.11$
$V_{1}$ | $30\,h^{-1}\mathrm{Mpc}$ | $0.3144\pm 0.0072$ | $0.6758\pm 0.0088$ | $3.04\pm 0.12$
$V_{1}$ | $38\,h^{-1}\mathrm{Mpc}$ | $0.3149\pm 0.0081$ | $0.6768\pm 0.0086$ | $3.03\pm 0.13$
$V_{5}$ | $22\,h^{-1}\mathrm{Mpc}$ | $0.3156\pm 0.0027$ | $0.6772\pm 0.0043$ | $3.044\pm 0.053$
$V_{5}$ | $30\,h^{-1}\mathrm{Mpc}$ | $0.3153\pm 0.0036$ | $0.6760\pm 0.0054$ | $3.014\pm 0.061$
$V_{5}$ | $38\,h^{-1}\mathrm{Mpc}$ | $0.3151\pm 0.0039$ | $0.6771\pm 0.0055$ | $2.999\pm 0.069$
Table 4: Constraints on our cosmological parameters when fitting our ABACUS
LRG sample for three minimum scales when using different methodologies, Full
Modeling methodology and ShapeFit methodology together with two different
volumes $V_{1}=8\,h^{-3}\mathrm{Mpc}^{3}$ and
$V_{5}=40\,h^{-3}\mathrm{Mpc}^{3}$.

####  5.1.2 Maximal and minimal freedom

We explore two different parameter space configurations of the galaxy bias
terms, that we refer to as maximal and minimal freedom. In the first, we let
all parameters vary simultaneously and it is considered as our baseline
configuration. While the second configuration assumes a co-evolution for which
non-local biases appear only through gravity evolution, and so are not present
in Lagrangian coordinates scheme, corresponding to fix $b_{s^{2}}=0$.

Figure 8 shows the mean value of the posterior distribution of the parameters
found when running our methodology using both maximal and minimal freedom
configurations and using both ShapeFit and Full Modeling methods. The
corresponding quantitative results are presented in Tables 3, 4, 5 and 6. We
notice that in both our it Full Modeling and it Shape Fit results, the central
values are fairly consistent with each other. However, the Minimal freedom
configuration leads to $\sim 10\%$ smaller error bars, this trend is expected
as we reduce the number of free parameters. We notice that the it Full
Modeling scenario shows an underestimation of the $A_{s}$ parameters in this
larger scale which is closer to its expected value at
$22\,h^{-1}\mathrm{Mpc}$.

On the other hand, if we explore the compressed parameter space for the
$22h^{-1}\mathrm{Mpc}$ scale specifically, we observe that the Maximal freedom
mean value predictions are closer to the true value of the simulation for
$f\sigma_{8}$, while the Minimal freedom case fares better with its
predictions for $\alpha_{\parallel}$. Additionally, we show that the error
bars are of similar size for all parameters except $m$, where the Maximal
freedom case leads to tighter constraints.

Full Modeling (Min Freedom).
---
$s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
$22\,h^{-1}\mathrm{Mpc}$ | $0.3150\pm 0.0072$ | $0.6755\pm 0.0072$ | $3.058\pm 0.065$
$30\,h^{-1}\mathrm{Mpc}$ | $0.3139\pm 0.0082$ | $0.6767\pm 0.0078$ | $2.989\pm 0.082$
$38\,h^{-1}\mathrm{Mpc}$ | $0.314\pm 0.010$ | $0.6786\pm 0.0080$ | $2.96\pm 0.11$
Table 5: Constraints to our parameters when fitting our ABACUS LRG sample for
three minimum scales when using different methodologies, Full Modeling
methodology. ShapeFit (Min Freedom).
---
$s_{\text{min}}$ | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | $m$
$22\,h^{-1}\mathrm{Mpc}$ | $0.453\pm 0.020$ | $1.009\pm 0.017$ | $0.9967\pm 0.0076$ | $-0.004^{+0.028}_{-0.018}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.442\pm 0.021$ | $1.005\pm 0.017$ | $0.9963\pm 0.0077$ | $-0.001\pm 0.027$
$38\,h^{-1}\mathrm{Mpc}$ | $0.440\pm 0.024$ | $1.005\pm 0.018$ | $0.9945\pm 0.0080$ | $0.004\pm 0.032$
Table 6: Constraints to our parameters when fitting our ABACUS LRG sample for
three minimum scales when using different methodologies for ShapeFit
methodology.

Figure 8: Error bar comparison for LRG using a volume of
$8\,h^{-3}\mathrm{Gpc}^{3}$. The red bars were obtained using maximal freedom,
while the red ones were obtained using minimal freedom. We show the model’s
behavior at various minimum scales, each one separated by 1 bin of 8
$h^{-1}\mathrm{Mpc}$. On the left side, we have the ShapeFit fits, while on
the right side, we have the Full Modeling fits.

In conclusion, our comparison between the predicted central values and the
expected values of our parameters reveals similar performances between our
Minimal and Maximal freedom configurations. This holds true for both the
compressed parameter space and the cosmological parameter space. We observe no
significant shift of the fitted parameters from their true value. However, our
baseline analysis with $s_{min}=22\,h^{-1}\mathrm{Mpc}$ indicates a
improvement in most parameters of approximately $\sim 10\%$ in the size of the
error bars when using the more restricted Minimal freedom configuration, this
is true regardless of whether we employ the compressed or Full Modeling
approach. However, throughout this work we opt for a conservative approach and
designate the Maximal freedom scenario as our baseline configuration.

###  5.2 Extensions to Baseline Analysis

We will present three extensions to the main analysis we presented above.
First we focus on the information we can extract from the hexadecapole of the
correlation function, then we explore the possibility to constrain the
spectral index $n_{s}$ and finally we allow the Dark Energy equation of state
evolve with the time using the $(w_{o},w_{a})$ prescription.

####  5.2.1 Effect of including hexadecapole

In this subsection, we explore the effect of adding the hexadecapole in our
baseline analysis. It is well known that the hexadecapole is considered to be
the multipole where PT at first-loop starts to struggle to reproduce the small
scales for both the 2-point correlation function and the power spectrum,
particularly at small scales (below $30\,h^{-1}\mathrm{Mpc}$ in configuration
space), where non-linear effects start to be important.

In Figure 9 we can observe the effect of including the hexadecapole in the
compressed parameter space and the cosmological parameter space for three
different ranges, $s_{min}=22$, $30$, and $38\,h^{-1}\mathrm{Mpc}$. The fits
following our baseline analysis are plotted in red, while fits adding
hexadecapole are shown in green. Due to the difficulties in small scale model
of the hexadecapole, we remove the bin $22\,h^{-1}\mathrm{Mpc}$ and maintain
the bin $s_{min}=22\,h^{-1}\mathrm{Mpc}$ for the monopole and the quadrupole.
This is done because when we include the smaller hexadecapole bin, we note
that our posterior distribution becomes bimodal due to the appearance of a
second region of the parameter space with high likelihood. However, this
region has already been discarded by Planck to more than $3\sigma$. Given that
we don’t trust our model at these smaller $r$ scales, we opt to exclude this
bin. The quantitative results are shown in Tables 3, 4, 7 and 8. In the
compressed parameters, we note slight variations within the 1$\sigma$ range
for $f\sigma_{8}$, $\alpha_{\parallel}$, $\alpha_{\perp}$, and $m$. Similarly,
the cosmological parameters exhibit consistent results, with $\Omega_{m}$,
$h$, and $A_{s}$ all falling within the 1$\sigma$ range.

On the other side, the error bars when we include the hexadecapole for the
three minimum scales have a small reduction in ShapeFit. For $f\sigma_{8}$ the
error bar is reduced by $0.26\%-0.07\%$ compared with the baseline analysis,
for $\alpha_{\parallel}$ this change is $0.26\%-0.30\%$, while
$\alpha_{\perp}$ remains in $0.06\%$, and finally for $m$ this change goes
from an increase of $0.04\%$ to a reduction of $0.08\%$. For Full Modeling we
have a change in the error bars for $\Omega_{m}$ that goes from an increment
of $0.13\%$ to a decrease of $0.17\%$, for $h$ we have a diminution in the
error bars of $0.03\%-0.14\%$, while $\mathrm{ln}(10^{10}A_{s})$ has a
reduction in the error bars of $0\%-0.21\%$.

In conclusion, adding the hexadecapole does not significantly improve the
results, even when removing specifically the lowest bin distance in order to
prevent the problems on quasi-linear scale model. The main effect is slightly
reducing the error on the Full Modeling analysis. However, we can notice that
all the constraints remain below 1$\sigma$ from the true value when we include
the hexadecapole.

Full Modeling (M+Q+H)
---
$s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
$22\,h^{-1}\mathrm{Mpc}$ | $0.3138^{+0.0099}_{-0.0085}$ | $0.6707\pm 0.0072$ | $3.098^{+0.066}_{-0.082}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.3094\pm 0.0082$ | $0.6711\pm 0.0071$ | $3.060\pm 0.095$
$38\,h^{-1}\mathrm{Mpc}$ | $0.3107\pm 0.0093$ | $0.6741\pm 0.0075$ | $3.01\pm 0.13$
Table 7: 1D constraints when fitting the ABACUS LRG mocks with (M+Q+H) and
without (M+Q) the hexadecapole for three minimum scales for Full Modeling.
ShapeFit (M+Q+H)
---
$s_{\text{min}}$ | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | $m$
$22\,h^{-1}\mathrm{Mpc}$ | $0.459\pm 0.017$ | $1.009^{+0.011}_{-0.012}$ | $1.0015^{+0.0074}_{-0.0065}$ | $-0.013^{+0.026}_{-0.023}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.453\pm 0.022$ | $1.010^{+0.012}_{-0.013}$ | $0.9984\pm 0.0074$ | $-0.017\pm 0.025$
$38\,h^{-1}\mathrm{Mpc}$ | $0.451\pm 0.026$ | $1.011\pm 0.014$ | $0.9965\pm 0.0078$ | $-0.007\pm 0.031$
$30\,h^{-1}\mathrm{Mpc}$ | $0.449\pm 0.025$ | $1.005\pm 0.018$ | $0.9965\pm 0.0079$ | $-0.001\pm 0.031$
Table 8: 1D constraints when fitting the ABACUS LRG mocks with (M+Q+H) and
without (M+Q) the hexadecapole for three minimum scales for ShapeFit
methodology.

Figure 9: 1D constraints as a function of minimum scale for the LRGs using a
volume of 8 $h^{-1}\mathrm{Mpc}$. The results for the baseline analysis
without the hexadecapole are shown in red while the results including the
hexadecapole are shown in green. On the left side, we have the ShapeFit fits,
while on the right side, we have the Full Modeling fits. The lower bin,
represented in blue, has the same settings, but the bin
$s_{min}=22\,h^{-1}\mathrm{Mpc}$ is excluded for the hexadecapole.

####  5.2.2 Exploring the impact of $n_{s}$ free

The scalar spectral index, denoted as $n_{s}$, serves as a key parameter that
characterizes the primordial power spectrum of density perturbations in the
early universe. As mentioned in Section 4.6, in our standard approach, we fix
this parameter to its Planck 2018 reported value of $n_{s}=0.9649$ [75].
However, in this section, we incorporate this parameter into our methodology
as one of the extensions of our baseline analysis. This extension enables a
more comprehensive and complex exploration of the parameter space. In this
section, we impose the following prior on $n_{s}$:

$\displaystyle n_{s}$ $\displaystyle:\,\mathcal{U}[0.5,1.5]$ (5.1)

Given the new degree of freedom included, the models become more complex, in
order to maintain the precision of our emulator, we expand the number of
multipoles utilized in the neural network to 90,000. Among these, 80,000 are
allocated for training the network, while 5,000 served as a validation set.
The remaining 5,000 are our test set dedicated to generating statistics for
evaluating the model’s performance.

Figure 10 shows the results of our analysis when varying minimum scales,
$s_{\text{min}}$, while treating $n_{s}$ as a free parameter for the LRG
sample and using Full Modeling. The results are presented for three distinct
minimum scales, and each one has a different performance in the precision at a
parameter is recovered. Our analysis consistently captures the expected values
of all parameters within the $1\sigma$ uncertainty contour, including the
extra parameter $n_{s}$. Specifically, for
$s_{\text{min}}=22\,h^{-1}\mathrm{Mpc}$, we achieve remarkable accuracy in
determining $h$ and $\log(10^{10}A_{s})$ with $0.13\sigma$ and $0.13\sigma$
deviations from the true value respectively. Similarly, when
$s_{\text{min}}=30\,h^{-1}\mathrm{Mpc}$, our precision improves notably for
$\Omega_{m}$, $h$, and $n_{s}$ with $0.13\sigma$, $0.19\sigma$ and
$0.05\sigma$ deviations respectively from the true value. Finally, when
$s_{\text{min}}=38\,h^{-1}\mathrm{Mpc}$, our analysis aligns exceptionally
well with $\Omega_{m}$ and $n_{s}$ within $0.07\sigma$ and $0.04\sigma$
deviations respectively. The resulting constraints on our parameters are
summarized in Table 9 (to be compared with baseline analysis in Table 4).

$n_{s}$ Free
---
$s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$ | $n_{s}$
$22\,h^{-1}\mathrm{Mpc}$ | $0.307\pm 0.014$ | $0.6746\pm 0.0076$ | $3.05^{+0.10}_{-0.12}$ | $0.975\pm 0.037$
$30\,h^{-1}\mathrm{Mpc}$ | $0.313^{+0.015}_{-0.017}$ | $0.6752\pm 0.0086$ | $2.99\pm 0.15$ | $0.967\pm 0.045$
$38\,h^{-1}\mathrm{Mpc}$ | $0.314\pm 0.018$ | $0.6773\pm 0.0089$ | $2.96\pm 0.17$ | $0.967\pm 0.050$

Table 9: Summary of the constraints obtained when fitting the ABACUS LRG
sample with our Full Modeling methodology for our three minimum scales and
when we allow $n_{s}$ to vary as a free parameter

.

This shows how the sensitivity with which we can recover parameters varies
with the minimum ($s_{\text{min}}$) scale when $n_{s}$ is free.The distinct
performance levels observed across different $s_{\text{min}}$ values
underscore the importance of carefully selecting this parameter in
cosmological analyses. We see that, $s_{\text{min}}=30\,h^{-1}\mathrm{Mpc}$ is
the scale exhibiting the most robust performance across all parameters, except
$\mathcal{A}_{s}$ which is still predicted to within less than one sigma from
the true value. Our results highlight the capability of our analysis to
consistently capture parameters within one sigma uncertainty, highlighting its
robustness in characterizing extended cosmological models.

The results presented in this section are built using our Full Modeling
methodology only.

Figure 10: The triangular plot displays the posterior one-sigma and two-sigma
contours of our baseline analysis when $n_{s}$ is treated as a free parameter,
using our Full Modeling methodology on the LRG sample. This plot illustrates
the variation in predictions when adjusting the minimum scales. It is
equivalent to the left panel of Figure 7 but with $n_{s}$ set as a free
parameter.

####  5.2.3 Exploring the impact of $w_{0}$ and $w_{a}$ as free parameters

The parameters $w_{0}$ and $w_{a}$ characterize the behavior of dark energy
within the Chevallier–Polarski–Linder model [105, 106]. This model provides a
framework to study an evolving equation of state for dark energy and is an
extension of the cosmological constant $\Lambda$CDM model. The function
$w(a)=w_{0}+w_{a}(1-a)$ describes the evolution of the equation of state
parameter over cosmic times, offering insights into the dynamic nature of dark
energy throughout the universe’s history.

In this subsection, we analyze two different extensions of our methodology. In
the first, we incorporate variations of the parameter $w_{0}$ only into our
baseline analysis. Then, in the second extension, we vary both $w_{0}$ and
$w_{a}$ and study their combined effect.

$w_{0}$CDM Extension
---
$s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$ | $w_{0}$ | $w_{a}$
$22\,h^{-1}\mathrm{Mpc}$ | $0.309^{+0.015}_{-0.021}$ | $0.681^{+0.022}_{-0.019}$ | $3.031\pm 0.090$ | $-1.031^{+0.077}_{-0.10}$ | $0$
$30\,h^{-1}\mathrm{Mpc}$ | $0.313^{+0.017}_{-0.023}$ | $0.679^{+0.025}_{-0.022}$ | $2.99\pm 0.12$ | $-1.011^{+0.085}_{-0.11}$ | $0$
$38\,h^{-1}\mathrm{Mpc}$ | $0.314^{+0.018}_{-0.024}$ | $0.681^{+0.026}_{-0.022}$ | $2.94\pm 0.13$ | $-1.014^{+0.090}_{-0.12}$ | $0$
$w_{0}w_{a}$CDM Extension
$22\,h^{-1}\mathrm{Mpc}$ | $0.314^{+0.020}_{-0.026}$ | $0.674^{+0.027}_{-0.024}$ | $3.13^{+0.11}_{-0.16}$ | $-0.88\pm 0.25$ | $-0.67\pm 0.97$
$30\,h^{-1}\mathrm{Mpc}$ | $0.319^{+0.022}_{-0.027}$ | $0.673\pm 0.026$ | $3.07^{+0.14}_{-0.19}$ | $-0.86^{+0.29}_{-0.26}$ | $-0.7\pm 1.0$
$38\,h^{-1}\mathrm{Mpc}$ | $0.322^{+0.023}_{-0.029}$ | $0.673\pm 0.027$ | $3.07^{+0.16}_{-0.23}$ | $-0.83^{+0.30}_{-0.27}$ | $-0.8\pm 1.1$
Table 10: Constraints when fitting the ABACUS LRG mock for three minimum
scales when we allow $w_{0}$ an $w_{a}$ to vary as a free parameter and when
using our Full Modeling methodology.

We expand our standard setup, to incorporate $w_{0}$ and $w_{a}$ as additional
parameters, and we use the resulting models to train a new neural network. To
ensure the precision of our emulator, we increase the number of multipoles
used to train and test our neural networks to 170,000 points. Of these,
150,000 are allocated for train the networks, with 10,000 designated for
validation and the remaining 10,000 are our test set that evaluate the model’s
performance. We also apply the following priors to our parameters:

$\displaystyle w_{0}$ $\displaystyle:\,\mathcal{U}[-2,0]$ (5.2) $\displaystyle
w_{a}$ $\displaystyle:\,\mathcal{U}[-3,3]$ (5.3)

Here, we only consider the extension of our Full Modeling methodology.

Let us first present our results for the extended methodology that only
includes $w_{0}$. Figure 11 shows the triangular plot of the constraints on
the set of cosmological parameters, where we can see that our analysis yields
parameter constraints that are all within 1$\sigma$ from their true value. The
$s_{\text{min}}=22\,h^{-1}\mathrm{Mpc}$ case shows better performance when
recovering cosmological parameters with smaller 1$\sigma$contours, while also
being more accurate when predicting $A_{s}$. The model presents a
$0.34\,\sigma$ deviation for the true $\Omega_{m}$ value, $0.41\sigma$ for
$h$, $0.05\,\sigma$ for $\log(10^{10}A_{s})$, and $0.09\,\sigma$ for $w_{0}$.
Our extended model successfully recovers all parameters, including $w_{0}$,
showcasing its robustness when recovering the cosmological parameters.

Figure 12 shows the results for our second extension when we also add $w_{a}$
as a free parameter and our extended model depends on both $w_{0}$ and
$w_{a}$. As with the first extension, our analysis recovers all parameter
constraints to within less than 1 $\sigma$ from the true value for all
parameters. The three models perform similarly well showing similar accuracy
and prediction in most parameters, we note that the model with a minimum scale
of $s_{\text{min}}=30\,h^{-1}\mathrm{Mpc}$ and
$s_{\text{min}}=38\,h^{-1}\mathrm{Mpc}$ slightly perform better than
$s_{\text{min}}=22\,h^{-1}\mathrm{Mpc}$ when recovering $A_{s}$ and
$\Omega_{m}$, but it still has the smaller error bars. The deviations for
$s_{\text{min}}=22\,h^{-1}\mathrm{Mpc}$ are within 1$\sigma$ for this model
and are $0.05\,\sigma$ for $\Omega_{m}$, $0.01\,\sigma$ for $h$,
$0.63\,\sigma$ for $\log(10^{10}A_{s})$, $0.47\,\sigma$ for $w_{0}$, and
$0.69\sigma$ for $w_{a}$.

Table 10 summarizes our final constraints for both extensions and all three
scales. Our extended model successfully recovers all parameters, including
$w_{0}$ and $w_{a}$, showing its ability to reliably recover the cosmological
parameters.

Figure 11: The triangular plot displays the posterior one-sigma and two-sigma
contours of our baseline analysis when $w_{0}$ is treated as a free parameter,
using our Full Modeling methodology on the LRG sample. This plot illustrates
the variation in predictions when adjusting the minimum scales. Figure 12:
This plot shows the error bar comparison for LRG using a volume of 8
$h^{-1}\mathrm{Mpc}$. The red bars were obtained using maximal freedom, while
the red ones were obtained using minimal freedom. We show the model’s behavior
at various minimum scales, each one separated by 1 bin of 8
$h^{-1}\mathrm{Mpc}$. On the left side, we have the ShapeFit fits, while on
the right side, we have the Full Modeling fits.

###  5.3 All DESI tracers: LRG, ELG, QSO mocks

In this section, we expand our analysis to include fits to the ELG and QSO
mocks introduced in Section 2. Additionally, we create a new dataset by
fitting all three samples together, this is done by merging all three distinct
data vectors. Each vector contributes to the total sum of the logarithm of the
likelihood, which is subsequently fitted using our model.

As with our previous sections, we fit three iterations of our model,
additionally, we perform all fits twice, once using our ShapeFit template and
once using the Full Modeling approach.

Tables 11 and 12 show the constraints on all ShapeFit and Full Modeling
parameters for all fits mentioned above. The tables demonstrate that, in
general, all constraints are tighter for our smaller minimum scale of
$s_{\text{min}}=22,h^{-1}\mathrm{Mpc}$. To visually compare these results more
intuitively, Figure 13 summarizes them as a set of error bar plots that also
include our previous LRG results from Section 5.1.1. It is worth noting that
all results in the figure are within one standard deviation from their
expected value, indicating good accuracy. Given that all results are
consistent with their expected value and that our measurements at the smaller
minimum scale are slightly more precise, we select this as our main result.

Figure 14 depicts the triangular plot of our complete shape model alongside
our findings for $s_{\text{min}}=22,h^{-1}\mathrm{Mpc}$. We observe that,
generally, the constraints from LRG and ELG exhibit similar precision, with
ELG being marginally more precise, while QSO constraints are comparatively
less precise overall. This aligns with our expectations based on the different
effective volumes of the samples, with ELG being the largest DESI sample and
QSO being the smallest. Additionally, we observe that our combined fit matches
the expected value of all parameters within less than one sigma, with a
$0.23\,\sigma$ deviation for $\Omega_{m}$, a $0.68\,\sigma$ deviation for $h$,
and a $0.26\,\sigma$ deviation for $\log(10^{10}A_{s})$.

DESI tracers
---
$s_{\text{min}}$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
ELG Full Modeling
$22\,h^{-1}\mathrm{Mpc}$ | $0.314\pm 0.011$ | $0.6779\pm 0.0066$ | $3.015^{+0.085}_{-0.10}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.314\pm 0.011$ | $0.6761\pm 0.0072$ | $3.014\pm 0.092$
$38\,h^{-1}\mathrm{Mpc}$ | $0.316^{+0.012}_{-0.014}$ | $0.6769\pm 0.0076$ | $3.00^{+0.11}_{-0.10}$
QSO Full Modeling
$22\,h^{-1}\mathrm{Mpc}$ | $0.309\pm 0.017$ | $0.681\pm 0.014$ | $3.02\pm 0.14$
$30\,h^{-1}\mathrm{Mpc}$ | $0.310\pm 0.019$ | $0.683\pm 0.015$ | $3.01\pm 0.16$
$38\,h^{-1}\mathrm{Mpc}$ | $0.3139\pm 0.0073$ | $0.6776\pm 0.0057$ | $2.983\pm 0.076$
COMBINED Full Modeling
$22\,h^{-1}\mathrm{Mpc}$ | $0.3139\pm 0.0057$ | $0.6770\pm 0.0050$ | $3.021\pm 0.060$
$30\,h^{-1}\mathrm{Mpc}$ | $0.3129\pm 0.0059$ | $0.6772\pm 0.0057$ | $3.002\pm 0.062$
$38\,h^{-1}\mathrm{Mpc}$ | $0.309^{+0.015}_{-0.021}$ | $0.681^{+0.022}_{-0.019}$ | $3.031\pm 0.090$
Table 11: Constraints to our Full shape (top) and ShapeFit (bottom) parameters
when fitting our ABACUS ELG, QSO and Combined samples for all three minimum
scales explored. DESI tracers
---
$s_{\text{min}}$ | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | m
ELG ShapeFit
$22\,h^{-1}\mathrm{Mpc}$ | $0.411\pm 0.019$ | $1.007\pm 0.014$ | $0.9971\pm 0.0076$ | $-0.001^{+0.023}_{-0.020}$
$30\,h^{-1}\mathrm{Mpc}$ | $0.415^{+0.014}_{-0.016}$ | $1.007\pm 0.014$ | $0.9960\pm 0.0070$ | $-0.004\pm 0.027$
$38\,h^{-1}\mathrm{Mpc}$ | $0.415^{+0.016}_{-0.018}$ | $1.007\pm 0.015$ | $0.9959^{+0.0066}_{-0.0075}$ | $-0.003\pm 0.030$
QSO ShapeFit
$22\,h^{-1}\mathrm{Mpc}$ | $0.385\pm 0.023$ | $1.006\pm 0.026$ | $0.995\pm 0.015$ | $-0.011\pm 0.027$
$30\,h^{-1}\mathrm{Mpc}$ | $0.385\pm 0.025$ | $1.004\pm 0.027$ | $0.994\pm 0.015$ | $-0.006\pm 0.035$
$38\,h^{-1}\mathrm{Mpc}$ | $0.384^{+0.025}_{-0.029}$ | $1.004\pm 0.027$ | $0.994\pm 0.015$ | $0.003\pm 0.038$
Table 12: Constraints to our Full shape (top) and ShapeFit (bottom) parameters
when fitting our ABACUS ELG, QSO, and Combined samples for all three minimum
scales explored.

Figure 13: Error bar plots on LRG (red), ELG (green), and QSO (blue)
constraints for three minimum scales and for Shape Fit (left) and Full
Modeling (right) parameters. The black dashed lines show the true value of the
parameters (note that the true value of $f\sigma_{8}$ is dependent on
redshift). Figure 14: We present results for our three tracers: LRG (green
dashed lines), ELG (red dashed lines) and QSO (salmon dashed lines). We also
include the resulting contour of the combined $LRG+ELG+QSO$ (blue solid lines
and contours), and the expected value of each parameter (gray faded lines).

###  5.4 ShapeFit and Full Modeling methodology comparison

In this section, we transform ShapeFit from compressed parameters to
cosmological parameters to be able to compare with Full Modeling in the same
parameter space. In Figure 15 we present a plot where we can observe small
differences for ShapeFit and Full Modeling when we fit the correlation
function using LRG for a volume of $V=8\,h^{-3}\mathrm{Gpc^{3}}$ for the
scales, $s_{min}=30\,h^{-1}\mathrm{Mpc}$. In particular, we observe a good
agreement for the Maximal Freedom, Minimal Freedom, the case when we include
the hexadecapole and the $n_{s}$ extension having a similar constraining
power.

First, we look at the baseline analysis, which shows good performance as we
can notice from Table 4, where for $\Omega_{m}$, we have nearly the same
results in both methodologies, with a central value difference in the range of
$10^{-3}$ and error bars of $0.68\sigma_{FM}<\sigma_{SF}<0.74\sigma_{FM}$
being ShapeFit more constrictive for the three minimum scales. For $h$, we
have a similar result obtaining a difference in the central value of around
$10^{-3}$ and error bars of $1.11\sigma_{FM}<\sigma_{SF}<1.0\sigma_{FM}$ this
time, shape fit shows less constrictive power. In the last parameter,
$\mathrm{ln}(10^{10}A_{s})$, the central value shows a difference of the order
of $10^{-3}$ with error bars $1.0\sigma_{FM}<\sigma_{SF}<1.18\sigma_{FM}$,
ShapeFit again is less constrictive.

The Minimal Freedom case shows similar performance, the statistic is detailed
in Table 13, where we obtain for $\Omega_{m}$ roughly the same results in the
central value with a difference of the order of $10^{-3}$ and error bars
$0.76\sigma_{FM}<\sigma_{SF}<0.80\sigma_{FM}$ being ShapeFit more
constrictive. For $h$, we have a central value which shows a difference around
$10^{-3}$ and error bars of $1.05\sigma_{FM}<\sigma_{SF}<1.1\sigma_{FM}$, with
ShapeFit with less constrictive power. Finally, $\mathrm{ln}(10^{10}A_{s})$
shows a central value difference of the order of $10^{-2}$ with error bars
$1.09\sigma_{FM}<\sigma_{SF}<1.33\sigma_{FM}$, again ShapeFit less
constrictive.

When we include the hexadecapole in the baseline analysis, we obtain similar
constraining power in both analyses. In Table 13, we can see that for
$\Omega_{m}$, we have a central value with a difference around $10^{-3}$ and
error bars of $0.71\sigma_{FM}<\sigma_{SF}<0.77\sigma_{FM}$ being ShapeFit
more constrictive for all the minimum scales. For the $h$ parameter, we have a
central value with a difference around $10^{-3}$ and error bars of
$1.00\sigma_{FM}<\sigma_{SF}<1.11\sigma_{FM}$, with ShapeFit less
constrictive. Lastly, $\mathrm{ln}(10^{10}A_{s})$, shows a central value
difference of the order of $10^{-2}$ with error bars
$1.0\sigma_{FM}<\sigma_{SF}<1.09\sigma_{FM}$, ShapeFit being less
constrictive.

On the other hand, when we include $n_{s}$ as a model extension, we implement
the extra parameter $n$, in the ShapeFit methodology, so we vary
$\alpha_{\perp}$, $\alpha_{\parallel}$, $f\sigma_{8}$, $m$ and $n$, then, we
convert these parameters into cosmological parameter using
$n=n_{s}-n_{\texttt{ref}}$ [94], these results are in Table 13. We notice that
for $\Omega_{m}$, we have a central value of the order of $10^{-3}$, with a
difference around $0.72\sigma_{FM}<\sigma_{SF}<0.93\sigma_{FM}$ being ShapeFit
more constrictive. In the case of $h$, we have a central value difference of
the order $10^{-3}$ and error bars which show a difference of
$0.99\sigma_{FM}<\sigma_{SF}<1.11\sigma_{FM}$. Then, we have the
$\mathrm{ln}(10^{10}A_{s})$ parameter, which shows a central value difference
of the order of $<10^{-2}$ with error bars
$0.88\sigma_{FM}<\sigma_{SF}<1.18\sigma_{FM}$. Finally, the $n_{s}$ parameter
shows a difference in the central value of the order of $10^{-3}$ while the
error bars are $1.32\sigma_{FM}<\sigma_{SF}<1.54\sigma_{FM}$ with ShapeFit
being less constrictive.

So far, we extensively explore two distinct methodologies for predicting the
recovery of cosmological parameters: ShapeFit and FullModeling. We demonstrate
the accuracy and efficacy of both approaches for the baseline analysis.

Figure 15: Comparison of the Full Modeling and ShapeFit methodologies in the
cosmological parameter space, the red bars and the blue bars represent
respectively the ShapeFit and the Full Modeling results obtained using the
mean of the 25 realizations of the ABACUS-SUMMIT simulations with a non-
rescaled covariance matrix when we fitted the monopole and quadrupole.
Cosmological Parameter Constraints
---
Min Freedom
T | $s$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$ |  |
FM | $22$ | $0.3150\pm 0.0072$ | $0.6755\pm 0.0072$ | $3.058\pm 0.065$ | - | -
FM | $30$ | $0.3139\pm 0.0082$ | $0.6767\pm 0.0078$ | $2.989\pm 0.082$ | - | -
FM | $38$ | $0.314\pm 0.010$ | $0.6786\pm 0.0080$ | $2.96\pm 0.11$ | - | -
SF | $22$ | $0.3144\pm 0.0058$ | $0.6743\pm 0.0080$ | $3.082\pm 0.087$ | - | -
SF | $30$ | $0.3146\pm 0.0062$ | $0.6758\pm 0.0082$ | $3.01\pm 0.10$ | - | -
SF | $38$ | $0.3155\pm 0.0077$ | $0.6783\pm 0.0085$ | $22.97\pm 0.12$ | - | -
Hexadecapole
FM | $22$ | $0.3138^{+0.0099}_{-0.0085}$ | $0.6707\pm 0.0072$ | $3.098^{+0.066}_{-0.082}$ | - |
FM | $30$ | $0.3094\pm 0.0082$ | $0.6711\pm 0.0071$ | $3.060\pm 0.095$ | - | -
FM | $38$ | $0.3107\pm 0.0093$ | $0.6741\pm 0.0075$ | $3.01\pm 0.13$ | - | -
SF | $22$ | $0.3142\pm 0.0065$ | $0.6684\pm 0.0072$ | $3.117\pm 0.081$ | - | -
SF | $30$ | $0.3113\pm 0.0060$ | $0.6709\pm 0.0076$ | $3.09\pm 0.11$ | - | -
SF | $38$ | $0.3128\pm 0.0072$ | $0.6730\pm 0.0083$ | $3.07\pm 0.13$ | - | -
$n_{s}$
T | $s$ | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$ | $n_{s}$ |
FM | $22$ | $0.307\pm 0.014$ | $0.6746\pm 0.0076$ | $3.05^{+0.10}_{-0.12}$ | $0.975\pm 0.037$ | -
FM | $30$ | $0.313^{+0.015}_{-0.017}$ | $0.6752\pm 0.0086$ | $2.99\pm 0.15$ | $0.967\pm 0.045$ | -
FM | $38$ | $0.314\pm 0.018$ | $0.6773\pm 0.0089$ | $2.96\pm 0.17$ | $0.967\pm 0.050$ | -
SF | $22$ | $0.326^{+0.033}_{-0.038}$ | $0.681\pm 0.013$ | $2.97^{+0.25}_{-0.29}$ | $0.93\pm 0.14$ | -
SF | $30$ | $0.319^{+0.031}_{-0.037}$ | $0.680^{+0.013}_{-0.015}$ | $2.99\pm 0.27$ | $0.96\pm 0.14$ | -
SF | $38$ | $0.317^{+0.031}_{-0.038}$ | $0.680^{+0.013}_{-0.015}$ | $3.01\pm 0.30$ | $0.97\pm 0.15$ | -
Table 13: Constraints on our parameters when fitting our ABACUS LRG sample for
three minimum scales when using Full Modeling methodology and ShapeFit
methodologies.

##  6 Comparisons with others Effective Field Theory models

As stated at the beginning of this manuscript, our EFT-GSM model is one of
four EFT theories utilized in the analysis of galaxy Full-Shape within DESI
year one data. In the companion paper [52], the three models in Fourier space
(FS) (namely, the Velocileptors model of [51], the Folps$\nu$ model of [48],
and the PyBird model of [46]) are compared to each other, this section extends
the comparison by showing how their results align with our configuration space
(CS) counterpart presented in this work. We highlight that our GSM-EFT
methodology is the only model fully developed in CS out of the four models
considered here.

We divide this comparison into two parts. The first section 6.1 shows two
results where the analysis of the three models we are comparing with is done
in FS. The first result is a direct comparison between our parameter
constraints in CS and the original FS results of Velocileptors, Folps$\nu$,
and PyBird. The second comparison involves using a different implementation of
these three codes, where we calculate their results using the
desilike111111https://github.com/cosmodesi/desilike/ tool. This tool is a
compendium of all the methodologies required for the likelihood analysis of
DESI. It’s designed to provide a common framework for all analyses. Therefore,
it allows us to execute all the important parts of our methodologies using the
same tool for all models. This includes the MCMC sampler, likelihood
estimations, methods for fabricating the input linear power spectrum and so
on. Therefore when using desilike, the only difference between the
methodologies used to analyze each of the four models is the method used to
compute the redshift space multipoles from the input linear power spectrum.

Working with desilike is not only useful because it allows us to standardize
the methodology between codes, but it also provides a standardized tool to
transform FS multipoles into the CS. Therefore, it permits us to reanalyze the
Velocileptors, Folps$\nu$, and PyBird results, this time in CS. Section (6.2)
shows how our EFT-GSM results compare with these new CS model constraints.

Throughout this whole section, we present all analyses in both Full Modeling
and ShapeFit configurations. Table 14 provides a short description of the
codes used in this comparison. For a full description of the codes and
individual performance, we refer the reader to the companion papers.

Model | Description
---|---
Folps-$\nu$ | EPT/EFT with beyond-EdS kernels [48]
PyBird | EPT /EFT [46]
Velocileptors | LPT/EFT with the possibility to run EPT/EFT [51].
Table 14: Brief description of the models compared in FS

The four models utilize distinct conventions for their different galaxy-halo
bias parameters. The translation between these conventions can sometimes be
challenging to reproduce across all models. Therefore, in this study and its
companion paper, we concentrate only on cases where this mapping is well-
defined. While [51] explores the two cases minimal and maximal freedom, all
comparisons presented here are done using only the minimal freedom
configuration. This is done because our maximal freedom GSM-EFT model does not
implement the tidal bias parameter of [52] , and therefore our maximal freedom
definition is not equivalent to one of the other models.

The definitions of minimal freedom for FS models are presented in [52],
equations 4.15 and 4.16. Our comparisons are conducted by fitting the mean of
the Abacus mocks using our covariance matrix without any re-scaling, which
corresponds to having a volume of V=8$h^{-1}\mathrm{Mpc}$.

###  6.1 Comparison with EFT models in Fourier Space

As stated above, throughout this section, we present a set of results from
analyses conducted in FS. All of these FS analyses select the same fiducial
settings as those utilized in the companion paper. Therefore, we fit both the
monopole and quadrupole up to $k_{\rm max}=0.25h^{-1}$Mpc. We begin this
section by introducing a comparison between our CS results and the original FS
results of Velocileptors, Folps$\nu$, and PyBird. These comparisons are done
for all three tracers: LRG, ELG, and QSO. We utilize the priors presented in
Tables 1, 2, and 3 of the companion paper [52].

Figure 16 shows a set of triangular plots with the constraints we get using
our Full Modeling (left panels) and ShapeFit (right panels) configurations,
for our LRG (top panels), ELG (intermediate panels), and QSO (bottom panels)
data respectively. It is important to notice that the posteriors of the FS
results are generated with the original codes and are the same as those
analyzed in the companion paper. We observe that the 1-$\sigma$ contours of
our GSM-EFT model align with those from the FS models showing very good
agreement, while our model gives slightly narrower constraints. The observed
difference between results may not be unexpected, considering that not all
analyses are conducted in the same space and are not done with one
standardized framework.

We observe that our CS constraints consistently yield smaller values compared
to the three FS results obtained from other codes. Additionally, this trend
holds true for all three tracers. As we mentioned, we have not standardized
our framework in terms of MCMC sampler, likelihood evaluation and the fitting
ranges between FS and CS are not exactly equivalent, which could explain the
differences observed between methodologies. With this in mind, we test how our
model compares with the rest when all the FS models are analyzed using the
common framework from desilike. Figure 17 shows the results from this analysis
and compares it with our GSM-EFT results. This comparison is complementary to
the analysis of the next sub-section, where we use the desilike framework to
redo all likelihood analysis in CS instead. These tests are intended as a
verification that the differences observed so far are still present when we
homogenize our analysis with desilike. Tables 15 and 16 present the
quantitative results of this FS desilike analysis. Before analysing desilike
results we verified that the results that the Full Modeling and ShapeFit
constraints and errors were consistent between the original pipeline and their
implementation within desilike. As we highlight in the next section, the
difference found between FS and CS is also found in the CS implementation of
the FS models, and is larger than the difference found between FS results with
different pipelines.

Figure 16: All models using their baseline analysis for all tracers. Need to
replace left plots with the Full Modeling cases, change the format of the
plot, i.e put in solid GSM-EFT

###  6.2 Comparison with EFT models in Configuration space with DesiLike

In this final section we discuss the performance of all four models in CS. It
is important to note that the Velocileptors, Folps$\nu$, and PyBird, models
work in FS by construction. The CS counterparts of the original FS code are
also implemented in desilike, the transformation between one space and the
other is done using the Hankel transform of the FS multipoles,

$\xi_{\ell}(s)=\frac{(-i)^{\ell}}{2\pi^{2}}\int
dkk^{2}P_{\ell}(k)j_{\ell}(ks)$ (6.1)

The Hankel transform is implemented in the
cosmoprimo121212https://github.com/cosmodesi/cosmoprimo/ package using the
fftlog algorithm. For the Hankel transform we used a $k_{min}=10^{-4}$ and
$k_{max}=0.6$. Beyond 0.6, there is a Gaussian damping.

Figure 18 shows the comparison between our GSM-EFT results and the desilike
configurations space versions of Velocileptors, Folps$\nu$, and PyBird. The
left panel shows the Full Modeling results and the right panel the ShapeFit
counterparts for the LRG sample and for different values of $s_{min}$. We
observe good agreement between the models at all scales, we also note that
there is slightly more dispersion than the one we observed in the FS analysis,
which may be related to the sensitivity of the model to the parameters used in
the Hankel transform. In order to getting a closer look, Figure 19 shows the
triangular plot for both methodologies for $s_{min}=30$. We observe that the
contours are all well aligned with each other, and most of them are centered
within less than $1\sigma$ from the true values. Tables 15 and 16 present the
quantitative results of the CS with desilike. We notice that the errors
obtained from the FS models in CS are consistently lower than their
counterparts in FS, confirming the conclusion of the past section that CS
provides slightly smaller error bars. The error bars of GSM-EFT are smaller
than the mean error by 5-10 percent depending on the parameter. We consider
the level of agreement to be reasonable given that the only model fully
developed in CS is GSM-EFT and the others are generated from transforming the
multipoles which may add a layer of potential numeric errors that require to
be explored by tuning the internal parameters of the Hankel transform some of
the parameters, such as the binning and the value of $k_{max}$. A careful
exploration of the dependence of the configuration space constraints on the
numerical precision on the Hankel transform and the exploration of the
differences between precision depending of the different spaces is left of a
future job.

Figure 17: All models using their baseline analysis for the minimal freedom
case with desilike

Figure 18: All models using their baseline analysis for the minimal freedom
case.

Figure 19: All models using their baseline analysis for the minimal freedom
case for the correlation function.  Comparison Models in CS Full Modeling with
desilike .
---
Model | Space | $\Omega_{m}$ | $h$ | $\mathrm{ln}(10^{10}A_{s})$
Velocileptors LPT | FS | $0.309\pm 0.012$ | $0.6739\pm 0.0089$ | $3.030^{+0.083}_{-0.097}$
Folps$\nu$ | FS | $0.309\pm 0.011$ | $0.6740\pm 0.0091$ | $3.051^{+0.083}_{-0.097}$
Pybird | FS | $0.311^{+0.012}_{-0.010}$ | $0.6737\pm 0.0090$ | $3.03\pm 0.11$
Velocileptors LPT | CS | $0.3117\pm 0.0086$ | $0.6757\pm 0.0076$ | $3.008\pm 0.089$
Folps$\nu$ | CS | $0.3138\pm 0.0093$ | $0.6782\pm 0.0077$ | $2.989\pm 0.093$
Pybird | CS | $0.3095\pm 0.0087$ | $0.6762\pm 0.0086$ | $3.01^{+0.10}_{-0.12}$
GSM-EFT | CS | $0.3139\pm 0.0082$ | $0.6767\pm 0.0078$ | $2.989\pm 0.082$
GSM-EFT-22 | CS | $0.3150\pm 0.0072$ | $0.6755\pm 0.0072$ | $3.058\pm 0.065$
Table 15: Constraints on Full modeling parameters when fitting the ABACUS-LRG
for different models and different spaces, minimal freedom, using desilike
compared with this work. Comparison Models in CS ShapeFit with desilike .
---
Model | Space | $f\sigma_{8}$ | $\alpha_{\parallel}$ | $\alpha_{\perp}$ | m
Velocileptors LPT | FS | $0.448\pm 0.024$ | $1.004\pm 0.021$ | $0.9987\pm 0.0099$ | $-0.023^{+0.039}_{-0.032}$
Folps$\nu$ | FS | $0.453\pm 0.024$ | $1.004\pm 0.021$ | $0.9989\pm 0.0099$ | $-0.024^{+0.037}_{-0.031}$
Pybird | FS | $0.451\pm 0.026$ | $1.001\pm 0.021$ | $0.9995\pm 0.010$ | $-0.024^{+0.043}_{-0.031}$
Velocileptors LPT | CS | $0.444\pm 0.022$ | $1.002\pm 0.017$ | $0.9972\pm 0.0078$ | $-0.009^{+0.031}_{-0.028}$
Folps$\nu$ | CS | $0.447\pm 0.023$ | $1.002\pm 0.017$ | $0.9982\pm 0.0076$ | $-0.012^{+0.034}_{-0.031}$
Pybird | CS | $0.451\pm 0.024$ | $1.012^{+0.020}_{-0.018}$ | $0.9957\pm 0.0074$ | $-0.046\pm 0.032$
GSM-EFT-30 | CS | $0.442\pm 0.021$ | $1.005\pm 0.017$ | $0.9963\pm 0.0077$ | $-0.001\pm 0.027$
GSM-EFT-22 | CS | $0.453\pm 0.020$ | $1.009\pm 0.017$ | $0.9966\pm 0.0074$ | $-0.004^{+0.027}_{-0.021}$
Table 16: Constraints on ShapeFit parameters when fitting the ABACUS-LRG for
different models and different spaces, minimal freedom using desilike compared
with this work.

##  7 Conclusions

In preparation for the analysis of the upcoming first batch of DESI data, it
is imperative to analyze and test the tools that will be used to extract
cosmological information. This includes testing different configurations and
degrees of freedom that models can have when predicting the full-shape of
galaxy two-point statistics. In particular, there is a lot of interest in
understanding the differences in accuracy and precision between compressed
methodologies like the Standard and Shape-Fit methods, which are more model
agnostic, and direct methods like the Full Modeling methodology where a
particular cosmological model is selected a priori.

Currently, four independent methodologies that model the redshift space
clustering statistics are being developed and tested within the collaboration
in preparation for the first data release. These are the PyBird, FOLPS$\nu$,
and velocileptors methodologies, all of which study the FS power spectrum of
tracers. The fourth method is our own EFT-GSM methodology, which is the only
approach working on CS.

Throughout this work, we have analyzed the performance of our GSM-EFT model,
along with some extensions that allow us to explore a wider parameter range,
in recovering the cosmological information from a new suite of 25 ABACUS-
SUMMIT mocks. These mocks have been independently populated with different
tracers of matter to generate samples resembling ELGs, LRGs, QSOs. Here, we
have analyzed all three samples independently, and we also examine a combined
sample. We have compared our constraints with those produced by the three
other FS models when transforming their predictions into CS. The analysis
presented here utilizes our optimized C-code tool introduced in this paper,
which computes the CS multipoles of our GSM-EFT model. Furthermore, we have
optimized our methodology by employing Neural Networks to construct surrogate
models of our GSM-EFT models, accelerating their evaluation and thus reducing
the convergence time required to obtain cosmological constraints.

We have analyzed the ability of our Full Modeling configuration to recover the
ABACUS -SUMMIT LRG cosmology under variations in the volume of the simulation.
This was achieved by running our analysis on volumes of
$8,h^{-3},\textrm{Gpc}^{3}$, $40,h^{-3},\textrm{Gpc}^{3}$, and
$200,h^{-3},\textrm{Gpc}^{3}$, respectively, which correspond to approximately
the first and fifth years of DESI observations, and a much larger volume than
DESI, respectively. We find that our methodology is relatively robust to
variations in volume, with all parameters being recovered to within $1\sigma$
of their original values, with the exception of the values of $h$ and
$\mathrm{ln}(10^{10}A_{s})$ for the larger volume, which are recovered to
within $1.14\sigma$ and $1.4\sigma$, respectively. Throughout this work, we
choose to report our results based on an $8,h^{-3},\textrm{Gpc}^{3}$ box, as
this approximates the size of the DESI year one volume.

We have also tested how our results vary in our baseline analysis of the LRG
sample under variations of the minimum scale utilized in our analysis,
employing both Full Shape and Shape fit methodologies. We explore three
possible values of the minimum scale, namely
$s_{\text{min}}=22,h^{-1}\mathrm{Mpc}$,
$s_{\text{min}}=30,h^{-1}\mathrm{Mpc}$, and
$s_{\text{min}}=38,h^{-1}\mathrm{Mpc}$. We found that both methodologies are
in good agreement with each other and with the true values of the simulation
for all minimum scales explored and for all parameters within $1\sigma$. Given
that the smaller scale corresponds to tighter constraints, we select
$s_{\text{min}}=22,h^{-1}\mathrm{Mpc}$ as the baseline value utilized
throughout the rest of this work. We have also tested how our methodology is
affected by switching between the minimal and maximal freedom scenarios, where
$b_{s^{2}}$ is either set free or fixed to zero, and we define maximal freedom
as our baseline configuration. We found that while all predictions are
consistent within $1\sigma$ for $V_{1}$, the minimal freedom case yields
constraints approximately $10\%$ tighter. This improvement in constraining
power is expected when one reduces the parameter space explored.

In this work, we have also analyzed several extensions to our standard
configuration. The first of these extensions is the inclusion of the
hexadecapole of the correlation function into our analysis, which is
considered as the multipole where first-loop PT starts to find difficulties
modeling our scales of interest. In this extension, we include a model of the
hexadecapole but only down to a minimum scale of
$s_{\text{min}}=30\,h^{-1}\mathrm{Mpc}$ as we do not trust our model at
smaller scales. We note that this extension does not significantly improve our
results, with its only advantage being a small improvement in the precision of
the Full Modeling estimates. We also highlight that all predictions recover
the true value of all parameters within $1\sigma$.

We have explored the inclusion of the spectral index $n_{s}$ into our
methodology as another extension to our model and analyzed how this affects
the fits on the LRG sample using the Full Shape configuration of our model. We
found that our model predicts all parameters, including $n_{s}$, to within
$1\sigma$. Additionally, we observe that the ability to recover the original
parameters improves over their baseline counterparts at all scales. It is
noteworthy that under this extension, our minimum scale of
$s_{\text{min}}=30\,h^{-1}\mathrm{Mpc}$ shows the best performance overall.

As another extension to our baseline analysis we have expanded our standard
setup to include the $w_{0}$ and $w_{a}$ parameters of the
Chevallier–Polarski–Linder model of dark energy. We run our fitting
methodology on our LRG sample and for two distinct configurations of these
extensions, first we allow for only $w_{0}$ to vary while keeping $w_{a}$
fixed at its true values. Then we allow both parameters to vary
simultaneously. We find that all constrains using our Full-Modeling
configuration are capable of recovering the ABACUS-SUMMIT parameters,
including $w_{0}$ and $w_{a}$ well bellow the $1\sigma$ range, we also note
that in both cases the $s_{\text{min}}=22\,h^{-1}\mathrm{Mpc}$ shows the
overall best performance. We believe the accuracy of our fits is in part a
result of working in CS, where volume effects play a smaller role on these new
parameters.

All the tests and extensions discussed so far have been analyzed using our LRG
sample. As a separate test, we have analyzed the accuracy of our methodology
when fitting the QSO and ELG samples from ABACUS-SUMMIT. We have conducted
fits for these two new samples using our standard configuration for three
different minimum scales. In general, we observe that the ELG fits are
slightly more accurate and precise than the LRG fits, while the QSO fits are
generally less precise. This discrepancy is expected due to the differences in
the number of tracers between the samples. Nevertheless, our methodology
successfully recovers all parameters within $1\sigma$ for all three samples.
Additionally, we find that using a minimum scale of
$s_{\text{min}}=22,h^{-1}\mathrm{Mpc}$ provides tighter constraints on our
parameter space while also predicting all parameters precisely.

We have also conducted tests comparing the performance of Full Modeling and
ShapeFit within cosmological parameter space. They involved examining the
impact of variations in certain settings, including bias settings and minimum
scale ranges. Additionally, we explored how these configurations were affected
by our three model extensions: hexadecapole and $n_{s}$ free. We observed that
the various configurations of our vanilla model generally exhibit agreement in
their estimates of the central values between both methodologies. There are
slight variations in the size of their error estimates, ranging from a few to
30 percent depending on the case studied and the specific parameter. On the
other hand, when considering extensions beyond the baseline that involve
varying additional cosmological parameters, we observed larger disparities in
both the central values and the size of the error bars. Specifically, when
allowing $n_{s}$ to vary freely, we found that the Full Modeling approach
imposes greater constraints, particularly on $n_{s}$ compared to ShapeFit.

The last tests we have presented analyze the performance of our CS
methodology, GSM-EFT, when compared to three FS codes: Velocileptors,
Folps$\nu$, and PyBird. This analysis serves as a complement to the main
comparison between RSD methodologies presented in a companion paper dedicated
to the consistency of all FS models. We found that our GSM-EFT results are in
good agreement with those from all FS models, with our results showing smaller
error bars for both Full Modeling and ShapeFit configurations. We have also
investigated the performance of models in CS through the common framework
desilike, confirming the aforementioned consistency of models in both spaces
and the reduction of error bars in CS.

##  8 Data Availability

The data will be made publicly available once the paper is accepted.

## Acknowledgments

SR, HN, SF, FR and MV are supported by Investigacion in Ciencia Basica
CONAHCYT grant No. A1-S-13051 and PAPIIT IN108321 and IN116024, and Proyecto
PIFF. SF and FR is supported by PAPIIT IN115424. HN and AA are supported by
Ciencia de Frontera grant No. 319359 and PAPIIT IG102123.

MIL recognizes the support of the National Research Foundation of Korea (NRF)
grant funded by the Korea government (MIST), Grant No. 2021R1A2C1013024

This material is based upon work supported by the U.S. Department of Energy
(DOE), Office of Science, Office of High-Energy Physics, under Contract No.
DE–AC02–05CH11231, and by the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility under the same contract.
Additional support for DESI was provided by the U.S. National Science
Foundation (NSF), Division of Astronomical Sciences under Contract No.
AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research
Laboratory; the Science and Technology Facilities Council of the United
Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation;
the French Alternative Energies and Atomic Energy Commission (CEA); the
National Council of Humanities, Science and Technology of Mexico (CONAHCYT);
the Ministry of Science and Innovation of Spain (MICINN), and by the DESI
Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the U. S. National Science Foundation, the U. S. Department of Energy, or
any of the listed funding agencies.

The authors are honored to be permitted to conduct scientific research on
Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the
Tohono O’odham Nation.

##  Appendix A Author Affiliations

.5cm1

1Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México
C.P. 04510, México

2Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu,
Daejeon 34055, Republic of Korea

3Instituto de Ciencias Físicas, Universidad Autónoma de México, Cuernavaca,
Morelos, 62210, (México)

4Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720,
USA

5Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215,
USA

6University of Michigan, Ann Arbor, MI 48109, USA

7Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK

8Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK

9Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK

10NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA

11Department of Physics & Astronomy and Pittsburgh Particle Physics,
Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941
O’Hara Street, Pittsburgh, PA 15260, USA

12Kavli Institute for Particle Astrophysics and Cosmology, Stanford
University, Menlo Park, CA 94305, USA

13SLAC National Accelerator Laboratory, Menlo Park, CA 94305, USA

14Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10,
Edificio Ip, CP 111711, Bogotá, Colombia

15Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10,
Edificio H, CP 111711 Bogotá, Colombia

16Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

17Institute of Cosmology and Gravitation, University of Portsmouth, Dennis
Sciama Building, Portsmouth, PO1 3FX, UK

18Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans
s/n, 08913 Bellaterra, Barcelona, Spain

19Departament de Física Quàntica i Astrofísica, Universitat de Barcelona,
Martí i Franquès 1, E08028 Barcelona, Spain

20Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), c.
Martí i Franquès, 1, 08028 Barcelona, Spain.

21Center for Cosmology and AstroParticle Physics, The Ohio State University,
191 West Woodruff Avenue, Columbus, OH 43210, USA

22Department of Physics, The Ohio State University, 191 West Woodruff Avenue,
Columbus, OH 43210, USA

23The Ohio State University, Columbus, 43210 OH, USA

24School of Mathematics and Physics, University of Queensland, 4072, Australia

25Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona,
08193 Bellaterra (Barcelona), Spain

26Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain

27University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA
94720, USA

28Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís
Companys, 23, 08010 Barcelona, Spain

29Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH,
U.K

30Department of Physics & Astronomy, University of Wyoming, 1000 E.
University, Dept. 3905, Laramie, WY 82071, USA

31National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun
Rd., Chaoyang District, Beijing, 100012, P.R. China

32Department of Physics and Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada

33Perimeter Institute for Theoretical Physics, 31 Caroline St. North,
Waterloo, ON N2L 2Y5, Canada

34Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave
W, Waterloo, ON N2L 3G1, Canada

35Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way,
Berkeley, CA 94720, USA

36Department of Physics, Kansas State University, 116 Cardwell Hall,
Manhattan, KS 66506, USA

37Department of Physics and Astronomy, Sejong University, Seoul, 143-747,
Korea

38CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

39Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

40Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA

41Physics Department, Stanford University, Stanford, CA 93405, USA