© The Authors 2024

1. Introduction

The large-scale distribution of galaxies is an extremely important source of cosmological information from the low-redshift Universe, complementing observations of the cosmic microwave background (CMB), such as those made using the Wilkinson Microwave Anisotropy Probe (WMAP; Hinshaw et al. 2013) and Planck (Planck Collaboration VI 2020). In the course of the past two decades, observations of the large-scale structure (LSS) from spectroscopic galaxy redshift surveys, such as the 2dF Galaxy Redshift Survey (2dFGRS; Colless et al. 2001), the 6dF Galaxy Survey (6dFGS; Jones et al. 2009), the VIMOS VLT Deep Survey (VVDS; Le Fèvre et al. 2013), the Sloan Digital Sky Survey (SDSS; York et al. 2000), the WiggleZ Dark Energy Survey (WiggleZ; Drinkwater et al. 2010), the VIMOS Public Extragalactic Redshift Survey (VIPERS; Guzzo et al. 2014), the Galaxy And Mass Assembly (GAMA; Driver et al. 2011), and the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2012) and its extension (eBOSS; Dawson et al. 2016), have provided a wealth of information on how gravitational instability shapes the large-scale matter distribution and on the relation between matter and galaxy density perturbations. At the same time, such observations have stood as a testing ground for what has ultimately emerged as the standard cosmological model.

In the next decade, this picture is going to be significantly enriched by analyses performed by Stage-IV spectroscopic surveys, such as the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016) and Euclid (Laureijs et al. 2011), which are going to explore a still relatively uncharted epoch at 1 ≲ z ≲ 2, when the Universe was only about half of its current age. In particular, Euclid is going to collect the redshift of millions of Hα-emitting galaxies across a total sky surface of 15 000deg², therefore increasing the statistical constraining power on the cosmological parameters to an unprecedented level for spectroscopic analyses in the low-redshift Universe. It comes with no surprise that the increase in statistical significance of the observations must necessarily be accompanied by an equivalent increase in the accuracy of the theoretical recipes used to analyse the data in order to keep systematic errors in the theory at a fraction of the statistical error budget. This becomes even more relevant in terms of the range of validity of the considered models, whose reach must be properly benchmarked against realistic mock samples.

The standard cosmological probe for galaxy clustering is the galaxy two-point correlation function (2PCF), or its Fourier transform, the galaxy power spectrum. Both statistics quantify the excess probability of finding galaxy pairs at a given separation with respect to the case of a purely random (Poissonian) distribution. These observables have been extensively used by recent experiments to place constraints on the cosmological parameters, either focusing on specific features such as baryon acoustic oscillation (BAO) and redshift-space distortions (RSDs) in the so-called template-fitting approach (Peacock et al. 2001; Tegmark et al. 2006; Guzzo et al. 2008; Blake et al. 2011; Reid et al. 2012; Beutler et al. 2012, 2016; Contreras et al. 2013; Howlett et al. 2015; Okumura et al. 2016; Alam et al. 2017; Pezzotta et al. 2017; Gil-Marín et al. 2018; Hou et al. 2018; Wang et al. 2018; Zhao et al. 2018), or in full-shape analyses (Sánchez et al. 2014, 2016; Grieb et al. 2017; Ivanov et al. 2020; Tröster et al. 2020; d’Amico et al. 2020; Semenaite et al. 2022, 2023; Chen et al. 2022; Philcox & Ivanov 2022; Carrilho et al. 2023; Moretti et al. 2023).

For both kinds of approaches, the inference of cosmological information is made intrinsically more difficult by the presence of three separate effects that build on linear theory predictions. These are the non-linear gravitational evolution of the dark matter density field (Bernardeau et al. 2002; Baumann et al. 2012; Carrasco et al. 2012), the relationship between the galaxy δ_g and the matter δ density fields, known as thr ‘galaxy bias’ (Kaiser 1984; Bardeen et al. 1986; see Desjacques et al. 2018 for a recent review), and finally the apparently anisotropic pattern in the distribution of galaxies due to the effect of their peculiar velocities on the observed redshift (RSDs, Kaiser 1987; Hamilton 1992; Fisher 1995; Scoccimarro et al. 1999b; Scoccimarro 2004; Taruya et al. 2010; Senatore & Zaldarriaga 2014; Perko et al. 2016). Each of these effects needs to be carefully modelled in order to infer accurate cosmological constraints from the full shape of the galaxy power spectrum/2PCF. This goal can be achieved in different ways: using numerical methods, such as N-body simulations (e.g. Kuhlen et al. 2012; Schneider et al. 2016; Springel et al. 2021), adopting analytical approaches based on a perturbative solution to the equations governing the evolution of the matter and galaxy density fields (e.g. Fry & Gaztanaga 1993; Bernardeau et al. 2002; McDonald & Roy 2009; Carrasco et al. 2012; Assassi et al. 2014; Senatore 2015; Desjacques et al. 2018), or resorting to hybrid methods that combine the previous two methodologies (e.g. Euclid Collaboration 2021; Angulo et al. 2021; Zennaro et al. 2023; Pellejero-Ibañez et al. 2023).

In terms of galaxy bias, it has become standard practice to adopt a perturbative expansion of the galaxy density field using Eulerian coordinates, which can be expressed as a sum of partial derivatives of the gravitational and velocity divergence potentials, each one weighted by a corresponding free parameter to be fitted against measurements. When considering the one-loop galaxy power spectrum, the list includes the linear bias, b₁, expressed in terms of the dark matter density field in the large-scale limit, as δ_g = b₁δ (Kaiser 1984; Bardeen et al. 1986; Cole & Kaiser 1989; Nusser & Davis 1994; Mo & White 1996; Sheth & Tormen 1999), plus the next-to-leading-order correction obtained from a spherically symmetric gravitational collapse via a power-law Taylor expansion δ_g = ∑_nb_n δⁿ/n! of the matter density field (Szalay 1988; Coles 1993; Fry & Gaztanaga 1993; Scoccimarro et al. 2001; Smith et al. 2007; Manera et al. 2010; Desjacques et al. 2010; Frusciante & Sheth 2012; Schmidt et al. 2013). This is further supplemented by the presence of non-local contributions generated by the cosmic tidal field (Bouchet et al. 1992; Catelan et al. 1998, 2000; McDonald & Roy 2009), which have been proven to be essential for a correct description of the clustering of dark matter halos (Manera & Gaztañaga 2011; Roth & Porciani 2011) and to secure consistency between the results from the analysis of two- and three-point correlation measurements (Pollack et al. 2012, 2014). These extra corrections were first detected in N-body simulations in Chan et al. (2012) and Baldauf & Seljak (2012), and have since then become a standard ingredient in the bias expansion. Additionally, the latter also takes into account the effects of short-range non-localities during the processes of galaxy formation, which lead to the presence of higher-than-second-order derivatives of the gravitational potential. At leading order in the power spectrum, higher derivatives appear with a term scaling as ∇²δ (Bardeen et al. 1986; Matsubara 1999; Desjacques 2008; Desjacques et al. 2010). Finally, the dependence on short-wavelength modes is included via an additional stochastic field ε_g(x) (Dekel & Lahav 1999; Sheth & Lemson 1999; Taruya & Soda 1999; Matsubara 1999; Bonoli & Pen 2009; Hamaus et al. 2010; Schmidt 2016; Ginzburg et al. 2017)¹, which is responsible for a shot-noise contribution to the power spectrum. This correction deviates from the predictions of a purely Poissonian distribution, and at the same time can introduce a scale dependence due to the physical scale at which two objects can be mistaken for a single one, similarly to the exclusion effect for dark matter halos (Scherrer & Weinberg 1998; Sheth & Lemson 1999; Cooray & Sheth 2002; Smith et al. 2007; Baldauf et al. 2013, 2016).

The high dimensionality of the parameter space for the model described above can be reduced by employing a set of physically motivated relations expressing a few higher-order bias parameters in terms of lower-order ones. A typical assumption is the conserved evolution of tracers (coevolution), which, from a local-in-matter-density expansion at the moment of formation, leads to the well-known local Lagrangian relations (Chan et al. 2012; Baldauf & Seljak 2012; Eggemeier et al. 2019). The latter have been adopted in the literature as a fairly conservative trade-off between sampling the whole set of bias parameters and fixing some of the model degrees of freedom, most notably in the analysis of the BOSS DR12 data release (Sánchez et al. 2016; Grieb et al. 2017) to improve the statistical constraints on the cosmological parameters obtained from the anisotropic 2PCF and power spectrum.

The standard bias expansion has been the subject of several tests in the literature, together with a validation of the coevolution relations mentioned in the previous paragraph. As an example, Saito et al. (2014) checked the consistency between the bias parameters fitted from the halo power spectrum and bispectrum (the Fourier transform of the three-point function) using a sample of measurements in different mass bins and at different redshifts, revealing an agreement between the two sets of bias measurements up to k ∼ 0.1 h Mpc⁻¹. The use of an irreducible bias basis, and also properly including a higher-derivative correction, was tested in Angulo et al. (2015), who showed that with this approach it is possible to extend the validity of the one-loop galaxy bias expansion up to k ∼ 0.3 h Mpc⁻¹ even at z = 0. More recently, Eggemeier et al. (2020) analysed the accuracy of this expansion at fixed cosmology using simulated halo occupation distribution (HOD) catalogues built to mimic the clustering properties of the SDSS Main Galaxy Sample (Strauss et al. 2002) and of the BOSS CMASS and LOWZ samples (Eisenstein et al. 2011; Dawson et al. 2012; Reid et al. 2016). The authors focused on the necessity to introduce both a higher-derivative term and a scale-dependent correction to shot noise while analysing the auto galaxy and cross galaxy-matter power spectrum. Findings from this study indicate that the standard one-loop bias expansion can be broken on scales k ∼ 0.2 h Mpc⁻¹ unless higher-order stochastic corrections are taken into consideration. Pezzotta et al. (2021) and Eggemeier et al. (2021) extended this analysis to include a determination of the cosmological parameters, and explored the additional constraining power coming from the one-loop galaxy bispectrum. These works show how fixing the quadratic tidal bias as a function of the linear bias provides accurate results up to k ∼ 0.35 h Mpc⁻¹ for the galaxy power spectrum. A similar analysis was carried out by Oddo et al. (2021), who assessed the constraining power of the galaxy bispectrum on the cosmological parameters, displaying a consistency up to k ∼ 0.3 h Mpc⁻¹ for the one-loop power spectrum and k ∼ 0.09 h Mpc⁻¹ for the tree-level bispectrum. Equivalent analyses in redshift space (Markovič et al. 2019; Bose et al. 2020; de la Bella et al. 2020; Gualdi et al. 2021; Rizzo et al. 2023; Nicola et al. 2024), or in terms of field-level comparisons (Schmittfull et al. 2019), have also appeared in recent years, leading to compatible scenarios.

On a partially different side, numerical simulations (see e.g. Kuhlen et al. 2012, for a review) have proven to be an optimal way to reproduce the evolution of the matter and galaxy density field deep into the non-linear regime, and their use in analyses of galaxy redshift surveys has therefore multiplied in recent years thanks to a large number of different suites, such as DEMNUni (Castorina et al. 2015), UNIT (Chuang et al. 2019), Quijote (Villaescusa-Navarro et al. 2020), Uchuu (Ishiyama et al. 2021), and AbacusSummit (Maksimova et al. 2021; Yuan et al. 2022) simulations. In quantitative terms, state-of-art N-body simulations can achieve an accuracy of better than 2% on the shape of the non-linear matter power spectrum down to scales of k ∼ 10 h Mpc⁻¹ (Schneider et al. 2016; Springel et al. 2021; Angulo et al. 2021), but unfortunately their application as a tool to infer cosmological information from real observations is limited by their extreme computational cost. However, in recent years, different methods have been proposed with the goal of increasing their range of applicability. This ranges from methods meant to speed up their production (e.g. Monaco et al. 2002; Tassev et al. 2013; Izard et al. 2016) to ones designed to find an optimal interpolation strategy among a limited pool of high-resolution simulations (Heitmann et al. 2013; Liu et al. 2018; Nishimichi et al. 2019; DeRose et al. 2019; Giblin et al. 2019; Wibking et al. 2018; Winther et al. 2019; Euclid Collaboration 2019, 2021). Among this second category of methods, we highlight baccoemu (Angulo et al. 2021), an emulator for the non-linear matter power spectrum that was recently extended to also include biased tracers in real (Zennaro et al. 2023) and redshift space (Pellejero-Ibañez et al. 2023), assuming a hybrid Lagrangian bias model, with the individual terms of the expansion directly emulated from high-resolution simulations, and a cosmology-rescaling technique meant to reduce the number of simulations required to train the emulator.

In this paper, we compare different models and test different scale cuts and bias relations on a sample of synthetic galaxy catalogues tailored to reproduce – to the best of our knowledge – the clustering signal of the Hα sample that will be targeted by Euclid. As we are interested in the relative performance of different theory models, we do not consider the presence of observational systematic uncertainties in this analysis. For example, effects such as purity and completeness of the sample will induce variations in the comoving number density considered in this work. At the same time, the sample purity, which is determined by the presence of line and noise interlopers, will also modify the overall shape of the measured n-point statistics. All of these effects are going to be investigated by a dedicated group in the Euclid Consortium, while a specific analysis on theory model selection with a more realistic analysis (including survey mask, selection effects, and combining multiple redshift bins) is going to be developed in a future paper (Euclid Collaboration: Moretti et al., in prep.).

Our goal is to test the range of validity of the one-loop galaxy bias expansion, which we quantify by means of three different performance metrics (Osato et al. 2019): the figure of bias, quantifying the accuracy of the model in terms of the recovery of the model parameters; the goodness of fit, measuring how well the best-fit model compares to the input data vectors; and the figure of merit, quantifying the statistical power of the model in constraining the cosmological parameters. These metrics are computed for each of the configurations we test, as a function of the maximum wave mode k_max included in the fit, exploring different bias relations meant to reduce the dimensionality of the parameter space. As we limit our attention to the real-space galaxy power spectrum alone, we focus on the recovery of the dimensionless Hubble parameter h, which is defined in terms of the Hubble constant H₀ as H₀ = 100 h km s⁻¹ Mpc⁻¹, and of the cold dark matter density parameter ω_c ≡ Ω_ch², where Ω_c is the corresponding fractional density parameter. At the same time, we avoid sampling the scalar amplitude A_s, as this would lead to a strong degeneracy with the linear bias parameter b₁. This degeneracy can be partially broken if only considering the additional constraining power from higher-order statistics, cross-correlation with the matter density field, or in a multi-tracer analysis, or by considering the apparently anisotropic clustering amplitude when also including RSDs.

