# Testing one-loop galaxy bias: cosmological constraints from the power
spectrum

Andrea Pezzotta<EMAIL_ADDRESS>Max-Planck-Institut für
extraterrestrische Physik, Postfach 1312, Giessenbachstr., 85741 Garching,
Germany  Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can
Magrans, s/n, 08193 Barcelona, Spain  Institut d’Estudis Espacials de
Catalunya (IEEC), 08034 Barcelona, Spain  Martin Crocce  Institute of Space
Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona,
Spain  Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona,
Spain  Alexander Eggemeier  Institute for Computational Cosmology, Department
of Physics, Durham University, South Road, Durham DH1 3LE, United Kingdom
Ariel G. Sánchez  Max-Planck-Institut für extraterrestrische Physik, Postfach
1312, Giessenbachstr., 85741 Garching, Germany  Román Scoccimarro  Center for
Cosmology and Particle Physics, Department of Physics, New York University, NY
10003, New York, USA

###### Abstract

We investigate the impact of different assumptions in the modeling of one-loop
galaxy bias on the recovery of cosmological parameters, as a follow up of the
analysis done in the first paper of the series at fixed cosmology. We use
three different synthetic galaxy samples whose clustering properties match the
ones of the CMASS and LOWZ catalogues of BOSS and the SDSS Main Galaxy Sample.
We investigate the relevance of allowing for either short range non-locality
or scale-dependent stochasticity by fitting the real-space galaxy auto power
spectrum or the combination of galaxy-galaxy and galaxy-matter power spectrum.
From a comparison among the goodness-of-fit ($\chi^{2}$), unbiasedness of
cosmological parameters (FoB), and figure-of-merit (FoM), we find that a four-
parameter model (linear, quadratic, cubic non-local bias, and constant shot-
noise) with fixed quadratic tidal bias provides a robust modelling choice for
the auto power spectrum of the three samples, up to $k_{\rm
max}=0.3\,h\,\mathrm{Mpc}^{-1}$ and for an effective volume of
$6\,h^{-3}\,\mathrm{Gpc}^{3}$. Instead, a joint analysis of the two
observables fails at larger scales, and a model extension with either higher
derivatives or scale-dependent shot-noise is necessary to reach a similar
$k_{\rm max}$, with the latter providing the most stable results. These
findings are obtained with three, either hybrid or perturbative, prescriptions
for the matter power spectrum, RESPRESSO, gRPT and EFT. In all cases, the
inclusion of scale-dependent shot-noise increases the range of validity of the
model in terms of FoB and $\chi^{2}$. Interestingly, these model extensions
with additional free parameters do not necessarily lead to an increase in the
maximally achievable FoM for the cosmological parameters
$\left(h,\,\Omega_{c}h^{2},\,A_{s}\right)$, which are generally consistent to
those of the simpler model at smaller $k_{\rm max}$.

###### pacs:

Valid PACS appear here

##  I Introduction

Over the past decades galaxy redshift surveys have provided a wealth of
information on the large-scale distribution of galaxies across the Universe.
Clustering measurements of two-point statistics – the galaxy power spectrum
$P(k)$ and the two-point correlation function $\xi(s)$ – from large data
samples can indeed provide precise measurements about the underlying
cosmological model Davis and Peebles (1983); Maddox  _et al._ (1990); Tegmark
_et al._ (2004); Cole  _et al._ (2005).

The inference of cosmological parameters from large-scale structure is made
intrinsically more difficult by the realisation that galaxies are a biased
tracer of the total matter density field Hauser and Peebles (1973); Kaiser
(1984); Davis  _et al._ (1985); Bardeen  _et al._ (1986); Tegmark  _et al._
(2004); Zehavi  _et al._ (2005); Desjacques  _et al._ (2018). This translates
into the necessity of having an accurate description of the relationship
between the galaxy and matter density fields, a phenomenon which is commonly
referred to as galaxy bias. Even though the latter can be partially understood
in a phenomenological way, e.g. using results from N-body simulations, the
complexity of dealing with tracers featuring different morphological
properties makes desirable to develop an analytic formulation that is based on
a more theoretical background. From this point of view, perturbative
approaches stand as a natural way of describing galaxy bias in a physically
motivated way.

The main idea behind this formulation is that the galaxy overdensity
$\delta_{\rm g}$ can be expressed in terms of a series of operators involving
spatial derivatives of the gravitational and velocity potentials. At leading
order, this relation is captured by a single multiplicative factor, i.e.
$\delta_{\rm g}=b_{1}\,\delta$, where $\delta$ is the matter density contrast
and $b_{1}$ is a multiplicative factor called linear bias Kaiser (1984).
Higher-order contributions become progressively more important on mildly non-
linear scales, as expected from a spherically-symmetric gravitational collapse
Mo and White (1996); Mo  _et al._ (1997), in a way that the expression for
$\delta_{\rm g}$ can be expanded to include higher powers of $\delta$, i.e.
$\delta_{\rm g}=\sum_{n}\,b_{n}/n!\,\delta^{n}$ Kaiser (1984); Frieman and
Gaztanaga (1994). At the same time it has been shown that anisotropies in the
process of gravitational collapse are responsible for the generation of non-
negligible tidal effects, which also contribute to the local distribution of
galaxies Chan  _et al._ (2012); Baldauf  _et al._ (2012). This finding
followed the realization that the local-in-matter-density bias model had
limitations in providing a proper description of the clustering of dark matter
halos Manera and Gaztañaga (2011); Roth and Porciani (2011) and was leading to
incompatible constraints on the quadratic bias $b_{2}$ from measurements of
the power spectrum and bispectrum Pollack  _et al._ (2012, 2014).

The most important aspect of the previously described model is that galaxy
bias is treated as a spatially local quantity. However, it is well known that
the formation of halos and galaxies is triggered by the gravitational collapse
of matter from a spatially finite region, and therefore the local assumption
is bound to fail when approaching scales that roughly correspond to the
Lagrangian radius $R$ of the host halos. In terms of $\delta_{\rm g}$, this
effect can be taken into account by considering not only its dependency on
density and tidal fields, but also on functionals of $\delta$ Desjacques  _et
al._ (2018). This is equivalent to introducing higher derivatives of the
matter density field, that at leading order provide a contribution to
$\delta_{g}$ of the form $R^{2}\,\nabla^{2}\delta$ McDonald and Roy (2009);
Desjacques (2008); Desjacques  _et al._ (2010). A further ingredient to the
galaxy-matter relation is represented by stochastic terms, that on large
scales behave as an additional contribution to Poissonian shot-noise Scherrer
and Weinberg (1998). Stochasticity is the direct result of small-scale
perturbations, which are not correlated over long distances under the
assumption of Gaussian initial conditions Dekel and Lahav (1999); Taruya and
Soda (1999); Matsubara (1999). At higher wave modes the halo-exclusion effect
Mo and White (1996); Sheth and Lemson (1999); Smith  _et al._ (2007) imprints
a scale dependence on the stochastic contributions, whose strength is
controlled by the Lagrangian radius Baldauf  _et al._ (2013), similarly to
higher derivatives. At next-to-leading order, this contribution scales as
$k^{2}$, and might become relevant even for clustering analyses based on two-
point statistics, as shown in Eggemeier  _et al._ (2020).

Galaxy bias models based on the one-loop perturbative expansion have been used
to extract cosmological constraints from big data collaborations, such as
BOSS, using clustering measurements both in configuration Sánchez  _et al._
(2017) and in Fourier space Gil-Marín  _et al._ (2017); Beutler  _et al._
(2017); Grieb  _et al._ (2017). More recently, the same data have been re-
analysed in a number of works Ivanov  _et al._ (2019); D’Amico  _et al._
(2019); Tröster  _et al._ (2020) with novel techniques that nevertheless
assume the same biasing scheme for the galaxy-matter relationship. The
majority of these analyses assume a different ansatz in terms of the degrees
of freedom of the model, by fixing some of the free bias parameters to some
physically motivated relations. For example Sánchez  _et al._ (2017)Grieb  _et
al._ (2017) and Tröster  _et al._ (2020) fix the cubic non-local bias
parameter to the Local Lagrangian (LL) relation, while Beutler  _et al._
(2017) also fix the quadratic non-local parameter. On the other side, Ivanov
_et al._ (2019) set the cubic term to zero and leave the quadratic tidal bias
completely free.

This paper is part of a set of works Eggemeier  _et al._ (2020)Eggemeier  _et
al._ (2021) in which we explore the impact of different bias modeling choices
using a set of three performance metrics, namely the goodness-of-fit, the
unbiasedness of sampling parameters, and the merit of the model. To do that,
we make use of a set of three different simulated galaxy samples with an
effective volume of $6\,h^{-3}\,\mathrm{Gpc}^{3}$, whose clustering properties
and number densities reproduce the ones of three real catalogues, the CMASS
and LOWZ samples of BOSS Eisenstein  _et al._ (2011); Dawson  _et al._ (2013);
Reid  _et al._ (2015), and the Main Galaxy Sample of SDSS Strauss  _et al._
(2002). In order to exclusively concentrate on the modeling of one-loop galaxy
bias, we carry out this analysis in real-space removing the impact of
redshift-space distortions, which would require an additional modeling layer.
In Eggemeier  _et al._ (2020) we performed a fixed cosmology analysis by using
the measured non-linear matter power spectrum. In that case, we used the
linear bias parameter as a proxy for the goodness of our model. Here, we
additionally model the impact of matter non-linear evolution and sample over
cosmological parameters to determine the impact of one-loop galaxy bias on the
recovery of such parameters. In the first part of the paper we use the hybrid
perturbative-simulated approach, RESPRESSO, while comparing results for
different PT-based dark matter models in a later section.

Our paper is organised as follows. In Section II we summarise the main
ingredients of our galaxy bias model including a description of all the terms
contributing at one-loop in perturbation theory. This includes also a short
review of the three non-linear matter predictions we employ in this work. In
Section III we describe the simulated galaxy samples we use to test one-loop
galaxy bias, along with fitting procedure, parameter priors and the
performance metrics we introduced above. Main results from fit of the auto
galaxy power spectrum and the combination between the latter and the galaxy-
matter cross power spectrum are given in Section IV. We finally draw our
conclusions in Section V.

##  II One-loop perturbation theory for biased tracers

The theory describing the evolution of the clustering of biased tracers on
mildly non-linear scales is a well established topic (for a review of galaxy
bias see Desjacques  _et al._ (2018)) that can be naturally described in the
framework of perturbation theory. This section stands as an overview of the
most important results obtained at one-loop in standard perturbation theory
(SPT) and derived approaches. Notice that for sake of readability we omit all
the dependences on redshift $z$.

###  II.1 Galaxy bias expansion

The general perturbative expansion of galaxy bias can be interpreted as a sum
of different operators that are functions of the gravitational potential
$\Phi$ and velocity potential $\Phi_{\rm v}$. Focusing on one-loop
contributions, the relationship between biased tracers and the underlying dark
matter density field can be described considering terms up to third order in
the matter perturbation $\delta$. Following the notation of Eggemeier  _et
al._ (2019) we can write

$\begin{split}\delta_{\rm
g}(\bm{x})=\;\,&\overline{b}_{1}\,\delta(\bm{x})+\frac{\overline{b}_{2}}{2}\,\delta^{2}(\bm{x})+\overline{\gamma}_{2}\,\mathcal{G}_{2}(\Phi_{\rm
v}|\,\bm{x})\\\
&+\overline{\gamma}_{21}\,\mathcal{G}_{2}(\varphi_{1},\varphi_{2}|\,\bm{x})+\overline{\beta}_{1}\,\nabla^{2}\delta(\bm{x})+\ldots\,,\end{split}$
(1)

where the first two terms on the right-hand side are part of the standard
local expansion in powers of $\delta$.

In the previous equation $\mathcal{G}_{2}$ is a Galilean invariant operator,
representing the tidal stress tensor generated by the velocity potential
$\Phi_{\rm v}$ , and it is given by

$\mathcal{G}_{2}(\Phi_{\rm v})=(\nabla_{ij}\Phi_{\rm
v})^{2}-(\nabla^{2}\Phi_{\rm v})^{2}.$ (2)

In Fourier space 111We adopt the usual convention for the Fourier transform,
$\delta(\bm{x})=\int_{\bm{k}}\exp{(-i\bm{k}\cdot\bm{x})}\,\delta(\bm{k})$, and
use the short-hand notation for the 3d integrals,
$\int_{\bm{k}_{1},\ldots,\bm{k}_{n}}\equiv\int\text{d}^{3}k_{1}/(2\pi)^{3}\cdots\text{d}^{3}k_{n}/(2\pi)^{3}$.
this translates into

$\mathcal{G}_{2}(\bm{k})=\\!\\!\int_{\bm{q}}\\!\left[\frac{\bm{q}\cdot(\bm{k}-\bm{q})}{q^{2}\,|\bm{k}-\bm{q}|^{2}}-1\right]\theta(\bm{q})\,\theta(\bm{k}-\bm{q})\,,$
(3)

where $\theta$ is the divergence of the matter velocity field, such that
$\nabla^{2}\Phi_{\rm v}\equiv\theta$.

Differently from the first two terms of Equation (1), terms involving
$\mathcal{G}_{2}$ incorporate non-local (in matter density $\delta$)
contributions that spontaneously arise from the non-linear evolution of the
matter density field. The second of these terms is obtained by expressing the
non-linear velocity potential up to second order (i.e. $\Phi_{\rm v}=\Phi_{\rm
v}^{(1)}+\Phi_{\rm v}^{(2)},\,\varphi_{1}=-\Phi_{\rm
v}^{(1)},\,\varphi_{2}=-\Phi_{\rm v}^{(2)}$), and keeping the next-to-leading
order correction Chan  _et al._ (2017), leading to

$\mathcal{G}_{2}(\varphi_{1},\varphi_{2})=\nabla_{ij}\varphi_{1}\nabla_{ij}\varphi_{2}-\nabla^{2}\varphi_{1}\nabla^{2}\varphi_{2},$
(4)

where $\nabla^{2}\varphi_{1}=-\theta$ is the linear velocity divergence field,
and $\nabla^{2}\varphi_{2}=-\mathcal{G}_{2}(\varphi_{1})$ is the next-to-
leading order. The net result of adding this higher-order correction is that
$\delta_{\rm g}$ collects contributions up to third order in the matter
perturbation $\delta$.

In addition to the standard expansion, any biased tracer also comes with a
physical scale which regulates the importance of higher-derivative operators
Bardeen  _et al._ (1986); Matsubara (1999); McDonald and Roy (2009);
Desjacques  _et al._ (2010); Musso and Sheth (2012), whose leading order
scales as $\nabla^{2}\delta$. This scale quantifies the size of the region in
which galaxy formation occurs, and it is therefore influenced by any short-
range gravitational effect and baryonic corrections. For halos, this scale is
close to the Lagrangian radius McDonald and Roy (2009); Lazeyras and Schmidt
(2019), while it differs for other kinds of tracers, such as galaxies and
quasars, depending on their type.

An important aspect of the previously described biasing scheme is the
renormalization of the parameters on which it is based. As a matter of fact,
the bare bias parameters listed in Equation (1) are sensitive to the UV cutoff
scale used to make one-loop integrals convergent. In particular, the current
basis carries a dependence on the variance of the matter density field
$\sigma^{2}=\mathcal{h}\delta^{2}(\bm{x})\mathcal{i}$, and more generally to
any higher-order correlator $\mathcal{h}\delta^{(n)}\mathcal{i}$. This
dependence is completely non-physical and can be reabsorbed by means of an
adequate renormalization of the bias parameters McDonald (2006); Assassi  _et
al._ (2014). The new basis (for convention this is denoted without the
overscript hat) can be identified using the peak-background split formalism
Schmidt  _et al._ (2013), and it is constructed in a way that its components
quantify the response of the cosmic mean abundance of tracers to a change in
the background density, with no dependence on the one-loop cutoff scale.

An alternative approach to renormalization can be obtained expanding
$\delta_{\rm g}$ in terms of the galaxy multi-point propagators, a topic which
was firstly described in Eggemeier  _et al._ (2019). Consistently with the
naming convention adopted in the context of renormalized perturbation theory
Crocce and Scoccimarro (2006a), these quantities are defined as the ensemble
averaged derivatives of $\delta_{\rm g}$ with respect to $\delta_{\rm L}$,

$\begin{split}\left<\frac{\partial\delta_{\rm g}(\bm{k})}{\partial\delta_{\rm
L}(\bm{k}_{1})\,\cdots\,\partial\delta_{\rm
L}(\bm{k}_{n})}\right>\equiv\;\,&(2\pi)^{3}\,\Gamma_{\rm
g}^{(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})\\\ &\times\,\delta_{\rm
D}(\bm{k}-\bm{k}_{1\cdots n})\,,\end{split}$ (5)

