<figure><img src="image_18.png" /><figcaption>**Fig. C.1.** Ratio between the no-wiggle and the linear matter power spectrum. Top: Comparison between the raw featureless power spectrum $P_{EH}$ and the no-wiggle power spectrum $P_{nw}$ computed using Eq. (C.1). 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_{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.</figcaption></figure>

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 non-linear 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 \text{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 non-linear scales, at $0.5 h \text{Mpc}^{-1} \le k \le 1 h \text{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 ΛCDM model, such as those including massive neutrinos. In these cases, depending on the magnitude of the deviation from ΛCDM, the broadband shape of $P_{EH}$ can deviate from that of $P_L$ by up to 10% (as quoted in Eisenstein & Hu 1998). We therefore adopt one of the other two methods for future analyses, such as for the upcoming analysis of the Flagship II simulation.

## Appendix D: Comparison between model implementation

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<sup>30</sup> (Hahn 2005) used to generate the training set of the COMET emulator (Eggemeier et al. 2022), and implemented in the CosmoSIS-gClust<sup>31</sup> code (courtesy of A. Moradinezhad Dizgah). Alternative methods take advantage of a Fast Fourier Transform approach, such as FastPT<sup>32</sup> (McEwen et al. 2016) – implemented in the PBJ code (Oddo et al. 2020, 2021; Rizzo et al. 2023; Moretti et al. 2023) – and FFTLog<sup>33</sup> (Hamilton 2000; Simonović et al. 2018) – used in the TNSToolBox<sup>34</sup> code (courtesy of S. de la Torre) and again in CosmoSIS-gClust.

In Fig. D.1 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. (B.14–B.19). 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 [0.00001, 100] $h \text{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. D.1. We observe a slightly worse concordance between different computations of the propagator-like terms, that is, $P_{\text{mm}}^{\text{IR-(LO+NLO)}}$, which contains the contribution $P_{13}$, and more importantly $P_{b_1 b_3}$, for which the discrepancy can become as large as ~0.05% on scales of $k \sim 1 h \text{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.

<footnote n="30"><sup>30</sup> http://www.feynarts.de/cuba/</footnote>

<footnote n="31"><sup>31</sup> To be made publicly available soon.</footnote>

<footnote n="32"><sup>32</sup> https://github.com/JoeMcEwen/FAST-PT/</footnote>

<footnote n="33"><sup>33</sup> https://jila.colorado.edu/~ajsh/FFTLog/</footnote>

<footnote n="34"><sup>34</sup> https://github.com/sdlt/TNS ToolBox/</footnote>
<figure><img src="image_19.png" /><figcaption>Fig. D.1. 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. (B.14–B.19).</figcaption></figure>

## Appendix E: Comparison between posterior averaged and maximum likelihood position

In the main body of this article, we determine the goodness of fit of a given model configuration in terms of the posterior-averaged $\chi^2$ statistics. For this purpose, we simply iterate through the sampled positions of the posterior distribution and define the weighted mean as

$$ \langle \chi^2 \rangle = \sum_i w_i \chi_i^2, \qquad (E.1) $$

where $\chi_i^2$ and $w_i$ correspond to the $\chi^2$ and corresponding weight for the i-th sampled position, respectively. As already mentioned in Sec. 4.1.2, this quantity can be estimated with less uncertainty from a sampled posterior distribution with respect to the maximum-likelihood position. However, the conclusions drawn in Sec. 5 are independent from which definition is used to quantify the goodness-of-fit statistics.

In Fig. E.1 we show the goodness-of-fit performance metric for the case of the Model 3 sample at $z = 1.2$, which is one of the few cases with a posterior-averaged $\chi^2$ getting above the $2\sigma$ in Fig. 9 and 10. Differently from those two plots, in this case we compare the trend when using the posterior-averaged $\chi^2$ (solid lines) and maximum-likelihood value (dash-dotted lines). As expected, the latter is constantly outperforming the former, with a relative improvement that becomes larger at low $k_{max}$ values. Moving to larger $k_{max}$, we note that for this specific case the maximum-likelihood position makes all the points consistent within $2\sigma$ of the considered $\chi^2$ distribution. However, the

<figure><img src="image_20.png" /><figcaption>Fig. E.1. Normalised $\chi^2$ values for the Model 3 sample at $z = 1.2$. Solid and dot-dashed lines corresponds to the posterior-averaged and minimum $\chi^2$, respectively. Different colours represent different bias modelling assumptions, as listed in the legend. The two grey shaded areas represent the 68% and 95% of the $\chi^2$ distribution with the same number of degrees of freedom of the considered model.</figcaption></figure>

marginal gains are not significantly changing the interpretations that can be made when using the posterior-averaged values.
<figure><img src="image_21.png" /><figcaption>**Fig. E.2.** Normalised $\chi^2$ distribution for the Model 3 sample at $z = 1.2$ and four different value of $k_{\text{max}}$, as shown in the upper right corner of the individual panels. The shaded grey areas mark the same confidence intervals as in Fig. E.1. The magenta vertical line marks the posterior-averaged value.</figcaption></figure>

As a complement, in Fig. E.2 we show the complete distribution of $\chi^2$ values, and we compare it with the 68th and 95th percentile of the reference $\chi^2$ distribution with the same numbers of degrees of freedom. Different panels correspond to different $k_{\text{max}}$ values, while the solid magenta line marks the position of the posterior-averaged value. The only configurations for which there is a partially significant gain in picking the maximum-likelihood value is clearly the one at $k_{\text{max}} = 0.4 h \, \text{Mpc}^{-1}$. Otherwise, our conclusions are not affected by the choice of which $\chi^2$ value is used to estimate the goodness-of-fit performance metric. For this reason, in Sec. 5 we show results employing the posterior-averaged $\chi^2$ statistics.

## Appendix F: Dependency on the selected sampler

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, that is, 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,<fnref n="35" /> which features an implementation of a Metropolis-Hastings algorithm coupled with the likelihood for the galaxy power spectrum described in Sect. 3.

### 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).<fnref n="36" />

### 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. 2009a, 2019b), which can be interfaced with Python using a dedicated wrapper module.<fnref n="37" />

### Comparison

The comparison between the marginalised two-dimensional distributions obtained with the different samplers is shown in Fig. F.1, using as data the galaxy power spectrum of the Model 3 HOD sample at $z = 0.9$ with $k_{\text{max}} = 0.2 h \, \text{Mpc}^{-1}$. We test two different cases, one where we vary the nuisance parameters ($b_1, b_2, b_{\Gamma_1}, c_0, \alpha_{P,1}$) while keeping the quadratic tidal bias $b_{G_2}$ fixed to its local Lagrangian relation (35), and one where the latter is also allowed to freely vary with the rest of the parameters.<fnref n="38" />

In the first case, we note an almost perfect match among the three different set of contours, with a statistically negligi-

<footnote n="35">³⁵ COMPASS, courtesy of Ariel G. Sánchez, Martín Crocce, and Román Scocimarro.</footnote>

<footnote n="36">³⁶ https://emcee.readthedocs.io/en/stable/</footnote>

<footnote n="37">³⁷ https://johannesbuchner.github.io/PyMultiNest/</footnote>

<footnote n="38">³⁸ As shown in the main body of this article, fixing $b_{G_2}$ is the only case for which it is possible to completely break the strong degeneracies of the parameter space.</footnote>
<figure><img src="image_22.png" /><figcaption>**Fig. F.1.** 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_{\tilde{G}_2}$ fixed to the local Lagrangian relation (Eq. 35), while in the right panel all the parameters are left free to vary. Different colours correspond to different samplers, as listed in the legend.</figcaption></figure>

ble 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_{\tilde{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_{\text{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.

## Appendix G: 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 = \{0.9, 1.2, 1.5, 1.8\}$. 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_{\tilde{G}_2}$ and $b_{\Gamma_3}$ are obtained from the excursion-set relation (37)

and the coevolution relation (36), 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} \text{ Gpc}^3$.

Figure G.1 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_{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 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_{\tilde{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_{\tilde{G}_2}$ degeneracy is still present. Consistently among the different redshifts, we note that 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_{\tilde{G}_2}$ is fixed as a function of $b_1$ is systematically better in terms of amplitude of the error bars 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_{\tilde{G}_2}$ and $b_{\tilde{G}_2}-b_{\Gamma_3}$ degeneracies, for which a clear example can be found in the left plot of Fig. F.1.
<figure><img src="image_23.png" /><figcaption>**Fig. G.1.** 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. (34), 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.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. G.2.** Same as in Fig. G.1, but using theory vectors displaced by Gaussian noise realisation from a box of size $3780 \ h^{-1} \ \text{Mpc}$.</figcaption></figure>

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 G.2 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_{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_{\text{max}} \sim 0.2 \ h \ \text{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_{\max}$ for the marginalised con-
straints of the configurations we are considering, that is, with
all the parameters free to vary, or with either $b_{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 vec-

tors, for which we find a sharp running in the $(b_1, \alpha_{P,1})$ plane
at $k_{\max} \gtrsim 0.2\,h\,\mathrm{Mpc}^{-1}$. We therefore conclude that a more real-
istic 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.