This work is the first installment in a series of Euclid preparation papers meant to validate the theoretical framework used to analyse the full shape of two- and three-point clustering measurements from the final data sample. Here, we focus on the real-space galaxy power spectrum, while the corresponding three-point equivalent for the real-space galaxy bispectrum is going to be presented in Euclid Collaboration: Eggemeier et al. (in prep.). Both of these papers describe tests conducted in real space, that is, using the true comoving positions of galaxies inside the box instead of the positions displaced because of RSDs. While this choice excludes one of the main observational probes of galaxy clustering, such analyses can provide an important testing ground for the model of galaxy bias. This includes calibration of optimal scale cuts for the model², as well as testing different ways to reduce the dimensionality of the parameter space, such as using the coevolution relations. In addition, real-space analyses can become relevant in the context of modelling 3 × 2 point statistics (photometric galaxy clustering, weak gravitational lensing, and galaxy–galaxy lensing), such as in the analysis performed by the Dark Energy Survey (see e.g. Pandey et al. 2022; Porredon et al. 2022, for cosmological inference that requires a proper modelling of photometric galaxy bias). On the other hand, future installments of this series will focus on the modelling of the redshift-space equivalents of the statistics adopted in these works.

This article is structured in the following way. In Sect. 2 we present the simulated galaxy samples and the power-spectrum measurements and covariances used throughout the paper. In Sect. 3 we describe the theoretical models that we employed for the analysis, while in Sect. 4 we describe the fitting procedure and the performance metrics used to quantify the goodness of the models for different configurations as a function of the maximum mode included in the fit. Finally, in Sect. 5 we present the results of the analysis, and we draw our conclusions in Sect. 6.

2. Data

2.1. Euclid simulations

In order to determine the performance of the selected theoretical models, we first need a set of simulated data samples spanning the same redshift range that will be observed by Euclid, and for which the input cosmology is known a priori.

In the following we make use of four comoving outputs, selected to cover the redshift range 0.9 < z < 1.8 of the Flagship I simulation³. The latter has been carried out on the supercomputer Piz Daint, which is hosted by the Swiss National Supercomputing Center (CSCS), using the PKDGRAV3 algorithm (Potter et al. 2017), and consists of a record-setting N-body run with two trillions dark matter particles moving under the effect of gravity within a box of size L = 3780 h⁻¹ Mpc. The mass resolution of the simulation (m_p ∼ 2.398 × 10⁹ h⁻¹ M_⊙ ) allows us to marginally resolve halos with a typical mass M_h of few 10¹⁰ h⁻¹ M_⊙, which host the majority of the Hα emission line galaxies that are going to be targeted by Euclid. The nominal flat ΛCDM cosmology adopted to run the simulation as stated in Potter et al. (2017) differs from the fiducial cosmology assumed in this paper in the value of the spectral index (n_s = 0.96 vs. n_s = 0.97). This choice has been motivated since, during the course of this study, we observed subtle yet significant inconsistencies between our models and the measurement in the Flagship I simulation. After contacting the team responsible for running the simulation, they confirmed that the nominal parameters of the simulation were previously wrongly communicated, and that they are in agreement with the ones we have identified. The latter are obtained by performing cosmological fits to high-resolution dark matter power spectrum measurements at various redshifts, and are summarized in Table 1. This procedure is detailed in Appendix A.

Each comoving snapshot has been populated with galaxies by firstly generating a friends-of-friends (FoF) halo catalogue with a linking length b = 0.2 and a minimum halo mass corresponding to ten dark matter particles⁴, where the halo mass M_h is defined as the sum of all the identified particles. Subsequently, halos have been populated with galaxies using a HOD algorithm to match the abundance and clustering of the Hα samples implemented in the main Flagship I lightcone catalog⁵. The latter, in turn, is meant to reproduce the number density and clustering properties corresponding to two different Hα profiles, labelled as Model 1 and Model 3 in Pozzetti et al. (2016). These samples are defined by different templates for the evolution of the luminosity function ϕ(L, z), from the use of a standard Schechter parametrization for Model 1 (Schechter 1976), to the direct fit to real observations for Model 3. The net effect in terms of number density is that the Model 1 sample has almost twice as many objects as Model 3, which is more conservative in the selection of Hα emitters, as shown in Fig. 4 of Pozzetti et al. (2016)⁶. A more detailed description of the Flagship Hα lightcone and the pipeline for its construction will be provided in Euclid Collaboration: Castander et al., (in prep).

In terms of the comoving snapshots, the HOD we implemented consists of a 8-parameter model, where the mean occupation numbers of central and satellite galaxies are defined as

$$\begin{array}{cc}\hfill \langle N_{\mathrm{cen}}\rangle (M_{\mathrm{h}})& =\frac{1}{2}{f}_{\mathrm{cen}}^{\max }[1+\mathrm{erf}\left(\frac{log{M}_{\mathrm{h}}-log{M}_{\min }}{{\sigma }_{logM}}\right)]\hfill \\ \hfill & \times [1-\frac{1-f_{\mathrm{cen}}^{\min }/f_{\mathrm{cen}}^{\max }}{1+{10}^{\frac{2}{k}(log{M}_{\mathrm{h}}-log{M}_{\mathrm{drop}})}}],\hfill \end{array}$$(1)

$$\begin{array}{c}\hfill \langle N_{\mathrm{sat}}\rangle (M_{\mathrm{h}})=N_{\mathrm{cen}}{\left(\frac{M_{\mathrm{h}}}{M_1}\right)}^{\alpha }.\end{array}$$(2)

Here, M_min is the typical minimum mass of halos hosting a central galaxy, σ_{log M} is the dispersion around M_min, and $f_{\mathrm{cen}}^{\max }$ is the amplitude of the central galaxy occupation. We include mass-dependency of ⟨N_cen⟩ above the transition scale M_min using three additional parameters, M_drop, k, and $f_{\mathrm{cen}}^{\min }$. Finally, the mean occupancy of satellite galaxies is regulated by M₁, which is a simple normalisation factor, and α, which corresponds to the slope of the power law distribution of satellites. In order to determine the distribution of galaxies inside halos, we employ a standard NFW profile (Navarro et al. 1996).

In Fig. 1 we show the fitted HOD profiles as a function of the host halo mass M_h for the eight samples: four redshifts times two different models. In all panels, both centrals and satellites profiles are shown, marked respectively with continuous and dashed lines. The mean occupation of central galaxies $f_{\mathrm{cen}}^{\max }$ inside dark matter halos does not converge to one, even for the most massive halos selected by the halo finder. This is a consequence of having selected a subsample (Hα in this case) from the whole population present in the ligthcone. The typical expection value for the occupancy of central galaxies for M_h > M_min slightly varies with the different samples, but is typically close to 0.2. This follows from the property of Hα emitters to be relatively young, blue, and star-forming galaxies, whereas, in massive halos, environmental effects such as galaxy-galaxy interactions, ram pressure stripping, and AGN feedback can suppress star formation in galaxies, effectively reducing the likelihood of finding actively star-forming central galaxies.

The total number of galaxies for each sample, their number density, and the scale k_sn at which the Poissonian shot-noise $P_{\mathrm{sn}}\equiv 1/\overline{n}$ becomes the dominant contribution in the data vectors, are listed in Table 2. A warning to be made is that, to determine the parameters for the HOD models, we selected galaxies from the lightcone assuming a significantly faint Hα flux limit, corresponding to f_Hα = 2 × 10⁻¹⁶ erg cm⁻² s⁻¹ (Scaramella et al. 2022), without assuming any realistic observational effect, such as target incompleteness, purity of the sample, and the impact of the angular footprint and radial selection function (Euclid Collaboration: Granett et al., in prep., Euclid Collaboration: Monaco et al., in prep.). This results in a sample with higher number density (see e.g. Bagley et al. 2020, for a recent forecast of Hα emitters from HST), with measured galaxy power spectra that are less affected by shot-noise. At the same time, the lack of line and noise interlopers allows us to neglect any extra contribution (Euclid Collaboration: Risso et al., in prep., Euclid Collaboration: Lee et al., in prep.) to the model galaxy power spectrum presented in Sect. 3. As a consequence, given the high precision assumed to validate the theory models, we believe that our tests should provide a conservative estimate of their range of validity. We leave to future Euclid analyses a more dedicated study of the impact of observational systematics.

2.2. Measurements and covariances

For each of the samples described above, we measured the real-space galaxy power spectrum P_gg(k) using the publicly available PowerI4 code⁷. The latter provides the functionality to compute the power spectrum from a particle distribution within a regular cubic box, using a variety of particle assignment schemes. For this analysis, we made use of a fourth-order interpolation scheme, otherwise known as piecewise cubic spline (PCS; see Sefusatti et al. 2016, for the exact form of the kernel), coupled with an interlacing method to reduce the aliasing contribution at high wave modes k.

We measured the power spectrum in the range defined by [k_F,k_Nyq], where k_F = 2π/L ∼ 0.0017 h Mpc⁻¹ is the fundamental frequency in a box of linear size L , and $k_{\mathrm{Nyq}}=\pi N_{\mathrm{grid}}/L$ is the Nyquist frequency corresponding to a density grid of linear size N_grid. We choose a grid resolution of N_grid = 1024 for the three dimensions of the box, to obtain measurements of the power spectrum up to a maximum wave mode of 0.8 h Mpc⁻¹, and we sampled the k range using a linear binning with step Δk = k_F⁸. The top and central panels of Fig. 2 show the power spectrum measurements for the Model 1 and 3 HOD samples respectively, with the redshift evolution over the available simulation snapshots marked by different lines in each panel. Differently from the evolution of the matter power spectrum, the galaxy power spectrum features an increasingly lower amplitude at lower redshifts. This can be explained by a larger linear galaxy bias at high redshift that overcomes the growth of matter fluctuations.

Since only a single realization is available for each of the HOD samples, we estimate the error covariance matrices associated to the data vectors using an analytical prediction in the Gaussian approximation, as explained in Grieb et al. (2016). This implies that the variance ${\sigma }_P^2(k)$ associated to each k bin is independent from the value of the galaxy power spectrum at different modes, and can be written as

$$\begin{array}{c}\hfill {\sigma }_P^2(k)=\frac{2}{N_k}{P}_{\mathrm{gg}}^2(k),\end{array}$$(3)

where N_k identifies the number of independent wave modes falling in the bin [k−Δk/2,k+Δk/2], while P_gg(k) is the theoretical non-linear galaxy power spectrum including shot-noise contributions. The latter has been obtained from a preliminary fit of the full non-linear model to the data vector of each sample, assuming, in a first iteration, an approximate but reasonable evaluation of the covariance itself. We expect the Gaussian approximation to be sufficient to our goals. The only additional contribution due to the galaxy trispectrum, here neglected, while noticeable at the relevant scales (Scoccimarro et al. 1999a; Sefusatti et al. 2006; Blot et al. 2015, 2016; Bertolini et al. 2016; Wadekar & Scoccimarro 2020) does not lead to significant differences (≲10%) on parameters constraints in the mildly non-linear regime (Blot et al. 2019; Wadekar et al. 2020). This is also supported by the goal of this analysis, which is testing the relative performance of different theory models rather than providing absolute values for the parameter uncertainties. However, we highlight how the Gaussian approximation is not expected to deliver completely realistic error bars (to the level of accuracy mentioned above), and that future analyses, both on simulated and real data will be integrated with a more complex model also including non-linear corrections (either analytical or from N-body simulations).

The bottom panel of Fig. 2 shows the standard deviation normalised by the corresponding galaxy power spectrum, where the latter is shot-noise-subtracted to highlight the different level of noise in our samples. For this reason, the eight cases exhibit a deviation from the linear relation in Eq. (3) (shown with a black dashed line) at small (k ≳ 0.1 h Mpc⁻¹) scales, where the shot-noise correction starts to become dominant with respect to the power spectrum signal. It should be noted how the relative error over this range of scales is well below the 1% level. Similarly, at large scales (k ≲ 0.005 h Mpc⁻¹) we find a departure from the linear relationship due to the small amplitude of P_gg, which can be clearly observed from the top and middle panel of Fig. 2. Finally, we note that this range of scales is also partially dominated by cosmic variance, due to the use of a single realisation of the Flagship simulation.

2.3. Volume rescaling

The main goal of our analysis is to carry out stringent tests to determine the range of validity of the standard one-loop galaxy bias model on the redshift range that will be explored by Euclid. For this we make use of a volume V_box corresponding to the full-box size of the Flagship comoving outputs, which is significantly larger than the volume that will be covered by Euclid. At the same time, we are interested in assessing the constraining power of the real-space galaxy power spectrum using a reference volume close to the one of an expected redshift bin of the full Euclid volume, V_shell. With this purpose in mind, we define new covariance matrices for the different samples presented in the previous sections, with an overall amplitude rescaled by the ratio between the volume of the Flagship comoving outputs and that of the Euclid-like shells,

$$\begin{array}{c}\hfill \eta =\frac{V_{\mathrm{box}}}{V_{\mathrm{shell}}},\end{array}$$(4)

such that the rescaled covariance C_shell can be expressed⁹ in terms of the original full-box covariance C_box as

$$\begin{array}{c}\hfill C_{\mathrm{shell}}=\eta C_{\mathrm{box}}.\end{array}$$(5)

We follow Euclid Collaboration (2020) and assume four non-overlapping redshift shells, centered at z = (0.9, 1.2, 1.5, 1.8), and with a depth of Δz = (0.2, 0.2, 0.2, 0.3), respectively, over a total projected area of 15 000 square degrees. With these values, we derive volume factors η for each of the considered redshift bins, shown in the last column of Table 2. We note that the mean values of the four redshift shells used in Euclid Collaboration (2020) do not match perfectly the redshifts of the four comoving snapshots used in this work. However, this is only marginally relevant, since we do not carry out a proper comparison to the Fisher forecasts obtained in that analysis. In fact, this will be a more suited aspect of investigation when considering the same observables, that is, the Legendre multipoles of the anisotropic galaxy power spectrum, and especially when considering more realistic number densities, as pointed out in Sect. 2.1.