and, being observables themselves, they are already normalized by
construction. For this reason, the multi-point propagators can be identified
with the scale-dependent bias parameters in a way that

$\delta_{\rm g}=\Gamma_{\rm g}^{(1)}\otimes{\cal H}_{1}+\Gamma_{\rm
g}^{(2)}\otimes{\cal H}_{2}+\ldots\,,$ (6)

where the ${\cal H}_{n}$ are the Wiener-Hermite functionals Eggemeier  _et
al._ (2019); Bernardeau  _et al._ (2002), and the operator $\otimes$ is
defined as

$\begin{split}\left[\Gamma_{\rm g}^{(n)}\otimes{\cal
H}_{n}\right](\bm{k})\equiv\;\,&(2\pi)^{3}\int_{\bm{k}_{1},\ldots,\bm{k}_{n}}\delta_{\rm
D}(\bm{k}-\bm{k}_{1\cdots n})\\\ &\times\,\Gamma_{\rm
g}^{(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})\,{\cal
H}_{n}(\bm{k}_{1},\ldots,\bm{k}_{n})\,,\end{split}$ (7)

Using this approach, we can define our observables in a fully consistent
framework. In this paper we are interested in the galaxy auto power spectrum
$P_{\rm gg}$ and galaxy-matter cross power spectrum $P_{\rm gm}$, which are
defined as the auto and cross-correlation of the two density fields
$\delta_{\rm g}$ and $\delta_{\rm m}$, as

$\displaystyle\mathcal{h}\delta_{\rm g}(\bm{k})\,\delta_{\rm
g}(\bm{k}^{\prime})\mathcal{i}$ $\displaystyle\equiv(2\pi)^{3}P_{\rm
gg}(k)\,\delta_{\rm D}(\bm{k}+\bm{k}^{\prime})\,,$ (8)
$\displaystyle\mathcal{h}\delta_{\rm g}(\bm{k})\,\delta_{\rm
m}(\bm{k}^{\prime})\mathcal{i}$ $\displaystyle\equiv(2\pi)^{3}P_{\rm
gm}(k)\,\delta_{\rm D}(\bm{k}+\bm{k}^{\prime})\,.$ (9)

Substituting Equation (1) in the previous set, we get to the full expressions
for the galaxy auto and cross power spectra at one-loop,

$\displaystyle\begin{split}P_{\rm gg}(k)=&\;\,b_{1}^{2}\,P_{\rm
mm}(k)+b_{1}b_{2}\,P_{b_{1}b_{2}}(k)+b_{1}\gamma_{2}\,P_{b_{1}\gamma_{2}}(k)\\\
&+b_{2}^{2}\,P_{b_{2}b_{2}}(k)+b_{2}\gamma_{2}\,P_{b_{2}\gamma_{2}}(k)+\gamma_{2}^{2}\,P_{\gamma_{2}\gamma_{2}}(k)\\\
&+b_{1}\gamma_{21}\,P_{b_{1}\gamma_{21}}(k)-2b_{1}\beta_{1}\,k^{2}P_{\rm
L}(k)\,,\end{split}$ (10) $\displaystyle\begin{split}P_{\rm
gm}(k)=&\;\,b_{1}\,P_{\rm
mm}(k)+\frac{1}{2}\bigg{[}b_{2}\,P_{b_{1}b_{2}}(k)+\gamma_{2}\,P_{b_{1}\gamma_{2}}(k)\\\
&+\gamma_{21}\,P_{b_{1}\gamma_{21}}(k)\bigg{]}-\beta_{1}\,k^{2}P_{\rm
L}(k)\,,\end{split}$ (11)

where $P_{\rm L}$ is the linear matter power spectrum, $P_{\rm mm}$ is the
non-linear matter power spectrum, and the one-loop bias corrections read

