# An analysis of parameter compression and Full-Modeling techniques with
Velocileptors for DESI 2024 and beyond

M. Maus  0000-0002-9020-911X S. Chen  0000-0002-5762-6405 M. White
0000-0001-9912-5070 J. Aguilar  S. Ahlen  0000-0001-6098-7247 A. Aviles
0000-0001-5998-3986 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 P. Doel  S. Ferraro  0000-0003-4992-7854 N. Findlay
0009-0007-0716-3477 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 C.
Hahn  0000-0003-1197-0902 K. Honscheid  C. Howlett  0000-0002-1081-9410 M.
Ishak  0000-0002-6024-466X S. Juneau  A. Kremin  0000-0001-6356-7424 Y. Lai
M. Landriau  0000-0003-1838-8528 M. E. Levi  0000-0003-1887-1018 M. Manera
0000-0003-4962-8934 R. Miquel  E. Mueller  A. D. Myers  S. Nadathur
0000-0001-9070-3102 J. Nie  0000-0001-6590-8122 H. E. Noriega
0000-0002-3397-3998 N. Palanque-Delabrouille  0000-0003-3188-784X W. J.
Percival  0000-0002-0644-5727 C. Poppett  S. Ramirez-Solano  M. Rezaie
0000-0001-5589-7116 A. Rocher  0000-0003-4349-6424 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 M. Vargas-Magaña
0000-0003-3841-1836 B. A. Weaver  S. Yuan  0000-0002-5992-7586 P. Zarrouk
0000-0002-7305-9578 H. Zhang  0000-0001-6847-5254 R. Zhou  0000-0001-5381-4372
H. Zou  0000-0002-6684-3997

###### Abstract

In anticipation of forthcoming data releases of current and future
spectroscopic surveys, we present the validation tests and analysis of
systematic effects within velocileptors modeling pipeline when fitting mock
data from the AbacusSummit N-body simulations. We compare the constraints
obtained from parameter compression methods to the direct fitting (Full-
Modeling) approaches of modeling the galaxy power spectra, and show that the
ShapeFit extension to the traditional template method is consistent with the
Full-Modeling method within the standard $\Lambda$CDM parameter space. We show
the dependence on scale cuts when fitting the different redshift bins using
the ShapeFit and Full-Modeling methods. We test the ability to jointly fit
data from multiple redshift bins as well as joint analysis of the pre-
reconstruction power spectrum with the post-reconstruction BAO correlation
function signal. We further demonstrate the behavior of the model when opening
up the parameter space beyond $\Lambda$CDM and also when combining likelihoods
with external datasets, namely the Planck CMB priors. Finally, we describe
different parametrization options for the galaxy bias, counterterm, and
stochastic parameters, and employ the halo model in order to physically
motivate suitable priors that are necessary to ensure the stability of the
perturbation theory.

##  1 Introduction

The large-scale structure (LSS) of the Universe is the observed, coherent
spatial distribution of material on scales larger than the typical galaxy or
halo scale, and provides a powerful observational tool for probing cosmic
evolution. LSS observations allow us to study 3D volumes of the sky that span
a long range of cosmic times, enabling us to study the initial conditions of
the primordial universe as well as its evolution at later times. [1, 2, 3, 4].

One of the primary methods of measuring the evolution of LSS is through galaxy
redshift surveys that aim to probe the clustering of matter on a wide range of
scales using galaxies as tracers. Spectroscopic galaxy surveys have had
significant success over the years in scanning large regions of the sky. These
include the 2dF [5], 6dF [6], GAMA [7], WiggleZ [8], and most recently the
completed Sloan Digital Sky Survey (SDSS), composed of data from SDSS, SDSS-II
[9], BOSS [10, 11, 12], and eBOSS [13, 14, 15]. The next telescope surveys to
further push the boundaries of LSS observations that have recently begun
operations are the Euclid Satellite [16, 17] and the ground-based Dark Energy
Spectroscopic Instrument (DESI) [18, 19, 20]. DESI aims to cover over 14,000
deg2 by the end of 5 years of observations, with target samples of stars from
the Milky Way Survey (MWS), bright galaxies from the Bright Galaxy Survey
(BGS, $0.0<z<0.4$), Luminous Red Galaxies (LRG, $0.4<z<1.1$), Emission Line
Galaxies (ELG, $1.1<z<1.6$), and Quasars (QSO, $1.6<z<2.1$). Altogether the
DESI survey will span an effective volume of about
$20\,(h^{-1}\mathrm{Gpc})^{3}$ by the end of its 5 years of observation [21].

In anticipation of the upcoming Year-1 data release of DESI [22, 23, 24, 25,
26, 27, 28, 29, 30] (as well as later releases along with Euclid), it is
important to characterize the performance of the current state-of-the-art
models for analyzing the observed galaxy clustering 2-point statistics and the
resultant cosmological constraints. The growth of large-scale structure is a
competition between gravity, the dominant force on large scales, and the
expansion of the universe. Models must also include several other effects:
First, galaxies are not perfect tracers of the underlying matter overdensity
field, and thus a ‘biasing’ scheme is needed in order to relate the matter
power spectrum to the observed galaxy spectrum (see ref. [31] for a recent
review). Second, since distances along the line-of-sight (LOS) are inferred
from redshifts, components of galaxy peculiar velocities in the LOS direction
influence the inferred distances and are a source of anisotropy in the
observed clustering signal [32, 33]. This latter effect is known as redshift
space distortions (RSD) and provides both a challenge to modeling while also
giving direct access to information about the growth rate of LSS. Finally,
nonlinear effects on small scales must be included. We use perturbation theory
to model the mildly non-linear regime, with additional parameters to account
for the small-scale physics such that the models are not sensitive to the
complicated processes e.g. involved with galaxy formation (sometimes known as
Effective Field Theory or EFT terms [34, 35, 36]). The model considered in
this work, velocileptors111https://github.com/sfschen/velocileptors/tree/2.0
[37, 38], is one of the models that will will be used for analyzing the full-
shape power spectra from the upcoming DESI survey data releases, the others
being the Fourier space Eulerian PT codes PyBird [39, 40, 41] and FOLPS$\nu$
[42] and the configuration space code EFT-GSM [43]. The purpose of this work
is to characterize the performance of velocileptors and understand any
systematic issues by comparing to a suite of simulated, or ‘mock’, data.
Similar tests are being performed with the other three models in addition to a
comparison between models, and will be reported in companion publications[44,
45, 46, 47]. While velocileptors has been tested previously on simulations
[38, 48, 49], here we focus on DESI-like galaxies and redshift ranges, and
also use the new AbacusSummit [50] suite of simulations produced for the DESI
collaboration that is also used to test the other theory models.

Within the framework of the model, there are still various approaches to
fitting data. One method, previously used by the BOSS and eBOSS
collaborations, involved choosing a fiducial template for the linear power
spectrum while compressing the observed power spectrum multipoles into three
parameters: the amplitude of the redshift-space anisotropy $f\sigma_{8}$, and
the two scaling parameters parallel and perpendicular to the line of sight,
i.e. $\alpha_{\parallel}$ and $\alpha_{\perp}$. This technique was meant to
encode the intuition that, for currently popular cosmological models, primary
CMB anisotropies fix the parameters determining the shape of the power
spectrum but late-time effects such as non-trivial dark energy evolution or
spatial curvature can affect the total growth and the distance-redshift
relation. These impacts are accounted for by the three parameters above and
redshift surveys can constrain them well. An extension to this standard
“template” fit is to include another compressed “ShapeFit” parameter to allow
a set of modifications to the shape of the linear power spectrum [51]. The
extra shape information of this method allows for tighter constraints on
cosmological parameters when interpreting the compressed statistics in light
of a given cosmological model without including CMB priors. This partially
bridges the gap in constraining power between the traditional template fit and
the direct fitting or “Full-Modeling” approach of directly varying the
parameters of a specific cosmological model. In this paper we compare these
three methods under a variety of conditions in order to better understand the
advantages and disadvantages of the methods. A comparison of the template and
Full-Modeling approaches was investigated in ref. [52] on the BOSS DR12
dataset, specifically focusing on shifts in $f\sigma_{8}$ constraints between
the two methods. Here we extend that analysis to include the ShapeFit method
and compare the three methods for the range of different settings,
parameterizations, and modeling choices.

This paper is organized as follows. We begin by describing the Abacus
simulations in Sect. 2 and give an overview of Lagrangian Perturbation Theory
(LPT) and velocileptors in Sect. 3. We describe the parameter compression and
Full-Modeling fitting methods in more detail in Sect. 4. The results of our
primary tests, namely the dependence on scale cuts, joint fitting of multiple
redshift bins, post-reconstruction statistics, $w$CDM models, CMB priors,
varying $n_{s}$, Lagrangian vs Eulerian (EPT) Perturbation Theory, and freeing
$\sigma_{8}$ are presented in Sect. 5. We conclude the paper in Sect. 6. We
also provide a brief discussion of our method for analytic marginalization
over the linear parameters in our model in Appendix A along with some further
tests, namely the dependence of $\omega_{b}$ prior, inclusion of cubic bias,
and inclusion of hexadecapole moment in Appendix D. In Appendix B we discuss
the issue of parameter projection effects and the dependence on priors within
our model, a problem that also arises in many other areas of cosmology. We
follow this up with a section dedicated to the halo model in Appendix C, which
allows us to estimate typical scales for stochastic parameters in our model
and provide physical motivation for our prior choices. Appendix E explains our
use of emulators based on Taylor series in order to speed up likelihood
evaluations, and we show that they perform consistently with the direct theory
predictions.

##  2 Mock data

Figure 1: Power spectrum monopole (left) and quadrupole (right) mock data for
the LRG, ELG, and QSO tracers. For each tracer, the mean of the 25 N-body
realizations is used. The error bars of the data correspond to the covariance
re-scaled by the number of realizations, which represents a survey volume of
200 $h^{-3}$Gpc3. The shaded regions show the error bars for a single cubic
box, of volume $V=8\,h^{-3}$Gpc3.

To test our theory model we make use of the AbacusSummit [50] suite of N-body
simulations in their native, cubic geometry. These simulations were run with
the Abacus [53] N-body code on the Summit supercomputer at the Oak Ridge
Leadership Computing Facility for use by the DESI collaboration. The
simulations relevant to this work use a fixed cosmology222The Abacus fiducial
cosmology has $h=0.6736$, $\omega_{b}=0.02237$, $\omega_{\rm cdm}=0.12$,
$A_{s}=2.0830\times 10^{-9}$, and $n_{s}=0.9649$, with a corresponding BAO
drag scale of $r_{d}=99.08\,h^{-1}$Mpc, with 25 boxes each with a different
random number seed for the initial conditions run in a
$(2\,h^{-1}\mathrm{Gpc})^{3}$ volume for a combined volume of 200$h^{-3}$Gpc3.
The mock galaxy catalogs have been produced for three types of tracers, each
produced at a different redshift: Luminous Red Galaxies (LRGs) at $z=0.8$,
Emission Line Galaxies (ELGs) at $z=1.1$, and Quasars (QSOs) at $z=1.4$. For
this study we ignore light-cone and evolution effects in order to better study
the non-linear dynamics and biasing models. The RSD power spectrum data for
each tracer is shown in Fig. 1.

The covariance we use for each tracer is calculated by Monte-Carlo from 1000
“effective Zeldovich approximation” (EZmock [54]) simulations of the same
cosmology333Since these computationally efficient simulalions make use of the
Zel’dovich approximation they may not be as accurate at small scales. As we
will show later, our models are able to obtain unbiased constraints up to
$k_{\rm max}=0.2\,h{\rm Mpc}^{-1}$ but analytic covariances may be desirable
in the future.. We compute this covariance numerically via:

$\displaystyle Cov[P(k)]_{ij}=\frac{1}{N-1}\sum_{n}^{N}[P_{n}(k_{i})-\langle
P(k_{i})\rangle][P_{n}(k_{j})-\langle P(k_{j})\rangle]^{\top}$ (2.1)

In principle, when using as data the mean of 25 cubic boxes the error bars of
the data should also be re-scaled to reflect the increase in volume because
$\sigma^{2}\propto V^{-1}$. A proper treatment of the mean of 25 realizations
would therefore involve re-scaling the covariance from the EZmocks by a factor
of 1/25. However, we must be careful in interpreting results when the error
bars of the data are so tight, as the “survey volume” of the simulations is
orders of magnitude larger than any realistic survey will ever be able to
achieve. For example, if we consider a future survey covering 18 000 deg2 with
tracers in a single redshift bin spanning $0.75<z<1.25$, then the comoving
volume of that data would be about 24 ($h^{-1}$Gpc)3, which is still much less
than the 200 ($h^{-1}$Gpc)3 volume of the simulations. The
$8(h^{-1}\mathrm{Gpc})^{3}$ volume of a single box in our simulations is much
closer to what we expect for any tracers/redshift bin by the end of five years
of DESI observations.

The motivation for the large simulation volume is to detect systematic errors
in the models relevant to the DESI Y5 data. If we define the detection of a
systematic error as being larger than twice the statistical error $\sigma_{\rm
sim}$ of the simulations and would like to keep systematic errors below some
fraction $1/n$ of the Y5 data errors ($\sigma_{Y5}$), then this implies that
we desire simulations with $\sigma_{\rm sim}\leq(2n)^{-1}\sigma_{\rm Y5}$. If
$\sigma\propto 1/\sqrt{V}$, then for $n=3$ and a DESI Y5 volume of
$~{}5\,(h^{-1}\mathrm{Gpc})^{3}$, we would require a simulation volume of
$~{}180\,(h^{-1}\mathrm{Gpc})^{3}$. The Abacus simulations fulfill this
requirement. However the above argument fails to account for the systematic
errors of the N-body simulations themselves. The fractional errors of the
Abacus mock LRG monopole data with 25 box covariance (re-scaled by $1/25$) are
roughly $0.15\%$ between $0.15<k<0.2\,\,h{\rm Mpc}^{-1}$. Ref. [55] compared
different cosmological N-body codes and found that RSD power spectra
multipoles differed by $\approx 0.5\%$ in the same $k-$range, i.e. the
simulations themselves do not agree to these levels of precision, even before
uncertainties from initial condition generation, halo finding and additional
physics are included [56]. In addition to this, the large volume also reflects
a level of precision that our models are not designed for, meaning that
contributions from, e.g., two-loop terms that we don’t include in our theory
can result in poor fits. For all of these reasons, we will primarily focus on
results using the un-rescaled covariance of the more reasonable single-box
volume in the analysis of this paper, while only commenting briefly on the 25
box covariance results when relevant. Finally, when computing the covariance
from a finite number of simulations, one should in principle include
corrections such as the Hartlap factor[57], which depends on the number of
bins in the data vector versus the number of independent mock data sets used.
Given the large number of EZmock simulations that we use, this factor is close
to 1 and we therefore do not observe any noticeable change in constraints when
including the correction. We also do not observe any significant bias in
constraints arising from the finite number of mocks and therefore neglect the
Hartlap correction in our analyses.

##  3 Theory and Model

The velocileptors code is based on the Lagrangian Perturbation Theory (LPT)
approach to large-scale structure. This approach treats dark matter as
collisionless particles whose mapping from initial (Lagrangian) positions,
$\boldsymbol{q}$, to their final observed coordinates, $\boldsymbol{x}$ is
given by $\boldsymbol{x}=\boldsymbol{q}+\boldsymbol{\Psi}(\boldsymbol{q})$,
where $\boldsymbol{\Psi}(\boldsymbol{q})$ is the displacement field. The
dynamical equation, based on Newtonian gravity in an expanding spacetime,
$\ddot{\Psi}+\mathcal{H}\dot{\Psi}=-\nabla_{\boldsymbol{x}}\Phi$, is
perturbatively expanded and solved as
$\boldsymbol{\Psi}=\boldsymbol{\Psi}^{(1)}+\boldsymbol{\Psi}^{(2)}+\boldsymbol{\Psi}^{(3)}+...$.
The observed galaxy overdensity is derived from number conservation, with the
inclusion of a bias functional in the initial conditions,
$F[\delta_{0}(\boldsymbol{q})]$, that relates the tracer overdensity field to
the linear matter field in the form of a Taylor series [37, 38]. In Fourier
space, this results in

$\displaystyle 1+\delta_{g}(\boldsymbol{k})=\int d^{3}\boldsymbol{q}\
F[\delta_{0}(\boldsymbol{q})]e^{-i\boldsymbol{k}\cdot(\boldsymbol{q}+\boldsymbol{\Psi}(\boldsymbol{q}))}$
$\displaystyle
F[\delta_{0}(\boldsymbol{q})]=1+b_{1}\delta_{0}+\frac{1}{2}b_{2}(\delta_{0}(\boldsymbol{q})^{2}-\left\langle\delta_{0}^{2}\right\rangle)+b_{s}(s_{0}^{2}(\boldsymbol{q})-\left\langle
s_{0}^{2}\right\rangle)+b_{3}\mathcal{O}_{3}(\boldsymbol{q}),$ (3.1)