A proper comparison between the results obtained using the full-box volume and the rescaled ones is presented in Sect. 5.5. In addition to the Euclid-like shells, we consider three additional volume rescalings, by dividing the range between V_box and V_shell into four evenly sized intervals. This leads to a total of five different sets of covariances, based on the volumes defined above.

3. Theoretical model

In this section we describe the theoretical framework of perturbation theory (PT), which is essential to understand the evolution of post-inflationary fluctuations in the matter density field δ into the current large-scale distribution of galaxies via gravitational instability. This description is expected to be accurate only down to the mildly non-linear regime, where the amplitude of the density contrast δ is small enough to be perturbatively expanded. In the strong non-linear regime we expect this model to fail, as gravitational collapse leads to the formation of bound structures beyond the regime of validity of perturbative approaches.

For convenience, in the rest of this article we use the following notation for the integration over the infinite volume of a loop variable q,

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\mathit{q}}}\equiv {\displaystyle \int \frac{{\mathrm{d}}^{3}q}{{(2\pi )}^3}},\end{array}$$(6)

and adopt the following convention for the direct and inverse Fourier transform of the density contrast,

$$\begin{array}{cc}& \delta (\mathit{k})\equiv {(2\pi )}^3{\displaystyle {\int }_{\mathit{x}}}{\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{x}}\delta (\mathit{x}),\hfill \end{array}$$(7)

$$\begin{array}{cc}& \delta (\mathit{x})\equiv {\displaystyle {\int }_{\mathit{k}}}{\mathrm{e}}^{\mathrm{i}\mathit{k}·\mathit{x}}\delta (\mathit{k}).\hfill \end{array}$$(8)

The three-dimensional Dirac function is represented with the standard notation ${\delta }_{\mathrm{D}}^{(3)}$. Finally, the power spectrum P_XX(k, z) of any component, matter or biased tracer, is defined as the auto-correlation of the corresponding density field δ_X, such that

$$\begin{array}{c}\hfill \langle {\delta }_{\mathrm{X}}(\mathit{k}){\delta }_{\mathrm{X}}({\mathit{k}}^{\prime })\rangle \equiv {(2\pi )}^3{P}_{\mathrm{XX}}(k){\delta }_{\mathrm{D}}^{(3)}(\mathit{k}+{\mathit{k}}^{\prime }),\end{array}$$(9)

where the presence of the Dirac function and the independence of the power spectrum from the orientation of the wave mode k reflect the underlying assumption of homogeneity and isotropy.

3.1. Eulerian framework and effective field theory

3.1.1. Modelling of the non-linear matter power spectrum

We begin by summarising the most relevant outcomes of standard perturbation theory (SPT; see e.g. Bernardeau et al. 2002, for a review on the subject). Its main assumption is that on large scales the dynamics of dark matter can be approximated as that of a perfectly pressureless fluid, with negligible effects from particle shell-crossing in multi-streaming regions. Under the so-called Einstein–de Sitter (EdS) approximation, we can write the matter density contrast using a perturbative expansion,

$$\begin{array}{c}\hfill \delta ={\delta }^{(1)}+{\delta }^{(2)}+{\delta }^{(3)}+\dots ,\end{array}$$(10)

where at each order n the individual contribution δ⁽ⁿ⁾ is a function of the linear density contrast δ_L¹⁰,

$$\begin{array}{cc}& {\delta }^{(n)}(\mathit{k})={\displaystyle {\int }_{{\mathit{q}}_1}}{\delta }_{\mathrm{D}}^{(3)}(\mathit{k}-{\mathit{q}}_{1\dots n})F_n({\mathit{q}}_1,\dots ,{\mathit{q}}_n)\hfill \\ \hfill & \times {\delta }_{\mathrm{L}}({\mathit{q}}_1)\dots {\delta }_{\mathrm{L}}({\mathit{q}}_n).\hfill \end{array}$$(11)

Here q_1… n ≡ q₁ + ⋯ + q_n, and F_n is the nth order symmetrised PT kernel describing the non-linear interaction among fluctuations at different wave modes q₁,  …, q_n. The vanishing argument of ${\delta }_{\mathrm{D}}^{(3)}$ reflects the translational invariance of the equations of motion in a spatially homogeneous universe.

Similarly, the non-linear matter power spectrum P_mm(k) can be expanded by combining Eqs. (9) and (10), leading to

$$\begin{array}{c}\hfill P_{\mathrm{mm}}(k)=P_{\mathrm{L}}(k)+P^{1\text{-}\mathrm{loop}}(k)+P^{2\text{-}\mathrm{loop}}(k)+\dots ,\end{array}$$(12)

where $P_{\mathrm{L}}\sim \langle {\delta }_{\mathrm{L}}^2\rangle$ corresponds to the linear matter power spectrum, and at one-loop the only non-vanishing contributions are

$$\begin{array}{cc}& P^{1\text{-}\mathrm{loop}}(k)=P_{22}(k)+P_{13}(k)\hfill \\ \hfill & =2{\displaystyle {\int }_{\mathit{q}}}{F}_2^2(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}\left(|\mathit{k}-\mathit{q}|\right)P_{\mathrm{L}}(q)\hfill \\ \hfill & +6P_{\mathrm{L}}(k){\displaystyle {\int }_{\mathit{q}}}{F}_3(\mathit{q},-\mathit{q},\mathit{k})P_{\mathrm{L}}(q).\hfill \end{array}$$(13)

For the sake of completeness, we report the expanded expressions for the second- and third-order symmetrised kernels, F₂(q₁, q₂) and F₃(q₁, q₂, q₃) in Appendix B.

The one-loop model in SPT, however, fails to accurately describe the non-linear damping of the acoustic oscillations due to bulk flow displacements (Eisenstein et al. 2007; Crocce & Scoccimarro 2008; Baldauf et al. 2015b). At first order, this effect can be reproduced in the theoretical model for P_mm(k) by a proper resummation of all infrared (IR) modes q < k, that is, of comoving separations larger than the one under consideration (see Crocce & Scoccimarro 2006, 2008 for a description of the BAO smearing in the context of renormalised perturbation theory).

A more standard procedure to include these corrections is based on the split of the linear power spectrum P_L as the sum of a smooth P_nw and wiggly P_w component (Seo et al. 2008; Baldauf et al. 2015b; Blas et al. 2016), that is

$$\begin{array}{c}\hfill P_{\mathrm{L}}(k)=P_{\mathrm{nw}}(k)+P_{\mathrm{w}}(k).\end{array}$$(14)

At leading order, it is possible to estimate the amplitude of the damping factor making use of the Zeldovich approximation (Zel’dovich 1970). This leads to an expression for the leading-order, IR-resummed power spectrum,

$$\begin{array}{c}\hfill P_{\mathrm{mm}}^{\mathrm{IR}\text{-}\mathrm{LO}}(k)=P_{\mathrm{nw}}(k)+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}(k),\end{array}$$(15)

where Σ², representing the variance of the relative displacement field (Eisenstein et al. 2007), is defined as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}^2=\frac{1}{6{\pi }^2}{\displaystyle {\int }_0^{k_{\mathrm{s}}}}{P}_{\mathrm{nw}}(q)[1-j_0\left(\frac{q}{k_{\mathrm{osc}}}\right)+2j_2\left(\frac{q}{k_{\mathrm{osc}}}\right)]\mathrm{d}q.\end{array}$$(16)

Here j_n is the nth order spherical Bessel function of the first kind, k_osc = 1/ℓ_osc is the wavelength corresponding to the BAO scale ℓ_osc = 110  h⁻¹ Mpc¹¹, and k_s = 0.2 h Mpc⁻¹ is the ultraviolet (UV) integration limit¹².

The next-to-leading order correction can be written by using the leading order term of Eq. (15) inside the expression for the one-loop corrections of Eq. (13). This leads to the final formulation for the non-linear IR-resummed power spectrum (Baldauf et al. 2015b; Blas et al. 2016),

$$\begin{array}{cc}\hfill P_{\mathrm{mm}}^{\mathrm{IR}\text{-}(\mathrm{LO}+\mathrm{NLO})}(k)& =P_{\mathrm{nw}}(k)+(1+k^2{\mathrm{\Sigma }}^2){\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}(k)\hfill \\ \hfill & +P^{1\text{-}\mathrm{loop}}[P_{\mathrm{nw}}+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}](k),\hfill \end{array}$$(17)

where the square brackets of the last term mean that the evaluation of the one-loop correction is carried out using the leading order IR-resummed power spectrum in place of the linear one.

Another partial failure of the model, which is equally shared by any recipe based on perturbative methods, is that its range of validity is limited to quasi-linear scales, where the assumption of a pressureless fluid is still justified. However, on scales approaching the non-linear scale k_NL¹³, dark matter particles experience shell-crossing, effectively introducing a non-zero pressure, which under more realistic conditions is further enhanced by the presence of baryonic processes, such as galaxy formation, ISM cooling, and AGN and supernovae feedback. These effects can be described in terms of a non-trivial stress-energy tensor which, at leading order, results in an additional contribution to the matter power spectrum (Pueblas & Scoccimarro 2009; Carrasco et al. 2012; Baumann et al. 2012),

$$\begin{array}{c}\hfill P_{\mathrm{ctr}}(k)=-2c_{\mathrm{s}}^2{k}^2{P}_{\mathrm{mm}}^{\mathrm{IR}\text{-}\mathrm{LO}}(k),\end{array}$$(18)

usually denoted as counterterm in the EFTofLSS framework. Here, the parameter c_s can be interpreted as an effective speed of sound (Baumann et al. 2012; Carrasco et al. 2014; Baldauf et al. 2015a), reflecting the influence of short-wavelength perturbations, but accounts as well for the complex physics behind galaxy formation (when considering biased tracers of the matter density field).

Summarising, we can write the final expression for the model of the non-linear matter power spectrum as

$$\begin{array}{c}\hfill P_{\mathrm{mm}}(k)=P_{\mathrm{mm}}^{\mathrm{IR}\text{-}(\mathrm{LO}+\mathrm{NLO})}(k)+P_{\mathrm{ctr}}(k),\end{array}$$(19)

which contains one free parameter, c_s, which must be treated as a nuisance parameter to be fitted against real or, in our case, simulated measurements.

3.1.2. Modelling of the non-linear galaxy power spectrum

The general perturbative expansion of the galaxy density field δ_g is based on the sum of all the individual operators that are a function of properties of the environment in which galaxies reside, such as the underlying matter density field and the large-scale tidal field. More precisely, this sum includes all those operators that are sourced by second derivatives of the gravitational potential Φ and the velocity potential Φ_v (see Desjacques et al. 2018, for a detailed review on the subject).

If we restrict our model to the one-loop prediction for the power spectrum, the relation between δ_g and δ can be described considering only terms up to third order in the perturbations. In detail, this relation can be written as

$$\begin{array}{cc}& {\delta }_{\mathrm{g}}(\mathit{x})=b_1\delta (\mathit{x})+b_{{\mathrm{\nabla }}^2\delta }{\mathrm{\nabla }}^2\delta (\mathit{x})+{\epsilon }_{\mathrm{g}}(\mathit{x})\hfill \\ \hfill & +\frac{b_2}{2}{\delta }^2(\mathit{x})+b_{{\mathcal{G}}_2}{\mathcal{G}}_2\left({\mathrm{\Phi }}_v|\mathit{x}\right)+b_{{\mathrm{\Gamma }}_3}{\mathrm{\Gamma }}_3(\mathit{x})+\dots ,\hfill \end{array}$$(20)

where each operator is multiplied by a free bias parameter that determines its overall amplitude¹⁴. The different terms in Eq. (20) can be summarised as follows.

• (i)

  At leading order, the shot-noise-corrected galaxy density field can be expressed using a linear and local relation in δ. This relation is characterised by a linear bias parameter, b₁, which simply rescales the underlying matter density contrast by a constant factor (Kaiser 1984).

• (ii)

  The effect of short-range non-localities during the process of galaxy formation is characterised by the presence of higher derivatives of the gravitational potential (Bardeen et al. 1986; McDonald & Roy 2009; Desjacques et al. 2010). At leading order, the only non-zero term scales with the Laplacian of the matter density field, ∇²δ, and has an amplitude regulated by the free parameter b_∇²δ. The formation of structures involves the collapse of matter from a finite region of space, which for dark matter halos is well approximated by their Lagrangian radius R. Since the estimation of the corresponding radius for a given galaxy sample can be cumbersome, here we absorb the value of R inside the definition of b_∇²δ.

• (iii)

  The impact of short-scale fluctuations on the galaxy density field at larger separations is determined by an additional stochastic field, ε_g, which, under the assumption of Gaussian initial conditions, is completely uncorrelated from large-scale perturbations. If galaxies are randomly distributed, the stochastic contribution to the galaxy power spectrum is purely represented by the Poisson limit, $1/\overline{n}$, with $\overline{n}$ corresponding to the mean number density of the selected sample.

• (iv)

  Moving to mildly non-linear scales, higher-order correlations of the density field appear (Coles 1993; Fry & Gaztanaga 1993), starting with a term proportional to δ², characterised by a quadratic local bias b₂. This factor is expected from a spherically symmetric gravitational collapse, in which higher powers of δ become more relevant at progressively smaller scales. The third power of the matter density field is not included in Eq. (20) since its effect on the one-loop galaxy power spectrum is an extra contribution to the large-scale limit, which can be absorbed in the renormalization of the linear bias.