$\displaystyle\begin{split}P_{b_{1}b_{2}}(k)=&\;\,2\\!\\!\int_{\bm{q}}\\!F_{2}(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(|\bm{k}-\bm{q}|)P_{\rm L}(q)\,,\end{split}$ (12)
$\displaystyle\begin{split}P_{b_{1}\gamma_{2}}(k)=&\;\,P_{b_{1}\gamma_{2}}^{\rm
mc}(k)+P_{b_{1}\gamma_{2}}^{\rm prop}(k)\\\
=&\;\,4\\!\\!\int_{\bm{q}}\\!F_{2}(\bm{k}-\bm{q},\bm{q})S(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(|\bm{k}-\bm{q}|)P_{\rm L}(q)\\\ &+8P_{\rm
L}(k)\\!\\!\int_{\bm{q}}\\!G_{2}(\bm{k},\bm{q})S(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(q)\,,\end{split}$ (13)
$\displaystyle\begin{split}P_{b_{2}b_{2}}(k)=&\;\,\frac{1}{2}\int_{\bm{q}}\\!\left[P_{\rm
L}(|\bm{k}-\bm{q}|)P_{\rm L}(q)-P_{\rm L}^{2}(q)\right]\,,\end{split}$ (14)
$\displaystyle\begin{split}P_{b_{2}\gamma_{2}}(k)=&\;\,2\\!\\!\int_{\bm{q}}\\!S(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(|\bm{k}-\bm{q}|)P_{\rm L}(q)\,,\end{split}$ (15)
$\displaystyle\begin{split}P_{\gamma_{2}\gamma_{2}}(k)=&\;\,2\\!\\!\int_{\bm{q}}\\!S^{2}(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(|\bm{k}-\bm{q}|)P_{\rm L}(q)\,,\end{split}$ (16)
$\displaystyle\begin{split}P_{b_{1}\gamma_{21}}(k)=&\;\,4P_{\rm
L}(k)\\!\\!\int_{\bm{q}}\\!S(\bm{k}-\bm{q},\bm{q})S(\bm{k},\bm{q})P_{\rm
L}(q)\,.\end{split}$ (17)

Here $F_{2}$ and $G_{2}$ are the symmetrised second-order mode-coupling
kernels Bernardeau  _et al._ (2002), and
$S(\bm{k}_{1},\bm{k}_{2})=\left(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}\right)^{2}-1$
is the Fourier transform of the kernel describing $\mathcal{G}_{2}(\Phi_{\rm
v})$, as shown in Equation (3).

The only two propagator-like contributions are perfectly degenerate with each
other, and follow the relation

$P_{b_{1}\gamma_{21}}(k)=-\frac{7}{6}P_{b_{1}\gamma_{2}}^{\rm prop}(k)\,.$
(18)

For this reason, in a real analysis, it is common practice to either neglect
one of the two tidal field related parameters (e.g. Ivanov  _et al._ (2019))
or to assume perfectly local-in-matter-density relations (see Section II.4) to
express one of them in terms of lower order local bias parameters (e.g.
Sánchez  _et al._ (2017); Grieb  _et al._ (2016)).

Since the $P_{b_{2}b_{2}}$ contribution does not asymptote to 0 in the large-
scale limit, we renormalise it as in Equation (14), and absorb the constant
low-$k$ amplitude as an additional contribution to the shot noise error, that
will be discussed in Section II.3.

###  II.2 Matter modeling

In this section we discuss the modeling options for the non-linear matter
power spectrum $P_{\rm mm}$ in Equations (10) and (11). As a matter of fact an
accurate modeling of $P_{\rm mm}$ is essential, as any systematic effect in
the description of the matter density field in the range of scales we are
considering might lead to invalid interpretations of the galaxy-matter bias
relationship.

Differently from Eggemeier  _et al._ (2020), where one-loop bias was
investigated at fixed cosmology and adopting the measured matter power
spectrum as reference, here the main goal is to assess the level of accuracy
of our model in terms of cosmological parameters. For this reason, we have to
explicitly assume a model for $P_{\rm mm}$. In the rest of this section we
provide a description of the three models we are going to test. In particular
we use 1) a refined RPT-derived model, based on the preservation of Galilean
invariance (dubbed gRPT), 2) an EFT-like approach based on BAO damping and a
non-trivial stress tensor, and 3) a mixed approach, RESPRESSO, based on
accurate N-body simulations and a perturbative expansion around the fiducial
cosmology.

####  II.2.1 Standard perturbation theory

The basic assumption of SPT is that dark matter behaves as a perfect
pressureless fluid on large enough scales, where matter is not subject to
shell crossing as in multi-streaming regions. Under this assumption, and after
having expanded the matter density contrast in a Taylor series, i.e.
$\delta=\delta^{(1)}+\delta^{(2)}+\delta^{(3)}+\ldots$, we find solutions at
every order in perturbation as Bernardeau  _et al._ (2002)

$\begin{split}\delta^{(n)}(\bm{k})=\int_{\bm{q}_{1}\ldots\bm{q}_{n}}&\delta_{\rm
D}(\bm{k}-\bm{q}_{1}-\ldots-\bm{q}_{n})\,F_{n}(\bm{q}_{1},\ldots,\bm{q}_{n})\\\
&\times\delta_{\rm L}(\bm{q}_{1})\ldots\delta_{\rm
L}(\bm{q}_{n})\,,\end{split}$ (19)

where $F_{n}$ is the $n$-th order symmetrized kernel describing the non-linear
mode coupling between fluctuations at different wavelengths.

Moving to two-point statistics, the expansion for the matter power spectrum
can be written as

$P_{\rm mm}(k)=P_{\rm L}(k)+P^{1\mbox{-}\mathrm{loop}}(k)+\ldots\,,$ (20)

where at one-loop the only non-vanishing contributions are

$\begin{split}P^{1\mbox{-}\mathrm{loop}}(k)=&\>\,P_{22}(k)+P_{13}(k)\\\
=&\>\,2\\!\\!\int_{\bm{q}}\\!F_{2}(\bm{k}-\bm{q},\bm{q})P_{\rm
L}(\bm{k}-\bm{q})P_{\rm L}(q)\\\ &+6\,P_{\rm
L}(k)\\!\\!\int_{\bm{q}}\\!F_{3}(\bm{q},-\bm{q},\bm{k})P_{\rm
L}(q)\,.\end{split}$ (21)

It is now well established that a SPT approach like the one described above
lead to significant residuals if compared to the output of numerical
simulations, even including higher-order corrections Blas  _et al._ (2013).
The source of this inaccuracy can be mostly identified in two separate
effects, whose description is the subject of the next two sections.

####  II.2.2 BAO damping from large-scale “infrared” modes

One of the most acknowledged deviations between one-loop SPT predictions and
the matter power spectrum measured from numerical simulations is the shape of
the BAO features. Since the characteristic scale of the BAO peak is much
larger than the scale at which non-linear contributions become important, one
may expect a standard perturbative approach to provide accurate predictions on
that scale. However, large-scale bulk motions produce a non-negligible effect
on the amplitude of the power spectrum at the BAO scale, the most significant
of which is a smearing of the BAO signal due to the large-scale relative
displacement field Crocce and Scoccimarro (2006b); Eisenstein  _et al._
(2007a, b); Matarrese and Pietroni (2008); Crocce and Scoccimarro (2008).

These corrections to the matter power spectrum were firstly resummed in the
context of RPT Crocce and Scoccimarro (2006b), and at leading order the net
effect is to apply a damping factor to the BAO wiggles. Practically, we can
express the matter power spectrum as the sum of a smooth ($P_{\rm nw}$) and
wiggly ($P_{\rm w}$) term Seo  _et al._ (2008),

$P(k)=P_{\rm nw}(k)+P_{\rm w}(k).$ (22)

The smooth-wiggle split can be realised adopting several different recipes,
and throughout this paper we follow the approach described in Vlah  _et al._
(2016) and Osato  _et al._ (2019), where the smooth component is defined as a
rescaling of the featureless spectrum firstly defined in Eisenstein and Hu
(1999) to account for broadband difference with the linear power spectrum.

At leading order, the damping factor can be calculated assuming the Zel’dovich
approximation, and subsequently applied to the wiggly component, so that

$P_{\rm LO}(k)=P_{\rm nw}(k)+e^{-k^{2}\Sigma^{2}}P_{\rm w}(k),$ (23)

where

$\Sigma^{2}=\\!\int_{0}^{k_{\rm S}}\\!P_{\rm
nw}(q)\left[1-j_{0}\left(\frac{q}{k_{\mathrm{BAO}}}\right)+2j_{2}\left(\frac{q}{k_{\mathrm{BAO}}}\right)\right]\frac{dq}{6\pi^{2}}\,,$
(24)

is the relative displacement field two-point function at the BAO scale
Eisenstein  _et al._ (2007a). Here $j_{n}$ is the $n$-th order spherical
Bessel function, $k_{\mathrm{BAO}}=\pi/l_{\mathrm{BAO}}$ is the wavemode
corresponding to the reference BAO scale
$l_{\mathrm{BAO}}=110\,h^{-1}\,\mathrm{Mpc}$, and $k_{\rm S}$ is the UV limit
of integration. To properly account for the resummation of IR modes at any
given scale $k$, one should integrate over all modes $q<k$, in a way that
$k_{\rm S}=k_{\rm S}(k)$. However, pushing the integration to significantly
large values of $k_{\rm S}$ would result in the breaking of the range of
validity of the pertubative IR expansion. At the same time, it can be shown
that the integrand of Equation (24) gives significant contributions only up to
$k_{\rm S}\sim 0.2\,h\,\mathrm{Mpc}^{-1}$ and therefore we restrict the
integration to the range $[0,\,0.2]$ when computing the value of $\Sigma$ Blas
_et al._ (2016); Ivanov  _et al._ (2019).

At next-to-leading order, the IR-resummed matter power spectrum can be written
as Blas  _et al._ (2016)

$\begin{split}P_{\rm NLO}(k)=&\,P_{\rm
nw}(k)+\left(1+k^{2}\Sigma^{2}\right)e^{-k^{2}\Sigma^{2}}P_{\rm w}(k)\\\
&+\left(P_{\rm nw}^{1\mbox{-}\mathrm{loop}}(k)+e^{-k^{2}\Sigma^{2}}P_{\rm
w}^{1\mbox{-}\mathrm{loop}}(k)\right),\end{split}$ (25)

where $P_{\rm nw}^{1\mbox{-}\mathrm{loop}}$ is the one-loop matter correction
defined in Equation (21) but evaluated using the smooth component $P_{\rm nw}$
rather than the full linear power spectrum $P_{\rm L}$, and $P_{\rm
w}^{1\mbox{-}\mathrm{loop}}=P^{1\mbox{-}\mathrm{loop}}-P_{\rm
nw}^{1\mbox{-}\mathrm{loop}}$.

####  II.2.3 Small-scale corrections: non-trivial stress tensor

Along with the resummation of infrared modes, we also have to consider the
impact of small-scale physics on long wavelength fluctuations. This happens
because the original assumption of a perfectly pressureless fluid is bound to
fail on non-linear scales, where dark matter experiences shell-crossing in
multistreaming regions Bernardeau  _et al._ (2002). Moreover, the effect of
baryonic physics, such as galaxy formation, cooling and feedback, also
contributes in generating a baryonic pressure that impacts the clustering of
dark matter on larger scales.

The net effect of UV scales on dark matter clustering is to generate a non-
zero stress tensor Pueblas and Scoccimarro (2009) whose leading contribution
to the matter power spectrum is to add a counter-term of the form Pueblas and
Scoccimarro (2009); Carrasco  _et al._ (2012); Baumann  _et al._ (2012)

$P_{\rm ctr}(k)=-2\,c_{1}\,k^{2}P_{\rm L}(k).$ (26)

Here $c_{1}$ can be treated as an effective speed of sound, that reflects the
influence of short wavelength perturbations, and in particular of the complex
physics beyond galaxy formation. Given the poor knowledge about these types of
processes, a standard assumption is to treat $c_{1}$ as a free parameter (see
e.g. Carrasco  _et al._ (2013); Baldauf  _et al._ (2015)) and marginalise over
it to obtain the posterior distribution of the parameters of interest.

By inspection of the individual terms contributing to one-loop formulas in
Equations (10) and (11), we can notice that the counter-term defined in
Equation (26) is completely degenerate with the leading order higher-
derivative contribution defined in Section II.1, as they both scale as
$k^{2}P_{\rm L}(k)$. In principle we may remove this degeneracy when jointly
fitting the auto and cross galaxy power spectra, as the non-vanishing stress
tensor affects only the one-loop matter power spectrum $P_{\rm mm}$. This
results in different imprints on $P_{\rm gg}$ and $P_{\rm gm}$,

$\begin{split}&P_{\rm
gg}(k)\supset-2b_{1}\left(b_{1}c_{1}+\beta_{1}\right)k^{2}P_{\rm L}(k)\,,\\\
&P_{\rm gm}(k)\supset-\left(2b_{1}c_{1}+\beta_{1}\right)k^{2}P_{\rm
L}(k)\,.\end{split}$ (27)

In practice, the sensitivity to UV modes might affect in different ways
$P_{\rm gg}$ and $P_{\rm gm}$, leading to inconsistent values of $c_{1}$
between the two observables. For this reason in the rest of the paper we will
employ two independent free parameters $\beta_{P}$ and $\beta_{P}^{\times}$
characterising the $k^{2}P_{\rm L}(k)$ contributions coming from $P_{\rm gg}$
and $P_{\rm gm}$, respectively.

####  II.2.4 Modeling of the non-linear matter power spectrum

In this section we briefly summarise the three different prescriptions we
adopt to model the non-linear matter power spectrum $P_{\rm mm}$ throughout
the rest of the paper.

The first of such models is based on a perturbative-simulated mixed approach,
which revolves around a fiducial high-resolution measurement of the non-linear
matter power spectrum from N-body simulations, and a two-loop perturbative
expansion in the cosmological parameter space for the response function
Bernardeau  _et al._ (2014); Nishimichi  _et al._ (2016). The latter
quantifies the variation of the non-linear power spectrum at scale $k$ induced
by a variation of the linear power spectrum at scale $q$, namely

$K(k,q)\equiv q\,\frac{\partial P_{\rm mm}(k)}{\partial P_{\rm L}(q)}.$ (28)

Under the assumption of having a reliable measurement of $P_{\rm mm}$ for a
fiducial cosmology $\bm{\theta}_{\mathrm{fid}}$, it is possible to predict the
same observable at a generic position $\bm{\theta}$ as

$\begin{split}P_{\rm mm}(k\,|\,\bm{\theta})=&\,P_{\rm
mm}(k\,|\,\bm{\theta}_{\mathrm{fid}})\int d(\log q)\,K(k,q)\\\
&\times\left[P_{\rm L}(q\,|\,\bm{\theta})-P_{\rm
L}(q\,|\,\bm{\theta}_{\mathrm{fid}})\right]\,.\end{split}$ (29)

The range of validity of this mixed approach becomes progressively less
accurate for cosmologies that are far way from $\bm{\theta}_{\mathrm{fid}}$,
but this issue can be overcome by employing a multi-step reconstruction
starting from the fiducial cosmology.

The RESPRESSO public package Nishimichi  _et al._ (2017) makes use of this
approach, starting from a fiducial measurements of $P_{\rm mm}$ from a set of
high-resolution N-body simulations with the Planck 2015 cosmology Planck
Collaboration (2016).

The second model we consider is based on the Effective Field Theory of Large
Scale Structure Carrasco  _et al._ (2012); Angulo  _et al._ (2015), and it is
close to what was recently used in the full shape analysis of the BOSS DR12
galaxy power spectrum Ivanov  _et al._ (2019). At one-loop in real-space, this
model is based on SPT results, but it also accounts for the effect of IR and
UV modes on the evolution of the matter power spectrum, as described in the
previous two sections. Namely, we can write a simple expression for the one-
loop matter power spectrum, that reads

$P_{\rm mm}(k)=P_{\rm NLO}(k)+P_{\rm ctr}(k)\,,$ (30)

where $P_{\rm NLO}(k)$ and $P_{\rm ctr}(k)$ are defined in Equations (25) and
(26), respectively.

The third model we consider is based on a particular flavour of Renormalised
Perturbation Theory (RPT) Crocce and Scoccimarro (2006a, b). In this kind of
approach, the non-linear matter power spectrum is typically separated into a
component that evolves the initial density contrast independently at each
wavelength, called propagator $G(k)$, and a mode-coupling term that accounts
for the mixing of scales due to non-linear evolution, so that we can write

$P(k)=G^{2}(k)P_{\rm L}(k)+P^{\rm MC}(k)\,.$ (31)

In a RPT-based approach (e.g. Bernardeau  _et al._ (2008); Taruya  _et al._
(2013)), the propagator is resummed while keeping the mode-coupling term at a
fixed order, leading to a breaking of the Galilean invariance (GI) of equal-
time correlators Scoccimarro and Frieman (1996). This translates into an
unphysical damping of the broadband power, which becomes mostly significant in
the UV regime. The RPT flavour we consider here, known as gRPT, effectively
attempts to resum the mode-coupling term in a way that is consistent with the
resummation of the propagator (see equation 17 of Eggemeier  _et al._ (2020)
for the explicit formula). At the same time, this approach naturally
incorporates IR resummation consistently with what is described in Section
II.2.3. We notice that, although even in this case the impact of the non-zero
stress tensor should be taken into account with the addition of an effective
speed of sound, the broadband predicted by gRPT is slightly suppressed with
respect to SPT predictions. This would in principle lead to an even smaller UV
counter-term, and for this reason we fix $c_{1}=0$ when modeling $P_{\rm mm}$
with gRPT.

Being partially calibrated on numerical measurements, RESPRESSO provides much
more accurate results on non-linear scales, as it intrinsically incorporates
higher-order corrections with respect to the previous two models, which
include only one-loop contributions to the matter density field. For this
reason, in the first part of this paper we fix the description of dark matter
non-linear evolution following this approach, and check the impact of using
different matter models only in Section IV.4. This choice is also motivated by
the analysis performed in Eggemeier  _et al._ (2020), that showed how
RESPRESSO is the model that most closely reproduces the performances of the
true matter power spectrum measured from simulations.

###  II.3 Stochasticity

One ingredient that is still missing from Equation (10) is the stochastic
contribution to the galaxy power spectrum. As a matter of fact, the previous
relations are completely deterministic, and assume a one-to-one correspondence
between the distribution of galaxies and the combined effect of the matter
density and tidal fields. However, galaxy formation is determined not only by
large-scales perturbations, but also by short wavelength modes. Under the
assumption of Gaussian initial conditions, these modes are completely
uncorrelated from large-scales fluctuations, and give birth to an additional
stochastic field $\varepsilon_{\rm g}$ which is dependent on the local
distribution of matter.

In practice, we can define the stochastic contribution $P_{\varepsilon_{\rm
g}\varepsilon_{\rm g}}$ to the power spectrum as Desjacques  _et al._ (2018)

$\mathcal{h}\varepsilon_{\rm g}(\bm{k})\,\varepsilon_{\rm
g}(\bm{k}^{\prime})\mathcal{i}=(2\pi)^{3}P_{\varepsilon_{\rm
g}\varepsilon_{\rm g}}(k)\,\delta_{\rm D}(\bm{k}+\bm{k}^{\prime}),$ (32)

and add this contribution to the one-loop galaxy power spectrum in Equation
(10). As for the modeling of $P_{\varepsilon_{\rm g}\varepsilon_{\rm g}}$, we
notice that the relation between the distribution of galaxies and the high-$k$
modes is not exactly local, as it depends on the distribution of matter within
a small finite region, similarly to the higher-derivatives of the galaxy
density field. For this reason we can write

$P_{\varepsilon_{\rm g}\varepsilon_{\rm
g}}(k)=\frac{1}{\overline{n}}\,\left(1+N_{0}+N_{2}\,k^{2}+\ldots\right)\,,$
(33)

where $\overline{n}$ is the mean number density of galaxies in the considered
volume. The constant contribution $N_{0}$ represents deviations from purely
Poissonian shot-noise ($1/\overline{n}$), while higher-order corrections, of
which the leading term scales as $k^{2}$, are generated to account for the
short-range non-locality described above.

$N_{0}$ is expected from the halo exclusion effect Sheth and Lemson (1999);
Smith  _et al._ (2007); Baldauf  _et al._ (2013) for which two different dark
matter halos cannot overlap (the same principle is what drives the effective
matter pressure in multi-streaming regions). This implies a deviation from
Poissonian shot-noise that can be either positive (super-Poisson) or negative
(sub-Poisson), depending on the considered tracer. Sub-Poissonian shot-noise
is more expected for central galaxies of massive halos, while super-Poissonian
values are more typical of galaxy populations with high satellite fractions
Baldauf  _et al._ (2013); Casas-Miranda  _et al._ (2002).

Notice that assuming Gaussian initial conditions all cross-correlators of the
form $\mathcal{h}\varepsilon_{\rm g}\,\delta\mathcal{i}$ are null by
construction, and therefore Equation (33) is the only stochastic contribution
to the galaxy auto power spectrum $P_{\rm gg}$. In reality, non-linear
gravitational evolution introduces a degree of correlation between long and
short wavelengths, so that $\mathcal{h}\varepsilon_{\rm
g}\,\delta\mathcal{i}\neq 0$ at later times, but all of these contributions
are subdominant in the case of the power spectrum, and can thus be neglected.

As anticipated in Section II.1, the quadratic term $P_{b_{2}b_{2}}$ needs to
be renormalized in order to provide a null contribution in the low-$k$ limit.
For this reason, we subtract from Equation (14) its large-scale asymptote
McDonald (2006); McDonald and Roy (2009), defined by

$P_{b_{2}b_{2}}^{\mathrm{noise}}=\frac{b_{2}}{2}\int_{\bm{q}}P_{\rm
L}^{2}(q)\,,$ (34)

and reabsorb it into the constant shot-noise parameter $N_{0}$.

Along with $P_{\rm gg}$ it can be shown that also $P_{\rm gm}$ requires an
additional stochastic component. In this case, stochasticity is not sourced by
$\varepsilon_{\rm g}$, as once again all correlators of the form
$\mathcal{h}\varepsilon_{\rm g}\,\delta\mathcal{i}$ vanish, but rather by the
matter density field itself via a new stochastic field $\varepsilon_{\rm m}$.
The reason of this is that dark matter ceases to behave as an ideal
pressureless fluid on small-scales, where the dynamics of gravitational
collapse is subject to shell crossing. This translates into an effective
pressure exerted by the matter density field, whose contribution to the
galaxy-matter cross power spectrum scales as $k^{2}$ in the low-$k$ limit
Akbar Abolhasani  _et al._ (2016). For this reason, it follows Desjacques  _et
al._ (2018)

$\mathcal{h}\varepsilon_{\rm g}(\bm{k})\,\varepsilon_{\rm
m}(\bm{k}^{\prime})\mathcal{i}=(2\pi)^{3}P_{\varepsilon_{\rm
g}\varepsilon_{\rm m}}(k)\,\delta_{\rm D}(\bm{k}+\bm{k}^{\prime})\,,$ (35)

where

$P_{\varepsilon_{\rm g}\varepsilon_{\rm
m}}(k)=\frac{1}{\overline{n}}\,(N_{2}^{\times}k^{2}+\ldots)\,.$ (36)

###  II.4 Co-evolution relations

Although the previously described galaxy bias expansion is complete at one-
loop, the large number of free bias parameters makes its applicability to real
datasets difficult, particularly when using the information contained only in
the two-point statistics (using additional constraints from e.g. the galaxy
bispectrum can break some of the degeneracies between parameters). For this
reason it is common practice to adopt empirical relations among bias
parameters in order to reduce the degrees of freedom of our model.

The most natural expressions are the so-called local Lagrangian relations
Catelan  _et al._ (1998, 2000); McDonald and Roy (2009); Matsubara (2011). The
latter are based on the assumption that galaxies (or more generally any biased
tracer of the matter density field) are formed instantaneously at an infinite
past time, that galaxy formation is driven exclusively by the local matter
density field, and that the number of tracers is conserved after their
formation. Under these assumptions, it is possible to describe the subsequent
evolution of the galaxy density field under the effect of gravity, which leads
to the appearance of higher-order non-local operators even in the presence of
a purely local relation at the time of formation.

It is possible to relax the Lagrangian local-in-matter-density assumption,
while still requiring the conservation of the total number of tracers. In this
way, we can express the Eulerian non-local parameters $\gamma_{2}$ and
$\gamma_{21}$ as a function of the corresponding Lagrangian counterparts,

$\displaystyle\gamma_{2}$
$\displaystyle=-\frac{2}{7}(b_{1}-1)+\gamma_{2,\mathcal{L}}\,,$ (37)
$\displaystyle\gamma_{21}$
$\displaystyle=-\frac{2}{21}(b_{1}-1)+\frac{6}{7}\gamma_{2}+\gamma_{21,\mathcal{L}}\,,$
(38)

where the subscript $\mathcal{L}$ indicates Lagrangian quantities, and the
remaining terms in the right-hand sides are the result of gravitational
evolution Chan  _et al._ (2012); Baldauf  _et al._ (2012); Eggemeier  _et al._
(2019).

Although local Lagrangian relations (i.e.
$\gamma_{2,\mathcal{L}}=\gamma_{21,\mathcal{L}}=0$) have been proven to be
much more accurate than simply neglecting non-local terms, recent measurements
of non-local bias parameters from a wide range of halo samples showed a slight
deviation with respect to observations Lazeyras and Schmidt (2018); Abidi and
Baldauf (2018). An alternative approach for the quadratic tidal parameter is
based on the excursion set theory Sheth  _et al._ (2013), for which
$\gamma_{2}$ can be predicted as a function of $b_{1}$ using a quadratic fit,

$\gamma_{2,\,\rm{ex}}(b_{1})=0.524-0.547\,b_{1}+0.046\,b_{1}^{2}\,.$ (39)

As shown in Figure 1 of Eggemeier  _et al._ (2020), this relation is more
accurate than the local Lagrangian approximation for tracers with
$b_{1}\gtrsim 1.3$. Therefore this relation should apply much better to our
datasets, which show a linear bias consistently larger than this value (see
Table 1). For this reason, we fix $\gamma_{2}$ to Equation (39) throughout the
rest of the paper.

##  III Data and methodology

###  III.1 Simulated galaxy samples

The robustness of our biasing scheme can be assessed only by validating the
model over a wide range of tracers, featuring different host halo masses,
galaxy bias, and redshifts. For this purpose, we make use of a set of three
different synthetic galaxy samples, whose main properties are summarised in
Table 1.

Table 1: Main properties of the three synthetic galaxy samples used in this work. The table shows the label we use to identify each sample, the N-body run on which it is based, the total number of independent realizations, redshift, galaxy number density, effective volume scaling ratio (see Section III.2), linear galaxy bias and large-scale deviation from Poisson shot-noise (in units of $1/\overline{n}$). In order to obtain fiducial values for the last two columns, we adopt the same strategy described in Eggemeier  _et al._ (2020) that makes use of the large-scale limit of the quantities $P_{\rm gg}$, $P_{\rm gm}$ and $P_{\rm mm}$. Identifier | Simulation | $N_{\rm R}$ | $z$ |  $\overline{n}$ $\left[\left(h/\mathrm{Mpc}\right)^{3}\right]$ | $\eta$ | $b_{1}$ |  $N_{0}$
---|---|---|---|---|---|---|---
MGS |  LasDamas Carmen | $40$ | $0.132$ | $1.1\times 10^{-3}$ | $7.04$ | $1.414\pm 0.003$ |  $-0.16\pm 0.07$
LOWZ |  LasDamas Oriana | $40$ | $0.342$ | $9.4\times 10^{-5}$ | $1$ | $2.235\pm 0.012$ | $-0.176\pm 0.018$
CMASS | Minerva | $100$ | $0.57$ | $4.0\times 10^{-4}$ | $2.47$ | $2.022\pm 0.003$ |  $-0.29\pm 0.02$
Table 2: Cosmological parameters of the dark matter simulations used to generate the galaxy samples analysed in this work. Columns show the matter, dark energy and baryon density, Hubble constant, scalar index, and rms density fluctuations within a sphere of radius $8\,h^{-1}\,\mathrm{Mpc}$. Simulation | $\Omega_{m}$ | $\Omega_{\Lambda}$ | $\Omega_{b}$ | $h$ | $n_{s}$ |  $\sigma_{8}$
---|---|---|---|---|---|---
LasDamas | 0.25 | 0.75 | 0.04 | 0.7 | 1.0 | 0.8
Minerva | 0.285 | 0.715 | 0.046 | 0.695 | 0.9632 | 0.828

The three catalogues were generated by populating dark matter halos with
galaxies using Halo Occupation Distribution (HOD) prescriptions. Since the HOD
parameters were calibrated to obtain number densities consistent with the ones
of pre-existing real observations, for sake of easiness we label our simulated
catalogues with the name of the corresponding data samples. Nevertheless, we
remind the reader that here we make use of the full comoving snapshot volume,
without considering any selection effect, as the goal is to investigate the
range of validity of 1-loop galaxy bias, leaving aside the impact of
observational systematics.

The CMASS sample is based on the Minerva simulations Grieb  _et al._ (2016),
which consist in a set of 100 realizations of a
$(1500\,h^{-1}\,\mathrm{Mpc})^{3}$ cubic box with periodic boundary
conditions. The simulations were run using the public Gadget code Springel
_et al._ (2001); Springel (2005), that regulated the motion of $1000^{3}$ dark
matter particles within the aforementioned volume. Initial conditions were set
up using the linear power spectrum obtained from CAMB Lewis and Bridle (2002)
and displacing particles according to second-order Lagrangian Perturbation
Theory (2LPT) Crocce  _et al._ (2006). The comoving snapshot has a redshift of
$z=0.57$ and it is meant to reproduce the properties of the BOSS CMASS galaxy
sample Alam  _et al._ (2017).

The LOWZ $(z=0.342)$ and MGS $(z=0.132)$ samples are based on the Oriana and
Carmen boxes of the LasDamas N-body simulations McBride  _et al._ (2009);
Sinha  _et al._ (2018). These are a set of 40 independent realizations with
periodic boundary conditions, with a volume of
$(2400\,h^{-1}\,\mathrm{Mpc})^{3}$ and $(1000\,h^{-1}\,\mathrm{Mpc})^{3}$, and
a mass resolution of $4.6\times 10^{11}\,h^{-1}\,M_{\odot}$ and $4.9\times
10^{10}\,h^{-1}\,M_{\odot}$, respectively. Initial conditions were also set up
using 2LPT, but in this case the initial power spectrum is computed using
CMBFast Seljak and Zaldarriaga (1996). The HOD parameters of these synthetic
catalogues were selected to reproduce the properties of the BOSS LOWZ and SDSS
Main Galaxy Sample (MGS) at $M_{r}<-21$.

The complete set of cosmological parameters for Minerva and LasDamas is
summarised in Table 2.

###  III.2 Measurements of power spectra and their covariances

We make use of the estimator described in Sefusatti  _et al._ (2016), based on
fourth order particle assignment scheme and interlacing optimization to reduce
the effect of large modes aliasing, to measure the galaxy auto power spectrum
$P_{\rm gg}$ and the galaxy-matter cross power spectrum $P_{\rm gm}$ for all
the tracers described in Table 1. We adopt a linear $k$ binning from $k_{\rm
min}=k_{f}$, where $k_{f}=2\pi/L$ is the fundamental frequency of the box of
size $L$, to $k_{\rm max}=k_{\rm Nyq}$, where $k_{\rm Nyq}=\pi N_{\rm grid}/L$
is the Nyquist frequency corresponding to a FFT grid of size $N_{\rm grid}$.
We use a linear binning with step $\Delta_{k}=k_{f}$ for CMASS and LOWZ, and
$\Delta_{k}=2k_{f}$ for MGS. Since we model the stochastic contribution to the
galaxy auto power spectrum in terms of deviations from Poisson shot-noise, we
correct our raw measurements of $P_{\rm gg}$ by subtracting the constant
factor $P^{\mathrm{noise}}=1/\overline{n}$.

Our final data vector consists in the ensemble average over the number of
independent realizations $N_{\rm R}$ (100 for CMASS, 40 for both LOWZ and
MGS). First, we define the averaged observable

$\overline{X}_{i}=\frac{1}{N_{\rm R}}\sum_{n=1}^{N_{\rm R}}X^{(n)}_{i}$ (40)

where $X\in\\{P_{\rm gg},\,P_{\rm gm}\\}$, $i$ is the index running through
the $k$ binning, and the superscript $n$ refers to the $n$-th realization of
the considered observable. From this definition we estimate the auto- and
cross- covariance matrices as

$C_{X\times Y,ij}=\frac{1}{N_{\rm R}}\sum_{n=1}^{N_{\rm
R}}\left(X_{i}^{(n)}-\overline{X}_{i}\right)\,\left(Y_{j}^{(n)}-\overline{Y}_{j}\right)\,,$
(41)

where once again $X,\,Y\in\\{P_{\rm gg},\,P_{\rm gm}\\}$ .

Similarly to what was done in Eggemeier  _et al._ (2020) we retain only the
diagonal entries of the previously defined covariance matrices. This
approximation is justified given the significant low number density of the
synthetic samples we consider (actual values are listed in Table 1). These
numbers translate in a significant shot-noise contribution to the power
spectrum, whose main effect is to boost the diagonal entries of the covariance
matrix, leading to sub-dominant off-diagonal terms. Therefore we approximate
our covariance matrices with a block-diagonal shape, and consider only the
auto- and cross-correlation at the same wavelength (i.e. $C_{ij}=0$ for
$k_{i}\neq k_{j}$). Additionally, in order to reduce noise due to the limited
number of independent realizations, we compare the raw covariance matrix to
Gaussian predictions Grieb  _et al._ (2016) for each $k$ bin, and retain the
maximum of the two values.

The error budgets of our galaxy samples are subsequently rescaled in order to
match the same _effective_ volume Feldman  _et al._ (1994); Tegmark (1997),
defined as

$V_{\mathrm{eff}}(k^{\star})=\left[\frac{\overline{n}\,P_{\rm
gg}(k^{\star})}{1+\overline{n}\,P_{\rm gg}(k^{\star})}\right]^{2}\,V\,,$ (42)

where $\overline{n}$ and $V$ are the mean number density and volume of the
considered sample, respectively, and we choose the reference scale to evaluate
the effective volume as $k^{\star}=0.1\,h\,\mathrm{Mpc}^{-1}$. In this way the
constraining power and the signal-to-noise ratio of the three galaxy
catalogues are artificially set to roughly match the same amplitude. We choose
to rescale the covariance of both CMASS and MGS to match the effective volume
of the LOWZ sample $\left(\approx 6\,(\mathrm{Gpc}/h)^{3}\right)$. Table 1
shows the rescaling factor for the three samples, that is simply defined as
the sample-to-LOWZ ratio between the respective effective volumes.

###  III.3 Fitting procedure and prior choices

We make use of the large number of independent realizations for each synthetic
galaxy sample to define the overall likelihood function as the product of the
individual $N_{\rm R}$ likelihoods Eggemeier  _et al._ (2020); Oddo  _et al._
(2020). In this way, under the assumption that the individual likelihoods are
well described by a multivariate Gaussian distribution, we can define the
overall likelihood function as

$\begin{split}-2\log{{\cal L}_{\mathrm{tot}}}&=-\frac{2}{N_{\rm
R}}\sum_{n=1}^{N_{\rm R}}\log{{\cal L}_{(n)}}=\frac{1}{N_{\rm
R}}\sum_{n=1}^{N_{\rm R}}\chi_{(n)}^{2}\\\ &=\frac{1}{N_{\rm
R}}\sum_{n=1}^{N_{\rm R}}\sum_{i,j=1}^{N_{\rm
b}}\left(X_{i}^{(n)}-\mu_{i}\right)\,C_{X,ij}^{-1}\,\left(X_{j}^{(n)}-\mu_{j}\right)\,,\end{split}$
(43)

where $X_{i}^{(n)}$ is the measurement from the $n$-th realization of the
$i$-th $k$ bin, $\mu_{i}$ is the corresponding model prediction, and $C_{ij}$
is the rescaled covariance matrix as described in Sec. III.2. With this
definition, our likelihood ends up coinciding with the likelihood of the mean
of the data with covariance $C_{ij}$. In practice, there is a constant factor
between the two definitions, that depends on the number of independent
realizations $N_{\rm R}$ and the number of data points $N_{\rm b}$. This
factor is taken into account when deriving the goodness of fit for each of the
tested configurations.

The inference of model parameters is carried out through a least-$\chi^{2}$
analysis based on a standard Metropolis-Hastings MCMC algorithm. The
likelihood incorporates a complete recipe for one-loop galaxy clustering and
an interface to CAMB. For each tested configuration of the MCMC, we first run
some preliminary chains that are needed to obtain a robust estimate of the
parameter covariance, which is essential for the good convergence of the
Markov chain in a highly dimensional parameter space. We iterate this process
twice, each time specifying the parameter covariance obtained at the previous
step, before running the final set of chains. These are terminated as soon as
the chains satisfy a standard Gelman-Rubin convergence criterium, i.e. as soon
as the between-chain and within-chain variances (we run 8 independent chains
for each case, changing the initial random seed) are in agreement within
$R<1.1$ Gelman and Rubin (1992). After removing the burn-in, these chains are
combined into a single object, which we pass to the software package getdist
Lewis (2019) to extract the parameter posteriors.

Table 3: Adopted priors for the full list of cosmological and nuisance
parameters of the fits. We impose a uniform prior ($\mathrm{U}$) for all the
sampling parameters. The scalar amplitude $A_{\rm s}$ is not sampled over when
fitting only $P_{\rm gg}$ but it is when adding the additional constraint from
$P_{\rm gm}$. In all cases the second-order tidal bias $\gamma_{2}$ is fixed
to the excursion set relation defined in Equation (39).
COSMOLOGY
---
$h$ | $\mathrm{U}[0.5,\,1]$
$\Omega_{\rm c}h^{2}$ | $\mathrm{U}[0.08,\,0.16]$
$10^{9}A_{\rm s}$ | $P_{\rm gg}:$ | $P_{\rm gg}+P_{\rm gm}:$
Fixed to fiducial | $\mathrm{U}[0.7,\,3.3]$
BIAS
$b_{1}\sigma_{8}$ | $\mathrm{U}[0.25,\,4]$
$b_{2}\sigma_{8}^{2}$ | $\mathrm{U}[-10,\,10]$
$\gamma_{2}$ | Fixed to $\gamma_{2,ex}(b_{1})$
$\gamma_{21}\sigma_{8}^{3}$ | $\mathrm{U}[-8,\,8]$
$\beta_{P}\,[k_{\mathrm{HD}}^{-2}]$ | $\mathrm{U}[-10,\,10]$
$\beta_{P}^{\times}\,[k_{\mathrm{HD}}^{-2}]$ | $\mathrm{U}[-10,\,10]$
SHOT-NOISE
$N_{0}\,[\overline{n}^{-1}]$ | $\mathrm{U}[-1,\,0.5]$
$N_{2}\,[k_{\mathrm{HD}}^{-2}\,\overline{n}^{-1}]$ | $\mathrm{U}[-10,\,10]$
$N_{2}^{\times}\,[k_{\mathrm{HD}}^{-2}\,\overline{n}^{-1}]$ | $\mathrm{U}[-10,\,10]$
COUNTER-TERMS
$c_{1}$ | $\mathrm{U}[0,\,10]$

The models described in Section II embed a large number of free parameters. In
this analysis we vary both cosmological and nuisance terms (the latter include
bias, shot-noise and matter counter-term). Since we are constraining two-point
statistics only via galaxy clustering, we fix both the baryon density
parameter $\Omega_{\rm b}h^{2}$ and the primordial spectral index $n_{\rm s}$
to their fiducial values (listed in Table 2), and vary the CDM density
parameter $\Omega_{\rm c}h^{2}$ and the Hubble constant $h$. We add the
primordial scalar amplitude $A_{\rm s}$, but only when fitting the combination
of $P_{\rm gg}$ and $P_{\rm gm}$, since the fit of the auto galaxy power
spectrum alone cannot break the strong degeneracy between $A_{\rm s}$ and
$b_{1}$. On the contrary, while both $P_{\rm gg}$ and $P_{\rm gm}$ have the
same dependency on $A_{\rm s}$, they scale differently with the linear bias,
i.e. $P_{\rm gg}\propto b_{1}^{2}$ and $P_{\rm gm}\propto b_{1}$ in the linear
limit. We choose a flat prior for the three cosmological parameters that is
large enough to fully contain their posterior distribution up to $2\sigma$
even for the least constraining configuration.

The same choice of uniform priors is chosen for the nuisance parameters. We
sample the combination of the bias parameters with the corresponding power of
$\sigma_{8}$ in order to avoid strong degeneracies between the one-loop bias
expansion and the intrinsic non-linearity of the matter power spectrum
$P_{mm}$, that at first order is well captured by $\sigma_{8}$. As anticipated
in Section II.4, we fix $\gamma_{2}$ to Equation (39), and vary all other bias
parameters freely. We quote higher-derivatives parameters with respect to an
arbitrary fiducial scale $k_{\mathrm{HD}}=0.4\,h\,\mathrm{Mpc}^{-1}$, while
all the shot-noise parameters can be naturally expressed in units of
$1/\overline{n}$. Since terms involving $N_{2}$ and $N_{2}^{\times}$ both
carry a $k^{2}$-dependency, we also express the latter in units of
$k_{\mathrm{HD}}$.

The full list of priors for the model parameters is shown in Table 3.

###  III.4 Rescaling of the input linear power spectrum

Figure 1: Ratio between the CMBFast and CAMB linear prediction at the
reference cosmology of the LasDamas simulations. The two shaded grey bands
mark the intrinsic $1\mbox{-}$sigma error of the LOWZ and MGS galaxy power
spectrum $P_{\rm gg}$.

As discussed in Section III, Minerva and LasDamas adopt different Boltzmann
solvers to obtain the power spectra for the initial particle displacements. In
particular, the Minerva runs assume initial conditions based on CAMB, while
the LasDamas simulations make use of the power spectrum computed with CMBFast.
Figure 1 shows the fractional deviation between the output of the two codes at
the reference cosmology of the LasDamas runs and at the redshift of the MGS
sample. This difference is compared to LOWZ and MGS intrinsic signal-to-noise
ratio, represented by the two shaded areas. The two spectra deviate at the
level of $\sim 0.5\%$ in the range of scales that are relevant to this
analysis, and this might represent an issue when modeling the galaxy-galaxy
and galaxy-matter power spectra. For this reason, when fitting LOWZ and MGS,
we employ a rescaling scheme based on the renormalization of the input power
spectrum of CAMB by the ratio shown in Figure 1. This approximation is valid
as long as the sampled cosmology is not far from the fiducial one, and this
becomes a completely legitimate assumption when hitting exactly the true
cosmology. However, we claim that even for cosmologies that are slightly off
with respect to the fiducial one, this rescaling scheme is already better than
directly using the CAMB prediction with no further correction.

In order to perform this rescaling, we first apply the units transformation
$h\,\mathrm{Mpc}^{-1}\rightarrow\mathrm{Mpc}^{-1}$. This is meant to assure
that the power spectrum ratio is obtained at the same set of scales, avoiding
the dependency on the Hubble parameter $h$ Sánchez (2020). At each step of the
Markov chain, the current CAMB prediction is firstly transformed into
$\mathrm{Mpc}^{-1}$ units, rescaled using the approach described above, and
finally transformed back into $h\,\mathrm{Mpc}^{-1}$ units.

###  III.5 Performance metrics

The validation of one-loop galaxy bias as described in Section II is mostly
based on its range of validity. In order to select the maximum mode at which
our model stops providing good performances, we run multiple MCMC chains
varying $k_{\rm max}$ in the interval $(0.1,\,0.3)\,h\,\mathrm{Mpc}^{-1}$ with
a step of $0.025\,h\,\mathrm{Mpc}^{-1}$, leading to 9 different fitting ranges
for each configuration. We stop at $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ as
pushing the model to even smaller scales would inevitably require accounting
for two-loop contributions in the description of galaxy bias and non-linear
matter evolution. Clearly, exploiting the small-scale information contained in
the non-linear regime naturally provides a better constraining power for the
model parameters, but at the cost of introducing theory systematics in their
posterior distribution. As described in Eggemeier  _et al._ (2020) (see Osato
_et al._ (2019) for a similar approach), we make use of a set of three
different performance metrics, whose goal is to provide a quantitative way to
compare the posteriors extracted from our Markov chains, and to determine when
a particular model no longer gives an accurate description of our datasets.

####  III.5.1 Figure of bias

The first quantity we are interested in evaluating is the ability of the model
to provide unbiased measurements of its free parameters. In particular, in
this work we focus mostly on the systematic error of the cosmological
parameters: the Hubble constant $h$, the cold dark matter density $\Omega_{\rm
c}h^{2}$, and, when considering the joint fit between $P_{\rm gg}$ and $P_{\rm
gm}$, the primordial scalar amplitude $A_{\rm s}$.

We define the Figure of Bias (FoB) of the considered model for a given
parameter set $\theta$ as

$\mathrm{FoB}(\bm{\theta})\equiv\left[\left(\,\overline{\bm{\theta}}-\bm{\theta}_{\mathrm{fid}}\right)^{\intercal}{\mathbb{S}}^{-1}\left(\,\overline{\bm{\theta}}-\bm{\theta}_{\mathrm{fid}}\right)\right]^{1/2},$
(44)

where $\overline{\bm{\theta}}$ represents the mean of the posterior
distribution, $\bm{\theta}_{\mathrm{fid}}$ represents the fiducial position,
and ${\mathbb{S}}$ is the parameter covariance expressed in matrix form. In
this way, we are simply quantifying the relative separation of the measured
parameter from its true value in terms of the variance of the posterior
distribution. If the FoB is evaluated only for one parameter, then the
standard $68\%-95\%$ percentile thresholds correspond to a FoB of $1-2$,
respectively. On the contrary, when more than one parameter is considered to
compute the FoB, the two percentiles can be calculated by direct integration
of a multivariate Gaussian distribution with the corresponding dimensionality.
For the case we consider $n=2\,(n=3)$ it follows that the $68\%-95\%$
percentiles correspond to a FoB of $1.52-2.49\,(1.88-2.83)$, respectively.

####  III.5.2 Goodness-of-fit

Along with finding unbiased constraints on the parameters of interest, we also
ask our model to provide a good description of the observables we use in the
fit. We quantify the goodness-of-fit in terms of the standard $\chi^{2}$
extracted from the Markov chains (Equation (43)), but after rescaling it to
account for the additional factor $\eta$ introduced in the covariance of the
data (see Section III.2). With this rescaling, we can compare the recovered
$\chi^{2}$ to the $68\%$ and $95\%$ percentiles of a $\chi^{2}$ distribution
with $N_{\mathrm{dof}}$ degrees of freedom, where

$N_{\mathrm{dof}}=N_{\rm R}\times N_{\rm b}-N_{\rm p}\,.$ (45)

Here $N_{\rm R}$ is the number of independent realizations for each galaxy
sample, $N_{\rm b}$ is the total number of $k$ bins included in the fit, and
$N_{\rm p}$ is the number of free model parameters.

####  III.5.3 Figure of merit

The final metric we employ is based on the merit of the considered model with
respect to the cosmological parameters that are varied in the fit. We define
the Figure of Merit (FoM) of a given model for a subset of parameters
$\bm{\theta}$ as the inverse of the determinant of the parameter covariance
matrix corresponding to the subset $\bm{\theta}$, i.e.

$\mathrm{FoM}(\bm{\theta})\equiv\frac{1}{\sqrt{\mathrm{det}\left[\,\overline{\mathbb{S}}(\bm{\theta})\right]}}\,.$
(46)

Here
$\overline{\mathbb{S}}(\bm{\theta})\equiv{\mathbb{S}}(\bm{\theta})/\bm{\theta}^{2}_{\mathrm{fid}}$
is the block corresponding to the parameter subset $\bm{\theta}$ rescaled by
their fiducial values. With this last normalization we are explicitly asking
for a relative FoM rather than an absolute quantity, that would be more
difficult to compare to the one of other parameters.

##  IV Testing one-loop galaxy bias

Table 4: List of all the considered combinations between the modeling of $P_{\rm mm}$ and the extensions to one-loop galaxy bias, as described in Section II. For each model, we denote with a check mark all the free parameters of the $P_{\rm gg}$ fit, while we use a double check mark for the free parameters that are exclusive to the joint $P_{\rm gg}+P_{\rm gm}$ fit. $P_{mm}$ | RESPRESSO | gRPT | EFT
---|---|---|---
Bias model | STD | k2N | HD | STD | k2N | HD | STD | k2N
$h$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$\Omega_{\rm c}h^{2}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$A_{\rm s}$ | ✓✓ | ✓✓ | ✓✓ | ✓✓ | ✓✓ | ✓✓ | ✓✓ | ✓✓
$b_{1}\sigma_{8}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$b_{2}\sigma_{8}^{2}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$\gamma_{21}\sigma_{8}^{3}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$\beta_{P}$ |  |  | ✓ |  |  | ✓ |  |
$\beta_{P}^{\times}$ |  |  | ✓✓ |  |  | ✓✓ |  |
$N_{0}$ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
$N_{2}$ |  | ✓ |  |  | ✓ |  |  | ✓
$N_{2}^{\times}$ |  | ✓✓ |  |  | ✓✓ |  |  | ✓✓
$c_{1}$ |  |  |  |  |  |  | ✓ | ✓

In this section we analyse parameter posteriors extracted from several Markov
chains, where we vary both the maximum fitting scale $k_{\rm max}$ and the
final expression we adopt to model one-loop galaxy power spectra. The complete
list of model configurations is summarised in Table 4, along with the total
number of free parameters in each case. Based on results from Eggemeier  _et
al._ (2020), we decide to employ RESPRESSO to model $P_{\rm mm}$ in this
section, and to leave the comparison between different matter models to
Section IV.4. In this way we can focus exclusively on assessing the validity
of one-loop bias in real-space leaving aside both the impact of non-linear
matter evolution and redshift-space distortions (that will be the topic of a
future paper).

For each configuration, we compute the three performance metrics defined in
Section III.5 and show their trends as a function of $k_{\rm max}$. At the
same time, we perform sanity checks on the linear bias $b_{1}\sigma_{8}$ and
the large-scale shot-noise parameter $N_{0}$, for which we have fiducial
values to compare with (see Table 1). As anticipated in Section II.4,
throughout the rest of the paper we fix the second-order non-local bias
parameter $\gamma_{2}$ to the excursion-set relation (Equation (39)), in order
to break the strong degeneracy between the former and $\gamma_{21}$. This
choice is also based on results from Eggemeier  _et al._ (2020), where it was
shown how this choice led to stable and overall accurate measurement of the
linear bias.

In the next two sections we analyse the posterior derived from fits of the
auto galaxy power spectrum $P_{\rm gg}$ alone, and from the combination of the
auto and cross galaxy power spectra, $P_{\rm gg}+P_{\rm gm}$. In the latter
case, the additional information contained in $P_{\rm gm}$ allows to break the
degeneracy between the linear bias and the amplitude of primordial
fluctuations, allowing us to also sample over $A_{\rm s}$.

In both cases, we adopt a criterion for defining the range of validity of a
given model, which is based on a combination of FoB and goodness-of-fit.
Following Eggemeier  _et al._ (2020), we define the model-breaking statistics
as

$\sigma_{\mathrm{MB}}(k)\equiv\frac{\mathrm{FoB}(k)}{\mathrm{FoB}_{68\%}}+\frac{\chi^{2}(k)-1}{\chi^{2}_{95\%}(k)-1},$
(47)

where the subscript percentages correspond to the percentiles of the
corresponding distribution (FoB and $\chi^{2}$). We say that a given model
breaks down at a scale $k^{\dagger}$ when the model-breaking statistics
assumes a critical value $\sigma_{\mathrm{crit}}$ at $k^{\dagger}$, i.e.

$\sigma_{\mathrm{MB}}(k^{\dagger})=\sigma_{\mathrm{crit}}\,.$ (48)

We arbitrary choose the critical threshold $\sigma_{\mathrm{crit}}=1.5$. In
this way we would accept models with a maximum FoB one and a half times larger
than the value corresponding to the $1$-sigma of the FoB distribution, but
only under a perfect recovery of the shape of the input dataset. Practically,
this case is unrealistic at high $k_{\rm max}$, as the $\chi^{2}$
progressively deviates from the number of degrees of freedom because of the
weakening of the model. For this reason $\sigma_{\mathrm{MB}}$ receives
individual contributions from the FoB and the goodness-of-fit.

###  IV.1 Validity of one-loop galaxy bias for the auto power spectrum

Figure 2: FoM, FoB and goodness-of-fit extracted from chains corresponding to
the fits of the galaxy auto power spectrum $P_{\rm gg}$. Performance metrics
refer to the combinations of $h$ and $\Omega_{\rm c}h^{2}$, which are the only
two free cosmological parameters of the model. The standard (STD),
$k^{2}$-dependent shot-noise (k2N), and higher derivatives (HD) models are
colour-coded as shown in the legend. In all three cases, the non-linear matter
power spectrum is modelled with RESPRESSO, and the second-order tidal bias
$\gamma_{2}$ is fixed to Equation (39). The two shaded grey regions mark the
$68$th and $95$th percentiles of the FoB and $\chi^{2}$ distribution.

In this section we analyse results from the fits to the auto power spectrum
$P_{\rm gg}$. We fix the description of the non-linear matter power spectrum
and test different extensions of one-loop galaxy bias. The standard (STD)
model is based on two free cosmological parameters $(h,\,\Omega_{\rm c}h^{2})$
plus four free nuisance parameters
$(b_{1}\sigma_{8},\,b_{2}\sigma_{8}^{2},\,\gamma_{21}\sigma_{8}^{3},\,N_{0})$.
The $k^{2}$-dependent shot-noise (k2N) and higher-derivatives (HD) extensions
have one additional free parameter each, i.e. $N_{2}$ and $\beta_{P}$,
respectively.

In Figure 2 we show the three performance metrics extracted from the three
different $P_{\rm gg}$ measurements (CMASS, LOWZ and MGS), plotted against the
maximum scale $k_{\rm max}$ included in the fit. In this case we show the FoB
and FoM corresponding to the combination $h,\,\Omega_{\rm c}h^{2})$, i.e. to
all the cosmological parameters that are varied in the model. The two grey
shaded areas in the panels corresponding to FoB and goodness-of-fit mark the
$68$th and $95$th percentiles of the corresponding quantity. As described in
Section III.5.1, the $68\%-95\%$ FoB limits do not correspond to values of FoB
$=1$ and FoB $=2$, but rather have to be computed integrating a multivariate
Gaussian distribution with $n=2$ number of dimensions.

