listed in Table 4.

Figure 5: Comparison between the marginalised constraints on the linear bias
parameter $b_{1}$ and the shot-noise parameter $\alpha_{P,1}$ obtained at
fixed cosmology, and the fiducial values listed in Table 4 obtained using only
the large scale-limit of Eq. (43). The first two and last two rows show
results for the Model 1 and Model 3 HOD samples, respectively. In both cases,
the upper panels show constraints on the linear bias $b_{1}$, while the bottom
ones show constraints on the constant shot-noise parameter $\alpha_{P,1}$.
Different colors correspond to different assumptions on the total number of
free bias parameters, as shown in the legend. Star symbols highlight the
position of the maximum-likelihood for the case with all bias parameters free
to vary. Dashed grey lines and shaded bands mark the fiducial value and
$1\sigma$ confidence interval from Table 4. Figure 6: Marginalised 2d
constraints on the pair $\left(b_{1},\alpha_{P,1}\right)$ for the case
corresponding to the Model 3 HOD sample at $z=1.2$. The three different panels
correspond to the case in which all the parameters are left free (left panel),
$b_{{\cal G}_{2}}$ is fixed to the excursion-set relation as in Eq. (33)
(middle panel), and $b_{\Gamma_{3}}$ is fixed to the coevolution relation as
in Eq. (34) (right panel). Different colours represent different $k_{\rm max}$
values, as listed in the legend. Fiducial $68\%$ confidence intervals for both
$b_{1}$ and $\alpha_{P,1}$ are shown with grey bands.

In Fig. 5 we show the marginalised constraints obtained fitting the full model
of Eq. (32) to the eight data vectors, against the fiducial values of
$\left\\{b_{1},\alpha_{P,1}\right\\}$ obtained from the large-scale limit as
in Eq. (43). We test four different model configurations, which differ by the
total number of bias parameters that are kept fixed to the relations presented
in Sect. 3.1.3.

  * (i)

All nuisance parameters are left free to vary while sampling the posterior
distribution, for a total of six free parameters – linear bias $b_{1}$, local
quadratic bias $b_{2}$, non-local quadratic bias $b_{{\cal G}_{2}}$, non-local
cubic bias $b_{\Gamma_{3}}$, matter counterterm and higher-derivative bias
$c_{0}$, and constant shot-noise parameter $\alpha_{P,1}$ .

  * (ii)

$b_{{\cal G}_{2}}$ is fixed to the excursion-set-based relation defined in Eq.
(35), for a total number of five parameters.

  * (iii)

$b_{\Gamma_{3}}$ is fixed to the coevolution relation defined in Eq. (34), for
a total number of five parameters.

  * (iv)

Both $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$ are fixed to the relations
assumed in (ii) and (iii), respectively, for a total number of four
parameters.

In all these cases, we keep the scale-dependent shot-noise parameter
$\alpha_{P,2}$ fixed to zero (we test the validity of this assumption in Sect.
5.3).

Overall, the most relaxed configuration – case (i) in the previous list, shown
with blue points and errorbars in Fig. 5 – with all the bias parameters free
to vary is capable of capturing the correct amplitude of both $b_{1}$ and
$\alpha_{P,1}$ for the majority of the tested $k_{\rm max}$ values and
redshifts, showing a mild running of the one-dimensional marginalised values
that becomes relevant only for the lowest redshift snapshot we consider, on
scales $k_{\rm max}>0.2\,h\,{\rm Mpc}^{-1}$. The same effect is partially
present for the $z=1.2$ snapshot, although less significant: as a matter of
fact, the marginalised constraints are consistent with their fiducial values
at better than $2\sigma$. Rather than only considering the mean posterior
distribution, it is instructive to also plot the maximum-likelihood point in
the parameter space under consideration. We estimate this quantity using as a
proxy the point in the sampled posterior distribution that maximises the
likelihood, even though the latter is partially affected by a certain degree
of stochasticity. In this case (star symbols in Fig. 5) we observe a shift
towards the fiducial values, even if not for all configurations. A discrepancy
between the maximum-likelihood point and the mean of the marginalised
posterior is a clear hint at the presence of projection effects, also known as
prior volume effects, due to the high dimensionality of the parameter space
and to nonlinear degeneracies among the model parameters.

Fixing either $b_{{\cal G}_{2}}$ – case (ii), orange points and errorbars – or
$b_{\Gamma_{3}}$ – case (iii), green points and errorbars – does not
significantly help in terms of accuracy of the marginalised constraints, with
systematic deviations that can still become larger than the $1\sigma$
confidence interval. However we find that, while fixing $b_{\Gamma_{3}}$
typically results in similar constraining power on both $b_{1}$ and
$\alpha_{P,1}$, imposing a relation on $b_{{\cal G}_{2}}$ leads to definitely
tighter posteriors. This is the result of breaking the strong degeneracy
between the two non-local bias parameters, and at the same time the one
between the quadratic biases, $b_{2}$ and $b_{{\cal G}_{2}}$, leaving the
remaining parameters to be more tightly constrained (a clear example of these
is displayed in the right panel of Fig. 17). The same clearly happens when
combining the two previous relations – case (iv), red points and errorbars –
since with this setup we completely break the degeneracies in the considered
parameter space. However, in this case we observe a deviation from the
fiducial values of $b_{1}$ which can reach more than $2\sigma$ for some of the
configurations, in particular at low redshift, hinting at a departure from the
assumption of conserved evolution.

The effect of the strong $b_{2}$-$b_{{\cal G}_{2}}$-$b_{\Gamma_{3}}$
degeneracy can be observed in a more direct way by inspection of the 2d
marginalised constraints in the $b_{1}$-$\alpha_{P,1}$ subspace. In Fig. 6 we
show such posterior distributions, taking as a reference the Model 3 sample at
$z=1.2$. The different panels correspond to different bias relations, from
left to right: the case with all the parameters free to vary, with $b_{{\cal
G}_{2}}$ fixed to the excursion-set-based relation, and with $b_{\Gamma_{3}}$
fixed to the coevolution relation. A first consideration to make is that there
is a non-trivial degeneracy between the two parameters, for which projection
effects might bias the 1d constraints without necessarily meaning that the
hyper-dimensional posterior distribution does not cover the fiducial values of
the parameters. Secondly, we can observe how the case with fixed $b_{{\cal
G}_{2}}$ gives the tightest constraints for both parameters, with an increase
in the merit of the constraints that is directly related to the maximum scale
adopted in the fit, up to $k_{\rm max}=0.4\,h\,{\rm Mpc}^{-1}$. Once more,
this trend can be easily explained by the effective breaking of the degeneracy
among the higher-order bias parameters.

For a limited number of cases, even when fixing one or more degrees of
freedom, we find that the final posterior distribution can still appear multi-
modal, leading to enlarged constraints when marginalising over the remaining
parameters. For this reason, some of the chains where both tidal bias
parameters $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$ are fixed feature
marginalised constraints that are larger than the ones with one additional
degree of freedom. This is clearly noticeable for the largest $k_{\rm max}$
bin of the Model 3 sample at $z=0.9$.

Figure 7: Number of effective parameters that can be properly constrained by
the input data vectors of the Model 1 (circles) and 3 (squares) HOD samples,
respectively, as a function of the maximum wave mode included in the analysis
and for different configurations of the bias model, as listed in the legend.
Different panels correspond to different redshifts, as shown in the
corresponding top left corner.

The self-consistency of the different models in terms of the number of model
parameters can be assessed using an additional statistics. In this context, we
are interested in determining the total number of parameters that can be
effectively constrained by the data vectors. A commonly employed statistics is
represented by the $p_{\rm V}$ value, defined as (Gelman et al. 2014)

$p_{\rm
V}=\frac{1}{2}\left\langle\left(\,\chi^{\,2}-\left\langle\,\chi^{\,2}\right\rangle\right)^{2}\right\rangle\,,$
(44)

that is, the variance of the corresponding $\chi^{\,2}$ distribution. This
number indirectly tracks the presence of degeneracies among the model
parameters, and only converges to the total number of free parameters for a
normal distribution. In order for a theory model to effectively constrain a
given number of parameters, the $p_{V}$ value is expected to reach that same
value, and can therefore be used as a proxy for the self-consistency of
different model configurations. In Fig. 7 we show this value as measured from
both sets of HOD samples and for different values of $k_{\rm max}$. In
practice, we notice how the model with all parameters free never reaches the
expected value of $p_{V}=6$, even for the largest value of $k_{\rm max}$, with
the exception of a couple of configurations. This shows how a six-parameter
model is most likely resulting in overfitting. On the contrary, fixing one of
the two tidal biases makes the $p_{V}$ reach the expected limit above some
$k_{\rm max}$, with a transition that typically happens sooner for the case
with fixed $b_{{\cal G}_{2}}$. This reinforces the conclusion that this
configuration is preferred with respect to the others under consideration.

We note that these results may be partially affected by the presence of cosmic
variance in the data vectors. For this reason, in Appendix F we explicitly
assess the impact of this extra contribution, using both a smooth and a noisy
realization of the data, generated using the theory code. This test shows how
most of the residuals observed in Fig. 5 can be explained by sample variance
affecting our data vectors.

Figure 8: Marginalised 1d constraints on the scale-dependent shot-noise
parameter $\alpha_{P,2}$ in the fits with fixed cosmological parameters.
Different rows correspond to different redshifts (top to bottom, low to high
redshift), while different columns correspond to different HOD samples. Within
each panel, the colour gradient marks different values of $k_{\rm max}$, as
detailed in the legend. Solid/dashed lines correspond to the configurations
with all the nuisance parameters free to vary, and with $b_{{\cal G}_{2}}$
fixed to the excursion-set relation as in Eq. (35), respectively.

###  5.3 Constraints on scale-dependent shot-noise

So far, the stochastic field $\varepsilon_{\rm g}$ entering the expression for
the galaxy density field in Eq. (18) has been assumed responsibile only for a
constant offset from Poissonian predictions, via the parameter $\alpha_{P,1}$.
An immediate check on the performance of the one-loop galaxy bias expansion
can be carried out by further extending the model parameter space to also
include the next-to-leading order correction to the stochastic field
$\varepsilon_{\rm g}$. As already mentioned in Sect. 3, this leads to the
presence of an additional $k^{2}$-dependent term in the galaxy power spectrum,
whose amplitude is regulated by the extra parameter $\alpha_{P,2}$.

Figure 8 shows the marginalised one-dimensional constraints on $\alpha_{P,2}$,
for both HOD samples, Model 1 on the left and Model 3 on the right,
respectively. Since the large-scale limit of the galaxy power spectrum does
not have enough constraining power on $\alpha_{P,2}$, we only consider values
of $k_{\rm max}$ above $0.35\,h\,{\rm Mpc}^{-1}$.212121For some of the
samples, this range of scales is already dominated by the Poissonian shot-
noise contribution, as can be observed from the top and middle panel of Fig.
2. We never observe a statistically significant detection of the
$\alpha_{P,2}$ parameter, with the majority of the marginalised constraints
being consistent with $\alpha_{P,2}=0$ well within the $2\sigma$ confidence
interval. The only configurations for which this does not happen are the ones
at high redshifts, specifically when considering high values of $k_{\rm max}$,
since these are the configurations for which the parameter $\alpha_{P,2}$ is
constrained with the highest precision. Performing the same test with one of
the tidal bias fixed 222222We choose to fix $b_{{\cal G}_{2}}$, motivated by
the results of Sect. 5.2 that suggested this is the configuration less
affected by degeneracies between the model parameters does not lead to
significantly different conclusions. This seems to suggest that the
$\alpha_{P,2}$ parameter might have a more important role over a range of
redshifts where the shot-noise correction is more relevant.

We postpone further tests to a next installment of this series of papers,
since a more careful check should be carried out adopting samples with
different number densities – as this value determines the range of scales
where the transition from signal to noise takes place – in particular
considering values that would represent in a more reliable way the expected
H$\alpha$ galaxy distribution detected by Euclid. At the same time, an
important test should be carried out using the redshift-space galaxy power
spectrum (Euclid Collaboration: Camacho et al., in prep.), for which extra
$k^{2}$-dependent noise corrections are required, as a function of the
orientation with respect to the line of sight (see, e.g., Philcox & Ivanov
2022; Carrilho et al. 2023; Moretti et al. 2023). Finally, we notice that the
findings of this analysis are in line with the conclusions from Pezzotta et
al. (2021), for which, in terms of constraints on the cosmological parameters
$\left\\{h,\omega_{\rm c}\right\\}$, a clear detection of scale-dependent
stochastic parameters happens only when considering the combined information
from the galaxy-galaxy and galaxy-matter power spectra.

###  5.4 Testing the EFT model: results on cosmological parameters