• (v)

  Even starting with a purely local-in-matter-density bias expansion at the time of formation, non-linear evolution is responsible for the generation of large-scale tidal fields (Chan et al. 2012; Baldauf & Seljak 2012). At leading order, the corrections given by the tidal stress tensor are represented by a non-local quadratic bias, $b_{{\mathcal{G}}_2}$, and by the second-order Galileon operator, 𝒢₂, defined as

  $$\begin{array}{c}\hfill {\mathcal{G}}_2(\mathrm{\Phi }|\mathit{x})\equiv {\left[{\mathrm{\nabla }}_{\mathit{ij}}\mathrm{\Phi }(\mathit{x})\right]}^2-{\left[{\mathrm{\nabla }}^2\mathrm{\Phi }(\mathit{x})\right]}^2.\end{array}$$(21)

  In Fourier space, Eq. (21) can be written as

  $$\begin{array}{c}\hfill {\mathcal{G}}_2(\mathit{k})={\displaystyle {\int }_{\mathit{q}}}S(\mathit{q},\mathit{k}-\mathit{q})\delta (\mathit{q})\delta (\mathit{k}-\mathit{q}),\end{array}$$(22)

  where

  $$\begin{array}{c}\hfill S({\mathit{k}}_1,{\mathit{k}}_2)\equiv \frac{{({\mathit{k}}_1·{\mathit{k}}_2)}^2}{k_1^2{k}_2^2}-1\end{array}$$(23)

  is the Fourier-space kernel corresponding to the second-order Galileon operator 𝒢₂.

• (vi)

  The next-to-leading-order correction to the tidal field can be obtained considering terms up to second-order in the potential of the displacement field (Chan et al. 2012). This contribution is represented by an additional non-local cubic bias, b_Γ3, and by the operator

  $$\begin{array}{c}\hfill {\mathrm{\Gamma }}_3(\mathit{x})\equiv {\mathcal{G}}_2(\mathrm{\Phi }|\mathit{x})-{\mathcal{G}}_2\left({\mathrm{\Phi }}_v|\mathit{x}\right),\end{array}$$(24)

  whose net effect inside Eq. (20) is to include terms up to third order in perturbations of δ.

All the terms giving a non-zero contribution to the one-loop galaxy power spectrum are listed in Eq. (20). The complete expression for P_gg then reads

$$\begin{array}{c}\hfill P_{\mathrm{gg}}(k)=P_{\mathrm{gg}}^{\mathrm{tree}}(k)+P_{\mathrm{gg}}^{1\text{-}\mathrm{loop}}(k)+P_{\mathrm{gg}}^{\mathrm{ctr}}(k)+P_{\mathrm{gg}}^{\mathrm{noise}}(k),\end{array}$$(25)

where the individual contributions can be written as

$$\begin{array}{c}\hfill P_{\mathrm{gg}}^{\mathrm{tree}}(k)=b_1^2{P}_{\mathrm{L}}(k),\end{array}$$(26)

$$\begin{array}{cc}\hfill P_{\mathrm{gg}}^{1\text{-loop}}(k)& =P_{\mathrm{gg},22}(k)+P_{\mathrm{gg},13}(k)\hfill \\ \hfill & =2{\displaystyle {\int }_{\mathit{q}}}{K}_2^2(\mathit{q},\mathit{k}-\mathit{q})P_{\mathrm{L}}\left(|\mathit{k}-\mathit{q}|\right)P_{\mathrm{L}}(q)\hfill \\ \hfill & +6b_1{P}_{\mathrm{L}}(k){\displaystyle {\int }_{\mathit{q}}}{K}_3(\mathit{q},-\mathit{q},\mathit{k})P_{\mathrm{L}}(q),\hfill \end{array}$$(27)

$$\begin{array}{cc}& P_{\mathrm{gg}}^{\mathrm{ctr}}(k)=-2b_1(b_1{c}_{\mathrm{s}}^2+b_{{\mathrm{\nabla }}^2\delta })k^2{P}_{\mathrm{L}}(k)\hfill \\ \hfill & \equiv -2c_0{k}^2{P}_{\mathrm{L}}(k),\hfill \end{array}$$(28)

$$\begin{array}{c}\hfill P_{\mathrm{gg}}^{\mathrm{n}oise}(k)=\frac{1}{\overline{n}}(1+{\alpha }_{P,1}+{\alpha }_{P,2}{k}^2).\end{array}$$(29)

For the sake of completeness, a complete list of the individual one-loop corrections can be found in Appendix B. In the previous expressions, K₂ and K₃ are the symmetrised mode-coupling kernels for a generic biased tracer of the matter density field that follows the parametrization given in Eq. (20). In detail, they read

$$\begin{array}{c}\hfill K_2({\mathit{k}}_1,{\mathit{k}}_2)=b_1{F}_2({\mathit{k}}_1,{\mathit{k}}_2)+\frac{1}{2}{b}_2+b_{{\mathcal{G}}_2}S({\mathit{k}}_1,{\mathit{k}}_2),\end{array}$$(30)

and

$$\begin{array}{cc}& K_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)=b_1{F}_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)+b_2{F}_2({\mathit{k}}_1,{\mathit{k}}_2)\hfill \\ \hfill & +2b_{{\mathcal{G}}_2}S({\mathit{k}}_1,{\mathit{k}}_{23})F_2({\mathit{k}}_2,{\mathit{k}}_3)\hfill \\ \hfill & +2b_{{\mathrm{\Gamma }}_3}S({\mathit{k}}_1,{\mathit{k}}_{23})[F_2({\mathit{k}}_2,{\mathit{k}}_3)-G_2({\mathit{k}}_2,{\mathit{k}}_3)],\hfill \end{array}$$(31)

where G₂(k₁, k₂) is the standard one-loop kernel for the non-linear evolution of the velocity divergence field, and Eq. (31) has to be symmetrised with respect to its arguments (k₁,k₂,k₃).

Inside Eq. (29), α_P, 1 is a free nuisance parameter that accounts for deviations from a purely Poissonian shot-noise¹⁵. In addition, it is also required as a way to reabsorb the otherwise non-zero low-k limit of one of the individual one-loop contributions, as explained in Appendix B. Similarly, α_P, 2 parametrises the next-to-leading order correction, which scale as k².

Since the leading-order higher-derivative correction is completely degenerate with the matter counterterm, as they are both proportional to the combination k² P_L(k), we define a new more suited parameter,

$$\begin{array}{c}\hfill c_0\equiv b_1(b_1{c}_{\mathrm{s}}^2+b_{{\mathrm{\nabla }}^2\delta }),\end{array}$$(32)

to avoid the presence of unnecessary degeneracies between the parameters of the model.

In the previous scheme we have deliberately omitted the resummation of infrared modes, but, similarly to the case discussed in Sect. 3.1.1, galaxy two-point clustering also has to be corrected for the effect of large-scale bulk motions. For this reason, we write the relations for the leading- and next-to-leading order IR-resummed galaxy power spectra (mimicking Eqs. (15) and (17)) as

$$\begin{array}{c}\hfill P_{\mathrm{gg}}^{\mathrm{IR}-\mathrm{LO}}(k)=b_1^2[P_{\mathrm{nw}}(k)+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}(k)]+\frac{1}{\overline{n}}(1+{\alpha }_{P,1}),\end{array}$$(33)

$$\begin{array}{cc}\hfill P_{\mathrm{gg}}^{\mathrm{IR}-(\mathrm{LO}+\mathrm{NLO})}(k)& =b_1^2[P_{\mathrm{nw}}(k)+(1+k^2{\mathrm{\Sigma }}^2){\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}(k)]\hfill \\ \hfill & +P_{\mathrm{gg}}^{1\text{-loop}}[P_{\mathrm{nw}}+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}](k)\hfill \\ \hfill & +P_{\mathrm{gg}}^{\mathrm{ctr}}[P_{\mathrm{nw}}+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2}{P}_{\mathrm{w}}](k)+P_{\mathrm{gg}}^{\mathrm{n}oise}(k),\hfill \end{array}$$(34)

where, once again, the square brackets of the second and third terms in Eq. (34) reflect how the evaluation of the one-loop and counterterm contributions is carried out sourcing the leading order IR-resummed matter power spectrum, $P_{\mathrm{mm}}^{\mathrm{IR}\text{-}\mathrm{LO}}$, in place of the linear power spectrum, P_L(k).

3.1.3. Coevolution relations

A significant fraction of the bias parameters that have been introduced in this section enters in the expression for P_gg(k) only at higher-order, as clearly pointed out by the presence of only the linear bias b₁ in the expression for the leading-order galaxy power spectrum Eq. (33). This is significantly different from higher-order correlators of the galaxy density field, such as the galaxy bispectrum, for which both the local and non-local quadratic biases, b₂ and $b_{{\mathcal{G}}_2}$, appear also in the expression for the leading-order term, and can therefore be constrained with much better accuracy (Oddo et al. 2021; Eggemeier et al. 2021).

Given the poor constraining power of the galaxy power spectrum alone, it has become standard practice in real-data analyses to fix some of them to some physically motivated values or relations. This is important not only to obtain a larger constraining power for the remaining parameters, but also to ensure that none of them experiences strong degeneracies such as the one exhibited by the $(b_{{\mathcal{G}}_2},b_{{\mathrm{\Gamma }}_3})$ pair (see Appendix B). In this work, we test two different relations, which are briefly summarised in the next paragraphs.

As already explained in Sect. 3.1.2, even starting with a purely local-in-matter-density expression, δ_g(δ), at the time of formation, non-linear gravitational evolution is responsible for the generation of a large-scale tidal field (Fry 1996; Chan et al. 2012). This means that, even expressing the initial galaxy density field assuming only a spherically symmetric gravitational collapse – and thus with only local bias parameters b_n ≠ 0 – tidal contributions appear at later times because of gravitational evolution, leading to the presence of non-negligible tidal biases. Assuming that the total number of objects is conserved in time, it is possible to find a relation between the late-time non-local parameters and lower-order bias parameters, such that

$$\begin{array}{cc}& b_{{\mathcal{G}}_2}^{\mathrm{coev}}=-\frac{2}{7}(b_1-1)+b_{{\mathcal{G}}_2}^{\mathcal{L}},\hfill \end{array}$$(35)

$$\begin{array}{cc}& b_{{\mathrm{\Gamma }}_3}^{\mathrm{coev}}=-\frac{1}{6}(b_1-1)-\frac{5}{2}{b}_{{\mathcal{G}}_2}+b_{{\mathrm{\Gamma }}_3}^{\mathcal{L}},\hfill \end{array}$$(36)

where the bias parameters with a superscript ℒ stand for the corresponding Lagrangian quantities, that is, at the time of formation. The previous relations are commonly referred to as coevolution, or local Lagrangian relations when setting to zero the Lagrangian bias, and have been extensively used in most real-data analyses to fix one or both non-local parameters (Feldman et al. 2001; Gil-Marín et al. 2015; Sánchez et al. 2016; Grieb et al. 2017). However, recent results (Lazeyras & Schmidt 2018; Abidi & Baldauf 2018) have indicated that measurements from numerical simulations seem to suggest lower values for $b_{{\mathcal{G}}_2}$ with respect to its local Lagrangian relation.

An alternative approach for fixing $b_{{\mathcal{G}}_2}$, found to be more accurate when compared to results from N-body simulations and derived using the excursion-set formalism, has been proposed by Sheth et al. (2013). In this case, it is possible to express $b_{{\mathcal{G}}_2}$ as a quadratic form in terms of the linear bias b₁, such that

$$\begin{array}{c}\hfill b_{{\mathcal{G}}_2}^{\mathrm{ex}-\mathrm{set}}=0.524-0.547b_1+0.046b_1^2.\end{array}$$(37)

Such expressions are based on theoretical considerations on halo bias that only take into account the halo mass. As a consequence, they neglect potentially important effects, such as assembly bias (see e.g. Lazeyras et al. 2021; Barreira et al. 2021, for an analysis carried out, respectively, on dark matter halos and galaxies from hydrodinamical simulations). This means that their applicability to a real-data analysis must be carefully assessed (see Eggemeier et al. 2020; Pezzotta et al. 2021, for recent applications). Nonetheless, their use in this analysis is well justified, since we focus on HOD samples for which the assignment of a galaxy into a host halo is only determined by the mass of the latter.

In Sect. 5 we carry out tests to determine whether the previously defined relations can be employed to analyse clustering measurements adopting Euclid requirements.

3.2. Hybrid Lagrangian bias expansion model

In the previous sections, the relationship between the galaxy and the matter density field has been described through an Eulerian-based framework. However, this is not the only description of the galaxy power spectrum in the quasi-linear regime. Other approaches are possible, often based to various degrees on results from numerical simulations. We consider here the so-called ‘hybrid Lagrangian’ models. They draw from Lagrangian perturbation theory for the bias expression connecting galaxy and matter overdensities, but rely on simulations to capture the development of non-linearities when converting Lagrangian quantities to the observable Eulerian quantities.

The Lagrangian bias expansion describes the clustering of biased tracers in terms of a superposition of Lagrangian operators advected to Eulerian coordinates. It was first developed at one-loop in PT by Matsubara (2008), while Modi et al. (2020) proposed to combine the perturbative approach on bias with measurements of the advected operators from N-body simulations. This hybrid approach potentially allows us to push the bias expansion formalism to smaller scales with respect to purely perturbative approaches (see e.g. Hadzhiyska et al. 2021a; DeRose et al. 2023, for recent applications of the model to DES Y1 data and numerical simulations). Recent implementations include the works of Kokron et al. (2021) and Zennaro et al. (2023) in real space, and of Pellejero-Ibañez et al. (2023) in redshift space. In the present work we consider the implementation in the code baccoemu (Zennaro et al. 2023)¹⁶. It describes the Eulerian galaxy overdensity in terms of a second-order expansion of the Lagrangian galaxy density field δ_g(q) where q is the Lagrangian position corresponding to the Eulerian position x = q + Ψ(q) with Ψ(q) being the displacement field-connecting initial and final positions. This means that the Eulerian overdensity is given by

$$\begin{array}{c}\hfill 1+{\delta }_{\mathrm{g}}(\mathit{x})={\displaystyle {\int }_{\mathit{q}}}w(\mathit{q}){\delta }_{\mathrm{D}}^{(3)}(\mathit{x}-\mathit{q}-\mathbf{\Psi }(\mathit{q})),\end{array}$$(38)

where w(q) expresses the weighting function that transforms the matter field into the galaxy field,

$$\begin{array}{cc}& w(\mathit{q})=1+b_1^{\mathcal{L}}\delta (\mathit{q})+b_2^{\mathcal{L}}({\delta }^2(\mathit{q})-\langle {\delta }^2\rangle )\hfill \\ \hfill & +b_{s^2}^{\mathcal{L}}[s^2(\mathit{q})-\langle s^2\rangle ]+b_{{\mathrm{\nabla }}^2\delta }^{\mathcal{L}}{\mathrm{\nabla }}^2\delta (\mathit{q}).\hfill \end{array}$$(39)