For the three galaxy samples, we find that the STD model (blue line) alone is
enough to provide a good description of the dataset in terms of goodness-of-
fit, as the reduced $\chi^{2}$ constantly stays within the $68\%$ (dark grey
shaded area) of the corresponding $\chi^{2}$ distribution. At the same time,
the combined constraint on the $(h,\,\Omega_{\rm c}h^{2})$ pair is unbiased,
showing a multivariate posterior distribution that is able to capture the
fiducial position at typically $\sim 0.5\,\sigma$. Adding either an additional
$k^{2}$-dependent shot-noise (red line) or the leading higher derivatives
correction (green line) does not significantly alter the two metrics, while
reducing the amplitude of the corresponding FoM ($\sim 10\%$ when adding the
stochastic contribution $N_{2}$, $\sim 20-30\%$ when adding the higher
derivatives parameter $\beta_{P}$). Even for the STD model, for which we
expect the model to fail earlier than the other two extensions) the model-
breaking scale $k^{\dagger}$ is higher than the maximum scale we probe
$(k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1})$. However, we decide to stop at this
scale without exploring smaller scales, even if the model-breaking criterion
were most likely satisfied at higher $k_{\rm max}$. As a matter of fact it
would be hard to separate the effective goodness of the model from the impact
of two-loop contributions to galaxy bias.