Figure 9: Performance metrics – figure of bias (top), goodness of fit
(middle), and figure of merit (bottom) – of the Model 1 HOD samples, for the
various configurations defined in Sect. 5.2, as a function of the maximum wave
mode $k_{\rm max}$, and for the four different redshifts of the samples.
Different colours correspond to different model configurations, as listed in
the legend. The black dashed line shows as a reference the case in which both
tidal bias parameters, $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$, are set to 0.
The grey bands in the FoB and $\chi^{2}$ panels represent the $68\%$ and
$95\%$ percentiles of the corresponding distributions. The FoM panels show the
figure of merit normalised to the one of the standard run – with all bias
parameters free to vary – at $k_{\rm max}=0.1\,h\,{\rm Mpc}^{-1}$. Figure 10:
Same as in Fig. 9, but for the Model 3 HOD samples. Figure 11: Residuals of
the maximum-likelihood best fits against the input galaxy power spectrum data
vectors, assuming the case with all bias parameters free to vary. Different
columns correspond to different redshifts, as shown on top of the
corresponding column, while different rows mark either the Model 1 (top row)
or Model 3 (bottom row) HOD sample. Different colour shades mark the best fits
obtained at different $k_{\rm max}$ values, from $0.1\,h\,{\rm Mpc}^{-1}$ up
to $0.45\,h\,{\rm Mpc}^{-1}$.

After having investigated the performance of the EFT model, we now turn our
attention to the study of how cosmological constraints can be affected by
different choices of model configuration. In detail, we assume the same
parameter space already used in Sect. 5.1, which also includes the Hubble
parameter $h$ and the cold dark matter density parameter $\omega_{\rm c}$,
while keeping the rest of the cosmological parameters fixed to their fiducial
values, as listed in Table 1. A standard full-shape analysis of the redshift-
space galaxy power spectrum would typically also include the scalar amplitude
of the power spectrum, $A_{\rm s}$, since the anisotropies introduced by
peculiar velocities make possible to break the strong $A_{\rm s}$-$\,b_{1}$
degeneracy that is otherwise present when considering real-space
coordinates.232323Specifically, the $A_{\rm s}$-$\,b_{1}$ degeneracy can be
broken thanks to the different impact that these two parameters have on the
amplitude of the leading-order power spectrum multipoles, $P^{\,(\ell)}\propto
A_{\rm s}\,b_{1}^{2-\ell/2}$. However, since this analysis revolves around the
real-space galaxy power spectrum, we ought to choose a more conservative
approach in order to obtain as least degenerate constraints as possible on the
rest of the cosmological parameters. The sampling of $A_{\rm s}$ will be
performed in the rest of the papers of this series, when considering the
additional information content of the galaxy bispectrum (Euclid Collaboration:
Eggemeier et al., in prep.) and RSD (Euclid Collaboration: Camacho et al., in
prep., Euclid Collaboration: Pardede et al., in prep).

Figures 9 and 10 show the three performance metrics defined in Sect. 4.1 for
the Model 1 and Model 3 HOD samples respectively, assuming the full-box volume
of the Flagship I simulation. In both cases, the FoB and FoM panels refer to
the combination between the two cosmological parameters we are sampling over.
In terms of FoB, we notice that the typical trend for each model configuration
is to provide an unbiased combined measurement of the cosmological parameters,
even well within the mildly nonlinear regime, at $k_{\rm max}\gtrsim
0.3\,h\,{\rm Mpc}^{-1}$. The only exception is represented by the
configuration in which both tidal bias parameters, $b_{{\cal G}_{2}}$ and
$b_{\Gamma_{3}}$, are simultaneously kept fixed to the excursion-set relation
(Eq. 35) and to the coevolution relation (Eq. 34), respectively. This comes as
no surprise, as we are reducing by two the total number of degrees of freedom
of the model. The bias on the cosmological parameters gets larger at lower
redshift, hinting to a departure from the coevolution relations when the
effects of nonlinear gravitational evolution become more important. This
deviation typically occurs at scales of approximately $k_{\rm
max}=0.3\,h\,{\rm Mpc}^{-1}$. However, we notice how this configuration
behaves surprisingly well for high-redshift snapshots, with a FoB that is well
within the $68\%$ confidence interval.

In terms of overall stability of the results, for the Model 3 sample we notice
a deviation at low $k_{\rm max}$ values. This is affecting more heavily
samples at the highest redshifts, for which the value of the FoB at $k_{\rm
max}\sim 0.1\,h\,{\rm Mpc}^{-1}$ is already larger than the corresponding
$68\%$ confidence level, and only gets below the threshold when including the
additional signal from smaller scales. This effect is mostly related to the
presence of projection effects, given the large dimensionality of the selected
parameter space. In detail, we find that all samples display a non-negligible
correlation between the cold dark matter density parameter, $\omega_{\rm c}$,
and the EFT counterterm $c_{0}$, leading to a systematic shift of $\omega_{\rm
c}$ for the lowest values of $k_{\rm max}$, where there is not enough
constraining power to properly constrain $c_{0}$. Fixing the cubic tidal bias
$b_{\Gamma_{3}}$ to the coevolution relation typically helps in bringing the
cosmological parameters back on top of their fiducial positions. This happens
because of the further degeneracy between $b_{\Gamma_{3}}$ and $c_{0}$ over
the mildly nonlinear regime. The relative importance of the $\omega_{\rm
c}$-$c_{0}$ degeneracy gets amplified only when considering the snapshots at
the highest redshifts. As a partial confirmation of this trend, Pezzotta et
al. (2021) did not report either any low-$k$ systematic effect when analysing
mock galaxies meant to reproduce the clustering properties of the BOSS – CMASS
and LOWZ – and SDSS MGS samples, since, in that case, the considered redshift
range was much lower ($0.1\lesssim z\lesssim 0.6$) than the one analysed in
this work. As a further cross-check, the same effect is partially present when
combining the full shape of the galaxy power spectrum and bispectrum in a
joint analysis (Euclid Collaboration: Eggemeier et al., in prep.), albeit with
a lower significance, due to the additional constraining power of higher-order
statistics.

Similarly to the case at fixed cosmology, we find that the goodness of fit for
the different models is consistent among the various model configurations,
with only a small departure of the case where both tidal biases, $b_{{\cal
G}_{2}}$ and $b_{\Gamma_{3}}$, are kept fixed to their corresponding relations
and for the largest $k_{\rm max}$ values we consider in this analysis. This is
visible in the case of the Model 1 sample at $z=0.9$, for which there is an
increase in the averaged $\chi^{2}$ value at $k_{\rm max}\gtrsim 0.3\,h\,{\rm
Mpc}^{-1}$, which also corresponds to the transition of the FoB to values
above the $68\%$ percentile value. Otherwise, we find that the different model
configurations provide a systematically consistent goodness of fit, with a
reduced average $\chi^{2}$ value that is typically well within the $95\%$
percentile of the corresponding $\chi^{2}$ distribution. When considering the
two different HOD models, the $\chi^{2}$ for the Model 3 HOD samples is
consistently larger than for the Model 1 case (see the middle panel of Fig.
10). The most significant deviation is affecting the high-redshift snapshots,
for which the average $\chi^{2}$ spuriously gets larger than the $95\%$
confidence interval for some of the selected $k_{\rm max}$ values. In
practice, this deviation is still consistent to better than the $3\sigma$
confidence interval. In addition, we remind that here we are using a single
noisy realization, meaning sample variance could partially be driving some of
the constraints. Moreover, we are analysing the data vectors with an extremely
high level of precision, due to choice of using the full volume of the
simulation and to the high number density of the HOD samples. As a further
evidence for the goodness of our fits, in Fig. 11 we show the residuals
between the maximum-likelihood theory vectors obtained at different values of
$k_{\rm max}$ against the input data vectors, assuming the most relaxed model
configuration. For all the samples that are under examination, we find that
the broadband of the input galaxy power spectrum is perfectly recovered, and
that the worst performance (in terms of goodness of fit) is only imputable to
the scatter of the noisy data vector around the best fit – with a significance
that is larger for some of the samples, such as for the Model 3 sample at
$z=1.8$. Finally, in Figs. 9 and 10 we show with a dashed black line the FoB
and averaged $\chi^{2}$ of the case where the values of the tidal biases are
set to 0. In this case, we observe a departure of the goodness of fit from the
other configurations, in a redshift-dependent way, which is also accompanied
by a breaking of the model in terms of FoB. This shows how the use of
coevolution relations can drastically improve the performance of the model,
with respect to simply set the non-local bias parameters to zero.

As expected, the FoM monotonically increases when including additional
information from more nonlinear scales, with a relative gain with respect to
the most relaxed configuration – i.e., all bias parameters free to vary, at
$k_{\rm max}=0.1\,h\,{\rm Mpc}^{-1}$ – that becomes larger moving towards
lower redshifts. In fact, extending the fitting range to the maximal value of
$k_{\rm max}=0.45\,h\,{\rm Mpc}^{-1}$, we find that the trend for the FoM of
the different bias configurations is approximately $10,15,20$ and $30$ times
larger than the reference at $z=1.8,1.5,1.2$ and $0.9$, respectively. The only
exception is represented by the high-redshift snapshots of the Model 3 sample,
for which the total galaxy power spectrum on mildly nonlinear scales becomes
dominated by the shot-noise correction earlier than for the rest of the
samples, and for which we observe that the FoM reaches a plateau at $k_{\rm
max}\gtrsim 0.3\,h\,{\rm Mpc}^{-1}$. As expected, the case with two less
degrees of freedom consistently gains more constraining power on the
cosmological parameters, leading to much tighter constraints in particular at
the largest $k_{\rm max}$ value we probe. However, we notice that these gains
are directly correlated with the breaking of the model in terms of FoB at
$z=0.9$, and might therefore lead to biased cosmological constraints if used
in a real-data analysis. Nevertheless, for most of the tested cases, combining
the two relations still leads to acceptable results up to the maximal scale we
are considering.

Overall, we find that fixing only the quadratic tidal bias $b_{{\cal G}_{2}}$
leads to the most stable results, with a FoB that is typically – except for
some spurious scale cut – well within the $68\%$ percentile of the
corresponding distribution, and with a FoM which is systematically larger than
in the case where all the parameters are free to vary. The performance of the
case with a fixed cubic tidal bias is also consistent, but with the caveat
that the underlying bias parameters experience a strong degeneracy among
themselves, as shown in Sect. 5.2. However, we find that this case typically
achieves a FoM larger than the one with fixed quadratic tidal bias, with the
latter catching up only at large enough values of $k_{\rm max}$. Also, in a
range of scales up to $k_{\rm max}\sim 0.3\,h\,{\rm Mpc}^{-1}$, the case with
relations applied to both tidal biases matches almost identically the case
with fixed cubic bias, highlighting again how this parameter has a much larger
impact when constraining the cosmological parameters considered in this
analysis.

###  5.5 Dependence on sample volume

Figure 12: Evolution of the FoB and FoM as a function of the different choice
for the rescaling of the reference volume, as shown in the legend. Solid and
dashed curves correspond to the case with $b_{\Gamma_{3}}$ fixed to the
coevolution relation and to the case with all the parameters free to vary,
respectively. All configurations of the top two rows correspond to the fits of
the Model 1 HOD samples, while the bottom two rows do the same for the Model 3
HOD samples. Grey bands in the FoM panels identify the scale at which the
Poisson shot-noise contribution assumes the same height of the underlying
clustering signal, marking the transition to the shot-noise-dominated scales.

As anticipated in Sect. 2.3, in addition to performing a model selection using
extremely high-precision measurements – with a Gaussian covariance matrix
derived assuming the full box volume of the Flagship I simulation – in this
section we test the performance of the Eulerian bias expansion for smaller
values of the galaxy sample volume, in order to provide more realistic
forecasts for the analysis of the Euclid spectroscopic data. Once again, we
remind the reader that for these tests we consider HOD catalogs with a number
density larger than the one expected for the real galaxy samples, and that
more realistic mocks will be used in forthcoming analyses that will also
consider observational systematic effects, such as target purity and
incompleteness, and observational effects such as the radial and angular
selection function.

In this section, we consider the same data vectors obtained from measurements
of the Flagship simulation snapshots already used in the previous ones. The
dependence on the volume is explored by rescaling the corresponding
covariance. We consider four different volumes corresponding to possible
Euclid-like shells (i.e. redshift bins), as explained in Sect. 2.3. For each
comoving snapshot – at $z=(0.9,1.2,1.5,1.8)$ – with reference volume $V_{\rm
box}$, we define the volume $V_{\rm shell}$ of a spectroscopic bin
corresponding to a total angular surface of $15\,000$ square degrees – and
with a depth of $\Delta z=(0.2,0.2,0.2,0.3)$ – following the choices made in
Euclid Collaboration: Blanchard et al. (2020). The three additional volume
rescalings are obtained by selecting the values that divide the interval
$[\eta,1]$ into four equi-partitioned subintervals, where $\eta=V_{\rm
box}/V_{\rm shell}$.

