© The Authors 2024

1. Introduction

The recently launched Euclid mission will measure the positions and redshifts of billions of galaxies, granting us an unprecedented picture of the large-scale structure of the Universe (Laureijs et al. 2011; Amendola et al. 2018; Euclid Collaboration 2024). This will serve as a fantastic experimental test of the standard model of cosmology, ΛCDM¹, and as a means to constrain or detect deviations from it. In particular, the Euclid mission will offer a unique opportunity to probe and constrain modifications in both the matter and gravitational sectors (Euclid Collaboration 2020).

Euclid can thus offer insight into a more fundamental description of gravity beyond General Relativity (GR), and perhaps into one that could account for the ongoing accelerated expansion of the Universe (Perlmutter 1999; Riess et al. 1998), one of the most poorly understood phenomena in Nature. It could also provide clues as to whether or not dark energy evolves over time, or whether or not it interacts with CDM, the latter of which can help solve emerging data-set tensions within ΛCDM (see Pourtsidou & Tram 2016, for example).

A key observable that Euclid will provide is the power spectrum – the Fourier transform of the 2-point correlation function – of the galaxy distribution in redshift space. Given the immense number of galaxies that Euclid will detect, the statistical uncertainty on the measured power spectrum will be tiny. In particular, Euclid is set to measure more than twenty million galaxy redshifts (Amendola et al. 2018)². Since there will be more galaxy pairs found at small physical separations, the statistical uncertainty will be smaller as we move to the mildly nonlinear, small-scale, regime. There is therefore, potentially, a wealth of cosmological and gravitational information to be extracted by going beyond the linear, large-scale regime (see, for example, Bernardeau et al. 2002; Nishimichi et al. 2007; Lacasa 2022; Bose et al. 2019).

In order to seize this opportunity, the precision of the data must be matched with the accuracy of the theoretical predictions, which becomes a challenge to maintain when nonlinear effects become important. Further, if modelling uncertainty is not well estimated, then a biased picture of the Universe will be inferred from the data (see, for example, Markovic et al. 2019; Martinelli et al. 2021).

There exist various approaches to modelling large-scale structure observables on scales where nonlinear effects become significant. N-body simulations are powerful tools for this purpose. However, they are not well-suited for data-theory inference analyses due to their high computational cost. Another alternative is to employ faster, semi-analytical methods based on standard perturbation theory (SPT, see Bernardeau et al. 2002, for a review), such as the effective field theory of large-scale structure (EFTofLSS; Baumann et al. 2012; Carrasco et al. 2012) or other related approaches (see for example Pueblas & Scoccimarro 2009; Pietroni et al. 2012). A procedure one can follow is to model all relevant physics within the chosen perturbative framework and then determine the range of scales where this approximate method remains valid – where the theoretical uncertainty is smaller than the statistical and known systematic uncertainties – by comparison with an N-body benchmark measurement (see, for example, Heitmann et al. 2016, 2019; Bose et al. 2019, 2017; Rossi et al. 2021).

In this paper, we aim to quantify the uncertainty of various approximations and approaches to the nonlinear redshift space galaxy clustering power spectrum in beyond-ΛCDM cosmologies. The goal is to identify which such choices are worthy of further assessment and eventual use in computationally demanding, statistical analyses.

In particular, we have considered four representative theories of gravity and dark energy that have been identified by the Euclid consortium as target theories to test (Amendola et al. 2018): the normal branch of the Dvali–Gabadadze–Porrati (DGP) model of gravity (Dvali et al. 2000), f(R) gravity (Carroll et al. 2004; Hu & Sawicki 2007), a time-varying dark energy component (Chevallier & Polarski 2001; Linder 2003), and an interacting dark matter – dark energy model with pure momentum exchange (Simpson 2010; Pourtsidou et al. 2013). Two gravity models provide examples of phenomenologically distinct types of screening³, these being the Vainshtein mechanism (Vainshtein 1972, in the case of DGP) and the chameleon mechanism (Khoury & Weltman 2004, in the case of f(R) gravity). These two gravity models also provide distinct effects on the linear growth of structure, with f(R) inducing a modification to growth that varies with physical scale while DGP induces a modification that is scale independent.

The dark sector models provide two phenomenological examples of possible deviations from a cosmological constant. The time-varying dark energy model is the simplest way to approximate and parameterise dark energy evolution as predicted by, for example, quintessence scalar field theories. The specific interacting dark energy model we have chosen represents more complex dynamics and has been shown to be able to address the tension between CMB and structure growth measurements (Pourtsidou & Tram 2016; Carrilho et al. 2023).

On top of a modified gravitational or dark sector, we also considered massive neutrinos, the total mass of which Euclid is expected to measure with a high significance (see for example Archidiacono et al. 2017). The effects of a massive neutrino species becomes significant at the scales of interest and have demonstrated a notable degeneracy with modifications to gravity at low redshift (He 2013; Baldi et al. 2014; García-Farieta et al. 2019; Wright et al. 2019; Hagstotz et al. 2019; Moretti et al. 2023).

An additional choice we considered is the redshift space distortion effect (Kaiser 1987). We looked at two prominent models: the Taruya–Nishimishi–Saito (TNS) model (Taruya et al. 2010) and the EFTofLSS model (Perko et al. 2016; Chudaykin et al. 2020; D’Amico et al. 2020). We also examined choices for how to resum infrared modes and model tracer bias.

This paper is structured as follows. In Sect. 2 we introduce the four theories of gravity and dark energy that we consider. In Sect. 3 we cover all the various nonlinear redshift-space power spectrum models and approximations. In Sect. 4 we present the N-body simulations which we then use to quantify the validity of the various theoretical choices in Sect. 5. We summarised these findings and our conclusions in Sect. 6.

2. Perturbation theory beyond-ΛCDM

We consider a perturbed flat Friedmann–Lemaître–Robertson–Walker (FLRW) universe with a background metric given by ds² = −c²dt² + a²(t)δ_ijdxⁱdx^j, where a represents the scale factor. The Hubble rate, H, is defined as H := ȧ/a, where a dot denotes the derivative with respect to cosmic time t. We focus on scalar perturbations, and we adopt the Newtonian gauge where the perturbed FLRW metric can be written as

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-c^2[1+2\mathrm{\Phi }(\mathit{x},t)]\mathrm{d}{t}^2+a^2(t)[1-2\mathrm{\Psi }(\mathit{x},t)]{\delta }_{\mathit{ij}}\mathrm{d}{x}^i\mathrm{d}{x}^j,\end{array}$$(1)

where the gravitational potential, Φ, appears in the time-time component of the metric. To express the continuity and Euler equations, we use the rescaled CDM velocity divergence, θ := ∇ ⋅ v/(a H), and the CDM density contrast, δ := δρ/ρ, where ρ and δρ are respectively the energy density background and perturbations.

SPT assumes that density and velocity perturbations are small and can be expanded as

$$\begin{array}{c}\hfill \delta (\mathit{k},a)={\displaystyle \sum _{n=1}^{\infty }}{\delta }_n(\mathit{k},a),\theta (\mathit{k},a)={\displaystyle \sum _{n=1}^{\infty }}{\theta }_n(\mathit{k},a),\end{array}$$(2)

where δ_n, θ_n ∼ δ₁ⁿ, with δ₁ being the linear theory solution. Explicitly we can write the n^th order perturbations in terms of scale and time dependent kernels F_n and G_n, defined implicitly as

$$\begin{array}{cc}\hfill {\delta }_n(\mathit{k},a)& =\frac{1}{{(2\pi )}^{3(n-1)}}{\displaystyle \int {\mathrm{d}}^3}{\mathit{k}}_1\cdots {\mathrm{d}}^3{\mathit{k}}_n{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{1\cdots n})\hfill \\ \hfill & \times F_n({\mathit{k}}_1,\dots ,{\mathit{k}}_n,a){\delta }_{1,\mathrm{i}}({\mathit{k}}_1)\cdots {\delta }_{1,\mathrm{i}}({\mathit{k}}_n),\hfill \\ \hfill {\theta }_n(\mathit{k},a)& =\frac{1}{{(2\pi )}^{3(n-1)}}{\displaystyle \int {\mathrm{d}}^3}{\mathit{k}}_1\cdots {\mathrm{d}}^3{\mathit{k}}_n{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{1\cdots n})\hfill \\ \hfill & \times G_n({\mathit{k}}_1,\cdots ,{\mathit{k}}_n,a){\delta }_{1,\mathrm{i}}({\mathit{k}}_1)\cdots {\delta }_{1,\mathrm{i}}({\mathit{k}}_n),\hfill \end{array}$$(3)

where k_1⋯n = k₁ + … + k_n, δ_D denotes the Dirac delta and a subscript ‘i’ denotes a quantity computed at some early initial time. We define the linear growth rate of structure as

$$\begin{array}{c}\hfill f:=q\frac{\mathrm{d}ln{\delta }_1}{\mathrm{d}lna}=-G_1/F_1,\end{array}$$(4)

where in the second equality we have used the linear version of Eq. (5). This can now clearly be both time and scale-dependent, as in f(R) gravity (see Sect. 2.4).

To solve for the kernels, F_n and G_n, we can write down generic energy and momentum conservation equations for the CDM density and rescaled velocity divergence perturbations, as well as the Poisson equation which relates the Newtonian gravitational potential (Φ) to the density perturbation⁴

$$\begin{array}{cc}& a\delta \prime (\mathit{k},a)+\theta (\mathit{k},a)\hfill \\ \hfill & =-{\displaystyle \int \frac{{\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2}{{(2\pi )}^3}}{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{12})\alpha ({\mathit{k}}_1,{\mathit{k}}_2)\theta ({\mathit{k}}_1,a)\delta ({\mathit{k}}_2,a),\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill & a\theta \prime (\mathit{k},a)+(2+\frac{aH\prime (a)}{H(a)}+\frac{\mathrm{\Xi }(a)}{H(a)})\theta (\mathit{k},a)-{\left(\frac{ck}{aH(a)}\right)}^2\mathrm{\Phi }(\mathit{k},a)\hfill \\ \hfill & =-\frac{1}{2}{\displaystyle \int \frac{{\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2}{{(2\pi )}^3}}{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{12})\beta ({\mathit{k}}_1,{\mathit{k}}_2)\theta ({\mathit{k}}_1,a)\theta ({\mathit{k}}_2,a),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill & =-\frac{1}{2}{\displaystyle \int \frac{{\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2}{{(2\pi )}^3}}{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{12})\beta ({\mathit{k}}_1,{\mathit{k}}_2)\theta ({\mathit{k}}_1,a)\theta ({\mathit{k}}_2,a),\hfill \\ \hfill & -{\left(\frac{ck}{aH(a)}\right)}^2\mathrm{\Phi }(\mathit{k},a)=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}(a)}{2}\mu (k,a)\delta (\mathit{k},a)+S(\mathit{k},a),\hfill \end{array}$$(7)

where a prime denotes a derivative with respect to the scale factor, Ω_m is the total matter density fraction and Ξ encodes the effects of a possible additional drag force, coming from, say, an interaction within the dark sector (see Sect. 2.2). The kernels α(k₁, k₂) and β(k₁, k₂) are given by

$$\begin{array}{cc}\hfill \alpha ({\mathit{k}}_1,{\mathit{k}}_2)& =1+\frac{{\mathit{k}}_1·{\mathit{k}}_2}{|{\mathit{k}}_1{|}^2},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \beta ({\mathit{k}}_1,{\mathit{k}}_2)& =\frac{({\mathit{k}}_1·{\mathit{k}}_2){|{\mathit{k}}_1+{\mathit{k}}_2|}^2}{|{\mathit{k}}_1{|}^2{|{\mathit{k}}_2|}^2}.\hfill \end{array}$$(9)

The function μ(k, a) expresses any linear modification to GR in the Poisson equation (Bean & Tangmatitham 2010; Silvestri et al. 2013; Pogosian et al. 2010) and S(k, a) captures higher-order modifications, which include screening effects. Up to third order, we have

$$\begin{array}{cc}\hfill S(\mathit{k},a)& ={\displaystyle \int \frac{{\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2}{{(2\pi )}^3}}{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{12}){\gamma }_2({\mathit{k}}_1,{\mathit{k}}_2,a)\delta ({\mathit{k}}_1,a)\delta ({\mathit{k}}_2,a)\hfill \\ \hfill & +{\displaystyle \int \frac{{\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2{\mathrm{d}}^3{\mathit{k}}_3}{{(2\pi )}^6}}{\delta }_{\mathrm{D}}(\mathit{k}-{\mathit{k}}_{123}){\gamma }_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,a)\hfill \\ \hfill & \times \delta ({\mathit{k}}_1,a)\delta ({\mathit{k}}_2,a)\delta ({\mathit{k}}_3,a),\hfill \end{array}$$(10)

where the functions γ₂ and γ₃ are additional kernels coming from higher-order scalar field interactions in the Klein–Gordon equation (see Koyama et al. 2009, for example). For ΛCDM, we have μ(k, a) = 1 and Ξ = γ₂ = γ₃ = 0. We now look at the forms of these equations for various theories beyond-ΛCDM.

2.1. w₀w_aCDM

For a dark energy fluid the equation of state parameter is defined as w = p_DE/(ρ_DEc²), where p_DE and ρ_DE are the pressure and the background energy density of the fluid. In ΛCDM, w = −1. In general, w can be a function of time. One such (phenomenological) form was proposed by Chevalier–Polarski–Linder (CPL, Chevallier & Polarski 2001; Linder 2003):

$$\begin{array}{c}\hfill w(a)=w_0+(1-a)w_a,\end{array}$$(11)

where w₀ and w_a are constants. This form is essentially the first order Taylor expansion of w(a) about the present time, where we have chosen the present day scale factor to be a = 1. Taking only the first order Taylor expansion also implies we assume the equation of state is slowly evolving. Eq. (11) effectively changes the background expansion H but leaves all other functional modifications to ΛCDM at their ΛCDM values, that is μ(k, a) = 1 and Ξ = γ₂ = γ₃ = 0.

2.2. Dark Scattering model

Dark Scattering is an interacting dark energy model in which dark energy exchanges momentum with dark matter without energy exchange (Simpson 2010; Pourtsidou et al. 2013; Pourtsidou & Tram 2016; Baldi & Simpson 2015, 2017; Bose et al. 2018a). The modifications to Eq. (10) are given by μ(k, a) = 1 and γ₂ = γ₃ = 0. On the other hand, we have a coupling between dark matter and dark energy perturbations which appears in the Euler equations for each species (Simpson 2010; Baldi & Simpson 2015).

In Eqs. (5), (6) we made the assumption that the fluctuations of dark energy propagate at the speed of light, ensuring that they remain significantly smaller than the fluctuations of other species. This allows the rescaled velocity divergence of the dark energy fluid, θ_DE, to be effectively neglected in the Euler equation. Then, the interaction introduces only an additional drag force on the left-hand side of Eq. (6), given by

$$\begin{array}{c}\hfill \mathrm{\Xi }(a)=[1+w(a)]\xi c{\rho }_{\mathrm{DE}}(a),\end{array}$$(12)

where ξ := σ_D/(m_cc²). Here m_c is the dark matter particle mass and σ_D is the cross-section of the interaction. ξ is treated as a free parameter of the theory. It is positive (or zero) and proportional to the strength of the modification to ΛCDM, with ΛCDM being recovered for ξ → 0 or w → −1.

2.3. Dvali–Gabadadze–Porrati gravity

The Dvali–Gabadadze–Porrati (DGP) model has been introduced by Dvali et al. (2000) and describes a model where matter lives on a four-dimensional brane which is embedded in a five-dimensional Minkowski background. This model is representative of the class of models exhibiting the Vainshtein mechanism, which screens additional scalar field-sourced forces at small physical scales (Vainshtein 1972).

We consider the normal branch of DGP, which is free of ghost instabilities (see Luty et al. 2003; Gorbunov et al. 2006; Charmousis et al. 2006). In this model dark matter particles follow geodesics, implying that Ξ = 0 in Eq. (6). The relation between Φ and δ is described by Eq. (7), with (see for example Bose & Koyama 2016; Bose et al. 2018b)

$$\begin{array}{cc}\hfill \mu (k,a)& =1+\frac{1}{3\beta (a)},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\gamma }_2({\mathit{k}}_1,{\mathit{k}}_2,a)& =-{\left[\frac{H_0}{H(a)}\right]}^2\frac{1}{24{\beta }^3(a){\mathrm{\Omega }}_{\mathrm{rc}}}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m},0}}{a^3}\right)}^2(1-{\mu }_{1,2}^2),\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\gamma }_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,a)& ={\left[\frac{H_0}{H(a)}\right]}^2\frac{1}{144{\beta }^5(a){\mathrm{\Omega }}_{\mathrm{rc}}^2}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m},0}}{a^3}\right)}^3\hfill \\ \hfill & \times (1-{\mu }_{2,3}^2)(1-{\mu }_{1,23}^2),\hfill \end{array}$$(15)