where $s_{0}=(\partial_{i}\partial_{j}/\partial^{2}-\delta_{ij}/3)\delta_{0}$
is the initial shear tensor. The Lagrangian biases $b_{O}$ describe the
response of galaxy formation to large-scale perturbations and are the free
parameters of the theory—absent a complete model of galaxy formation at small
scales their values must be measured directly from large-scale observables
like the power spectrum, though rough estimates for their sizes can be made
through toy models like halo occupation distributions. At 1-loop order there
is only one non-degenerate cubic bias contribution which we include
schematically as $\mathcal{O}_{3}$. Note that the Lagrangian bias parameters
here are not equivalent to the Eulerian ones (for example the standard linear
bias is $b=1+b_{1}$) but equivalent under a set of linear transformations (see
e.g. ref. [38]). Throughout most of this paper we will set $b_{3}=0$ under the
assumption that the cubic nonlinearities in galaxy clustering are consistent
with those from dynamical contributions alone [58]. We test this assumption in
Appendix D.

The modeling of observed galaxy clustering statistics is complicated by the
peculiar velocities of the galaxies, whose line-of-sight components introduce
anisotropies in the clustering signal, an effect known as Redshift Space
Distortions (RSD). In LPT, the transformation into redshift space amounts to a
boost along the LOS direction, $\hat{n}$ so that the redshift space
displacement field is

$\displaystyle\boldsymbol{\Psi}_{s}=\boldsymbol{\Psi}+\dot{\boldsymbol{\Psi}}=\boldsymbol{\Psi}+\frac{\hat{n}(\textbf{v}\cdot\hat{n})}{\mathcal{H}},$
(3.2)

where v is the galaxy peculiar velocity and $\mathcal{H}$ is the conformal
Hubble parameter. We can simplify this relation with the Einstein-deSitter
Approximation (EdS), such that

$\displaystyle\boldsymbol{\Psi}_{s}^{(n)}=\boldsymbol{\Psi}^{(n)}+nf(\hat{n}\cdot\boldsymbol{\Psi}^{(n)}),$
(3.3)

where $f$ is the linear growth rate. This can be expressed as a rotation of
the real space field via the matrix
$R^{(n)}=\delta_{ij}+nf\hat{n_{i}}\hat{n_{j}}$ such that
$\boldsymbol{\Psi}^{s,(n)}=R^{(n)}\boldsymbol{\Psi}^{(n)}$. Defining the
pairwise displacement field in redshift space as
$\Delta_{s}=\Psi_{s}(\boldsymbol{q}_{1})-\Psi_{s}(\boldsymbol{q}_{2})$, the
redshift-space galaxy power spectrum can be obtained from the cumulant
expansion of

$\displaystyle P_{s,g}(\boldsymbol{k})=\int d^{3}\boldsymbol{q}\left\langle
e^{i\boldsymbol{k}\cdot(\boldsymbol{q}+\Delta_{s})}F(\boldsymbol{q}_{1})F(\boldsymbol{q}_{2})\right\rangle_{\boldsymbol{q}=\boldsymbol{q}_{1}-\boldsymbol{q}_{2}}.$
(3.4)

In order to accurately capture the effects of long-wavelength (IR) linear
displacements on the power spectrum, particularly with respect to their
smearing of the BAO, it is necessary to include their effects beyond 1-loop
order in perturbation theory [59, 60, 61, 62]. This class of techniques is
known in the literature as “IR resummation”: in our scheme the linear piece,
i.e. the $A_{ij}^{s,(11)}$ component of
$A_{ij}^{s}=\left<\Delta^{s}_{i}\Delta^{s}_{j}\right>$, is split into long-
and short- wavelength components, $A_{ij}^{s,\rm
lin}=A_{ij}^{s,<}+A_{ij}^{s,>}$, with a cutoff scale $k_{\rm IR}$, and we keep
the $A_{ij}^{s,<}$ piece exponentiated while expanding all other contributions
to 1-loop order. Due to the matrix transformation between the real and
redshift space displacements,
$\boldsymbol{\Psi}^{s,(n)}=R^{(n)}\boldsymbol{\Psi}^{(n)}$, both velocities
and displacements contribute to the resummed $A_{ij}^{s}$. The expression for
the power spectrum becomes [38]

$\displaystyle P_{s,g}^{PT}(\boldsymbol{k})$ $\displaystyle=\int
d^{3}\boldsymbol{q}\
e^{i\boldsymbol{k}\cdot\boldsymbol{q}}e^{-\frac{1}{2}k_{i}k_{j}A_{ij}^{s,<}}\left\\{1-\frac{1}{2}k_{i}k_{j}A_{ij}^{s,>}+\frac{1}{8}k_{i}k_{j}k_{k}k_{l}A_{ij}^{s,>}A_{kl}^{s,>}\right.$
$\displaystyle-\frac{1}{2}k_{i}k_{j}A_{ij}^{s,\text{ loop
}}+\frac{i}{6}k_{i}k_{j}k_{k}W_{ijk}^{s}+2ib_{1}k_{i}\left(1-\frac{1}{2}k_{i}k_{j}A_{ij}^{s,>}\right)U_{i}^{s}-b_{1}k_{i}k_{j}A_{ij}^{s,10}$
$\displaystyle+b_{1}^{2}\left(1-\frac{1}{2}k_{i}k_{j}A_{ij}^{s,>}\right)\xi_{\text{lin
}}+ib_{1}^{2}k_{i}U_{i}^{s,11}-b_{1}^{2}k_{i}k_{j}U_{i}^{s,\operatorname{lin}}U_{j}^{s,\operatorname{lin}}$
$\displaystyle+\frac{1}{2}b_{2}^{2}\xi_{\text{lin
}}^{2}+2ib_{1}b_{2}\xi_{\operatorname{lin}}k_{i}U_{i}^{s,\operatorname{lin}}-b_{2}k_{i}k_{j}U_{i}^{s,\operatorname{lin}}U_{j}^{s,\operatorname{lin}}+ib_{2}k_{i}U_{i}^{s,20}$
$\displaystyle\left.+b_{s}\left(-k_{i}k_{j}\Upsilon_{ij}^{s}+2ik_{i}V_{i}^{s,10}\right)+2ik_{i}b_{1}b_{s}V_{i}^{s,12}+b_{2}b_{s}\chi+b_{s}^{2}\zeta+2ib_{3}k_{i}U_{b_{3},i}^{s}+2b_{1}b_{3}\theta+\ldots\right\\}.$
(3.5)

The other correlators appearing above ($\xi$, $W$, $V$, $U$, etc.) are defined
in [59, 63, 37, 38].

We account for the sensitivity to small scales by introducing counterterms
with coefficients, $\alpha_{n}$, that multiply the tree-level power spectrum.
These coefficients describe couplings with short-wavelength modes whose sizes
are not directly specified by perturbation theory. While their exact values
(or even signs) are not known, we can put reasonable priors on them based on
the size of gravitational nonlinearities seen in N-body simulations and
expected nonlocalities induced by galaxy formation and baryonic physics, all
of which contribute additively to the $\alpha_{n}$. Equivalently, the expected
contribution of these effects dictates the scales on which our perturbative
model is valid. We therefore put Gaussian priors on each counterterm centered
at zero with widths set such that their corrections are perturbative at our
chosen $k_{\rm max}$. We similarly include stochastic contributions which we
parametrize with SN${}_{0}=R_{h}^{3}$, SN${}_{2}=R_{h}^{3}\sigma_{2}$, and
SN${}_{4}=R_{h}^{3}\sigma_{4}$, where $R_{h}^{3}$ is the typical galaxy or
halo formation scale and the $\sigma_{n}$ arise from correlations of
stochastic modes in densities and velocities, (e.g. $\langle\delta
v\rangle,\langle v^{2}\rangle$, etc.). These stochastic terms again account
for the small-scale modes missing in perturbation theory, whose signs and
exact values are unknown, but whose rough size can be estimated based on our
understanding of the small-scale distribution and velocities of galaxies in
halos (see §4.2 and Appendix C and also Ref. [64]). These contributions are
added to the 1-loop power spectrum, $P_{s,g}^{PT}(\boldsymbol{k})$, above to
give our final LPT prediction

$\displaystyle P_{s,g}(\boldsymbol{k})=P_{s,g}^{PT}(\boldsymbol{k})$
$\displaystyle+(b+f\mu^{2})(b\alpha_{0}+f\alpha_{2}\mu^{2}+f\alpha_{4}\mu^{4})k^{2}P_{{\rm
s},b_{1}^{2}}(\boldsymbol{k})$
$\displaystyle+(\text{SN}_{0}+\text{SN}_{2}k^{2}\mu^{2}+\text{SN}_{4}k^{4}\mu^{4}),$
(3.6)

where $P_{{\rm s},b_{1}^{2}}$ is the term containing $b_{1}^{2}\xi_{\rm lin}$
in Eq. 3.5 evaluated to linear order outside of the exponential. This
parameterization of the counterterms differs slightly from previous works
using velocileptors . While giving consistent results, it makes it easier to
interpret the counterterms as “fractional corrections” to the linear theory
multipoles and motivates our choice of prior width on these parameters. For
example, a value of $\alpha_{n}=12.5\,h^{-2}\mathrm{Mpc}^{2}$ corresponds to a
$50\%$ correction to the $n^{\rm th}$ moment at $k_{\rm max}=0.20\,h{\rm
Mpc}^{-1}$.

In computing the observed power spectrum, we assume a fiducial cosmology to
convert $\mathbf{\theta}$ and $z$ to 3D distances using the fiducial distance-
redshift relation. We need to account for distortions in $P(k)$ between
assumed and true coordinates, the “Alcock-Paczynski (AP) effect” [65], in our
modeling. We do this by rescaling the theoretical power spectrum in true
cosmological coordinates to the observed coordinates by:

$\displaystyle P^{\rm obs}_{s}(\boldsymbol{k}_{\rm
obs})=q_{\perp}^{-2}q_{\parallel}^{-1}P_{s}(\boldsymbol{k})\quad,\quad k^{\rm
obs}_{\parallel,\perp}=q_{\parallel,\perp}\ k_{\parallel,\perp},\quad$ (3.7)

with the scaling parameters above are defined by444Previously in BOSS
analyses(e.g. [48, 52]) we have used the notation
$\alpha_{\parallel,\perp},\tilde{\alpha}_{\parallel,\perp}$ in place of
$q,\alpha_{\parallel,\perp}$ but in this paper we use the latter in order to
be consistent with the conventions of other DESI papers.:

$\displaystyle q_{\parallel}=\frac{H^{\rm ref}(z)}{H(z)}\quad,\quad
q_{\perp}=\frac{D_{A}(z)}{D^{\rm ref}_{A}(z)}\quad.$ (3.8)

$D_{A}(z)$ is the comoving angular diameter distance and the “ref” superscript
labels the values from the fiducial cosmology.

Finally, we use a Legendre transformation to compute the predicted power
spectrum as multipoles,

$\displaystyle P_{\ell}(k_{\rm obs})=\frac{(2\ell+1)}{2}\int_{-1}^{1}d\mu\
P(k,\mu_{\rm obs})\mathcal{L}_{\ell}(\mu)$ (3.9)

where ${\mathcal{L}}_{\ell}(\mu)$ is the Legendre polynomial of order $\ell$.

##  4 Fitting methods

###  4.1 Standard template and ShapeFit

The traditional parameter compression method used originally by the BOSS/eBOSS
collaborations involves choosing a reference cosmology,
$\boldsymbol{\Theta}^{\rm ref}$, and keeping the resultant linear power
spectrum, and by extension, the dependence on early-universe physics, fixed.
The “compressed” parameters being varied are then the amplitude,
$f\sigma_{s8}$ and the distance scalings transverse and along the line-of-
sight, $\alpha_{\perp},\alpha_{\parallel}$; all of which are only dependent on
late-time dynamics. The quantity $f\sigma_{s8}$, which controls the ratio of
monopole-to-quadrupole amplitudes, is a product of the growth rate,
$f\simeq\Omega_{m}^{0.55}$ and the total amplitude, $\sigma_{s8}$, at
$R=s\cdot 8\,h^{-1}$Mpc scales. Here $s=r_{\rm d}/r_{\rm d}^{\rm fid}$ with
$r_{\rm d}$ being the BAO scale at the drag epoch. We will comment on the $s$
scaling further below. The two distance scaling parameters are defined by,

$\alpha_{\parallel}=\frac{H^{\rm ref}(z)}{H(z)}\left(\frac{r_{\rm d}^{\rm
ref}}{r_{\rm d}}\right)=q_{\parallel}\left(\frac{r_{\rm d}^{\rm ref}}{r_{\rm
d}}\right)=\frac{q_{\parallel}}{s}\quad,\quad\alpha_{\perp}=\frac{D_{A}(z)}{D^{\rm
ref}_{A}(z)}\left(\frac{r_{\rm d}^{\rm ref}}{r_{\rm
d}}\right)=q_{\perp}\left(\frac{r_{\rm d}^{\rm ref}}{r_{\rm
d}}\right)=\frac{q_{\perp}}{s}\quad,$ (4.1)

We highlight that these parameters used in the template fitting are different
from the scaling parameters defined in eq. 3.8 by a factor of $\left(r_{\rm
d}^{\rm ref}/r_{\rm d}\right)$555Technically, this “ref” is not necessarily
the same as the “ref” in the definitions of $q_{\parallel,\perp}$. The one in
$\left(r_{\rm d}^{\rm ref}/r_{\rm d}\right)$ refers to the reference template
used in the standard template and ShapeFit fits, whereas in
$q_{\parallel,\perp}$ it refers to the fiducial cosmology assumed when
converting angles and redshift coordinates to physical distances when
measuring the power spectrum. However, in practice it is simplest to choose
the same cosmology for the template as was used for measuring the power
spectrum from the data, so this distinction is not important.. This is because
in the template method we assume that most information comes from the BAO
feature, and thus we account for the fact that both changes in $r_{d}$ and
$q_{\parallel,\perp}$ induce stretching in the observed BAO signal.666See
discussion in Appendix C of ref. [66], where however the pure AP parameters
are referred to as $\alpha$ and the BAO-rescaled ones are called
$\tilde{\alpha}$. In contrast, with a fitting method in which the underlying
cosmology is directly being varied (see next subsection), the changes to
$r_{d}$ affecting the BAO signal are automatically included in the linear
power spectrum which is self-consistently varied. We must also emphasize that
by including the factors of $s$ in our $\alpha$ scaling parameters we are
implicitly assuming distances in units of the BAO scale, which motivates our
use of the notation $f\sigma_{s8}$. This subtlety is discussed in detail in §
3 of Ref. [51].

Despite sacrificing constraining power through the lack of sensitivity to the
early universe (the shape of the transfer function is held fixed by the
reference cosmology), this “template” fitting method was sufficient at a time
when the tightest constraints on early-time physics came from the CMB and LSS
data was too noisy for direct fitting methods to be feasible without
significant priors from Planck. The advantages of the template fitting method
include the model-independence that allows for mapping the compressed
parameter constraints to a cosmological model of one’s choosing. Furthermore,
computing the linear power spectrum using a Boltzmann code such as CLASS or
CAMB at every step of a Markov Chain Monte Carlo (MCMC) sampler, in addition
to calculating nonlinear perturbation theory (PT) corrections, is
computationally very expensive. Fixing the linear power spectrum avoids this
step, allowing for a faster fitting procedure without needing to train an
emulator.

The “ShapeFit” method is an extension to the standard template-fit
compression, and was conceived as a way to partially bridge the gap in
constraining power between the standard template and direct/full modeling
methods, while preserving some of the model-independence of the former
technique [51]. This is achieved by allowing modifications to the shape of the
linear power spectrum via a multiplicative factor,

$\displaystyle P^{\prime}_{\rm lin}(\boldsymbol{k})=P_{\rm lin}^{\rm
ref}(\boldsymbol{k})\
\exp\left\\{\frac{m}{a}\tanh\left[a\ln\left(\frac{k}{k_{p}}\right)\right]+n\ln\left(\frac{k}{k_{p}}\right)\right\\},$
(4.2)