In Fig. 12 we show the trends of the FoB and FoM for the previously defined
samples, with Model 1 and 3 in the top and bottom two rows, respectively.
Thick solid lines correspond to the case with $b_{\Gamma_{3}}$ fixed to the
coevolution relation, which we selected as one of the best performing model
among the ones that we have tested in Sect. 5.4. On the contrary, thin dashed
lines correspond to the case with all the bias parameters free to vary.
Different colors identify the different covariance matrices used in the fits,
from the already shown full-box case – in light blue/orange – to the case
corresponding to the Euclid-like shells – in dark blue/orange. As expected, we
find that in both cases the FoB for the two cosmological parameters
$\left(h,\omega_{\rm c}\right)$ becomes progressively smaller, reflecting the
increasing amplitude of the covariance matrices used in the fits of the input
data vectors. At the same time, we observe how this trend is tightly
correlated to a decrease in the FoM, with the severity of the drop being
almost proportional to the factor between the original and the rescaled
volume, $\eta$.

It comes with no surprise that the reference bias model is well-performing up
to the highest value of $k_{\rm max}$ we consider, even under realistic
assumptions. Additionally, for the case with fixed $b_{\Gamma_{3}}$, we can
observe how the precision on the cosmological constraints reaches a plateau
above a typical threshold that corresponds to the transition scale between the
regimes dominated by the signal and by the shot-noise contribution,
respectively. For this reason, it is possible to gain additional constraining
power by pushing the analysis to high values of $k_{\rm max}$ at $z=0.9$, for
which the number density of the sample is significantly larger than the one of
the high-redshift samples. In contrast, when all the parameters are free to
vary, the trend for the FoM curves is to gain additional constraints from
smaller scales, even above the scale of transition, possibly pointing to a
further breaking of parameter degeneracies that are no longer present when
fixing the value of $b_{\Gamma_{3}}$ to the coevolution relation.

As a final consideration, we notice how the Eulerian bias expansion is
performing significantly well for all the considered rescalings of the
covariance matrix. The constraining power of the EFT model in terms of the
combination $\left\\{h,\omega_{\rm c}\right\\}$ can be enhanced employing one
of the coevolution relations described in Sect. 3.1.3 without the appearance
of systematic errors, even when considering samples with a number density
significantly larger than the one expected from the real Euclid data. This
analysis, limited to real space, motivates further tests including higher-
order statistics, such as the galaxy bispectrum, and taking into account
redshift-space distortions. These topics will be properly explored in the next
entries of this series.

##  6 Conclusions

In this paper we carried out an analysis meant to assess the performance of
state-of-art models for one-loop galaxy bias over a redshift range that is
well-representative of the spectroscopic galaxy sample that will be one of the
main target of Euclid. We employed a set of four FoF halo catalogs from
comoving snapshots of the Flagship I simulation, at $z=(0.9,1.2,1.5,1.8)$,
which have been subsequently populated with H$\alpha$ galaxies using HOD
prescriptions based on the Model 1 and 3 from Pozzetti et al. (2016). Each
snapshot features an outstanding volume of $(3780\,h^{-1}\,{\rm Mpc})^{3}$ and
a high comoving number density (from $\sim 10^{-4}\,h^{3}\,{\rm Mpc}^{-3}$ to
$\sim 10^{-3}\,h^{3}\,{\rm Mpc}^{-3}$), corresponding to a flux limit of
$f_{{\rm H}\alpha}=2\times 10^{-16}\,{\rm erg\,cm^{-2}\,s^{-1}}$. They can
therefore be used to assess the accuracy of the perturbative bias expansion at
a high level of precision.

We tested two galaxy bias models for the full shape of the real-space galaxy
power spectrum. The first one adopts an Eulerian bias expansion, and it is
based on the recently developed EFTofLSS modelling, in which the impact of
small-scale physics, as well as the integration of ultraviolet modes in SPT,
can be captured by a set of counterterms, which reduce to a single one when
considering real-space coordinates. The final parameter space consists of two
cosmological parameters, the Hubble parameter $h$ and the cold dark matter
density parameter $\omega_{\rm c}$, plus a set of six nuisance parameters,
consisting of the linear bias $b_{1}$, the quadratic bias $b_{2}$, the tidal
quadratic and cubic biases, $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$, the
matter counterterm and higher-derivative bias $c_{0}$, and two extra
parameters representing deviations from Poissonian shot-noise: a constant
offset, $\alpha_{P,1}$, and a scale-dependent term, $\alpha_{P,2}$. The second
model adopts a similar one-loop expansion of the galaxy power spectrum,
although using Lagrangian coordinates, and is based on the emulation of the
individual terms of such expansion starting from a limited set of high-
resolution N-body simulations. This is achieved thanks to the cosmology-
rescaling technique presented in Angulo et al. (2021). This model features, in
addition to the two cosmological parameters, the linear bias $b_{1}$, the
quadratic bias $b_{2}$, the tidal quadratic bias $b_{s^{2}}$, the Laplacian
bias $b_{\nabla^{2}\delta}$, and a stochastic parameter representing a
constant deviation from Poissonian shot-noise, $\alpha_{P,1}$. In this work we
made use of the two implementations available in the public codes COMET
(Eggemeier et al. 2022) for the EFT model – with accuracy tests carried out
against external benchmarks, as shown in Appendix D – and baccoemu (Zennaro et
al. 2023) for the hybrid model.

In the main section of this work we tested the relative performance of these
two galaxy bias models, while in the next sections we determined the range of
validity of the EFT model and tested the impact of fixing one or more
parameters of the Eulerian bias expansion to some physically-motivated
relations, in a way that allows us to break strong parameter degeneracies and
better constrain the cosmological parameters. In all the cases we determined
the range of validity of a given bias relation and scale cut by means of three
different performance metrics, the goodness of fit, the figure of bias, and
the figure of merit. The last two metrics are computed on the
$\left\\{h,\omega_{\rm c}\right\\}$ combination, in order to quantify how
accurate and precise the model is in terms of these two parameters.

We compare the performance of the Eulerian and hybrid Lagrangian bias models
using a rescaled covariance to match the size of Euclid-like redshift shells
(assuming a full-sky area of $15\,000\deg^{2}$). Our results highlight how
both models are capable of providing unbiased measurements of the cosmological
parameters up to $k_{\rm max}=0.45\,h\,{\rm Mpc}^{-1}$ for all four redshift
snapshots, consistently with the $1\sigma$ confidence interval for the FoB
distribution. In terms of FoM, baccoemu reaches the same amplitude of the
maximally-achievable FoM of the EFT case already at a lower $k_{\rm max}$
value, most likely due to the absence of the cubic tidal bias parameter that
is instead free to vary in the EFT chains. As expected, when fixing this
parameter the EFT model performs similarly to baccoemu, and for most of the
configurations it achieves a slightly larger FoM. When considering the
covariance matrix corresponding to the full-box of the comoving snapshots
(from three to seven times larger than the Euclid-like shells, depending on
the considered redshift) we find that the EFT model manages to recover the
cosmological parameters consistently with the $1\sigma$ confidence interval.
The exception are some spurious cases at low values of $k_{\rm max}$ that are
affected by projection effects, which can be alleviated by fixing the value of
$b_{\Gamma_{3}}$. On the other hand, with this level of precision we hit the
intrinsic emulation error of baccoemu, which leads to bias in the inferred
parameters when considering high $k_{\rm max}$ values at low redshift.
Including an extra component to the covariance matrix – corresponding to
$0.5\%$ of the amplitude of the galaxy power spectrum, based on the combined
systematic error from the emulation and from the measurements used to train
the emulator – brings the FoB of baccoemu back within the $1\sigma$ confidence
interval. Also in this case, on scales for which the intrinsic error of
baccoemu is not the dominant contribution, $k_{\rm max}\lesssim 0.2\,h\,{\rm
Mpc}^{-1}$, its FoM is consistent with the one of the EFT model. Not including
the $0.5\%$ extra error, we notice how the FoM is consistent with – in some
cases even larger than – the corresponding EFT results, showing the potential
gain that one can obtain by improving the accuracy of next-generation
emulators of the full-shape of the galaxy power spectrum.

We then focus exclusively on the EFT model: we first perform fits at fixed
cosmology to check the self-consistency of the one-loop galaxy bias expansion
in terms of the linear bias $b_{1}$ and the scale-independent shot-noise
parameter $\alpha_{P,1}$. The fiducial values for these parameters have been
fitted from the measured galaxy to matter power spectrum ratio, assuming a
leading-order recipe on scales $k<0.08\,h\,{\rm Mpc}^{-1}$. The result is
that, when leaving all parameters free to vary, it is possible to recover at
better than $2\sigma$ the value of both parameters for the majority of the
samples – two HOD models times four different redshifts – and scale cuts, up
to $k_{\rm max}=0.45\,h\,{\rm Mpc}^{-1}$. The only significant deviation takes
place at the lowest redshift we consider, $z=0.9$, for which we observe a
departure from the fiducial values soon after $k_{\rm max}=0.2\,h\,{\rm
Mpc}^{-1}$. The latter is however consistent with sample variance
expectations, as observed in a set of ten different noisy realizations of a
synthetic theory data vector (see Appendix F). The systematic errors are
partially alleviated when considering the position in parameter space
corresponding to the maximum of the likelihood, showing how the deviations
might be imputable to projection effects. Fixing one of the two tidal biases
to either an excursion-set-derived relation or to the coevolution relation
still results in constraints that are consistent with the $2\sigma$ confidence
interval. In particular, the former is preferred in terms of constraints of
the model parameters because of the simultaneous breaking of the strong
degeneracies with both the quadratic bias $b_{2}$ and the cubic tidal bias
$b_{\Gamma_{3}}$. We find that fixing both parameters at the same time still
works extremely well (with a typical recovery of $b_{1}$ and $\alpha_{P,1}$
within the 68% confidence interval), even for the largest $k_{\rm max}$
values, when considering the high-redshift snapshots. However, at low redshift
this choice can lead to deviations of more than $3\sigma$ for some of the
configurations we test, especially when the maximum scale included in the fit
is above a typical scale of $k_{\rm max}\sim 0.3\,h\,{\rm Mpc}^{-1}$.

We explicitly check if a next-to-leading order correction to the shot-noise
contribution ($\alpha_{P,2}$) can improve the model performance. In all cases
considered, the marginalised posterior distribution for $\alpha_{P,2}$ is
consistent with zero within $2\sigma$, suggesting this additional parameter is
not needed, at least for the description of these galaxy samples.
Additionally, the scales that can constrain this parameter soon become
dominated by the underlying shot-noise correction, effectively breaking the
perturbative description of the latter in a Taylor expansion. A more
significant test should be carried out considering more realistic galaxy
samples, in terms of galaxy number density, and also including RSD and
observational systematic effects.

When the parameter space is extended to also include the two cosmological
parameters, we notice a good recovery of the fiducial values across the whole
range of separations we test. The FoB exhibits an increasing trend moving
towards high redshifts and low value of $k_{\rm max}$, due to projection
effects when marginalising over all the nuisance parameters. This trend can be
partially corrected by fixing the cubic tidal bias $b_{\Gamma_{3}}$ to the
coevolution relation, and indeed with this configuration it is possible to
consistently recover unbiased – within the $68\%$ confidence interval –
constraints on the $\left(h,\omega_{\rm c}\right)$ pair. Also in this case,
fixing both tidal biases at the same time can lead to biased cosmological
constraints, with the amplitude of the systematic errors increasing towards
lower redshifts, pointing to a premature breaking of the tested relations. In
terms of goodness of fit, we do not observe a significant change in the
average $\chi^{2}$ when fixing some of the model parameters to the relations
presented in Sect. 3.1.3. Finally, the FoM of the cosmological parameters
clearly increases when reducing the degrees of freedom of the model. However,
since the configuration with both tidal biases fixed results in biased
constraints, we find that the case with only $b_{\Gamma_{3}}$ fixed is the one
mostly preferred by our data vectors, with a relative gain in FoM that ranges
from 1.5 to 2 times, with respect to the case with all the parameters free.
Relative gains in the FoM are more concentrated at $k_{\rm max}\lesssim
0.3\,h\,{\rm Mpc}^{-1}$, where the model with fixed $b_{\Gamma_{3}}$
experiences a steep increase that is after followed by a more modest growth.
The configuration with $b_{{\cal G}_{2}}$ fixed also displays a FoM larger
than the case with all parameters free to vary, but with a steady slope that
manages to catch up with the other configuration only for the largest $k_{\rm
max}$ values that we tested. Overall, we find that the one-loop galaxy bias
expansion is sufficiently accurate on the redshift range that we are
exploring, $1\lesssim z\lesssim 2$, even deep within the mildly nonlinear
regime, at $k_{\rm max}\sim 0.4\,h\,{\rm Mpc}^{-1}$, with a statistical
significance on the cosmological parameters that can be enhanced by fixing
some of the degrees of freedom of the model.