where H₀ and Ω_m, 0 are the Hubble parameter and total matter density fraction today,

$$\begin{array}{c}\hfill \beta (a):=q1+\frac{H(a)}{H_0}\frac{1}{\sqrt{{\mathrm{\Omega }}_{\mathrm{rc}}}}(1+\frac{aH\prime (a)}{3H(a)}).\end{array}$$(16)

The third order kernel, γ₃, needs to be symmetrised and ${\mu }_{i,j}={\widehat{\mathit{k}}}_i·{\widehat{\mathit{k}}}_j$ is the cosine of the angle between k_i and k_j (recall k_ij = k_i + k_j). Here Ω_rc := c²/(4r_c²H₀²), where r_c represents the scale above which gravity deviates from GR, meaning r_c → ∞, or Ω_rc → 0, is the GR-limit of the theory (see Dvali et al. 2000). In this work we assume that the background expansion, H, is exactly the one of a flat ΛCDM model (Schmidt 2009).

2.4. f(R) gravity

In this work we consider the specific form for f(R) proposed by Hu & Sawicki (2007). This model is representative of the class of models exhibiting the Chameleon mechanism (Khoury & Weltman 2004), which screens in a phenomenologically distinct way to Vainshtein screening.

As with DGP, in this theory we have Ξ = 0, and have the following modifications to the Poisson equation (Koyama et al. 2009; Taruya et al. 2014; Bose & Koyama 2016)

$$\begin{array}{cc}\hfill \mu (k,a)& =1+{\left(\frac{k}{a}\right)}^2\frac{1}{3\mathrm{\Pi }(k,a)},\hfill \end{array}$$(17) $$\begin{array}{cc}\hfill {\gamma }_2({\mathit{k}}_1,{\mathit{k}}_2,a)& =-\frac{3}{16}{\left(\frac{k{H}_0}{aH(a)}\right)}^2{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m},0}}{a^3}\right)}^2\frac{{\mathrm{{\rm Y}}}^5(a)}{f_0^2{(3{\mathrm{\Omega }}_{\mathrm{m},0}-4)}^4}\hfill \\ \hfill & \times \frac{1}{\mathrm{\Pi }(k,a)\mathrm{\Pi }(k_1,a)\mathrm{\Pi }(k_2,a)},\hfill \end{array}$$(18)

and

$$\begin{array}{cc}& {\gamma }_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,a)\hfill \\ \hfill & =\frac{1}{32}{\left(\frac{k{H}_0}{aH(a)}\right)}^2{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m},0}}{a^3}\right)}^3\frac{1}{\mathrm{\Pi }(k,a)\mathrm{\Pi }(k_1,a)\mathrm{\Pi }(k_2,a)\mathrm{\Pi }(k_3,a)}\hfill \\ \hfill & \times \{-5\frac{{\mathrm{{\rm Y}}}^7(a)}{f_0^3{(3{\mathrm{\Omega }}_{\mathrm{m},0}-4)}^6}+\frac{9}{2}\frac{1}{\mathrm{\Pi }(k_{23},a)}{\left[\frac{{\mathrm{{\rm Y}}}^5(a)}{f_0^2{(3{\mathrm{\Omega }}_{\mathrm{m},0}-4)}^4}\right]}^2\},\hfill \end{array}$$(19)

where we have implicitly assumed the constraint coming from the 1-loop integral: k = k₁ + k₂ for γ₂ and k = k₁ + k₂ + k₃ for γ₃. Again, we note that γ₃ must be symmetrised. The functions Π and Υ are given by

$$\begin{array}{cc}\hfill \mathrm{\Pi }(k,a)& ={\left(\frac{k}{a}\right)}^2+\frac{{\mathrm{{\rm Y}}}^3(a)}{2f_0{(3{\mathrm{\Omega }}_{\mathrm{m},0}-4)}^2},\hfill \end{array}$$(20) $$\begin{array}{cc}\hfill \mathrm{{\rm Y}}(a)& =\frac{{\mathrm{\Omega }}_{\mathrm{m},0}+4a^3(1-{\mathrm{\Omega }}_{\mathrm{m},0})}{a^3}.\hfill \end{array}$$(21)

Here we set f₀ = c²|f_R0|/H₀², where f_R0 is the value of f_R = df(R)/dR today, which gauges the strength of the modification to GR. The GR-limit is given by f_R0 → 0. As with the DGP case, we also assume a ΛCDM background expansion, H.

3. Nonlinear modelling

When mapping the underlying theories described in the previous section to the observable quantity of interest, there are several choices one can make. This range of choices expands significantly when we wish to model scales where nonlinearities become important. In this section, we present a number of such choices which we test to varying degrees against N-body simulations in Sect. 5.

3.1. The quasi-nonlinear power spectrum

Our goal is to calculate the power spectrum of the dark matter density and velocity fluctuations. This can be defined as

$$\begin{array}{c}\hfill {(2\pi )}^3{\delta }_{\mathrm{D}}(\mathit{k}+\mathit{k}\prime )P^{\mathit{ab}}(k,a):=⟨{\phi }^a(\mathit{k},a){\phi }^b(\mathit{k}\prime ,a)⟩,\end{array}$$(22)

where we have introduced the duplet φ^a := (δ, θ), and ⟨ ⋅ ⟩ denotes an ensemble average. By expanding δ and θ using Eq. (2), we can perturbatively express the power spectrum as a leading-order term plus a next-to-leading-order (nlo) term

$$\begin{array}{c}\hfill P^{\mathit{ab}}(k,a)=P_{11}^{\mathit{ab}}(k,a)+P_{\mathrm{nlo}}^{\mathit{ab}}(k,a)+\dots ,\end{array}$$(23)

where the leading term, $P_{11}^{\mathit{ab}}$, is the ‘tree-level’ or linear power spectrum. The nlo or ‘1-loop’ term can be expressed as the sum of three quantities,

$$\begin{array}{c}\hfill P_{\mathrm{nlo}}^{\mathit{ab}}(k,a)=P_{22}^{\mathit{ab}}(k,a)+P_{13}^{\mathit{ab}}(k,a)+P_{31}^{\mathit{ab}}(k,a),\end{array}$$(24)

where we have used the notation

$$\begin{array}{c}\hfill {(2\pi )}^3{\delta }_{\mathrm{D}}(\mathit{k}+\mathit{k}\prime )P_{\mathit{ij}}^{\mathit{ab}}(k,a)=⟨{\phi }_i^a(\mathit{k},a){\phi }_j^b(\mathit{k}\prime ,a)⟩.\end{array}$$(25)

The ellipses in Eq. (23) represent terms of higher-order, including contributions from 2-loop order and beyond. We define the ‘1-loop’ power spectrum as

$$\begin{array}{c}\hfill P_{1-\mathrm{loop}}^{\mathit{ab}}:=P_{11}^{\mathit{ab}}(k,a)+P_{\mathrm{nlo}}^{\mathit{ab}}(k,a).\end{array}$$(26)

3.2. Redshift space

So far we have worked in ‘real’ space coordinates, r, or its Fourier equivalent. But astronomical observations are composed of angular positions and redshifts, which contains all depth information. This means we must convert our theoretical predictions to redshift space coordinates, s. This involves peculiar velocity information which combines with the Hubble velocity to give the total measured redshift. In terms of positional coordinates, this is expressed as

$$\begin{array}{c}\hfill \mathit{s}=\mathit{r}+\frac{v_z(\mathit{r},a)}{aH(a)}\widehat{\mathit{z}},\end{array}$$(27)

where we have taken the line-of-sight direction to be along the real space z-axis, and v_z is the projection of the peculiar velocity along that axis. Conservation of mass then implies the translation

$$\begin{array}{c}\hfill {\delta }^s(\mathit{s},a)={\left|\frac{\partial \mathit{s}}{\partial \mathit{r}}\right|}^{-1}[1+\delta (\mathit{r},a)]-1,\end{array}$$(28)

where a superscript ‘s’ denotes a redshift-space quantity. If we Fourier transform the redshift-space density field and take the ensemble average with itself at two points, we can arrive at the following expression for the redshift-space power spectrum at 1-loop order (Matsubara 2008; Heavens et al. 1998; Taruya et al. 2010),

$$\begin{array}{cc}\hfill P_{1-\mathrm{loop}}^s(k,\mu )& =P_{1-\mathrm{loop}}^{\delta \delta }(k)+2{\mu }^2{P}_{1-\mathrm{loop}}^{\delta \theta }(k)+{\mu }^4{P}_{1-\mathrm{loop}}^{\theta \theta }(k)\hfill \\ \hfill & +A(k,\mu )+B(k,\mu )+C(k,\mu )\hfill \\ \hfill & -k^2{\mu }^2{\tilde{\sigma }}_{\mathrm{v}}^2{[F_1(k)+G_1(k){\mu }^2]}^2{P}_{11,\mathrm{i}}^{\delta \delta }(k),\hfill \end{array}$$(29)

where we have dropped the explicit time dependence of all functions for compactness, except for $P_{11,\mathrm{i}}^{\delta \delta }$ which is computed at some fixed early time. Here μ is the angle between the line of sight direction $\widehat{\mathit{z}}$ and the wave mode k⁵. The linear theory estimate for the velocity dispersion, ${\stackrel{\sim }{\sigma }}_{\mathrm{v}}$, is given by

$$\begin{array}{c}\hfill {\tilde{\sigma }}_{\mathrm{v}}^2=\frac{1}{6{\pi }^2}{\displaystyle \int \mathrm{d}}q{G}_1^2(q)P_{11,\mathrm{i}}^{\delta \delta }(q).\end{array}$$(30)

We refer the reader to Bose & Koyama (2016), Taruya et al. (2010) for expressions for A(k, μ), B(k, μ), and C(k, μ). These terms come from higher-order interactions between the density and velocity perturbations. We give their explicit forms in terms of the kernels defined in Eq. (3) in Appendix A.

Equation (29) performs poorly when modelling nonlinear redshift-space effects. In particular, it does not capture the fingers-of-god (FoG) effect, caused by the large positional distortions of measured objects towards the centres of potential wells, coming from their high peculiar velocities. Higher-order perturbation theory (higher than third order), has been shown to have poor convergence at the level of the power spectrum (see, for example, Carlson et al. 2009), and is also very computationally expensive. To move further into the nonlinear regime, a number of proposals have been made in the past decade which modify Eq. (29) by introducing new degrees of freedom, or nuisance parameters, quantifying our uncertainty on nonlinear effects.

We consider two prominent such proposals, both of which have been applied in the biggest galaxy survey to date, the BOSS survey (Beutler et al. 2017; Ivanov et al. 2020; D’Amico et al. 2020; Carrilho et al. 2023). Both of these can also be computed with a number of fast public codes, discussed next.

3.2.1. Codes

In this work we use a number of codes. In particular, we adopt:

• The publicly available MG-Copter code (Bose & Koyama 2016), recently absorbed by the ReACT code⁶ (Bose et al. 2020b, 2021, 2023) which solves the first, second and third order continuity and Euler equations simultaneously and numerically for each external momentum;

• The publicly available PyBird code⁷ (D’Amico et al. 2021b), which uses a Fast Fourier Transform (FFT) decomposition and has been used in recent BOSS analyses (for instance see Zhang et al. 2022; D’Amico et al. 2021b; Piga et al. 2023, for ΛCDM, wCDM, and nDGP respectively);

• The PBJ code (Moretti et al. 2023; Oddo et al. 2020, 2021; Rizzo et al. 2023), which implements the model of Ivanov et al. (2020) taking advantage of the FAST-PT algorithm (McEwen et al. 2016; Fang et al. 2017), and has also been used in recent BOSS analyses for Dark Scattering (Carrilho et al. 2023). This code is planned to be made public soon.

We summarise the codes in Table 1 along with relevant implementations. MG-Copter is by far the slowest, needing to solve eight sets of Eqs. (5), (6) for each external mode k to find the perturbative kernels up to third order. More quantitatively, to compute the first two multipoles of Eq. (29) in ΛCDM at fifty values of k in the range [0.01, 0.3] h Mpc⁻¹, it takes roughly 45 s on a MacBook Pro 2018 model. By comparison, PBJ and PyBird only need to compute a set of products on a fixed grid by using the FFT method, which allows them to produce the same result in roughly 20 ms. Of course, this assumes a specific scale dependence of the perturbative kernels, which may not be known or even analytically available. This is not true for MG-Copter which only requires the specification of H, μ, γ₂, γ₃ and Ξ.

We have performed validation of these codes both in real (see Eq. 26) and redshift space (see Eq. 29) to sub-percent precision for ΛCDM, w₀w_aCDM, nDGP, and Dark Scattering (see Appendix C for some of these comparisons). Additionally, MG-Copter, PBJ and PyBird have been validated independently against large N-body simulations in a number of papers: see Bose & Koyama (2016), Bose et al. (2017, 2018a) for MG-Copter, Oddo et al. (2021), Carrilho et al. (2021), Tsedrik et al. (2023) for PBJ and Nishimichi et al. (2020), D’Amico et al. (2020) for PyBird.

Within the Euclid collaboration, both PBJ and PyBird are being used to fit measurements from the Flagship simulation (Potter et al. 2017), with results to be presented in a series of papers. Finally, PBJ has been ported to the Euclid likelihood code, named CLOE, and will be used to perform the official analysis of the spectroscopic sample for ΛCDM.

3.2.2. Effective Field Theory of Large-Scale Structure

The effective field theory of large-scale structure (EFTofLSS) prescription for the redshift-space matter power spectrum can be written as (D’Amico et al. 2020; Ivanov et al. 2020; Chudaykin et al. 2020)

$$\begin{array}{cc}\hfill P_{\mathrm{EFT}}^s(k,\mu )& =P_{1-\mathrm{loop}}^s(k,\mu )+{\tilde{P}}_{\mathrm{ctr}}^s(k,\mu ),\hfill \end{array}$$(31)

where the second term, ${\stackrel{\sim }{P}}_{\mathrm{ctr}}^s(k,\mu )$, contains ‘counterterms’ which are used to model nonlinear effects, and contains a number of free constants that must be fit to the data.

In PBJ the EFTofLSS counterterms are given by

$$\begin{array}{c}\hfill {\tilde{P}}_{\mathrm{ctr}}^s(k,\mu )=-2k^2{P}_{11}^{\delta \delta }(k)({\tilde{c}}_0+{\tilde{c}}_{2}f{\mu }^2+{\tilde{c}}_4{f}^2{\mu }^4).\end{array}$$(32)

This model has therefore a total of three nuisance parameters at the dark matter level, $\{{\stackrel{\sim }{c}}_0,{\stackrel{\sim }{c}}_2,{\stackrel{\sim }{c}}_4\}$, with dimensions h⁻² Mpc².

In PyBird a slightly different, dimensionless, basis is used for the counterterms (see D’Amico et al. 2021b; Piga et al. 2023, for example). One has