Figure 3: Marginalised constraints on the linear bias $b_{1}\sigma_{8}$ (top)
and the constant stochastic term $N_{0}$ (bottom) as a function of
$k_{\mathrm{max}}$, for the fits of the galaxy auto power spectrum $P_{\rm
gg}$. In all cases dark matter non-linear evolution is modelled using
RESPRESSO. Fiducial values of the linear bias and constant shot-noise
contribution are reproduced with black dashed lines. Uncertainties on the
fiducial values are marked with shaded grey bands (marking $1$-and
$2$-$\sigma$ confidence intervals). Fiducial measurements are listed in Table
1.

We notice how, similarly to what was obtained in Eggemeier  _et al._ (2020)
when focusing on the linear galaxy bias, there is a tendency of the FoM to
flatten above an approximate scale of $k_{\rm
max}=0.25\,h\,\mathrm{Mpc}^{-1}$. This is an indication that pushing the non-
linear model to increasingly higher wave modes does not necessarily imply a
more stringent measurement of the cosmological parameters, as most of the
additional constraining power is absorbed by nuisance parameters such as
higher-order galaxy bias contributions.

In Figure 3 we show the dependence on $k_{\rm max}$ of the linear bias
parameter $b_{1}\sigma_{8}$ and the constant stochastic contribution $N_{0}$,
for which we have fiducial values obtained exploiting the large-scale limit of
the measured galaxy and matter power spectra. Although the extension to either
the k2N or the HD model enlarges the size of the errorbars, we notice how the
former model is the only one capable of providing a simultaneous unbiased
measurements of the two parameters, for the three galaxy samples we are
considering. In particular the LOWZ sample shows a clear detection of $N_{2}$,
that, if not accounted for, can lead to a significant $(>2$-$\sigma)$
systematic effect on $N_{0}$ above $k_{\rm max}\sim
0.25\,h\,\mathrm{Mpc}^{-1}$, while partially recovering the true amplitude of
$P_{\rm gg}$ with an underestimation of $b_{1}\sigma_{8}$. The only exception
is represented by the MGS sample, whose linear bias is constantly
overestimated by a $\sim 2\%$ factor for all the three modeling assumptions,
and that is consistent with the fiducial value only up to $2$-$\sigma$.

###  IV.2 Consistency between auto and cross power spectra

Figure 4: Same as Figure 2 but for the combinations of the galaxy auto and
cross power spectrum, $P_{\rm gg}$ and $P_{\rm gm}$. In this case, performance
metrics (FoM and FoB) refer to the combinations of the three cosmological
parameters $(h,\,\Omega_{\rm c}h^{2},\,A_{\rm s})$. The standard (STD),
$k^{2}$-dependent shot-noise (k2N), and higher derivatives (HD) models are
colour-coded as shown in the legend. The solid-to-dashed transition mark the
model-breaking scale $k^{\dagger}$. In all three cases, the non-linear matter
power spectrum is modelled with RESPRESSO. The two shaded grey regions mark
the $68$th and $95$th percentiles of the FoB and $\chi^{2}$ distribution.

In this section we focus on a much more stringent test, corresponding to the
simultaneous fit of both the galaxy auto power spectrum $P_{\rm gg}$ and the
galaxy-matter cross power spectrum $P_{\rm gm}$. As a matter of fact, a good
accuracy for this combined statistics might become crucial in analyses that
aim to exploit the whole information contained in galaxy clustering and
galaxy-galaxy weak lensing, as most of the upcoming large observational
projects are going to do ($3\times 2$-point analysis, see e.g. Abbott  _et
al._ (2018); Heymans  _et al._ (2020); Pandey  _et al._ (2020)).

The different scaling of $P_{\rm gg}$ and $P_{\rm gm}$ on the linear bias
$b_{1}$ leads to the breaking of the strong $b_{1}-A_{\rm s}$ degeneracy.
Therefore, we additionally sample the scalar amplitude $A_{\rm s}$,
effectively extending by one the number of degrees of freedom of our model.
Also in this case, we concentrate on the one-loop biasing scheme, fixing the
description of the matter clustering to the output of RESPRESSO. We test the
standard (STD) model against the $k^{2}$-dependent shot-noise (k2N) and higher
derivatives (HD) templates, for which we introduce two more degrees of freedom
(one free parameter for $P_{\rm gg}$ and $P_{\rm gm}$, each). In summary, the
STD model is based on three free cosmological parameters $(h\,,\Omega_{\rm
c}h^{2},\,A_{\rm s}$) plus four free nuisance parameters
$(b_{1}\sigma_{8},\,b_{2}\sigma_{8}^{2},\,\gamma_{21}\sigma_{8}^{3}\,,N_{0})$.
The k2N and HD models have two additional free parameters each, i.e.
$(N_{2},\,N_{2}^{\times})$ and $(\beta_{P},\,\beta_{P}^{\times})$,
respectively. Note that that the second-order tidal bias $\gamma_{2}$ is not a
free parameter, but it is fixed to the excursion-set relation defined in
Equation (39) for all the cases we consider.

In Figure 4 we show the three performance metrics in a similar way as we did
for the $P_{\rm gg}$-only fits. In this case both FoM and FoB are computed
from the combinations of the three cosmological parameters $(h,\,\Omega_{\rm
c}h^{2},\,A_{\rm s})$. We first notice how the use of the STD model is no
longer sufficient to provide an accurate recovery of the cosmological
parameters on the overall range of scales we are considering. Indeed, the FoB
of both CMASS and LOWZ quickly increases above $2$-$\sigma$ at an approximate
scale of $k_{\rm max}\sim 0.23\,h\,\mathrm{Mpc}^{-1}$ and keeps getting worse
at higher wave modes. The same trend is observed in panels referring to the
goodness-of-fit, as the reduced $\chi^{2}$ starts to fall outside the $95\%$
confidence region at approximately the same scale cut. According to the model-
breaking criterion, the range of validity of the STD model is limited to
scales below $k\sim 0.2\,h\,\mathrm{Mpc}^{-1}$ (this is reflected in the plot
by the solid-to-dashed line transition).