In order to understand the impact of a different statistical precision on the
input data vectors, we rescaled the Gaussian covariance matrix used in the
fitting procedure to match the volume of a Euclid-like spectroscopic bin, with
three additional intermediate volume choices selected between the Euclid-like
bin and the original volume of the comoving box. A smaller volume thus
corresponds to a reduced amplitude in the covariance matrix, resulting in a
decrease in both FoB and FoM in a way that is proportional to the fraction of
lost volume. We therefore confirm that these models of galaxy bias can be
eventually used to analyse the real spectroscopic data collected by Euclid.

This paper stands as a first installment of a series of works meant to
validate the theoretical framework that will be used to analyse the large-
scale galaxy distribution as observed in the actual measurements of Euclid.
Here we focus on the modelling of the real-space galaxy power spectrum of the
spectroscopic sample, something that stands as an important test for the
complementary photometric analysis that is going to be carried out by Euclid,
in the shape of the popular 3$\times$2-point data combination. On a parallel
track, in Euclid Collaboration: Eggemeier et al. (in prep.) we consider the
joint analysis of the real-space galaxy power spectrum and bispectrum,
exploring a consistent description of non-linear bias in both observables. Two
additional installments of the series (Euclid Collabpration: Camacho et al.,
in prep., Euclid Collaboration: Pardede et al., in prep.) will extend the
modelling tests to redshift space. In parallel, a different set of papers will
be devoted to a similar analysis of configuration-space statistics (Euclid
Collaboration: Guidi et al., in prep., Euclid Collaboration: Kärcher et al.,
in prep., Euclid Collaboration: Pugno et al., in prep.).

###### Acknowledgements.

The Euclid Consortium acknowledges the European Space Agency and a number of
agencies and institutes that have supported the development of Euclid, in
particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian
Science Policy, the Canadian Euclid Consortium, the French Centre National
d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish
Space Research Institute, the Fundação para a Ciência e a Tecnologia, the
Ministerio de Ciencia e Innovación, the National Aeronautics and Space
Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org). The majority of the analysis carried out in this
manuscript has been produced by a joint effort among several Euclid members
and work centers. A special acknowledgement goes to the Max Planck Computing
and Data Facility (MPCDF) in Garching (Germany), where most of the results
presented in this article have been obtained. We acknowledge the hospitality
of the Institute for Fundamental Physics of the Universe (IFPU) of Trieste for
the group meeting held there in November 2022. C.M.’s research for this
project was supported by a UK Research and Innovation Future Leaders
Fellowship [grant MR/S016066/2]. C.M.’s work is supported by the Fondazione
ICSC, Spoke 3 Astrophysics and Cosmos Observations, National Recovery and
Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID
CN_00000013 “Italian Research Center on High-Performance Computing, Big Data
and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4:
Potenziamento strutture di ricerca e creazione di ”campioni nazionali di R&S
(M4C2-19 )” - Next Generation EU (NGEU). A.E. is supported at the Argelander
Institut für Astronomie by an Argelander Fellowship. A.B. 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. R.E.A. acknowledges the support of
the ERC-StG number 716151 (BACCO), and the project PID2021-128338NB-I00 from
the Spanish Ministry of Science. M.B. is supported by the Programma Nazionale
della Ricerca (PNR) grant J95F21002830001 with title “FAIR-by-design”. M.K. is
funded by the Excellence Initiative of Aix-Marseille University - A*MIDEX, a
French “Investissements d’Aveni” programme (AMX-19-IET-008 - IPhU). This
research made use of matplotlib, a Python library for publication quality
graphics (Hunter 2007).