$$\begin{array}{c}\hfill P_{\mathrm{ctr}}^s(k,\mu )=2P_{11}^{\delta \delta }(k){\left(\frac{k}{k_M}\right)}^2(b_1+f{\mu }^2)(c_{\mathrm{ctr}}+c_{r,1}{\mu }^2+c_{r,2}{\mu }^4),\end{array}$$(33)

where k_M = 0.7 h⁻¹ Mpc⁻¹ parametrises the inverse spatial extension scale of galaxies. The total number of nuisance parameters is therefore the same: {c_ctr, c_r, 1, c_r, 2}.

When we only consider the monopole and quadrupole of Eq. (31), as is done in the majority of this paper, we can set ${\stackrel{\sim }{c}}_4=c_{r,2}=0$. This yields the following mapping between the two different counterterm bases

$$\begin{array}{c}\hfill c_{\mathrm{ctr}}=-\frac{k_M^2(35b_1{\tilde{c}}_0+30{\tilde{c}}_{0}f+3{\tilde{c}}_2{f}^2)}{35b_1^2+30b_{1}f+3f^2},\end{array}$$(34) $$\begin{array}{c}\hfill c_{r,1}=-\frac{35k_M^2(b_1{\tilde{c}}_2-{\tilde{c}}_0)}{35b_1^2+30b_{1}f+3f^2}.\end{array}$$(35)

This mapping was obtained after imposing that the monopole and the quadrupole of the two different expressions used for the counterterm, Eqs. (32), (33), have to be the same respectively.

We remind the reader that in theories such as f(R), the growth rate, f = −G₁(k, a)/F₁(k, a), is both time and scale-dependent, and the linear power spectrum has additional scale dependencies coming from F₁(k, a): $P_{11}^{\delta \delta }=F_1^2(k,a)P_{11,\mathrm{i}}^{\delta \delta }$.

3.2.3. Taruya–Nishimichi–Saito

The TNS prescription is given by (Taruya et al. 2010)

$$\begin{array}{cc}\hfill P_{\mathrm{TNS}}(k,\mu )& =D_{\mathrm{FoG}}({\mu }^2{k}^2{\sigma }_{\mathrm{v}}^2)[P_{1-\mathrm{loop}}^{\delta \delta }(k)+2{\mu }^2{P}_{1-\mathrm{loop}}^{\delta \theta }(k)\hfill \\ \hfill & +{\mu }^4{P}_{1-\mathrm{loop}}^{\theta \theta }(k)+A(k,\mu )+B(k,\mu )],\hfill \end{array}$$(36)

where A and B terms are as in Eq. (29), and the C term has been omitted as it can be largely absorbed in D_FoG (Taruya et al. 2010). The damping factor D_FoG encodes a phenomenological description of the FoG damping. This is a function of the velocity dispersion, σ_v, which is taken to be a free nuisance parameter. Here we assume this function takes a Lorentzian form so that

$$\begin{array}{c}\hfill D_{\mathrm{FoG}}(k^2{\mu }^2{\sigma }_{\mathrm{v}}^2):=\frac{1}{1+(k^2{\mu }^2{\sigma }_{\mathrm{v}}^2)/2}.\end{array}$$(37)

This form has been shown to be good at describing nonlinear effects (see for example Bose et al. 2020a). This model then has a single nuisance parameter at the dark matter level, {σ_v}, with units h⁻¹ Mpc, which represents an effective velocity dispersion, not to be confused with the linear estimate given in Eq. (30).

3.3. Tracer bias

The next step to connect theory with the spectroscopic observable is to translate the dark matter spectrum to that of galaxies, which are biased tracers of the CDM distribution. As a proxy, one can consider the CDM halo power spectrum. This can be modelled incorporating an appropriate bias model.

3.3.1. Eulerian bias expansion

A general Eulerian perturbative bias expansion can be used to relate the halo density field to the dark matter density field, with the (scale-independent) coefficients of each term being treated as free parameters to be fit to observations (McDonald & Roy 2009; Chan et al. 2012; Assassi et al. 2014; Senatore 2015; Mirbabayi et al. 2015; Desjacques et al. 2018; Fujita et al. 2020). The terms in the bias expansion that are relevant for the calculation of the 1-loop power spectrum of halos are

$$\begin{array}{c}\hfill {\delta }_{\mathrm{halo}}=b_1\delta +\frac{b_2}{2}{\delta }^2+b_{{\mathcal{G}}_2}{\mathcal{G}}_2\left(\tilde{\mathrm{\Phi }}|\mathbf{x}\right)+b_{{\mathrm{\Gamma }}_3}{\mathrm{\Gamma }}_3+b_{{\mathrm{\nabla }}^2\delta }{\mathrm{\nabla }}^2\delta +\mathrm{noise},\end{array}$$(38)

where b₁ and b₂ are the linear and quadratic bias parameters, ∇²δ is a higher-derivative operator and ‘noise’ denotes the stochastic contributions uncorrelated with δ and with zero means, encoding shot-noise. The ∇²δ term is completely degenerate with the effects of the first counterterm (see for example Perko et al. 2016, and Eq. (32)) and so we do not consider it further here, as was also done in previous EFTofLSS-based clustering data analyses (Ivanov et al. 2020; D’Amico et al. 2020). 𝒢₂ and Γ₃ are non-local operators which take into account the large-scale tidal fields at leading and next-to-leading order. In configuration space, they are defined as:

$$\begin{array}{c}\hfill {\mathcal{G}}_2\left(\tilde{\mathrm{\Phi }}|\mathbf{x}\right):={\left[{\mathrm{\nabla }}_{\mathit{ij}}\tilde{\mathrm{\Phi }}(\mathbf{x})\right]}^2-{\left[{\mathrm{\nabla }}^2\tilde{\mathrm{\Phi }}(\mathbf{x})\right]}^2,\end{array}$$(39)

and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_3(\mathbf{x}):={\mathcal{G}}_2\left(\tilde{\mathrm{\Phi }}|\mathbf{x}\right)-{\mathcal{G}}_2\left({\tilde{\mathrm{\Phi }}}_{\mathrm{v}}|\mathbf{x}\right),\end{array}$$(40)

where the potentials $\stackrel{\sim }{\mathrm{\Phi }}$ and ${\stackrel{\sim }{\mathrm{\Phi }}}_{\mathrm{v}}$ are defined as $\stackrel{\sim }{\mathrm{\Phi }}={\mathrm{\nabla }}^{-2}\delta$, ${\stackrel{\sim }{\mathrm{\Phi }}}_{\mathrm{v}}=-{\mathrm{\nabla }}^{-2}(\theta /f)$. We have also defined ∇_ij := ∇_i∇_j and ∇² := ∇_i∇_i.

Moving to redshift space, we obtain the following halo redshift-space power spectrum (see Senatore & Zaldarriaga 2014, for example)

$$\begin{array}{cc}\hfill P_g^s(k,\mu )& =D_{\mathrm{FoG}}(k^2{\mu }^2{\sigma }_{\mathrm{v}}^2)[P_{\mathrm{g},1-\mathrm{loop}}^s(k,\mu )+P_{\mathrm{g},\mathrm{noise}}^s(k,\mu )]\hfill \\ \hfill & +P_{\mathrm{ctr}}^s(k,\mu ).\hfill \end{array}$$(41)

The term in square brackets is what we get from taking the ensemble average of the redshift space halo density perturbation in Eq. (38).

We have also included the nonlinear corrections coming from the EFTofLSS and TNS. In the EFTofLSS, D_FoG(k²μ²σ_v²) = 1 and the counterterms $P_{\mathrm{ctr}}^s(k,\mu )$ are given by Eqs. (32), (33) for PBJ and PyBird respectively. In TNS, $P_{\mathrm{ctr}}^s(k,\mu )=0$ while D_FoG(k²μ²σ_v²) is given by Eq. (37). Further, in the TNS case we do not include the last two terms of Eq. (29), in accordance with Taruya et al. (2010).

The power spectrum of the stochastic part of Eq. (38), and as appearing in PBJ, is given by

$$\begin{array}{cc}\hfill P_{\mathrm{g},noise}^s(k,\mu )& =N+\frac{1}{\overline{n}}({ϵ}_0{k}^2+{ϵ}_2{\mu }^2{k}^2),\hfill \end{array}$$(42)

where $\overline{n}$ is the galaxy number density in units of h³ Mpc⁻³ and N, ϵ₀, ϵ₂ are constants with dimensions h⁻³ Mpc³, h⁻² Mpc² and h⁻² Mpc², respectively. Pybird uses the same modelling, with dimensionless parameters $c_{ϵ,0}=N\overline{n}$, c_ϵ, 1 = ϵ₀k_M² and c_ϵ, 2 = ϵ₂k_M²/f to give

$$\begin{array}{cc}\hfill P_{\mathrm{g},noise}^s(k,\mu )& =\frac{1}{\overline{n}}(c_{ϵ,0}+c_{ϵ,1}\frac{k^2}{k_M^2}+c_{ϵ,2}f{\mu }^2\frac{k^2}{k_M^2}).\hfill \end{array}$$(43)

In this work we set ϵ₀ = ϵ₂ = 0 and only consider the constant stochastic term, N, as done in the flagship BOSS analysis (Beutler et al. 2017) and in Ivanov et al. (2020). We have checked that the addition of ϵ₀ does not significantly improve the fits to the simulations. ϵ₂ however was found to improve the modelling significantly and has been included in some EFTofLSS-BOSS analyses (for example D’Amico et al. 2020). We have indeed checked that the goodness-of-fit to selected simulations does improve by including ϵ₂. Despite this, we choose to omit it in order to keep the complexity of the model minimal, which should not affect relative differences between ΛCDM and beyond-ΛCDM scenarios.

In the notation of Eq. (29), the 1-loop part can be expanded as

$$\begin{array}{cc}\hfill P_{\mathrm{g},1-\mathrm{loop}}^s(k,\mu )& =b_1^2{P}_{1-\mathrm{loop}}^{\delta \delta }(k)+2b_1{\mu }^2{P}_{1-\mathrm{loop}}^{\delta \theta }(k)+{\mu }^4{P}_{1-\mathrm{loop}}^{\theta \theta }(k)\hfill \\ \hfill & +A(k,\mu ,b_1)+B(k,\mu ,b_1)+C(k,\mu ,b_1)\hfill \\ \hfill & -k^2{\mu }^2{\tilde{\sigma }}_{\mathrm{v}}^2{[b_1{F}_1(k)+G_1(k){\mu }^2]}^2{P}_{11,\mathrm{i}}^{\delta \delta }(k)\hfill \\ \hfill & +B_{\mathrm{g},22}^s(k,\mu )+B_{\mathrm{g},13}^s(k,\mu ),\hfill \end{array}$$(44)

where we now have included the linear bias contributions explicitly for the 1-loop spectra and implicitly for the A, B, and C correction terms (see Appendix A).

The last two terms of Eq. (44) include the contributions from the higher-order bias. In this work we treat these terms under the Einstein–de Sitter approximation outlined in Sect. 3.6. We make this approximation explicit in the relevant expressions to follow. Our reasoning is that any additional scale dependencies coming from beyond-ΛCDM physics can be absorbed by the higher-order bias coefficients, and indeed a full derivation of these terms for general cosmologies and theories of gravity has proven not to be necessary under the considerations of this paper (see Sect. 5).

Under this assumption, the higher-order bias terms can be expressed as (see Appendix A of Ivanov et al. 2020, for example)

$$\begin{array}{cc}\hfill B_{\mathrm{g},22}^s(k,\mu )& =2{\displaystyle \int {\mathrm{d}}^3}q{P}_{11}^{\delta \delta }(q)P_{11}^{\delta \delta }(|\mathbf{k}-\mathbf{q}|)\hfill \\ \hfill & \times Z_{b,2}(\mathbf{q},\mathbf{k}-\mathbf{q})[\frac{b_2}{2}+b_{{\mathcal{G}}_2}{S}^b(\mathit{q},\mathit{k}-\mathit{q})],\hfill \end{array}$$(45) $$\begin{array}{cc}\hfill B_{g,13}^s(k,\mu )& =6Z_{b,1}(\mu )P_{11}^{\delta \delta }(k){\displaystyle \int {\mathrm{d}}^3}q{Z}_{b,3}^{\mathrm{sym}}(\mathbf{q},-\mathbf{q},\mathbf{k})P_{11}^{\delta \delta }(q),\hfill \end{array}$$(46)

where S^b(k₁, k₂) = (k₁ ⋅ k₂)²/k₁²k₂² − 1. We note here that we have included the linear growth implicitly in the linear spectra and have assumed it is calculated at the external wave mode, for example $P_{11}^{\delta \delta }(q)=F_1^2(k,a)P_{11,\mathrm{i}}^{\delta \delta }(q)$. The Z_b, n kernels can be written as (see Scoccimarro et al. 1999; Bernardeau et al. 2002, see)

$$\begin{array}{cc}\hfill Z_{b,1}(\mu )& =b_1+f{\mu }^2,\hfill \end{array}$$(47) $$\begin{array}{cc}\hfill Z_{b,2}({\mathbf{k}}_1,{\mathbf{k}}_2)& =2b_1{F}_2^{\mathrm{EdS}}({\mathbf{k}}_1,{\mathbf{k}}_2)+2f{\mu }^2{G}_2^{\mathrm{EdS}}({\mathbf{k}}_1,{\mathbf{k}}_2)\hfill \\ \hfill & +f\mu k[\frac{{\mu }_1}{k_1}(b_1+f{\mu }_2^2)+\frac{{\mu }_2}{k_2}(b_1+f{\mu }_1^2)]\hfill \\ \hfill & +\frac{b_2}{2}+b_{{\mathcal{G}}_2}{S}^b({\mathbf{k}}_1,{\mathbf{k}}_2),\hfill \end{array}$$(48) $$\begin{array}{cc}\hfill Z_{b,3}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)& =\frac{b_2}{2}f\mu k\frac{{\mu }_1}{k_1}+b_{{\mathcal{G}}_2}f\mu k\frac{{\mu }_1}{k_1}{S}^b({\mathbf{k}}_2,{\mathbf{k}}_3)\hfill \\ \hfill & +b_2{F}_2^{\mathrm{EdS}}({\mathbf{k}}_1,{\mathbf{k}}_2)+2b_{{\mathcal{G}}_2}{S}^b({\mathbf{k}}_1,{\mathbf{k}}_{23})F_2^{\mathrm{EdS}}({\mathbf{k}}_2,{\mathbf{k}}_3)\hfill \\ \hfill & +2b_{{\mathrm{\Gamma }}_3}{S}^b({\mathbf{k}}_1,{\mathbf{k}}_{23})[F_2^{\mathrm{EdS}}({\mathbf{k}}_2,{\mathbf{k}}_3)-G_2^{\mathrm{EdS}}({\mathbf{k}}_2,{\mathbf{k}}_3)],\hfill \end{array}$$(49)

where ‘EdS’ means the kernel is an Einstein–de Sitter expression. We note that f = −G₁(k, a)/F₁(k, a) is calculated exactly for the theory of gravity or dark energy under consideration at the external wave mode, k. Further, the Z_b, 3 kernel needs to be symmetrised to get $Z_{b,3}^{\mathrm{sym}}$ and variables k and μ are the amplitude of the sum of the momenta k = k₁ + k₂ + k₃ and cosine of angle between $\widehat{k}$ and the line of sight $\widehat{z}$.

Now we report here the mapping among the different bias bases used in the analyses of this work for both PBJ and Pybird codes. The bias basis used by PBJ is presented in Eq. (38), and uses [b₁, b₂, b_𝒢2, b_Γ3], while Pybird adopts the basis presented in D’Amico et al. (2020), that we indicate with $[{\widehat{b}}_1,{\widehat{b}}_2,{\widehat{b}}_3,{\widehat{b}}_4]$. There is a one-to-one correspondence between the two bases:

$$\begin{array}{cc}\hfill b_1& ={\widehat{b}}_1,\hfill \end{array}$$(50) $$\begin{array}{cc}\hfill b_2& =-2({\widehat{b}}_1-{\widehat{b}}_2-{\widehat{b}}_4),\hfill \end{array}$$(51) $$\begin{array}{cc}\hfill b_{{\mathcal{G}}_2}& =\frac{2}{7}[{\widehat{b}}_2+\frac{7}{4}(a_{\gamma }-2){\widehat{b}}_1],\hfill \end{array}$$(52) $$\begin{array}{cc}\hfill b_{{\mathrm{\Gamma }}_3}& =(\frac{a_{\gamma }^2}{2}-\frac{a_{\gamma a}}{2}+1){\widehat{b}}_1+\frac{2}{7}{a}_{\gamma }{\widehat{b}}_2-\frac{2}{21}{\widehat{b}}_3.\hfill \end{array}$$(53)

Here a_γ and a_γa are the bootstrap time-dependent functions introduced in D’Amico et al. (2021a) and are usually fixed to their EdS values a_γ^EdS = 10/7 and $a_{\gamma a}^{\mathrm{EdS}}=3/7$. We fix b_Γ3 = 0 because it is degenerate with b_𝒢2, as well as b_∇²δ = 0 as it is degenerate with a counterterm. We vary all the other parameters.

Further, using the TNS model, we wish to test the Local Lagrangian relationship (Sheth et al. 2013; Baldauf et al. 2012; Saito et al. 2014) for beyond-ΛCDM cases, which is known to hold well for ΛCDM. Again, we adopt this choice for the TNS model only, as was applied in the BOSS survey (Beutler et al. 2017). This leaves us with only the following bias parameters in the TNS case {b₁, b₂, N}, which are those used in Beutler et al. (2017) and can be related to the basis used in McDonald & Roy (2009) through Saito et al. (2014). To map these to the basis used in Assassi et al. (2014), which is used in our EFTofLSS models, see Desjacques et al. (2018).

3.3.2. Scale-dependent Q-bias

For the gravity and dark energy models that produce a scale-independent modification to the linear growth factor F₁(a), one does not expect any significant additional scale-dependencies entering the bias expansion (see, for example, Valogiannis et al. 2020; Valogiannis & Bean 2019). However, this is not the case in f(R) which yields a scale-dependent growth factor F₁(k, a) and growth rate f(k, a). To address possible scale-dependent modifications to the bias parameters, we also tested the phenomenological ‘Q-bias’ prescription for the linear bias of Song et al. (2015),

$$\begin{array}{c}\hfill b_1(k)=b_0\frac{1+A_2{k}^2}{1+A_{1}k},\end{array}$$(54)

which was shown to match N-body halo measurements well. This introduces two new free parameters {A₁, A₂} with dimensions h⁻¹ Mpc and h⁻² Mpc² respectively, characterising possible scale-dependencies of the linear bias. To keep the degrees of freedom comparable to the Eulerian bias expansion, when adopting this model, we set all higher-order bias terms in Eq. (44) to zero, and simply replace b₁ with Eq. (54).

3.4. Infrared-Resummation

A further complication comes from the baryonic acoustic oscillations whose imprint on the power spectrum should be treated properly to avoid percent-level oscillatory inaccuracies (Carrasco et al. 2014; Blas et al. 2016). To treat this insufficient damping of the baryon acoustic oscillations, various resummation methods have been proposed. This is particularly relevant for the EFTofLSS and we do not adopt an infrared-resummation scheme in the TNS case. We describe the method implemented by each EFTofLSS code below.

3.4.1. Wiggle-no-Wiggle decomposition

PBJ and MG-Copter adopt a “Wiggle-no-Wiggle” (WnW) decomposition approach, both under the Einstein–de Sitter approximation (see Sect. 3.6). This resummation method is based on a splitting of the linear power spectrum into a “wiggle” part (containing the baryon acoustic oscillations features) and a “no-wiggle”, broadband part (Baldauf et al. 2015; Vlah et al. 2015; de la Bella et al. 2017):

$$\begin{array}{c}\hfill P_{11}^{\delta \delta }=P_{11}^{\mathrm{w}}+P_{11}^{\mathrm{nw}}.\end{array}$$(55)

The loop integrals of Eq. (44) are then computed separately on the two components, and resummed with a damping applied to the wiggle component. Specifically, PBJ computes the smooth component by convolving the Eisenstein and Hu fit for the broadband linear spectrum (Eisenstein & Hu 1997) with a Gaussian filter:

$$\begin{array}{c}\hfill P_{11}^{\mathrm{nw}}(k)=P_{\mathrm{EH}}(k)[\mathcal{F}\ast R](k),\end{array}$$(56)

where P_EH is the Eisenstein and Hu prescription for the linear matter power spectrum, $R=P_{11}^{\delta \delta }(k)/P_{\mathrm{EH}}(k)$ and ℱ is a 1D Gaussian filter, so that the convolution reads

$$\begin{array}{c}\hfill [\mathcal{F}\ast R](k)=\frac{1}{\sqrt{2\pi }\lambda }{\displaystyle {\int }_{q_{\min }}^{q_{\max }}}\frac{\mathrm{d}q}{q}R(q)exp\{-\frac{{(lnk/q)}^2}{2{\lambda }^2}\},\end{array}$$(57)

with λ being the dimensionless width of the Gaussian filter.

We then compute the loop integrals using the linear ($P_{1-\mathrm{loop}}^s$) and smooth power spectra ($P_{1-\mathrm{loop}}^{s,\mathrm{nw}}$), with their difference being the loop corrections to the wiggle component ($P_{1-\mathrm{loop}}^{s,\mathrm{w}}=P_{1-\mathrm{loop}}^s-P_{1-\mathrm{loop}}^{s,\mathrm{nw}}$). Finally, we construct the infrared-resummed 1-loop power spectrum as

$$\begin{array}{cc}\hfill P_{\mathrm{NL}}^{s,\mathrm{IR}-\mathrm{resum}}(k)& ={(b_1+f{\mu }^2)}^2\{P_{11}^{\mathrm{nw}}(k)\hfill \\ \hfill & +{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2(\mathit{k},{\mathit{\ell }}_{\mathrm{osc}})}{P}_{11}^{\mathrm{w}}(k)[1+k^2{\mathrm{\Sigma }}^2(\mathit{k},{\mathit{\ell }}_{\mathrm{osc}})]\}\hfill \\ \hfill & +P_{1-\mathrm{loop}}^{s,\mathrm{nw}}(k)+{\mathrm{e}}^{-k^2{\mathrm{\Sigma }}^2(\mathit{k},{\mathit{\ell }}_{\mathrm{osc}})}{P}_{1-\mathrm{loop}}^{s,\mathrm{w}}(k),\hfill \end{array}$$(58)

where Σ²(k, ℓ_osc) is the RSD damping function computed as

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}^2(\mathit{k},\mathit{\ell })& =\frac{1}{2}\{{\mathrm{\Xi }}_0(\ell )(1+2f{\mu }^2+f^2{\mu }^2)\hfill \\ \hfill & +{\mathrm{\Xi }}_2(\ell )[{(\widehat{\mathit{k}}·\widehat{\mathit{\ell }})}^2+2f\mu {\mu }_{\ell }(\widehat{\mathit{k}}·\widehat{\mathit{\ell }})+f^2{\mu }^2{\mu }_{\ell }^2]\},\hfill \end{array}$$(59)

and we have defined

$$\begin{array}{cc}& {\mathrm{\Xi }}_0(\ell )=\frac{2}{3}{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}p}{2{\pi }^2}{\mathrm{e}}^{-p^2/{\mathrm{\Lambda }}_{\mathrm{IR}}^2}{P}_{11}^{\delta \delta }(p)[1-j_0(p\ell )-j_2(p\ell )],\hfill \end{array}$$(60) $$\begin{array}{cc}& {\mathrm{\Xi }}_2(\ell )=2{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}p}{2{\pi }^2}{\mathrm{e}}^{-p^2/{\mathrm{\Lambda }}_{\mathrm{I}R}^2}{P}_{11}^{\delta \delta }(p)j_2(p\ell ).\hfill \end{array}$$(61)

Here j_n are the spherical Bessel functions of order n and μ_ℓ is the angle between ℓ and the line of sight. Λ_IR is a cutoff scale. In PBJ we do not include the exponential cutoff but rather truncate the upper bound of the integral to k_s. PyBird also uses these functions (see Sect. 3.4.2), but instead of a cutoff, it sets ${\mathrm{\Lambda }}_{\mathrm{IR}}=0.2h{\mathrm{Mpc}}^{-1}$. ℓ_osc = |ℓ_osc| is set to be the baryon acoustic oscillations scale.

For PBJ, we adopt the following choice of parameters: q_min = k e^−4λ, q_max = k e^4λ, k_s = 0.2 h Mpc⁻¹, λ = 0.25 and ℓ_osc = 102.707 h⁻¹ Mpc with ${\widehat{\mathit{\ell }}}_{\mathrm{osc}}$ and $\widehat{\mathit{k}}$ having the same orientation. We note that ℓ_osc is in principle degenerate with cosmology, and in a full cosmological analysis it should also be varied. This being said, the impact of varying this parameter in such an analysis has been checked internally within the Euclid Collaboration and has been found to be minimal.

3.4.2. Lagrangian resummation

Pybird adopts a different approach to resummation, explicitly developed in Senatore & Zaldarriaga (2014). The details of the numerical implementation are explained in D’Amico et al. (2021b).

The Lagrangian resummation starts from the expression of the overdensity as a functional of the displacement field ψ(y, t):=l(y, t)−y, where l(y, t) is the final position of the particle and y is the initial position. The overdensity is given as

$$\begin{array}{cc}\hfill 1+\delta (\mathit{x},t)& ={\displaystyle \int {\mathrm{d}}^3}y{\delta }_{\mathrm{D}}^{(3)}[\mathit{x}-\mathit{y}-\mathit{\psi }(\mathit{y},t)]\hfill \\ \hfill & ={\displaystyle \int {\mathrm{d}}^3}y\int \frac{{\mathrm{d}}^{3}k}{{(2\pi )}^3}{\mathrm{e}}^{\mathrm{i}\mathit{k}·[\mathit{x}-\mathit{y}-\mathit{\psi }(\mathit{y},t)]}.\hfill \end{array}$$(62)

In redshift space, one needs to separate the components perpendicular and parallel to the line of sight

$$\begin{array}{cc}\hfill 1+{\delta }^s(\mathit{s},t)& ={\displaystyle \int {\mathrm{d}}^3}y{\delta }_{\mathrm{D}}^{(2)}[{\mathit{s}}_{\perp }-{\mathit{y}}_{\perp }-{\mathit{\psi }}_{\perp }^s(\mathit{y},t)]\hfill \\ \hfill & \times {\delta }_{\mathrm{D}}^{(1)}[s_{\Vert }-y_{\Vert }-{\psi }_{\Vert }^s(\mathit{y},t)-{\dot{\psi }}_{\Vert }^s(\mathit{y},t)/H]\hfill \\ \hfill & ={\displaystyle \int {\mathrm{d}}^3}y{\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }}{{(2\pi )}^2}\int \frac{\mathrm{d}{k}_{\Vert }}{2\pi }}{\mathrm{e}}^{\mathrm{i}{\mathit{k}}_{\perp }·[{\mathit{s}}_{\perp }-{\mathit{y}}_{\perp }-{\mathit{\psi }}_{\perp }^s(\mathit{y},t)]}\hfill \\ \hfill & \times {\mathrm{e}}^{\mathrm{i}{k}_{\Vert }[s_{\Vert }-y_{\Vert }-{\psi }_{\Vert }^s(\mathit{y},t)-{\dot{\psi }}_{\Vert }^s(\mathit{y},t)/H]}.\hfill \end{array}$$(63)

Here ${\mathit{\psi }}_{\Vert }^s(\mathit{y},t)=({\mathit{\psi }}^s(\mathit{y},t)·\widehat{z})\widehat{z}$ and ψ_⊥^s(y, t) = ψ^s(y, t)−ψ_∥^s(y, t). From Eq. (63), it is easy to derive the power spectrum

$$\begin{array}{c}\hfill P({\mathit{k}}_{\perp },k_{\Vert })={\displaystyle \int {\mathrm{d}}^3}y{\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{y}}{K}_r(\mathit{k},\mathit{y}),\end{array}$$(64)

where

$$\begin{array}{cc}\hfill K_r(\mathit{k},\mathit{y})& =\langle {\mathrm{e}}^{\mathrm{i}{\mathit{k}}_{\perp }·[{\mathit{\psi }}_{\perp }^s(\mathit{y})-{\mathit{\psi }}_{\perp }^s(\mathbf{0})]}{\mathrm{e}}^{\mathrm{i}{k}_{\Vert }[{\psi }_{\Vert }^s(\mathit{y})-{\psi }_{\Vert }^s(\mathbf{0})+\frac{{\dot{\psi }}_{\Vert }^s(\mathit{y},t)-{\dot{\psi }}_{\Vert }^s(\mathbf{0},t)}{H}]}\rangle \hfill \\ \hfill & =exp{\displaystyle \sum _{N=1}^{\infty }}\frac{1}{N!}\langle {\{\mathrm{i}{\mathit{k}}_{\perp }·{\mathbf{\Delta }}_{\perp }(\mathit{y})+i{k}_{\Vert }[{\mathrm{\Delta }}_{\Vert }(\mathit{y})+{\dot{\mathrm{\Delta }}}_{\Vert }(\mathit{y})/H]\}}^N\rangle ,\hfill \end{array}$$(65)

having defined Δ(y):=ψ(y)−ψ(0). Now, one only resums the linear displacement field on large scales, while expanding perturbatively the short-scale displacements and the density field. In practice, one cuts off the displacement integrals in the exponential, as done in Eqs. (60), (61). After some lengthy algebra, one arrives at the formula

$$\begin{array}{c}\hfill {P^{\ell }(k)|}_N={\displaystyle \sum _{j=0}^N\sum _{{\ell }^{\prime }}}4\pi {(-\mathrm{i})}^{{\ell }^{\prime }}{\displaystyle {\int }_0^{\infty }}\mathrm{d}r{r}^2{Q}_{||N-j}^{\ell {\ell }^{\prime }}(k,r){\xi }_j^{{\ell }^{\prime }}(r),\end{array}$$(66)

where P^ℓ(k)|_N is the resummed (ℓ-th multipole of the) power spectrum up to order N, ξ_j^ℓ(r) is the j-th loop order term in Eulerian perturbation theory of the (ℓ-th multipole of the) correlation function, and $Q_{\left|\right|N-j}^{\mathcal{l}{\mathcal{l}}^{\prime }}(k,r)$ is given by

$$\begin{array}{cc}\hfill Q_{||N-j}^{\ell {\ell }^{\prime }}(k,r)& ={\mathrm{i}}^{{\ell }^{\prime }}\frac{2\ell +1}{2}\hfill \\ \hfill & \times {\displaystyle {\int }_{-1}^1}\mathrm{d}{\mu }_k{\displaystyle \int \frac{{\mathrm{d}}^2\widehat{r}}{4\pi }}{\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{r}}F{||}_{N-j}(\mathit{k},\mathit{r}){\mathcal{L}}_{\ell }({\mu }_k){\mathcal{L}}_{{\ell }^{\prime }}({\mu }_r),\hfill \end{array}$$(67) $$\begin{array}{cc}\hfill {F||}_{N-j}(\mathit{k},\mathit{r})& =T_{0,r}(\mathit{k},\mathit{r})T_{0,r}^{-1}{||}_{N-j}(\mathit{k},\mathit{r}),\hfill \end{array}$$(68) $$\begin{array}{cc}\hfill T_{0,r}(\mathit{k},\mathit{r})& =exp\{-k^2{\mathrm{\Sigma }}^2(\mathit{k},\mathit{r})\},\hfill \end{array}$$(69)