where $P_{\rm lin}^{\rm ref}(\boldsymbol{k})$ is the template power spectrum
produced by CLASS and is fixed throughout the fit. The modified power spectrum
$P^{\prime}_{\rm lin}(\boldsymbol{k})$ is what we provide to velocileptors to
produce the full 1-loop prediction for a given
$(f\sigma_{8},\alpha_{\parallel},\alpha_{\perp},m)$. For simplicity we keep
fixed the second shape parameter, $n=0$. Allowing this parameter to vary
accounts variations of the template emulating a spectral index effect, which
in this paper we do not consider. Following the original ShapeFit paper [51]
we choose for $a$ and $k_{p}$ their proposed values, $a=0.6$ and
$k_{p}=0.03\,h\,{\rm Mpc}^{-1}$.

###  4.2 Full modeling: $\Lambda$CDM and extensions

The alternative modeling technique to parameter compression is a more
conventional forward-modeling approach that involves directly varying the
underlying parameters of a cosmological model and making a theoretical
prediction for the observed quantities. While the $\Lambda$CDM model depends
on six parameters: ($\omega_{b}$, $\omega_{cdm}$, $H_{0}$,
$\log(10^{10}A_{s})$, $M_{\nu}$ and $n_{\rm s}$), some of these parameters are
not constrained by galaxy clustering analyses independently. For these
quantities we use priors derived from e.g. Big-Bang Nucleosynthesis (BBN)
and/or CMB anisotropies. We initially fix the spectral tilt and neutrino mass
to the Abacus fiducial values of ($n_{s}$,$M_{\nu}$) = (0.9649,0.06) – though
see Section 5.7. For the baryon abundance we adopt a narrow gaussian BBN prior
of $\mathcal{N}[\mu=0.02237,\sigma=0.00037]$ [67] (though see discussion in
Appendix D). Within these constraints, in this “Full-Modeling” approach the
shape of the linear power spectrum is able to change at each step of the MCMC
as the shape of the transfer function is dependent on the $\Lambda$CDM
parameters being varied. If done directly, this method is more computationally
expensive because the linear power spectrum must be calculated using a
Boltzmann code such as CLASS or CAMB in addition to the Velocileptors PT
corrections. However, through the use of an emulator we can efficiently and
accurately approximate the predictions for a given set of $\Lambda$CDM
parameters. Under the assumption that the predicted power spectrum multipoles
are a smooth function of the underlying parameters when close to some
reasonably chosen values, we can use an emulator based on a Taylor series
expansion in the relevant parameter space [41, 48].777In the event that the
data require a significantly different parameter space the analysis can be
iterated with the Taylor series recomputed closer to the best fit, assuming
the data are sufficiently constraining. We find that the emulator agrees well
with the direct LPT prediction when going to fourth order in the Taylor
expansion. After employing such an emulator both for the Full-Modeling and
template/ShapeFit methods, the MCMC chains converge (Gelman-Rubin
$|R-1|<0.01$) within roughly 1-2 hours888This is when using 8 parallel chains
on a single node. By analytically marginalizing of stochastic and counterterm
contributions (see Appendix A), the MCMC converges in 5-10 minutes for all
methods. Therefore, the improved computational efficiency of a compression is
no longer relevant in our setup.

The advantage of the Full-Modeling approach is that it is sensitive to both
the early-universe physics that determines the shape of the transfer function,
as well as late-time dynamics/geometry. Parameters such as $\omega_{b}$,
$\omega_{cdm}$, and $H_{0}$ affect both the early- and late- universe
dynamics, and are thus expected to be more tightly constrained in the Full-
Modeling approach, when compared to the methods employing a template that
fixes the early-universe dependence. On the other hand, the Full-Modeling
approach requires choosing a specific cosmological model from the start, and a
new MCMC fit is needed for any other model being employed. The parameter
compression methods, however, only require one fit, and afterwards the results
can be reused and mapped to any model of choice, though the model of choice
must be sufficiently close to the template cosmology unlike in the Full-
Modeling approach which does not suffer from this requirement.

We show in Table 1 the parameters and priors used for the Full-Modeling and
ShapeFit methods. We show the priors on bias parameters for three
parametrizations. The standard setting in this paper is the “intermediate”
freedom case for which the cubic bias is fixed to zero while
$(1+b_{1})\sigma_{8}$, $b_{2}\sigma_{8}^{2}$, and $b_{s}\sigma_{8}^{2}$ are
varied with Gaussian priors applied to the latter two. The other parameter
choices are discussed in Appendix D. We analytically marginalize over the
parameters controlling the stochastic and counterterm contributions, and refer
readers to Appendix A for further details and validation of this method.

Finally we remark that in order to make contact with earlier work, and in
particular with our companion papers, we use $\log(10^{10}A_{\mathrm{s}})$ as
the “normalization” of the power spectrum throughout. This choice, being the
normalization of the curvature power spectrum at $k=0.05\,\mathrm{Mpc}^{-1}$,
is actually better motivated for CMB surveys than galaxy redshift surveys.
Most of the constraining power of our data comes from quasi-linear scales and
we better constrain the matter power spectrum than the curvature (or
potential) power spectrum. In this respect a better choice for normalization
may be $\sigma_{8}$. We will discuss constraints on $\sigma_{8}$ later. We
also reiterate that the Full-modeling method does not require any re-scaling
of distances by $s=r_{\rm d}/r^{\rm ref}_{\rm d}$, and therefore the amplitude
being constrained here is $\sigma_{8}$ not $\sigma_{s8}$.

Full-Modeling | ShapeFit | Bias | Stoch/Counter
---|---|---|---
|  | Min. F. | Int. F.* | Max. F. |
H0 | $f\sigma_{8}$ | $(1+b_{1})\sigma_{8}$ | $\tilde{\alpha}_{0}$
$\mathcal{U}[55,79]$ | $\mathcal{U}[0,2]$ | $\mathcal{U}[0.5,3.0]$ | $\mathcal{N}[0,12.5]$
$\omega_{\mathrm{b}}$ | $\alpha_{\parallel}$ | $b_{2}\sigma_{8}^{2}$ | $\tilde{\alpha}_{2}$
$\mathcal{N}[0.02237,0.00037]$ | $\mathcal{U}[0.5,1.5]$ | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,12.5]$
$\omega_{\mathrm{cdm}}$ | $\alpha_{\perp}$ | $b_{s}\sigma_{8}^{2}$ | SN0
$\mathcal{U}[0.08,0.16]$ | $\mathcal{U}[0.5,1.5]$ | 0 | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,\mathcal{O}(1/\bar{n}_{g})]$
$\log(10^{10}A_{\mathrm{s}})$ | $m$ | $b_{3}\sigma_{8}^{3}$ | SN2
$\mathcal{U}[2.03,4.03]$ | $\mathcal{U}[-3.0,3.0]$ | 0 | 0 | $\mathcal{N}[0,5]$ | $\mathcal{N}[0,\mathcal{O}(f_{\rm sat}\sigma_{v}^{2}/\bar{n}_{g})]$
Table 1: Velocileptors LPT priors on parameters used in the Full-Modeling ($\Lambda$CDM) and ShapeFit fitting methods. The $\Lambda$CDM model involves H0, $\Omega_{\mathrm{b}}$, $\omega_{\mathrm{cdm}}$, $\log(10^{10}A_{\mathrm{s}})$ and all of the bias, stochastic, and counterterms. The ShapeFit method fits $f\sigma_{8}$, $\alpha_{\parallel}$, $\alpha_{\perp}$, $m$ as well as the same bias, stochastic and counterterms. The entries $\mathcal{U}[{\rm min,max}]$ and $\mathcal{N}[\mu,\sigma]$ refer to uniform and Gaussian normal distributions, respectively. For the bias terms we show both minimal, intermediate (standard), and maximal freedom cases, defined in Appendix D. For the two counterterms we report the priors within the parameterization for which the counterterms scale relative to the linear theory multipoles. The priors on the stochastic terms are given in Table 2 and discussed in the text. Tracer | $z_{\rm eff}$ | $1/\bar{n}_{g}$ | $f_{\rm sat}$ | $\log_{10}\bar{M}_{h}$ | $\sigma_{v}^{\rm est.}$ | SN0 | SN2 | SN4
---|---|---|---|---|---|---|---|---
LRG | 0.8 | 1000 | 0.1 | 13.3 | 7.8 | 2000 | $5.0\times 10^{4}$ | $1.0\times 10^{6}$
ELG | 1.1 | 300 | 0.1 | 11.9 | 2.9 | 1000 | $2500$ | $2.5\times 10^{4}$
QSO | 1.4 | 8000 | 0.03 | 12.7 | 5.7 | $1.5\times 10^{4}$ | $5.0\times 10^{4}$ | $1.0\times 10^{6}$
Table 2: Relevant quantities used for the prior widths of stochastic
parameters (see text). The typical halo mass, $\log_{10}\bar{M}_{h}$, per
galaxy is expressed in units of $h^{-1}M_{\odot}$ and $1/\bar{n}_{g}$ is
expressed in $h^{-3}\,\mathrm{Mpc}^{3}$. The characteristic velocities,
$\sigma_{v}^{\rm est.}$ are in $h^{-1}$Mpc. Motivated by these numbers, the
last three cloumns show the widths of the Gaussian priors (centered on 0 and
in $h^{-1}$Mpc units) that are used in this paper for each stochastic
parameter within each redshift bin. The results do not depend upon the precise
values chosen.

###  4.3 Cosmological inference from compressed statistics

In order to interpret the ShapeFit and standard template results, we must do
so in the context of a chosen cosmological model such as $\Lambda$CDM. While
it is simple to take a set of $\Lambda$CDM parameters and compute the
distances, $H(z)$, $D_{\rm A}(z)$, and $r_{\rm d}$ using CLASS or CAMB, in
order to compute compressed parameters assuming a certain fiducial cosmology,
it is more tricky in reverse [51]. Instead we must fit $\Lambda$CDM parameters
to the results of a fixed template fit with another MCMC. We take the chains
in the compressed parameters that were obtained from the initial template
fits, and compute the parameter mean vector and covariance matrix, i.e.
$\bar{\bf\Theta}=(\bar{f\sigma_{8}},\bar{\alpha_{\parallel}},\bar{\alpha_{\perp}},\bar{m})$
and C4×4. Treating $\bar{\bf\Theta}$ and C4×4 as a “data” vector and
associated covariance, we can now sample in $\Lambda$CDM parameters so that
for each proposed set of ($\omega_{b}$, $\omega_{cdm}$, $h$, $\log A_{s}$) we
compute the corresponding vector $\boldsymbol{\Theta_{\rm
thy}}=(f\sigma_{8},\alpha_{\parallel},\alpha_{\perp},m)_{\rm thy}$. Assuming
all compressed parameters are Gaussian, we then use an MCMC to sample from the
likelihood,

$\mathcal{L}\propto\exp\left\\{-\frac{1}{2}(\boldsymbol{{\Theta}_{\rm
thy}}-\boldsymbol{\bar{\Theta}})^{T}\mathbf{C}_{4\times
4}^{-1}(\boldsymbol{{\Theta}_{\rm thy}}-\boldsymbol{\bar{\Theta}})\right\\}.$
(4.3)

When inferring cosmological constraints from the ShapeFit parameters, care
must be taken in interpreting the amplitude $f\sigma_{s8}$ appropriately, as
the slope rescaling via the $m$ parameter also changes $\sigma_{s8}$. As noted
in refs. [51, 68], the parameter $f$ that is varied in ShapeFit analyses is
actually $fA\equiv f(A_{sp}/A_{sp}^{\rm ref})^{1/2}$, where
$A_{sp}=s^{-3}P_{\rm no-wiggle}^{\rm lin}(k_{p}/s,\boldsymbol{\Theta})$ is the
amplitude of the no-wiggle power spectrum at the pivot scale, $k_{p}\simeq
0.03h$Mpc-1. The parameter $s$ describes the scaling of lengths relative to
the BAO and is defined to be the ratio $r_{\rm d}/r^{\rm ref}_{\rm d}$. In
order to generate the model 1-loop power spectrum multipoles, we must provide
velocileptors with the linear power spectrum $P^{\prime}_{\rm
lin}(\boldsymbol{k})$ from Eq. 4.2 and the growth factor $f$. Defining LPT_RSD
as the function that produces the power spectrum multipoles, the nearly exact
degeneracy between $f$ and the power spectrum amplitude (see § 5.9) implies
that

$\displaystyle\texttt{LPT\\_RSD}\left[f\times\left(\frac{A_{sp}}{A_{sp}^{\rm
ref}}\right)^{1/2};P^{\prime}_{\rm
lin}(\boldsymbol{k})\right]\leftrightarrow\texttt{LPT\\_RSD}\left[f;\left(\frac{A_{sp}}{A_{sp}^{\rm
ref}}\right)\times P^{\prime}_{\rm lin}(\boldsymbol{k})\right],$ (4.4)

and thus the true $f\sigma_{s8}$ is given by

$\displaystyle f\sigma_{s8}$
$\displaystyle=f\times\left[\left(\frac{A_{sp}}{A_{sp}^{\rm
ref}}\right)\int\frac{dk}{2\pi^{2}}k^{2}\tilde{W}_{R}^{2}(kR)P^{\prime}_{\rm
lin}(\boldsymbol{k})\right]^{1/2}$
$\displaystyle=f\times\left(\frac{A_{sp}}{A_{sp}^{\rm
ref}}\right)^{1/2}\left[\int\frac{dk}{2\pi^{2}}k^{2}\tilde{W}_{R}^{2}(kR)P_{\rm
lin}(\boldsymbol{k})\exp\left\\{\frac{m}{a}\tanh\left[a\ln\left(\frac{k}{k_{p}}\right)\right]\right\\}\right]^{1/2}$
(4.5) $\displaystyle\simeq\frac{(f\sigma_{s8})^{\rm ref}}{(fA_{sp}^{1/2})^{\rm
ref}}fA_{sp}^{1/2}\times\exp\left\\{\frac{m}{2a}\tanh\left[a\ln\left(\frac{r_{d}^{\rm
fid}}{\rm R}\right)\right]\right\\}.$ (4.6)

Here $R$ is the smoothing scale of the amplitude parameter $\sigma_{R}$ and is
chosen to be $R=8\,h^{-1}$Mpc by convention. There are now two ways in which
one could use ShapeFit chains in order to infer about cosmological parameters:
one can use the above equations (either the exact or approximate forms) to
transform the sampled $fA$ chain into $f\sigma_{s8}$, and then use CLASS to
compute $f\sigma_{s8}$ for every set of $\Lambda$CDM parameters at the
interpretation step; or one can directly perform the interpretation on $fA$ by
always computing $f$ and $A_{sp}$ while sampling in $\Lambda$CDM parameters.
We find that the two approaches give consistent constraints in the
$\Lambda$CDM parameter space.

Finally, the $m$ parameter in ShapeFit that controls the shape of the linear
power spectrum can be computed from $\Lambda$CDM parameters through the ratio
[51]

$m=\left.\frac{d}{dk}\left(\ln\left[\frac{T(k_{p}/s,\boldsymbol{\Theta})}{T(k_{p},\boldsymbol{\Theta}^{\rm
ref})}\right]\right)\right|_{k=k_{p}}\quad,\quad
T(\boldsymbol{\Theta},k)=\frac{P_{\rm no-wiggle}^{\rm
lin}(k,\boldsymbol{\Theta})}{\mathcal{P}_{\mathcal{R}}(k,\boldsymbol{\Theta})},$
(4.7)

with primordial power spectrum $\mathcal{P}_{\mathcal{R}}$.

##  5 Results

Figure 2: 1D posteriors from the Full-Modeling fit as the covariance volume is
varied from that of a single box (8 $[\,h^{-1}{\rm Gpc}]^{3}$) to 15 boxes
((120 $[\,h^{-1}{\rm Gpc})^{3}$])

Before we present the results from the various systematic tests of
velocileptors and the different modeling methods, we first revisit the issue
of covariance volume. In Fig. 2 we present 1D posterior constraints from the
Full-Modeling fit to LRG mock data as a function of covariance volume, i.e.
multiples of the single-box volume such that the covariance is rescaled by
$1/n,\,n=1,3,5,\cdots,15$. We show results for fits using two different
$k$-ranges, $0.02\leq k\,[\,h{\rm Mpc}^{-1}]\leq 0.18$ and $0.02\leq
k\,[\,h{\rm Mpc}^{-1}]\leq 0.20$ (which will be our ‘standard’ range). We find
that as the volume is increased, the constraints in $\Omega_{\rm m}$ shift
towards the truth as the error bars tighten, which is indicative of a prior
volume effect. For $H_{0}$ and $\log(10^{10}A_{s})$ the constraints remain
mostly stable as the volume is increased, with small shifts increasing with
volume that likely relate to the increasing sensitivity to two-loop effects
that are not included in the model. For similar reasons, we observe a
divergence in constraints between $k_{\rm max}=0.18$ and $k_{\rm
max}=0.20\,h{\rm Mpc}^{-1}$ that grows as the volume is increased. This shows
that when using an ultra-tight covariance such as that of the $200\,h^{-1}{\rm
Gpc}$ simulation volume, one can expect $1-2\sigma$ offsets in constraints
arising purely from theoretical errors due to the limited number of terms
included in the 1-loop power spectrum model. In addition, as mentioned
earlier, the N-body simulations themselves have systematic errors that become
important at these volumes and can contribute to the shifts we observe.