## References

  * Abidi & Baldauf (2018) Abidi, M. M. & Baldauf, T. 2018, J. Cosmology Astropart. Phys., 07(2018), 029
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617
  * Angulo et al. (2015) Angulo, R., Fasiello, M., Senatore, L., & Vlah, Z. 2015, J. Cosmology Astropart. Phys., 09(2015), 029
  * Angulo & White (2010) Angulo, R. E. & White, S. D. M. 2010, MNRAS, 405, 143
  * Angulo et al. (2021) Angulo, R. E., Zennaro, M., Contreras, S., et al. 2021, MNRAS, 507, 5869
  * Assassi et al. (2014) Assassi, V., Baumann, D., Green, D., & Zaldarriaga, M. 2014, J. Cosmology Astropart. Phys., 08(2014), 056
  * Baldauf et al. (2016) Baldauf, T., Codis, S., Desjacques, V., & Pichon, C. 2016, MNRAS, 456, 3985
  * Baldauf et al. (2015a) Baldauf, T., Mercolli, L., & Zaldarriaga, M. 2015a, Phys. Rev. D, 92, 123007
  * Baldauf et al. (2015b) Baldauf, T., Mirbabayi, M., Simonović, M., & Zaldarriaga, M. 2015b, Phys. Rev. D, 92, 043514
  * Baldauf et al. (2013) Baldauf, T., Seljak, U., Smith, R. E., Hamaus, N., & Desjacques, V. 2013, Phys. Rev. D, 88, 083507
  * Baldauf et al. (2012) Baldauf, T., Seljak, U. c. v., Desjacques, V., & McDonald, P. 2012, Phys. Rev. D, 86, 083540
  * Bardeen et al. (1986) Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15
  * Barreira et al. (2021) Barreira, A., Lazeyras, T., & Schmidt, F. 2021, J. Cosmology Astropart. Phys., 08(2021), 029
  * Baumann et al. (2012) Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, J. Cosmology Astropart. Phys., 07(2021), 051
  * Bernardeau et al. (2002) Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep, 367, 1
  * Bertolini et al. (2016) Bertolini, D., Schutz, K., Solon, M. P., Walsh, J. R., & Zurek, K. M. 2016, Phys. Rev. D, 93, 123505
  * Beutler et al. (2012) Beutler, F., Blake, C., Colless, M., et al. 2012, MNRAS, 423, 3430
  * Beutler et al. (2016) Beutler, F., Seo, H.-J., Ross, A. J., et al. 2016, MNRAS, 464, 3409
  * Blake et al. (2011) Blake, C., Brough, S., Colless, M., et al. 2011, MNRAS, 415, 2876
  * Blas et al. (2016) Blas, D., Garny, M., Ivanov, M. M., & Sibiryakov, S. 2016, J. Cosmology Astropart. Phys., 07(2016), 028
  * Blot et al. (2015) Blot, L., Corasaniti, P. S., Alimi, J.-M., Reverdy, V., & Rasera, Y. 2015, MNRAS, 446, 1756
  * Blot et al. (2016) Blot, L., Corasaniti, P. S., Amendola, L., & Kitching, T. D. 2016, MNRAS, 458, 4462
  * Blot et al. (2019) Blot, L., Crocce, M., Sefusatti, E., et al. 2019, MNRAS, 485, 2806
  * Bonoli & Pen (2009) Bonoli, S. & Pen, U. L. 2009, MNRAS, 396, 1610
  * Bose et al. (2020) Bose, B., Pourtsidou, A., Markovič, K., & Beutler, F. 2020, MNRAS, 493, 5301
  * Bouchet et al. (1992) Bouchet, F. R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Carrasco et al. (2014) Carrasco, J. J. M., Foreman, S., Green, D., & Senatore, L. 2014, J. Cosmology Astropart. Phys., 07(2014), 057
  * Carrasco et al. (2012) Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, J. High Energy. Phys., 09(2012), 82
  * Carrilho et al. (2023) Carrilho, P., Moretti, C., & Pourtsidou, A. 2023, J. Cosmology Astropart. Phys., 01(2023), 028
  * Catelan et al. (1998) Catelan, P., Lucchin, F., Matarrese, S., & Porciani, C. 1998, MNRAS, 297, 692
  * Catelan et al. (2000) Catelan, P., Porciani, C., & Kamionkowski, M. 2000, MNRAS, 318, L39
  * Chan et al. (2012) Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509
  * Chartier et al. (2021) Chartier, N., Wandelt, B., Akrami, Y., & Villaescusa-Navarro, F. 2021, MNRAS, 503, 1897
  * Chen et al. (2022) Chen, S.-F., Vlah, Z., & White, M. 2022, J. Cosmology Astropart. Phys., 02(2022), 008
  * Chudaykin et al. (2020) Chudaykin, A., Ivanov, M. M., Philcox, O. H., & Simonović, M. 2020, Phys. Rev. D, 102, 063533
  * Cole & Kaiser (1989) Cole, S. & Kaiser, N. 1989, MNRAS, 237, 1127
  * Coles (1993) Coles, P. 1993, MNRAS, 262, 1065
  * Colless et al. (2001) Colless, M., Dalton, G., Maddox, S., et al. 2001, MNRAS, 328, 1039
  * Contreras et al. (2013) Contreras, C., Blake, C., Poole, G. B., et al. 2013, MNRAS, 430, 924
  * Contreras et al. (2020) Contreras, S., Angulo, R. E., Zennaro, M., Aricò, G., & Pellejero-Ibañez, M. 2020, MNRAS, 499, 4905
  * Cooray & Sheth (2002) Cooray, A. & Sheth, R. 2002, Phys. Rep, 372, 1
  * Crocce & Scoccimarro (2006) Crocce, M. & Scoccimarro, R. 2006, Phys. Rev. D, 73, 063520
  * Crocce & Scoccimarro (2008) Crocce, M. & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533
  * d’Amico et al. (2020) d’Amico, G., Gleyzes, J., Kokron, N., et al. 2020, J. Cosmology Astropart. Phys., 05(2020), 005
  * Dawson et al. (2016) Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44
  * Dawson et al. (2012) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2012, AJ, 145, 10
  * de la Bella et al. (2020) de la Bella, L. F., Regan, D., Seery, D., & Parkinson, D. 2020, J. Cosmology Astropart. Phys., 07(2020), 011
  * Dekel & Lahav (1999) Dekel, A. & Lahav, O. 1999, ApJ, 520, 24
  * DeRose et al. (2019) DeRose, J., Wechsler, R. H., Tinker, J. L., et al. 2019, ApJ, 875, 69
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv:1611.00036
  * Desjacques (2008) Desjacques, V. 2008, Phys. Rev. D, 78, 103503
  * Desjacques et al. (2010) Desjacques, V., Crocce, M., Scoccimarro, R., & Sheth, R. K. 2010, Phys. Rev. D, 82, 103529
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep, 733, 1
  * Donath & Senatore (2020) Donath, Y. & Senatore, L. 2020, J. Cosmology Astropart. Phys., 10(2020), 039
  * Drinkwater et al. (2010) Drinkwater, M. J., Jurek, R. J., Blake, C., et al. 2010, MNRAS, 401, 1429
  * Driver et al. (2011) Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971
  * Eggemeier et al. (2022) Eggemeier, A., Camacho-Quevedo, B., Pezzotta, A., et al. 2022, MNRAS, 519, 2962
  * Eggemeier et al. (2020) Eggemeier, A., Scoccimarro, R., Crocce, M., Pezzotta, A., & Sánchez, A. G. 2020, Phys. Rev. D, 102, 103530
  * Eggemeier et al. (2019) Eggemeier, A., Scoccimarro, R., & Smith, R. E. 2019, Phys. Rev. D, 99, 123514
  * Eggemeier et al. (2021) Eggemeier, A., Scoccimarro, R., Smith, R. E., et al. 2021, Phys. Rev. D, 103, 123550
  * Eisenstein & Hu (1998) Eisenstein, D. J. & Hu, W. 1998, ApJ, 496, 605
  * Eisenstein et al. (2007) Eisenstein, D. J., Seo, H.-J., & White, M. 2007, ApJ, 664, 660
  * Eisenstein et al. (2011) Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72
  * Euclid Collaboration: Blanchard et al. (2020) Euclid Collaboration: Blanchard, A., Camera, S., Carbone, C., et al. 2020, A&A, 642, A191
  * Euclid Collaboration: Knabenhans et al. (2019) Euclid Collaboration: Knabenhans, M., Stadel, J., Marelli, S., et al. 2019, MNRAS, 484, 5509
  * Euclid Collaboration: Knabenhans et al. (2021) Euclid Collaboration: Knabenhans, M., Stadel, J., Potter, D., et al. 2021, MNRAS, 505, 2840
  * Feldman et al. (2001) Feldman, H. A., Frieman, J. A., Fry, J. N., & Scoccimarro, R. 2001, Phys. Rev. Lett., 86, 1434
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10
  * Fisher (1995) Fisher, K. B. 1995, ApJ, 448, 494
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Frusciante & Sheth (2012) Frusciante, N. & Sheth, R. K. 2012, J. Cosmology Astropart. Phys., 11(2012), 016
  * Fry (1996) Fry, J. N. 1996, ApJ, 461, L65
  * Fry & Gaztanaga (1993) Fry, J. N. & Gaztanaga, E. 1993, ApJ, 413, 447
  * Gelman et al. (2014) Gelman, A., Carlin, J. B., Stern, H. S., et al. 2014, Bayesian Data Analysis (Chapman & Hall/CRC)
  * Giblin et al. (2019) Giblin, B., Cataneo, M., Moews, B., & Heymans, C. 2019, MNRAS, 490, 4826
  * Gil-Marín et al. (2018) Gil-Marín, H., Guy, J., Zarrouk, P., et al. 2018, MNRAS, 477, 1604
  * Gil-Marín et al. (2015) Gil-Marín, H., Noreña, J., Verde, L., et al. 2015, MNRAS, 451, 539
  * Ginzburg et al. (2017) Ginzburg, D., Desjacques, V., & Chan, K. C. 2017, Phys. Rev. D, 96, 083528
  * Goodman & Weare (2010) Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Goroff et al. (1986) Goroff, M. H., Grinstein, B., Rey, S. J., & Wise, M. B. 1986, ApJ, 311, 6
  * Grieb et al. (2016) Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., & Dalla Vecchia, C. 2016, MNRAS, 457, 1577
  * Grieb et al. (2017) Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., et al. 2017, MNRAS, 467, 2085
  * Gualdi et al. (2021) Gualdi, D., Gil-Marín, H., & Verde, L. 2021, J. Cosmology Astropart. Phys., 07(2021), 008
  * Guzzo et al. (2008) Guzzo, L., Pierleoni, M., Meneux, B., et al. 2008, Nature, 451, 541
  * Guzzo et al. (2014) Guzzo, L., Scodeggio, M., Garilli, B., et al. 2014, A&A, 566, A108
  * Hahn (2005) Hahn, T. 2005, Comp. Phys. Comm., 168, 78
  * Hamann et al. (2010) Hamann, J., Hannestad, S., Lesgourgues, J., Rampf, C., & Wong, Y. Y. 2010, J. Cosmology Astropart. Phys., 07(2010), 022
  * Hamaus et al. (2010) Hamaus, N., Seljak, U., Desjacques, V., Smith, R. E., & Baldauf, T. 2010, Phys. Rev. D, 82, 043515
  * Hamilton (1992) Hamilton, A. J. S. 1992, ApJ, 385, L5
  * Hamilton (2000) Hamilton, A. J. S. 2000, MNRAS, 312, 257
  * Hastings (1970) Hastings, W. K. 1970, Biometrika, 57, 97
  * Heitmann et al. (2013) Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2013, ApJ, 780, 111
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19
  * Hou et al. (2018) Hou, J., Sánchez, A. G., Scoccimarro, R., et al. 2018, MNRAS, 480, 2521
  * Howlett et al. (2015) Howlett, C., Ross, A. J., Samushia, L., Percival, W. J., & Manera, M. 2015, MNRAS, 449, 848
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Ivanov et al. (2020) Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, J. Cosmology Astropart. Phys., 05(2020), 042
  * Izard et al. (2016) Izard, A., Crocce, M., & Fosalba, P. 2016, MNRAS, 459, 2327
  * Jain & Bertschinger (1994) Jain, B. & Bertschinger, E. 1994, ApJ, 431, 495
  * Jones et al. (2009) Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683
  * Kaiser (1984) Kaiser, N. 1984, ApJ, 284, L9
  * Kaiser (1987) Kaiser, N. 1987, MNRAS, 227, 1
  * Kokron et al. (2022) Kokron, N., Chen, S.-F., White, M., DeRose, J., & Maus, M. 2022, J. Cosmology Astropart. Phys., 09(2022), 059
  * Kokron et al. (2021) Kokron, N., DeRose, J., Chen, S.-F., White, M., & Wechsler, R. H. 2021, MNRAS, 505, 1422
  * Kuhlen et al. (2012) Kuhlen, M., Vogelsberger, M., & Angulo, R. 2012, Phys. Dark Univ., 1, 50
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193
  * Lazeyras et al. (2021) Lazeyras, T., Barreira, A., & Schmidt, F. 2021, J. Cosmology Astropart. Phys., 10(2021), 063
  * Lazeyras & Schmidt (2018) Lazeyras, T. & Schmidt, F. 2018, J. Cosmology Astropart. Phys., 09(2018), 008
  * Le Fèvre et al. (2013) Le Fèvre, Cassata, P., Cucciati, O., et al. 2013, A&A, 559, A14
  * Liu et al. (2018) Liu, J., Bird, S., Zorrilla Matilla, J. M., et al. 2018, J. Cosmology Astropart. Phys., 03(2018), 049
  * Manera & Gaztañaga (2011) Manera, M. & Gaztañaga, E. 2011, MNRAS, 415, 383
  * Manera et al. (2010) Manera, M., Sheth, R. K., & Scoccimarro, R. 2010, MNRAS, 402, 589
  * Markovič et al. (2019) Markovič, K., Bose, B., & Pourtsidou, A. 2019, Open J. Astrophys., 2, 13
  * Matsubara (1999) Matsubara, T. 1999, ApJ, 525, 543
  * Matsubara (2008) Matsubara, T. 2008, Phys. Rev. D, 77, 063530
  * McDonald (2006) McDonald, P. 2006, Phys. Rev. D, 74, 103512
  * McDonald & Roy (2009) McDonald, P. & Roy, A. 2009, J. Cosmology Astropart. Phys., 08(2009), 020
  * McEwen et al. (2016) McEwen, J. E., Fang, X., Hirata, C. M., & Blazek, J. A. 2016, J. Cosmology Astropart. Phys., 09(2016), 015
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087
  * Mo & White (1996) Mo, H. J. & White, S. D. M. 1996, MNRAS, 282, 347
  * Modi et al. (2020) Modi, C., Chen, S.-F., & White, M. 2020, MNRAS, 492, 5754
  * Monaco et al. (2002) Monaco, P., Theuns, T., & Taffoni, G. 2002, MNRAS, 331, 587
  * Moradinezhad Dizgah et al. (2021) Moradinezhad Dizgah, A., Biagetti, M., Sefusatti, E., Desjacques, V., & Noreña, J. 2021, J. Cosmology Astropart. Phys., 05(2021), 015
  * Moretti et al. (2023) Moretti, C., Tsedrik, M., Carrilho, P., & Pourtsidou, A. 2023, arXiv:2306.09275
  * Nicola et al. (2023) Nicola, A., Hadzhiyska, B., Findlay, N., et al. 2023, arXiv:2307.03226
  * Nishimichi et al. (2019) Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29
  * Nusser & Davis (1994) Nusser, A. & Davis, M. 1994, ApJ, 421, L1
  * Oddo et al. (2021) Oddo, A., Rizzo, F., Sefusatti, E., Porciani, C., & Monaco, P. 2021, J. Cosmology Astropart. Phys., 11(2021), 038
  * Oddo et al. (2020) Oddo, A., Sefusatti, E., Porciani, C., Monaco, P., & Sánchez, A. G. 2020, J. Cosmology Astropart. Phys., 03(2020), 056
  * Okumura et al. (2016) Okumura, T., Hikage, C., Totani, T., et al. 2016, PASJ, 68, 38
  * Osato et al. (2019) Osato, K., Nishimichi, T., Bernardeau, F., & Taruya, A. 2019, Phys. Rev. D, 99, 063530
  * Pandey et al. (2022) Pandey, S., Krause, E., DeRose, J., et al. 2022, Phys. Rev. D, 106, 043520
  * Peacock et al. (2001) Peacock, J. A., Cole, S., Norberg, P., et al. 2001, Nature, 410, 169
  * Pellejero-Ibañez et al. (2023) Pellejero-Ibañez, M., Angulo, R. E., Zennaro, M., et al. 2023, MNRAS, 520, 3725
  * Perko et al. (2016) Perko, A., Senatore, L., Jennings, E., & Wechsler, R. H. 2016, arXiv:1610.09321
  * Pezzotta et al. (2021) Pezzotta, A., Crocce, M., Eggemeier, A., Sánchez, A. G., & Scoccimarro, R. 2021, Phys. Rev. D, 104, 043531
  * Pezzotta et al. (2017) Pezzotta, A., de la Torre, S., Bel, J., et al. 2017, A&A, 604, A33
  * Philcox & Ivanov (2022) Philcox, O. H. & Ivanov, M. M. 2022, Phys. Rev. D, 105, 043517
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Pollack et al. (2012) Pollack, J. E., Smith, R. E., & Porciani, C. 2012, MNRAS, 420, 3469
  * Pollack et al. (2014) Pollack, J. E., Smith, R. E., & Porciani, C. 2014, MNRAS, 440, 555
  * Porredon et al. (2022) Porredon, A., Crocce, M., Elvin-Poole, J., et al. 2022, Phys. Rev. D, 106, 103530
  * Potter et al. (2017) Potter, D., Stadel, J., & Teyssier, R. 2017, CompAC, 4, 2
  * Pozzetti et al. (2016) Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3
  * Pueblas & Scoccimarro (2009) Pueblas, S. & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504
  * Reid et al. (2016) Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553
  * Reid et al. (2012) Reid, B. A., Samushia, L., White, M., et al. 2012, MNRAS, 426, 2719
  * Rizzo et al. (2023) Rizzo, F., Moretti, C., Pardede, K., et al. 2023, J. Cosmology Astropart. Phys., 01(2023), 031
  * Roth & Porciani (2011) Roth, N. & Porciani, C. 2011, MNRAS, 415, 829
  * Saito et al. (2014) Saito, S., Baldauf, T., Vlah, Z., et al. 2014, Phys. Rev. D, 90, 123522
  * Sánchez et al. (2014) Sánchez, A. G., Montesano, F., Kazin, E. A., et al. 2014, MNRAS, 440, 2692
  * Sánchez et al. (2016) Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2016, MNRAS, 464, 1640
  * Scaramella et al. (2022) Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297
  * Scherrer & Weinberg (1998) Scherrer, R. J. & Weinberg, D. H. 1998, ApJ, 504, 607
  * Schmidt (2016) Schmidt, F. 2016, Phys. Rev. D, 93, 063512
  * Schmidt et al. (2013) Schmidt, F., Jeong, D., & Desjacques, V. 2013, Phys. Rev. D, 88, 023515
  * Schmittfull et al. (2019) Schmittfull, M., Simonović, M., Assassi, V., & Zaldarriaga, M. 2019, Phys. Rev. D, 100, 043514
  * Schneider et al. (2016) Schneider, A., Teyssier, R., Potter, D., et al. 2016, J. Cosmology Astropart. Phys., 04(2016), 047
  * Scoccimarro (2004) Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007
  * Scoccimarro et al. (1999) Scoccimarro, R., Couchman, H. M. P., & Frieman, J. A. 1999, ApJ, 517, 531
  * Scoccimarro et al. (2001) Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20
  * Scoccimarro et al. (1999) Scoccimarro, R., Zaldarriaga, M., & Hui, L. 1999, ApJ, 527, 1
  * Sefusatti et al. (2006) Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, Phys. Rev. D, 74, 023522
  * Sefusatti et al. (2016) Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624
  * Semenaite et al. (2023) Semenaite, A., Sánchez, A. G., Pezzotta, A., et al. 2023, MNRAS, 521, 5013
  * Semenaite et al. (2022) Semenaite, A., Sánchez, A. G., Pezzotta, A., et al. 2022, MNRAS, 512, 5657
  * Senatore (2015) Senatore, L. 2015, J. Cosmology Astropart. Phys., 11(2015), 007
  * Senatore & Zaldarriaga (2014) Senatore, L. & Zaldarriaga, M. 2014, arXiv:1409.1225
  * Seo et al. (2008) Seo, H.-J., Siegel, E. R., Eisenstein, D. J., & White, M. 2008, ApJ, 686, 13
  * Sheth et al. (2013) Sheth, R. K., Chan, K. C., & Scoccimarro, R. 2013, Phys. Rev. D, 87, 083002
  * Sheth & Lemson (1999) Sheth, R. K. & Lemson, G. 1999, MNRAS, 304, 767
  * Sheth & Tormen (1999) Sheth, R. K. & Tormen, G. 1999, MNRAS, 308, 119
  * Simonović et al. (2018) Simonović , M., Baldauf, T., Zaldarriaga, M., Carrasco, J. J., & Kollmeier, J. A. 2018, J. Cosmology Astropart. Phys., 04(2018), 030
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833
  * Smith et al. (2007) Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys. Rev. D, 75, 063512
  * Springel et al. (2021) Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871
  * Strauss et al. (2002) Strauss, M. A., Weinberg, D. H., Lupton, R. H., et al. 2002, AJ, 124, 1810
  * Szalay (1988) Szalay, A. S. 1988, ApJ, 333, 21
  * Taruya et al. (2010) Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522
  * Taruya & Soda (1999) Taruya, A. & Soda, J. 1999, ApJ, 522, 46
  * Tassev et al. (2013) Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, J. Cosmology Astropart. Phys., 06(2013), 036
  * Tegmark et al. (2006) Tegmark, M., Eisenstein, D. J., Strauss, M. A., et al. 2006, Phys. Rev. D, 74, 123507
  * Tröster et al. (2020) Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10
  * Vlah et al. (2016) Vlah, Z., Seljak, U., Yat Chu, M., & Feng, Y. 2016, J. Cosmology Astropart. Phys., 03(2016), 057
  * Wadekar et al. (2020) Wadekar, D., Ivanov, M. M., & Scoccimarro, R. 2020, Phys. Rev. D, 102, 123521
  * Wadekar & Scoccimarro (2020) Wadekar, D. & Scoccimarro, R. 2020, Phys. Rev. D, 102, 123517
  * Wang (2008) Wang, Y. 2008, Phys. Rev. D, 77, 123525
  * Wang et al. (2018) Wang, Y., Zhao, G.-B., Chuang, C.-H., et al. 2018, MNRAS, 481, 3160
  * Wibking et al. (2018) Wibking, B. D., Salcedo, A. N., Weinberg, D. H., et al. 2018, MNRAS, 484, 989
  * Winther et al. (2019) Winther, H. A., Casas, S., Baldi, M., et al. 2019, Phys. Rev. D, 100, 123540
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84
  * Zennaro et al. (2019) Zennaro, M., Angulo, R. E., Aricò, G., Contreras, S., & Pellejero-Ibáñez, M. 2019, MNRAS, 489, 5938
  * Zennaro et al. (2023) Zennaro, M., Angulo, R. E., Pellejero-Ibáñez, M., et al. 2023, MNRAS, 524, 2407
  * Zhao et al. (2018) Zhao, G.-B., Wang, Y., Saito, S., et al. 2018, MNRAS, 482, 3497