where ${\mathrm{d}}^2\widehat{r}$ is an angular integration and μ = μ_k is now an integration variable. Here $T_{0,r}^{-1}{||}_{N-j}(\mathit{k},\mathit{y})$ means that we need to perturbatively expand (in powers of $P_{11}^{\delta \delta }$) $T_{0,r}^{-1}$ up to loop order N − j.

If in PyBird the option OPTIRESUM is chosen, the correlation function is split in a smooth part and a baryon acoustic oscillations peak, and Eq. (66) is only applied to the baryon acoustic oscillations peak. The result is then added to the smooth power spectrum. For all PyBird results shown in this paper we have used OPTIRESUM.

We find that OPTIRESUM and the WnW decomposition are highly consistent, but are both significantly discrepant with the Lagrangian resummation predictions. The reader is directed to Appendix C for more details on these comparisons.

3.4.3. Resummation in scale-dependent models

We discuss now the extension of the infrared-resummation procedure to models where the linear growth is scale-dependent, as in the f(R) model of Sect. 2.4 or when massive neutrinos are considered (see next section). Both infrared-resummation schemes discussed above consider contributions of the form (we discuss only real space, for simplicity)

$$\begin{array}{cc}& -\frac{k^2}{3}{\displaystyle \int \frac{{\mathrm{d}}^{3}q}{{(2\pi )}^3}}{\mathrm{e}}^{-q^2/{\mathrm{\Lambda }}_{\mathrm{IR}}^2}\frac{P_{11}^{\delta \delta }(q)}{q^2}[P(k)-P(|\mathit{k}-\mathit{q}|)],\hfill \end{array}$$(70)

where the power spectra inside square brackets are computed at a given order in perturbation theory. The cutoff ensures that q ≪ k in most of the integration domain. If P(k) contains a new scale, as is the case for ℓ_osc in P^w(k), the difference in parentheses is either of order $P_{11}^{\mathrm{w}}(k)q^2{\ell }_{\mathrm{osc}}^2$ (for q ℓ_osc ≪ 1) or $P_{11}^{\mathrm{w}}(k)$ (for q ℓ_osc ≫ 1, as the oscillations in the second term inside parentheses average out in this limit). Both regimes give a contribution to Eq. (70) which is k² enhanced with respect to other perturbative contributions of the same order, and therefore should be resummed.

If, on the other hand, no oscillatory feature is present in P(k), the difference inside parentheses can be approximated as

$$\begin{array}{cc}& [P(k)-P(|\mathit{k}-\mathit{q}|)]\simeq -\frac{1}{2}\frac{{\mathrm{d}}^{2}P(k)}{\mathrm{d}{k}^2}{q}^2{x}^2,\hfill \end{array}$$(71)

with

$$\begin{array}{c}\hfill x=\frac{\mathit{k}·\mathit{q}}{kq}.\end{array}$$(72)

For an approximately scale-invariant power spectrum, the second derivative gives a contribution of order P(k)/k² and, consequently, Eq. (70) is not k²-enhanced. Therefore, in resummation schemes which focus on the wiggly part of the power spectrum, as in WnW or in the OPTIRESUM option of PyBird, only leading contributions are resummed. In other schemes, as the full Lagrangian one, also subleading contributions are included, The difference between the two methods is of order 2-loop terms. The wiggle-no-wiggle method has been derived as an approximation of the Lagrangian resummation in Lewandowski & Senatore (2020).

In scale-dependent models the separation between leading and next-to-leading contributions could be, potentially, complicated by the presence of a new scale. However, considering for definiteness f(R), we can show that it is not the case. Indeed, even assuming that the whole scale-dependence of the function μ(k, a) of Eq. (18) is inherited by the power spectrum would imply a contribution to Eq. (71) proportional to

$$\begin{array}{c}\hfill \frac{{\partial }^{2}ln\mu (k,a)}{\partial k^2}\ll {\ell }_{\mathrm{osc}}^2,\end{array}$$(73)

where the inequality holds in the whole relevant k-range for the f(R) models considered in this work. We conclude that the extra-scale-dependence induced in f(R) is always subdominant with respect to the one induced by baryon acoustic oscillations, and does not require any modification of the infrared-resummation schemes with respect to those for scale-independent models. Pybird implements this scheme, and computes the linear power spectrum in Eqs. (60), (61) by taking the scale-dependent linear growth into account.

Since the scale-dependence in f(R) is much stronger than that induced by massive neutrinos, the above conclusions extend also to models with nonvanishing neutrino masses.

3.5. Massive neutrinos

With the measurements of flavour oscillations (Fukuda et al. 1998; Ahmed et al. 2004), massive neutrinos entered the realm of standard physics. These have been shown to have a significant effect on cosmological observables (Lesgourgues & Pastor 2006; Wright et al. 2019; Bose et al. 2021). We look to include these effects in Eq. (41).

Since galaxies are biased tracers of the CDM plus baryon field only (see Villaescusa-Navarro et al. 2014; Castorina et al. 2014; Costanzi et al. 2013, for example), we only consider the effects of massive neutrinos on the CDM+baryon (cb) spectrum, and opt to include them only through the linear cb power spectrum $P_{11}^{(\mathrm{cb})}(k,a)$. We do not include higher order massive neutrino effects in the kernels F_n^app and G_n^app (‘app’ stands for approximate), which are calculated assuming a single CDM matter fluid. This differs slightly from the prescription of Wright et al. (2019), which was shown to work very well by comparing to simulations. In particular, we do not include the effects of massive neutrinos on the growth rate f. We have checked that this effect is sub-percent for neutrino masses considered here.

In practice, we use the following initial linear spectrum in the 1-loop integrals

$$\begin{array}{c}\hfill P_{11,\mathrm{i}}^{(\mathrm{cb})}(k)=\frac{P_{11}^{(\mathrm{cb})}(k,a)}{F_1^{(\mathrm{m})2}(k,a)},\end{array}$$(74)

where F₁^(m)(k, a) assumes no massive neutrino contribution but the same total matter contribution to the background evolution and Poisson equation, Ω_m. We have checked that using Ω_cb instead makes negligible difference on the final RSD power spectrum predictions. On the other hand, $P_{11}^{(\mathrm{cb})}(k,a)$ has the full massive neutrino and modified gravity or dark energy dependence, which we calculate using the Boltzmann code MGCAMB (Zucca et al. 2019). The goal is to check for any deterioration of this approximation between the ΛCDM and beyond-ΛCDM comparisons with the simulation measurements.

3.6. Approximations

The efficient calculation of the 1-loop power spectrum can be challenging, particularly when the modifications include some sort of scale-dependence, for example in f(R) gravity. The numerical approach, as adopted by the MG-Copter code, is both time consuming and prone to numerical inaccuracies making it ill-suited for fast and comprehensive statistical analyses. On the other hand, the FFT approach, as adopted by PyBird and PBJ is highly computationally efficient and free from numerical instabilities. The challenge is extending the FFT to non-trivial kernels as introduced by modifications to gravity. Alternatively, one can apply various approximations to the computation which are easily implemented in current FFT frameworks.

Motivated by these issues, we examine various approximations that can alleviate computational cost while maintaining the desired degree of accuracy for the halo power spectrum in redshift space. These approximations, relevant for the loop integrals implicit in the first two lines of Eq. (44), are computed by rescaling the approximate second and third order kernels with the appropriate linear growth factors as follows,