###  5.1 Baseline Comparison

Figure 3: (left): Comparison of constraints of compressed parameters for the
standard template method vs. ShapeFit. (right): Comparison of constraints on
$\Lambda$CDM parameters for the standard template, ShapeFit, and full modeling
methods. The single-box covariance is used for these results, with our
‘standard’ $k_{\rm max}=0.2\,h\,\mathrm{Mpc}^{-1}$ (see §5.2 for the
discussion of $k_{\rm max}$ dependence).

We begin with comparing constraints in the compressed parameter space between
the standard template and ShapeFit approaches, using the single-box
covariance, as shown in the left panel of Fig. 3. We see that the posterior
means of the two methods agree very closely, with slightly smaller contours
for the standard template due to varying fewer parameters. Since the reference
template used in these fits is the true Abacus cosmology, we expect
$\alpha_{\parallel,\perp}=1$ and $m=0$. In both cases, the means of all
parameters are within 1$\sigma$ of the expected values. When interpreting
these results in terms of a $\Lambda$CDM cosmology, however, we see a
significant difference in the constraints from the two compression methods
(right panel of Fig. 3). While both methods give unbiased constraints on
$\Lambda$CDM parameters (within 1$\sigma$ of truth) the error bars for all
parameters are significantly larger for the template case due to the lack of
information from the power spectrum shape in the template approach. This is
expected, as the template method was traditionally combined with external data
sets under the assumption that the parameters determining the shape are not as
well constrained from LSS data than e.g. CMB anisotropies, but in our setup we
rely purely on the LSS data alone(but see §5.6). Meanwhile, when comparing the
constraints between the ShapeFit and Full-Modeling methods, we find a very
close agreement in the shape and orientations of the contours, showing that
the ShapeFit method is able to match the constraining power of direct model
fitting, at least for the $\Lambda$CDM case for which it was designed. We do
observe mild differences in the tightness of constraints between the ShapeFit
and Full-Modeling methods. These could be due to a combination of various
approximations in the ShapeFit method, such as controlling the shape of
$P_{\rm lin}$ with only one parameter and assuming the compressed parameters
to be perfectly Gaussian in the interpretation step.

###  5.2 Dependence on $k_{\rm max}$

Figure 4: Dependence on $k_{\rm max}$ of the Full-Modeling (upper panel) and
ShapeFit (lower panel) methods for the three tracer types: LRG ($z=0.8$). ELG
($z=1.1$), and QSO ($z=1.4$) (points slightly offset for clarity). The single-
box covariance is used for all of these fits.

We next test the dependence on scale cuts of our model, for the different
methods. In all cases we fix the lower bound of the $k$-range to $0.02h$Mpc-1.
This is fully in the linear regime so the the stability of the theory is not
affected by the specific value chosen, but this choice simply removes points
too close to the fundamental mode of the cubic box ($k=0.003\ \,h{\rm
Mpc}^{-1}$). We then run our fits with upper bounds of $k_{\rm
max}=0.16-0.26\,h{\rm Mpc}^{-1}$. The results are shown for Full-Modeling and
ShapeFit in Fig. 4 for the LRG, ELG, and QSO tracers. The higher $k$-modes,
above $\sim 0.2\,h\,\mathrm{Mpc}^{-1}$, correspond to smaller scales which are
more sensitive to nonlinear effects and galaxy/halo formation physics, which
are not well-understood and therefore difficult to model. Our model includes
non-linearities only at the 1-loop level and bias only up to cubic order. We
therefore expect biases to worsen as higher $k$-modes are included in the fit.
For the single-box volume we find the two methods to remain relatively stable
as $k_{\rm max}$ is increased as the observational errors match or exceed the
theoretical or modeling errors, however we do observe $\gtrsim 1\sigma$
offsets in the $\sigma_{8}$ constraints for LRG and ELG tracers in the Full-
modeling method when $k_{\rm max}\geq 0.22\ \,h{\rm Mpc}^{-1}$. In Fig. 5 we
repeat this test for the LRG tracers but using the 25 box covariance. We show
constraints in the $\Lambda$CDM as well as ShapeFit parameter spaces. In this
case we obtain significantly biased constraints when $k_{\rm max}>0.2\,h{\rm
Mpc}^{-1}$. In the $\Lambda$CDM parameters, we find a mild improvement in
constraining power of Full-Modeling at k${}_{\rm max}=0.18\,h{\rm Mpc}^{-1}$
versus our usual setting of $k_{\rm max}=0.20\,h{\rm Mpc}^{-1}$. This
worsening of constraints when $k_{\rm max}$ is increased is likely due to a
sensitivity to higher-order effects that our theory does not adequately
describe, and which become increasingly important with increasing $k$. When
using an extremely tight covariance, the additional high-$k$ points push the
fit towards incorrect models and away from the constraints coming from low-$k$
data points. In the compressed parameter space we observe slightly more
significant offsets ($\gtrsim 1.\sigma$) in the $\alpha_{\perp}$ and
$f\sigma_{s8}$ constraints for ShapeFit at k${}_{\rm max}=0.18\,h{\rm
Mpc}^{-1}$. When deriving summary statistics from the Full-Modeling
constraints, the $\alpha_{\parallel}$, and $\alpha_{\perp}$ parameters are
significantly more tightly constrained than in the ShapeFit and Template
methods because the $\Lambda$CDM priors in Full-Modeling restrict the
allowable values that the scaling parameters can take [52]. We use the results
from Figs. 4-5 to motivate a choice of $k_{\rm max}=0.20\ \,h{\rm Mpc}^{-1}$
as our baseline analysis setting, as this is the largest $k_{\rm max}$ for
which all three modeling methods are acceptably close to truth ($\lesssim
1\sigma$ offsets) in the $\Lambda$CDM parameter space in both the single-box
and full covariance volume cases.

Figure 5: Constraints on $\Lambda$CDM and compressed parameters from the three
modeling methods with varying $k_{\rm max}$ between $0.16\leq k_{\rm max}\leq
0.26\ \,h{\rm Mpc}^{-1}$(points slightly offset for clarity). These results
are obtained with a covariance appropriate to the 25 box volume, fitting to
the LRG cubic mock data.

As we proceed to the remainder of tests presented in this paper, we refer
readers to Fig. 6 for a summary figure of 1D constraints on $\Omega_{\rm
m},H_{0},$ and $\log 10^{10}A_{s}$ obtained from each of the tests.

###  5.3 Joint fitting of LRG, ELG, and QSO mocks

Figure 6: Constraints on $\Lambda$CDM parameters for the three modeling
methods for a variety of different fit settings and data sets to be discussed
in the text. The results are obtained with the covariance for a single box, 8
$(\,h^{-1}{\rm Gpc})^{3}$, volume and $k_{\rm max}=0.2\,h\,\mathrm{Mpc}^{-1}$.
In many of the cases $\Omega_{\rm m}$ appears slightly below the truth, which
is in part due to projection effects. Here “standard” refers to our baseline
result on the LRG mocks.

We now turn to the joint fitting of data samples from different tracers and
redshift bins. The three tracers are Luminous Red Galaxies (LRG, $z=0.8$),
Emission Line Galaxies (ELG, $z=1.1$), and Quasars (QSO, $z=1.4$). For the
Full-Modeling case, we still sample in $\Lambda$CDM parameters as usual but
compute separate $P_{\ell}(k)$ models for each redshift bin and the likelihood
is computed from all data sets, i.e. the data vector becomes
$\boldsymbol{d}=(P_{0}^{LRG},P_{2}^{LRG},P_{0}^{ELG},P_{2}^{ELG},P_{0}^{QSO},P_{2}^{QSO})$.
This results in a total effective volume of 600 $(\,h^{-1}{\rm Gpc})^{3}$. We
do not assume any correlation between tracers at different redshifts999The
mean data vectors for the LRG, ELG, and QSO tracers actually came from the
same 25 realizations and therefore share initial conditions. In principle this
means that the redshift bins are not truly uncorrelated, but we assume so in
this work for simplicity., so the total joint covariance matrix has zeros in
the indices corresponding to cross correlations between different tracers.
This ensures that contributions to the log-likelihood such as $\Delta
P_{i}^{LRG}C_{ij}^{-1}\Delta P_{j}^{ELG}=0$. We use a separate set of nuisance
parameters for each type of tracer. For the standard template and ShapeFit
fits, the free parameters ($f\sigma_{s8}$,$\alpha_{\parallel,\perp},m$) are in
general redshift dependent. While in principle one could use a single
$f\sigma_{s8}(z=0)$ as a free parameter and then rescale by the fiducial
growth factor $D(z,\Omega_{m}=\Omega_{m}^{\rm fid})$ in order to get the
corresponding parameter for the different samples, the redshift dependence of
the $\alpha$’s and $m$ parameters is not as obvious. Instead we perform the
parameter compression separately for the LRG, ELG, and QSO samples and obtain
three sets of $(f\sigma_{s8}$,$\alpha_{\parallel,\perp},m)_{z}$ to be used as
“summary statistics” of each tracer sample. It is in the cosmological
interpretation step that we can either infer $\Lambda$CDM parameters from a
single sample or from the combination of
$(f\sigma_{s8}$,$\alpha_{\parallel,\perp},m)_{z}$ sets of multiple tracer
samples.

In the three panels of Fig. 7 we show a comparison between results of fitting
a single sample versus joint fits of multiple tracers, for the standard
template, ShapeFit, and Full-Modeling methods respectively. We observe that in
each method, the ELG data is significantly more constraining than the LRG
sample, and thus the joint fitting constraints appear to be dominated by the
ELG sample. The QSO mocks are the least constraining data set, due to the
lower number density of Quasars from which the power spectrum is measured.
Therefore the error bars at each Fourier mode are larger than those of the ELG
and LRG data, resulting in significantly poorer constraints in the model
parameters governing the power spectrum shape, i.e. $\Omega_{\rm m}$ and
$H_{0}$. Meanwhile, the amplitude parameter $\log A_{s}$ is not as sensitive
to the type of tracer and we observe smaller differences in constraint between
the tracer types. Overall, the tightest constraints on all parameters are
obtained in the joint analysis of LRG+ELG+QSO, but with an almost negligible
improvement coming from the inclusion of QSO data.

(a) Standard Template

(b) ShapeFit

(c) Full-Modeling

(d) All Methods

Figure 7: Comparison of cosmological constraints of different tracers (LRG,
$z=0.8$; ELG, $z=1.1$; QSO, $z=1.4$) for the three fit methods. Here we only
show the results with the covariance of the single-box volume of
$(2\,h^{-1}{\rm Gpc})^{3}$ and $k_{\rm max}=0.2\,h\,\mathrm{Mpc}^{-1}$.

###  5.4 Full-shape + BAO Reconstruction

In addition to fitting the full-shape power spectra using our model, we can
gain extra constraining power through a joint analysis with the reconstructed
BAO correlation function. The BAO reconstruction procedure aims to undo some
of the damping of the BAO signal due to nonlinear structure growth in order to
sharpen its peak, allowing for a better measurement of the cosmological
distance-redshift relation via the well-defined drag horizon scale (see e.g.
refs. [69, 70, 71, 72, 73, 74]). This procedure begins by smoothing the
observed clustering signal by a Gaussian filter $S(k)=\exp(-(kR_{s})^{2}/2)$,
which serves to filter out small-scale modes. Next, we use this smoothed
density to estimate the smoothed Zel’dovich displacement,
$\boldsymbol{\chi}\approx S(k)\boldsymbol{\Psi}_{\rm Zel}$, which we subtract
from the observed galaxy field as well as from a random matter density field
in order to preserve large-scale power. The reconstructed galaxy density field
is then $\delta_{\rm rec}=\delta_{\rm d}-\delta_{\rm s}$, with $\delta_{\rm
d}$ and $\delta_{\rm s}$ being the displaced galaxy and shifted random fields,
respectively. Moving to redshift space once again amounts to a rotation of the
real-space field, $\boldsymbol{\chi}_{s}=R\boldsymbol{\chi}$ with matrix $R$
defined in Sect. 3. In the literature one commonly encounters two methods for
reconstructions in redshift space: RecSym [72] and RecIso [69, 71]. The first
applies the transformation into redshift space equally to both $\delta_{\rm
d}$ and $\delta_{\rm s}$, whereas the latter method keeps the shifted field in
real-space (see ref. [74] for further discussion). For the DESI simulations
considered in this work, the RecSym procedure is applied to produce the post-
reconstruction mock data.

We model the damping of the BAO feature in the reconstructed power spectrum,
$P_{\rm rec}=P_{\rm dd}+P_{\rm ss}-2P_{\rm ds}$ within the Zel’dovich
approximation by splitting the linear theory predictions into the wiggle and
no-wiggle components101010There are numerous methods for performing this
split. Here we use the method described in Appendix D of ref. [75] that uses a
sine transform to identify the BAO feature in real space and subtracts it
before transforming back to Fourier space to produce a wigge-free power
spectrum. and apply an exponential damping factor111111Previous works studying
BAO reconstruction have sometimes derived different damping factors for
$P_{dd}$, $P_{ds}$ and $P_{ss}$. This results from a $1^{\rm st}$ order
approximation in LPT, and a more consistent approach has the randoms damped by
the same factor. This subtlety is described in detail in ref. [74], as well as
in ref. [76] for a slightly different reconstruction scheme. However, we find
that the difference between the old and new methods results in negligible
effects to the fit posteriors. to the wiggle part [74]

$\displaystyle
P(k,\mu)=(b+f\mu^{2})^{2}\left(P^{NW}(k)+e^{-\frac{1}{2}k^{2}\Sigma^{2}(\mu)}P_{ab}^{W}(k)\right),$
(5.1)

where the $\Sigma^{2}$ in the damping factor is the isotropic component of the
linear pairwise displacement
$A_{ij}^{dd}=\left<\Delta_{i}^{dd}\Delta_{j}^{dd}\right>$, of the displaced
density field at $|\boldsymbol{q}|=r_{d}$, i.e.

$\displaystyle\Sigma^{2}(\mu)$
$\displaystyle=\left.\frac{1}{3}\delta_{ij}A_{ij}(\boldsymbol{q})\right|_{q=r_{d}}$
(5.2)
$\displaystyle=\left[1+f(2+f)\mu^{2}\right]\left[2\tilde{\Sigma}^{2}(0)-2\tilde{\Sigma}^{2}(r_{d})\right]$
(5.3) $\displaystyle\tilde{\Sigma}^{2}(q)$
$\displaystyle=\frac{1}{3}\int\frac{dk}{2\pi^{2}}(1-\mathcal{S})^{2}j_{0}(kq)P_{\rm
lin}(k)\quad.$ (5.4)

Finally, after generating the reconstructed power spectrum, we use a Fourier
transform to obtain the reconstructed correlation function. We limit our model
to linear bias as it has been found in previous works that the IR damping of
the BAO feature dominates over other nonlinear effects such as mode-coupling
which are largely cancelled by reconstruction. Following Ref. [74] we employ a
new method for modeling the broadband that is not degenerate with the BAO
signal, which in Fourier space involves using a basis of cubic splines. When
fitting the correlation function in configuration space this is equivalent to
setting a minimum scale, $r_{\rm min}$, with the exception of two Hankel
transformed basis functions that are included in the quadrupole:

$\displaystyle\mathcal{Q}_{2,n}(r\Delta)=\frac{i^{2}}{2\pi^{2}}\int
dkk^{2}W_{3}\left(\frac{k}{\Delta}-n\right)j_{2}(kr),\quad n=0,1$ (5.5)

where $W_{3}$ is the piecewise cubic spline kernel [77, 78], $j_{2}$ is a
$\nu=2$ spherical Bessel function, and we choose $\Delta=0.06\,h{\rm
Mpc}^{-1}$ for the separation scale of the splines. We additionally include a
template of polynomials in even powers of $r$ for the monopole and quadrupole
moments, truncated at quadratic order, to marginalize over contamination by
large-scale systematics below some $k_{\rm min}$. The broadband model in
configuration space is thus [74]:

$\displaystyle\mathcal{B}_{0}(r)=a_{0,0}+a_{0,1}\left(\frac{rk_{\rm
min}}{2\pi}\right)^{2}$
$\displaystyle\mathcal{B}_{2}(r)=a_{2,0}+a_{2,1}\left(\frac{rk_{\rm
min}}{2\pi}\right)^{2}+\Delta^{3}(a_{2,2}\mathcal{Q}_{2,0}(r\Delta)+a_{2,3}\mathcal{Q}_{2,1}(r\Delta))$
(5.6)

where $k_{\rm min}=0.02\,h{\rm Mpc}^{-1}$ and the parameters
{$a_{0,0},a_{0,1},a_{2,0},a_{2,1},a_{2,2},a_{2,3}$} can be analytically
marginalized over. We use broad Gaussian priors centered at 0 with widths of
$5\times 10^{5}$ for all of these broadband parameters. Finally, we note that
one should also include some more flexibility in the damping factor by
introducing parameters $\Sigma_{\parallel,\perp}$ in the exponent in Eq. 5.1
to marginalize over the effects of nonlinearities. However, we did not find
this necessary in the tests presented here, and so the damping factors vary
only as $f,P_{\rm lin},$ and $r_{\rm d}$ change in Full-modeling and likewise
with ShapeFit through the $f\sigma_{s8}$ and $m$ parameters.

The joint covariance matrix is computed numerically using the reconstructed
correlation function realizations of the EZmock simulations. So the joint data
vector is now $\textbf{d}=\\{P_{0}^{\rm pre}(k),P_{2}^{\rm
pre}(k),\xi_{0}^{\rm post}(r),\xi_{2}^{\rm post}(r)\\}$ with cross-
correlations between $P_{\ell}^{\rm pre}(k)$ and $\xi_{\ell}^{\rm post}(r)$
accounted for as nonzero off-diagonal elements in the joint covariance matrix.
(e.g. see Fig. 3 of [48])

We show in Fig. 8 comparisons of the cosmological constraints pre/post BAO
reconstruction. We find that for all three modeling methods there is
significant improvement in constraints when joint-fitting with the post-recon
correlation function, most significantly in $H_{0}$ as the cleaner measurement
of BAO scale from the sharpened peak allows for better calibration of the
distance-redshift relation that constrains Hubble’s constant. When comparing
all methods we find consistent constraints between ShapeFit and Full-Modeling
that are both tighter than those of the standard template.

(a) Standard Template

(b) ShapeFit

(c) Full-Modeling

(d) Post Recon (all methods)

Figure 8: Comparison of cosmological constraints with and without BAO
reconstruction for each modeling method. The bottom right panel compares the
post-reconstruction constraints of the three methods. For all plots above, we
present results using the covariance appropriate to the single-box volume.

###  5.5 Beyond $\Lambda$CDM: $w$CDM model

Figure 9: On the left we show a comparison of cosmological constraints on the
$w$CDM parameters for the three modeling methods. On the right we show $w$CDM
Full-Modeling fits with and without the inclusion of post-reconstruction BAO
data. In both plots we show the results with the single-box $(2\,h^{-1}{\rm
Gpc})^{3}$ volume.

With the expected improvement in cosmological parameter estimation from future
galaxy redshift surveys, we hope to place better constraints on parameters not
just underlying the standard $\Lambda$CDM model, but also departures from it.
From the Friedmann equations, the energy density of a specific component of
the Universe is related to the scale factor, $a$, by

$\displaystyle\rho\propto a^{-3(1+w)}$ (5.7)

where $w=p/\rho$ is the equation of state parameter. One of the simplest
extensions to $\Lambda$CDM involves allowing the dark energy equation of state
to differ from the value of $-1$ that corresponds to a cosmological constant
($\Lambda$) as the energy density is constant in that case. On the other hand,
“quintessence” models have $w\neq-1$ such that dark energy is a dynamic
quantity in the Universe.121212If dark energy is described by a scalar field,
$\phi$, with a canonical kinetic term then the equation of state can be
interpreted in terms of kinetic and potential energies via,  $\displaystyle
w=\frac{\frac{1}{2}\dot{\phi}^{2}-V(\phi)}{\frac{1}{2}\dot{\phi}^{2}+V(\phi)}\quad.$
(5.8) Under this assumption the equation of state is usually expected to lie
between $-1<w<1$, with values $w<-1/3$ leading to cosmic acceleration.
However, more exotic models exist that do allow for negative kinetic energies.

Fig. 9 shows in the left panel the constraints on $w$CDM parameters obtained
from each of the three modeling methods, for the covariance of the single-box
volume. Since the Abacus cosmology assumes a cosmological constant for dark
energy, the expected value is $-1$. We find that the ShapeFit and Full-
Modeling methods both give constraints on $w$ that are within $1\sigma$ of the
expected equation of state. Meanwhile the parameters in the template method
are very poorly constrained when $w$ is varied. When changing the properties
of dark energy away from the cosmological constant the universe’s expansion
history and geometry are significantly altered, thus affecting the
$\alpha_{\parallel,\perp}$ parameters and $f\sigma_{s8}$. This results in the
observed degeneracies between $w$ and the other parameters (which also
determine $f\sigma_{s8}$, and $\alpha_{\parallel,\perp}$). If those three
parameters are the only information we have from the data, as is the case in
the template fit, then this results in very poor constraints. However, moving
far along those degeneracy directions also significantly affects the shape of
the power spectrum, which the ShapeFit and Full-Modeling methods _are_
sensitive to. Therefore these two methods do not suffer from the degeneracies
as much as the template fit. Comparing ShapeFit to Full-Modeling, we find that
the constraints on parameters from the ShapeFit method are a bit wider than in
Full-Modeling. This is likely because all of the shape information is
contained in a single parameter, which then needs to be interpreted as
constraints on three different cosmological parameters ($w,\omega_{\rm m},$
and $H_{0}$), as these all control the shape of the power spectrum. Thus, a
poorer measurement of $m$ results in more sensitivity to the degeneracies in
shape that the template fit also suffered from. Finally, we also note that
projection effects (see Appendix B) in Full-Modeling cause close-to 1$\sigma$
offsets in the $w,\Omega_{\rm m},$ and $H_{0}$ parameters. While these shifts
are not huge for this dataset, we also are interested into what extent
including more data can mitigate projection effects. We show in the right
panel of Fig. 9 a comparison of Full-Modeling fits with and without the
inclusion of reconstructed BAO data. We find that including BAO results in
noticeable improvements in the constraints by shifting the posteriors closer
to the truth. These projection effects are not as significant in the ShapeFit
method, which suggests that the extra information that Full-modeling obtains
w.r.t. ShapeFit may come from regions of the power spectrum that are
degenerate with counterterm and/or stochastic parameters. A similar effect was
observed and reported in Ref. [52] when comparing $f\sigma_{s8}$ constraints
between Full-Modeling and standard template methods in BOSS data.

###  5.6 Priors from CMB

(a) Full Covariance

(b) Single Box Covariance

Figure 10: Comparison of $\Lambda$CDM constraints from the template fit(blue),
Full-Modeling (green), and ShapeFit (red) with Planck priors, on the Abacus
LRG ($z=0.8$) mock data. On the left we show the fits with the covariance for
the full ($25\times(2\,h^{-1}{\rm Gpc})^{3}$) volume while the right figure
shows the results with the single-box $(2\,h^{-1}{\rm Gpc})^{3}$ volume. In
the $\Omega_{\rm cdm}h^{2}-\Omega_{\rm b}h^{2}$ panels we include black lines
showing the Planck prior.

The ‘standard’ template method was conceived at a time when the data from
galaxy redshift surveys was not constraining enough on early-universe physics
to be competitive with constraints from probes such as Planck that modeled CMB
anisotropies. In particular, data from CMB anisotropies tightly constrain the
$\Lambda$CDM parameters that determine the shape of the power spectrum [79],
and this shape is left unaltered by late-time physics such as dark energy or
spatial curvature. These constraints are tighter than those from the galaxy
survey themselves. In such a scenario, the primary degrees of freedom to be
constrained by galaxy surveys are late time growth and the late-time distance-
redshift relation. The template method was intended to be used in conjunction
with the other probes, such that most of the information on $P(k)$ shape came
from strong priors using results from e.g. Planck. To demonstrate this, we
repeat the cosmological inference of the template results, but including an
additional likelihood derived from the Planck 2018 results [80]. We do this by
taking the chains from the baseline model of the Planck Legacy Archive, “base
plikHM TT lowl lowE”, and compute the covariance matrix, $C_{\rm pl}$, from
the ($\omega_{\rm b}$, $\omega_{\rm cdm}$) samples. We do not apply a prior on
$A_{s}$ or $H_{0}$ as we are interested in how information from galaxy
clustering constrains the late-time growth compared to Planck. When we sample
in these $\Lambda$CDM parameters we now include the additional likelihood

$\displaystyle\mathcal{L}_{\rm
pl}\propto\exp\\{-\frac{1}{2}\boldsymbol{\Delta\Theta}^{T}C_{\rm
pl}^{-1}\boldsymbol{\Delta\Theta}\\},$ (5.9)

where $\boldsymbol{\Delta\Theta}$ is the difference between the sampled
($\omega_{\rm b}$, $\omega_{\rm cdm}$) and the values in the Abacus cosmology.
Because we are including the CMB prior on $\omega_{b}$, we remove the BBN
prior that we usually use in our standard analyses. We show these results,
comparing the template, Shapefit, and Full-Modeling methods with Planck
priors, in Fig. 10, using the LRG ($z=0.8$) mock data within the standard
$\Lambda$CDM model. We see that the inclusion of Planck priors significantly
tightens the constraints on $\Omega_{\rm m}$. Despite us not applying any
prior on $H_{0}$ and $\log A_{s}$, we still observe a shift to the truth and
tightening in those parameter constraints for all three methods, with the
$\log A_{s}$ posterior slightly narrower for the Full-Modeling approach.
Overall, all three methods agree very closely in all of the parameters when
including these priors, suggesting that the difference in constraining power
of these methods is almost entirely due to shape information (which is better
determined by the CMB than the galaxy survey).

###  5.7 Varying $n_{s}$

(a) LRG ($0.4<z<0.6$)

(b) All Y1 bins

Figure 11: The left panel shows Full-Modeling fits to the mean of Cutsky mocks
in the DESI Y1, LRG $0.4<z<0.6$ redshift bin, with different priors on
$n_{s}$. We use the notation $\mathcal{U}[\mathrm{min},\mathrm{max}]$ and
$\mathcal{N}[\mu,\sigma]$ to denote uniform and Gaussian priors, respectively.
In the right panel we compare Full-Modeling $n_{s}$ fixed versus free using
the synthetic mocks created with velocileptors simulating the the full Y1
footprint: BGS ($0.1<z<0.4$), LRG ($0.4<z<0.6$), LRG ($0.6<z<0.8$), LRG
($0.8<z<1.1$), ELG ($0.8<z<1.1$), ELG ($1.1<z<1.6$), and QSO ($0.8<z<2.1$). We
show constraints with uniform priors on $n_{s}$ in the free case.

For previous fullshape analyses from spectroscopic surveys, it was
common/necessary to fix (or impose tight priors) on several of the
$\Lambda$CDM parameters such as $\omega_{b}$, $M_{\nu}$, and $n_{s}$, using
information from the CMB and BBN. With the increasing constraining power of
DESI and future surveys it is of interest to see how much we can untangle
fullshape analyses from other probes. While a tight prior on $\omega_{b}$ (see
Appendix D) is still necessary, the improved constraining power of DESI may
allow us to free $n_{s}$ and/or $M_{\nu}$131313In this paper we only perform
tests with $n_{\rm s}$ free and refer readers to Ref. [44] for a discussion on
varying $M_{\nu}$.. To investigate the impact of uncertainty in $n_{s}$ on our
analysis given the statistical uncertainties in Y1, we chose mock data from
one of the DESI Y1 redshift bins (LRG; $0.4<z<0.6$) with an appropriate
analytic covariance. We compare constraints on $\Lambda$CDM parameters with
various prior choices on $n_{s}$, including a uniform prior, Gaussian with
widths of 10$\times$ and 5$\times$ Planck 2018 constraints
($\sigma_{ns}=0.004$)[80], and with $n_{s}$ fixed. These results are shown in
the left panel of Fig. 11 for the Full-Modeling method. We find that for both
the 10$\times$ and 5$\times$ priors on $n_{s}$ the constraints on $\omega_{\rm
cdm}$, $H_{0}$, and $\log A_{s}$ are identical to those when $n_{s}$ is fixed,
suggesting that the Full-Modeling constraints on $\Lambda$CDM parameters are
robust even if the $n_{s}$ constraints from the CMB are systematically off by
10$\sigma$. In order to see how well $n_{s}$ can be constrained completely
independently from Planck we additionally fit to noiseless synthetic mock data
vectors simulating all seven DESI Y1 redshift bins: BGS ($0.1<z<0.4$), LRG
($0.4<z<0.6$), LRG ($0.6<z<0.8$), LRG ($0.8<z<1.1$), ELG ($0.8<z<1.1$), ELG
($1.1<z<1.6$), and QSO ($0.8<z<2.1$) using the appropriate Y1 analytic
covariance for each redshift bin. We compare the case with uniform priors on
$n_{s}$ to the case with $n_{s}$ fixed. These results are shown in the right
panel of Fig. 11. We find that despite the slight degradation in $\Omega_{\rm
m}$ constraint with the flat prior on $n_{s}$, we are able to measure $n_{s}$
to a 3% precision.

###  5.8 Comparison of LPT and EPT

In addition to the LPT model that we primarily focus on in this paper,
velocileptors also has an Eulerian perturbation theory module. The EPT kernels
are constructed from the Lagrangian kernels while setting the IR resummation
scale, $k_{IR}$, to zero. The Eulerian and Lagrangian theories differ in their
treatment of cold dark matter, the first describing dark matter as a perfect
pressureless fluid, and the latter describing it as collisionless particles.
The overdensities derived from both theories agree order-by-order except when
particle trajectories cross. The EPT model in velocileptors employs the galaxy
bias scheme described in Ref. [81]. The mapping between the Lagrangian and
Eulerian bias bases can be achieved within velocileptors via the
transformations [82]:

$\displaystyle b_{1}^{E}$ $\displaystyle=1+b_{1}^{L}$ $\displaystyle
b_{2}^{E}$ $\displaystyle=b_{2}^{L}+\frac{8}{21}b_{1}^{L}\quad,\quad
b_{s}^{E}=b_{s}^{L}-\frac{2}{7}b_{1}^{L}$ $\displaystyle b_{3}^{E}$
$\displaystyle=3b_{3}^{L}+b_{1}^{L}.$ (5.10)

Lastly, the IR resummation in EPT is performed by splitting the wiggle and no-
wiggle parts of the power spectrum, using the same method as is employed in
modeling the poste-reconstruction BAO correlation function (§ 5.4) and
applying a damping factor to the wiggle component. We refer readers to Ref.
[82] for full details of the Eulerian model and how it compares to LPT. We
show in Fig. 12 a comparison of Full-Modeling constraints when fitting the LRG
cubic mocks using LPT and EPT. We see that the constraints agree to within
fractions of a $\sigma$. A more detailed comparison between the two models,
including fits to the ELG and QSO mocks for ShapeFit and Full-Modeling, is
presented in Ref. [47] along with comparisons to other EFT models on the
market.

Figure 12: Full-Modeling fits to the mean of LRG Cubic mocks using the LPT and
EPT models within velocileptors

###  5.9 Varying $f$ and $\sigma_{8}$ separately

The “standard” method of compression involves varying $f$ while keeping
$\sigma_{s8}$ fixed to the fiducial value $\sigma_{8}^{\rm ref}$, and then
reporting the product as $f^{\rm true}\sigma_{s8}^{\rm true}$. In principle,
one should be able to vary $f$ and $\sigma_{s8}$ independently and present the
result as $f^{\rm true}\sigma_{8}^{\rm true}$. This is because the degeneracy
between $f$ and $\sigma_{s8}$ is broken in the 1-loop terms of the power
spectrum. In order to test the ability to constrain $\sigma_{s8}$, we run a
fit in which $\sigma_{s8}(z=0)$ is a free parameter in addition to $f(z)$ and
the other compressed parameters. We vary $\sigma_{s8}(z=0)$ by re-scaling the
linear power spectrum by:

$\displaystyle P_{\rm
lin}^{\prime}(k)=\left(\frac{\sigma_{s8}}{\sigma_{8}^{\rm
fid}}\right)^{2}\times P_{\rm lin}(k),$ (5.11)