##  Appendix A Matter power spectrum fits and the fiducial cosmology

In the main body of this article, we carried out tests meant to assess the
level of accuracy of different models for the one-loop galaxy power spectrum,
using measurements coming from the Flagship I simulation. In addition to this
analysis, we also tested our model for the matter power spectrum on a set of
measurements obtained on the fly while running the PKDGRAV3 code in the
redshift range $\left[0.7,2.4\right]$. Each measurement consists of the matter
power spectrum, measured in the range of wave modes defined by the interval
$\left[0.01,4\right]\,\,h\,{\rm Mpc}^{-1}$, using 18 linearly-spaced bins up
to $k\sim 0.03\,h\,{\rm Mpc}^{-1}$ and other 84 logarithmically-spaced bins
after.

We run two independent analyses, the first using the next-to-leading order
matter power spectrum obtained in the EFTofLSS framework (Eq. 17) and the
second using the nonlinear matter power spectrum from the baccoemu emulator.
We use an analytical Gaussian covariance matrix to describe the error on the
matter power spectrum measurements, assuming the full volume of the simulation
box, i.e., $(3780\,h^{-1}\,{\rm Mpc})^{3}$. For both models, we limit the
maximum mode of the fit to $k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1}$, within the
expected PT range of validity for the relevant redshifts.

Figure 13: Marginalised one-dimensional constraints as a function of redshift
obtained by fitting the measured matter power spectrum with the EFT model
(blue) and the baccoemu emulator (orange). Solid lines and shaded bands mark
the mean and the standard deviation of the posterior distribution,
respectively. For both models, the fit is carried out up to maximum wave mode
$k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1}$, using a Gaussian covariance matrix
corresponding to the full box volume of the Flagship I simulation. Dashed
lines denote the nominal fiducial values of the parameters
$\left\\{h,\omega_{\rm c},A_{\rm s}\right\\}$. In the bottom panel we show the
marginalised constraints on the $c_{0}$ EFT counterterm parameter.

In Fig. 13 we show the marginalised 1d posterior distributions for
$\left\\{h,\omega_{\rm c},A_{\rm s}\right\\}$ as a function of redshift. For
this test we limit the redshift range to the first four snapshots, i.e.,
$z<1.3$, since the results are sufficient to draw conclusions on the agreement
between data vector and theory models. Also shown, with dashed lines, are the
nominal values of the parameters provided in Potter et al. (2017). We find
that all cosmological parameters are obtained with a bias of 3–4$\sigma$ from
their fiducial values, in a way that is consistent across redshifts, as
highlighted by the almost constant trends in each of the panels. In addition,
the EFT model and baccoemu are consistent at better than 1$\sigma$ with each
other, pointing to a systematic effect that cannot be attributed to the
particular model used to describe the data vectors.

Figure 14: Ratio between the Flagship matter power spectrum measurements and
predictions obtained using baccoemu, considering both the the nominal
cosmology (upper panel) and the fiducial one (bottom panel). In both cases,
different colours correspond to a different redshift, as shown in the legend.
The grey shaded band mark the intrinsic error of the power spectrum
measurements to which an additional $1\%$ contribution has been added to
include the error contribution of the emulator, as explained in Sect. 5.1.

In the upper panel of Fig. 14 we show the ratio of the nonlinear matter power
spectrum from the baccoemu emulator assuming the nominal cosmology to the
measured one for the different redshifts. We notice how the discrepancy among
the two sets of curves is apparently due to a different tilt in the full-shape
of the matter power spectrum, here corresponding to the nominal value of
$n_{\rm s}=0.96$. Indeed we notice that the disagreement between nominal and
recovered cosmology can be alleviated for all parameters by assuming a
different value of $n_{\rm s}$, while keeping all the remaining parameters
fixed to their nominal value (with the exception of the scalar amplitude
$A_{s}$, since the latter has to be modified to recover the nominal value of
$\sigma_{8}$). This is also supported by the matter transfer function files
that were generated to set-up the initial conditions for the Flagship I
simulation, which confirmed how all the parameters affecting the shape of the
transfer function $\left\\{\omega_{\rm c},\,\omega_{\hbox{\hglue 0.5pt}\rm
b},\,M_{\nu}\right\\}$ are consistent with the nominal values listed in Table
1.

We therefore performed new fits to the Flagship I matter power spectra keeping
all parameters fixed (including the Hubble parameter $h$) to the nominal
values, and only sampling those controlling the primordial matter power
spectrum, i.e., the scalar amplitude $A_{\rm s}$ and index $n_{\rm s}$. The
best-fit value that we found for the spectral index is close to the value of
$n_{\rm s}=0.97$. This leads to consistent results across the redshift range
we are considering, as well as for a broad range of scales, up to $k\sim
4\,h\,{\rm Mpc}^{-1}$, as shown in the bottom panel of Fig. 14.

The new values of the spectral index and scalar amplitude, along with nominal
values for the other parameters constitute what we refer to as the fiducial
cosmology. This is the one given in Table 1 and adopted for all predictions
throughout this paper.

##  Appendix B Standard perturbation theory

In this section we report the full expressions for several quantities that are
defined in Sect. 3 and used throughout the rest of the paper. For a more
detailed and exhaustive description of the framework on which cosmological PT
is based, we refer the reader to the comprehensive review by Bernardeau et al.
(2002).

The main idea behind cosmological PT is that the generic solution to the
growth of nonlinear density and velocity fluctuations – $\delta(\bm{k})$ and
$\theta(\bm{k})$ – in an expanding universe can be expressed in terms of
linear theory solutions, $\delta_{\rm L}(\bm{k})$. Assuming an EdS universe,
it is possible to perfectly separate the time- and space-dependence of
$\delta(\bm{k})$ and $\theta(\bm{k})$ (Goroff et al. 1986; Jain & Bertschinger
1994). In an arbitrary $\Lambda$CDM cosmology, it can be shown that
approximate solutions can be found with the same separation (Donath & Senatore
2020), especially in a range of redshift for which the EdS cosmology is still
a valid approximation. In this case we can write the density and velocity
divergence fields using an expansion of the form

$\displaystyle\delta\,(\bm{k},\tau)=\sum_{n=1}^{\infty}D_{n}(\tau)\,\delta^{\,(n)}(\bm{k})\,,$
(45)
$\displaystyle\theta\,(\bm{k},\tau)=-\,\mathcal{H}(\tau)f(\tau)\sum_{n=1}^{\infty}D_{n}(\tau)\,\theta^{\,(n)}(\bm{k})\,,$
(46)

where $\tau$ is the conformal time defined via
$\mathrm{d}t=a(\tau)\,\mathrm{d}\tau$, $a(\tau)$ is the cosmic scale factor,
$\mathcal{H}(\tau)\equiv\mathrm{d}\ln{a(\tau)}/\mathrm{d}\tau$ is the
conformal Hubble expansion factor, and $f(\tau)\equiv\mathrm{d}\ln
D_{1}(\tau)/\mathrm{d}\ln a(\tau)$ is the growth rate. The $n$-th order growth
factor $D_{n}$ characterises the time-dependence of the density and velocity
field, and reduces to $D_{n}=a^{\,n}$ in the EdS limit. Assuming the
conservation of mass, momentum, and the Poisson equation, the individual
$n$-th order corrections to the density and velocity fields can be written as

$\displaystyle\begin{split}\delta^{\,(n)}(\bm{k})=&\int_{\bm{q}_{1}}\ldots\int_{\bm{q}_{n}}\delta_{\rm
D}^{\,(3)}(\bm{k}-\bm{q}_{1\ldots 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}$ (47)
$\displaystyle\begin{split}\theta^{\,(n)}(\bm{k})=&\int_{\bm{q}_{1}}\ldots\int_{\bm{q}_{n}}\delta_{\rm
D}^{\,(3)}(\bm{k}-\bm{q}_{1\ldots n})\,G_{n}(\bm{q}_{1},\ldots,\bm{q}_{n})\\\
&\times\delta_{\rm L}(\bm{q}_{1})\ldots\delta_{\rm
L}(\bm{q}_{n})\,,\end{split}$ (48)

where the $n$-th order PT kernels $F_{n}$ and $G_{n}$ are homogeneous
functions of the wave vectors $\left(\bm{q}_{1},\ldots,\bm{q}_{n}\right)$, and
are built starting from the fundamental mode-coupling functions,

$\displaystyle\alpha\,(\bm{k}_{1},\bm{k}_{2})\equiv\frac{\bm{k}_{12}\cdot\bm{k}_{1}}{k_{1}^{\,2}}\,,$
(49)
$\displaystyle\beta\,(\bm{k}_{1},\bm{k}_{2})\equiv\frac{k_{12}^{\,2}(\bm{k}_{1}\cdot\bm{k}_{2})}{2\,k_{1}^{\,2}k_{2}^{\,2}}\,,$
(50)

with $\bm{k}_{12}=\bm{k}_{1}+\bm{k}_{2}$. At linear order these quantities
clearly become unity, $F_{1}=G_{1}=1$, so to recover linear theory
predictions, e.g., $\delta^{\,(1)}(\bm{k})=\theta^{\,(1)}(\bm{k})=\delta_{\rm
L}(\bm{k})$. At higher order, these kernels can be derived using recursive
relations, which read

$\displaystyle F_{n}\,\left(\bm{q}_{1},\ldots,\bm{q}_{n}\right)=$
$\displaystyle\sum_{m=1}^{n-1}\frac{G_{m}(\bm{q}_{1},\ldots,\bm{q}_{m})}{(2n+3)\,(n-1)}$
$\displaystyle\times\Big{[}(2n+1)\,\alpha(\bm{k}_{1},\bm{k}_{2})\,F_{n-m}\,(\bm{q}_{m+1},\ldots,\bm{q}_{n})$
$\displaystyle\hskip
11.0pt+2\,\beta(\bm{k}_{1},\bm{k}_{2})\,G_{n-m}(\bm{q}_{m+1},\ldots,\bm{q}_{n})\Big{]}\,,$
(51) $\displaystyle G_{n}\,\left(\bm{q}_{1},\ldots,\bm{q}_{n}\right)=$
$\displaystyle\sum_{m=1}^{n-1}\frac{G_{m}(\bm{q}_{1},\ldots,\bm{q}_{m})}{(2n+3)\,(n-1)}$
$\displaystyle\times\Big{[}3\,\alpha(\bm{k}_{1},\bm{k}_{2})\,F_{n-m}(\bm{q}_{m+1},\ldots,\bm{q}_{n})$
$\displaystyle\hskip
11.0pt+2n\,\beta(\bm{k}_{1},\bm{k}_{2})\,G_{n-m}(\bm{q}_{m+1},\ldots,\bm{q}_{n})\Big{]}\,,$
(52)

where $\bm{k}_{1}=\bm{q}_{1}+\ldots+\bm{q}_{m}$ and
$\bm{k}_{2}=\bm{q}_{m+1}+\ldots+\bm{q}_{n}$. As a classical example, needed
for the calculation of the one-loop galaxy power spectrum, the second-order PT
kernels for the nonlinear evolution of the matter density and velocity fields
are defined as

$\displaystyle
F_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{5}{7}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)+\frac{2}{7}\frac{\left(\bm{k}_{1}\cdot\bm{k}_{2}\right)^{2}}{k_{1}^{\,2}\,k_{2}^{\,2}}\,,$
(53) $\displaystyle
G_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{3}{7}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)+\frac{4}{7}\frac{\left(\bm{k}_{1}\cdot\bm{k}_{2}\right)^{2}}{k_{1}^{\,2}\,k_{2}^{\,2}}\,,$
(54)