$$\begin{array}{cc}& F_2({\mathit{k}}_1,{\mathit{k}}_2)=\frac{F_1(k_1)F_1(k_2)}{F_1^{\mathrm{app}}(k_1)F_1^{\mathrm{app}}(k_2)}{F}_2^{\mathrm{app}}({\mathit{k}}_1,{\mathit{k}}_2),\hfill \\ \hfill & G_2({\mathit{k}}_1,{\mathit{k}}_2)=\frac{f(k_{12})F_1(k_1)F_1(k_2)}{f^{\mathrm{app}}(k_{12})F_1^{\mathrm{app}}(k_1)F_1^{\mathrm{app}}(k_2)}{G}_2^{\mathrm{app}}({\mathit{k}}_1,{\mathit{k}}_2),\hfill \\ \hfill & F_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)=\frac{F_1(k_1)F_1(k_2)F_1(k_3)}{F_1^{\mathrm{app}}(k_1)F_1^{\mathrm{app}}(k_2)F_1^{\mathrm{app}}(k_3)}{F}_3^{\mathrm{app}}({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3),\hfill \\ \hfill & G_3({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)=\frac{f(k_{123})F_1(k_1)F_1(k_2)F_1(k_3)}{f^{\mathrm{app}}(k_{123})F_1^{\mathrm{app}}(k_1)F_1^{\mathrm{app}}(k_2)F_1^{\mathrm{app}}(k_3)}\hfill \\ \hfill & \times G_3^{\mathrm{app}}({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3),\hfill \end{array}$$(75)

where we have dropped the time dependence for compactness and we remind the reader that f(k) = − G₁(k)/F₁(k) is the growth rate. Recall that ‘app’ stands for approximate. The approximations we consider for the kernels are

1. The Einstein–de Sitter approximation (EdS): F₁^app(k, a) = F₁(k, a) and G₁^app(k, a) = G₁(k, a), but the higher-order kernel’s scale dependence are given by the standard Einstein–de Sitter universe (Ω_m(a) = 1) expressions, F_n^EdS and $G_n^{\mathrm{EdS}}$ (see Bernardeau et al. 2002). The time dependence is then given by the appropriate factors of F₁ and G₁, for example F₂^app(k₁, k₂) = F₁(k₁, a) F₁(k₂, a) F₂^EdS(k₁, k₂).

2. The unscreened approximation (USA): F₁^app(k, a) = F₁(k, a) and G₁^app(k, a) = G₁(k, a) and the higher-order kernels are computed by solving the perturbation evolution equations using S(k) = 0 (see Eq. 7), that is we ignore the effects of higher-order mode-coupling terms responsible for screening.

3. The ΛCDM-screened approximation (ΛCDM-scr): F₁^app(k) and G₁^app(k) are the ΛCDM growth factors which are scale-independent to a very good approximation: F₁^ΛCDM(k, a)≈F₁^ΛCDM(a) and G₁^ΛCDM(k, a)≈G₁^ΛCDM(a). The higher-order kernels are computed using the approximate expression for S(k) as given in Appendix B. This expression does not depend on the integrated momentum mode and hence allows it to be decomposed using a FFT approach.

4. Simulations

Before discussing our results, we list the various simulations we consider in Table 2. We describe each of these in detail below.

4.1. The ELEPHANT simulations

The Extended LEnsing PHysics using ANalaytic ray Tracing (ELEPHANT) simulations (Cautun et al. 2018; Hernández-Aguayo et al. 2019) are a suite of five independent realisations of the GR (ΛCDM), F6 (f_R0 = −10⁻⁶), F5 (f_R0 = −10⁻⁵), N5 (H₀r_c/c = 5) and N1 (H₀r_c/c = 1) models, where F stands for an f(R) model and N for a DGP model. The simulations were run with the ECOSMOG adaptive mesh refinement code (Li et al. 2012, 2013) and they follow the evolution of 1024³ dark-matter particles in a cubical box of length 1024 h⁻¹Mpc, giving a mass resolution of m_p = 7.78 × 10¹⁰ h⁻¹ M_⊙. Their initial conditions were generated at z = 49 using the Zel’dovich approximation with the MPgraphic code (Prunet et al. 2008) and the cosmological parameters are consistent with those of the Wilkinson Microwave Anisotropy Probe (WMAP) data release nine collaboration (Hinshaw et al. 2013),

$$\begin{array}{c}\hfill \{{\mathrm{\Omega }}_{\mathrm{b},0},{\mathrm{\Omega }}_{\mathrm{m},0},h,n_{\mathrm{s}},A_{\mathrm{s}}\}=\{0.046,0.281,0.697,0.971,2.297\times {10}^{-9}\},\end{array}$$

where A_s is the amplitude of primordial scalar perturbations, n_s is the scalar spectral index and Ω_b, 0 is the baryon density fraction today. The halo catalogues were constructed with the rockstar halo finder (Behroozi et al. 2013), where we chose the M_200c halo mass definition, which is the mass enclosed within a sphere of radius r_200c with 200 times the critical density of the universe.

These simulations offer one of the largest suites of modified gravity simulations available. There are larger volume simulations (see for example Arnold et al. 2019), but single realisation, and so the effective volume the ELEPHANT simulations provide is deemed optimal. The unavailability of simulations with a volume comparable to Euclid largely stems from the fact that running large volume modified gravity simulations is currently extremely computationally challenging (see Arnold et al. 2022, and references therein).

4.2. The DAKAR simulations

We made use of the DAKAR simulations presented in Baldi & Simpson (2017) which were run using a modified version of the GADGET-2N-body code (Springel 2005) which consistently implements the effects of the momentum exchange within the dark sector. The simulations consisted of 1024³ dark matter particles in a periodic cosmological box of length 1 h⁻¹ Gpc, evolved from a starting redshift of z_i = 99. The resulting CDM particle mass is m_c = 8 × 10¹⁰ h⁻¹ M − ⊙ and the spatial resolution is ϵ = 24 h⁻¹ kpc. The cosmological parameters are

$$\begin{array}{c}\hfill \{{\mathrm{\Omega }}_{\mathrm{b},0},{\mathrm{\Omega }}_{\mathrm{m},0},h,n_{\mathrm{s}},A_{\mathrm{s}}\}=\{0.048,0.308,0.678,0.966,2.115\times {10}^{-9}\}.\end{array}$$

We refer the interested reader to Baldi & Simpson (2017) for a more extended description of the simulations and of the modified N-body code.

We note that these simulations have a very limited volume, far smaller than that to be probed by Euclid. As will be seen in Sections 5.2 and 5.3, the approximations considered for the Dark Scattering model are extremely good even without introducing bias degrees of freedom, and more so when considering recent lensing and clustering constraints on this model (Carrilho et al. 2023; Carrion et al. 2024).

4.3. The DEMNUni simulations

The “Dark Energy and Massive Neutrino Universe” (DEMNUni) simulations (Carbone et al. 2016; Parimbelli et al. 2022) have been produced with the aim of investigating large-scale structures in the presence of massive neutrinos and dynamical dark energy, and they were conceived for the nonlinear analysis and modelling of different probes, including dark matter, halo, and galaxy clustering (see Castorina et al. 2015; Moresco et al. 2017; Zennaro et al. 2018; Ruggeri et al. 2018; Bel et al. 2019; Parimbelli et al. 2021, 2022; Guidi et al. 2023; Baratta et al. 2023; Gouyou Beauchamps et al. 2023, and Carella et al., in prep.), weak lensing, CMB lensing, SZ and ISW effects (Roncarelli et al. 2015; Carbone et al. 2016; Fabbian et al. 2018; Hernández-Molinero et al. 2024), cosmic void statistics (Kreisch et al. 2019; Schuster et al. 2019; Verza et al. 2019, 2022, 2023; Vielzeuf et al. 2023), and cross-correlations among these probes (Cuozzo et al. 2024). The DEMNUni simulations combine a good mass resolution with a large volume to include perturbations both at large and small scales. They are characterised by a softening length ε = 20 h⁻¹ kpc, a comoving volume of 8 h⁻³ Gpc³ filled with 2048³ dark matter particles and, when present, 2048³ neutrino particles. The simulations are initialised at z_i = 99 with Zeldovich initial conditions. The initial power spectrum is rescaled to the initial redshift via the rescaling method developed in Zennaro et al. (2017). Initial conditions are then generated with a modified version of the N-GenIC software, assuming Rayleigh random amplitudes and uniform random phases. The DEMNUni simulations were run using the tree particle mesh-smoothed particle hydrodynamics (TreePM-SPH) code P-Gadget3 (Springel 2005), specifically modified as in Viel et al. (2010) to account for the presence of massive neutrinos. This modified version of P-Gadget3 follows the evolution of CDM and neutrino particles, treating them as two separated collisionless components.

The reference cosmological parameters are chosen to be close to the baseline Planck 2013 cosmology (Planck Collaboration XVI 2014)

$$\begin{array}{c}\hfill \{{\mathrm{\Omega }}_{\mathrm{b},0},{\mathrm{\Omega }}_{\mathrm{m},0},h,n_{\mathrm{s}},A_{\mathrm{s}}\}=\{0.05,0.32,0.67,0.96,2.127\times {10}^{-9}\}.\end{array}$$

Given these values, the reference (i.e., the massless neutrino case) CDM-particle mass resolution is $m_{\mathrm{CDM}}^{\mathrm{p}}=8.27\times {10}^{10}{h}^{-1}{M}_{\odot }$ and is decreased according to the mass of neutrino particles, in order to keep the same Ω_m, 0 among all the DEMNUni simulations. In fact, massive neutrinos are assumed to come as a particle component in a three mass-degenerate scenario. Therefore, to keep Ω_m, 0 fixed, an increase in the massive neutrino density fraction yields a decrease in the CDM density fraction.

5. Results

5.1. Setup

We tested the various beyond-ΛCDM modelling approaches against measurements from N-body simulations, with the ΛCDM case serving as a performance benchmark. We considered the redshift space power spectrum monopole and quadrupole as they contain most of the cosmological and gravitational information, and only considered z = 1, as this will be one of the lowest targeted redshifts of the Euclid clustering probe (see, for example, Euclid Collaboration 2020; Amendola et al. 2018; Laureijs et al. 2011), with the higher redshifts relying less on the accuracy of nonlinear modelling, and where modified gravity or dark energy effects on the power spectra are less pronounced.

The analyses we performed are limited by the speed of MG-Copter which is capable of providing the spectra predictions without employing any of the approximations outlined in Sect. 3.6. Because of this, we restricted our analysis to the calculation of the minimum χ² statistic between the theory predictions and simulation measurements, comparing this statistic between the various modelling choices. This gave us a basic measure of the applicability of each choice without having to perform more expensive Bayesian parameter inference analyses. Our results aim at informing and optimising such future parameter posterior analyses (for example, D’Amico et al., in prep.).

We performed two separate analyses, one for CDM and one for CDM halos. The first tested the validity of the approximations outlined in Sect. 3.6 as at the level of CDM we have a very reduced nuisance parameter set with which to fit the simulations. The second analysis tested the robustness of all the modelling approaches detailed in Sect. 3; the various kernel approximations of Sect. 3.6, the two RSD models, the two bias schemes, and massive neutrino modelling as described in Sect. 3.

The χ² statistic is given by

$$\begin{array}{cc}\hfill {\chi }^2(k_{\max })=& \frac{1}{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}}{\displaystyle \sum _{k=k_{\min }}^{k_{\max }}\sum _{\ell ,{\ell }^{\prime }=0,2}}[P_{\ell ,\mathrm{data}}^s(k)-P_{\ell ,\mathrm{model}}^s(k)]\hfill \\ \hfill & \times {\mathrm{Cov}}_{\ell ,{\ell }^{\prime }}^{-1}(k)[P_{{\ell }^{\prime },\mathrm{data}}^s(k)-P_{{\ell }^{\prime },\mathrm{model}}^s(k)],\hfill \end{array}$$(76)

where P_ℓ^s is the ℓ^th multipole of the redshift space CDM or halo power spectrum and Cov_{ℓ, ℓ′} is the Gaussian covariance matrix between the different multipoles. In the case of a Gaussian covariance, the number of degrees of freedom are given by N_dof = 2 N_bins − N_params, where N_bins is the number of k-bins summed over and N_params is the number of free parameters in the theoretical model. In this case, N_params is just be the number of nuisance parameters (bias and RSD) as we fixed the cosmological ones to the fiducial values.

We used a Gaussian, linear covariance between the multipoles (see Appendix C of Taruya et al. 2010, for details). This has been shown to reproduce N-body results up to k ≤ 0.3 h Mpc⁻¹ at z = 1 (Taruya et al. 2010). We used the simulation volume (see Table 2) to calculate the covariance, and in the case of halos, we assumed b₁ as measured from the simulations. No supersample covariance was included. The effects of this was shown to be minimal in past surveys (Wadekar et al. 2020), and will be investigated in the context of Euclid in an upcoming Euclid paper. As we mostly use simulations with a finite size, no supersample contribution should be included in these cases.

Further, no shot noise term was considered. We have found that including a shot noise term based on the simulation-measured halo number density yields unreasonably good fits by our theoretical prescriptions to the simulations. More specifically, we found χ² ≈ 1 for k_max ≤ 0.6 h Mpc⁻¹ and beyond in some cases. This is likely because the models are efficient at fitting shot noise above a certain scale. We have checked that the relative behaviour of the χ² between ΛCDM and beyond-ΛCDM models is unchanged around χ² ∼ 1 with and without the shot noise contribution. This gives us confidence that the χ² ∼ 1 criteria we used is still meaningful to compare the relative goodness-of-fit between the various approximations considered.

On this note, a few things should be mentioned. The actual measurements from Euclid will encode a larger number density of tracers than the simulations considered making us more sensitive to inaccuracies in the modelling. For example, the ELEPHANT halo measurements at z = 1 have $\overline{n}\approx 3.1\times {10}^{-4}{h}^3{\mathrm{Mpc}}^{-3}$ while Euclid is expected to have over twice this density of tracers (Euclid Collaboration 2020). Further, Euclid will measure galaxies not halos, which come with additional bias modelling considerations. Lastly, despite having fixed the cosmology and gravity to their fiducial values, such goodness-of-fit comparisons do not ensure an absence of bias in the recovered parameters in a full posterior analysis, which can emerge from complex degeneracies between nuisance and cosmological parameters, and the shape of the full parameter space posterior distribution.

The simulations also all have smaller volumes than the expected Euclid survey volume, V(0.9 ≤ z ≤ 1.1) = 7.94 h⁻³ Gpc³ (Euclid Collaboration 2020), except for DEMNUni, which has a slightly larger volume. This limits our capacity to make robust statements for Euclid based on the simulation measurements alone.

To address these issues, we performed some analyses on mock data produced using the theoretical prescriptions. In this case we could employ a larger volume in the covariance as well as the Euclid-estimated shot noise contribution, $\overline{n}=6.86\times {10}^{-4}{h}^3{\mathrm{Mpc}}^{-3}$ (Euclid Collaboration 2020). In what follows, this number density and V = 8.8 h⁻³ Gpc³ is assumed in the covariance whenever we refer to a Euclid-like setting. This volume is larger than the expected volume of Euclid at this redshift bin, which accounts partially for any additional information coming from higher redshift bin observations.

On this point, we remark that in all our analyses we only vary the redshift space and bias model nuisance parameters, which will be fit at every redshift bin in the Euclid spectroscopic analysis. This means our results will not likely degrade with the inclusion of other redshift bins. It is also true that all considered beyond-ΛCDM effects are diminished at higher redshifts, and so the tested approximations will perform better at all other Euclid redshift bins (see Table. 1 of Euclid Collaboration 2020, for the specific binning). It should however be noted that the lessening of modifications to ΛCDM can be compensated by the larger volumes probed at higher redshifts.

We use the largest value of k_max, the maximum applicable wave mode, at which theory gives a good fit to the simulation measurements, as our performance metric for testing the approximations and modelling approaches. To determine this, we followed the procedure of Bose et al. (2020a):

1. We performed a least-squares fit to the N-body data by varying the set of nuisance parameters but fixing all cosmological and gravitational parameters to their fiducial values. We repeated this for k-bins within 0.125 h Mpc⁻¹ ≤ k_max ≤ 0.300 h Mpc⁻¹.

2. We then calculated the 1σ confidence intervals (Δχ²) on a χ² distribution using N_dof degrees of freedom.

3. The model validity scale k_max was then chosen as the largest value of k_max which has [χ²(k_max)−Δχ²(k_max)] ≤ 1. A χ² value of one indicates a good fit to the data.

This procedure provides a fair estimate of the range of validity of the models in which we can recover the fiducial cosmology with a 1σ criterion. This procedure was validated in ΛCDM in Bose et al. (2020a), Markovic et al. (2019) using a Markov Chain Monte Carlo analysis.

We imposed the physically motivated priors b₁, σ_v > 0 in the fits. Otherwise all fits assume an extremely wide prior on all other nuisance parameters, representing minimal prior information. We made this choice, as it is now known from the results of Carrilho et al. (2023), Simon et al. (2023), that the choice of narrow priors on nuisance parameters can lead to differences in the inferred marginal posteriors of cosmological parameters. While also shown in the aforementioned references that prior volume effects are important and increase with the size of the priors⁸, our results are not affected, since we are maximising the posterior, which is this case of large priors is equivalent to a maximum likelihood estimation of the parameters.

5.2. Approximation selection

We begun by performing direct spectra comparisons for the approximations outlined in Sect. 3.6 in order to evaluate the effectiveness of each approximation. To do this, we compared the monopole and quadrupole as predicted using the TNS model (see Eq. 36) under the various approximations against the full, most rigorous numerical calculation. We considered the TNS as opposed to the 1-loop SPT prediction (see Eq. 29) as it can probe smaller scales while not invoking too many degrees of freedom which may be degenerate with the approximation’s effects. We fixed the free parameter, σ_v, for all the predictions, to the exact value which is given in the figures. This value was chosen to be the best-fit value found by comparing to the respective dark matter P₀ and P₂ simulation measurements, where χ²(k_max)∼1 (see Eq. (76)). The hexadecapole comparisons are also shown to check the approximation’s accuracy at modelling even greater nonlinear RSD effects.

We note that for f(R) and DGP, the EdS and USA approximations have been examined in real space, and to a lesser extent in redshift space, in Bose & Koyama (2016), Bose et al. (2017, 2018b), Aviles et al. (2021). On the other hand, the ΛCDM-scr approximation has not been tested at all, and is presented for the first time in this work. For Dark Scattering and w₀w_aCDM, the EdS approximation was mildly tested at the level of the redshift space power spectrum in Bose et al. (2018a), Carrilho et al. (2021). Given that PyBird has recently been extended to include the exact calculation for both DGP (Piga et al. 2023) and w₀w_aCDM (D’Amico et al. 2021b), we restrict ourselves to f(R) and Dark Scattering in this section.

Despite Dark Scattering having been tested independently in other works, we also performed a limited test of the EdS approximation. We show this in Fig. 1 along with multipole errors coming from a Gaussian covariance with our Euclid-like volume at z = 1. In particular we considered the largest modification to ΛCDM, the CPL2 case (see Table 2). We see that the EdS deviates by more than 2% in the quadrupole at small scales, which may signal that this approximation is inadequate for Euclid. We have checked that lowering the coupling parameter, ξ, produces better consistency but, more importantly, changing the RSD degree of freedom can account for the entire deviation down to less than 0.5% for monopole and quadrupole, and brings all multipole ratios well within the monopole error assuming the Euclid-like volume at this redshift bin. Given this, we find the EdS to be a good approximation for low interaction strengths within the Dark Scattering model. We investigate this further in the next section. Further, we note that in the limit ξ → 0 we recover w₀w_aCDM and so we confirm that for mild w₀w_aCDM modifications, the EdS approximation is very good.

In Fig. 2 we show the f(R) case. Here, scale dependencies in the linear growth factor and rate cause the EdS and USA approximations to break down completely, with deviations significantly larger than the Euclid-like measurement errors. This has already been shown in the literature (see, for example, Aviles et al. 2021). Remarkably though, the ΛCDM-scr approximation is sub-percent consistent with the exact kernel calculation for all scales considered and well within the Euclid-like error bands. This provides a promising implementation for analyses pipelines which can use the FFT approach with this approximation. We investigate this further in the next section.

To summarise, the EdS approximation works very well for models inducing a scale-independent modification to the linear growth of structure, such as DGP, w₀w_aCDM or Dark Scattering. In the case of the former, the absence of a scalar field mass term and any environment dependence of the modification, leaves its impact on scales below the crossover scale r_c uniform. Similarly, the latter two models give no additional scale dependencies in the Poisson equation nor the Euler or continuity equations Eqs. (5), (6), (10). This ensures deviations to the scale dependencies of the standard EdS kernels are minimal. On the other hand, f(R) gravity comes with a scalar field potential term and associated mass, which induces a scale-dependent modification of the linear growth of structure. It also has the more complex, environment-dependent, Chameleon screening mechanism. This breaks the EdS approximation as well as makes the omission of screening more noticeable on the resulting power spectrum. The ΛCDM-scr approximation attempts to match the large and small scale behaviour of γ₂ and γ₃ without introducing additional scale dependencies from integrated wave modes, allowing the perturbative kernels to remain accurate (see Appendix B for more details).

5.3. Dark matter χ² tests

This analysis validated the necessary approximations currently needed by the FFT-based fast codes described in Table 1. In particular, we tested the EdS approximation for Dark Scattering, the EdS and USA for DGP and f(R), and the ΛCDM-scr approximation for f(R). We restricted all these tests to the TNS model of RSD as only this implementation in MG-Copter uses the exact kernel calculations. If the approximation does sufficiently well using the TNS, we expect that it should do equally well or better in the EFTofLSS case which makes use of more RSD fitting parameters. Further, our ultimate goal is to test these approximations at the level of CDM halos, which we do in the next section. As we are only considering CDM, to compute the minimum χ² as a function of k_max, we only vary the TNS RSD nuisance parameter, σ_v, in order to minimise the χ² statistic given in Eq. (76).

In Fig. 3 we show the results for all cases. The top panel show the ELEPHANT simulation DGP model results. In this case, the low modification (N5, exact) does equally well as the ΛCDM case. The high modification (N1, exact) does slightly worse, which is expected as nonlinearities are amplified by the significantly enhanced growth, causing SPT to break down at smaller k. In both cases, the EdS combined with the USA approximation does very well, following the goodness of fit of the exact prediction.

In the middle panel, we show the ELEPHANT simulation f(R) results. Here the high modification case (F5, exact) does equally well as ΛCDM, while we see a much better fit in the low modification case (F6, exact). This is in part due to the degeneracy between damping effects induced by modified gravity and the TNS damping parameter σ_v. The added damping of the quadrupole due to f(R) increases the effectiveness of σ_v in damping both multipoles appropriately. This only works to an extent as illustrated by the F5 case. As the damping from f(R) is also redshift dependent, we do not expect this result to hold for all redshifts.

Reassuringly, the ΛCDM-scr approximation does only slightly worse than the exact computation. On the other hand, the EdS combined with the USA approximation does significantly worse than the exact computation in the F5 case. This was expected given Fig. 2. Despite this poor performance, it is yet to be seen if this approximation is sufficient (or necessary) at the level of CDM halos, when we introduce a number of bias parameters.

Finally, the bottom panel shows the Dark Scattering simulation results. In this case, we find the EdS approximation is a good approximation in all cases, even in the high interaction case (CPL2), with all cases following the ΛCDM goodness of fit very closely. This was expected given the top right panel of Fig. 1. Based on this, we do not consider Dark Scattering in the next section.

Finally, we comment on the particularly small χ² in the ELEPHANT simulations at small k_max. This is likely due to the very small scatter in these measurements as they are the averages of five realisations, yielding small χ² despite the smaller error expressed in the covariance. In the absence of this scatter at linear scales, the nonlinear models applied here, with their additional degree of freedom, will over-fit the data as seen. Moving to nonlinear scales, the size of the errors (larger in DAKAR), begins to play a bigger role as the measurements from both simulations are very low in scatter. This is also seen to some extent in the next section.

5.4. Halo χ² tests

Having validated the various approximations in the TNS case at the level of CDM, we move to compare the goodness of fit in the beyond-ΛCDM models to the ΛCDM case at the level of CDM halos. These are biased tracers of the CDM distribution and so serve as a proxy for galaxies, and a means to test the modelling prescriptions in the presence of bias degrees of freedom.

In the TNS case, we primarily adopted the model used in the flagship BOSS survey analysis (Beutler et al. 2017), which has two bias degrees of freedom b₁ and b₂, as well as a shot noise parameter N. We also considered a Q-bias model independently for the f(R) cases (see Eq. 54) which has three degrees of freedom, on top of which we also include the shot-noise parameter N. These are in addition to the TNS RSD parameter, σ_v.

In the EFTofLSS analyses for biased tracers we adopted a model with three perturbative biases, b₁, b₂ and b_𝒢2, one shot-noise parameter, N, and two counterterms, ${\stackrel{\sim }{c}}_0$ and ${\stackrel{\sim }{c}}_2$. Adding more shot-noise or counterterm parameters gradually improves the reached k_max but, at the same time, could erase the information about the specific model analysed, since many of these EFTofLSS parameters may be degenerate with the beyond-ΛCDM ones. This will be investigated in an upcoming Euclid paper (D’Amico et al., in prep.). Our main goal here is to compare the ΛCDM to beyond-ΛCDM scenarios, which should not depend on the inclusion of the additional shot-noise parameters. Since in EFTofLSS analyses these parameters are fixed and/or marginalised analytically, we fixed all the higher-order nuisance parameters to zero.

We summarise the best-fit results for all models with Eulerian bias in Table 3. The TNS model with the Q-bias prescription fits are shown in Table 4.

5.4.1. f(R) and DGP

We present the dependence on the χ² as a function of k_max for the ELEPHANT simulations in Fig. 4. Here we show both EFTofLSS and TNS model predictions for DGP (top two panels) and f(R) (bottom three panels). We find that the Eulerian bias expansion, keeping scale-independent bias coefficients, is effective for the DGP case. Both levels of modification follow the ΛCDM trend very closely, with the high modification doing slightly better at the mid-range of scales presented. Further, the EdS and USA are very good approximations, following the exact computation almost exactly as was seen in Fig. 3. This holds for both TNS and EFTofLSS cases.

In the f(R) case, we tested two different bias models and both EdS and USA as well as the ΛCDM-scr approximations. The bias models we examine are the Q-bias and Eulerian bias (see Sect. 3.3). For the Eulerian bias, we assume EdS kernels for the higher-order bias terms and scale-independent bias coefficients. Interestingly, when the bias degrees of freedom were introduced, the EdS and USA approximations seem to become much more applicable, with the inaccuracy it incurs likely being absorbed by the bias degrees of freedom. We also note that there is little difference between the TNS Eulerian bias prescription and the Q-bias model, despite the Q-bias having one additional parameter.

It is also interesting to note that all f(R) cases do better in terms of their fit to the simulations (a lower χ²) than the ΛCDM fits. This could be because of an increased efficiency of the bias and RSD degrees of freedom to damp both multipoles appropriately in the presence of f(R) effects, as was noted in the CDM case. We also see this enhanced fit in the EFTofLSS case which applies the EdS and USA approximation. This improved efficiency of the RSD models to capture the full shape of the power spectrum when f(R) effects are present does not necessarily mean we will gain more information as it also suggests there are no distinct f(R) features at these scales, but rather just an overall damping of the spectrum, highly degenerate with σ_v or counterterms.

This being said, we do note that an improved fit to the full shape does not mean there are no additional features in the multipoles coming from f(R), specifically ones that are not pure damping in nature. These effects may become relevant with smaller error bars. Moreover, both the damping and these possible additional features are crucially dependent on the modification strength, which we can see clearly in Fig. 4 with the F6 case fitting the data marginally better than the ΛCDM case, while the F5 case fitting the multipole measurements significantly better.

An important point here is that we have fit the simulation data using a linear covariance that assumes the simulation volume, which is much less than that which Euclid will be probing. We do however omit shot noise which would improve the fits at small scales. We have found that including shot noise contributions, based on number densities from the halo catalogue measurements, enables the TNS model to fit the data to well beyond k > 0.5 h Mpc⁻¹. This is likely due to the model simply fitting the shot noise which it can do efficiently.

To check whether or not the approximations are good enough for Euclid, we checked the multipoles using the best-fit parameters determined from the simulations. In particular, we investigated the case which is expected to perform worst: the USA combined with EdS for the F5 case, which we saw in Fig. 2 performed terribly for dark matter. In Fig. 5 we show the ratio of the monopole and quadrupole for the TNS EdS and USA approximation for F5 to the exact solution with the ELEPHANT-based simulation best-fit parameters found in Table 3. We also overlay as bands the expected cosmic variance and shot noise errors expected from the Euclid survey at the z = 1 redshift bin (Euclid Collaboration 2020). As dotted lines we show the EdS and USA approximation prediction but without refitting the nuisance parameters and using the same fits as the exact solution. This highlights how the nuisance parameters are effectively accounting for the approximation, even to within Euclid precision. Without this freedom the approximation fails, even before the shot noise begins to dominate.

To further investigate whether or not the approximations applied here would be valid in a Euclid-like setting, we prepared a mock data vector using the TNS model and exact SPT kernel predictions⁹, and attached to it a covariance computed with Euclid-like shot noise and volume (V = 8.8 h⁻³ Gpc³ and $\overline{n}=6.86\times {10}^{-4}{h}^3{\mathrm{Mpc}}^{-3}$, see Sect. 5.1). We further added scatter to the mock data vector using this covariance. Then, we fitted the theoretical prescriptions using approximations to this data vector and checked the goodness of fit as a function of k_max. The fits are shown in Fig. 6. All approximations do as well as the exact calculation (which the mock data vector was produced with), and give a χ² ≈ 1 for a large range of k_max, exceeding that determined by fitting the simulations.

We also fitted this mock data vector in the f(R) cases using a pure ΛCDM modelling (μ(k, a) = 1 and S(k) = 0), which also shows an excellent fit to the mock data. This suggests that the nuisance parameters are completely degenerate with f(R) effects, at least for |f_R0|≤10⁻⁵ and over the range of scales that our models are valid within, even in a Euclid-like setup.

One way of checking this hypothesis is to observe the behaviour of the nuisance parameters as a function of k_max. We find, in the TNS case, that the best-fit linear bias b₁ and σ_v are ∼5% and ∼10% consistent between all approximations and the exact results over the full range of k_max, with b₁ also being consistent with the simulation measurement. Interestingly, the best-fit values of σ_v in the ΛCDM-scr approximation match extremely well with the exact solution, which confirms it being a very good approximation for the SPT kernels. Differing values of σ_v for the other approximations and Q-bias model indicate that the effects of f(R) are being absorbed by a change in this parameter. Similarly, b₂ and N also vary significantly, both as a function of k_max and between the various approximations. Such degeneracies will likely be able to be broken by using a combination of redshift bins, giving far more information about the evolution of structure, as well as the bispectrum for example, which will yield more information on bias parameters, particularly if the modelling prescription for both statistics is consistent (Hashimoto et al. 2017; Philcox et al. 2022; Tsedrik et al. 2023; Ivanov et al. 2022). We leave further investigation of this to future work which should employ large-volume simulations.

5.4.2. w₀w_aCDM and massive neutrinos

In Fig. 7 we show the χ² as a function of k_max for the DEMNUni simulations (see Table 2). The TNS and EFTofLSS model predictions for the massless neutrino cases are shown in the top two panels, while the bottom two panels show the massive neutrino cases (∑m_ν = 0.16 eV). We find that the EdS approximation does equally well to the exact kernel calculations for all cases, at least until the χ² ≈ 1. Further, all beyond-ΛCDM cases follow very similar trends to the ΛCDM case. This suggests that the EdS approximation is sufficiently accurate for w₀w_aCDM cases, and our linear-only treatment of massive neutrinos (see Sect. 3.5) is equally good in all cases.

As with the DGP and f(R) cases, we find that the TNS model shows slightly lower χ² over the full range of k_max considered than EFTofLSS, but both RSD models have roughly a k_max(χ² ≈ 1)≈0.19 h Mpc⁻¹. This is similar to the k_max for most of the ELEPHANT simulations (except F5), despite the DEMNUni simulations having a larger volume. This is likely compensated by the slightly lower σ₈ of these simulations, and hence slightly lower levels of nonlinearity, leading to a better performance of SPT.

The volume of the DEMNUni simulations is comparable to the volume Euclid will probe at z = 1 and so we did not perform any additional tests of the EdS approximation against SPT mock data. On the other hand, we did test the importance of including massive neutrinos. We created a mock data vector using the exact calculation of the SPT kernels and including massive neutrino effects assuming ∑m_ν = 0.16 eV in the linear P_cb power spectrum (see Sect. 3.5), again using the best-fit nuisance parameters found in Table 3. We then created a covariance with the Euclid-like volume and tracer number density (V = 8.8 h⁻³ Gpc³ and $\overline{n}=6.86\times {10}^{-4}{h}^3{\mathrm{Mpc}}^{-3}$) and a linear bias as given in Table 3. Using this covariance, we added scatter to the mock data. We then fitted this mock data vector using the massless neutrino modelling. The fits of the various theoretical predictions to the mock data vectors are shown in Fig. 8.

Surprisingly, we find all massless neutrino fits, both in ΛCDM and w₀w_aCDM, follow the massive neutrino modelling fits almost perfectly, with χ² ≈ 1 at all k_max considered. This further suggests that the effects of massive neutrinos with masses ∑m_ν ≤ 0.16 eV are degenerate with nuisance parameters, at least at a fixed redshift bin and only with the clustering power spectrum. This is expected as the FoG damping effect, controlled by the free parameter σ_v in the TNS model or counterterms in EFTofLSS, are likely highly degenerate with the effects of massive neutrinos, which act to damp power on small scales.

We await further tests to be performed in D’Amico et al., in prep. to further explore these results. Unlike screening, the inclusion of massive neutrinos at the linear level is not significantly computationally expensive and so no approximation needs to be made in principle.

6. Conclusions

In this paper we have investigated various theoretical approximations in beyond-ΛCDM scenarios, necessary for the computational demands of forthcoming Euclid galaxy clustering analyses, as well as different RSD models. Before listing our results, as a reference for the reader, we provide a summary of previous related works and their main conclusions.

6.1. Previous work

The Einstein–de Sitter (EdS) approximation has been investigated in a number of works. Bose et al. (2018b) compare the EdS approximation to simulations and the exact 1-loop real space matter power spectrum in DGP using the EFTofLSS. They find this is a very good approximation, especially given the additional degrees of freedom of the EFTofLSS. Bose et al. (2017) perform Markov Chain Monte Carlo analyses on DGP simulations using the EdS TNS redshift space dark matter multipoles with and without screening effects, finding a 2σ bias is incurred on the growth rate, f, when screening is omitted at z = 1 for a survey of V_s = 20 h⁻³Gpc³, but not for V_s = 10 h⁻³Gpc³. Further, Carrilho et al. (2021) check the EdS approximation for the Dark Scattering real space spectra against the exact solutions and find sub-percent agreement for z = 1, for scales less than k ≈ 0.25 h Mpc⁻¹ and for ξ comparable to the values considered in this paper.

Regarding tracer bias schemes, D’Amico et al. (2021a) show that the scale-independent bias parameters are likely a good approximation for models of gravity inducing a scale-independent growth factor and rate of structure formation, without explicitly comparing to simulations. A number of other works have developed exact calculations for redshift space correlations for dark matter and biased tracers, but have not explicitly investigated the need for this by comparing, directly or through best-fit analyses, to simulations or mock data (Aviles & Cervantes-Cota 2017; Aviles et al. 2021; Valogiannis et al. 2020; Valogiannis & Bean 2019; Bose & Koyama 2016, 2017; Bose et al. 2018a), specifically when higher order nuisance parameters are included. We do however note that Aviles et al. (2021) find the local Lagrangian approximation to be valid for f(R) gravity.

6.2. This work

In this work we have considered the TNS and EFTofLSS models for RSD with an Eulerian bias expansion. We further considered the local Lagrangian bias relation, known to hold well for ΛCDM and f(R) (Aviles et al. 2021). For the scale-dependent theory, f(R), we also investigate the phenomenological Q-bias model (see Eq. (54)).

These approximations were checked by comparing the monopole and quadrupole of both dark matter and halo clustering to high quality numerical simulations as well as standard perturbation theory (SPT) based mock data generated with a minimal amount of approximations. For all comparisons, we have fixed cosmology and only varied the model nuisance parameters related to RSD, nonlinear clustering and tracer bias. Our main conclusions are summarised below.

6.2.1. Dark Matter

The EdS approximation without screening predictions gives an equally good χ² fit to the simulation dark matter multipoles as the predictions without these approximations in DGP. The EdS approximation is generally excellent for DGP, w₀w_aCDM, and Dark Scattering models, both at the level of dark matter and halos. It does however give a significantly worse fit for f(R). On the other hand, the ΛCDM-screened approximation (see Appendix B) for f(R) gravity, developed in this work, gives an equally good fit to the simulation dark matter multipoles as the exact SPT predictions. This approximation is suitable for a fast Fourier transform implementation.

6.2.2. Halos

The EdS approximation without screening predictions gives an equally good χ² fit to the simulation halo multipoles as the predictions without these approximations, both in f(R) gravity and DGP. Further, both Eulerian and Q-bias prescriptions give equally good fits to the simulation data, with the local Lagrangian bias relation seeming to be valid for all beyond-ΛCDM scenarios, with the χ² of these scenarios being comparable to that of ΛCDM.

We found that not modelling any f(R) effects, both in the bias terms and in the perturbation theory kernels, gives an equally good fit to the SPT-f(R) mock data as the exact kernel computation. This suggests such effects can efficiently be absorbed into higher-order bias and RSD nuisance parameters. In fact, we have noted a significant change in the RSD and higher order bias parameters from the fiducial mock values in these fits. Particularly prominent shifts were noted for the largest approximations, for example using the Einstein–de Sitter and unscreened approximations for f(R) gravity. These shifts also grew with the maximum scale, k_max, included in the fit.

Similarly, not modelling any massive neutrino effects gives an equally good fit to the SPT-ΛCDM and -w₀w_aCDM mock data as the modelling including massive neutrinos, for ∑m_ν ≤ 0.16 eV. This suggests such effects can also be efficiently absorbed into higher-order bias and RSD nuisance parameters. We note that changes in these parameters away from the fiducial mock values was also noted in the fits, in support of this claim.

Lastly, TNS and EFTofLSS both give a similar value of k_max such that |χ²(k_max)−1|≲0.15 when fit to the simulations, with the TNS yielding a slightly larger k_max systematically over all beyond-ΛCDM scenarios¹⁰. This range is roughly 0.19 h Mpc⁻¹ ≤ k_max ≤ 0.21 h Mpc⁻¹, with the exception being the strong f(R) modification where we find the TNS model can fit the data well up to k_max ≈ 0.26 h Mpc⁻¹. This is likely due to an enhanced efficiency of the phenomenological fingers-of-god RSD damping parameter σ_v when strong f(R) effects are at play.

We note that we do not vary cosmological parameters, and so these results, in particular the use of the various approximations considered for the perturbative kernels and bias, does not ensure an unbiased estimation of cosmology. They do however strongly suggest large degeneracies between screening and massive neutrino effects and higher-order bias as well as RSD effects. If validated, this would mean that a power spectrum monopole and quadrupole analysis at a single redshift will offer no significant information on modified gravity. This result was known in the case of massive neutrinos due to their degeneracies with bias parameters (see Marulli et al. 2011; Villaescusa-Navarro et al. 2018; Hahn et al. 2020, for example), but not in the case of modified gravity, which is a key finding of this work.

This result may change when including higher-order multipoles or statistics (García-Farieta et al. 2021), such as the hexadecapole or bispectrum, galaxy lensing, or informative priors on nuisance parameters, which will help to break clear degeneracies between beyond-ΛCDM effects and bias parameters (Chan & Blot 2017; Bose & Taruya 2018; Bose et al. 2020a; Markovic et al. 2019; Bose et al. 2019; Tsedrik et al. 2023; Kacprzak & Fluri 2022). Such informative priors must be chosen carefully though, as they can lead to biases on the inferred cosmology or gravitational model (Carrilho et al. 2023; Simon et al. 2023). If such priors can however be identified, or if we can identify a more efficient nuisance parameter set, through for example principal component analyses (see for example Eifler et al. 2015) or machine learning methods (Piras & Lombriser 2023), then the theoretical approximations used, particularly in the case of f(R) or scale-dependent models, must be revised. In this case a more accurate prescription, such as the ΛCDM-screened approximation, can be applied. Conversely, we can use the methods presented here to gain information on the appropriate width of nuisance parameter priors. For example, we find that σ_v changes by up to 10% when employing approximate perturbative kernels. So, any prior must be at least so wide about some central value if we choose to employ such an approximation in a real analysis.

In light of these findings, we do not advocate additional effort in either improving computational efficiency to calculate the exact beyond-ΛCDM SPT kernels or theoretically developing new, more accurate kernel approximations. We also deem the Eulerian bias expansion, with constant bias coefficients to be sufficient for the considered beyond-ΛCDM cases. The RSD models perform similarly and so we do not advocate one over the other. This being said, tools and theoretical prescriptions are already available which can compute the clustering multipoles highly accurately and efficiently for DGP, Dark Scattering, w₀w_aCDM, and massive neutrinos (D’Amico et al. 2021b; Noriega et al. 2022; Piga et al. 2023; Carrilho et al. 2023). A forthcoming Euclid paper (D’Amico et al., in prep.) will provide a more detailed Markov Chain Monte Carlo-based analyses of both RSD models.

In summary, the main contribution of this work to Euclid is to have identified useful approximations for determining scale cuts in beyond-ΛCDM scenarios, which will require both high-volume simulation measurements and computationally expensive Markov Chain Monte Carlo analyses. On this note, this work has instructed the minimal extension of current internal pipelines for beyond-ΛCDM theories, in particular for the TNS model. We have also developed a useful theoretical approximation (ΛCDM-screened) in the case where standard approximations in f(R)-gravity break under the combined power of the Euclid spectroscopic probe. This approximation can be integrated easily into codes employing the Fast Fourier transform method such as PyBird and PBJ.

Finally, we have also validated the clustering predictions from various SPT based codes being employed within the Euclid consortium. In particular, we have found the predictions from PBJ, Pybird and MG-Copter (now part of ReACT) (see Table 1) to be in sub-percent agreement for k ≤ 0.5 h Mpc⁻¹ (see Appendix C).

---

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 Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft-und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, 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). BB was supported by a UK Research and Innovation Stephen Hawking Fellowship (EP/W005654/2). AP is a UK Research and Innovation Future Leaders Fellow [grant MR/S016066/2]. CM and PC’s research for this project was supported by a UK Research and Innovation Future Leaders Fellowship [grant MR/S016066/2]. CM’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). MP acknowledges support by the MIUR ‘Progetti di Ricerca di Rilevante Interesse Nazionale’ (PRIN) Bando 2022 – grant 20228RMX4A. FV acknowledges partial support by the ANR Project COLSS (ANR-21-CE31-0029). NF is supported by the Italian Ministry of University and Research (MUR) through the Rita Levi Montalcini project “Tests of gravity on cosmic scales” with reference PGR19ILFGP and she also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021. CC acknowledges a generous CPU and storage allocation by the Italian Super-Computing Resource Allocation (ISCRA) as well as from the coordination of the “Accordo Quadro MoU per lo svolgimento di attivitá congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, together with storage from INFN-CNAF and INAF-IA2. FP acknowledges partial support from the INFN grant InDark, the Departments of Excellence grant L.232/2016 of the Italian Ministry of Education, University and Research (MUR) and from the European Union – Next Generation EU. FP also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021. We acknowledge the hospitality of the Institute for Fundamental Physics of the Universe (IFPU) of Trieste for the group meeting held there in December 2022. For the purpose of open access, the author(s) has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