Adding $k^{2}$-dependent stochasticity increases both the accuracy in the
recovery of the cosmological parameters and the overall likelihood between
data vectors and best-fitting model. Indeed, the k2N model shows a FoB which
is constantly within $1$-$\sigma$ for the three galaxy samples. Moreover, the
goodness-of-fit is significantly improved with respect to the STD model, with
the $\chi^{2}$ measured at the highest $k_{\rm max}$ still being consistent
within $1$-$\sigma$. The combination of FoB and goodness-of-fit can be
observed via the model-breaking criterion on the FoM, for which the k2N model
is the only one managing to be accepted up to $k_{\rm
max}=0.3\,h\,\mathrm{Mpc}^{-1}$.

The HD model also provides a better recovery (compared to STD) of the
cosmological parameters, within $1$-$\sigma$ from the fiducial position, but
it fails earlier than the k2N model to provide the required goodness-of-fit
for both CMASS and LOWZ, hinting for a stronger necessity of adding
$k^{2}$-dependent stochasticity with respect to short-range non-localities.

Figure 5: Same as Figure 3 but for the combined fits of the galaxy auto power
spectrum $P_{\rm gg}$ and the galaxy-matter cross power spectrum $P_{\rm gm}$.
In all cases dark matter non-linear evolution is modelled using RESPRESSO.
Fiducial values of the linear bias and constant shot-noise contribution are
reproduced with black dashed lines. Uncertainties on the fiducial values are
marked with shaded grey bands (marking $1\mbox{-}$ and $2\mbox{-}\sigma$
confidence interval).

In terms of FoM amplitude, the k2N and HD models behave in a mostly similar
way, significantly reducing the statistical constraints of the cosmological
parameters at any given $k_{\rm max}$ cut. However, we notice how according to
the model-breaking criterion both models allow to push the analysis to higher
$k_{\rm max}$ with respect to the STD case. In other words, a proper
comparison of the merit of the three models should be done considering the
individual scale $k^{\dagger}$ after which the model breaks down. Considering
the CMASS sample, the STD model provides a FoM consistent with the one
obtained from the HD case, and $\sim 30\%$ larger than the one of the k2N
model. However, the latter two models have yet to reach the maximum scale
$k^{\dagger}$ corresponding to the breaking of the model. A simple exercise of
extending the k2N model to slightly larger modes values shows that the latter
turns out to be competitive against the STD model at a scale of $k_{\rm
max}\sim 0.35\,h\,\mathrm{Mpc}^{-1}$. As for LOWZ, the HD case is
significantly less constraining than the STD model ($\sim 50\%$), mostly
because its range of validity does not run all the way up to $k_{\rm
max}=0.30\,h\,\mathrm{Mpc}^{-1}$ but stops earlier. In contrast, the k2N model
outperforms the standard model by a $\sim 10\%$ factor, although also in this
case the model breaking scale $k^{\dagger}$ has yet to be reached.

A particular case is represented by MGS, that differently from the previous
two samples shows a stable and accurate recovery of cosmological parameters
all the way up to $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ even with the STD
model. The reason is most likely related to our choice of employing RESPRESSO
to model the non-linear matter evolution. This ultimately provides an accurate
description of $P_{\rm mm}$ even at the low redshift of MGS, leaving only the
galaxy-matter relation to be modelled. Since MGS features the lowest large-
scale bias among the samples we consider, it is likely that the one-loop
standard expansion is enough to provide both a good $\chi^{2}$ and an accurate
parameter recovery on the full range of scales up to $k_{\rm
max}=0.3\,h\,\mathrm{Mpc}^{-1}$.

Differently from the fit of the galaxy auto-power spectrum, in this case the
FoM for the three different galaxy samples do not show a flattening on the
range of scales we are considering. Checking the independent contributions
coming from the individual cosmological parameters, we notice how the FoM does
actually flatten for both $h$ and $A_{\rm s}$, whereas it monotonically
increases for $\Omega_{\rm c}h^{2}$. This is a hint showing that the matter
density parameter still benefits of additional informations contained in the
mildly non-linear regime, when combining the two galaxy power spectra.

In Figure 5 we show the dependence of the linear bias $b_{1}\sigma_{8}$ and
the constant shot-noise contribution on $k_{\rm max}$, similarly to what we
showed in Figure 3 for $P_{\rm gg}$ alone. The more stringent requirement of
simultaneously fitting $P_{\rm gg}$ and $P_{\rm gm}$ makes the standard model
fail in providing accurate measurements of both parameters. In particular,
$b_{1}\sigma_{8}$ starts to be biased at approximately the same scale for
which the model also provides incorrect measurements of the cosmological
parameters. Adding either $k^{2}$-dependent shot noise terms or higher
derivatives parameters help in recovering the true values of the bias, with a
slight preference for the former. More generally the k2N model is the only one
that is capable of simultaneously providing unbiased measurements of the
cosmological parameters, the linear bias and the large-scale stochastic
contribution to the galaxy power spectrum. This is also in agreement with the
results obtained in Eggemeier  _et al._ (2020) at fixed cosmology.

###  IV.3 Constraints on stochasticity and higher-derivative parameters

Figure 6: Marginalised posterior of the $k^{2}$-dependent stochastic
parameters $N_{2}$ and $N_{2}^{\times}$ (top), and the higher derivatives
parameters $\beta_{P}$ and $\beta_{P}^{\times}$ (bottom), as a function of
$k_{\mathrm{max}}$, for the combination of $P_{\rm gg}+P_{\rm gm}$. In all
cases dark matter non-linear evolution is modelled using RESPRESSO. For
consistency check we also show the $1$-$\sigma$ standard deviations of the
parameters $N_{2}$ and $\beta_{P}$ obtained from the fit of $P_{\rm gg}$ only
as a shaded grey area.

After having focused on the marginalised posteriors of the cosmological
parameters, linear bias and constant shot noise, here we try to put internal
constraints on both the $k^{2}$-dependent stochastic parameters and the higher
derivatives parameters. In particular, we want to check whether these terms
can be detected with statistical significance under the assumption of using a
reference effective volume like the one described in Section III.

Figure 6 shows marginalised constraints for the two sets of parameters
(extracted from the $P_{\rm gg}+P_{\rm gm}$ chains) as a function of $k_{\rm
max}$. For consistency check, in the panels referring to parameters that enter
the modeling of the galaxy auto power spectrum, $N_{2}$ and $\beta_{P}$, we
also show less stringent constraints extracted from the chains corresponding
to the fit of $P_{\rm gg}$ (shaded grey band).

The consistency of the $k^{2}$-dependent shot-noise model is additionally
reinforced by the trend observed in the top panels of Figure 6. As a matter of
fact, the relation between $N_{2}/N_{2}^{\times}$ and $k_{\rm max}$ can be
well described by a straight line with zero slope, this fact intrinsically
hinting at a correct behaviour for both parameters, as their profiles remain
stable when adding constraints at higher wave modes. Moreover, there is a
clear detection of the $k^{2}$-dependent shot-noise parameters, with a
statistical significance that further reinforce the impact of adding $N_{2}$
and $N_{2}^{\times}$ in terms of goodness-of-fit. Indeed $N_{2}$ and
$N_{2}^{\times}$ are detected at $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ with
a statistical significance of $3.5$-$\sigma$ and $2.5$-$\sigma$, respectively,
for the CMASS sample. With LOWZ the detection becomes even larger, with
numbers that are close to $6$-$\sigma$ and $4$-$\sigma$, respectively. The MGS
sample is the one showing the least significant detection of both parameters,
with a possible modification of the recovered value of $N_{2}^{\times}$ at
$k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$. However, we remind that the STD model
performs surprisingly well on this dataset, and therefore the weaker detection
of the two parameters is partially expected.

The higher derivatives parameters also show stable results as a function of
$k_{\rm max}$, with the highest detection represented by $\beta_{P}^{\times}$
for LOWZ. However on these scales the HD model is already broken, as shown in
Figure 4. For the other cases, the typical significance of the detection is
set to $\sim 2$-$\sigma$ and $\sim 4$-$\sigma$ for $\beta_{P}$ and
$\beta_{P}^{\times}$, respectively.

Differently from Eggemeier  _et al._ (2020), we do not observe incompatible
results between the marginalised posteriors of $N_{2}$ and $\beta_{P}$ from
fits of $P_{\rm gg}$ and $P_{\rm gg}+P_{\rm gm}$. However in this case we
significantly extended the dimensionality of the parameter space leading to
weaker constraints on the model parameters, so that the deviation may be
hidden by the larger statistical noise.

###  IV.4 Results for alternative matter models

Figure 7: FoM, FoB and goodness-of-fit extracted from fits of $P_{\rm
gg}+P_{\rm gm}$. Performance metrics are separated into each individual
component, i.e. $h$, $\Omega_{\rm c}h^{2}$ and $A_{\rm s}$. Different dark
matter models and different one-loop modeling assumptions are color- and
style-coded as shown in the legend. The two shaded grey regions mark the
$68$th and $95$th percentiles of the FoB and $\chi^{2}$ distribution.

So far we have fixed the modeling of the non-linear matter power spectrum
using a hybrid approach such as RESPRESSO. Another way of proceeding is to
employ a perturbative description not only for one-loop bias but also for
$P_{\rm mm}$. In this section, we analyse the impact of a different modeling
of $P_{\rm mm}$ inside Equation (10) and (11), and for this purpose we make
use of the two additional models described in Section II.2.4.

We remind the reader that these models correspond to 1) an EFT-like model
based on resummation of IR modes and non-trivial stress tensor that at leading
order in real-space can be described by a single counterterm, and 2) a RPT-
based approach with a tighter requirement of preserving the Galilean
invariance of equal time correlators (gRPT). For the three different matter
models, we show results for the STD model and for the extension including
higher-order stochastic contributions, as in the previous sections we assessed
how this extension is preferred from the three data samples we are
considering. Differently from the other matter models, the EFT-based approach
requires an additional free parameter, represented by the counterterm $c_{1}$,
that is completely degenerate with leading order higher derivatives terms. For
this reason, in this case the STD model already includes a higher-derivative
correction. The full list of configurations, together with free model
parameters, is given in Table 4.

Figure 7 shows the performance metrics for the three different dark matter
models when fitting the combination of $P_{\rm gg}$ and $P_{\rm gm}$. In this
case, we show the individual contributions to FoM and FoB from the three
cosmological parameters, $h$, $\Omega_{\rm c}h^{2}$ and $A_{\rm s}$. In
addition, for the sake of simplicity, we identify the model breaking scale
$k^{\dagger}$ with the end of the corresponding line in the FoM panels.

We first notice how the recovered goodness-of-fit is severely worsened when
adopting gRPT to describe the non-linear matter power spectrum, particularly
when considering the STD model. We separately tested the impact of extending
this model to include higher derivatives, and while this extension slightly
improves the range of validity of the gRPT-based model, it still fails in
being accepted up to $k_{\rm max}=0.30\,h\,\mathrm{Mpc}^{-1}$. The
corresponding k2N model definitely improves the recovered trends but the net
result is a worse performance respect to RESPRESSO and EFT. The STD model of
EFT breaks at higher values of $k_{\rm max}$, and this is somewhat expected
given the presence of the additional free parameter $c_{1}$. Nevertheless in
order to provide a $\chi^{2}$ consistent within $95\%$ confidence interval at
all scales we still have to include additional stochastic terms.

The FoB is overall consistent among the three different parameters, showing a
strong bias on the marginalised posterior ($>2$-$\sigma$ at $k_{\rm
max}>0.25\,h\,\mathrm{Mpc}^{-1}$) when adopting the STD for either gRPT or
RESPRESSO. Similarly to the goodness-of-fit, the STD case combined with EFT
provides better results (consistent within $1$-$\sigma$) since it already
incorporates a $k^{2}P_{\rm L}(k)$ correction which is absent from the other
matter models. The overall interpretation is that all three matter models
strongly hint for the necessity of adding $k^{2}$-dependent stochastic
contributions, both to produce an accurate description of the joint data
vector $P_{\rm gg}+P_{\rm gm}$ and to provide unbiased measurements of the
cosmological parameters.

As for the FoM, we notice how its dependency on $k_{\rm max}$ is significantly
different between the $(h,\,A_{\rm s})$ pair and $\Omega_{\rm c}h^{2}$. For
the former, the difference in performance between STD and extended models is
tiny $(\sim 10\%)$, while, for the latter, adding $k^{2}$-dependent
stochasticity (or higher derivatives) results in the suppression of the FoM at
a reference scale of $k_{\rm max}=0.2\,h\,\mathrm{Mpc}^{-1}$ by a factor of
$\sim 50\%$ and $\sim 25\%$, respectively.

Focusing on the maximum FoM achieved by each model in the range of scales we
consider, we obtain different results. The usage of gRPT typically leads to
the breaking of the model at lower wave modes, $k_{\rm max}\sim
0.2\,h\,\mathrm{Mpc}^{-1}$, when using the STD model (slightly larger with the
HD extension). Substituting gRPT with RESPRESSO helps in extending the range
of validity of the HD model, as shown in Section IV.2, while leaving almost
untouched the model breaking scale of the STD case. However, this last
configuration, i.e. RESPRESSO with the STD biasing scheme systematically
results in one of the best performing models in terms of maximally achievable
FoM. Adding $k^{2}$-dependent stochastic terms to both matter models extends
the range of validity up to $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$, with a
value of FoM which is typically consistent with the one of the STD case at the
corresponding $k^{\dagger}$. Finally, the range of validity of the EFT model
is hitting our maximum value of $0.3\,h\,\mathrm{Mpc}^{-1}$ already with the
STD model for both CMASS and MGS, while being limited to intermediate values
$(k_{\rm max}\sim 0.25\,h\,\mathrm{Mpc}^{-1})$ for LOWZ. The latter can be
also be adjusted by considering $k^{2}$-dependent shot-noise parameters.

###  IV.5 Impact of $P_{\rm gm}$ on the Figure of Merit

Figure 8: Comparison between the posterior parameter distribution ($68\%$ and
$95\%$ confidence intervals) obtained by fitting the galaxy auto power
spectrum $P_{\rm gg}$ (blue) and the combination of galaxy-galaxy and galaxy
matter power spectrum $P_{\rm gg}+P_{\rm gm}$ (red). We show constraints
obtained at $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ for the cosmological
parameters $h$, $\Omega_{\rm c}h^{2}$, $A_{\rm s}$, and the linear bias
$b_{1}\sigma_{8}$. In all cases we assume non-linear dark-matter evolution to
be described by RESPRESSO. Dashed solid lines correspond to the fiducial
values of the corresponding parameter. Grey solid bands mark the $1$-$\sigma$
and $2$-$\sigma$ error on the fiducial linear bias parameter.

After having assessed the impact of different modeling assumptions we can now
check the statistical significance of adding the galaxy-matter cross power
spectrum to a full shape measurement fit like the one we carried out in this
analysis. Indeed, as we already pointed out in Section IV.2 $3\times 2$-point
analysis will be the source of several cosmological constraints in next-
generation galaxy surveys. Although in this analysis we are excluding
redshift-space distortions, in order to separate the additional model
systematics from the ones of one-loop galaxy bias, it is anyhow interesting to
quantify the improvement in FoM when including constraints from $P_{\rm gm}$.

In Figure 8 we compare the marginalised posterior distribution of the
cosmological parameters and linear bias obtained when fitting only $P_{gg}$
and the combination $P_{\rm gg}+P_{\rm gm}$. For this comparison we model
$P_{\rm mm}$ with RESPRESSO, include $k^{2}$-dependent stochastic
contributions, and include all scales up to $k_{\rm
max}=0.3\,h\,\mathrm{Mpc}^{-1}$. We remind that the comparison is not properly
fair, as we fix the scalar amplitude $A_{\rm s}$ in the fits of the auto power
spectrum alone, while this is treated as a free parameter in the fits of the
combined statistics.

The observed trend when including $P_{\rm gm}$ in the fit is a significant
reduction of the statistical uncertainty of all parameters. In particular we
measure a $1$-$\sigma$ standard deviation for $h$ and $\Omega_{\rm c}h^{2}$
that is $1.4$, $2$, and $1.2$ times smaller for CMASS, LOWZ and MGS,
respectively. On the other hand, the linear bias absorbs most of the
additional constraining power, leading to more precise measurements of
$b_{1}\sigma_{8}$ by a factor $3$ for all three samples.

##  V Conclusions