while the explicit expression for the third-order kernel of the matter density
field, $F_{3}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$, can be found in Goroff et
al. (1986), and it is not reported here for practical purposes.

When considering the galaxy power spectrum $P_{\rm gg}$, we need to evaluate
the usual two-point statistics defined by the ensemble average of the galaxy
density field times itself,

$\left\langle\delta_{\rm g}(\bm{k})\,\delta_{\rm
g}(\bm{k}^{\prime})\right\rangle=(2\pi)^{3}\,\delta_{\rm
D}^{\,(3)}(\bm{k}+\bm{k}^{\prime})\,P_{\rm gg}(k)\,.$ (55)

At third order in the perturbations of $\delta$, we obtain the one-loop
expression for the power spectrum presented in Eq. (25), where all the next-
to-leading order corrections are grouped into mode-coupling and propagator-
like contributions. In the former, the loop integrand is proportional to
$P_{\rm L}(\left|\bm{k}-\bm{q}\right|)\,P_{\rm L}(q)$, reflecting the mixing
of modes due to nonlinear evolution, while in the latter the integrals are
carried out on the factor $P_{\rm L}(k)\,P_{\rm L}(q)$, corresponding to a
time propagation of the initial density field. Separating these two groups
into individual corrections, each one multiplied by a given combination of
bias parameters, we end up with the following scheme,

$\begin{split}P_{\rm gg}^{\,\rm 1\mbox{-}loop}(k)=\>&b_{1}^{\,2}\,P^{\,\rm
1\mbox{-}loop}(k)\\\ &+b_{1}b_{2}\,P_{b_{1}b_{2}}(k)+b_{1}b_{{\cal
G}_{2}}\,P_{b_{1}b_{{\cal G}_{2}}}(k)\\\
&+b_{1}b_{\Gamma_{3}}\,P_{b_{1}b_{\Gamma_{3}}}(k)+b_{2}^{\,2}\,P_{b_{2}b_{2}}(k)\\\
&+b_{2}b_{{\cal G}_{2}}\,P_{b_{2}b_{{\cal G}_{2}}}(k)+b_{{\cal
G}_{2}}^{\,2}\,P_{b_{{\cal G}_{2}}b_{{\cal G}_{2}}}(k)\,,\end{split}$ (56)

where all the previous terms can be represented as loop integrals,