Here the total list of operators built on the matter density field δ(q) consists of 𝒪 = {1, δ, δ², s², ∇²δ}, and the individual entries correspond to the fully non-linear matter distribution, not weighted (1) and weighted (δ) by the linear overdensity field, the squared linear overdensity field δ², the squared traceless tidal field s²¹⁷, and the Laplacian of the linear overdensity field ∇²δ, respectively. Note that unlike the Eulerian bias basis presented before, the expansion in Eq. (39) does not include the next-to-leading-order correction to the tidal field, captured by the operator Γ₃¹⁸. This implies that the two bases are only equivalent under the assumption of coevolution for the Eulerian parameter b_Γ3 (see Eq. 36).

The final model depends on four free parameters, the linear bias $b_1^{\mathcal{L}}$, the local quadratic bias $b_2^{\mathcal{L}}$, the tidal quadratic bias $b_{s^2}^{\mathcal{L}}$, and the higher-derivative bias $b_{{\mathrm{\nabla }}^2\delta }^{\mathcal{L}}$, to which we add the extra free parameter α_P, 1 to account (at first order) for non-Poissonian shot-noise, in the same way as it is done in the Eulerian PT model. We use a different notation for the quadratic tidal bias, since the definition of the tidal field operator is slightly different from the one presented in Sect. 3.1. The same is true for the Laplacian bias, which in this case only models higher-derivative corrections, but could also (partially) absorb unmodelled non-local effects coming from higher orders, extra physics, such as baryonic effects, or the smoothing of the density field performed in Lagrangian space.

The galaxy power spectrum can then be expressed as

$$\begin{array}{c}\hfill P_{\mathrm{gg}}(k)={\displaystyle \sum _{i,j}}{b}_i^{\mathcal{L}}{b}_j^{\mathcal{L}}{P}_{\mathit{ij}}(k)+\frac{1+{\alpha }_{P,1}}{{\overline{n}}_{\mathrm{g}}},\end{array}$$(40)

where P_ij(k) are the 15 cross-spectra of the five previously defined advected operators. To compute the P_ij terms, baccoemu has been trained with high-resolution P_ij measurements from 800 combinations of cosmologies and redshifts, obtained applying the cosmology-rescaling technique to four main N-body simulations (Angulo & White 2010; Zennaro et al. 2019; Contreras et al. 2020).

As a final remark, we note that, even if it is possible to find a relation between the Lagrangian and Eulerian bias parameters, the two sets do not exactly correspond to the same physical quantities. This happens because in the purely perturbative Eulerian framework they properly represent the response of galaxy formation to large-scale perturbations, whereas in the hybrid Lagrangian one this physical meaning is lost due to the extrapolation of the individual operators to the non-linear regime (since the individual operators of the bias expansion are measured from N-body simulations and thus contain higher-order contributions).

4. Model selection and fitting procedure

In this section we describe the methodology used to determine the best combination between different models, scale cuts, and bias configurations. In addition we list the details of the fitting procedure and the priors of the selected parameter spaces.

4.1. Performance metrics

In the context of model selection, the most relevant aspects to take into consideration are the range of validity of a given model and the precision and accuracy of the constraints on the parameters of interest. The procedure that we adopt is based on the selection of the maximum wave mode k_max up to which the model is still capable of providing a good description of the data vectors, while still recovering the correct input parameters. This can be quantified by means of three different performance metrics (employed in e.g. Osato et al. 2019; Eggemeier et al. 2020, 2021; Pezzotta et al. 2021), which are described in the next subsections.

4.1.1. Figure of bias

One of the main requirement that the theoretical model has to satisfy is that its fit to the data return unbiased model parameters. The parameters controlling bias, shot-noise, and counterterms can be effectively treated as free nuisance parameters, to be marginalised over after sampling the joint posterior distribution. The set of parameters of interest is therefore restricted to the cosmological parameters, in our case θ ≡ {h, ω_c}.

We quantify the unbiasedness of the model in the recovery of θ in terms of the figure of bias (FoB) defined as

$$\begin{array}{c}\hfill \mathrm{FoB}(\mathit{\theta })\equiv {[(\langle \mathit{\theta }\rangle -{\mathit{\theta }}_{\mathrm{f}id}{)}^{⊺}{S}^{-1}(\mathit{\theta })(\langle \mathit{\theta }\rangle -{\mathit{\theta }}_{\mathrm{fid}})]}^{\frac{1}{2}},\end{array}$$(41)

where ⟨θ⟩ and θ_fid represent the mean of the posterior distribution of the selected parameters and their fiducial values, respectively, and S(θ) is a square matrix containing the auto- and cross-covariance among all the entries of the vector θ¹⁹. The meaning of Eq. (41) is straightforward: we are quantifying the deviation of the posterior distribution from the fiducial values of the corresponding parameters, and expressing this information in terms of the intrinsic error of those parameters. In the case where θ consists of only one parameter, the FoB simply expresses how far the posterior is from the fiducial value in units of the standard deviation of the parameter, with the 68% and 95% percentiles corresponding to values of FoB of 1 and 2, respectively. Note that when considering more than one parameter these values change, as they need to be computed by directly integrating a multivariate normal distribution with the corresponding number of dimensions. For n = 2, we evaluate that the new thresholds for the 68% and 95% percentiles are 1.52 and 2.49, respectively.

4.1.2. Goodness of fit

The goodness of fit quantifies the consistency of the theoretical model P^th with the input data vector P^data. We consider the standard χ² test, corresponding to

$$\begin{array}{c}\hfill {\chi }^2(\mathit{\theta })={\displaystyle \sum _{i=1}^{N_{\mathrm{bins}}}\sum _{j=i}^{N_{\mathrm{bins}}}}[P_i^{\mathrm{th}}(\mathit{\theta })-P_i^{\mathrm{data}}]C_{\mathit{ij}}^{-1}(\mathit{\theta })[P_j^{\mathrm{th}}(\mathit{\theta })-P_j^{\mathrm{data}}].\end{array}$$(42)

This results in a distribution of χ² values across the sampled parameter space. Instead of picking the χ² corresponding to the maximum-likelihood position, whose estimation from the sampled posterior distribution is subject to noise²⁰, we compute the posterior-averaged value, ⟨ χ²⟩, from a weighted average over all sampled parameter combinations, which is instead a more stable quantity (see Appendix E for a comparison between the two approaches). The posterior-averaged χ² is then compared to the predictions from the 68% and 95% percentiles of the χ² distribution with the corresponding number of degrees of freedom. The latter is simply defined as N_dof = N_bins − N_pars, where N_bins is the total number of independent wave mode bins up to the selected k_max, and N_pars is the total number of free parameters of the model.

4.1.3. Figure of merit

Finally, each configuration of the model – that is a given scale cut and bias assumptions – is inspected to determine its statistical power in constraining the parameters θ. For this purpose, similarly to what is done for the figure of bias, we define a figure of merit (FoM) for a given set of model parameters θ as (Wang 2008)

$$\begin{array}{c}\hfill \mathrm{FoM}(\mathit{\theta })={[det(S(\mathit{\theta }))]}^{-1/2},\end{array}$$(43)

where S(θ) is once again the covariance matrix of the parameters θ, and det(S) its determinant. The meaning of this quantity can be more clearly understood assuming a flat posterior distribution with null correlation between the entries of θ. In this case, det(S) represents the volume of the hyper-rectangle over which the posterior distribution of θ is distributed. Similarly, for non-zero parameter correlations, det(S) represents the hyper-volume contained in the hyper-surface defined by the covariance matrix S. Therefore, a high value of the FoM corresponds to a more statistically significant constraint of the model parameters.

In order to visualise how much can be gained by pushing the model to higher k_max values, in the next section we plot the FoM of each individual configuration normalised to that of a reference case, corresponding to the configuration with the EFTofLSS model at k_max = 0.1 h Mpc⁻¹ with all nuisance parameters sampled as free parameters.

4.2. Fitting procedure

In order to properly sample the posterior distribution we need to compute the galaxy power spectrum and the likelihood for a large number of points in parameter space. To achieve convergence while keeping the number of evaluations as low as possible, an efficient sampling algorithm is needed.

All the results presented in this work have been obtained using a nested sampling approach (Skilling 2006), which differs from a standard Metropolis–Hastings (Metropolis et al. 1953; Hastings 1970) Markov-chain sampler in a number of ways. The main difference is that, using nested sampling, the whole hyper-dimensional parameter space is explored within the specified priors by means of a given number of live points, which are subsequently modified to track the posterior distribution of the parameters according to the value of the evidence. In this analysis we make use of the public code PyMultiNest (Buchner et al. 2014) with a total number of 1800 live points, after having checked that the output posterior distribution has properly converged with this number. Further details, together with a comparison of different samplers, are presented in Appendix F.

We adopt the approach of a full-shape analysis. This means that we directly sample the cosmological parameter space, with a model galaxy power spectrum that is generated at each step. A single evaluation of the theory models presented in Sect. 3 can take up to few seconds, since it combines a call to the Boltzmann solver to obtain linear theory predictions, and a call to the routines responsible for computing the non-linear corrections. Since the typical number of model evaluations for a single Markov chain can reach order of 𝒪 (10⁶), the final running time necessary to obtain a converged posterior distribution can take up to several days.

In order to speed up the model evaluation, we make use of the publicly available COMET package (Eggemeier et al. 2022)²¹ to emulate the EFTofLSS model, providing an evaluation of the full one-loop prediction in about 𝒪 (10 ms). The code has been validated against a set of 1500 theory data vectors in a range of redshifts that covers the one we explore in this analysis, showing an averaged 0.1% systematic error for the final P_gg(k) model, and it is therefore suited to be used for this analysis²². The evaluation of the hybrid Lagrangian-bias-based model is instead carried out using the public emulator baccoemu, as mentioned in Sect. 3.2.

In all the cases, we assume a Gaussian likelihood function defined as

$$\begin{array}{c}\hfill -2ln\mathcal{L}(\mathit{\theta })=\sum _{i=1}^{N_{\mathrm{bins}}}\sum _{j=i}^{N_{\mathrm{bins}}}[P_i^{\mathrm{th}}(\mathit{\theta })-P_i^{\mathrm{data}}]C_{\mathit{ij}}^{-1}(\mathit{\theta })[P_j^{\mathrm{th}}(\mathit{\theta })-P_j^{\mathrm{data}}],\end{array}$$(44)

which is computed at each point in parameter space explored by the sampler, and whose value is used to determine whether to assign to the current point one of the live points. The final output, which is saved to external files ready to be post-processed, consists of a list of points in parameter space together with the corresponding value of the log-likelihood.

4.3. Parameter priors

Our parameter space consists of both cosmological and nuisance parameters. Sampled cosmological parameters comprise the Hubble parameter h and the cold dark matter density parameter ω_c. The latter can be constrained only through the full-shape of the galaxy power spectrum, especially via the position of the matter-radiation equality k_eq, since geometric distance information is lost due to the fact that we conduct our analysis using real-space coordinates. For the same reason, h can be artificially constrained because we fix all the other parameters affecting the amplitude of the matter power spectrum (Sánchez et al. 2022)²³. We keep fixed the baryon density parameter ω_b and the scalar spectral index n_s, since galaxy clustering measurements on their own are not able to constrain them with the same level of precision of CMB data. At the same time, since in real space the primordial scalar amplitude A_s is strongly degenerate with the linear bias parameter b₁, at least on sufficiently large scales²⁴, we keep A_s fixed to its fiducial value, along with the rest of the cosmological parameters, to the values shown in Table 1.

The nuisance parameters are split into two sets, depending on the considered model. The parameters of the Eulerian bias expansion are composed of a mixture of bias parameters, $\{b_1,b_2,b_{{\mathcal{G}}_2},b_{{\mathrm{\Gamma }}_3}\}$, counterterms, {c₀}, and shot-noise parameters, {α_P, 1,α_P, 2}. All of them enter in the final expression for the galaxy power spectrum as shown in Sect. 3.1. When testing the bias relations presented in Sect. 3.1.3, the parameters subject to the bias relations are not sampled over, but computed at each step in the chain as a function of the lower-order bias parameters. The scale-dependent noise parameter α_P, 2 is kept fixed to 0 for the majority of the runs we carry out, except for the ones presented in Sect. 5.3, where we explicitly test the constraining power of the EFTofLSS model on this parameter in the range of redshifts that we are considering.

For the hybrid Lagrangian model we sample over a different set of bias parameters, $\{b_1^{\mathcal{L}},b_2^{\mathcal{L}},b_{s^2}^{\mathcal{L}},b_{{\mathrm{\nabla }}^2\delta }^{\mathcal{L}}\}$, and shot-noise, {α_P, 1}. In this case we do not consider relations among bias parameters, but every run will assume the full set.

When not mentioned otherwise, we adopt a completely agnostic approach, setting an uninformative flat prior for all the parameters, as shown in Table 3. The size of the prior for the two cosmological parameters and for most of the nuisance parameters has been selected to prevent the posterior distribution from becoming dominated by the imposed prior.

5. Results

In this section we present the results obtained by fitting the data samples presented in Sect. 2 with the two theoretical models described in Sect. 3. We start off with a comparison between the performance of these models in Sect. 5.1, and leave to later sections a more detailed description of the model selection carried out for the EFTofLSS model in terms of scale cuts and bias relations. For compactness, we call this model simply the EFT model.

5.1. Performance of Eulerian and hybrid Lagrangian bias expansion

In this section we carry out a comparison between the Eulerian expansion and the hybrid Lagrangian bias approach implemented in baccoemu in terms of the three performance metrics previously defined in Sect. 4.1. For this goal, we focus on fitting the galaxy power spectra of the Model 3 HOD sample at the redshifts of the four comoving snapshots, using both the rescaled Euclid-like covariances, and the ones from the full simulation box. For each case, we run multiple chains to assess the stability of the results as a function of the maximum wave mode k_max. The latter is selected in the range [0.10,0.45] h Mpc⁻¹ using a linear spacing of Δk = 0.05 h Mpc⁻¹, for a total of eight different cases.