Where $\sigma_{8}^{\rm fid}=0.8076$ for the Abacus fiducial cosmology. The
reported $f\sigma_{s8}$ is then $f\sigma_{s8}=f(z)\sigma_{s8}D(z)$ where the
growth factor $D(z)$ is computed from the fiducial value of $\Omega_{\rm
m}=0.315$. We show these results in Fig. 13. We observe that even though
$f\sigma_{s8}$ agrees with that obtained from the standard method, the
$\sigma_{s8}$ constraint of $0.570\pm 0.087$ is significantly below the true
value of $0.8076$. This implies a growth rate $f(z)\sim[\Omega_{\rm
m}(z)]^{0.55}>1$ which is unphysical. While it is unfortunate that the 1-loop
corrections to the power spectrum can not sufficiently constrain $f$ and
$\sigma_{s8}$ independently, we reiterate that our constraint on
$f\sigma_{s8}$ remains robust. We also note a slight degeneracy between
$\sigma_{s8}$ and $m$. While $m$ is designed to change the shape of the power
spectrum, $\sigma_{s8}$ is an integrated quantity that is also mildly affected
by changes in the shape.

Figure 13: ShapeFit constraints, where in the red data we allow
$\sigma_{s8}(z=0)$ to vary independently of $f(z)$. This test is performed
using the covariance for the 25 box volume of 200 $(\,h^{-1}{\rm Gpc})^{3}$.

##  6 Conclusion