$\begin{split}P^{\,\rm 1\mbox{-}loop}(k)=\>&P^{\,\rm
1\mbox{-}loop,\,MC}(k)+P^{\,\rm 1\mbox{-}loop,\,Prop}(k)\\\
=\>&2\int_{\bm{q}}F_{2}^{\;2}(\bm{k}-\bm{q},\bm{q})\,P_{\rm
L}(\left|\bm{k}-\bm{q}\right|)\,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}$
(57) $P_{b_{1}b_{2}}(k)=2\int_{\bm{q}}F_{2}(\bm{k}-\bm{q},\bm{q})\,P_{\rm
L}(\left|\bm{k}-\bm{q}\right|)\,P_{\rm L}(q)\,,$ (58)
$\begin{split}P_{b_{1}b_{{\cal G}_{2}}}(k)=\>&P_{b_{1}b_{{\cal G}_{2}}}^{\,\rm
MC}(k)+P_{b_{1}b_{{\cal G}_{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}(\left|\bm{k}-\bm{q}\right|)\,P_{\rm L}(q)\\\ &+8\,P_{\rm
L}(k)\int_{\bm{q}}F_{2}(\bm{k},-\bm{q})\,S(\bm{k}-\bm{q},\bm{q})\,P_{\rm
L}(q)\,,\end{split}$ (59) $P_{b_{1}b_{\Gamma_{3}}}(k)=-\frac{16}{7}P_{\rm
L}(k)\int_{\bm{q}}S(\bm{k}-\bm{q},\bm{q})\,S(\bm{k},\bm{q})\,P_{\rm L}(q)\,,$
(60) $P_{b_{2}b_{2}}(k)=\frac{1}{2}\int_{\bm{q}}P_{\rm
L}(\left|\bm{k}-\bm{q}\right|)\,P_{\rm L}(q)\,,$ (61) $P_{b_{2}b_{{\cal
G}_{2}}}(k)=2\int_{\bm{q}}S(\bm{k}-\bm{q},\bm{q})\,P_{\rm
L}(\bm{k}-\bm{q},\bm{q})\,P_{\rm L}(q)\,,$ (62) $P_{b_{{\cal G}_{2}}b_{{\cal
G}_{2}}}(k)=2\int_{\bm{q}}S^{2}(\bm{k}-\bm{q},\bm{q})\,P_{\rm
L}(\bm{k}-\bm{q},\bm{q})\,P_{\rm L}(q)\,.$ (63)

In the previous set of equations, the $P_{b_{1}b_{\Gamma_{3}}}$ contribution
is characterised by a single propagator-like term, which is perfectly
degenerate with the second addend contributing to $P_{b_{1}b_{{\cal G}_{2}}}$.
Following the expansion that has been adopted in this paper, the degeneracy
results from the equality

$P_{b_{1}b_{\Gamma_{3}}}(k)=\frac{2}{5}P_{b_{1}b_{{\cal G}_{2}}}^{\,\rm
MC}(k)\,.$ (64)

As shown in the main body of this article, breaking this degeneracy with the
information coming from the galaxy power spectrum alone is not possible, and
for this reason we often keep one of the two tidal bias parameters fixed to
some physically-motivated relation, such as the excursion-set-based relation
(Eq. 35) for $b_{{\cal G}_{2}}$, or the coevolution relations for both
$b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$ (Eqs. 33 – 34).

The behaviour of all the loop integrals listed above is to consistently
converge to zero at infrared modes, since the nonlinear kernel
$F_{n}(\bm{k}_{1},\ldots,\bm{k}_{n})$ scales as $k^{2}$ when
$\bm{k}\equiv\bm{k}_{1}+\ldots+\bm{k}_{n}$ goes to zero, reflecting the range
of validity of linear theory predictions. The only exception is represented by
the $P_{b_{2}b_{2}}$ term, which features a non-zero asymptote for
$k\rightarrow 0$. This limit can be manually set to zero via a redefinition of
the loop integral, such that

$P_{b_{2}b_{2}}(k)=\frac{1}{2}\int_{\bm{q}}P_{\rm L}(q)\,\left[P_{\rm
L}(\left|\bm{k}-\bm{q}\right|)-P_{\rm L}(q)\right]\,.$ (65)

In turn, the extra contribution

$P_{b_{2}b_{2}}^{\,\rm noise}(k)=\int_{\bm{q}}P_{\rm L}^{\,2}(q)$ (66)

can be absorbed by the constant shot-noise parameter, $\alpha_{P,1}$, that we
have defined in Sect. 3.1.2, as they both correspond to constant shift in the
amplitude of the galaxy power spectrum.

##  Appendix C Implementation of the wiggle vs no-wiggle split

In this section, we investigate the prescriptions used to obtain a smooth
template $P_{\rm nw}$ starting from the linear matter power spectrum. This is
an important aspect of the theoretical recipe that we adopt, as the wiggle-no
wiggle split is essential for the correct implementation of IR-resummation, as
shown in Sect. 3.1. While several different algorithms can be found in the
literature, here we test three different methods.

The first one is based on a one-dimensional Gaussian smoothing (GS1D), and
consists in a rescaling of the original formula for the featureless matter
power spectrum $P_{\rm EH}$, originally presented in Eisenstein & Hu (1998),
to match the broadband amplitude of the linear matter power spectrum. In
practice, we follow the approach of Vlah et al. (2016), who defines the smooth
component of the linear matter power spectrum as

$P_{\rm nw}(k)=P_{\rm EH}(k)\,{\cal F}\left[\frac{P_{\rm L}(k)}{P_{\rm
EH}(k)}\right]\,,$ (67)

where ${\cal F}$ is meant to filter out the broadband difference between
$P_{\rm L}$ and $P_{\rm EH}$. We choose a functional form for ${\cal F}$
corresponding to a Gaussian filter, that is,

${\cal F}\left[f(k)\right]=\;\frac{\operatorname{log_{10}}({\rm
e})}{\sqrt{2\pi}\lambda}\int\mathrm{d}q\,\frac{f(q)}{q}\exp\left[-\frac{1}{2\lambda^{2}}\log^{\,2}_{10}\left(\frac{k}{q}\right)\right]\,,$
(68)

where $\lambda$ determines the variance of the Gaussian filter used to rescale
the ratio of the linear to the featureless power spectrum. In this analysis,
we fix its value to $\lambda=0.25$.

The result of this approach is presented in the top panel of Fig. 15, where we
show a comparison between the shape of the original Eisenstein & Hu (1998)
smooth function and the one presented in Eq. (67) using the fiducial cosmology
from Table 1 at $z=0$. From the comparison, it is clear that the broadband of
$P_{\rm EH}$ features a non-negligible tilt with respect to the one of the
linear matter power spectrum, reaching deviations of up to $2\%$ across the
whole BAO wave mode interval. This difference can be corrected using the
filtering strategy, whose output is completely consistent with the broadband
shape of $P_{\rm L}$, with the only exception of a small residual of
${\sim}\,0.5\%$ peaking at $0.5\,h\,{\rm Mpc}^{-1}\lesssim k\lesssim
1\,h\,{\rm Mpc}^{-1}$.

Figure 15: _Top_ : comparison between the raw featureless power spectrum
$P_{\rm EH}$ and the no-wiggle power spectrum $P_{\rm nw}$ computed using Eq.
(67). In both cases, the power spectra are computed at $z=0$ using the
fiducial cosmology from Table 1. _Bottom_ : comparison between the three
different methods we tested to obtain the no-wiggle power spectrum $P_{\rm
nw}$. The thick orange line corresponds to the method we selected, i.e., the
convolution with a Gaussian smoothing function, while the other two lines
represent a Discrete Sine Transform (red) and a basis spline (blue). In both
panels, the grey shaded bands represent the $1\%$ (dark grey) and $2\%$ (light
grey) thresholds.

The second approach makes use of a discrete sine transform (DST), and has been
originally proposed in Hamann et al. (2010) (see Chudaykin et al. 2020 and
Ivanov et al. 2020 for an application to real data). This consists in a fast
Fourier transformation of the input matter power spectrum, and in the removal
of the bump corresponding to the BAO peak. This step is carried out based on
the value of the second derivative of the sine transform, and the generated
gap is subsequently filled using a cubic-spline interpolation. The new
function is finally transformed back into Fourier space to deliver a power
spectrum shape deprived from BAO oscillations.

The third and final approach is based on the approximation of the BAO wiggles
with a basis spline (B-spline) curve, starting from a set of knots
$\\{k_{i},P_{i}\\}$, and subsequently finding the spline coefficients that
maximise the likelihood with the original power spectrum (see, e.g., Vlah et
al. 2016, for a similar implementation).

The bottom panel of Fig. 15 shows a comparison between the three methods
summarised above. We notice a non-negligible difference, with discrepancies on
the BAO scales that can reach a fraction of percent, depending on the
considered approach (see Moradinezhad Dizgah et al. 2021, for similar
conclusions). Nevertheless, at linear order these differences are diluted by
the recombination of the smooth and wiggling component at a later stage, so
that the net result on the IR-resummed nonlinear galaxy power spectrum is
going to be much smaller than the values exhibited in this plot. We highlight
a major discrepancy between the DST method and the other two that reaches its
maximum at a scale $k\sim 0.01\,h\,{\rm Mpc}^{-1}$, which roughly coincides
with the position of the turnaround in the matter power spectrum. Because of
this major discrepancy, and since all three methods behave in a slightly
different way on the BAO scales, we decided to adopt the GS1D method
throughout the analysis presented in this paper. We observe a marginal
residuals of this method on mildly nonlinear scales, at $0.5\,h\,{\rm
Mpc}^{-1}\lesssim k\lesssim 1\,h\,{\rm Mpc}^{-1}$, which, once again, is not
likely to significantly bias our results

At the same time we highlight how this approach is most likely going to
perform worse when considering cosmologies beyond the vanilla $\Lambda$CDM
model, such as those including massive neutrinos. In these cases, depending on
the magnitude of the deviation from $\Lambda$CDM, the broadband shape of
$P_{\rm EH}$ can deviate from that of $P_{\rm L}$ by up to $10\%$ (as quoted
in Eisenstein & Hu 1998). We will therefore adopt one of the other two methods
for future analyses, e.g.,the upcoming analysis of the Flagship II simulation.

##  Appendix D Comparison between model implementation

Figure 16: Comparison among the different methods to obtain the one-loop bias
integrals, colour-coded according to the individual code/method, as listed in
the legend. The reference code for these plots is the COMET emulator.
Different panels corresponds to different diagrams: the first column shows the
leading order IR-resummed matter power spectrum – for which no actual
integration is required – and the next-to-leading order correction. The
remaining columns show the bias diagrams from Eqs. (58 – 63).

In order to estimate the systematic error budget due to the implementation of
the algorithm for the one-loop model presented in Sect. 3.1 and Appendix B, we
make use of four individual implementations, based on four independent codes
provided by several group members. Each code features a different way to
compute the loop corrections presented in Sect. B.

These include a two-dimensional integration implemented within the Cuba
library 242424http://www.feynarts.de/cuba/ (Hahn 2005) used to generate the
training set of the COMET emulator (Eggemeier et al. 2022), and implemented in
the CosmoSIS-gClust 252525To be made publicly available soon. code (courtesy
of A. Moradinezhad Dizgah). Alternative methods take advantage of a Fast
Fourier Transform approach, such as FastPT
262626https://github.com/JoeMcEwen/FAST-PT/ (McEwen et al. 2016) – implemented
in the PBJ code (Oddo et al. 2020, 2021; Rizzo et al. 2023; Moretti et al.
2023) – and FFTLog 272727https://jila.colorado.edu/~ajsh/FFTLog/ (Hamilton
2000; Simonović et al. 2018) – used in the TNSToolBox
282828https://github.com/sdlt/TNS_ToolBox/ code (courtesy of S. de la Torre)
and again in CosmoSIS-gClust.

In Fig. 16 we show the systematic deviation between different computations of
the same terms, including the leading- and next-to-leading order correction to
the IR-resummed matter power spectrum, and the six one-loop bias corrections
defined in Eqs. (58 – 63). For this exercise, we compute a common data vector
using the fiducial Flagship cosmology from Table 1 at $z=0$. We assume as a
reference measurement the one performed using a direct integration of the wave
number $q$ within the range of scales $\left[0.00001,100\right]\,h\,{\rm
Mpc}^{-1}$, which proved to be a sensible choice to achieve the convergence of
the different loop integrals. Finally, we compare the different terms
evaluated using the four codes described in the previous paragraph. We find an
overall optimal consistency among the different theory implementations, with
all terms being in agreement at better than $0.01\%$ on the overall range of
scales shown in Fig. 16. We observe a slightly worse concordance between
different computations of the propagator-like terms, that is, $P_{\rm
mm}^{\,{\rm IR\mbox{-}NLO}}$, which contains the contribution $P_{13}$, and
more importantly $P_{b_{1}b_{\Gamma_{3}}}$, for which the discrepancy can
become as large as ${\sim}\,0.05\%$ on scales of $k\sim 1\,h\,{\rm Mpc}^{-1}$.
However, we consider this difference completely negligible, since this value
is much smaller than the statistical precision of the data vectors, and since
this scale is completely outside of the range of scales considered in this
work.

##  Appendix E Comparison of different samplers

As anticipated in Sect. 4, in order to obtain a posterior distribution for
each of the model configurations that we test, we need to select a robust
algorithm to sample the large multi-dimensional parameter space of the models
described in Sect. 3. Indeed, the choice of the sampler is critical in
presence of multi-modal distributions. As already thoroughly discussed in the
main body of this paper, this situation is prevalent when we focus solely on
the full shape of the galaxy power spectrum, using only two-point statistics,
due to the strong degeneracies between model parameters. In these situations
traditional algorithms, which are meant to explore smaller and more well-
behaved – Gaussian-like – random variables, may be under-performing, and
therefore affecting the efficiency of the sampling.

In this analysis we compare three different samplers, testing them against a
reference subset of the galaxy power spectrum data vectors already used to
validate the theoretical model for one-loop galaxy bias. The properties of
these sampling algorithms are listed hereafter.

###  Metropolis–Hastings sampler

This approach (Metropolis et al. 1953; Hastings 1970) is a Markov-Chain Monte
Carlo (MCMC) method based on the construction of a random walk inside the
parameter space. Subsequent points in the Markov chain are determined based on
a proposal function that has to be specified as a free parameter of the model.
For each new point, the likelihood of the model determines whether the
candidate state is accepted or discarded, i.e., whether to move to the
candidate state or stay in the current state. For this standard algorithm we
make use of a non-public code,292929COMPASS, courtesy of Ariel G. Sánchez,
Martín Crocce, and Román Scoccimarro. which features an implementation of a
Metropolis-Hastings algorithm coupled with the likelihood for the galaxy power
spectrum described in Sect. 3.

Figure 17: Comparison between runs at fixed cosmology obtained using different
samplers, carried out using a reference galaxy power spectrum (Model 3 HOD
sample at $z=0.9$). The left panel shows the posterior distribution of the
nuisance parameters with $b_{{\cal G}_{2}}$ fixed to the local Lagrangian
relation (Eq. 33), while in the right panel all the parameters are left free
to vary. Different colours correspond to different samplers, as listed in the
legend.

###  Affine-invariant sampler

This MCMC approach (Goodman & Weare 2010) is still based on a random walk
across the parameter space, but with an improved ensemble sampler, in which a
large number of walkers interact with each other in a way that reduces the
dependence of the sampling on the aspect ratio of the particular posterior
distribution under consideration. We use the affine-invariant sampler
implemented in the public Python package emcee (Foreman-Mackey et al.
2013).303030https://emcee.readthedocs.io/en/stable/

### Nested sampling

This approach (Skilling 2006) relies on iteratively refining a set of live
points, initially drawn from a prior distribution, to explore the parameter
space. At each iteration, the point with the lowest likelihood is selected
from the set and replaced with a new point. The latter is drawn from a
constrained region of the prior distribution, where the likelihood must be
higher than that of the replaced point. The selected point is then assigned a
weight, based on its likelihood and the proportion of the prior volume it
represents: several iterations will allow us to thus construct the posterior
distribution. These algorithms are particularly suited for multi-modal
distributions, given they are not subject to getting stuck in local minima of
the loglikelihood, as it commonly happens with standard algorithms based on
Markov chains. We use the public package Multinest (Feroz et al. 2009, 2019),
which can be interfaced with Python using a dedicated wrapper
module.313131https://johannesbuchner.github.io/PyMultiNest/

### Comparison

The comparison between the marginalised two-dimensional distributions obtained
with the different samplers is shown in Fig. 17, using as data the galaxy
power spectrum of the Model 3 HOD sample at $z=0.9$ with $k_{\rm
max}=0.2\,h\,{\rm Mpc}^{-1}$. We test two different cases, one where we vary
the nuisance parameters
$\left(b_{1},b_{2},b_{\Gamma_{3}},c_{0},\alpha_{P,1}\right)$ while keeping the
quadratic tidal bias $b_{{\cal G}_{2}}$ fixed to its local Lagrangian relation
(33), and one where the latter is also allowed to freely vary with the rest of
the parameters.323232As shown in the main body of this article, fixing
$b_{{\cal G}_{2}}$ is the only case for which it is possible to completely
break the strong degeneracies of the parameter space.

In the first case, we notice an almost perfect match among the three different
set of contours, with a statistically-negligible difference only appreciable
at the tails of the posterior distributions, which can be anyhow partially
explained by the intrinsic variance of the individual realizations of the
posterior distributions. The second case shows a slightly larger discrepancy,
mostly due to the presence of a multi-modal profile driven by the strong
degeneracy between $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$. In this case, the
peak of the posterior distribution is less consistent when using different
samplers, and we find that the size of the discrepancy is a very strong
function of $k_{\rm max}$. In details, the main differences arise because of
the non-homogeneous sampling of different peaks, for which traditional Markov
chain algorithms may get stuck in a particular minimum of the distribution,
therefore affecting the overall convergence of the chain. Given the purpose
that nested sampling algorithms were initially developed for, we decided to
employ the latter as our baseline sampler, and used it to sample the parameter
space for all the results presented in this work.

Figure 18: Marginalised 1d constraints on the linear bias $b_{1}$ and the
shot-noise parameter $\alpha_{P,1}$ obtained from a set of synthetic theory
vectors created using the same recipe from Eq. (32), for the four different
redshifts already explored with the Flagship data vectors. Different colours
correspond to different assumptions on the total number of degrees of freedom
of the model, as listed in the legend. Figure 19: Same as in Fig. 18, but
using theory vectors displaced by Gaussian noise realisation from a box of
size $3780\,h^{-1}\,{\rm Mpc}$.

##  Appendix F Sample variance effects

Even though the volume of the Flagship I simulation is quite large, most of
our results are, to some extent, affected by the sample variance of the single
N-body realization available.

In this section we quantify this effect reproducing some of the parameters
fits starting from noiseless, synthetic data vectors obtained as the real-
space galaxy power spectrum predictions from the COMET emulator. These are
generated adopting the reference cosmological parameters from Table 1 at the
four different redshifts of the comoving snapshots,
$z=\left\\{0.9,1.2,1.5,1.8\right\\}$. For each of the four samples, the
nuisance parameters are derived as follows: $b_{1}$ and $\alpha_{P,1}$ are
fixed to the best-fit value listed in Table 4 for the Model 1 HOD samples,
$b_{2}$ is computed using a $b_{2}(b_{1})$ relation derived from the
corresponding HOD model, $b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$ are obtained
from the excursion-set relation (35) and the coevolution relation (34),
respectively, $c_{0}$ is set to unity, and $\alpha_{P,2}$ is set to zero. The
covariance matrix for each sample is computed assuming only the Gaussian
component from the full Flagship I box, and therefore corresponds to a volume
of about $58\,h^{-3}\,{\rm Gpc}^{3}$.

Figure 18 shows the marginalised one-dimensional constraints of the $b_{1}$
and $\alpha_{P,1}$, similarly to Fig. 5, for fits at fixed cosmology of the
synthetic data vectors. Even though the data vectors are smooth – no noise
component has been added to any of the $P_{\rm gg}(k_{i})$ bins – and also
perfectly consistent with the theory model used to fit them, we observe some
discrepancy of the constraints with the fiducial values, both with all
parameters free as for the the case with $b_{\Gamma_{3}}$ fixed to the
coevolution relation. In both cases, the statistical constraints are less
tight than the ones obtained with other bias relations, with a partial break
of the degeneracies between parameters at higher redshift then the coevolution
relation is applied. Such deviations are induced by projection – or prior
volume – effects in the marginalised posteriors, due to the strong
degeneracies between parameters. In particular, these affect $b_{2}$,
$b_{{\cal G}_{2}}$ and $b_{\Gamma_{3}}$, and cannot be broken completely even
when reducing the dimensionality of the parameter space by fixing
$b_{\Gamma_{3}}$, as the $b_{2}$-$b_{{\cal G}_{2}}$ degeneracy is still
present. Consistently among the different redshifts, we notice how there is a
trend for $b_{1}$ and $\alpha_{P,1}$ to be under- and overestimated,
respectively, similarly to what we observe for the real Flagship I data
vectors.

On the contrary, the case where $b_{{\cal G}_{2}}$ is fixed as a function of
$b_{1}$ is systematically better in terms of amplitude of the errorbars and
accuracy in the recovery of the parameters, because the posterior distribution
of this case is closer to a Gaussian distribution. This is due to the
simultaneous breaking of both the $b_{2}$-$b_{{\cal G}_{2}}$ and $b_{{\cal
G}_{2}}$-$b_{\Gamma_{3}}$ degeneracies, for which a clear example can be found
in the left plot of Fig. 17.

As a follow-up test, we add to the synthetic data vectors a Gaussian noise
consistent with the covariance assumed for the fits to the simulation. Figure
19 shows the new marginalised constraints on $b_{1}$ and $\alpha_{P,1}$. In
this case it is possible to observe a much larger discrepancy from the
fiducial values, present also for the configurations in which $b_{{\cal
G}_{2}}$ is fixed as a function of $b_{1}$. We remark how the case at $z=1.2$
features a large fluctuation at intermediate scales, $k_{\rm max}\sim
0.2\,h\,{\rm Mpc}^{-1}$, which is symptomatic of the particular realization of
the Gaussian noise. We repeated this exercise, on the sample at $z=1.2$, with
ten different noise realizations, and in all cases we observe a slightly
different trend as a function of $k_{\rm max}$ for the marginalised
constraints of the configurations we are considering, i.e., with all the
parameters free to vary, or with either $b_{{\cal G}_{2}}$ or $b_{\Gamma_{3}}$
fixed in terms of $b_{1}$. For some of the configurations we observe trends as
the one seen in the $z=0.9$ Flagship I data vectors, for which we find a sharp
running in the $(b_{1},\alpha_{P,1})$ plane at $k_{\rm max}\gtrsim
0.2\,h\,{\rm Mpc}^{-1}$. We therefore conclude that a more realistic analysis,
also including observational effects and a proper survey window function
should be conducted using a wider set of simulations, in order to reduce the
overall impact of cosmic variance.