References

Appendix A: Redshift space correction terms

In this appendix we present the general forms for the A, B and C corrections terms in Eq. (29) as well as their linear bias dependence. These can be matched to the EdS expressions in, for example, Scoccimarro et al. (1999), Bernardeau et al. (2002) by expanding Eq. (44). We note that these are the forms that are also found in MG-Copter, which can be found https://github.com/nebblu/ACTio-ReACTio/blob/master/reactions/src/SPT.cpp. Following Bose & Koyama (2016) we have

$$\begin{array}{cc}\hfill A(k,\mu )& =-(k\mu ){\displaystyle \int {\mathrm{d}}^3}q[\frac{{\mu }_q}{q}{B}_{\sigma }(\mathit{q},\mathit{k}-\mathit{q},-\mathit{k})\hfill \\ \hfill & +\frac{k\mu -q{\mu }_q}{{|\mathit{k}-\mathit{q}|}^2}{B}_{\sigma }(\mathit{k}-\mathit{q},\mathit{q},-\mathit{k})],\hfill \end{array}$$(A.1)

where ${\mu }_q=\widehat{\mathit{q}}·\widehat{\mathit{z}}$ and the cross bispectrum is given in terms of the general kernels and linear bias as

$$\begin{array}{cc}& B_{\sigma }({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)=2\hfill \\ \hfill & \times \{[b_1{F}_1(k_2)-{\mu }_{k_2}^2{G}_1(k_2)][b_1{F}_1(k_3)-{\mu }_{k_3}^2{G}_1(k_3)]G_2({\mathit{k}}_2,{\mathit{k}}_3)\hfill \\ \hfill & \times P_{11,\mathrm{i}}(k_2)P_{11,\mathrm{i}}(k_3)\hfill \\ \hfill & +[b_1{F}_1(k_3)-{\mu }_{k_3}^2{G}_1(k_3)][b_1{F}_2({\mathit{k}}_1,{\mathit{k}}_3)-{\mu }_{k_2}^2{G}_2({\mathit{k}}_1,{\mathit{k}}_3)]G_1(k_1)\hfill \\ \hfill & \times P_{11,\mathrm{i}}(k_1)P_{11,\mathrm{i}}(k_3)\hfill \\ \hfill & +[b_1{F}_1(k_2)-{\mu }_{k_2}^2{G}_1(k_2)][b_1{F}_2({\mathit{k}}_1,{\mathit{k}}_2)-{\mu }_{k_3}^2{G}_2({\mathit{k}}_1,{\mathit{k}}_2)]G_1(k_1)\hfill \\ \hfill & \times P_{11,\mathrm{i}}(k_1)P_{11,\mathrm{i}}(k_2)\}.\hfill \end{array}$$(A.2)

The B term is given as

$$\begin{array}{c}\hfill B(k,\mu )={(k\mu )}^2{\displaystyle \int {\mathrm{d}}^3}qF(\mathit{q})F(\mathit{k}-\mathit{q}),\end{array}$$(A.3)

with

$$\begin{array}{c}\hfill F(\mathit{q})=\frac{{\mu }_q}{q}{G}_1(q)[b_1{F}_1(q)-{\mu }_q^2{G}_1(q)]P_{11,\mathrm{i}}(q).\end{array}$$(A.4)

Finally, the C term is given as

$$\begin{array}{cc}\hfill C(k,\mu )& ={(k\mu )}^2{\displaystyle \int {\mathrm{d}}^3}p{\mathrm{d}}^{3}q{\delta }_{\mathrm{D}}(\mathit{k}-\mathit{p}-\mathit{q})P_{11,\mathrm{i}}(p)P_{11,\mathrm{i}}(q)\hfill \\ \hfill & \times \frac{{\mu }_p^2}{p^2}{G}_1^2(p){[b_1{F}_1(q)-{\mu }_q^2{G}_1(q)]}^2.\hfill \end{array}$$(A.5)

Appendix B: ΛCDM-screened for f(R)

In this approximation we have the following modifications to the Poisson equation

$$\begin{array}{cc}\hfill \mu (k,a)& =1+{\left(\frac{k}{a}\right)}^2\frac{1}{3\mathrm{\Pi }(k,a)},\hfill \end{array}$$(B.1) $$\begin{array}{cc}\hfill {\gamma }_2^{\mathrm{app}}(\mathit{k}-\mathit{q},\mathit{q},a)& ={\gamma }_2(\mathit{k},0,a)\frac{\mathrm{\Pi }(0,a_{\mathrm{f}})}{\mathrm{\Pi }(q,a_{\mathrm{f}})}\{\frac{\mathrm{\Pi }(0,a_{\mathrm{f}})}{\mathrm{\Pi }(q,a_{\mathrm{f}})}\hfill \\ \hfill & +\frac{[\mathrm{\Pi }(k,a_{\mathrm{f}})-\mathrm{\Pi }(0,a_{\mathrm{f}})][\mathrm{\Pi }(q,a_{\mathrm{f}})-\mathrm{\Pi }(0,a_{\mathrm{f}})]}{{\mathrm{\Pi }}^2(q,a_{\mathrm{f}})}\},\hfill \\ \hfill {\gamma }_2^{\mathrm{app}}(\mathit{k},-\mathit{q},a)& ={\gamma }_2(\mathit{k},0,a)\frac{\mathrm{\Pi }(0,a_{\mathrm{f}})}{\mathrm{\Pi }(q,a_{\mathrm{f}})},\hfill \\ \hfill {\gamma }_2^{\mathrm{app}}(\mathit{k},\mathit{q}-\mathit{k},a)& =\frac{q^2}{k^2}{\gamma }_2^{\mathrm{app}}(\mathit{k}-\mathit{q},\mathit{q},a),\hfill \end{array}$$(B.2)

where a_f is the final scale factor at which the power spectrum is computed. For γ₃ we use the following expressions for their cyclic permutations

$$\begin{array}{cc}& {\gamma }_3^{\mathrm{app}}(\mathit{q},\mathit{k},-\mathit{q},a)={\gamma }_3^{\mathrm{app}}(\mathit{q},-\mathit{q},\mathit{k},a)={\gamma }_3(0,\mathit{k},0,a)\frac{{\mathrm{\Pi }}^2(0,a_{\mathrm{f}})}{{\mathrm{\Pi }}^2(q,a_{\mathrm{f}})},\hfill \\ \hfill & {\gamma }_3^{\mathrm{app}}(\mathit{k},\mathit{q},-\mathit{q},a)={\gamma }_3(\mathit{k},0,0,a)\frac{{\mathrm{\Pi }}^2(0,a_{\mathrm{f}})}{{\mathrm{\Pi }}^2(q,a_{\mathrm{f}})}.\hfill \end{array}$$(B.3)

The terms depending on the loop momentum q in these expressions are evaluated at a fixed value for the scale factor, a = a_f, allowing a factorization of scale and time integrations. Moreover, the momentum dependence is only through powers of

$$\begin{array}{c}\hfill a_{\mathrm{f}}^2{\mathrm{\Pi }}^2(q,a_{\mathrm{f}})=q^2+m^2(a_{\mathrm{f}}),\end{array}$$(B.4)

where, from Eq. (21),

$$\begin{array}{c}\hfill m^2(a_{\mathrm{f}}):=\frac{a_{\mathrm{f}}^2{\mathrm{{\rm Y}}}^3(a_{\mathrm{f}})}{2f_0{(3{\mathrm{\Omega }}_{\mathrm{m},0}-4)}^2},\end{array}$$(B.5)

which are well suited for an extension of the FFT approach for the fast evaluation of loop integrals.

Appendix C: Code and implementation validation

We have performed some basic validation tests between the three perturbation theory codes used in this work. In particular, we have compared predictions for the 1-loop SPT dark matter monopole, quadrupole and hexadecapole (see Eq. 29) which involves no free nuisance parameters. We have done this for ΛCDM and DGP at z = 0.5 where nonlinearities are larger compared to z = 1, providing a better validation test. We have also compared the different resummation implementations described in Sect. 3.4. Additionally, we note that the χ² obtained from the different codes when comparing the 1-loop SPT prediction to the ΛCDM ELEPHANT simulations at z = 1 are extremely similar over a large range of k_max (see Eq. 76).

C.1. Code validation

In Fig. C.1 we show the ratio of the 1-loop SPT multipoles between the various codes for ΛCDM (left) and DGP (right) for k ≤ 0.5 h Mpc⁻¹ at z = 0.5, without infrared-resummation applied. For all cases we used the MG-Copter prediction as our reference spectrum. We have also included error bands from a ΛCDM Gaussian covariance assuming V = 3.8 h⁻³ Gpc³ at z = 0.5, estimated using Laureijs et al. (2011).¹¹ In the ΛCDM case, we find sub-0.5% agreement between MG-Copter, PBJ and an additional code by Atsushi Taruya (AT).¹² We note that PBJ and AT’s code both implement the expression for the 1-loop redshift space power spectrum as given by Matsubara (2008) while MG-Copter uses the expression given in Eq. (29) (see also Eq. 23 of Taruya et al. 2010). These comparisons highlight the equivalence of these expressions. Further, the oscillations seen in the top left panel are caused by the interpolation over k. These three predictions have all made use of the Einstein–de Sitter (EdS) approximation.

In the bottom left panel we compare the exact kernel predictions between MG-Copter and Pybird. We find sub-1% agreement in this case, with the monopole being sub-0.5%. Interestingly, we observe almost the same discrepancy between these two codes when applying the EdS approximation. We have found that these predictions can be brought into better agreement by retuning the Fast Fourier Transform (FFT) bias parameter in Pybird. We note that the spike observed in the monopole is due to numerical instabilities in MG-Copter.

We observe similar agreement in the DGP scenario. In this case we compare MG-copter to PBJ, both using the EdS and unscreened approximations as this is the only current implementation in PBJ. Again we observe sub-0.5% agreement. For the comparison with Pybird we use the exact time and screening predictions for both codes and observe sub-1% agreement between the predictions, with some oscillations caused by interpolation.

Finally, we find all comparisons are within the estimated measurement errors down to scales well beyond the expected validity regime of the predictions, which at this lower redshift will be lower than the rough estimates found in Table 3 (k_max ≤ 0.2 h Mpc⁻¹).

C.2. Resummation

In Fig. C.2 we move on to compare the different infrared-resummation implementations, considering only ΛCDM, and again for k ≤ 0.5 h Mpc⁻¹ at z = 0.5. In this case, PBJ and MG-Copter only have the wiggle-no-wiggle (WnW) decomposition approach (see Sect. 3.4), both within the EdS approximation, while Pybird has the option of both Lagrangian and the approximate OPTIRESUM method. The former codes also have slightly different WnW implementations, but despite these differences, we still observe sub-1% agreement. This suggests that the infrared-resummation prescription will have little impact on inferred parameter posteriors from upcoming data analyses.

On the other hand, applying the Lagrangian resummation method leads to significant differences, far outside the estimated measurement errors within the predictions validity range and for all multipoles. As discussed in Sect. 3.4.3 the presence of the baryon acoustic oscillations scale, ℓ_osc, provides a clear criterion to identify the leading contributions in k² at each order in perturbation theory. Both WnW and PyBird’s OPTIRESUM exploit this criterion and then resum the same class of contributions at each order. The full Lagrangian resummation, on the other hand, is applied to the full power spectrum, and therefore also includes subleading contributions like those in Eq. (71). These are genuine higher-order contributions but, on the other hand, there is no guiding principle on why to include them while neglecting other contributions of the same order.

The difference between WnW and the Lagragian resummation may be degenerate with counterterms, which is left to be seen. A full comparison between PBJ and Pybird at the posterior level will be conducted in a forthcoming Euclid paper (D’Amico et al., in prep.).

All Tables

All Figures