Observations are probing the Universe and its evolution with unprecedented
precision, allowing for significant improvements in measurements of
fundamental parameters. The increased constraining power of these data also
increases the sensitivity of our results to systematic effects present in
models and analysis methods. The largest galaxy redshift survey to date, the
Dark Energy Spectroscopic Instrument (DESI), is currently under way with its
first year of Fullshape data being unblinded in the spring of 2024. To prepare
for unblinding we must have a detailed understanding of the sources of
systematic and theoretical error when fitting observations, the flexibility
and limitations of our models, and the performance of different analysis
methods. In this paper we presented tests of these effects using the public
effective-perturbation-theory code velocileptors, fitting data from the the
AbacusSummit suite of simulations. Our focus will be on cosmological
constraints using the Lagrangian Perturbation Theory (LPT) module in
velocileptors, though we also explore fits using its Eulerian Perturbation
Theory (EPT) counterpart. In particular, we fit LRG, ELG, and QSO mock data at
effective redshifts of $z=0.8,1.1,1.4$ respectively, consisting of clustering
measurements from 25 cubic boxes of 8 ($h^{-1}$Gpc)3 each for a total volume
of 200 ($h^{-1}$Gpc)3 for each tracer type. Companion papers to this one,
using other effective perturbation theory codes Folps$\nu$ and PyBird, are
scheduled to appear concurrently (Refs. [44, 45], including in addition a
comparison paper (Ref. [47]) showing that all three effective-theory pipelines
and models behave very similarly when the underlying assumptions and settings
are consistent.

In this paper we discussed three modeling methods: (1) the standard Template
fit, the default method used in previous BOSS and eBOSS analyses, that
compresses observed multipoles into three summary statistics,
($f\sigma_{s8}$,$\alpha_{\parallel}$,$\alpha_{\perp}$) while keeping the
linear power spectrum fixed; (2) the ShapeFit method which introduces an
additional compressed parameter $m$ to the standard Template that modulates
the shape of the linear template power spectrum which depends on early
universe physics; and (3) the Full-Modeling method which directly samples in
the parameter space of a cosmological model in order to fit the data. The
first two methods are model-agnostic and so the compression only needs to be
performed once, after which the obtained summary statistics can be mapped to
any cosmological model ($\Lambda$CDM or extensions) of ones choosing. Despite
the Full-Modeling method technically requiring a Boltzmann code to compute the
linear power spectrum at every step of an MCMC, the use of Taylor series
expansion emulators make the difference in computational cost/time negligible
when compared to the compressed analyses.

We showed throughout the paper that the increased information from the shape
of the linear power spectrum results in significant improvements in
cosmological constraints in ShapeFit when compared to the standard Template
analysis, when CMB data are not included. Compared to the Full-Modeling
approach, ShapeFit provides consistent results on $\Lambda$CDM (and $w$CDM)
parameters with minimal loss in constraining power. In varying the upper bound
of the fitting range, we found that the models give unbiased constraints for
scale cuts up to $k_{\rm max}\leq 0.2\,h{\rm Mpc}^{-1}$. When including priors
from Planck in order to constrain early universe information, all three
methods give consistent results. Since the upcoming data will include tracers
from different redshifts, we tested the ability of our pipelines in fitting
simultaneously the tracers from three redshift bins, finding the joint
analysis to improve the constraints without any noticeable systematic effects.

Because one of the most powerful sources of cosmological information in LSS
that DESI can detect is the Baryon Acoustic Oscillation (BAO) signal, whose
well-defined scale can be used as a standard ruler to constrain the distance-
redshift relation, we combined our fullshape analyses with post-reconstruction
BAO correlation function, finding significant improvements in constraints for
each modeling method. Finally, we also show how each method performs when
extending the parameter space beyond the standard $\Lambda$CDM model by
varying the dark energy equation of state parameter $w$. The ShapeFit and
Full-Modeling methods are both able to obtain consistent and unbiased
constraints within the wCDM model, whereas the standard template suffers
greatly from degeneracies that can not be broken without shape information.

In addition to the velocileptors LPT model, the pipeline also has a module
based on Eulerian perturbation theory (EPT). We show that these two
theoretical frameworks provide consistent constraints, in agreement with the
more extensive comparisons along with other PT pipelines, FOLPS$\nu$ and
PyBird, presented in Ref. [47].

We conclude by summarizing the optimal setup for velocileptors for DESI Y1
fullshape analyses. The scaling of the biases with $\sigma_{8}$ appears to be
a more natural choice of parameterization that is closer to the constraints
from the data and can ameliorate shifts to lower $\sigma_{8}$ in the
posteriors when the data is not sufficiently constraining. We recommend
against the use of the partial Jeffrey’s prior in attempts to reduce
projection effects, due to it being a highly informative prior in the
cosmological parameters. Our counterterm parameterization that scales relative
to linear theory allows for a more intuitive choice of priors on the
$\alpha_{n}$ parameters as “fractional corrections to linear theory”. When
fitting the hexadecapole we strongly suggest restricting the $k-$range in
$P_{4}$ to a $k_{4,\rm max}\sim 0.1\ \,h{\rm Mpc}^{-1}$ as this minimizes the
model’s sensitivity to higher orders in perturbation theory and non-linear
effects such as Fingers of God. For the monopole and quadrupole a scale cut of
$k_{\rm max}=0.20\ \,h{\rm Mpc}^{-1}$ has been found to perform well. Finally,
we also suggest the use of physically motivated Gaussian priors on the
stochastic parameters that can be justified based on the characteristic
physical scales in the system (as captured, for example, in the halo model).

##  7 Data availability

Data from the plots in this paper are available on Zenodo as part of DESI’s
Data Management Plan (DOI: 10.5281/zenodo.10951714). The data used in this
analysis will be made public along the Data Release 1 (details in
https://data.desi.lbl.gov/doc/releases/)

## Acknowledgements

We thank Arnaud de Mattia, Pat McDonald, and other members of the Galaxy and
Quasar Clustering working group within DESI for helpful discussions pertaining
to this work. SC thanks Misha Ivanov and Matias Zaldarriaga for useful
discussions on velocity stochasticities. MM and MW are supported by the DOE.
SC acknowledges the support of the National Science Foundation at the
Institute for Advanced Study. 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 Analytic Marginalization

We can substantially speed up our MCMC fits by analytically marginalizing over
the linear nuisance parameters in our model, i.e. the parameters of the
stochastic and counterterm contributions ($\alpha_{0}$, $\alpha_{2}$,
$\alpha_{4}$, $SN_{0}$, $SN_{2}$, $SN_{4}$). By reducing the number of sampled
parameters our chains are able to converge in under 10 minutes instead of an
hour or two. The procedure for marginalizing over the linear parameters
$b_{i}$ involves splitting the theoretical prediction, into the piece
dependent on the nonlinear parameters $\boldsymbol{a}$ that we sample in and
the “template” piece that is multiplied by the linear parameters:
$\boldsymbol{\Psi}=\boldsymbol{\Psi}_{0}(\boldsymbol{a})+\sum_{i}\theta_{i}\boldsymbol{\Psi}_{t,i}$.
The likelihood distribution marginalized over the linear nuisance parameters
is given by[83, 84]

Figure 14: Comparison of Full-Modeling constraints with and without analytic
marginalization of linear parameters: $\alpha_{0}$, $\alpha_{2}$, SN0, SN2.
Here we show results for the full $25\times(2\,h^{-1}{\rm Gpc})^{3}$ volume to
show that our analytical marginalization method is robust even at very large
volumes.

$\displaystyle
P(\boldsymbol{\Psi}_{d}|\boldsymbol{\Psi}_{0},\boldsymbol{\Psi}_{t},\sigma_{\theta})=\int
d\boldsymbol{\theta}\
\mathcal{L}(\boldsymbol{\Psi}_{d}|\boldsymbol{\Psi}_{0},\boldsymbol{\Psi}_{t},\boldsymbol{\theta})P(\boldsymbol{\theta}),$
(A.1)

where $\boldsymbol{\Psi}_{d}$ is the data and $P(\boldsymbol{\theta})$ denotes
the priors on parameters $\theta_{i}$, which we choose to be Gaussian
(centered at zero) with widths $\sigma_{\theta,i}$:

$\displaystyle
P(\theta_{i}|\sigma_{\theta,i})=\frac{1}{\sqrt{2\pi\sigma_{\theta,i}^{2}}}\exp\left(-\frac{\theta_{i}^{2}}{2\sigma_{\theta,i}^{2}}\right)$
(A.2)

The model likelihood in the integrand is

$\displaystyle\mathcal{L}(\boldsymbol{\Psi}_{d}|\boldsymbol{\Psi}_{0},\boldsymbol{\Psi}_{t},\boldsymbol{\theta})$
$\displaystyle=(2\pi)^{-n/2}\left|\mathcal{C}^{-1}\right|$
$\displaystyle\times
e^{-\frac{1}{2}\left[\boldsymbol{\Psi}_{d}-\left(\boldsymbol{\Psi}_{0}+\sum_{i}\theta_{i}\boldsymbol{\Psi}_{t,i}\right)\right]^{\rm
T}\mathcal{C}^{-1}\left[\boldsymbol{\Psi}_{d}-\left(\boldsymbol{\Psi}_{0}+\sum_{i}\theta_{i}\boldsymbol{\Psi}_{t,i}\right)\right]}.$
(A.3)

Defining $\boldsymbol{\Delta}=\boldsymbol{\Psi}_{d}-\boldsymbol{\Psi}_{0}$ and
$\log\mathcal{L}_{0}=-\frac{1}{2}\boldsymbol{\Delta}^{\rm
T}\mathcal{C}^{-1}\boldsymbol{\Delta}$ we get

$\displaystyle
P(\boldsymbol{\Psi}_{d}|\boldsymbol{\Psi}_{0},\boldsymbol{\Psi}_{t},\sigma_{\theta})$
$\displaystyle\propto\mathcal{L}_{0}\int d\boldsymbol{\theta}\
e^{-\frac{1}{2}\sum_{i,j}\theta_{i}\theta_{j}\left(\boldsymbol{\Psi}_{t,i}^{\rm
T}\mathcal{C}^{-1}\boldsymbol{\Psi}_{t,j}+\frac{1}{\sigma_{i}\sigma_{j}}\delta_{ij}\right)+\sum_{i}\boldsymbol{\Delta}^{\rm
T}\mathcal{C}^{-1}\theta_{i}\boldsymbol{\Psi}_{t,i}}$
$\displaystyle=\mathcal{L}_{0}\int d\boldsymbol{\theta}\
\exp\left[-\frac{1}{2}\left(\boldsymbol{\theta}^{\rm
T}L\boldsymbol{\theta}-V^{\rm T}L^{-1}V\right)\right]$
$\displaystyle=\frac{(2\pi)^{n/2}}{\sqrt{|L|}}\mathcal{L}_{0}e^{\frac{1}{2}V^{\rm
T}L^{-1}V},$ (A.4)

where we completed the square in the second line and defined the matrices
$L_{ij}=\boldsymbol{\Psi}_{t,i}^{\rm
T}\mathcal{C}^{-1}\boldsymbol{\Psi}_{t,j}+\delta_{ij}/(\sigma_{i}\sigma_{j})$
and $V_{i}=\boldsymbol{\Psi}_{t,i}^{\rm T}\mathcal{C}^{-1}\boldsymbol{\Delta}$
before taking the multivariate Gaussian integral. So then the log-likelihood
consists of the four terms

$\displaystyle\log P=\log\mathcal{L}_{0}+\frac{1}{2}V^{\rm
T}L^{-1}V-\frac{1}{2}\log|L|+\frac{n}{2}\log(2\pi).$ (A.5)

Despite analytically marginalizing over the linear parameters, we can always
recover their distribution using the chain containing non-linear parameters.
At each step of the chain, the nonlinear parameters are fixed and the
likelihood is a Gaussian function of the linear parameters with known mean and
variance, i.e. for step $n$ in the MCMC, the likelihood depends on linear
parameter $\theta_{i}$ like:

$\log\mathcal{L}_{n,i}=(\theta_{i}-\bar{\theta}_{i})^{\rm
T}\mathcal{N}^{-1}(\theta_{i}-\bar{\theta}_{i})+\mathrm{const}$ (A.6)

with variance $\mathcal{N}$ and the mean $\bar{\theta}_{i}$ determined by the
(fixed) non-linear parameters. Reconstructing the distribution of parameter
$\theta_{i}$ simply amounts to averaging over all of these Gaussians. This
allows us to still be able to e.g. check the effects of our priors or to
identify any degeneracies between linear parameters and others in the model
that could be driving projection effects. We show in Fig. 14 a comparison of
constraints from the Full Modeling method with and without analytic
marginalization of the linear parameters. For the parameters that are being
sampled in both cases, we find consistent behavior in the contours as
expected. In order to make sure that the analytic marginalization is also
correctly handling the parameters that we marginalize over, we maximize the
first two terms in A.5 (the latter terms describe the volume/width of the
likelihood surface). This gives us the best-fitting values for the nonlinear
parameters. From the maximized posterior, the corresponding best-fit points of
the analytically marginalized parameters can then be directly calculated:

$\displaystyle\theta_{j}^{\rm bf}=\sum_{i}V_{i}L_{ij}^{-1}.$ (A.7)

Once we have found the best-fitting nonlinear parameters and by extension
$\boldsymbol{\Psi}^{\rm bf}=\boldsymbol{\Psi}_{0}^{\rm
bf}+\sum_{i}\theta_{i}^{bf}\boldsymbol{\Psi}_{t,i}^{\rm bf}$, the maximum log-
likelihood is just:

$\displaystyle\log P^{\rm max}=-\frac{1}{2}[\boldsymbol{\Psi}^{\rm bf}]^{\rm
T}\mathcal{C}^{-1}\boldsymbol{\Psi}^{\rm
bf}+\log|\mathcal{C}^{-1}|-\frac{n}{2}\log(2\pi).$ (A.8)

In Table 3 we show the best-fitting parameter values from Full-Modeling fits
with and without analytic marginalization. We see that the parameters that we
marginalize over are well behaved and on the same order as they take when
being sampled.

We also note the third term in Eq. A.5, $-(1/2)\log|L|$, which is the log of
the determinant of the (linear parameter) part of the Fisher matrix. One prior
choice that one can very easily implement is a “partial Jeffrey’s prior” which
removes this term from the likelihood. This prior can cause significant shifts
in constraints in cases where parameter projection effects are noticeable, as
the Jeffrey’s prior removes some of the phase space volume from the
likelihood. We discuss the implications of such a prior in Appendix B.

Non-linear | Params | FM Standard ($\sigma$) | FM Analytic Marg ($\sigma$)
---|---|---|---
$H_{0}$ | 67.67 (0.35) | 67.63 (0.34)
$\Omega_{\rm m}$ | 0.3139 (0.0023) | 0.3143 (${}^{+0.0026}_{-0.0023}$)
$\log(10^{10}A_{\mathrm{s}})$ | 2.998 (${}^{+0.017}_{-0.023}$) | 3.001 (${}^{+0.017}_{-0.026}$)
$b\sigma_{8}$ | 1.642 (${}^{+0.019}_{-0.013}$) | 1.644 (${}^{+0.022}_{-0.013}$)
$b_{2}$ | 0.8982 (${}^{+0.49}_{-0.32}$) | 0.8705 (${}^{+0.51}_{-0.33}$)
| $b_{s}$ | -0.7607 (${}^{+0.55}_{-0.87}$) | -0.8512 (${}^{+0.55}_{-0.95}$)
Linear | $\alpha_{0}$ | 0.6987 (6.1) | 2.468
$\alpha_{2}$ | -11.69 (5.7) | -13.08
$SN_{0}$ | -890.3 (420) | -962.4
$SN_{2}$ | -1.919e4 (4300) | -1.911e4
Table 3: Comparison of Full-Modeling best-fit parameters with and without
analytic marginalization. Uncertainties of the posterior distributions are
given in parentheses for all sampled parameters.

##  Appendix B Parameter projection effects and the role of priors

Figure 15: Toy model examples of information loss in projected posteriors. The
left panel shows the posteriors from sampling from a likelihood distribution
that is constructed out of the sum of a small Gaussian and a Rosenbrock
function in 2D. The dashed lines label the maximum likelihood values of the
two parameters $x_{1}$ and $x_{2}$. The ‘truth’, and likelihood maximum,
appears to be in the tail of the 1D posteriors due to the large volume (area)
at only slightly lower likelihood near $x_{1}\approx x_{2}\approx 0$. The
right panel shows posteriors from two different “data sets” (different
likelihood distributions). Data 1 is constructed from a Rosenbrock function
and Data 2 is a Gaussian distribution. While in the full space (2D) it is
clear the posteriors disagree, in projection (here 1D) the posteriors appear
consistent.

In this section we discuss the role of priors on the parameters of our model
and the effect they can have on parameter projection effects – defined here as
shifts in the marginal posteriors away from the maximum likelihood regions due
to a non-Gaussian posterior surface. These effects frequently arise when there
are several parameters in the model that are poorly constrained or partially
degenerate. If there are degeneracies between parameters in the model, regions
of the parameter space far from the maximum likelihood point may have very
little likelihood penalty compared to the best fit. In spaces with large
numbers of dimensions the “parameter volume” in such regions can be large, and
integration over a subset of these parameters can shift the peaks or means of
the marginal posterior distribution significantly away from the maximum
likelihood values or the “input cosmology” in our tests. In addition, when the
data are not sufficiently powerful the constraints on the cosmological
parameters can depend on the choice of priors and the parameterization.

It is notoriously difficult to visualize complex probability distributions in
high-dimensional spaces, and unfortunately projections _necessarily_ remove
information even if they are given from many viewpoints. For this reason
marginal likelihoods can appear consistent (i.e. overlap in projection) when
they are not and they can appear inconsistent when they are actually
consistent. Even linear changes of the projection axes can change the
appearance of concordance. Such issues are by no means specific to our models:
projection effects in high-dimensional parameter spaces have been encountered
in many areas of cosmology and have been widely discussed in the literature
(see e.g. refs. [85, 86, 87, 88, 89] for recent discussions).

In Fig. 15 we show two toy model examples of projections, where the left plot
is inspired by Fig. 1 of Ref. [87] and the right plot is inspired by Fig. 1 of
Ref. [86]. For the first example, we construct a fake likelihood distribution
by adding a Rosenbrock function,
$f(x_{1},x_{2})=(1.0-x_{1})^{2}+0.5(x_{2}-x_{1}^{2})^{2}$, and a sharp 2D
Gaussian centered at $(\bar{x}_{1}=2.5,\bar{x}_{2}=6)$ with a width of
$\sigma=0.25$ along both parameter directions. The maximum of the total
likelihood distribution is very close to the center of the Gaussian, and is
labeled with grey dashed lines in the figure. However, the contribution of the
Rosenbrock function peaks at $(1.0,1.0)$ but in a much more gradual way. The
result is more likelihood “volume” for the MCMC to explore near $(1.0,1.0)$
than near the true maximum of the whole likelihood. As a result, the marginal
posterior distributions for parameters $x_{1}$ and $x_{2}$ are significantly
offset from the true best-fitting points.

The second cautionary example of projections is presented in the right panel
of Fig. 15 and shows posteriors from two “data sets”, which we simulate by
constructing two different fake likelihood distributions. For Data 1 we again
use a Rosenbrock function,
$f(x_{1},x_{2})=(1.0-x_{1})^{2}+10(x_{2}-x_{1}^{2})^{2}$ and for Data 2 we use
a Gaussian with means $(\bar{x}_{1}=1.5,\bar{x}_{2}=0.0)$ and widths of 0.2.
In this example we demonstrate how the constraints on $x_{1}$ and $x_{2}$
appear to agree for the two data sets when looking at the 1D posteriors, but
in the 2D panel the two data sets are clearly in tension. This serves as a
cautionary tale about interpreting constraints from a multi-dimensional
posterior surface when looking at the projections onto lower dimensions. It is
naturally difficult to visualize an N-dimensional volume, but looking only at
1D or 2D projections of the full distributions might lead one to misinterpret
results.

Finally, as an honorable mention, we refer readers to Fig. 7 of Ref. [85] in
which the authors show a toy model of posteriors from two different data sets
with three sampled parameters, $x$, $y$, $z$. The posteriors for these three
parameters are consistent between data sets. However, after performing a
linear transformation to new coordinates, ($x+y-z$, $x+z-y$, $y+z-x$) one
finds discrepant constraints on $x+y-z$. This shows that tensions can be
hidden due to particular choices of parameterization, and that appropriate
coordinate-independent metrics are necessary to measure the consistency
between data sets or results.

###  B.1 Projection effects for DESI

Figure 16: Full-Modeling $w$CDM fits to synthetic data created in each of the
DESI Y1 redshift bins with corresponding analytic covariances. Left:
Intermediate fits with the real Y1 volumes (black) compared to fits with the
covariances rescaled by a factor of 1/5. Right: Minimal freedom fits with
three different prior choices on the analytically marginalized counter and
stochastic terms. The black contours correspond to a partial Jeffrey’s prior
(only on linear parameters) discussed in the text, the red contours show the
fit with our usual Gaussian priors described in Table 1, and the green
contours correspond to ‘infinite’ priors. The stars in the 2D panels and solid
vertical lines in the diagonal (1D posterior) panels denote the best-fit
models obtained by running a minimizer starting at the MAP values of the
chains.

To demonstrate the impact of projection effects in the specific case of DESI
data with covariances similar to those expected from the first year we turn to
synthetic data created with velocileptors for each of the seven DESI Y1
redshift bins: BGS ($0.1<z<0.4$), LRG ($0.4<z<0.6$, $0.6<z<0.8$, $0.8<z<1.1$),
ELG ($0.8<z<1.1$, $1.1<z<1.6$), and QSO ($0.8<z<2.1$). Since the data we are
fitting to have been generated from the model, with no noise added, the best-
fit point occurs at “truth” and has $\chi^{2}=0$. However $\chi^{2}$ may rise
slowly along some directions which have significant volume, shifting the
marginalized posteriors away from the best-fit point. While the $\Lambda$CDM
(with and without fixing $n_{s}$) and k$\Lambda$CDM models do not exhibit
significant projection effects, we do observe them for wCDM. We show the wCDM
joint fits to the seven Y1 redshift bins in Fig. 16. Note that the marginal
posteriors on several parameters (black lines in the left hand panels of Fig.
16) peak way from the input model, even though the model is, by construction,
a good fit to the (mock) data and the maximum likelihood point is (again by
construction) at the true values of the parameters. As the data become more
constraining these projection effects are reduced – shown as the red contours
in the same figure where the errors have been scaled down by a factor of 5.
Note that some projection effects are still visible in the red contours. The
posterior for $\Omega_{m}$ is still offset by a non-trivial fraction of its
“new” error bar, but the absolute value of the offset is reduced. As we
continue to reduce the error bars the contours shrink to eventually be
$\delta$-functions at the true values. It is also worth noting another feature
of these projection effects. They typically occur when there are many
parameters, some of which are partially degenerate. They also tend to lead to
shifts that are $\mathcal{O}(1\,\sigma)$. This is because the likelihood falls
as $\exp[-\chi^{2}/2]$ moving away from the best-fit point, while the volume
in parameter space grows as a power of the “parameter distance”. Eventually
the Gaussian overcomes the impact of the volume. In the right panel of Fig. 16
we show wCDM constraints to the same synthetic data using three choices of
priors on the linear parameters ($\alpha_{0}$,$\alpha_{2}$,SN0,SN2): infinite
uniform, Gaussian, and the (partial) Jeffrey’s prior. The stars and solid
vertical lines denote the best-fit values obtained from running a minimizer,
and demonstrate that the shift between marginal posteriors and maximum
likelihood values are due to projection effects. We find that these projection
effects are slightly reduced when switching from the flat to Gaussian prior,
showing that the Gaussian priors on the linear parameters are not entirely
uninformative. The projection effects are more significantly reduced when
applying the Jeffrey’s prior and we discuss the implications of using such a
prior in the next section.

###  B.2 Jeffrey’s prior and reparameterizations

In addition to shifts in the posteriors such that they peak away from the
‘true’ values, insufficiently constraining data in a high-dimensional
parameter space can lead to increased sensitivity to priors and choice of
parameterizations. This is another manifestation of the likelihood not
dominating the posterior and is a generic feature of inference in high
dimensions. If we had firm theoretical reasons to prefer one model
parameterization over another this would not be a problem, but in practice
there are several choices between which there is little theoretical
preference. We discuss some of these implications here – first discussing the
choice of parameters and then the Jeffrey’s prior.

A natural141414This is not the only choice. One could imagine choosing e.g.
log priors in the mass scale of the halos hosting the galaxies, or linear
deviations from the peak-background split prediction (where the $b_{n>1}$ are
non-linear functions of $b_{1}$), or many other choices. set of parameters for
the model would be the cosmological parameters (e.g. $\sigma_{8}$) and the
bias parameters and counterterms ($b_{i}$ and $\alpha_{i}$). However some of
these are at least partially degenerate. Lowering $A_{s}$ or $\sigma_{8}$
while raising $\alpha$ can leave $\alpha\,k^{2}P$ unchanged, and a similar
upward adjustment of $b_{i}$ can reduce much of the impact from the other
terms so that $\chi^{2}$ changes little. Since, for linear priors on $b_{i}$
and $\alpha_{i}$, there is more “volume” at large values than small there is a
natural tendency to shift the posterior to lower $\sigma_{8}$. The quantities
best-constrained from observation are the power spectrum multipoles, and in
particular the monopole. For this reason we use parameters that are closer to
the data space, i.e. $b\sigma_{8}$ rather than $b$ (see Table 1). While this
is a natural choice, in terms of the $b_{i}$ it corresponds to a prior that
rises with $\sigma_{8}$ [90]. For example, the Jacobian translating between
$(b,\sigma_{8})$ and $(b\sigma_{8},\sigma_{8})$ is simply $\sigma_{8}$.
Inference using the second set of parameters is thus equivalent to inference
using the first, plus a prior $P(\sigma_{8})\propto\sigma_{8}$. When
$\sigma_{8}$ is not well constrained by the data, this prior choice will shift
the marginal posterior. Similar comments hold for the other parameters of
course.

Figure 17: A 2D slice through the partial Jeffrey’s prior for $\Lambda$CDM. We
show the variation of the prior in the $\Omega_{m}$ and $A_{s}$ directions,
with all other parameters fixed to their best fit. The dashed grey lines in
the lower left panel show the values of $\log(10^{10}A_{s})$ and $\Omega_{m}$
along which $\det F$ was evaluated in each of the other panels. The dashed
grey line in the upper left panel shows a power law $\propto
A_{s}^{3}\propto\sigma_{8}^{6}$. Note the strong dependence of this prior on
the cosmological parameters (see text).

A method that is sometimes used in the statistics literature to reduce the
impact of parameter changes is to include a “Jeffrey’s prior”. This
corresponds to the square root of the determinant of the Fisher matrix, and
has the same role as the familiar $\sqrt{-g}\,d^{4}x$ in General Relativity.
If implemented consistently, this removes the Jacobian from transformations of
variables and so is sometimes termed151515While common, this nomenclature is
incorrect. A much better term would be “reparametersation invariant” since in
general – and in our case – the prior is “informative” from the point of view
of inference. “uninformative”. There are some concerns about taking this
approach in our situation however161616The Jeffrey’s prior and problems with
it are also discussed in ref. [91], including an example from ref. [92]..
First, we do not believe that the physics indicates that e.g.
$(\ln[1+b^{10}],\cosh\sigma_{8})$ is as good a parameter set as
$(b\sigma_{8},\sigma_{8})$ for example. Our parameters have at least some
theoretical justification that we’d like to include as “prior information” in
our model specification. Secondly, as usually implemented, the Jeffrey’s prior
is a strong function of several key cosmological parameters.

To see this, let us consider the partial Jeffrey’s prior that is sometimes
introduced. This involves computing $\sqrt{\det F}$ for only those parameters
that enter the model linearly (if all parameters enter linearly, then this is
the “full” Jeffrey’s prior, however in that limit the likelihood is Gaussian
so the issue of projection effects does not arise). The calculation in the
previous appendix shows that introducing such a prior is equivalent to
dropping the $\log||L||$ term in Eq. (A.5) (see also ref. [88]), making this a
very easy change to make. That this prior is a strong function of the
underlying cosmological parameters is most easily seen by again considering
$\sigma_{8}$. The Fisher matrix has the form

$F\sim\frac{\partial(\mathrm{theory})}{\partial(\mathrm{param})}C^{-1}\frac{\partial(\mathrm{theory})}{\partial(\mathrm{param})}\sim(\mathrm{template})\,C^{-1}\,(\mathrm{template})$
(B.1)

where in the second step we have used the fact that for parameters entering
linearly the derivative is just some linear-parameter-independent template –
e.g. for $\alpha k^{2}P$ it would be $k^{2}P$. In the case of our perturbative
model, each of these ‘templates’ is $P_{\rm lin}$ or some integral over one or
more powers of $P_{\rm lin}$ and thus we expect the template to scale as a
power of $A_{s}$ or $\sigma_{8}$. The Fisher matrix is thus also a (high)
power of $A_{s}$ or $\sigma_{8}$ and so including such a prior has the effect
of shifting the marginal posterior to higher $\sigma_{8}$.

Fig. 17 shows a 2D slice through this (high-dimensional) prior to illustrate
the previous points. We have chosen to show the variation in the $\Omega_{m}$
and $A_{s}$ directions with all of the other parameters held fixed at their
best-fit points. The strong dependence on $A_{s}$ is clear ($\propto
A_{s}^{3}\propto\sigma_{8}^{6}$), and has been described above. The
$\Omega_{m}$ dependence can be understood similarly. Raising $\Omega_{m}$,
with all other parameters fixed, changes the shape of $P_{\rm lin}$ with more
power on the quasi-linear scales of relevance to DESI (and less power at large
scales). The increase in the amplitude of $P_{\rm lin}$ increases $\det F$ in
the same manner as for $A_{s}$ or $\sigma_{8}$. The dependence on each of the
other parameters can be similarly computed and understood, though they are not
shown here for simplicity. The introduction of such a prior is thus
“informative” or “strongly informative” in the sense of introducing non-
negligible shifts in the marginal posteriors given the size of the
uncertainties. We note that in making Fig. 17 we used the more traditional
form for the counterterms, e.g. $\alpha k^{2}P_{\rm lin}$ instead of the
parameterization of Eq. 3.6, since it is in that context that (partial)
Jeffrey’s priors have typically been discussed. For most of this paper we have
chosen parameters scaling like $\alpha\sigma_{8}^{2}$, meaning that the
“template” is closer to $k^{2}P_{\rm lin}/\sigma_{8}^{2}$ and is therefore
largely independent of $\sigma_{8}$. Indeed, we find that in our preferred
parameterization the (partial) Jeffrey’s prior scales much more weakly with
$\sigma_{8}$ than what is usually encountered. However, the strong dependence
on $\Omega_{\rm m}$ and other cosmological parameters is unaffected by this
particular reparameterization.

There are two things to note about these examples. First, in each case the
shift in the marginal posterior was accomplished by the introduction of a what
is effectively a prior, and not by any change in the model or the data. It
relies on the fact that the data are not sufficiently constraining such that
such prior or parameterization choices are relevant. Second, the two
approaches change the prior through different parts of the theory. In the
first case we modified the biases while in the second we introduced a prior
through the counterterms.

Luckily the existing theoretical models are sufficiently accurate to model
much more constraining data than DESI Y1 without the need to introduce
additional free parameters (see the main body of the paper and refs. [38,
47]). As the data become more constraining the impact of parameter choices and
priors is expected to reduce, as shown earlier. Combining the DESI data with
other datasets that can break degeneracies is also expected to reduce the
impact of these effects. In this sense, the Y1 data may well be a “worst case”
scenario.

##  Appendix C Connection to the halo model

It is sometimes helpful to establish the expected sizes of the terms in the
theoretical model. This can be done through arguments of self-consistency (see
main text), and by comparing to other models. In this appendix we compare the
PT approach to a simplified, analytical halo model [93, 94] with the goal of
understanding the expected size of the stochastic terms (see also the
discussion in ref. [66]). Since our goal is to gain insight, we shall deal
with an analytically tractable version of the halo model in which galaxies
reside in spherical, self-similar halos whose centers are distributed
according to biased linear theory with scale-independent bias. If $n(M)$ is
the volume density of halos per unit mass, and each halo has a Fourier-space
density profile $u(k,m,z)$, normalized to unity as $k\to 0$, then the power
spectrum is (see e.g. ref. [95] for a recent, pedagogical discussion with
references to the original literature)

$P_{g}(k,\mu,z)=P_{g}^{\rm 2-halo}(k,\mu,z)+P_{g}^{\rm
1-halo}(k,\mu,z)+P_{g}^{\rm shot}\quad.$ (C.1)

If $N_{\text{cen}}(m)$ and $N_{\text{sat}}(m)$ denote the mean number of
centrals and satellites in a halo of mass $m$ the mean number density of
galaxies is simply $\bar{n}_{g}=\int dm\
n(m)\left[N_{\text{cen}}(m)+N_{\text{sat}}(m)\right]$. To compute the
clustering we need to know the statistics of the galaxy occupation, and we
shall follow standard practice in assuming the centrals are Bernoulli
distributed while the satellites are Poisson distributed.

Under the above assumptions the 2-halo term in the power spectrum is given by:

$P_{g}^{\text{2-halo}}=\left(b_{g}+F\mu^{2}\right)^{2}\;P_{\text{lin}},$ (C.2)

where the bias is

$b_{g}(k,\mu,z)\equiv\frac{1}{\bar{n}_{g}}\int dm\
n(m)\;b(m)\;\left[N_{\text{cen}}+N_{\text{sat}}u(k,m,z)e^{-k^{2}\mu^{2}\sigma_{d}^{2}(m)/2}\right]$
(C.3)

and the effective growth rate of structure is

$F(k,\mu,z)\equiv f\;\int
dm\;n(m)\;\left(\frac{m}{\bar{\rho}}\right)u(k,m)e^{-k^{2}\mu^{2}\sigma_{d}^{2}(m)/2}$
(C.4)

which tends to $f$ as $k\to 0$. In the above we have written the (linear) bias
of a halo of mass $m$ as $b(m)$ and the mean matter density in the Universe as
$\bar{\rho}$. We have also used the fact that in going into redshift-space,
the density profile acquires a damping factor from the virial motions in
halos:

$\displaystyle u(k,m,z)\rightarrow u(k,m,z)e^{-k^{2}\mu^{2}\sigma_{d}^{2}/2},$
(C.5)

where $\sigma_{d}^{2}(m)$ is the velocity dispersion of such a halo in
distance units. The 1-halo term has in its integrand the term $\langle
N(N-1)\rangle$ which, when expanded is:

$\displaystyle\langle N(N-1)\rangle$
$\displaystyle=\langle(N_{\text{cen}}+N_{\text{sat}})(N_{\text{cen}}+N_{\text{sat}}-1)\rangle$
$\displaystyle=\langle
N_{\text{cen}}^{2}-N_{\text{cen}}+2N_{\text{cen}}N_{\text{sat}}+N_{\text{sat}}(N_{\text{sat}}-1)\rangle$
$\displaystyle=2\langle N_{\text{cen}}N_{\text{sat}}\rangle+\langle
N_{\text{sat}}(N_{\text{sat}}-1)\rangle$ $\displaystyle=2\langle
N_{\text{cen}}\rangle\langle N_{\text{sat}}\rangle+\langle
N_{\text{sat}}\rangle$ (C.6)

where in going from the second to third line we used that
$N_{\text{cen}}=0,1\rightarrow\langle N_{\text{cen}}^{2}\rangle=\langle
N_{\text{cen}}\rangle$. We obtain the last equality by assuming that the
centrals and satellites are uncorrelated and that $N_{\text{sat}}$ follows a
Poisson distribution, such that $\langle N_{\text{sat}}^{2}\rangle=\langle
N_{\text{sat}}\rangle^{2}+\langle N_{\text{sat}}\rangle$. Using this, the
1-halo term becomes (dropping the $\langle\rangle$’s for simplicity):

$P_{g}^{\text{1-halo}}=\frac{1}{\bar{n}_{g}^{2}}\int dm\;n(m)\
\left[N_{\text{sat}}^{2}\left|u(k,m,z)\right|^{2}e^{-k^{2}\mu^{2}\sigma_{d}^{2}(m)}+2N_{\text{sat}}N_{\text{cen}}u(k,m,z)e^{-k^{2}\mu^{2}\sigma_{d}^{2}(m)/2}\right].$
(C.7)

Finally, the shot noise power spectrum is simply
$P^{\text{shot}}_{g}=\bar{n}_{g}^{-1}$ if we assume Poisson fluctuations for
the galaxies and halos.

Our perturbative model should be able to describe any ‘complete’ model of
galaxy clustering, whether or not that model is correct in detail. We can make
the connection by considering the low-$k$ limit of the halo model. To make our
expressions slightly simpler we shall make an additional approximation that
$u(k,m,z)\approx 1$ on the scales of interest, which corresponds to assuming
that $kr_{\text{vir}}\ll 1$. We shall further assume that
$\sigma_{d}>r_{\text{vir}}$ so that the impact of virial velocities is more
important than the fact that the satellites do not sit at the halo center.
Under these approximations, and for small $k$,

$\displaystyle b_{g}(k,\mu,z)$ $\displaystyle\simeq\frac{1}{\bar{n}_{g}}\int
dm\
n(m)\;b(m)\left[N_{\text{cen}}+N_{\text{sat}}\left(1-\frac{1}{2}k^{2}\mu^{2}\sigma_{d}^{2}(m)\right)\right]$
(C.8) $\displaystyle=\frac{1}{\bar{n}_{g}}\int dm\
n(m)\;b(m)N_{\text{gal}}-\frac{1}{2}k^{2}\mu^{2}\ \frac{1}{\bar{n}_{g}}\int
dm\ n(m)\;b(m)N_{\text{sat}}\sigma_{d}^{2}(m)$ (C.9) $\displaystyle=b_{\rm
eff}\left(1-\frac{1}{2}k^{2}\mu^{2}\sigma_{\rm 2,eff}^{2}\right)$ (C.10)

The $k^{2}\mu^{2}$ term above, combined with the $b$ or $f\mu^{2}$ term from
the other power of $b_{g}$ in Eq. (C.2) contributes to the counterterms,
$\alpha_{i}$.

Since the mass-integral in $F$ extends all the way to $m=0$, the
$k^{2}\mu^{2}$ correction is smaller than for the bias and we shall neglect
it, taking $F\to f$ henceforth. The 1-halo term becomes

$\displaystyle P_{g}^{\text{1-halo}}$
$\displaystyle\simeq\frac{1}{\bar{n}_{g}^{2}}\int dm\;n(m)\ N_{\rm
sat}e^{-k^{2}\mu^{2}\sigma_{d}^{2}/2}\left[N_{\rm cen}+N_{\rm
sat}e^{-k^{2}\mu^{2}\sigma_{d}^{2}/2}\right]$ (C.11)
$\displaystyle\approx\frac{f_{\rm
sat}}{\bar{n}_{g}}\left(\cdots-\frac{1}{2}k^{2}\mu^{2}\sigma_{\rm
1,eff}^{2}+\cdots\right)$ (C.12)

Thus we see that the halo model predicts that the stochastic terms are of
order $\mathrm{SN}_{0}\sim\bar{n}_{g}^{-1}$ (from $P_{g}^{\rm shot}$ in Eq.
C.1) and $\mathrm{SN}_{2}\sim f_{\rm sat}\sigma_{1,{\rm eff}}^{2}/\bar{n}_{g}$
(from Eq. C.12) as described in the main text. Here $f_{\rm sat}$ is the
satellite fraction such that $f_{\rm sat}\sigma^{2}_{1,\rm eff}$ is the mean
velocity dispersion of halos weighted by $N_{\rm cen}N_{\rm sat}$, such that
roughly speaking $\sigma^{2}_{1,\rm eff}$ is the mean velocity dispersion of
the satellites in question. We often refer to $f_{\rm sat}^{1/2}\sigma_{1,{\rm
eff}}$ as a “characteristic halo velocity” for simplicity.

The simple derivation above neglects several physical effects, including halo
compensation and exclusion, correlations between the halo density and velocity
profiles and between local environment and profile, correlations between mass
bins in the halo shot noise, etc. It is sufficient for order of magnitude
estimates, since most of the neglected effects also have characteristic size
set by the mean inter-galaxy separation or the virial or infall velocity of
the halo but it should not be taken as a ‘complete’ model of clustering. As a
single example of an effect missed by this simple treatment, let us further
consider the effect of virial motions in Eq. C.5. Another way to account for
the effect of FoG in the galaxy power spectrum is to introduce a random
velocity field $\epsilon_{i}(\boldsymbol{s})$ to each galaxy, such that the
observed position is $\boldsymbol{s}+\hat{n}\cdot\epsilon_{\hat{n}}$. In this
case the galaxy 2-point function with these additional velocities is [96, 97]

$\displaystyle P(k,\mu)$ $\displaystyle=\int d^{3}\boldsymbol{s}\
e^{i\boldsymbol{k}\cdot\boldsymbol{s}}\Big{\langle}e^{ik\mu(\epsilon_{\hat{n}}(\boldsymbol{s})-\epsilon_{\hat{n}}(\bf{0}))}\big{(}1+\delta_{g}(\boldsymbol{s})\big{)}\big{(}1+\delta_{g}(\bf{0})\big{)}\Big{\rangle}$
$\displaystyle\approx\int d^{3}\boldsymbol{s}\
e^{i\boldsymbol{k}\cdot\boldsymbol{s}}\Big{\langle}e^{ik\mu(\epsilon_{\hat{n}}(\boldsymbol{s})-\epsilon_{\hat{n}}(\bf{0}))}\Big{\rangle}\Big{(}1+\xi_{g}(\boldsymbol{s})\Big{)}$
(C.13)

where in the second line we have made the (unphysical) assumption that the
virial motions and galaxy densities are uncorrelated in order to isolate the
pure effect of virial velocities usually called FoGs (in the literature models
making this approximation are frequently referred to as “dispersion” or
“streaming” models). The expectation value of the exponential can be expanded
in powers of $k\mu$ as

$\ln\Big{\langle}e^{ik\mu(\epsilon_{\hat{n}}(\boldsymbol{s})-\epsilon_{\hat{n}}(\bf{0}))}\Big{\rangle}=1-k^{2}\mu^{2}\left[\sigma_{v}^{2}-\xi_{\epsilon}(\boldsymbol{s})\right]+\mathcal{O}(k^{3}\mu^{3}).$
(C.14)

where $\xi_{\epsilon}$ is the correlation function of the virial velocities
projected along the line of sight. Since it describes virial motions, this
correlation must fall rapidly to zero outside of the halo radius, $R_{h}$, and
asymptote to the mean square velocity, $\sigma_{v}^{2}$, as $s\rightarrow 0$.
Expanding this cumulant to first order we see that, in addition to the damping
of the profile coming from $-k^{2}\mu^{2}\sigma^{2}_{v}$ in Eq. C.14 we also
gain the contribution

$P(k,\mu)\supset k^{2}\mu^{2}\int d^{3}\boldsymbol{s}\
e^{i\boldsymbol{k}\cdot\boldsymbol{s}}\
\xi_{\epsilon}(\boldsymbol{s})\big{(}1+\xi_{g}(\boldsymbol{s})\big{)}\approx
k^{2}\mu^{2}(1+\sigma^{2}_{g})\int d^{3}\boldsymbol{s}\
e^{i\boldsymbol{k}\cdot\boldsymbol{s}}\ \xi_{\epsilon}(\boldsymbol{s})$ (C.15)

where we have used that the linear galaxy density is smooth compared to the
support of $\xi_{\epsilon}$ and $\sigma^{2}_{g}$ is its the mean on the halo
scale. The integral in the final expression is simply the noise spectrum of
the virial motions, which we expect to be positive and white on large
($>R_{h}$) scales and of order $\sim\sigma_{v}^{2}R_{h}^{3}$. In order to
differentiate between satellites and centrals we can simply set $\epsilon=0$
for central galaxies such that the cumulant in Equation C.14 is instead simply
unity for the central-central correlation and
$1-\frac{1}{2}k^{2}\mu^{2}\sigma^{2}_{v}$ for the central-satellite cross
correlation. This gives the FoG prescription in the ‘analytic halo model’,
derived above, with the addition of a positive, scale-dependent noise along
the line of sight.171717We thank Misha Ivanov for pointing out that the sign
of this effect in N-body simulations is often positive.

We reiterate that our aim here was to motivate the scale of stochastic
contributions and not to make claims about what numerical value (or even sign)
they will take. We see that the term discussed above, while missed by the halo
model, did scale in the same manner as the included terms as we stated above.
Other allowed parameter combinations, such as $R_{h}^{4}\sigma_{v}$ for the
stochastic piece, should be subdominant.

##  Appendix D Further tests

###  D.1 Dependence on $\omega_{\rm b}$ prior

(a) Standard Template

(b) ShapeFit

(c) Full-Modeling

(d) All Methods

Figure 18: Comparison of constraints when loosening the prior on $\omega_{b}$
from $\sigma=0.00037$ to $\sigma=0.001$. In all cases we use the single box
covariance. The bottom right plot shows a comparison of the three modeling
methods while using a $\sigma=0.001$ prior on $\omega_{b}$

We next test the dependence of our constraints on the prior set on
$\omega_{\rm b}$. The standard setting that we choose is a Gaussian prior
centered on $\omega_{\rm b}^{\rm true}=0.02237$ with a width of
$\sigma=0.00037$, which is based on the recentmost Big-Bang Nucleosynthesis
(BBN) constraints on primordial deuterium abundance [67] which places
stringent constraints on $\omega_{\rm b}$. We test the dependence on the prior
by loosening it to $\sigma=0.001$. The results are shown in Fig. 18. Within
each individual method we show results for the covariance appropriate to the
single-box volume. We find that for all three methods, $H_{0}$ becomes
significantly less constrained. Meanwhile the $\Omega_{\rm m}$ constraints
remain unchanged in all methods.

In the Full-Modeling analysis, the measurement of $\Omega_{\rm m}$ is
extracted from the shape of the power spectrum and scale of matter-radiation
equality keq, and these depend on the full matter abundance rather than
$\omega_{b}$ and $\omega_{cdm}$ separately. We thus do not see a degradation
in the $\Omega_{\rm m}$ constraint when the prior on $\omega_{b}$ is relaxed.
In the template and ShapeFit analyses $\Omega_{\rm m}$ is inferred from the
compressed parameters, and because $f\simeq\Omega_{\rm m}^{0.55}$ we can
extract a measurement of $\Omega_{\rm m}$ from the compressed amplitude
parameter without any dependence on $\omega_{b}$ prior. In the ShapeFit case,
additional constraining power on $\Omega_{\rm m}$ comes from the shape
parameter $m$, but just like in the Full-Modeling case this power spectrum
shape information translates to a measurement $\Omega_{\rm m}$ without any
reliance on $\omega_{b}$ specifically.

For the $H_{0}$ measurement we do observe a significant degradation in
constraining power when the prior on $\omega_{b}$ is relaxed. In the template
analysis, information about cosmological distances is extracted from the BAO
feature and thus constrains $H(z)r_{d}$ and $D_{A}(z)/r_{d}$. Breaking the
degeneracy between $H_{0}$ and $r_{d}$ requires a physical (dimensionful)
length scale for the distance-redshift relation beyond just the angular size
of the BAO feature [98]. This is accomplished with knowledge about
$\omega_{b}$ (which determines $r_{d}$) from either BBN or CMB and then leads
to a direct measurement of $H_{0}$. Therefore, relaxing the prior on
$\omega_{b}$ worsens the constraint on $H_{0}$. The inclusion of the shape
parameter $m$, while in general improving constraints when compared to the
standard template, does not compensate for the changes in $\omega_{b}$
information and therefore ShapeFit also experiences worse $H_{0}$ constraint.
The Full-Modeling method can in principle constrain $\omega_{b}$ (and by
extension $r_{d}$) in the absence of an external prior because the amplitude
of BAO wiggles depend on $\omega_{b}$ and $\omega_{\rm cdm}$ and can be
modulated in Full-Modeling analyses, but this is still a much weaker
constraint than what can be accomplished with a BBN prior [99].