In this work we performed a full-shape analysis of the galaxy power spectrum
meant to assess the robustness of one-loop galaxy bias models, extending the
investigation carried out in Eggemeier  _et al._ (2020) to include also the
sampling over cosmological parameters. Since we deal with galaxy clustering in
real space, we decided to fix both the baryon density $\Omega_{\rm b}\,h^{2}$
and the spectral index $n_{s}$ while leaving as free parameters the cold dark
matter density $\Omega_{\rm c}\,h^{2}$, the Hubble constant $h$, and the
scalar amplitude $A_{\rm s}$. The latter is kept as a free parameter only when
considering joint fits to the galaxy-galaxy and galaxy-matter power spectra.

We measured both observables from a set of three different synthetic galaxy
samples, whose clustering properties are meant to reproduce the ones of three
real data catalogues, i.e. CMASS and LOWZ from BOSS, and MGS from SDSS. We
rescaled the statistical uncertainties of each of these samples to match an
effective volume of $6\,\left(h^{-1}\,\mathrm{Gpc}\right)^{3}$ (i.e. the one
of the LOWZ sample), which is representative of the volume of tomographic bins
from next-generation galaxy surveys (e.g. Euclid Blanchard  _et al._ (2020)).

The analytical recipes we adopted to model these observables are based on a
standard one-loop expansion of the galaxy density field on the matter density
field, collecting terms related to both spherically-symmetric gravitational
collapse and tidal fields up to third-order in perturbations. In the first
part of the paper, we fixed the description of the non-linear matter power
spectrum using a hybrid perturbative-simulated approach, i.e. RESPRESSO, which
was already shown in Eggemeier  _et al._ (2020) to provide the closest
behaviour to the true measured matter power spectrum on the overall range of
scales and redshifts that we consider. In a later section we explicitly tested
the impact of different matter modeling assumptions, by employing an EFT-like
model and a RPT-derived model (gRPT). In all the considered cases we fixed the
quadratic tidal bias $\gamma_{2}$ using an excursion-set-derived relation, in
order to break the strong degeneracy between the latter and the cubic non-
local parameter $\gamma_{21}$.

Our main interest was in understanding the range of validity of the standard
one-loop expansion, and if our synthetic galaxy sample would hint for the
necessity of adding additional terms to the standard recipe. For this, we
considered two extensions of the standard model that take into account either
higher-derivatives of the gravitational potential or scale-dependent
stochasticity, as expected from short-range non-locality due to galaxy
formation and the halo-exclusion effect, respectively.

In order to quantify the goodness of each model we make use of a set of three
different performance metrics, which are 1) the Figure of Bias, to assess the
level of bias introduced in the recovery of cosmological parameters, 2) the
goodness-of-fit, to quantify the likelihood between the input data vector and
the best fit model, and 3) the Figure of Merit, to compare the constraining
power of each model at any given $k_{\rm max}$. All of these quantities are
easily obtained by post-processing the MCMC chains that we run for each
combination of matter modeling, one-loop bias extensions, and maximum mode
$k_{\rm max}$.

Our results can be summarised as follows:

  1. (i)

When using the fiducial description of the non-linear matter power spectrum
(i.e RESPRESSO), we find that a standard one-loop bias model (linear bias,
quadratic bias, cubic non-local bias, constant shot-noise parameter) with
fixed quadratic tidal bias can provide a good description of the galaxy power
spectrum for all our samples up to $k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$
while also returning unbiased value of the cosmological parameter set
$\left(h,\,\Omega_{\rm c}h^{2}\right)$.

  2. (ii)

Similarly to the fixed cosmology analysis of Eggemeier  _et al._ (2020) we
notice that when we consider the additional information from the galaxy-matter
cross power spectrum, the standard model is no longer capable of providing a
good performance on the overall range of scales we consider for all the three
samples. In particular, this model breaks at a scale of $k_{\rm
max}=0.25\,h\,\mathrm{Mpc}^{-1}$ and $k_{\rm max}=0.2\,h\,\mathrm{Mpc}^{-1}$
for CMASS and LOWZ, respectively, while being accepted all the way up to
$k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ for MGS. This might be representative
of the level of non-linearities in the galaxy-matter relationship, which at
first order can be captured by the parameter $b_{1}\sigma_{8}$ (that is $\sim
1.52$ for LOWZ, $\sim 1.26$ for CMASS and $\sim 1.06$ for MGS).

  3. (iii)

Extending the standard model to account also for the presence of either higher
derivatives or scale-dependent stochasticity provides a better performance
both in terms of Figure of Bias (in this case for the parameter set
$(h,\,\Omega_{\rm c}h^{2},\,A_{\rm s})$) and goodness-of-fit. Adding scale-
dependent stochastic terms restores the range of validity of the model up to
$k_{\rm max}=0.3\,h\,\mathrm{Mpc}^{-1}$ for both CMASS and LOWZ, while the
importance of higher-derivatives is limited by a worse $\chi^{2}$ for $k_{\rm
max}\gtrsim 0.2\,h\,\mathrm{Mpc}^{-1}$, and slightly worse recovery of the
linear bias parameter and the constant shot-noise correction. As an additional
check, we verify that the marginalised posterior distribution of both
$k^{2}$-dependent terms and higher derivatives parameters is stable as a
function of $k_{\rm max}$. However we notice how the maximally achievable FoM
of the extended models is consistent with the one of the standard
configuration at its lower $k_{\rm max}$, leading to equivalent statistical
constraints on the cosmological parameters.

  4. (iv)

Changing the description of the non-linear matter power spectrum to either EFT
or gRPT induces modifications to the range of validity of one-loop galaxy bias
models. Overall, we find that, as highlighted in Eggemeier  _et al._ (2020),
RESPRESSO is the best performing model, followed by EFT and gRPT. The EFT-
based model is penalised in terms of Figure of Merit because of the presence
of an additional degree of freedom (two in case of the combined fits), but it
compensates this with an extended range of validity. On the contrary, gRPT
breaks sooner than the other models, at a typical scale of $k_{\rm max}\sim
0.2\,h\,\mathrm{Mpc}^{-1}$, which is consistent with the range of validity
observed in direct fits of the matter power spectrum. This indicates that this
model would also benefit of the presence of an additional free parameter
representing a non-zero stress tensor such as in the case of the EFT-based
model, as expected. Although the gRPT-based model breaks earlier than the
other two models, its maximally achievable FoM is consistent with the one
recovered using RESPRESSO.

  5. (v)

The additional constraints coming from the galaxy-matter cross power spectrum
results in an improved statistical precision on measurements of cosmological
parameters. Although we fix the scalar amplitude $A_{\rm s}$ when fitting only
the galaxy auto power spectrum, we still find that in the combined case we
achieve better constraints on both $h$ and $\Omega_{\rm c}h^{2}$ by a factor
$1.4$, $2$ and $1.2$ for CMASS, LOWZ and MGS, respectively. At the same time
we notice how a significant fraction of the additional constraining power can
be absorbed by nuisance parameters, such as the linear bias, for which
constraints are tighter by a factor of $3$ for all the three samples.

In this work we analyse three samples with an effective volume which is
rescaled to a value of $6\,h^{-3}\,\mathrm{Gpc}^{3}$. Although the results we
obtained are partially based on this choice, as shown in Eggemeier  _et al._
(2020), this volume can be well representative of individual tomographic
redshift bins that will be adopted by next-generation galaxy surveys. In
addition we explore the additional constraining power brought by adding the
galaxy-mass cross power spectrum to galaxy power spectrum fits, which is a
relevant study in the context of $3\times 2$pt analyses of next-generation
imaging surveys. Therefore, these results can provide useful insights when
adopting one-loop perturbation theory to describe the relationship between the
galaxy and matter density fields in that context Pandey  _et al._ (2020). In
order to completely concentrate on galaxy bias, we performed this analysis
using real-space coordinates, avoiding a further modeling layer for redshift-
space distortions. We leave to future works the exploration of this effect,
both using the power spectrum alone, and combining it with the galaxy
bispectrum.

###### Acknowledgements.