We select two different configurations of the EFT model. The first one corresponds to the case in which all the nuisance parameters $\{b_1,b_2,b_{{\mathcal{G}}_2},b_{{\mathrm{\Gamma }}_3},c_0,{\alpha }_{P,1}\}$ are free to vary, with the only exception of the scale-dependent shot-noise parameter α_P, 2 which we set to zero. In the second one we additionally fix the cubic tidal bias b_Γ3 to its coevolution relation (Eq. (36)). The latter case is chosen in order to provide an alternative model based on the Eulerian expansion of Sect. 3.1, with the same assumptions on galaxy bias as in baccoemu (see discussion in Sect. 3.2). In addition, this is one of the best configurations when considering the performance metrics on the combination {h,ω_c}, as we properly validate in Sect. 5.4.

As for baccoemu, we leave all bias parameters $\{b_1^{\mathcal{L}},b_2^{\mathcal{L}},b_{s^2}^{\mathcal{L}},b_{{\mathrm{\nabla }}^2\delta }^{\mathcal{L}}\}$ free to vary, with the addition of the parameter controlling the amplitude of the non-Poissonian stochastic noise, α_P, 1. Since baccoemu is an emulator based on N-body simulations, it is affected by two sources of noise: first, the emulation error, that is the noise introduced by the accuracy of the trained neural network itself; second, the training set error, that is the inaccuracies already present in the data used for training. The former is a scale-dependent quantity, which becomes progressively larger at small scales and caps at a maximum 0.5% of the galaxy power spectrum signal at k ∼ 0.7 h Mpc⁻¹ for ΛCDM cosmologies well within the allowed parameter space; it can get to the order of 𝒪 (1%) of the power spectrum signal for cosmologies closer to the limits of the emulator parameter space (Zennaro et al. 2023). On the other hand, the intrinsic error of the training set is induced by the cosmology-rescaling technique employed during its construction; it once again depends on scale, and is subpercent in the case of ΛCDM cosmologies, but could reach percent levels when also massive neutrinos and dynamical dark energy are considered (Contreras et al. 2020; Zennaro et al. 2023). To account for these combined effects, we consider two cases for the chains run with baccoemu. In the first one we employ the same covariance matrix used to analyse the data galaxy power spectra as in the EFT chains, while in the second one we add in quadrature a theory error corresponding to 0.5% of the galaxy power spectrum signal²⁵.

In Fig. 3 we show the performance metrics (FoB, FoM) in the case of the realistic Euclid-like volume. Note that we do not show the averaged reduced χ² in this case, since the rescaled covariance matrix does not describe the fluctuations in the data vector, and therefore the collection of ∑(P^th − P^data)²/σ² values deviates from a χ² distribution. At all redshifts the fits obtained with both models display a FoB within the 68% confidence interval up to k_max = 0.45 h Mpc⁻¹, with only a partial preference for the EFT framework when considering the value of the FoB, which is anyway consistent to 1σ for all the cases.

In the figure, the FoM is normalised by the value obtained with the EFT configuration at k_max = 0.1 h Mpc⁻¹, to show relative gains. As expected, the FoM of the EFT model is larger with a fixed b_Γ3, because of the smaller number of free parameters. Similarly, we note that the combined constraining power on (h, ω_c) of baccoemu is degraded when including theory errors in the data covariance matrix.

When comparing the two different models, we note that baccoemu without including theory errors (orange line) reaches its maximum FoM value, comparable to the maximum value achieved by the Eulerian model (blue line), already at a lower k_max, of about 0.25–0.3 h Mpc⁻¹. This is a consequence of the extra parameter, b_Γ3, present in the Eulerian bias model. As further evidence, the Eulerian bias model display higher FoM values when considering a fixed b_Γ3 (light-blue line), in particular on scales k_max < 0.2 h Mpc⁻¹. Above this threshold, we note that the EFT configuration features a slightly larger FoM than the one of the hybrid model, with the only exception of the z = 0.9 snapshot, for which the two curves have a similar amplitude at all scales (light blue vs. orange).

Except for the main EFT case, including smaller scales does not seem to increase the FoM beyond a scale of about 0.3 h Mpc⁻¹. Since on these scales the theory error associated to baccoemu is of similar magnitude as the data covariance, we note that, including the extra 0.5% contribution (red line), the FoM starts flattening at a slightly lower k_max ∼ 0.25 h Mpc⁻¹. This is mostly noticeable for the z = 0.9 snapshot, for which shot-noise becomes the dominant contribution at a much larger k_max.

A plot similar to the one in Fig. 3 is shown in Fig. 4, this time considering the covariance matrix corresponding to the full simulation volume of about 54 h⁻³ Gpc³. Since now the covariance matrix correctly represents the statistical fluctuations in the data vectors, we additionally show the χ² averaged over the chain and normalised to the numbers of degrees of freedom. In this case, it is clear that not accounting for the theory error of baccoemu can lead to a bias in the cosmological parameters, most notably at low redshift. On the contrary, including the reference 0.5% theory error is enough to recover unbiased results, with the sole exception of the case at k_max = 0.45 h Mpc⁻¹ and z = 0.9.

The EFT model also returns unbiased measurements, with some spurious configurations outside the 1σ confidence interval for low k_max values at z = 1.5. The main reason for this effect is likely imputable to the presence of projection effects when marginalising the posterior distribution in the {h,ω_c} plane, as we explain later in Sect. 5.4. The averaged χ² behaves in a consistent way between the two models, displaying an amplitude that is constantly lower than the 95th percentile of the corresponding χ² distribution for both sets of curves, with the only exception of the largest k_max values of the z = 1.8 snapshot.

In terms of goodness of fit, we note that, due to our choice of reporting the posterior-averaged χ² value instead of the maximum-likelihood value, the normalised χ² can start off with values larger than 1 at low k_max values. This is mostly caused by the non-gaussianity of the sampled posterior distribution when the data vectors cannot properly constrain the whole set of sampled parameters. In this case the averaged posterior can increase the value of the χ² and making it appear artificially larger. In Appendix E we include an example using the maximum-likelihood χ², showing how in this case the goodness of fit typically assumes values consistent with 1 on those scales.

While the constraining power of baccoemu is in this case limited by the theory error being of similar order as the statistical error of the synthetic data considered, it is highly competitive with the Eulerian approach on scales that are free from this limitation, at k_max ≲ 0.2 h Mpc⁻¹. On the one hand, the full-volume test considered here leads to very conservative results: the errors associated with the full volume of the Flagship simulation are roughly a factor 2 smaller than the scaled errors considered in this work, and these scaled errors for Model 3, in turn, are roughly another factor 2 smaller than the errors expected assuming the volumes and number densities for typical redshift bins of the spectroscopic sample described in the forecasts of Euclid Collaboration (2020). In addition, all covariance matrices are computed in the Gaussian approximation, which might underestimate the amplitude of the errors at small scales. On the other hand, these results provide a motivation to further reduce the noise associated with emulators – for example through larger training sets, and employing Zeldovich control variates (Chartier et al. 2021; Kokron et al. 2022). This is key for the design of the next generation of emulators. In fact, even if the configuration without theory errors shows the limitation of currently available codes, the corresponding FoM curve highlights the potential gain achievable with a more accurate version of the emulator.

5.2. Testing the EFT model: Fixed cosmology

In the rest of this section we focus on testing the range of validity of the EFT model using different scale cuts, bias relations, and reference volumes. This test is limited only to the EFT model because, as shown in the previous subsection, we cannot run fits using baccoemu without accounting for the extra contribution from theory errors in the covariance matrix, especially when considering the extremely large precision of the full-box covariance. In order to assess the level of accuracy of the EFT model and determine its range of validity, we first carry out fits at fixed cosmology assuming the full volume of the simulation box. In this way, we focus exclusively on the performance of one-loop galaxy bias prediction with highly precise measurements, testing which scale cuts and bias relations lead to the best agreement between the theory model and the input data vectors.

The validation of the model includes an accuracy test consisting in recovering fiducial values for the linear bias b₁ and the shot-noise parameter α_P, 1 determined from the large-scale limit of the ratio between the measurements of the galaxy and of the matter power spectrum. This will reduce the effect of cosmic variance on the linear bias estimate (the cross galaxy-matter power spectrum is unfortunately not available). At large scales, we can assume a simple two-parameter, linear model given by

$$\begin{array}{c}\hfill P_{\mathrm{gg}}(k)=b_1^2{P}_{\mathrm{mm}}(k)+\frac{1+{\alpha }_{P,1}}{\overline{n}},\end{array}$$(45)

to be fit on scales k_max < 0.08 h Mpc⁻¹. As a reference, the marginalised mean posterior values of both b₁ and α_P, 1 are listed in Table 4.

In Fig. 5 we show the marginalised constraints obtained fitting the full model of Eq. (34) to the eight data vectors, against the fiducial values of {b₁,α_P, 1} obtained from the large-scale limit as in Eq. (45). 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₁, local quadratic bias b₂, non-local quadratic bias $b_{{\mathcal{G}}_2}$, non-local cubic bias b_Γ3, matter counterterm and higher-derivative bias c₀, and constant shot-noise parameter α_P, 1.

• (ii)

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

• (iii)

  b_Γ3 is fixed to the coevolution relation defined in Eq. (36), for a total number of five parameters.

• (iv)

  Both $b_{{\mathcal{G}}_2}$ and b_Γ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 α_P, 2 fixed to zero (we test the validity of this assumption in Sect. 5.3).

Overall, the configuration with the largest number of free parameters – case (i) in the previous list, shown with blue points and error bars in Fig. 5 – is capable of capturing the correct amplitude of both b₁ and α_P, 1 for the majority of the tested k_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_max > 0.2 h Mpc⁻¹. 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σ. 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 non-linear degeneracies among the model parameters.

Fixing either $b_{{\mathcal{G}}_2}$ – case (ii), orange points and error bars – or b_Γ3 – case (iii), green points and error bars – does not significantly help in terms of accuracy of the marginalised constraints, with systematic deviations that can still become larger than the 1σ confidence interval. However we find that, while fixing b_Γ3 typically results in similar constraining power on both b₁ and α_P, 1, imposing a relation on $b_{{\mathcal{G}}_2}$ leads to definitely tighter posteriors. This is the result of breaking the strong degeneracy between the two non-local bias parameters, $b_{{\mathcal{G}}_2}$ and b_Γ3, and at the same time the one between the quadratic biases, b₂ and $b_{{\mathcal{G}}_2}$, leaving the remaining parameters to be more tightly constrained (a clear example of these is displayed in the right panel of Fig. F.1). The same clearly happens when combining the two previous relations – case (iv), red points and error bars – 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₁ which can reach more than 2σ 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₂-$b_{{\mathcal{G}}_2}$-b_Γ3 degeneracy can be observed in a more direct way by inspection of the 2d marginalised constraints in the b₁-α_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_{{\mathcal{G}}_2}$ fixed to the excursion-set-based relation, and with b_Γ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_{{\mathcal{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_max = 0.4 h Mpc⁻¹. 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_{{\mathcal{G}}_2}$ and b_Γ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_max bin of the Model 3 sample at z = 0.9.

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_V value, defined as (Gelman et al. 2014)

$$\begin{array}{c}\hfill p_{\mathrm{V}}=\frac{1}{2}\langle {({\chi }^2-\langle {\chi }^2\rangle )}^2\rangle ,\end{array}$$(46)

that is, the variance of the corresponding χ² 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_max. In practice, we observe that the model with all parameters free never reaches the expected value of p_V = 6, even for the largest value of k_max, with the exception of a couple of configurations. This shows that 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_max, with a transition that typically happens sooner for the case with fixed $b_{{\mathcal{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 G 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 that most of the residuals observed in Fig. 5 can be explained by sample variance affecting our data vectors.

5.3. Constraints on scale-dependent shot noise

So far, the stochastic field ε_g entering the expression for the galaxy density field in Eq. (20) has been assumed responsibile only for a constant offset from Poissonian predictions, via the parameter α_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 ε_g. As already mentioned in Sect. 3, this leads to the presence of an additional k²-dependent term in the galaxy power spectrum, whose amplitude is regulated by the extra parameter α_P, 2.

Figure 8 shows the marginalised one-dimensional constraints on α_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 α_P, 2, we only consider values of k_max above 0.35 h Mpc⁻¹²⁶. We never observe a statistically significant detection of the α_P, 2 parameter, with the majority of the marginalised constraints being consistent with α_P, 2 = 0 well within the 2σ confidence interval. The only configurations for which this does not happen are the ones at high redshifts, specifically when considering high values of k_max, since these are the configurations for which the parameter α_P, 2 is constrained with the highest precision. Performing the same test with one of the tidal bias fixed²⁷ does not lead to significantly different conclusions. This seems to suggest that the α_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α 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²-dependent noise corrections are required, as a function of the orientation with respect to the line of sight (Philcox & Ivanov 2022; Carrilho et al. 2023; Moretti et al. 2023). Finally, we note 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 {h,ω_c}, 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

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. Specifically, 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 ω_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_s, since the anisotropies introduced by peculiar velocities make possible to break the strong A_s- b₁ degeneracy that is otherwise present when considering real-space coordinates²⁸. 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_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 RSDs (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.

5.4.1. Figure of bias

In terms of FoB, we observe a consistent trend across each model configuration, indicating an unbiased combined measurement of the cosmological parameters even well within the mildly non-linear regime, at k_max ≳ 0.3 h Mpc⁻¹. The only exception is represented by the configuration in which both tidal bias parameters, $b_{{\mathcal{G}}_2}$ and b_Γ3, are simultaneously kept fixed to the excursion-set relation (Eq. (37)) and to the coevolution relation (Eq. (36)), respectively. This outcome is unsurprising, as we are reducing by two the total number of degrees of freedom of the model. Notably, the bias on the cosmological parameters gets larger at lower redshift, hinting to a departure from the coevolution relations as non-linear gravitational effects become more pronounced. Typically, this deviation occurs at scales of approximately k_max = 0.3 h Mpc⁻¹. However, we observe that this configuration behaves surprisingly well for high-redshift snapshots, exhibiting a FoB well within the 68% confidence interval.

In terms of overall stability of the results, we observe a deviation at low k_max values for the Model 3 sample. This is more strongly affecting high-redshift snapshots, for which the value of the FoB at k_max ∼ 0.1 h Mpc⁻¹ already exceeds the corresponding 68% confidence level, and only gets below the threshold when including additional signals from smaller scales. This effect is primarily attributed to the presence of projection effects, owing to the large dimensionality of the selected parameter space. Specifically, we find that all samples display a non-negligible correlation between the cold dark matter density parameter, ω_c, and the EFT counterterm c₀, resulting in a systematic shift of ω_c for the lowest values of k_max, where there is insufficient constraining power to accurately constrain c₀. Fixing the cubic tidal bias b_Γ3 to the coevolution relation typically helps to restore the cosmological parameters to their fiducial positions. This happens due the further degeneracy between b_Γ3 and c₀ over the mildly non-linear regime. The relative importance of the ω_c-c₀ 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 ≲ z ≲ 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.

5.4.2. Goodness of fit

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_{{\mathcal{G}}_2}$ and b_Γ3, are kept fixed to their corresponding relations and for the largest k_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 χ² value at k_max ≳ 0.3 h Mpc⁻¹, 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 χ² value that is typically well within the 95% percentile of the corresponding χ² distribution.

When considering the two different HOD models, the χ² 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 χ² spuriously gets larger than the 95% confidence interval for some of the selected k_max values. In practice, this deviation is still consistent to better than the 3σ 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_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 χ² 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 that 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.

5.4.3. Figure of merit

As expected, the FoM monotonically increases when including additional information from more non-linear scales, with a relative gain with respect to the most relaxed configuration (all bias parameters free to vary, at k_max = 0.1 h Mpc⁻¹) that becomes larger moving towards lower redshifts. In fact, extending the fitting range to the maximal value of k_max = 0.45 h Mpc⁻¹, 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 non-linear 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_max ≳ 0.3 h Mpc⁻¹. 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_max value we probe. However, we note 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.

5.4.4. Summary

Overall, we find that fixing only the quadratic tidal bias $b_{{\mathcal{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_max. Also, in a range of scales up to k_max ∼ 0.3 h Mpc⁻¹, 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.

Nevertheless, as already mentioned in Sec. 3.1.3, we argue that the applicability of the excursion set and coevolution relations may partially fail when adopting them with real data, given the lack of any realistic dependency on other quantities different from the halo mass, such as assembly bias (Croton et al. 2007; Barreira et al. 2021; Hadzhiyska et al. 2021b; Lazeyras et al. 2023) or the scatter in the ⟨N⟩(M_h) relation (Behroozi et al. 2010; Zehavi et al. 2011). For this reason, extended studies on more realistic simulated samples will be needed to properly benchmark these relations when applying them to real Euclid observations.

As a final remark, in Nicola et al. (2024) the authors carried out a comparison similar to the one presented in this paper, by forecasting the impact of different bias models in the analysis of the Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) photometric observations. Differently from this work, their analysis showed that the use of minimal bias models with 1/2 less degrees of freedom could bias the recovery of cosmological parameters at more than 3σ also when restricting the analysis to relatively large scales (k_max ∼ 0.15 h Mpc⁻¹). We argue that this apparent inconsistency is mainly imputable to the different observables selected in the two analyses. While in this work we make use of the real-space galaxy power spectrum at some fixed redshifts z ≳ 1, their analysis is focused on the 3 × 2-point data combination in tomographic bins that are for the most part at lower redshifts (0 < z < 1), for which the assumption of a negligible tidal or higher-derivative bias is no longer valid. As a further evidence, the dashed black lines in Figs. 9 and 10 show how setting both tidal biases to 0 makes the Eulerian bias model break sooner when considering the snapshots at lower redshifts, especially the one at z = 0.9.

5.5. Dependence on sample volume

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, as explained in Sect. 2.3. For each comoving snapshot, at z = (0.9, 1.2, 1.5, 1.8), with reference volume V_box, we define the volume V_shell of a spectroscopic bin corresponding to a total angular surface of 15 000 square degrees – and with a depth of Δz = (0.2, 0.2, 0.2, 0.3) – following the choices made in Euclid Collaboration (2020). The three additional volume rescalings are obtained by selecting the values that divide the interval [η, 1] into four equi-partitioned subintervals, where η = V_box/V_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_Γ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 colours 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 (h,ω_c) 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, η.

It comes with no surprise that the reference bias model is well-performing up to the highest value of k_max we consider, even under realistic assumptions. Additionally, for the case with fixed b_Γ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_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_Γ3 to the coevolution relation.

As a final consideration, we note that 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 {h,ω_c} 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 targets of Euclid. We employed a set of four FoF halo catalogues from comoving snapshots of the Flagship I simulation at z = (0.9, 1.2, 1.5, 1.8), which were subsequently populated with Hα 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⁻¹ Mpc)³ and a high comoving number density (from ∼10⁻⁴ h³ Mpc⁻³ to ∼10⁻³ h³ Mpc⁻³), which corresponds to a flux limit of f_Hα = 2 × 10⁻¹⁶ erg cm⁻² s⁻¹. These snapshots 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 a Eulerian bias expansion, and 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 counter terms, 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 ω_c, plus a set of six nuisance parameters, consisting of the linear bias b₁, the quadratic bias b₂, the tidal quadratic and cubic biases, $b_{{\mathcal{G}}_2}$ and b_Γ3, the matter counter term and higher-derivative bias c₀, and two extra parameters representing deviations from Poissonian shot noise: a constant offset, α_P, 1, and a scale-dependent term, α_P, 2. The second model adopts a similar one-loop expansion of the galaxy power spectrum, but this time using Lagrangian coordinates, and is based on the emulation of the individual terms of the 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). In addition to the two cosmological parameters, this model features the linear bias b₁, the quadratic bias b₂, the tidal quadratic bias b_s², the Laplacian bias b_∇²δ, and a stochastic parameter representing a constant deviation from Poissonian shot noise, α_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 present tests of the relative performance of these two galaxy bias models, while in the subsequent sections we show how 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 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 latter two metrics are computed on the {h,ω_c} combination in order to quantify the accuracy and precision of the model 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 000deg²). Our results highlight how both models are capable of providing unbiased measurements of the cosmological parameters up to k_max = 0.45 h Mpc⁻¹ for all four redshift snapshots, consistently with the 1σ 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_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σ confidence interval. However, there are some spurious cases at low values of k_max that are affected by projection effects, which can be alleviated by fixing the value of b_Γ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_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σ confidence interval. Also in this case, on scales for which the intrinsic error of baccoemu is not the dominant contribution, that is, k_max ≲ 0.2 h Mpc⁻¹, its FoM is consistent with that of the EFT model. Not including the 0.5% extra error, we note that the FoM is consistent with – and in some cases even larger than – the corresponding EFT results, demonstrating the potential gain to be obtained 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 describe the fits we performed at fixed cosmology to check the self-consistency of the one-loop galaxy bias expansion in terms of the linear bias b₁ and the scale-independent shot-noise parameter α_P, 1. The fiducial values for these parameters were fitted from the measured galaxy-to-matter-power-spectrum ratio, assuming a leading-order recipe on scales of k < 0.08 h Mpc⁻¹. The result is that, when leaving all parameters free to vary, it is possible to recover – at better than 2σ – the value of both parameters for the majority of the samples – two HOD models times four different redshifts – and scale cuts up to k_max = 0.45 h Mpc⁻¹. 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_max = 0.2 h Mpc⁻¹. The latter is however consistent with sample variance expectations, as observed in a set of ten different noisy realisations of a synthetic theory data vector (see Appendix G). 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σ 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₂ and the cubic tidal bias b_Γ3. We find that fixing both parameters at the same time still works extremely well (with a typical recovery of b₁ and α_P, 1 within the 68% confidence interval) even for the largest k_max values when considering the high-redshift snapshots. However, at low redshift, this choice can lead to deviations of more than 3σ for some of the configurations we test, especially when the maximum scale included in the fit is above a typical scale of k_max ∼ 0.3 h Mpc⁻¹.

We explicitly checked whether or not a next-to-leading-order correction to the shot-noise contribution (α_P, 2) can improve the model performance. In all cases considered, the marginalised posterior distribution for α_P, 2 is consistent with zero within 2σ, 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 RSDs and observational systematic effects.

When the parameter space is extended to also include the two cosmological parameters, we note 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 values of k_max, which is due to projection effects when marginalising over all the nuisance parameters. This trend can be partially corrected by fixing the cubic tidal bias b_Γ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 (h,ω_c) 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, which suggests a premature breaking of the tested relations. In terms of goodness of fit, we do not observe a significant change in the average χ² 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, the configuration with both tidal biases fixed results in biased constraints for the low-redshift samples, while the case with only b_Γ3 fixed can be employed down to smaller scales. This configuration displays 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_max ≲ 0.3 h Mpc⁻¹, where the model with fixed b_Γ3 experiences a steep increase that is then followed by a more modest growth. The configuration with $b_{{\mathcal{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_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 ≲ z ≲ 2, even deep within the mildly non-linear regime, at k_max ∼ 0.4 h Mpc⁻¹, 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 therefore 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 the 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, which 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 × 2-point data combination. In parallel, 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 papers in the series (Euclid Collaboration: 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.).

---

Acknowledgments

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’Avenir” programme (AMX-19-IET-008 – IPhU). This research made use of matplotlib, a Python library for publication quality graphics (Hunter 2007).

References

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 [0.7,2.4]. Each measurement consists of the matter power spectrum, measured in the range of wave modes defined by the interval [0.01,4]  h Mpc⁻¹, using 18 linearly spaced bins up to k ∼ 0.03 h Mpc⁻¹ 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. 19) and the second using the non-linear 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, that is, (3780 h⁻¹ Mpc)³. For both models, we limit the maximum mode of the fit to k_max = 0.25 h Mpc⁻¹, within the expected PT range of validity for the relevant redshifts.

In Fig. A.1 we show the marginalised 1d posterior distributions for {h,ω_c,A_s} as a function of redshift. For this test we limit the redshift range to the first four snapshots (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σ 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σ with each other, pointing to a systematic effect that cannot be attributed to the particular model used to describe the data vectors.

In the upper panel of Fig. A.2 we show the ratio of the non-linear matter power spectrum from the baccoemu emulator assuming the nominal cosmology to the measured one for the different redshifts. We note that 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_s = 0.96. Indeed we note that the disagreement between nominal and recovered cosmology can be alleviated for all parameters by assuming a different value of n_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 σ₈). 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 {ω_c, ω_b, M_ν} 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, that is, the scalar amplitude A_s and index n_s. The best-fit value that we found for the spectral index is close to the value of n_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 ∼ 4 h Mpc⁻¹, as shown in the bottom panel of Fig. A.2.

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 non-linear density and velocity fluctuations – δ(k) and θ(k) – in an expanding universe can be expressed in terms of linear theory solutions, δ_L(k). Assuming an EdS universe, it is possible to perfectly separate the time- and space-dependence of δ(k) and θ(k) (Goroff et al. 1986; Jain & Bertschinger 1994). In an arbitrary Λ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

$$\begin{array}{cc}& \delta (\mathit{k},\tau )=\sum _{n=1}^{\infty }{D}_n(\tau ){\delta }^{(n)}(\mathit{k}),\hfill \end{array}$$(B.1)

$$\begin{array}{cc}& \theta (\mathit{k},\tau )=-\mathcal{H}(\tau )f(\tau )\sum _{n=1}^{\infty }{D}_n(\tau ){\theta }^{(n)}(\mathit{k}),\hfill \end{array}$$(B.2)

where τ is the conformal time defined via dt = a(τ) dτ, a(τ) is the cosmic scale factor, ℋ(τ)≡dlna(τ)/dτ is the conformal Hubble expansion factor, and f(τ)≡dlnD₁(τ)/dlna(τ) 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ⁿ 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

$$\begin{array}{cc}\hfill {\delta }^{(n)}(\mathit{k})=& {\displaystyle {\int }_{{\mathit{q}}_1}}\dots {\displaystyle {\int }_{{\mathit{q}}_n}}{\delta }_{\mathrm{D}}^{(3)}(\mathit{k}-{\mathit{q}}_{1\dots n})F_n({\mathit{q}}_1,\dots ,{\mathit{q}}_n)\hfill \\ \hfill & \times {\delta }_{\mathrm{L}}({\mathit{q}}_1)\dots {\delta }_{\mathrm{L}}({\mathit{q}}_n),\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill {\theta }^{(n)}(\mathit{k})=& {\displaystyle {\int }_{{\mathit{q}}_1}}\dots {\displaystyle {\int }_{{\mathit{q}}_n}}{\delta }_{\mathrm{D}}^{(3)}(\mathit{k}-{\mathit{q}}_{1\dots n})G_n({\mathit{q}}_1,\dots ,{\mathit{q}}_n)\hfill \\ \hfill & \times {\delta }_{\mathrm{L}}({\mathit{q}}_1)\dots {\delta }_{\mathrm{L}}({\mathit{q}}_n),\hfill \end{array}$$(B.4)

where the n-th order PT kernels F_n and G_n are homogeneous functions of the wave vectors (q₁,…,q_n), and are built starting from the fundamental mode-coupling functions,

$$\begin{array}{cc}& \alpha ({\mathit{k}}_1,{\mathit{k}}_2)\equiv \frac{{\mathit{k}}_{12}·{\mathit{k}}_1}{k_1^2},\hfill \end{array}$$(B.5)

$$\begin{array}{cc}& \beta ({\mathit{k}}_1,{\mathit{k}}_2)\equiv \frac{k_{12}^2({\mathit{k}}_1·{\mathit{k}}_2)}{2k_1^2{k}_2^2},\hfill \end{array}$$(B.6)

with k₁₂ = k₁ + k₂. At linear order these quantities clearly become unity, F₁ = G₁ = 1, so to recover linear theory predictions, that is, δ⁽¹⁾(k)=θ⁽¹⁾(k)=δ_L(k). At higher order, these kernels can be derived using recursive relations, which read

$$\begin{array}{cc}\hfill F_n({\mathit{q}}_1,\dots ,{\mathit{q}}_n)=& \sum _{m=1}^{n-1}\frac{G_m({\mathit{q}}_1,\dots ,{\mathit{q}}_m)}{(2n+3)(n-1)}\hfill \\ \hfill & \times [(2n+1)\alpha ({\mathit{k}}_1,{\mathit{k}}_2)F_{n-m}({\mathit{q}}_{m+1},\dots ,{\mathit{q}}_n)\hfill \\ \hfill & +2\beta ({\mathit{k}}_1,{\mathit{k}}_2)G_{n-m}({\mathit{q}}_{m+1},\dots ,{\mathit{q}}_n)],\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill G_n({\mathit{q}}_1,\dots ,{\mathit{q}}_n)=& \sum _{m=1}^{n-1}\frac{G_m({\mathit{q}}_1,\dots ,{\mathit{q}}_m)}{(2n+3)(n-1)}\hfill \\ \hfill & \times [3\alpha ({\mathit{k}}_1,{\mathit{k}}_2)F_{n-m}({\mathit{q}}_{m+1},\dots ,{\mathit{q}}_n)\hfill \\ \hfill & +2n\beta ({\mathit{k}}_1,{\mathit{k}}_2)G_{n-m}({\mathit{q}}_{m+1},\dots ,{\mathit{q}}_n)],\hfill \end{array}$$(B.8)

where k₁ = q₁ + … + q_m and k₂ = q_m + 1 + … + q_n. As a classical example, needed for the calculation of the one-loop galaxy power spectrum, the second-order PT kernels for the non-linear evolution of the matter density and velocity fields are defined as

$$\begin{array}{cc}& F_2({\mathit{k}}_1,{\mathit{k}}_2)=\frac{5}{7}+\frac{1}{2}\frac{{\mathit{k}}_1·{\mathit{k}}_2}{k_1{k}_2}(\frac{k_1}{k_2}+\frac{k_2}{k_1})+\frac{2}{7}\frac{{({\mathit{k}}_1·{\mathit{k}}_2)}^2}{k_1^2{k}_2^2},\hfill \end{array}$$(B.9)

$$\begin{array}{cc}& G_2({\mathit{k}}_1,{\mathit{k}}_2)=\frac{3}{7}+\frac{1}{2}\frac{{\mathit{k}}_1·{\mathit{k}}_2}{k_1{k}_2}(\frac{k_1}{k_2}+\frac{k_2}{k_1})+\frac{4}{7}\frac{{({\mathit{k}}_1·{\mathit{k}}_2)}^2}{k_1^2{k}_2^2},\hfill \end{array}$$(B.10)

while the explicit expression for the third-order kernel of the matter density field, F₃(k₁, k₂, k₃), can be found in Goroff et al. (1986), and it is not reported here for practical purposes.

When considering the galaxy power spectrum P_gg, we need to evaluate the usual two-point statistics defined by the ensemble average of the galaxy density field times itself,

$$\begin{array}{c}\hfill \langle {\delta }_{\mathrm{g}}(\mathit{k}){\delta }_{\mathrm{g}}({\mathit{k}}^{\prime })\rangle ={(2\pi )}^3{\delta }_{\mathrm{D}}^{(3)}(\mathit{k}+{\mathit{k}}^{\prime })P_{\mathrm{gg}}(k).\end{array}$$(B.11)

At third order in the perturbations of δ, we obtain the one-loop expression for the power spectrum presented in Eq. (27), 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_L(|k−q|) P_L(q), reflecting the mixing of modes due to non-linear evolution, while in the latter the integrals are carried out on the factor P_L(k) P_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{array}{c}\hfill \begin{array}{cc}\hfill P_{\mathrm{gg}}^{1\text{-loop}}(k)=& b_1^2{P}^{1\text{-}\mathrm{loop}}(k)\hfill \\ \hfill & +b_1{b}_2{P}_{b_1{b}_2}(k)+b_1{b}_{{\mathcal{G}}_2}{P}_{b_1{b}_{{\mathcal{G}}_2}}(k)\hfill \\ \hfill & +b_1{b}_{{\mathrm{\Gamma }}_3}{P}_{b_1{b}_{{\mathrm{\Gamma }}_3}}(k)+b_2^2{P}_{b_2{b}_2}(k)\hfill \\ \hfill & +b_2{b}_{{\mathcal{G}}_2}{P}_{b_2{b}_{{\mathcal{G}}_2}}(k)+b_{{\mathcal{G}}_2}^2{P}_{b_{{\mathcal{G}}_2}{b}_{{\mathcal{G}}_2}}(k),\hfill \end{array}\end{array}$$(B.12)

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

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P^{1\text{-}\mathrm{loop}}(k)=& P^{1\text{-loop},MC}(k)+P^{1\text{-loop},Prop}(k)\hfill \\ \hfill =& 2{\displaystyle {\int }_{\mathit{q}}}{F}_2^2(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(|\mathit{k}-\mathit{q}|)P_{\mathrm{L}}(q)\hfill \\ \hfill & +6P_{\mathrm{L}}(k){\displaystyle {\int }_{\mathit{q}}}{F}_3(\mathit{q},-\mathit{q},\mathit{k})P_{\mathrm{L}}(q),\hfill \end{array}\end{array}$$(B.13)

$$\begin{array}{c}\hfill P_{b_1{b}_2}(k)=2{\displaystyle {\int }_{\mathit{q}}}{F}_2(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(|\mathit{k}-\mathit{q}|)P_{\mathrm{L}}(q),\end{array}$$(B.14)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_{b_1{b}_{{\mathcal{G}}_2}}(k)=& P_{b_1{b}_{{\mathcal{G}}_2}}^{\mathrm{M}C}(k)+P_{b_1{b}_{{\mathcal{G}}_2}}^{\mathrm{P}rop}(k)\hfill \\ \hfill =& 4{\displaystyle {\int }_{\mathit{q}}}{F}_2(\mathit{k}-\mathit{q},\mathit{q})S(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(|\mathit{k}-\mathit{q}|)P_{\mathrm{L}}(q)\hfill \\ \hfill & +8P_{\mathrm{L}}(k){\displaystyle {\int }_{\mathit{q}}}{F}_2(\mathit{k},-\mathit{q})S(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(q),\hfill \end{array}\end{array}$$(B.15)

$$\begin{array}{c}\hfill P_{b_1{b}_{{\mathrm{\Gamma }}_3}}(k)=-\frac{16}{7}{P}_{\mathrm{L}}(k){\displaystyle {\int }_{\mathit{q}}}S(\mathit{k}-\mathit{q},\mathit{q})S(\mathit{k},\mathit{q})P_{\mathrm{L}}(q),\end{array}$$(B.16)

$$\begin{array}{c}\hfill P_{b_2{b}_2}(k)=\frac{1}{2}{\displaystyle {\int }_{\mathit{q}}}{P}_{\mathrm{L}}(|\mathit{k}-\mathit{q}|)P_{\mathrm{L}}(q),\end{array}$$(B.17)

$$\begin{array}{c}\hfill P_{b_2{b}_{{\mathcal{G}}_2}}(k)=2{\displaystyle {\int }_{\mathit{q}}}S(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(q),\end{array}$$(B.18)

$$\begin{array}{c}\hfill P_{b_{{\mathcal{G}}_2}{b}_{{\mathcal{G}}_2}}(k)=2{\displaystyle {\int }_{\mathit{q}}}{S}^2(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(\mathit{k}-\mathit{q},\mathit{q})P_{\mathrm{L}}(q).\end{array}$$(B.19)

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

$$\begin{array}{c}\hfill P_{b_1{b}_{{\mathrm{\Gamma }}_3}}(k)=\frac{2}{5}{P}_{b_1{b}_{{\mathcal{G}}_2}}^{\mathrm{M}C}(k).\end{array}$$(B.20)

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. 37) for $b_{{\mathcal{G}}_2}$, or the coevolution relations for both $b_{{\mathcal{G}}_2}$ and b_Γ3 (Eqs. 35 – 36).

The behaviour of all the loop integrals listed above is to consistently converge to zero at infrared modes, since the non-linear kernel F_n(k₁, …, k_n) scales as k² when k ≡ k₁ + … + k_n goes to zero, reflecting the range of validity of linear theory predictions. The only exception is represented by the P_b2b2 term, which features a non-zero asymptote for k → 0. This limit can be manually set to zero via a redefinition of the loop integral, such that

$$\begin{array}{c}\hfill P_{b_2{b}_2}(k)=\frac{1}{2}{\displaystyle {\int }_{\mathit{q}}}{P}_{\mathrm{L}}(q)[P_{\mathrm{L}}(|\mathit{k}-\mathit{q}|)-P_{\mathrm{L}}(q)].\end{array}$$(B.21)

In turn, the extra contribution

$$\begin{array}{c}\hfill P_{b_2{b}_2}^{\mathrm{n}oise}(k)={\displaystyle {\int }_{\mathit{q}}}{P}_{\mathrm{L}}^2(q)\end{array}$$(B.22)

can be absorbed by the constant shot-noise parameter, α_P, 1, which 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_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 of a rescaling of the original formula for the featureless matter power spectrum P_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

$$\begin{array}{c}\hfill P_{\mathrm{nw}}(k)=P_{\mathrm{EH}}(k)\mathcal{F}\left[\frac{P_{\mathrm{L}}(k)}{P_{\mathrm{EH}}(k)}\right],\end{array}$$(C.1)

where ℱ is meant to filter out the broadband difference between P_L and P_EH. We choose a functional form for ℱ corresponding to a Gaussian filter, that is,

$$\begin{array}{c}\hfill \mathcal{F}[f(k)]=\frac{{\log }_{10}(\mathrm{e})}{\sqrt{2\pi }\lambda }\int \mathrm{d}q\frac{f(q)}{q}exp[-\frac{1}{2{\lambda }^2}{log}_{10}^2\left(\frac{k}{q}\right)],\end{array}$$(C.2)

where λ 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 λ = 0.25.

The result of this approach is presented in the top panel of Fig. C.1, where we show a comparison between the shape of the original Eisenstein & Hu (1998) smooth function and the one presented in Eq. (C.1) using the fiducial cosmology from Table 1 at z = 0. From the comparison, it is clear that the broadband of P_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_L, with the only exception of a small residual of ∼ 0.5% peaking at 0.5 h Mpc⁻¹ ≲ k ≲ 1 h Mpc⁻¹.

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 is based on 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. C.1 shows a comparison between the three methods summarised above. We note 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 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 ∼ 0.01 h Mpc⁻¹, 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 Mpc⁻¹ ≲ k ≲ 1 h Mpc⁻¹, 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 ³⁰ (Hahn 2005) used to generate the training set of the COMET emulator (Eggemeier et al. 2022), and implemented in the CosmoSIS-gClust³¹ code (courtesy of A. Moradinezhad Dizgah). Alternative methods take advantage of a Fast Fourier Transform approach, such as FastPT³² (McEwen et al. 2016) – implemented in the PBJ code (Oddo et al. 2020, 2021; Rizzo et al. 2023; Moretti et al. 2023) – and FFTLog³³ (Hamilton 2000; Simonović et al. 2018) – used in the TNSToolBox³⁴ 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 Mpc⁻¹, 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_{\mathrm{mm}}^{\mathrm{IR}\text{-}(\mathrm{LO}+\mathrm{NLO})}$, which contains the contribution P₁₃, and more importantly P_b1bΓ3, for which the discrepancy can become as large as ∼ 0.05% on scales of k ∼ 1 h Mpc⁻¹. 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 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 χ² statistics. For this purpose, we simply iterate through the sampled positions of the posterior distribution and define the weighted mean as

$$\begin{array}{c}\hfill \langle {\chi }^2\rangle =\sum _i{w}_i{\chi }_i^2,\end{array}$$(E.1)

where ${\chi }_i^2$ and w_i correspond to the χ² 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 χ² getting above the 2σ in Fig. 9 and 10. Differently from those two plots, in this case we compare the trend when using the posterior-averaged χ² (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σ of the considered χ² distribution. However, the marginal gains are not significantly changing the interpretations that can be made when using the posterior-averaged values.

As a complement, in Fig. E.2 we show the complete distribution of χ² values, and we compare it with the 68th and 95th prercentile of the reference χ² distribution with the same numbers of degrees of freedom. Different panels correspond to different k_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_max = 0.4 h Mpc⁻¹. Otherwise, our conclusions are not affected by the choice of which χ² value is used to estimate the goodness-of-fit performance metric. For this reason, in Sec. 5 we show results employing the posterior-averaged χ² 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,³⁵ 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).³⁶

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.³⁷

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_max = 0.2 h Mpc⁻¹. We test two different cases, one where we vary the nuisance parameters (b₁,b₂,b_Γ3,c₀,α_P, 1) while keeping the quadratic tidal bias $b_{{\mathcal{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.³⁸

In the first case, we note 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_{{\mathcal{G}}_2}$ and b_Γ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_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₁ and α_P, 1 are fixed to the best-fit value listed in Table 4 for the Model 1 HOD samples, b₂ is computed using a b₂(b₁) relation derived from the corresponding HOD model, $b_{{\mathcal{G}}_2}$ and b_Γ3 are obtained from the excursion-set relation (37) and the coevolution relation (36), respectively, c₀ is set to unity, and α_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⁻³ Gpc³.

Figure G.1 shows the marginalised one-dimensional constraints of the b₁ and α_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_Γ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₂, $b_{{\mathcal{G}}_2}$ and b_Γ3, and cannot be broken completely even when reducing the dimensionality of the parameter space by fixing b_Γ3, as the b₂-$b_{{\mathcal{G}}_2}$ degeneracy is still present. Consistently among the different redshifts, we note that there is a trend for b₁ and α_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_{{\mathcal{G}}_2}$ is fixed as a function of b₁ 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₂-$b_{{\mathcal{G}}_2}$ and $b_{{\mathcal{G}}_2}$-b_Γ3 degeneracies, for which a clear example can be found in the left plot of Fig. F.1.

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₁ and α_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_{{\mathcal{G}}_2}$ is fixed as a function of b₁. We remark how the case at z = 1.2 features a large fluctuation at intermediate scales, k_max ∼ 0.2 h Mpc⁻¹, 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 constraints of the configurations we are considering, that is, with all the parameters free to vary, or with either $b_{{\mathcal{G}}_2}$ or b_Γ3 fixed in terms of b₁. 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₁, α_P, 1) plane at k_max ≳ 0.2 h Mpc⁻¹. 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.

All Tables

All Figures