AP and MC acknowledge support from the Spanish Ministry of Science and
Innovation through grants PGC2018-102021-B-100 and ESP2017-89838-C3-1-R, and
EU grants LACEGAL 734374 and EWC 776247 with ERDF funds. IEEC is funded by the
CERCA program of the Generalitat de Catalunya. AE acknowledges support from
the European Research Council (grant num- ber ERC-StG-716532-PUNCA. AGS
acknowledges support from the Excellence Cluster ORIGINS, which is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311 AP and AGS would like to
thank Daniel Farrow, Jiamin Hou, Martha Lippich and Agne Semenaite for their
help and useful discussions. The analysis presented here has been performed on
the high-performance computing resources of the Max Planck Computing and Data
Facility (MPCDF) in Garching. This research made use of matplotlib, a Python
library for publication quality graphics Hunter (2007).

## References

  * Davis and Peebles (1983) M. Davis and P. J. E. Peebles, ApJ 267, 465 (1983).
  * Maddox  _et al._ (1990) S. J. Maddox, G. Efstathiou, W. J. Sutherland, and J. Loveday, MNRAS 242, 43 (1990).
  * Tegmark  _et al._ (2004) M. Tegmark, M. R. Blanton, M. A. Strauss, F. Hoyle, D. Schlegel, R. Scoccimarro, M. S. Vogeley, D. H. Weinberg, I. Zehavi, A. Berlind, and et al., The Astrophysical Journal 606, 702–740 (2004).
  * Cole  _et al._ (2005) S. Cole, W. J. Percival, J. A. Peacock, P. Norberg, C. M. Baugh, C. S. Frenk, I. Baldry, J. Bland-Hawthorn, T. Bridges, R. Cannon, and et al., Monthly Notices of the Royal Astronomical Society 362, 505–534 (2005).
  * Hauser and Peebles (1973) M. G. Hauser and P. J. E. Peebles, ApJ 185, 757 (1973).
  * Kaiser (1984) N. Kaiser, ApJ 284, L9 (1984).
  * Davis  _et al._ (1985) M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, ApJ 292, 371 (1985).
  * Bardeen  _et al._ (1986) J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, ApJ 304, 15 (1986).
  * Tegmark  _et al._ (2004) M. Tegmark, M. Blanton, M. Strauss, F. Hoyle, D. Schlegel, R. Scoccimarro, M. Vogeley, D. Weinberg, I. Zehavi, A. Berlind, T. Budavari, A. Connolly, D. Eisenstein, D. Finkbeiner, J. Frieman, J. Gunn, A. Hamilton, L. Hui, B. Jain, D. Johnston, S. Kent, H. Lin, R. Nakajima, R. Nichol, J. Ostriker, A. Pope, R. Scranton, U. Seljak, R. Sheth, A. Stebbins, A. Szalay, I. Szapudi, L. Verde, Y. Xu, J. Annis, N. Bahcall, J. Brinkmann, S. Burles, F. Castander, I. Csabai, J. Loveday, M. Doi, M. Fukugita, J. Gott, G. Hennessy, D. Hogg, Ž. Ivezić, G. Knapp, D. Lamb, B. Lee, R. Lupton, T. McKay, P. Kunszt, J. Munn, L. O’Connell, J. Peoples, J. Pier, M. Richmond, C. Rockosi, D. Schneider, C. Stoughton, D. Tucker, D. Vanden Berk, B. Yanny, and D. York, ApJ 606, 702 (2004).
  * Zehavi  _et al._ (2005) I. Zehavi, Z. Zheng, D. H. Weinberg, J. A. Frieman, A. A. Berlind, M. R. Blanton, R. Scoccimarro, R. K. Sheth, M. A. Strauss, I. Kayo, Y. Suto, M. Fukugita, O. Nakamura, N. A. Bahcall, J. Brinkmann, J. E. Gunn, G. S. Hennessy, Ž. Ivezić, G. R. Knapp, J. Loveday, A. Meiksin, D. J. Schlegel, D. P. Schneider, I. Szapudi, M. Tegmark, M. S. Vogeley, and D. G. York, ApJ 630, 1 (2005).
  * Desjacques  _et al._ (2018) V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rep. 733, 1 (2018).
  * Mo and White (1996) H. J. Mo and S. D. M. White, MNRAS 282, 347 (1996), arXiv:astro-ph/9512127 .
  * Mo  _et al._ (1997) H. J. Mo, Y. P. Jing, and S. D. M. White, MNRAS 284, 189 (1997).
  * Frieman and Gaztanaga (1994) J. Frieman and E. Gaztanaga, ApJ 425, 392 (1994).
  * Chan  _et al._ (2012) K. C. Chan, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 85, 083509 (2012), arXiv:1201.3614 [astro-ph.CO] .
  * Baldauf  _et al._ (2012) T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Phys. Rev. D 86, 083540 (2012), arXiv:1201.4827 [astro-ph.CO] .
  * Manera and Gaztañaga (2011) M. Manera and E. Gaztañaga, MNRAS 415, 383 (2011), arXiv:0912.0446 [astro-ph.CO] .
  * Roth and Porciani (2011) N. Roth and C. Porciani, MNRAS 415, 829 (2011), arXiv:1101.1520 [astro-ph.CO] .
  * Pollack  _et al._ (2012) J. E. Pollack, R. E. Smith, and C. Porciani, MNRAS 420, 3469 (2012), arXiv:1109.3458 .
  * Pollack  _et al._ (2014) J. E. Pollack, R. E. Smith, and C. Porciani, MNRAS 440, 555 (2014), arXiv:1309.0504 [astro-ph.CO] .
  * McDonald and Roy (2009) P. McDonald and A. Roy, J. Cosmology Astropart. Phys 8, 020 (2009), arXiv:0902.0991 [astro-ph.CO] .
  * Desjacques (2008) V. Desjacques, Phys. Rev. D 78, 103503 (2008), arXiv:0806.0007 [astro-ph] .
  * Desjacques  _et al._ (2010) V. Desjacques, M. Crocce, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 82, 103529 (2010), arXiv:1009.3449 [astro-ph.CO] .
  * Scherrer and Weinberg (1998) R. J. Scherrer and D. H. Weinberg, ApJ 504, 607 (1998), arXiv:astro-ph/9712192 [astro-ph] .
  * Dekel and Lahav (1999) A. Dekel and O. Lahav, ApJ 520, 24 (1999), arXiv:astro-ph/9806193 [astro-ph] .
  * Taruya and Soda (1999) A. Taruya and J. Soda, ApJ 522, 46 (1999), arXiv:astro-ph/9809204 [astro-ph] .
  * Matsubara (1999) T. Matsubara, ApJ 525, 543 (1999), arXiv:astro-ph/9906029 [astro-ph] .
  * Sheth and Lemson (1999) R. Sheth and G. Lemson, MNRAS 304, 767 (1999).
  * Smith  _et al._ (2007) R. E. Smith, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 75, 063512 (2007), arXiv:astro-ph/0609547 .
  * Baldauf  _et al._ (2013) T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Phys. Rev. D 88, 083507 (2013), arXiv:1305.2917 [astro-ph.CO] .
  * Eggemeier  _et al._ (2020) A. Eggemeier, R. Scoccimarro, M. Crocce, A. Pezzotta, and A. G. Sánchez, Phys. Rev. D 102, 103530 (2020).
  * Sánchez  _et al._ (2017) A. G. Sánchez, R. Scoccimarro, M. Crocce, J. N. Grieb, S. Salazar-Albornoz, C. Dalla Vecchia, M. Lippich, F. Beutler, J. R. Brownstein, C.-H. Chuang, D. J. Eisenstein, F.-S. Kitaura, M. D. Olmstead, W. J. Percival, F. Prada, S. Rodríguez-Torres, A. J. Ross, L. Samushia, H.-J. Seo, J. Tinker, R. Tojeiro, M. Vargas-Magaña, Y. Wang, and G.-B. Zhao, MNRAS 464, 1640 (2017), arXiv:1607.03147 .
  * Gil-Marín  _et al._ (2017) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, MNRAS 465, 1757 (2017), arXiv:1606.00439 [astro-ph.CO] .
  * Beutler  _et al._ (2017) F. Beutler, H.-J. Seo, S. Saito, C.-H. Chuang, A. J. Cuesta, D. J. Eisenstein, H. Gil-Marín, J. N. Grieb, N. Hand, F.-S. Kitaura, C. Modi, R. C. Nichol, M. D. Olmstead, W. J. Percival, F. Prada, A. G. Sánchez, S. Rodriguez-Torres, A. J. Ross, N. P. Ross, D. P. Schneider, J. Tinker, R. Tojeiro, and M. Vargas-Magaña, MNRAS 466, 2242 (2017), arXiv:1607.03150 [astro-ph.CO] .
  * Grieb  _et al._ (2017) J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz, R. Scoccimarro, M. Crocce, C. Dalla Vecchia, F. Montesano, H. Gil-Marín, A. J. Ross, F. Beutler, S. Rodríguez-Torres, C.-H. Chuang, F. Prada, F.-S. Kitaura, A. J. Cuesta, D. J. Eisenstein, W. J. Percival, M. Vargas-Magaña, J. L. Tinker, R. Tojeiro, J. R. Brownstein, C. Maraston, R. C. Nichol, M. D. Olmstead, L. Samushia, H.-J. Seo, A. Streblyanska, and G.-b. Zhao, MNRAS 467, 2085 (2017), arXiv:1607.03143 .
  * Ivanov  _et al._ (2019) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, arXiv e-prints , arXiv:1909.05277 (2019), arXiv:1909.05277 [astro-ph.CO] .
  * D’Amico  _et al._ (2019) G. D’Amico, J. Gleyzes, N. Kokron, D. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, arXiv e-prints , arXiv:1909.05271 (2019), arXiv:1909.05271 [astro-ph.CO] .
  * Tröster  _et al._ (2020) T. Tröster, A. G. Sánchez, M. Asgari, C. Blake, M. Crocce, C. Heymans, H. Hildebrandt, B. Joachimi, S. Joudaki, A. Kannawadi, C.-A. Lin, and A. Wright, A&A 633, L10 (2020), arXiv:1909.11006 [astro-ph.CO] .
  * Eggemeier  _et al._ (2021) A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta, and A. G. Sánchez, “Testing one-loop galaxy bias: joint analysis of power spectrum and bispectrum,” (2021), arXiv:2102.06902 [astro-ph.CO] .
  * Eisenstein  _et al._ (2011) D. J. Eisenstein, D. H. Weinberg, E. Agol, H. Aihara, C. Allende Prieto, S. F. Anderson, J. A. Arns, E. Aubourg, S. Bailey, E. Balbinot, and et al., The Astronomical Journal 142, 72 (2011).
  * Dawson  _et al._ (2013) K. S. Dawson, D. J. Schlegel, and et al., AJ 145, 10 (2013), arXiv:1208.0022 [astro-ph.CO] .
  * Reid  _et al._ (2015) B. Reid, S. Ho, N. Padmanabhan, W. J. Percival, J. Tinker, R. Tojeiro, M. White, D. J. Eisenstein, C. Maraston, and et al., (2015), arXiv:1509.06529 [astro-ph.CO] .
  * Strauss  _et al._ (2002) M. A. Strauss, D. H. Weinberg, R. H. Lupton, V. K. Narayanan, and et al., AJ 124, 1810 (2002), arXiv:astro-ph/0206225 [astro-ph] .
  * Eggemeier  _et al._ (2019) A. Eggemeier, R. Scoccimarro, and R. E. Smith, Phys. Rev. D 99, 123514 (2019), arXiv:1812.03208 [astro-ph.CO] .
  * Chan  _et al._ (2017) K. C. Chan, R. K. Sheth, and R. Scoccimarro, Phys. Rev. D 96, 103543 (2017), arXiv:1511.01909 .
  * Musso and Sheth (2012) M. Musso and R. K. Sheth, MNRAS 423, L102 (2012), arXiv:1201.3876 [astro-ph.CO] .
  * Lazeyras and Schmidt (2019) T. Lazeyras and F. Schmidt, J. Cosmology Astropart. Phys 2019, 041 (2019), arXiv:1904.11294 [astro-ph.CO] .
  * McDonald (2006) P. McDonald, Phys. Rev. D 74, 103512 (2006), astro-ph/0609413 .
  * Assassi  _et al._ (2014) V. Assassi, D. Baumann, D. Green, and M. Zaldarriaga, J. Cosmology Astropart. Phys 8, 056 (2014), arXiv:1402.5916 .
  * Schmidt  _et al._ (2013) F. Schmidt, D. Jeong, and V. Desjacques, Phys. Rev. D 88, 023515 (2013), arXiv:1212.0868 [astro-ph.CO] .
  * Crocce and Scoccimarro (2006a) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063519 (2006a), arXiv:astro-ph/0509418 .
  * Bernardeau  _et al._ (2002) F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. 367, 1 (2002), astro-ph/0112551 .
  * Grieb  _et al._ (2016) J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz, and C. Dalla Vecchia, MNRAS 457, 1577 (2016), arXiv:1509.04293 [astro-ph.CO] .
  * Blas  _et al._ (2013) D. Blas, M. Garny, and T. Konstandin, ArXiv e-prints (2013), arXiv:1309.3308 [astro-ph.CO] .
  * Crocce and Scoccimarro (2006b) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063520 (2006b), arXiv:astro-ph/0509419 .
  * Eisenstein  _et al._ (2007a) D. J. Eisenstein, H.-J. Seo, and M. White, ApJ 664, 660 (2007a), arXiv:astro-ph/0604361 .
  * Eisenstein  _et al._ (2007b) D. J. Eisenstein, H.-J. Seo, E. Sirko, and D. N. Spergel, ApJ 664, 675 (2007b), arXiv:astro-ph/0604362 [astro-ph] .
  * Matarrese and Pietroni (2008) S. Matarrese and M. Pietroni, Modern Physics Letters A 23, 25 (2008), arXiv:astro-ph/0702653 .
  * Crocce and Scoccimarro (2008) M. Crocce and R. Scoccimarro, Phys. Rev. D 77, 023533 (2008), arXiv:0704.2783 [astro-ph] .
  * Seo  _et al._ (2008) H.-J. Seo, E. R. Siegel, D. J. Eisenstein, and M. White, ApJ 686, 13 (2008), arXiv:0805.0117 .
  * Vlah  _et al._ (2016) Z. Vlah, U. Seljak, M. Yat Chu, and Y. Feng, J. Cosmology Astropart. Phys 2016, 057 (2016), arXiv:1509.02120 [astro-ph.CO] .
  * Osato  _et al._ (2019) K. Osato, T. Nishimichi, F. Bernardeau, and A. Taruya, Phys. Rev. D 99, 063530 (2019), arXiv:1810.10104 [astro-ph.CO] .
  * Eisenstein and Hu (1999) D. J. Eisenstein and W. Hu, The Astrophysical Journal 511, 5 (1999).
  * Blas  _et al._ (2016) D. Blas, M. Garny, M. M. Ivanov, and S. Sibiryakov, J. Cosmology Astropart. Phys 2016, 028 (2016), arXiv:1605.02149 [astro-ph.CO] .
  * Pueblas and Scoccimarro (2009) S. Pueblas and R. Scoccimarro, Phys. Rev. D 80, 043504 (2009), arXiv:0809.4606 .
  * Carrasco  _et al._ (2012) J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore, ArXiv e-prints (2012), arXiv:1206.2926 [astro-ph.CO] .
  * Baumann  _et al._ (2012) D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, J. Cosmology Astropart. Phys 7, 051 (2012), arXiv:1004.2488 [astro-ph.CO] .
  * Carrasco  _et al._ (2013) J. J. M. Carrasco, S. Foreman, D. Green, and L. Senatore, ArXiv e-prints (2013), arXiv:1310.0464 [astro-ph.CO] .
  * Baldauf  _et al._ (2015) T. Baldauf, L. Mercolli, and M. Zaldarriaga, Phys. Rev. D 92, 123007 (2015), arXiv:1507.02256 [astro-ph.CO] .
  * Bernardeau  _et al._ (2014) F. Bernardeau, A. Taruya, and T. Nishimichi, Phys. Rev. D 89, 023502 (2014).
  * Nishimichi  _et al._ (2016) T. Nishimichi, F. Bernardeau, and A. Taruya, Physics Letters B 762, 247 (2016), arXiv:1411.2970 [astro-ph.CO] .
  * Nishimichi  _et al._ (2017) T. Nishimichi, F. Bernardeau, and A. Taruya, Phys. Rev. D 96, 123515 (2017), arXiv:1708.08946 [astro-ph.CO] .
  * Planck Collaboration (2016) Planck Collaboration, A&A 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Angulo  _et al._ (2015) R. Angulo, M. Fasiello, L. Senatore, and Z. Vlah, Journal of Cosmology and Astro-Particle Physics 2015, 029 (2015), arXiv:1503.08826 [astro-ph.CO] .
  * Bernardeau  _et al._ (2008) F. Bernardeau, M. Crocce, and R. Scoccimarro, Phys. Rev. D 78, 103521 (2008), arXiv:0806.2334 .
  * Taruya  _et al._ (2013) A. Taruya, T. Nishimichi, and F. Bernardeau, Phys. Rev. D 87, 083509 (2013), arXiv:1301.3624 [astro-ph.CO] .
  * Scoccimarro and Frieman (1996) R. Scoccimarro and J. Frieman, ApJS 105, 37 (1996), arXiv:astro-ph/9509047 .
  * Casas-Miranda  _et al._ (2002) R. Casas-Miranda, H. J. Mo, R. K. Sheth, and G. Boerner, MNRAS 333, 730 (2002), arXiv:astro-ph/0105008 [astro-ph] .
  * Akbar Abolhasani  _et al._ (2016) A. Akbar Abolhasani, M. Mirbabayi, and E. Pajer, J. Cosmology Astropart. Phys 2016, 063 (2016), arXiv:1509.07886 [hep-th] .
  * Catelan  _et al._ (1998) P. Catelan, F. Lucchin, S. Matarrese, and C. Porciani, MNRAS 297, 692 (1998), arXiv:astro-ph/9708067 .
  * Catelan  _et al._ (2000) P. Catelan, C. Porciani, and M. Kamionkowski, MNRAS 318, L39 (2000), astro-ph/0005544 .
  * Matsubara (2011) T. Matsubara, Phys. Rev. D 83, 083518 (2011), arXiv:1102.4619 [astro-ph.CO] .
  * Lazeyras and Schmidt (2018) T. Lazeyras and F. Schmidt, Journal of Cosmology and Astro-Particle Physics 2018, 008 (2018), arXiv:1712.07531 [astro-ph.CO] .
  * Abidi and Baldauf (2018) M. M. Abidi and T. Baldauf, J. Cosmology Astropart. Phys 7, 029 (2018), arXiv:1802.07622 .
  * Sheth  _et al._ (2013) R. K. Sheth, K. C. Chan, and R. Scoccimarro, Phys. Rev. D 87, 083002 (2013), arXiv:1207.7117 [astro-ph.CO] .
  * Springel  _et al._ (2001) V. Springel, N. Yoshida, and S. D. White, New Astronomy 6, 79–117 (2001).
  * Springel (2005) V. Springel, Monthly Notices of the Royal Astronomical Society 364, 1105–1134 (2005).
  * Lewis and Bridle (2002) A. Lewis and S. Bridle, Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 [astro-ph] .
  * Crocce  _et al._ (2006) M. Crocce, S. Pueblas, and R. Scoccimarro, MNRAS 373, 369 (2006), arXiv:astro-ph/0606505 .
  * Alam  _et al._ (2017) S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, A. S. Bolton, J. R. Brownstein, A. Burden, C.-H. Chuang, J. Comparat, A. J. Cuesta, K. S. Dawson, D. J. Eisenstein, S. Escoffier, H. Gil-Marín, J. N. Grieb, N. Hand, S. Ho, K. Kinemuchi, D. Kirkby, F. Kitaura, E. Malanushenko, V. Malanushenko, C. Maraston, C. K. McBride, R. C. Nichol, M. D. Olmstead, D. Oravetz, N. Padmanabhan, N. Palanque-Delabrouille, K. Pan, M. Pellejero-Ibanez, W. J. Percival, P. Petitjean, F. Prada, A. M. Price-Whelan, B. A. Reid, S. A. Rodríguez-Torres, N. A. Roe, A. J. Ross, N. P. Ross, G. Rossi, J. A. Rubiño-Martín, S. Saito, S. Salazar-Albornoz, L. Samushia, A. G. Sánchez, S. Satpathy, D. J. Schlegel, D. P. Schneider, C. G. Scóccola, H.-J. Seo, E. S. Sheldon, A. Simmons, A. Slosar, M. A. Strauss, M. E. C. Swanson, D. Thomas, J. L. Tinker, R. Tojeiro, M. V. Magaña, J. A. Vazquez, L. Verde, D. A. Wake, Y. Wang, D. H. Weinberg, M. White, W. M. Wood-Vasey, C. Yèche, I. Zehavi, Z. Zhai, and G.-B. Zhao, MNRAS 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * McBride  _et al._ (2009) C. McBride, A. Berlind, R. Scoccimarro, R. Wechsler, M. Busha, J. Gardner, and F. van den Bosch, in  _American Astronomical Society Meeting Abstracts #213_ , American Astronomical Society Meeting Abstracts, Vol. 213 (2009) p. 425.06.
  * Sinha  _et al._ (2018) M. Sinha, A. A. Berlind, C. K. McBride, R. Scoccimarro, J. A. Piscionere, and B. D. Wibking, MNRAS 478, 1042 (2018), arXiv:1708.04892 [astro-ph.CO] .
  * Seljak and Zaldarriaga (1996) U. Seljak and M. Zaldarriaga, ApJ 469, 437 (1996).
  * Sefusatti  _et al._ (2016) E. Sefusatti, M. Crocce, R. Scoccimarro, and H. M. P. Couchman, MNRAS 460, 3624 (2016), arXiv:1512.07295 [astro-ph.CO] .
  * Feldman  _et al._ (1994) H. A. Feldman, N. Kaiser, and J. A. Peacock, ApJ 426, 23 (1994), arXiv:astro-ph/9304022 .
  * Tegmark (1997) M. Tegmark, Phys. Rev. Lett. 79, 3806 (1997), arXiv:astro-ph/9706198 [astro-ph] .
  * Oddo  _et al._ (2020) A. Oddo, E. Sefusatti, C. Porciani, P. Monaco, and A. G. Sánchez, J. Cosmology Astropart. Phys 2020, 056 (2020), arXiv:1908.01774 [astro-ph.CO] .
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statistical Science 7, 457 (1992).
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Sánchez (2020) A. G. Sánchez, Phys. Rev. D 102, 123511 (2020).
  * Abbott  _et al._ (2018) T. Abbott, F. Abdalla, A. Alarcon, J. Aleksić, S. Allam, S. Allen, A. Amara, J. Annis, J. Asorey, S. Avila, and et al., Physical Review D 98 (2018), 10.1103/physrevd.98.043526.
  * Heymans  _et al._ (2020) C. Heymans, T. Tröster, M. Asgari, C. Blake, H. Hildebrandt, B. Joachimi, K. Kuijken, C.-A. Lin, A. G. Sánchez, J. L. van den Busch, A. H. Wright, A. Amon, M. Bilicki, J. de Jong, M. Crocce, A. Dvornik, T. Erben, M. C. Fortuna, F. Getman, B. Giblin, K. Glazebrook, H. Hoekstra, S. Joudaki, A. Kannawadi, F. Köhlinger, C. Lidman, L. Miller, N. R. Napolitano, D. Parkinson, P. Schneider, H. Shan, E. Valentijn, G. V. Kleijn, and C. Wolf, (2020), arXiv:2007.15632 [astro-ph.CO] .
  * Pandey  _et al._ (2020) S. Pandey, E. Krause, B. Jain, N. MacCrann, J. Blazek, M. Crocce, J. DeRose, X. Fang, I. Ferrero, O. Friedrich, M. Aguena, S. Allam, J. Annis, S. Avila, G. M. Bernstein, D. Brooks, D. L. Burke, A. Carnero Rosell, M. Carrasco Kind, J. Carretero, M. Costanzi, L. N. da Costa, J. De Vicente, S. Desai, J. Elvin-Poole, S. Everett, P. Fosalba, J. Frieman, J. García-Bellido, D. Gruen, R. A. Gruendl, J. Gschwend, G. Gutierrez, K. Honscheid, K. Kuehn, N. Kuropatkin, M. A. G. Maia, J. L. Marshall, F. Menanteau, R. Miquel, A. Palmese, F. Paz-Chinchón, A. A. Plazas, A. Roodman, E. Sanchez, V. Scarpine, M. Schubnell, S. Serrano, I. Sevilla-Noarbe, M. Smith, M. Soares-Santos, E. Suchyta, M. E. C. Swanson, G. Tarle, and J. Weller (DES Collaboration), Phys. Rev. D 102, 123522 (2020).
  * Blanchard  _et al._ (2020) A. Blanchard, S. Camera, C. Carbone, V. F. Cardone, S. Casas, S. Clesse, S. Ilić, M. Kilbinger, T. Kitching, and et al., A&A 642, A191 (2020).
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).