# Euclid preparation. Modelling spectroscopic clustering on mildly nonlinear scales in beyond-ΛCDM models

B Bose, P Carrilho, M Marinucci, C Moretti, M Pietroni, E Carella, L Piga, B.S Wright, F Vernizzi, C Carbone, et al.

► To cite this version:

B Bose, P Carrilho, M Marinucci, C Moretti, M Pietroni, et al. Euclid preparation. Modelling spectroscopic clustering on mildly nonlinear scales in beyond-ΛCDM models. Astron.Astrophys., 2024, 689, pp.A275. 10.1051/0004-6361/202348784. hal-04308845

HAL Id: hal-04308845

https://hal.science/hal-04308845v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Euclid preparation

XLIV. Modelling spectroscopic clustering on mildly nonlinear scales in
beyond-ΛCDM models

Euclid Collaboration: B. Bose¹,*, P. Carrilho¹, M. Marinucci²,³, C. Moretti⁴,¹,⁵,⁶, M. Pietroni²,⁷, E. Carella⁸,⁹,

L. Piga²,⁷, B. S. Wright¹⁰, F. Vernizzi¹¹, C. Carbone⁹, S. Casas¹², G. D'Amico²,⁷, N. Frusciante¹³,

K. Koyama¹⁴, F. Pace¹⁵,¹⁶,¹⁷, A. Pourtsidou¹,¹⁸, M. Baldi¹⁹,²⁰,²¹, L. F. de la Bella¹⁴, B. Fiorini¹⁴,¹⁰,

C. Giocoli²⁰,²², L. Lombriser²³, N. Aghanim²⁴, A. Amara¹⁴,²⁵, S. Andreon²⁶, N. Auricchio²⁰, S. Bardelli²⁰,

C. Bodendorf²⁷, D. Bonino¹⁷, E. Branchini²⁸,²⁹,²⁶, M. Brescia¹³,³⁰,³¹, J. Brinchmann³², S. Camera¹⁵,¹⁶,¹⁷,

V. Capobianco¹⁷, V. F. Cardone³³,³⁴, J. Carretero³⁵,³⁶, M. Castellano³³, S. Cavuoti³⁰,³¹, A. Cimatti³⁷, G. Congedo¹,

C. J. Conselice³⁸, L. Conversi³⁹,⁴⁰, Y. Copin⁴¹, A. Costille⁴², F. Courbin⁴³, H. M. Courtois⁴⁴, A. Da Silva⁴⁵,⁴⁶,

H. Degaudenzi⁴⁷, A. M. Di Giorgio⁴⁸, F. Dubath⁴⁷, C. A. J. Duncan³⁸,⁴⁹, X. Dupac⁴⁰, S. Dusini⁵⁰, M. Farina⁴⁸,

S. Farrens⁵¹, S. Ferriol⁴¹, P. Fosalba⁵²,⁵³, M. Frailis⁵⁴, E. Franceschi²⁰, S. Galeotta⁵⁴, B. Garilli⁹, B. Gillis¹,

A. Grazian⁵⁵, F. Grupp²⁷,⁵⁶, L. Guzzo⁸,²⁶,⁵⁷, S. V. H. Haugan⁵⁸, F. Hormuth⁵⁹, A. Hornstrup⁶⁰,⁶¹, K. Jahnke⁶²,

B. Joachimi⁶³, E. Keihänen⁶⁴, S. Kermiche⁶⁵, A. Kiessling⁶⁶, M. Kilbinger⁶⁷, T. Kitching⁶⁸, M. Kunz²³,

H. Kurki-Suonio⁶⁹,⁷⁰, S. Ligori¹⁷, P. B. Lilje⁵⁸, V. Lindholm⁶⁹,⁷⁰, I. Lloro⁷¹, D. Maino⁸,⁹,⁵⁷, E. Maiorano²⁰,

O. Mansutti⁵⁴, O. Marggraf⁷², K. Markovic⁶⁶, N. Martinet⁴², F. Marulli⁷³,²⁰,²¹, R. Massey⁷⁴, E. Medinaceli²⁰,

M. Meneghetti²⁰,²¹, G. Meylan⁴³, M. Moresco⁷³,²⁰, L. Moscardini⁷³,²⁰,²¹, D. F. Mota⁵⁸, E. Munari⁵⁴, S.-M. Niemi⁷⁵,

C. Padilla³⁵, S. Paltani⁴⁷, F. Pasian⁵⁴, K. Pedersen⁷⁶, W. J. Percival⁷⁷,⁷⁸,⁷⁹, V. Pettorino⁸⁰, S. Pires⁵¹, G. Polenta⁸¹,

M. Poncet⁸², L. A. Popa⁸³, L. Pozzetti²⁰, F. Raison²⁷, A. Renzi⁸⁴,⁵⁰, J. Rhodes⁶⁶, G. Riccio³⁰, E. Romelli⁵⁴,

M. Roncarelli²⁰, R. Saglia⁸⁵,²⁷, D. Sapone⁸⁶, B. Sartoris⁸⁵,⁵⁴, P. Schneider⁷², A. Secroun⁶⁵, G. Seidel⁶², M. Seiffert⁶⁶,

S. Serrano⁵²,⁸⁷,⁸⁸, C. Sirignano⁸⁴,⁵⁰, G. Sirri²¹, L. Stanco⁵⁰, J.-L. Starck⁶⁷, P. Tallada-Crespi⁸⁹,³⁶, A. N. Taylor¹,

I. Tereno⁴⁵,⁹⁰, R. Toledo-Moreo⁹¹, F. Torradeflot³⁶,⁸⁹, I. Tutusaus⁹², E. A. Valentijn⁹³, L. Valenziano²⁰,⁹⁴,

T. Vassallo⁸⁵,⁴, A. Veropalumbo²⁶,²⁹, Y. Wang⁹⁵, J. Weller⁸⁵,²⁷, G. Zamorani²⁰, J. Zoubian⁶⁵, E. Zucca²⁰,

A. Biviano⁵⁴,⁶, E. Bozzo⁴⁷, C. Burigana⁹⁶,⁹⁴, C. Colodro-Conde⁹⁷, D. Di Ferdinando²¹, J. Graciá-Carpio²⁷,

N. Mauri³⁷,²¹, C. Neissner³⁵,³⁶, Z. Sakr⁹⁸,⁹²,⁹⁹, V. Scottez¹⁰⁰,¹⁰¹, M. Tenti²¹, M. Viel⁶,⁵⁴,⁴,¹⁰²,⁵, M. Wiesmann⁵⁸,

Y. Akrami¹⁰³,¹⁰⁴, V. Allevato³⁰, S. Anselmi⁸⁴,⁵⁰,¹⁰⁵, M. Ballardini¹⁰⁶,¹⁰⁷,²⁰, F. Bernardeau¹¹,¹⁰⁸, S. Borgani¹⁰⁹,⁶,⁵⁴,¹⁰²,

S. Bruton¹¹⁰, R. Cabanac⁹², A. Cappi²⁰,¹¹¹, C. S. Carvalho⁹⁰, G. Castignani⁷³,²⁰, T. Castro⁵⁴,¹⁰²,⁶,⁵,

G. Cañas-Herrera⁷⁵,¹¹², K. C. Chambers¹¹³, A. R. Cooray¹¹⁴, J. Coupon⁴⁷, S. Davini²⁹, S. de la Torre⁴²,

G. De Lucia⁵⁴, G. Desprez¹¹⁵, S. Di Domizio²⁸,²⁹, H. Dole²⁴, A. Díaz-Sánchez¹¹⁶, J. A. Escartin Vigo²⁷,

S. Escoffier⁶⁵, P. G. Ferreira⁴⁹, I. Ferrero⁵⁸, F. Finelli²⁰,⁹⁴, L. Gabarra⁸⁴,⁵⁰, K. Ganga¹¹⁷, J. García-Bellido¹⁰³,

F. Giacomini²¹, G. Gozalias¹¹⁸,⁶⁹, D. Guinet⁴¹, A. Hall¹, S. Joudaki¹⁴,⁷⁷,⁷⁸, J. J. E. Kajava¹¹⁹,¹²⁰, V. Kansal¹²¹,¹²²,¹²³,

D. Karagiannis¹²⁴, C. C. Kirkpatrick⁶⁴, L. Legrand²³, A. Loureiro¹²⁵,¹²⁶, J. Macias-Perez¹²⁷, M. Magliocchetti⁴⁸,

R. Maoli¹²⁸,³³, M. Martinelli³³,³⁴, C. J. A. P. Martins¹²⁹,³², S. Matthew¹, M. Maturi⁹⁸,¹³⁰, L. Maurin²⁴,

R. B. Metcalf⁷³,²⁰, M. Migliaccio¹³¹,¹³², P. Monaco¹⁰⁹,⁵⁴,¹⁰²,⁶, G. Morgante²⁰, S. Nadathur¹⁴, Nicholas A. Walton¹³³,

L. Patrizii²¹, A. Pezzotta²⁷, V. Popa⁸³, C. Porciani⁷², D. Potter¹³⁴, M. Pöntinen⁶⁹, P. Reimberg¹⁰⁰, P.-F. Rocci²⁴,

A. G. Sánchez²⁷, A. Schneider¹³⁴, E. Sefusatti⁵⁴,¹⁰²,⁶, M. Sereno²⁰,²¹, A. Silvestri¹¹², A. Spurio Mancini⁶⁸,

J. Steinwagner²⁷, G. Testera²⁹, R. Teyssier¹³⁵, S. Toft⁶¹,¹³⁶, S. Tosi²⁸,²⁹,²⁶, A. Troja⁸⁴,⁵⁰, M. Tucci⁴⁷, J. Valiviita⁶⁹,⁷⁰,
and D. Vergani²⁰

(Affiliations can be found after the references)

Received 29 November 2023 / Accepted 24 May 2024

<footnote>* Corresponding author; ben.bose@ed.ac.uk</footnote>
ABSTRACT

Context. The *Euclid* space satellite mission will measure the large-scale clustering of galaxies at an unprecedented precision, providing a unique probe of modifications to the ΛCDM model.

Aims. We investigated the approximations needed to efficiently predict the large-scale clustering of matter and dark matter halos in the context of modified gravity and exotic dark energy scenarios. We examined the normal branch of the Dvali–Gabadadze–Porrati model, the Hu–Sawicki $f(R)$ model, a slowly evolving dark energy model, an interacting dark energy model, and massive neutrinos. For each, we tested approximations for the perturbative kernel calculations, including the omission of screening terms and the use of perturbative kernels based on the Einstein–de Sitter universe; we explored different infrared-resummation schemes, tracer bias models and a linear treatment of massive neutrinos; we investigated various approaches for dealing with redshift-space distortions and modelling the mildly nonlinear scales, namely the Taruya–Nishimishi–Saito prescription and the effective field theory of large-scale structure. This work provides a first validation of the various codes being considered by *Euclid* for the spectroscopic clustering probe in beyond-ΛCDM scenarios.

*Methods.* We calculated and compared the $\chi^2$ statistic to assess the different modelling choices. This was done by fitting the spectroscopic clustering predictions to measurements from numerical simulations and perturbation theory-based mock data. We compared the behaviour of this statistic in the beyond-ΛCDM cases, as a function of the maximum scale included in the fit, to the baseline ΛCDM case.

*Results.* We find that the Einstein-de Sitter approximation without screening is surprisingly accurate for the modified gravity cases when comparing to the halo clustering monopole and quadrupole obtained from simulations and mock data. Further, we find the same goodness-of-fit for both cases – the one including and the one omitting non-standard physics in the predictions. Our results suggest that the inclusion of multiple redshift bins, higher-order multipoles, higher-order clustering statistics (such as the bispectrum), and photometric probes such as weak lensing, will be essential to extract information on massive neutrinos, modified gravity and dark energy. Additionally, we show that the three codes used in our analysis, namely, PBJ, Pybird and MG-Copter, exhibit sub-percent agreement for $k \le 0.5 \hbar \text{ Mpc}^{-1}$ across all the models. This consistency underscores their value as reliable tools.

**Key words.** gravitation – cosmology: theory – dark energy – large-scale structure of Universe

# 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<fnref n="1" />, 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)<fnref n="2" />. 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<fnref n="3" />, these being the Vainshtein mechanism (Vainshtein

<footnote n="1">¹ This model assumes the dominating energy components of the Universe are a constant dark energy, Λ, and cold dark matter (CDM), with gravity described by Einstein's general relativity (GR).</footnote>

<footnote n="2">² See also https://www.euclid-ec.org/public/core-science/</footnote>

<footnote n="3">³ Screening is a theoretical mechanism used to suppress additional forces coming from modifications to gravity in the Solar System (Brax et al. 2021), where experiments show theory must accord with GR's predictions (Will 2018).</footnote>
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^2 = -c^2dt^2 + a^2(t)\delta_{ij}dx^i dx^j$, where $a$ represents the scale factor. The Hubble rate, $H$, is defined as $H := \dot{a}/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

$$ds^2 = -c^2 [1 + 2\Phi(\mathbf{x}, t)] dt^2 + a^2(t) [1 - 2\Psi(\mathbf{x}, t)] \delta_{ij} dx^i dx^j, \quad (1)$$

where the gravitational potential, $\Phi$, appears in the time-time component of the metric. To express the continuity and Euler equations, we use the rescaled CDM velocity divergence, $\theta := \nabla \cdot \mathbf{v}/(aH)$, and the CDM density contrast, $\delta := \delta\rho/\rho$, where $\rho$ and $\delta\rho$ are respectively the energy density background and perturbations.

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

$$\delta(\mathbf{k}, a) = \sum_{n=1}^{\infty} \delta_n(\mathbf{k}, a), \quad \theta(\mathbf{k}, a) = \sum_{n=1}^{\infty} \theta_n(\mathbf{k}, a), \quad (2)$$

where $\delta_n, \theta_n \sim \delta_1^n$, with $\delta_1$ 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{align}
\delta_n(\mathbf{k}, a) &= \frac{1}{(2\pi)^{3(n-1)}} \int d^3\mathbf{k}_1 \cdots d^3\mathbf{k}_n \delta_D(\mathbf{k} - \mathbf{k}_{1\cdots n}) \nonumber \\
&\quad \times F_n(\mathbf{k}_1, \dots, \mathbf{k}_n, a) \delta_{1,i}(\mathbf{k}_1) \cdots \delta_{1,i}(\mathbf{k}_n), \nonumber \\
\theta_n(\mathbf{k}, a) &= \frac{1}{(2\pi)^{3(n-1)}} \int d^3\mathbf{k}_1 \cdots d^3\mathbf{k}_n \delta_D(\mathbf{k} - \mathbf{k}_{1\cdots n}) \nonumber \\
&\quad \times G_n(\mathbf{k}_1, \dots, \mathbf{k}_n, a) \delta_{1,i}(\mathbf{k}_1) \cdots \delta_{1,i}(\mathbf{k}_n), \tag{3}
\end{align}$$

where $\mathbf{k}_{1\cdots n} = \mathbf{k}_1 + \cdots + \mathbf{k}_n$, $\delta_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

$$f := q \frac{d \ln \delta_1}{d \ln a} = -G_1/F_1, \quad (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 ($\Phi$) to the density perturbation$^4$

$$\begin{align}
a \delta'(\mathbf{k}, a) + \theta(\mathbf{k}, a) \nonumber \\
&= - \int \frac{d^3 \mathbf{k}_1 d^3 \mathbf{k}_2}{(2\pi)^3} \delta_D(\mathbf{k} - \mathbf{k}_{12}) \alpha(\mathbf{k}_1, \mathbf{k}_2) \theta(\mathbf{k}_1, a) \delta(\mathbf{k}_2, a), \tag{5}
\end{align}$$

$$\begin{align}
a \theta'(\mathbf{k}, a) + \left(2 + \frac{aH'(a)}{H(a)} + \frac{\Xi(a)}{H(a)}\right) \theta(\mathbf{k}, a) - \left(\frac{ck}{aH(a)}\right)^2 \Phi(\mathbf{k}, a) \nonumber \\
= -\frac{1}{2} \int \frac{d^3 \mathbf{k}_1 d^3 \mathbf{k}_2}{(2\pi)^3} \delta_D(\mathbf{k} - \mathbf{k}_{12}) \beta(\mathbf{k}_1, \mathbf{k}_2) \theta(\mathbf{k}_1, a) \theta(\mathbf{k}_2, a), \tag{6}
\end{align}$$

$$- \left( \frac{ck}{aH(a)} \right)^2 \Phi(\mathbf{k}, a) = \frac{3\Omega_m(a)}{2} \mu(\mathbf{k}, a) \delta(\mathbf{k}, a) + S(\mathbf{k}, a), \quad (7)$$

where a prime denotes a derivative with respect to the scale factor, $\Omega_m$ is the total matter density fraction and $\Xi$ encodes the effects of a possible additional drag force, coming from, say, an interaction within the dark sector (see Sect. 2.2). The kernels $\alpha(\mathbf{k}_1, \mathbf{k}_2)$ and $\beta(\mathbf{k}_1, \mathbf{k}_2)$ are given by

$$\alpha(\mathbf{k}_1, \mathbf{k}_2) = 1 + \frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{|\mathbf{k}_1|^2}, \quad (8)$$

$$\beta(\mathbf{k}_1, \mathbf{k}_2) = \frac{(\mathbf{k}_1 \cdot \mathbf{k}_2) |\mathbf{k}_1 + \mathbf{k}_2|^2}{|\mathbf{k}_1|^2 |\mathbf{k}_2|^2}. \quad (9)$$

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

$$\begin{align}
S(\mathbf{k}, a) = \int \frac{d^3 \mathbf{k}_1 d^3 \mathbf{k}_2}{(2\pi)^3} \delta_D(\mathbf{k} - \mathbf{k}_{12}) \gamma_2(\mathbf{k}_1, \mathbf{k}_2, a) \delta(\mathbf{k}_1, a) \delta(\mathbf{k}_2, a) \nonumber \\
+ \int \frac{d^3 \mathbf{k}_1 d^3 \mathbf{k}_2 d^3 \mathbf{k}_3}{(2\pi)^6} \delta_D(\mathbf{k} - \mathbf{k}_{123}) \gamma_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, a) \nonumber \\
\times \delta(\mathbf{k}_1, a) \delta(\mathbf{k}_2, a) \delta(\mathbf{k}_3, a), \tag{10}
\end{align}$$

<footnote>$^4$ The exact kernels can be obtained by solving these equations numerically as described in Taruya (2016), Bose & Koyama (2016) or by using a Fast Fourier Transform decomposition.</footnote>
where the functions $\gamma_2$ and $\gamma_3$ are additional kernels coming from higher-order scalar field interactions in the Klein-Gordon equation (see Koyama et al. 2009, for example). For $\Lambda$CDM, we have $\mu(k, a) = 1$ and $\Xi = \gamma_2 = \gamma_3 = 0$. We now look at the forms of these equations for various theories beyond-$\Lambda$CDM.

## 2.1. $w_0w_a$CDM

For a dark energy fluid the equation of state parameter is defined as $w = p_{DE}/(\rho_{DE}c^2)$, where $p_{DE}$ and $\rho_{DE}$ are the pressure and the background energy density of the fluid. In $\Lambda$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):

$$w(a) = w_0 + (1 - a) w_a, \qquad (11)$$

where $w_0$ 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 $\Lambda$CDM at their $\Lambda$CDM values, that is $\mu(k, a) = 1$ and $\Xi = \gamma_2 = \gamma_3 = 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 $\mu(k, a) = 1$ and $\gamma_2 = \gamma_3 = 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, $\theta_{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

$$\Xi(a) = [1 + w(a)] \xi c \rho_{DE}(a), \qquad (12)$$

where $\xi := \sigma_D/(m_c c^2)$. Here $m_c$ is the dark matter particle mass and $\sigma_D$ is the cross-section of the interaction. $\xi$ is treated as a free parameter of the theory. It is positive (or zero) and proportional to the strength of the modification to $\Lambda$CDM, with $\Lambda$CDM being recovered for $\xi \to 0$ or $w \to -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 $\Xi = 0$ in Eq. (6). The relation between $\Phi$ and $\delta$ is described by Eq. (7), with (see for example Bose & Koyama 2016; Bose et al. 2018b)

$$\mu(k, a) = 1 + \frac{1}{3\beta(a)}, \qquad (13)$$

$$\gamma_2(\mathbf{k}_1, \mathbf{k}_2, a) = - \left[ \frac{H_0}{H(a)} \right]^2 \frac{1}{24\beta^3(a)\Omega_{\text{rc}}} \left( \frac{\Omega_{\text{m},0}}{a^3} \right)^2 (1 - \mu_{1,2}^2), \qquad (14)$$

$$\gamma_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, a) = \left[ \frac{H_0}{H(a)} \right]^2 \frac{1}{144\beta^5(a)\Omega_{\text{rc}}^2} \left( \frac{\Omega_{\text{m},0}}{a^3} \right)^3 \times (1 - \mu_{2,3}^2)(1 - \mu_{1,23}^2), \qquad (15)$$

where $H_0$ and $\Omega_{\text{m},0}$ are the Hubble parameter and total matter density fraction today,

$$\beta(a) := q1 + \frac{H(a)}{H_0} \frac{1}{\sqrt{\Omega_{\text{rc}}}} \left( 1 + \frac{a H'(a)}{3H(a)} \right). \qquad (16)$$

The third order kernel, $\gamma_3$, needs to be symmetrised and $\mu_{i,j} = \hat{k}_i \cdot \hat{k}_j$ is the cosine of the angle between $\mathbf{k}_i$ and $\mathbf{k}_j$ (recall $\mathbf{k}_i = \mathbf{k}_i + \mathbf{k}_j$). Here $\Omega_{\text{rc}} := c^2/(4r_c^2 H_0^2)$, where $r_c$ represents the scale above which gravity deviates from GR, meaning $r_c \to \infty$, or $\Omega_{\text{rc}} \to 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 $\Lambda$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 $\Xi = 0$, and have the following modifications to the Poisson equation (Koyama et al. 2009; Taruya et al. 2014; Bose & Koyama 2016)

$$\mu(k, a) = 1 + \left(\frac{k}{a}\right)^2 \frac{1}{3\Pi(k, a)}, \qquad (17)$$

$$\gamma_2(\mathbf{k}_1, \mathbf{k}_2, a) = -\frac{3}{16} \left(\frac{k H_0}{a H(a)}\right)^2 \left(\frac{\Omega_{\text{m},0}}{a^3}\right)^2 \frac{\Upsilon^5(a)}{f_0^2 (3\Omega_{\text{m},0} - 4)^4} \times \frac{1}{\Pi(k, a) \Pi(k_1, a) \Pi(k_2, a)}, \qquad (18)$$

and

$$\gamma_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, a) = \frac{1}{32} \left(\frac{k H_0}{a H(a)}\right)^2 \left(\frac{\Omega_{\text{m},0}}{a^3}\right)^3 \frac{1}{\Pi(k, a) \Pi(k_1, a) \Pi(k_2, a) \Pi(k_3, a)} \times \left\{ -5 \frac{\Upsilon^7(a)}{f_0^3 (3\Omega_{\text{m},0} - 4)^6} + \frac{9}{2} \frac{1}{\Pi(k_{23}, a)} \left[ \frac{\Upsilon^5(a)}{f_0^2 (3\Omega_{\text{m},0} - 4)^4} \right]^2 \right\}, \qquad (19)$$

where we have implicitly assumed the constraint coming from the 1-loop integral: $\mathbf{k} = \mathbf{k}_1 + \mathbf{k}_2$ for $\gamma_2$ and $\mathbf{k} = \mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3$ for $\gamma_3$. Again, we note that $\gamma_3$ must be symmetrised. The functions
Π and $\Upsilon$ are given by

$$ \Pi(k, a) = \left(\frac{k}{a}\right)^2 + \frac{\Upsilon^3(a)}{2f_0(3\Omega_{m,0} - 4)^2}, \qquad (20) $$

$$ \Upsilon(a) = \frac{\Omega_{m,0} + 4a^3(1 - \Omega_{m,0})}{a^3}. \qquad (21) $$

Here we set $f_0 = c^2|f_{R0}|/H_0^2$, 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} \to 0$. As with the DGP case, we also assume a $\Lambda$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

$$ (2\pi)^3 \delta_D(\mathbf{k} + \mathbf{k}') P^{ab}(k, a) := \langle \varphi^a(\mathbf{k}, a) \varphi^b(\mathbf{k}', a) \rangle, \qquad (22) $$

where we have introduced the duplet $\varphi^a := (\delta, \theta)$, and $\langle \cdot \rangle$ denotes an ensemble average. By expanding $\delta$ and $\theta$ using Eq. (2), we can perturbatively express the power spectrum as a leading-order term plus a next-to-leading-order (nlo) term

$$ P^{ab}(k, a) = P_{11}^{ab}(k, a) + P_{\text{nlo}}^{ab}(k, a) + \dots, \qquad (23) $$

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

$$ P_{\text{nlo}}^{ab}(k, a) = P_{22}^{ab}(k, a) + P_{13}^{ab}(k, a) + P_{31}^{ab}(k, a), \qquad (24) $$

where we have used the notation

$$ (2\pi)^3 \delta_D(\mathbf{k} + \mathbf{k}') P_{ij}^{ab}(k, a) = \langle \varphi_i^a(\mathbf{k}, a) \varphi_j^b(\mathbf{k}', a) \rangle. \qquad (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

$$ P_{1-\text{loop}}^{ab} := P_{11}^{ab}(k, a) + P_{\text{nlo}}^{ab}(k, a). \qquad (26) $$

#### 3.2. Redshift space

So far we have worked in ‘real’ space coordinates, $\mathbf{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, $\mathbf{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

$$ \mathbf{s} = \mathbf{r} + \frac{v_z(\mathbf{r}, a)}{a H(a)} \hat{\mathbf{z}}, \qquad (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

$$ \delta^s(s, a) = \left| \frac{\partial s}{\partial \mathbf{r}} \right|^{-1} [1 + \delta(\mathbf{r}, a)] - 1, \qquad (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),

$$ P_{1-\text{loop}}^{s}(k, \mu) = P_{1-\text{loop}}^{\delta\delta}(k) + 2\mu^2 P_{1-\text{loop}}^{\delta\theta}(k) + \mu^4 P_{1-\text{loop}}^{\theta\theta}(k) \\ + A(k, \mu) + B(k, \mu) + C(k, \mu) \\ - k^2 \mu^2 \tilde{\sigma}_{\mathbf{v}}^2 [F_1(k) + G_1(k) \mu^2]^2 P_{11,i}^{\delta\delta}(k), \qquad (29) $$

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

$$ \tilde{\sigma}_{\mathbf{v}}^2 = \frac{1}{6\pi^2} \int dq G_1^2(q) P_{11,i}^{\delta\delta}(q). \qquad (30) $$

We refer the reader to Bose & Koyama (2016), Taruya et al. (2010) for expressions for $A(k, \mu)$, $B(k, \mu)$, and $C(k, \mu)$. 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<fnref n="6" /> (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;

<footnote n="5">⁵ Not to be confused with the linear Poisson equation modification, $\mu(k, a)$, which always appears with its arguments.</footnote>

<footnote n="6">⁶ Download ReACT: https://github.com/nebblu/ACTio-ReACTio</footnote>
- The publicly available PyBird code<fnref n="7" /> (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, \mu, \gamma_2, \gamma_3$ and $\Xi$.

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_0w_a$CDM, 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)

$$P_{\text{EFT}}^{s}(k, \mu) = P_{1-\text{loop}}^{s}(k, \mu) + \tilde{P}_{\text{ctr}}^{s}(k, \mu), \quad (31)$$

where the second term, $\tilde{P}_{\text{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

$$\tilde{P}_{\text{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). \quad (32)$$

This model has therefore a total of three nuisance parameters at the dark matter level, {$\tilde{c}_0, \tilde{c}_2, \tilde{c}_4$}, with dimensions $h^{-2}$ 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

$$P_{\text{ctr}}^{s}(k, \mu) = 2 P_{11}^{\delta\delta}(k) \left(\frac{k}{k_M}\right)^2 (b_1 + f \mu^2) (c_{\text{ctr}} + c_{r,1}\mu^2 + c_{r,2}\mu^4), \quad (33)$$

where $k_M = 0.7 h \text{ Mpc}^{-1}$ parametrises the inverse spatial extension scale of galaxies. The total number of nuisance parameters is therefore the same: {$c_{\text{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 $\tilde{c}_4 = c_{r,2} = 0$. This yields the following mapping between the two different counterterm bases

$$c_{\text{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}, \quad (34)$$

$$c_{r,1} = - \frac{35 k_M^2 (b_1 \tilde{c}_2 - \tilde{c}_0)}{35 b_1^2 + 30 b_1 f + 3 f^2}. \quad (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_1(k, a)/F_1(k, a)$, is both time and scale-dependent, and the linear power spectrum has additional scale dependencies coming from $F_1(k, a)$: $P_{11}^{\delta\delta} = F_1^2(k, a)P_{11,i}^{\delta\delta}$.

### 3.2.3. Taruya–Nishimichi–Saito

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

$$P_{\text{TNS}}(k, \mu) = D_{\text{FoG}}(\mu^2 k^2 \sigma_v^2) [P_{1-\text{loop}}^{\delta\delta}(k) + 2\mu^2 P_{1-\text{loop}}^{\delta\theta}(k) + \mu^4 P_{1-\text{loop}}^{\theta\theta}(k) + A(k, \mu) + B(k, \mu)], \quad (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_{\text{FoG}}$ (Taruya et al. 2010). The damping factor $D_{\text{FoG}}$ encodes a phenomenological description of the FoG damping. This is a function of the velocity dispersion, $\sigma_v$, which is taken to be a free nuisance parameter. Here we assume this function takes a Lorentzian form so that

$$D_{\text{FoG}}(k^2 \mu^2 \sigma_v^2) := \frac{1}{1 + (k^2 \mu^2 \sigma_v^2)/2}. \quad (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, {$\sigma_v$}, with units $h^{-1}$ 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.

<footnote n="7">⁷ Download PyBird: https://github.com/pierrexyz/pybird.</footnote>
<table><caption>**Table 1.** Summary of codes used in this work.</caption><thead><tr><th>Name</th><th>Solving method</th><th>Speed [ms]</th><th>Beyond-ΛCDM Models</th><th>RSD Model</th><th>Resummation</th></tr></thead><tbody><tr><td>PBJ</td><td>FFT</td><td>O(10)</td><td>DGP<sub>USA,EdS</sub>, DS<sub>EdS</sub>, w<sub>0</sub>w<sub>a</sub>CDM<sub>EdS</sub></td><td>TNS, EFTofLSS</td><td>WnW</td></tr><tr><td>PyBird</td><td>FFT</td><td>O(10)</td><td>DGP, w<sub>0</sub>w<sub>a</sub>CDM</td><td>EFTofLSS</td><td>Lagrangian</td></tr><tr><td>MG-Copter</td><td>Solve Eqs. (5), (6) per k</td><td>O(10000)</td><td>f(R), DGP, DS, w<sub>0</sub>w<sub>a</sub>CDM</td><td>TNS, EFTofLSS</td><td>WnW</td></tr></tbody></table>

**Notes.** We show each code's method of solving for the SPT kernels (see main text for more details), the computational speed (in order of magnitude in ms) and available modelling. 'EdS' and 'USA' subscripts indicate only the Einstein-de Sitter and Unscreened approximations are available for that model (see Sect. 3.6). Here 'DS' stands for Dark Scattering and 'WnW' stands for Wiggle-no-Wiggle (see Sect. 3.4). We note that MG-Copter only supports an EdS implementation of the WnW-EFTofLSS 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; Assassini 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

$$ \delta_{\text{halo}} = b_1 \delta + \frac{b_2}{2} \delta^2 + b_{G_2} \mathcal{G}_2(\tilde{\Phi}|\mathbf{x}) + b_{\Gamma_3} \Gamma_3 + b_{\nabla^2 \delta} \nabla^2 \delta + \text{noise}, \quad (38) $$

where $b_1$ and $b_2$ are the linear and quadratic bias parameters, $\nabla^2 \delta$ is a higher-derivative operator and 'noise' denotes the stochastic contributions uncorrelated with $\delta$ and with zero means, encoding shot-noise. The $\nabla^2 \delta$ 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). $\mathcal{G}_2$ and $\Gamma_3$ 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:

$$ \mathcal{G}_2(\tilde{\Phi}|\mathbf{x}) := [\nabla_{ij} \tilde{\Phi}(\mathbf{x})]^2 - [\nabla^2 \tilde{\Phi}(\mathbf{x})]^2, \quad (39) $$

and

$$ \Gamma_3(\mathbf{x}) := \mathcal{G}_2(\tilde{\Phi}|\mathbf{x}) - \mathcal{G}_2(\tilde{\Phi}_v|\mathbf{x}), \quad (40) $$

where the potentials $\tilde{\Phi}$ and $\tilde{\Phi}_v$ are defined as $\tilde{\Phi} = \nabla^{-2}\delta$, $\tilde{\Phi}_v = -\nabla^{-2}(\theta/f)$. We have also defined $\nabla_{ij} := \nabla_i \nabla_j$ and $\nabla^2 := \nabla_i \nabla_i$.

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

$$ P_g^s(k, \mu) = D_{\text{FoG}}(k^2 \mu^2 \sigma_v^2) [P_{g,1\text{-loop}}^s(k, \mu) + P_{g,\text{noise}}^s(k, \mu)] + P_{\text{ctr}}^s(k, \mu). \quad (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_{\text{FoG}}(k^2 \mu^2 \sigma_v^2) = 1$ and the counterterms $P_{\text{ctr}}^s(k, \mu)$ are given by Eqs. (32), (33) for PBJ and PyBird respectively. In TNS, $P_{\text{ctr}}^s(k, \mu) = 0$ while $D_{\text{FoG}}(k^2 \mu^2 \sigma_v^2)$ 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

$$ P_{g,\text{noise}}^s(k, \mu) = N + \frac{1}{\bar{n}} (\epsilon_0 k^2 + \epsilon_2 \mu^2 k^2), \quad (42) $$

where $\bar{n}$ is the galaxy number density in units of $h^3 \text{Mpc}^{-3}$ and $N, \epsilon_0, \epsilon_2$ are constants with dimensions $h^{-3} \text{Mpc}^3, h^{-2} \text{Mpc}^2$ and $h^{-2} \text{Mpc}^2$, respectively. Pybird uses the same modelling, with dimensionless parameters $c_{\epsilon,0} = N\bar{n}$, $c_{\epsilon,1} = \epsilon_0 k_M^2$ and $c_{\epsilon,2} = \epsilon_2 k_M^2/f$ to give

$$ P_{g,\text{noise}}^s(k, \mu) = \frac{1}{\bar{n}} \left( c_{\epsilon,0} + c_{\epsilon,1} \frac{k_M^2}{k_M^2} + c_{\epsilon,2} f \frac{\mu^2}{k_M^2} \right). \quad (43) $$

In this work we set $\epsilon_0 = \epsilon_2 = 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 $\epsilon_0$ does not significantly improve the fits to the simulations. $\epsilon_2$ 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 $\epsilon_2$. 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

$$ P_{g,1\text{-loop}}^s(k, \mu) = b_1^2 P_{1\text{-loop}}^{\delta\delta}(k) + 2b_1\mu^2 P_{1\text{-loop}}^{\delta\theta}(k) + \mu^4 P_{1\text{-loop}}^{\theta\theta}(k) \\ + A(k, \mu, b_1) + B(k, \mu, b_1) + C(k, \mu, b_1) \\ - k^2 \mu^2 \sigma_v^2 [b_1 F_1(k) + G_1(k) \mu^2]^2 P_{11,i}^{\delta\delta}(k) \\ + B_{g,22}^s(k, \mu) + B_{g,13}^s(k, \mu), \quad (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)

$$B_{g,22}^{s}(k, \mu) = 2 \int d^3 q P_{11}^{\delta\delta}(q) P_{11}^{\delta\delta}(|\mathbf{k} - \mathbf{q}|) \\ \times Z_{b,2}(\mathbf{q}, \mathbf{k} - \mathbf{q}) \left[ \frac{b_2}{2} + b_{G_2} S^b(\mathbf{q}, \mathbf{k} - \mathbf{q}) \right], \quad (45)$$

$$B_{g,13}^{s}(k, \mu) = 6 Z_{b,1}(\mu) P_{11}^{\delta\delta}(k) \int d^3 q Z_{b,3}^{\text{sym}}(\mathbf{q}, -\mathbf{q}, \mathbf{k}) P_{11}^{\delta\delta}(q), \quad (46)$$

where $S^b(\mathbf{k}_1, \mathbf{k}_2) = (\mathbf{k}_1 \cdot \mathbf{k}_2)^2/k_1^2 k_2^2 - 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,i}^{\delta\delta}(q)$. The $Z_{b,n}$ kernels can be written as (see Scoccimarro et al. 1999; Bernardeau et al. 2002, see)

$$Z_{b,1}(\mu) = b_1 + f \mu^2, \quad (47)$$

$$Z_{b,2}(\mathbf{k}_1, \mathbf{k}_2) = 2b_1 F_2^{\text{EdS}}(\mathbf{k}_1, \mathbf{k}_2) + 2f\mu^2 G_2^{\text{EdS}}(\mathbf{k}_1, \mathbf{k}_2) \\ + f \mu k \left[ \frac{\mu_1}{k_1}(b_1 + f\mu_2^2) + \frac{\mu_2}{k_2}(b_1 + f\mu_1^2) \right] \\ + \frac{b_2}{2} + b_{G_2} S^b(\mathbf{k}_1, \mathbf{k}_2), \quad (48)$$

$$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_{G_2} f \mu k \frac{\mu_1}{k_1} S^b(\mathbf{k}_2, \mathbf{k}_3) \\ + b_2 F_2^{\text{EdS}}(\mathbf{k}_1, \mathbf{k}_2) + 2b_{G_2} S^b(\mathbf{k}_1, \mathbf{k}_{23}) F_2^{\text{EdS}}(\mathbf{k}_2, \mathbf{k}_3) \\ + 2b_{\Gamma_3} S^b(\mathbf{k}_1, \mathbf{k}_{23}) \left[ F_2^{\text{EdS}}(\mathbf{k}_2, \mathbf{k}_3) - G_2^{\text{EdS}}(\mathbf{k}_2, \mathbf{k}_3) \right], \quad (49)$$

where ‘EdS’ means the kernel is an Einstein–de Sitter expression. We note that $f = -G_1(k, a)/F_1(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}^{\text{sym}}$ and variables $k$ and $\mu$ are the amplitude of the sum of the momenta $\mathbf{k} = \mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3$ and cosine of angle between $\hat{k}$ and the line of sight $\hat{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_1, b_2, b_{G_2}, b_{\Gamma_3}]$, while Pybird adopts the basis presented in D'Amico et al. (2020), that we indicate with $[\hat{b}_1, \hat{b}_2, \hat{b}_3, \hat{b}_4]$. There is a one-to-one correspondence between the two bases:

$$b_1 = \hat{b}_1, \quad (50)$$

$$b_2 = -2(\hat{b}_1 - \hat{b}_2 - \hat{b}_4), \quad (51)$$

$$b_{G_2} = \frac{2}{7} \left[ \hat{b}_2 + \frac{7}{4}(a_\gamma - 2)\hat{b}_1 \right], \quad (52)$$

$$b_{\Gamma_3} = \left( \frac{a_{\gamma}^2}{2} - \frac{a_{\gamma a}}{2} + 1 \right) \hat{b}_1 + \frac{2}{7} a_{\gamma} \hat{b}_2 - \frac{2}{21} \hat{b}_3. \quad (53)$$

Here $a_\gamma$ and $a_{\gamma a}$ are the bootstrap time-dependent functions introduced in D'Amico et al. (2021a) and are usually fixed to their EdS values $a_\gamma^{\text{EdS}} = 10/7$ and $a_{\gamma a}^{\text{EdS}} = 3/7$. We fix $b_{\Gamma_3} = 0$ because it is degenerate with $b_{G_2}$, as well as $b_{\gamma z_\delta} = 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_1, b_2, 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_1(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_1(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),

$$b_1(k) = b_0 \frac{1 + A_2 k^2}{1 + A_1 k}, \quad (54)$$

which was shown to match N-body halo measurements well. This introduces two new free parameters {$A_1, A_2$} with dimensions $h^{-1}$ Mpc and $h^{-2}$ 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_1$ 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):

$$P_{11}^{\delta\delta} = P_{11}^{w} + P_{11}^{nw}. \quad (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:

$$P_{11}^{nw}(k) = P_{EH}(k) [\mathcal{F} * R](k), \quad (56)$$

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

$$[\mathcal{F} * R](k) = \frac{1}{\sqrt{2\pi\lambda}} \int_{q_{\min}}^{q_{\max}} \frac{dq}{q} R(q) \exp \left\{ -\frac{(\ln k/q)^2}{2\lambda^2} \right\}, \quad (57)$$

with $\lambda$ being the dimensionless width of the Gaussian filter.

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

$$P_{\text{NL}}^{s,\text{IR-resum}}(k) = (b_1 + f \mu^2)^2 \left\{ P_{11}^{\text{nw}}(k) \right. \\ \left. \qquad + e^{-k^2 \Sigma^2(k, \ell_{\text{osc}})} P_{11}^{\text{w}}(k) [1 + k^2 \Sigma^2(k, \ell_{\text{osc}})] \right\} \\ + P_{1-\text{loop}}^{s,nw}(k) + e^{-k^2 \Sigma^2(k, \ell_{\text{osc}})} P_{1-\text{loop}}^{s,w}(k), \quad (58)$$

where $\Sigma^2(k, \ell_{\text{osc}})$ is the RSD damping function computed as

$$\Sigma^2(k, \ell) = \frac{1}{2} \left\{ \Xi_0(\ell) (1 + 2f\mu^2 + f^2\mu^2) \right. \\ \left. \qquad + \Xi_2(\ell) [(\hat{\mathbf{k}} \cdot \hat{\mathbf{l}})^2 + 2f\mu\mu_\ell (\hat{\mathbf{k}} \cdot \hat{\mathbf{l}}) + f^2\mu^2\mu_\ell^2] \right\}, \quad (59)$$

and we have defined

$$\Xi_0(\ell) = \frac{2}{3} \int_0^\infty \frac{dp}{2\pi^2} e^{-p^2/\Lambda_{\text{IR}}^2} P_{11}^{\delta\delta}(p) [1 - j_0(p\ell) - j_2(p\ell)], \quad (60)$$

$$\Xi_2(\ell) = 2 \int_0^\infty \frac{dp}{2\pi^2} e^{-p^2/\Lambda_{\text{IR}}^2} P_{11}^{\delta\delta}(p) j_2(p\ell). \quad (61)$$

Here $j_n$ are the spherical Bessel functions of order $n$ and $\mu_\ell$ is the angle between $\ell$ and the line of sight. $\Lambda_{\text{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 $\Lambda_{\text{IR}} = 0.2 \, \text{Mpc}^{-1}$. $\ell_{\text{osc}} = |\ell_{\text{osc}}|$ is set to be the baryon acoustic oscillations scale.

For PBJ, we adopt the following choice of parameters: $q_{\min} = k e^{-4\lambda}$, $q_{\max} = k e^{4\lambda}$, $k_s = 0.2 \, \text{Mpc}^{-1}$, $\lambda = 0.25$ and $\ell_{\text{osc}} = 102.707 \, h^{-1} \, \text{Mpc}$ with $\ell_{\text{osc}}$ and $\hat{k}$ having the same orientation. We note that $\ell_{\text{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 $\psi(y, t) := l(y, t) - \bar{y}$, where $l(y, t)$ is the final position of the particle and $\bar{y}$ is the initial position. The overdensity is given as

$$1 + \delta(\mathbf{x}, t) = \int d^3 y \delta_D^{(3)}[\mathbf{x} - \mathbf{y} - \psi(\mathbf{y}, t)] \\ = \int d^3 y \int \frac{d^3 k}{(2\pi)^3} e^{i k \cdot [\mathbf{x} - \mathbf{y} - \psi(\mathbf{y}, t)]}. \quad (62)$$

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

$$1 + \delta^s(\mathbf{s}, t) = \int d^3 y \delta_D^{(2)}[\mathbf{s}_{\perp} - \mathbf{y}_{\perp} - \psi_{\perp}^s(\mathbf{y}, t)] \\ \times \delta_D^{(1)}[\mathbf{s}_{\parallel} - \mathbf{y}_{\parallel} - \psi_{\parallel}^s(\mathbf{y}, t) - \psi_{\parallel}^s(\mathbf{y}, t)/H] \\ = \int d^3 y \int \frac{d^2 k_{\perp}}{(2\pi)^2} \int \frac{d k_{\parallel}}{2\pi} e^{i \mathbf{k}_{\perp} \cdot [\mathbf{s}_{\parallel} - \mathbf{y}_{\perp} - \psi_{\perp}^s(\mathbf{y}, t)]} \\ \times e^{i \mathbf{k}_{\parallel} [\mathbf{s}_{\parallel} - \mathbf{y}_{\parallel} - \psi_{\parallel}^s(\mathbf{y}, t) - \psi_{\parallel}^s(\mathbf{y}, t)/H]}. \quad (63)$$

Here $\psi_{\parallel}^s(\mathbf{y}, t) = (\psi^s(\mathbf{y}, t) \cdot \hat{z}) \hat{z}$ and $\psi_{\perp}^s(\mathbf{y}, t) = \psi^s(\mathbf{y}, t) - \psi_{\parallel}^s(\mathbf{y}, t)$. From Eq. (63), it is easy to derive the power spectrum

$$P(\mathbf{k}_{\perp}, \mathbf{k}_{\parallel}) = \int d^3 y e^{-i \mathbf{k} \cdot \mathbf{y}} K_r(\mathbf{k}, \mathbf{y}), \quad (64)$$

where

$$K_r(\mathbf{k}, \mathbf{y}) = \left\langle e^{i \mathbf{k}_{\perp} \cdot [\psi_{\perp}^{s}(\mathbf{y}) - \psi_{\parallel}^{s}(\mathbf{0})]} e^{i \mathbf{k}_{\parallel} \cdot \left[ \psi_{\parallel}^{s}(\mathbf{y}) - \psi_{\perp}^{s}(\mathbf{0}) + \frac{\psi_{\parallel}^{s}(\mathbf{y}, t) - \psi_{\perp}^{s}(\mathbf{0}, t)}{H} \right]} \right\rangle \\ = \exp \sum_{N=1}^{\infty} \frac{1}{N!} \left\langle \left\{ i \mathbf{k}_{\perp} \cdot \Delta_{\perp}(\mathbf{y}) + i \mathbf{k}_{\parallel} \cdot [\Delta_{\parallel}(\mathbf{y}) + \dot{\Delta}_{\parallel}(\mathbf{y})/H] \right\}^N \right\rangle, \quad (65)$$

having defined $\Delta(\mathbf{y}) := \psi(\mathbf{y}) - \psi(\mathbf{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

$$P^{\ell}(k)|_{N} = \sum_{j=0}^{N} \sum_{\ell'} 4\pi(-i)^{\ell-\ell'} \int_{0}^{\infty} dr r^{2} Q_{\|N-j}^{\ell\ell'}(k, r) \xi_{j}^{\ell'}(r), \quad (66)$$

where $P^{\ell}(k)|_N$ is the resummed ($\ell$-th multipole of the) power spectrum up to order $N$, $\xi_j^{\ell}(r)$ is the $j$-th loop order term in Eulerian perturbation theory of the ($\ell$-th multipole of the) correlation function, and $Q_{\|N-j}^{\ell\ell'}(k, r)$ is given by

$$Q_{\|N-j}^{\ell\ell'}(k, r) = i^{\ell'} \frac{2\ell + 1}{2} \times \int_{-1}^{1} d\mu_k \int \frac{d^2 \hat{r}}{4\pi} e^{-i\mathbf{k} \cdot \mathbf{r}} F_{\|N-j}(\mathbf{k}, \mathbf{r}) \mathcal{L}_{\ell}(\mu_k) \mathcal{L}_{\ell'}(\mu_r), \quad (67)$$

$$F_{\|N-j}(\mathbf{k}, \mathbf{r}) = T_{0,r}(\mathbf{k}, \mathbf{r}) T_{0,r}^{-1} \|_{N-j}(\mathbf{k}, \mathbf{r}), \quad (68)$$

$$T_{0,r}(\mathbf{k}, \mathbf{r}) = \exp \left\{ -k^2 \Sigma^2(\mathbf{k}, \mathbf{r}) \right\}, \quad (69)$$

where $d^2\hat{r}$ is an angular integration and $\mu = \mu_k$ is now an integration variable. Here $T_{0,r}^{-1}\|_{N-j}(\mathbf{k}, \mathbf{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)

$$-\frac{k^2}{3} \int \frac{d^3 q}{(2\pi)^3} e^{-q^2/\Lambda_{\text{IR}}^2} \frac{P_{11}^{\delta\delta}(q)}{q^2} [P(k) - P(|k - q|)], \quad (70)$$

where the power spectra inside square brackets are computed at a given order in perturbation theory. The cutoff ensures that $q \ll k$ in most of the integration domain. If $P(k)$ contains a new scale, as is the case for $\ell_{\text{osc}}$ in $P^w(k)$, the difference in parentheses is either of order $P_{11}^w(k) q^2 \ell_{\text{osc}}^2$ (for $q \ell_{\text{osc}} \ll 1$) or $P_{11}^w(k)$ (for $q \ell_{\text{osc}} \gg 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^2$ 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

$$[P(k) - P(|k - q|)] \approx -\frac{1}{2} \frac{d^2 P(k)}{dk^2} q^2 x^2, \quad (71)$$

with

$$x = \frac{k \cdot q}{k q}. \quad (72)$$

For an approximately scale-invariant power spectrum, the second derivative gives a contribution of order $P(k)/k^2$ and, consequently, Eq. (70) is not $k^2$-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 $\mu(k, a)$ of Eq. (18) is inherited by the power spectrum would imply a contribution to Eq. (71) proportional to

$$\frac{\partial^2 \ln \mu(k, a)}{\partial k^2} \ll \ell_{\text{osc}}^2, \quad (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}^{(cb)}(k, a)$. We do not include higher order massive neutrino effects in the kernels $F_n^{\text{app}}$ and $G_n^{\text{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

$$P_{11,i}^{(cb)}(k) = \frac{P_{11}^{(cb)}(k, a)}{F_1^{(m)2}(k, a)}, \quad (74)$$

where $F_1^{(m)}(k, a)$ assumes no massive neutrino contribution but the same total matter contribution to the background evolution and Poisson equation, $\Omega_m$. We have checked that using $\Omega_{cb}$ instead makes negligible difference on the final RSD power spectrum predictions. On the other hand, $P_{11}^{(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 $\Lambda$CDM and beyond-$\Lambda$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 MGCopter 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,

$$F_2(\mathbf{k}_1, \mathbf{k}_2) = \frac{F_1(\mathbf{k}_1) F_1(\mathbf{k}_2)}{F_1^{\text{app}}(\mathbf{k}_1) F_1^{\text{app}}(\mathbf{k}_2)} F_2^{\text{app}}(\mathbf{k}_1, \mathbf{k}_2), \\ G_2(\mathbf{k}_1, \mathbf{k}_2) = \frac{f(\mathbf{k}_{12}) F_1(\mathbf{k}_1) F_1(\mathbf{k}_2)}{f^{\text{app}}(\mathbf{k}_{12}) F_1^{\text{app}}(\mathbf{k}_1) F_1^{\text{app}}(\mathbf{k}_2)} G_2^{\text{app}}(\mathbf{k}_1, \mathbf{k}_2),$$
$$ 
\begin{aligned}
F_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) &= \frac{F_1(\mathbf{k}_1) F_1(\mathbf{k}_2) F_1(\mathbf{k}_3)}{F_1^{\text{app}}(\mathbf{k}_1) F_1^{\text{app}}(\mathbf{k}_2) F_1^{\text{app}}(\mathbf{k}_3)} F_3^{\text{app}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3), \\
G_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) &= \frac{f(\mathbf{k}_{123}) F_1(\mathbf{k}_1) F_1(\mathbf{k}_2) F_1(\mathbf{k}_3)}{f^{\text{app}}(\mathbf{k}_{123}) F_1^{\text{app}}(\mathbf{k}_1) F_1^{\text{app}}(\mathbf{k}_2) F_1^{\text{app}}(\mathbf{k}_3)} \\
&\quad \times G_3^{\text{app}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3),
\end{aligned}
\quad (75) $$

where we have dropped the time dependence for compactness and we remind the reader that $f(k) = -G_1(k)/F_1(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_1^{\text{app}}(k, a) = F_1(k, a)$ and $G_1^{\text{app}}(k, a) = G_1(k, a)$, but the higher-order kernel’s scale dependence are given by the standard Einstein-de Sitter universe ($\Omega_m(a) = 1$) expressions, $F_n^{\text{EdS}}$ and $G_n^{\text{EdS}}$ (see Bernardeau et al. 2002). The time dependence is then given by the appropriate factors of $F_1$ and $G_1$, for example $F_2^{\text{app}}(\mathbf{k}_1, \mathbf{k}_2) = F_1(\mathbf{k}_1, a) F_1(\mathbf{k}_2, a) F_2^{\text{EdS}}(\mathbf{k}_1, \mathbf{k}_2)$.

2. The unscreened approximation (USA): $F_1^{\text{app}}(k, a) = F_1(k, a)$ and $G_1^{\text{app}}(k, a) = G_1(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_1^{\text{app}}(k)$ and $G_1^{\text{app}}(k)$ are the ΛCDM growth factors which are scale-independent to a very good approximation: $F_1^{\text{ΛCDM}}(k, a) \approx F_1^{\text{ΛCDM}}(a)$ and $G_1^{\text{ΛCDM}}(k, a) \approx G_1^{\text{Λ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^{-6}$), F5 ($f_{R0} = -10^{-5}$), N5 ($H_0r_c/c = 5$) and N1 ($H_0r_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^3$ dark-matter particles in a cubical box of length $1024 h^{-1}$Mpc, giving a mass resolution of $m_p = 7.78 \times 10^{10} h^{-1} M_\odot$. 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),

$$ \{\Omega_{b,0}, \Omega_{m,0}, h, n_s, A_s\} = \{0.046, 0.281, 0.697, 0.971, 2.297 \times 10^{-9}\}, $$

where $A_s$ is the amplitude of primordial scalar perturbations, $n_s$ is the scalar spectral index and $\Omega_{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-2 N-body code (Springel 2005) which consistently implements the effects of the momentum exchange within the dark sector. The simulations consisted of $1024^3$ dark matter particles in a periodic cosmological box of length $1 h^{-1}$ Gpc, evolved from a starting redshift of $z_i = 99$. The resulting CDM particle mass is $m_c = 8 \times 10^{10} h^{-1} M_\odot$ and the spatial resolution is $\epsilon = 24 h^{-1}$ kpc. The cosmological parameters are

$$ \{\Omega_{b,0}, \Omega_{m,0}, h, n_s, A_s\} = \{0.048, 0.308, 0.678, 0.966, 2.115 \times 10^{-9}\}. $$

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” (DEMNUuni) 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 DEMNUuni 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 $\epsilon = 20 h^{-1}$ kpc, a comoving volume of $8 h^{-3}$ Gpc$^3$ filled with $2048^3$ dark matter particles and, when present, $2048^3$ neutrino particles. The simulations are initialised at $z_i = 99$ with Zel'dovich 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 DEMNUuni simulations were run using the tree
<table><caption>**Table 2.** Simulations and associated cosmological or gravitational models considered in this work.</caption><thead><tr><th colspan="3">ELEPHANT<br>V = 5 × 1 h<sup>-3</sup> Gpc<sup>3</sup></th><th colspan="4">DAKAR<br>V = 1 h<sup>-3</sup> Gpc<sup>3</sup></th><th colspan="4">DEMNUni<br>V = 8 h<sup>-3</sup> Gpc<sup>3</sup></th></tr><tr><th>Code</th><th>|f<sub>Ro</sub>|</th><th>Ω<sub>rc</sub></th><th>Code</th><th>w<sub>0</sub></th><th>w<sub>a</sub></th><th>ξ [bn GeV<sup>-1</sup>]</th><th>Code</th><th>w<sub>0</sub></th><th>w<sub>a</sub></th><th>Σ m<sub>v</sub> [eV]</th></tr></thead><tbody><tr><td>F5</td><td>10<sup>-5</sup></td><td>–</td><td>w09</td><td>-0.9</td><td>0</td><td>10</td><td>CPL3</td><td>-1.1</td><td>-0.3</td><td>0.00</td></tr><tr><td>F6</td><td>10<sup>-6</sup></td><td>–</td><td>w11</td><td>-1.1</td><td>0</td><td>10</td><td>CPL3-16</td><td>-1.1</td><td>-0.3</td><td>0.16</td></tr><tr><td>N1</td><td>–</td><td>0.25</td><td>CPL2</td><td>-1.1</td><td>0.3</td><td>50</td><td>CPL4</td><td>-0.9</td><td>0.3</td><td>0.00</td></tr><tr><td>N5</td><td>–</td><td>0.01</td><td>–</td><td>–</td><td>–</td><td>–</td><td>CPL4-16</td><td>-0.9</td><td>0.3</td><td>0.16</td></tr></tbody></table>

**Notes.** We show the volume of each simulation box, noting that for ELEPHANT we have five realisations. Further, each simulation set has a ΛCDM simulation with the same initial seeds as the beyond-ΛCDM simulations.

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)

$$ \{\Omega_{b,0}, \Omega_{m,0}, h, n_s, A_s\} = \{0.05, 0.32, 0.67, 0.96, 2.127 \times 10^{-9}\}. $$

Given these values, the reference (i.e., the massless neutrino case) CDM-particle mass resolution is $m_{\text{CDM}}^{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 $\Omega_{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 $\Omega_{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 $\chi^2$ 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 approxima-

tions 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 $\chi^2$ statistic is given by

$$ \chi^2(k_{\text{max}}) = \frac{1}{N_{\text{d.o.f.}}} \sum_{k=\text{min}}^{k_{\text{max}}} \sum_{\ell, \ell'=0,2} [P_{\ell,\text{data}}^s(k) - P_{\ell,\text{model}}^s(k)] \times \text{Cov}_{\ell, \ell'}^{-1}(k) [P_{\ell',\text{data}}^s(k) - P_{\ell',\text{model}}^s(k)], \quad (76) $$

where $P_{\ell}^s$ is the $\ell^{th}$ multipole of the redshift space CDM or halo power spectrum and $\text{Cov}_{\ell, \ell'}$ 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_{\text{dof}} = 2N_{\text{bins}} - N_{\text{params}}$, where $N_{\text{bins}}$ is the number of k-bins summed over and $N_{\text{params}}$ is the number of free parameters in the theoretical model. In this case, $N_{\text{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 \le 0.3 h \text{ Mpc}^{-1}$ 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_1$ 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 $\chi^2 \approx 1$ for $k_{\text{max}} \le 0.6 h \text{ Mpc}^{-1}$ 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 $\chi^2$ between ΛCDM and beyond-ΛCDM models is unchanged around $\chi^2 \sim 1$ with and without the shot noise contribution. This gives us confidence that the $\chi^2 \sim 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 $\bar{n} \approx 3.1 \times 10^{-4} h^3 \text{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 \le z \le 1.1) = 7.94 h^{-3} \text{Gpc}^3$ (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, $\bar{n} = 6.86 \times 10^{-4} h^3 \text{Mpc}^{-3}$ (Euclid Collaboration 2020). In what follows, this number density and $V = 8.8 h^{-3} \text{Gpc}^3$ 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_{\text{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 \text{Mpc}^{-1} \le k_{\text{max}} \le 0.300 h \text{Mpc}^{-1}$.

2. We then calculated the $1\sigma$ confidence intervals ($\Delta\chi^2$) on a $\chi^2$ distribution using $N_{\text{dof}}$ degrees of freedom.

3. The model validity scale $k_{\text{max}}$ was then chosen as the largest value of $k_{\text{max}}$ which has $|\chi^2(k_{\text{max}}) - \Delta\chi^2(k_{\text{max}})| \le 1$. A $\chi^2$ 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\sigma$ 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_1, \sigma_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 nar-

row 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<fnref n="8" />, 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, $\sigma_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_0$ and $P_2$ simulation measurements, where $\chi^2(k_{\text{max}}) \sim 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_0w_a$CDM, 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_0w_a$CDM (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, $\xi$, 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 $\xi \to 0$ we recover $w_0w_a$CDM and so we confirm that for mild $w_0w_a$CDM 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

<footnote n="8">⁸ Solutions to this issue are currently being investigated in the literature, including using Jeffreys priors (Donald-McCann et al. 2023) or profile likelihoods (Moretti et al. 2023; Holm et al. 2023).</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Ratio of the EdS approximation to the exact calculation of the redshift space matter power spectrum multipoles, computed with the TNS prescription, for the Dark Scattering model with equation of state of dark energy described by the CPL parametrization (see Eq. 11) with $w_0 = -1.1$ and $w_a = 0.3$ (denoted CPL2 in the text, see Table 2), and at redshift $z = 1$. Blue solid, orange dashed and green dotted lines denote the monopole, the quadrupole and the hexadecapole ratios, respectively. In the top-left, bottom-left and bottom-right panels, we use $\sigma_v = 5.2 h^{-1} \text{ Mpc}$ and $\xi = 50 \text{ bn GeV}^{-1}$, $\xi = 30 \text{ bn GeV}^{-1}$ and $\xi = 10 \text{ bn GeV}^{-1}$, respectively. In the top-right panel, we use $\sigma_v = 5.14 h^{-1} \text{ Mpc}$ in the EdS prediction, which is the value obtained by refitting to the exact calculation which uses $\sigma_v = 5.2 h^{-1} \text{ Mpc}$. Blue and beige bands indicate errors on the monopole and quadrupole assuming a $\Lambda$CDM Gaussian covariance with $V = 8.8 h^{-3} \text{ Gpc}^3$ and no shot noise contribution (we compare dark matter multipoles). We note that the hexadecapole error fills the plot and so we have omitted it.</figcaption></figure>

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 $\Lambda$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_0w_a$CDM 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 $\Lambda$CDM-scr approximation attempts to match the large and small scale behaviour of $\gamma_2$ and $\gamma_3$ 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 $\chi^2$ 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 $\Lambda$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 $\chi^2$ as a function of $k_{max}$, we only vary the TNS RSD nuisance parameter, $\sigma_v$, in order to minimise the $\chi^2$ 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 $\Lambda$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 $\Lambda$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 $\sigma_v$. The added damping of
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Ratio of the three approximations considered in this paper to the exact calculation of the redshift space matter power spectrum multipoles, computed with the TNS prescription, for two $f(R)$ models: $|f_{\text{R0}}| = 10^{-5}$ and $\sigma_v = 5.8 \, h^{-1} \, \text{Mpc}$ (left panels), and $|f_{\text{R0}}| = 10^{-6}$ and $\sigma_v = 5.6 \, h^{-1} \, \text{Mpc}$ (right panels), at redshift $z = 1$. The top, middle and bottom panels show the USA, the USA and EdS, and the ΛCDM-scr approximations, respectively. Blue and beige bands indicate errors on the monopole and quadrupole assuming a ΛCDM Gaussian covariance with $V = 8.8 \, h^{-1} \, \text{Gpc}^3$ and no shot noise contribution. We note that the hexadecapole error fills the plot and so we have omitted it.</figcaption></figure>

the quadrupole due to $f(R)$ increases the effectiveness of $\sigma_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 $\chi^2$ in the ELEPHANT simulations at small $k_{\text{max}}$. This is likely due to the very small scatter in these measurements as they are the averages of five realisations, yielding small $\chi^2$ 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 $\chi^2$ 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_1$ and $b_2$, 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, $\sigma_v$.

In the EFTofLSS analyses for biased tracers we adopted a model with three perturbative biases, $b_1, b_2$ and $b_{\tilde{G}_2}$, one shot-noise parameter, $N$, and two counterterms, $\tilde{c}_0$ and $\tilde{c}_2$. Adding more shot-noise or counterterm parameters gradually improves the reached $k_{\text{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 $\chi^2$ as a function of $k_{\text{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
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Reduced $\chi^2$ for the fit of the redshift space matter power spectrum multipoles, computed with the TNS prescription, to the one from numerical simulations, as a function of $k_{\text{max}}$ and at redshift $z = 1$. In the top panel we consider the DGP ELEPHANT simulations, in the middle panel the $f(R)$ ELEPHANT simulations, and in the bottom panel the Dark Scattering DAKAR simulations. Crosses indicate the exact computation while circles indicate the EdS and USA approximations. For the $f(R)$ models, the $\Lambda$CDM-scr approximation is shown as pentagons. The $\Lambda$CDM case is shown as grey squares. The error bars are the $1\sigma$ errors on the $\chi^2$ statistic with $N_{\text{d.o.f.}} = 2N_k - 1$ degrees of freedom, where $N_k$ are the number of wave modes used in calculating $\chi^2$. We subtract one because the TNS model has one degree of freedom.</figcaption></figure>

DGP case. Both levels of modification follow the $\Lambda$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 $\Lambda$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 $\chi^2$) than the $\Lambda$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 $\sigma_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 mul-

tipoles 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 $\Lambda$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 \, \text{Mpc}^{-1}$. 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*.
<table><caption>**Table 3.** Summary of TNS and EFTofLSS fits to the simulations.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Approx.</th><th colspan="4">TNS</th><th colspan="4">EFTofLSS</th></tr><tr><th>k<sub>max</sub></th><th>χ<sup>2</sup></th><th>{b<sub>1</sub>, b<sub>2</sub>, N, σ<sub>v</sub>}</th><th>k<sub>max</sub></th><th>χ<sup>2</sup></th><th>{b<sub>1</sub>, b<sub>2</sub>, b<sub>G2</sub>, N, c̃<sub>0</sub>, c̃<sub>2</sub>}</th></tr></thead><tbody><tr><td rowspan="2">ΛCDM</td><td>Exact</td><td>0.21</td><td>1.14</td><td>{2.22, 0.94, -627, 4.34}</td><td>0.18</td><td>1.04</td><td>{2.20, 0.11, 0.05, 733, 5.39, 33.6}</td></tr><tr><td>EdS</td><td>0.21</td><td>1.16</td><td>{2.22, 0.77, -535, 4.32}</td><td>0.18</td><td>1.08</td><td>{2.21, 0.12, 0.03, 654, 4.98, 39.92}</td></tr><tr><td rowspan="2">N5</td><td>Exact</td><td>0.19</td><td>1.12</td><td>{2.17, -0.32, 306, 4.20}</td><td>0.18</td><td>1.09</td><td>{2.19, -0.02, 0.01, 400, 1.47, 34.11}</td></tr><tr><td>EdS and USA</td><td>0.19</td><td>1.12</td><td>{2.17, -0.39, 365, 4.18}</td><td>0.18</td><td>1.13</td><td>{2.19, -0.02, -0.01, 309, 0.92, 40.57}</td></tr><tr><td rowspan="2">N1</td><td>Exact</td><td>0.22</td><td>1.17</td><td>{2.12, 0.56, -553, 4.84}</td><td>0.18</td><td>0.95</td><td>{2.12, -0.10, -0.03, 0, -2.21, 38.16}</td></tr><tr><td>EdS and USA</td><td>0.22</td><td>1.19</td><td>{2.12, 0.39, -443, 4.84}</td><td>0.18</td><td>0.98</td><td>{2.12, -0.01, -0.05, -297, -4.38, 46.82}</td></tr><tr><td rowspan="3">F6</td><td>Exact</td><td>0.21</td><td>1.12</td><td>{2.20, 1.00, -637, 4.27}</td><td>—</td><td>—</td><td>—</td></tr><tr><td>EdS and USA</td><td>0.21</td><td>1.15</td><td>{2.21, 0.96, -692, 4.49}</td><td>0.18</td><td>1.00</td><td>{2.19, -0.04, 0.04, 0.16, -0.64, -5.45}</td></tr><tr><td>ΛCDM-scr</td><td>0.21</td><td>1.12</td><td>{2.20, 0.83, -566, 4.33}</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="3">F5</td><td>Exact</td><td>0.26</td><td>0.99</td><td>{2.11, 1.07, -769, 4.62}</td><td>—</td><td>—</td><td>—</td></tr><tr><td>EdS and USA</td><td>0.26</td><td>1.04</td><td>{2.12, 1.18, -1030, 5.21}</td><td>0.21</td><td>1.18</td><td>{2.07, 0.37, -0.14, 0.17, -4.21, -7.54}</td></tr><tr><td>ΛCDM-scr</td><td>0.26</td><td>0.97</td><td>{2.10, 0.90, -672, 4.62}</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="2">ΛCDM</td><td>Exact</td><td>0.20</td><td>1.12</td><td>{2.12, -1.38, 1118, 3.12}</td><td>0.19</td><td>0.98</td><td>{2.17, -0.09, -0.41, 1666, -4.85, 32.91}</td></tr><tr><td>EdS</td><td>0.20</td><td>1.12</td><td>{2.12, -1.36, 1108, 3.12}</td><td>0.19</td><td>1.02</td><td>{2.17, -0.07, -0.43, 1592, -4.87, 39.18}</td></tr><tr><td rowspan="2">CPL3</td><td>Exact</td><td>0.19</td><td>0.97</td><td>{2.08, -0.92, 747, 3.53}</td><td>0.19</td><td>1.02</td><td>{2.16, -0.93, -0.82, 2200, -0.65, 13.17}</td></tr><tr><td>EdS</td><td>0.19</td><td>0.97</td><td>{2.08, -0.93, 751, 3.53}</td><td>0.19</td><td>1.07</td><td>{2.11, -0.22, -0.46, 1403, -5.86, 41.34}</td></tr><tr><td rowspan="2">CPL4</td><td>Exact</td><td>0.20</td><td>1.11</td><td>{2.31, 2.96, -603, 3.09}</td><td>0.19</td><td>1.02</td><td>{2.16, -0.32, -0.66, 1262, -4.53, 16.55}</td></tr><tr><td>EdS</td><td>0.20</td><td>1.11</td><td>{2.31, 3.00, -597, 3.06}</td><td>0.19</td><td>1.04</td><td>{2.31, 0.20, -0.54, 1720, -2.99, 39.77}</td></tr><tr><td rowspan="2">ΛCDM-16</td><td>Exact</td><td>0.22</td><td>1.10</td><td>{2.21, -1.95, 1326, 3.00}</td><td>0.21</td><td>1.11</td><td>{2.23, 0.98, -0.46, 1433, -2.65, 50.17}</td></tr><tr><td>EdS</td><td>0.22</td><td>1.10</td><td>{2.21, -1.95, 1326, 3.00}</td><td>0.20</td><td>1.09</td><td>{2.23, 0.58, -0.44, 1578, -3.03, 51.02}</td></tr><tr><td rowspan="2">CPL3-16</td><td>Exact</td><td>0.19</td><td>0.97</td><td>{2.20, -0.53, 522, 3.45}</td><td>0.19</td><td>1.00</td><td>{2.22, -1.83, -0.55, 3182, -14.14, -2.22}</td></tr><tr><td>EdS</td><td>0.19</td><td>0.85</td><td>{2.23, 1.75, -560, 3.61}</td><td>0.19</td><td>0.94</td><td>{2.18, 0.16, -0.35, 1366, -10.98, 45.87}</td></tr><tr><td rowspan="2">CPL4-16</td><td>Exact</td><td>0.22</td><td>1.07</td><td>{2.37, -1.94, 1622, 2.42}</td><td>0.19</td><td>1.10</td><td>{2.24, -0.31, -0.49, 1239, -15.73, 13.38}</td></tr><tr><td>EdS</td><td>0.22</td><td>1.10</td><td>{2.39, 3.82, -146, 2.61}</td><td>0.19</td><td>1.06</td><td>{2.38, 0.01, -0.34, 2283, -7.25, 35.05}</td></tr></tbody></table>

**Notes.** We show the approximations employed, the highest $k_{max}$ with its associated value of the reduced $\chi^2$ and best-fitting model parameter. The upper table shows the ELEPHANT fits, the lower table the DEMNUni fits. We note that $b_1$ and $N$ are normalised differently in the different models, so they are not expected to have similar values. $b_1$, $b_2$ and $b_{G2}$ are dimensionless. $N$ has dimensions $h^{-3}$ Mpc$^3$, $\sigma_v$ and $\{\tilde{c}_0, \tilde{c}_2\}$ have dimensions $h^{-1}$ Mpc and $h^{-2}$ Mpc$^2$ respectively.

<table><caption>**Table 4.** Q-bias TNS fits for $f(R)$ ELEPHANT models.</caption><thead><tr><th>Model</th><th>Approximation</th><th>k<sub>max</sub></th><th>χ<sup>2</sup></th><th>{b<sub>0</sub>, A<sub>1</sub>, A<sub>2</sub>, N, σ<sub>v</sub>}</th></tr></thead><tbody><tr><td rowspan="3">F6</td><td>Exact</td><td>0.23</td><td>1.10</td><td>{2.17, 0.02, -2.00, 912, 3.48}</td></tr><tr><td>EdS and USA</td><td>0.23</td><td>1.10</td><td>{2.18, 0.04, -2.00, 886, 3.71}</td></tr><tr><td>ΛCDM-scr</td><td>0.23</td><td>1.10</td><td>{2.18, 0.06, -2.00, 940, 3.52}</td></tr><tr><td rowspan="3">F5</td><td>Exact</td><td>0.26</td><td>0.96</td><td>{2.07, -0.12, -1.40, 322, 4.10}</td></tr><tr><td>EdS and USA</td><td>0.26</td><td>1.00</td><td>{2.08, -0.11, -1.15, 137, 4.74}</td></tr><tr><td>ΛCDM-scr</td><td>0.26</td><td>0.98</td><td>{2.07, -0.13, -1.04, 196, 4.26}</td></tr></tbody></table>

**Notes.** $A_1$ and $A_2$ have units $h^{-1}$ Mpc and $h^{-2}$ Mpc$^2$ respectively. $N$ has dimensions $h^{-3}$ Mpc$^3$ and $\sigma_v$ has units $h^{-1}$ Mpc.

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<fnref n="9" />, and attached to it a covariance computed with Euclid-like shot noise and volume ($V = 8.8 h^{-3}$ Gpc$^3$ and $\bar{n} = 6.86 \times 10^{-4} h^3$ 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 $\chi^2 \approx 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 ($\mu(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}| \le 10^{-5}$ 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_1$ and $\sigma_v$ are ~5% and ~10% consistent between all approximations and

<footnote n="9">⁹ We used the values of the nuisance parameters found at the $k_{max}$ given in Table 3.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Reduced $\chi^2$ for the fit of the redshift space halo power spectrum multipoles, computed with the TNS and the EFTofLSS prescriptions, to the measured multipoles from the DGP and $f(R)$ ELEPHANT simulations, as a function of $k_{\text{max}}$ and at redshift $z = 1$. From top to bottom, TNS fit to DGP, EFTofLSS fit to DGP, TNS fit to $f(R)$ with $|f_{R0}| = 10^{-5}$, TNS fit to $f(R)$ with $|f_{R0}| = 10^{-6}$, and EFTofLSS fit to $f(R)$. Crosses, circles and pentagons indicate the exact computation, the EdS and USA approximation and the ΛCDM-scr approximation, respectively. The ΛCDM cases are shown as black (TNS) and grey (EFTofLSS) squares. The error bars are the $1\sigma$ errors on the $\chi^2$ statistic with $N_{\text{def}} = 2N_k - N_\chi$ degrees of freedom, where $N_k$ are the number of wave modes used in calculating $\chi^2$. We use $N_\chi = 4$ for the TNS model using Eulerian bias and $N_\chi = 5$ for EFTofLSS and the TNS model with Q-bias. All EFTofLSS exact and $f(R)$ calculations are performed using PyBird, while the DGP EdS and USA calculations are performed using PBJ. The ΛCDM predictions are made using the exact kernel calculations.</figcaption></figure>

the exact results over the full range of $k_{\text{max}}$, with $b_1$ also being consistent with the simulation measurement. Interestingly, the best-fit values of $\sigma_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 $\sigma_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_2$ and $N$ also vary significantly, both as a function of $k_{\text{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_0w_a$ CDM and massive neutrinos

In Fig. 7 we show the $\chi^2$ as a function of $k_{\text{max}}$ for the DEMUni 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 ($\sum m_\nu = 0.16$ eV). We find that the EdS approximation does equally well to the exact kernel calculations for all cases, at least until the $\chi^2 \approx 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_0w_a$ CDM 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 $\chi^2$ over the full range of $k_{\text{max}}$ considered than EFTofLSS, but both RSD models have roughly a $k_{\text{max}}(\chi^2 \approx 1) \approx 0.19 \, h \, \text{Mpc}^{-1}$. This is similar to the $k_{\text{max}}$ for most of the ELEPHANT simulations (except F5), despite the DEMUni simulations having
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Ratio of the USA and EdS approximation to the exact calculation of the redshift space matter power spectrum multipoles, computed with the TNS prescription, for the $|f_{R0}| = 10^{-5}$ (F5) model at redshift $z = 1$. The models assume the fits as found in Table 3. The dotted lines show a USA and EdS prediction where the nuisance parameters are not refit to the simulation data, and take the same values as the exact computation. Blue and beige bands indicate errors on the monopole and quadrupole assuming a Euclid-like Gaussian covariance with $V = 8.8 h^{-3} \text{ Gpc}^3$ and $\bar{n} = 6.86 \times 10^{-4} h^3 \text{ Mpc}^{-3}$.</figcaption></figure>

a larger volume. This is likely compensated by the slightly lower $\sigma_8$ of these simulations, and hence slightly lower levels of non-linearity, leading to a better performance of SPT.

The volume of the DEMUni 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 $\sum m_\nu = 0.16 \text{ 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^{-3} \text{ Gpc}^3$ and $\bar{n} = 6.86 \times 10^{-4} h^3 \text{ 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 $\Lambda$CDM and $w_0w_a$CDM, follow the massive neutrino modelling fits almost perfectly, with $\chi^2 \approx 1$ at all $k_{\text{max}}$ considered. This further suggests that the effects of massive neutrinos with masses $\sum m_\nu \le 0.16 \text{ 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 $\sigma_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-$\Lambda$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\sigma$ bias is incurred on the growth rate, $f$, when screening is omitted at $z = 1$ for a survey of $V_s = 20 h^{-3} \text{Gpc}^3$, but not for $V_s = 10 h^{-3} \text{Gpc}^3$. 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 \approx 0.25 h \text{ Mpc}^{-1}$ and for $\xi$ 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 $\Lambda$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 $\chi^2$ fit to the simulation dark matter multipoles as the predictions without these approximations in DGP. The EdS approximation is generally excellent for DGP, $w_0w_a$CDM, and Dark Scattering models, both at the level of dark matter and
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as the TNS panels of Fig. 4 but fitting to a fiducial SPT-mock data vector using the exact prediction for the SPT kernels and the Eulerian bias (red crosses, and blue crosses for the N5 model). The model parameters are those found in Table 3. The mock data vector is given scatter using a Euclid-like Gaussian covariance, using $V = 8.8 h^{-3} \text{Gpc}^3$ and $\bar{n} = 6.86 \times 10^{-4} h^3 \text{Mpc}^{-3}$ and the same linear bias as the mock data. The fits are also made using the same covariance. We delimit the $k_{\text{max}}$ which was used to choose the best-fit bias and RSD nuisance parameters when comparing to the ELEPHANT simulations, roughly corresponding to the values found in Table 3. The black triangles indicate a modelling where no $f(R)$ effects are accounted for, i.e. a pure ΛCDM modelling.</figcaption></figure>

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 $\chi^2$ 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 $\chi^2$ 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_{\text{max}}$, included in the fit.

Similarly, not modelling any massive neutrino effects gives an equally good fit to the SPT-ΛCDM and $-w_0w_a$CDM mock data as the modelling including massive neutrinos, for $\sum m_\nu \le$

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_{\text{max}}$ such that $|\chi^2(k_{\text{max}}) - 1| \le 0.15$ when fit to the simulations, with the TNS yielding a slightly larger $k_{\text{max}}$ systematically over all beyond-ΛCDM scenarios$^{10}$. This range is roughly $0.19 h \text{Mpc}^{-1} \le k_{\text{max}} \le 0.21 h \text{Mpc}^{-1}$, with the exception being the strong $f(R)$ modification where we find the TNS model can fit the data well up to $k_{\text{max}} \approx 0.26 h \text{Mpc}^{-1}$. This is likely due to an enhanced efficiency of the phenomenological fingers-of-god RSD damping parameter $\sigma_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

<footnote n="10">10 Note that we have not included the EFTofLSS ~$\mu^2 k^2$ shot-noise term, which improves the fit.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Fig. 4 but for the DEMNUri simulation models. We show the cases without massive neutrinos (top two panels: TNS upper, EFTofLSS lower) and with massive neutrinos (bottom two panels: TNS upper, EFTofLSS lower).</figcaption></figure>

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 $\sigma_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_0w_a$CDM, 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
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 6 but for the $w_0w_a$ CDM DEMNUni cosmologies (see Table 2). The mock data vector in this case includes massive neutrinos ($\sum m_\nu = 0.16$ eV), and the black squares use the same modelling as the mock data. The crosses assume $\sum m_\nu = 0.0$ eV in the modelling.</figcaption></figure>

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 \le 0.5 h \text{ Mpc}^{-1}$ (see Appendix C).

**Acknowledgements.** The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the 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 Nederlandse Onderzoeksschool 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 allo-

cation 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

* Ahmed, S. N., Agharmim, B., Anthony, A. E., et al. 2004, Phys. Rev. Lett., 92, 181301

* Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Rev. Rel., 21, 2

* Archidiacono, M., Brinckmann, T., Lesgourgues, J., & Poulin, V. 2017, in Prospects in Neutrino Physics, eds. F. Di Lodovico, S. Pascoli, & P. Ballett

* Arnold, C., Fosalba, P., Springel, V., Puchwein, E., & Blot, L. 2019, MNRAS, 483, 790

* Arnold, C., Li, B., Giblin, B., Harnois-Déraps, J., & Cai, Y.-C. 2022, MNRAS, 515, 4161

* Assassi, V., Baumann, D., Green, D., & Zaldarriaga, M. 2014, JCAP, 08, 056

* Aviles, A., & Cervantes-Cota, J. L. 2017, Phys. Rev. D, 96, 123526

* Aviles, A., Valogiannis, G., Rodriguez-Meza, M. A., et al. 2021, JCAP, 04, 039

* Baldauf, T., Seljak, U., Desjacques, V., & McDonald, P. 2012, Phys. Rev. D, 86, 083540

* Baldauf, T., Mirbabayi, M., Simonović, M., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 043514

* Baldi, M., & Simpson, F. 2015, MNRAS, 449, 2239

* Baldi, M., & Simpson, F. 2017, MNRAS, 465, 653

* Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75

* Baratta, P., Bel, J., Gouyou Beauchamps, S., & Carbone, C. 2023, A&A, 673, A1

* Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, JCAP, 07, 051

* Bean, R., & Tangmatitham, M. 2010, Phys. Rev. D, 81, 083534

* Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109

* Bel, J., Pezzotta, A., Carbone, C., Sefusatti, E., & Guzzo, L. 2019, A&A, 622, A109

* Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, Phys. Rept., 367, 1
Beutler, F., Seo, H.-J., Saito, S., et al. 2017, MNRAS, 466, 2242

Blas, D., Garny, M., Ivanov, M. M., & Sibiryakov, S. 2016, JCAP, 07, 028

Bose, B., & Koyama, K. 2016, JCAP, 08, 032

Bose, B., & Koyama, K. 2017, JCAP, 08, 029

Bose, B., & Taruya, A. 2018, JCAP, 10, 019

Bose, B., Koyama, K., Hellwing, W. A., Zhao, G.-B., & Winther, H. A. 2017, Phys. Rev. D, 96, 023519

Bose, B., Baldi, M., & Pourtsidou, A. 2018a, JCAP, 04, 032

Bose, B., Koyama, K., Lewandowski, M., Vernizzi, F., & Winther, H. A. 2018b, JCAP, 04, 063

Bose, B., Koyama, K., & Winther, H. A. 2019, JCAP, 10, 021

Bose, B., Pourtsidou, A., Markovič, K., & Beutler, F. 2020a, MNRAS, 493, 5301

Bose, B., Cataneo, M., Tröster, T., et al. 2020b, MNRAS, 498, 4650

Bose, B., Wright, B. S., Cataneo, M., et al. 2021, MNRAS, 508, 2479

Bose, B., Tsedrik, M., Kennedy, J., et al. 2022, MNRAS, 519, 4780

Brax, P., Casas, S., Desmond, H., & Elder, B. 2021, Universe, 8, 11

Carbone, C., Petkova, M., & Dolag, K. 2016, JCAP, 07, 034

Carlson, J., White, M., & Padmanabhan, N. 2009, Phys. Rev. D, 80, 043531

Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, JHEP, 09, 082

Carrasco, J. J. M., Foreman, S., Green, D., & Senatore, L. 2014, JCAP, 07, 057

Carrilho, P., Moretti, C., Bose, B., Markovič, K., & Pourtsidou, A. 2021, JCAP, 10, 004

Carrilho, P., Moretti, C., & Pourtsidou, A. 2023, JCAP, 01, 028

Carrion, K., Carrilho, P., Spurio Mancini, A., Pourtsidou, A., & Hidalgo, J. C. 2024, Arxiv e-prints [arXiv:2402.18562]

Carroll, S. M., Duvvuri, V., Trodden, M., & Turner, M. S. 2004, Phys. Rev. D, 70, 043528

Castorina, E., Sefusatti, E., Sheth, R. K., Villaescusa-Navarro, F., & Viel, M. 2014, JCAP, 02, 049

Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, JCAP, 07, 043

Cautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195

Chan, K. C., & Blot, L. 2017, Phys. Rev. D, 96, 023528

Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509

Charmousis, C., Gregory, R., Kaloper, N., & Padilla, A. 2006, JHEP, 10, 066

Chevallier, M., & Polarski, D. 2001, IJMPD, 10, 213

Chudaykin, A., Ivanov, M. M., Philcox, O. H. E., & Simonovič, M. 2020, Phys. Rev. D, 102, 063533

Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, JCAP, 12, 012

Cuozzo, V., Carbone, C., Calabrese, M., Carella, E., & Migliaccio, M. 2024, JCAP, 04, 073

D'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 05, 005

D'Amico, G., Marinucci, M., Pietroni, M., & Vernizzi, F. 2021a, JCAP, 10, 069

D'Amico, G., Senatore, L., & Zhang, P. 2021b, JCAP, 01, 006

de la Bella, L. F., Regan, D., Seery, D., & Hotchkiss, S. 2017, JCAP, 11, 039

Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rept., 733, 1

Donald-McCann, J., Gsponer, R., Zhao, R., Koyama, K., & Beutler, F. 2023, MNRAS, 526, 3461

Dvali, G., Gabadadze, G., & Porrati, M. 2000, Phys. Rev. Lett., B485, 208

Eifler, T., Krause, E., Dodelson, S., et al. 2015, MNRAS, 454, 2451

Eisenstein, D. J., & Hu, W. 1997, ApJ, 511, 5

Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191

Euclid Collaboration (Mellier, Y., et al.) 2024, A&A, submitted, [arXiv:2405.13491]

Fabbian, G., Calabrese, M., & Carbone, C. 2018, JCAP, 02, 050

Fang, X., Blazek, J. A., McEwen, J. E., & Hirata, C. M. 2017, JCAP, 02, 030

Fujita, T., Mauuerhofer, V., Senatore, L., Vlah, Z., & Angulo, R. 2020, JCAP, 01, 009

Fukuda, Y., et al. 1998, Phys. Rev. Lett., 81, 1158, [Erratum: Phys. Rev. Lett. 81, 4279 (1998)]

García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987

García-Farieta, J. E., Hellwing, W. A., Gupta, S., & Bilicki, M. 2021, Phys. Rev. D, 103, 103524

Gorbunov, D., Koyama, K., & Sibiryakov, S. 2006, Phys. Rev. D, 73, 044016

Gouyou Beauchamps, S., Baratta, P., Escoffier, S., et al. 2023, A&A, submitted [arXiv:2306.05988]

Guidi, M., Veropalumbo, A., Branchini, E., Eggemeier, A., & Carbone, C. 2023, JCAP, 08, 066

Hagstotz, S., Gronke, M., Mota, D., & Baldi, M. 2019, A&A, 629, A46

Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro, R. 2020, JCAP, 03, 040

Hashimoto, I., Rasera, Y., & Taruya, A. 2017, Phys. Rev. D, 96, 043526

He, J.-H. 2013, Phys. Rev. D, 88, 103523

Heavens, A. F., Matarrese, S., & Verde, L. 1998, MNRAS, 301, 797

Heitmann, K., Bingham, D., Lawrence, E., et al. 2016, ApJ, 820, 108

Heitmann, K., Finkel, H., Pope, A., et al. 2019, ApJS, 245, 16

Hernández-Aguayo, C., Hou, J., Li, B., Baugh, C. M., & Sánchez, A. G. 2019, MNRAS, 485, 2194

Hernández-Molinero, B., Carbone, C., Jimenez, R., & Peña Garay, C. 2024, JCAP, 01, 006

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Holm, E. B., Herold, L., Simon, T., et al. 2023, Phys. Rev. D, 108, 123514

Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004

Ivanov, M. M., Simonovič, M., & Zaldarriaga, M. 2020, JCAP, 05, 042

Ivanov, M. M., Philcox, O. H. E., Nishimichi, T., et al. 2022, Phys. Rev. D, 105, 063512

Kacprzak, T., & Fluri, J. 2022, Phys. Rev. X, 12, 031029

Kaiser, N. 1987, MNRAS, 227, 1

Khoury, J., & Weltman, A. 2004, Phys. Rev. D, 69, 044026

Koyama, K., Taruya, A., & Hiramatsu, T. 2009, Phys. Rev. D, 79, 123512

Kreisch, C. D., Pisani, A., Carbone, C., et al. 2019, MNRAS, 488, 4413

Lacasa, F. 2022, A&A, 661, A70

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

Lesgourgues, J., & Pastor, S. 2006, Phys. Rept., 429, 307

Lewandowski, M., & Senatore, L. 2020, JCAP, 03, 018

Li, B., Zhao, G.-B., Teyssier, R., & Koyama, K. 2012, JCAP, 01, 051

Li, B., Zhao, G.-B., & Koyama, K. 2013, JCAP, 05, 023

Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301

Luty, M. A., Porrati, M., & Rattazzi, R. 2003, JHEP, 09, 029

Markovic, K., Bose, B., & Pourtsidou, A. 2019, Open J. Astrophys., 2, 13

Martinelli, M., Tutusaus, I., Archidiacono, M., et al. 2021, A&A, 649, A100

Marulli, F., Carbone, C., Viel, M., Moscardini, L., & Cimatti, A. 2011, MNRAS, 418, 346

Matsubara, T. 2008, Phys. Rev. D, 77, 063530

McDonald, P., & Roy, A. 2009, JCAP, 08, 020

McEwen, J. E., Fang, X., Hirata, C. M., & Blazek, J. A. 2016, JCAP, 09, 015

Mirbabayi, M., Schmidt, F., & Zaldarriaga, M. 2015, JCAP, 07, 030

Moresco, M., Marulli, F., Moscardini, L., et al. 2017, A&A, 604, A133

Moretti, C., Tsedrik, M., Carrilho, P., & Pourtsidou, A. 2023, JCAP, 12, 025

Nishimichi, T., Ohmuro, H., Nakamichi, M., et al. 2007, PASJ, 59, 1049

Nishimichi, T., D'Amico, G., Ivanov, M. M., et al. 2020, Phys. Rev. D, 102, 123541

Noriega, H. E., Aviles, A., Fromenteau, S., & Vargas-Magaña, M. 2022, JCAP, 11, 038

Oddo, A., Sefusatti, E., Porciani, C., Monaco, P., & Sánchez, A. G. 2020, JCAP, 03, 056

Oddo, A., Rizzo, F., Sefusatti, E., Porciani, C., & Monaco, P. 2021, JCAP, 11, 038

Parimbelli, G., Anselmi, S., Viel, M., et al. 2021, JCAP, 01, 009

Parimbelli, G., Carbone, C., Bel, J., et al. 2022, JCAP, 11, 041

Perko, A., Senatore, L., Jennings, E., & Wechsler, R. H. 2016, ArXiv e-prints [arXiv:1610.09321]

Perlmutter, S., et al. 1999, ApJ, 517, 565

Philcox, O. H. E., Ivanov, M. M., Cabass, G., et al. 2022, Phys. Rev. D, 106, 043530

Pietroni, M., Mangano, G., Saviano, N., & Viel, M. 2012, JCAP, 01, 019

Piga, L., Marinucci, M., D'Amico, G., et al. 2023, JCAP, 04, 038

Piras, D., & Lombriser, L. 2023, ArXiv e-prints [arXiv:2310.10717]

Planck Collaboration XVI. 2014, A&A, 571, A16

Pogosian, L., Silvestri, A., Koyama, K., & Zhao, G.-B. 2010, Phys. Rev. D, 81, 104023

Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2

Pourtsidou, A., & Tram, T. 2016, Phys. Rev. D, 94, 043518

Pourtsidou, A., Skordis, C., & Copeland, E. J. 2013, Phys. Rev. D, 88, 083505

Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179

Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Rizzo, F., Moretti, C., Pardede, K., et al. 2023, JCAP, 01, 031

Roncarelli, M., Carbone, C., & Moscardini, L. 2015, MNRAS, 447, 1761

Rossi, G., Choi, P. D., Moon, J., et al. 2021, MNRAS, 505, 377

Ruggeri, R., Castorina, E., Carbone, C., & Sefusatti, E. 2018, JCAP, 03, 003

Saito, S., Baldauf, T., Vlah, Z., et al. 2014, Phys. Rev. D, 90, 123522

Schmidt, F. 2009, Phys. Rev. D, 80, 123003

Schuster, N., Hamaus, N., Pisani, A., et al. 2019, JCAP, 12, 055

Scoccimarro, R., Couchman, H. M. P., & Frieman, J. A. 1999, ApJ, 517, 531

Senatore, L. 2015, JCAP, 11, 007

Senatore, L., & Zaldarriaga, M. 2014, ArXiv e-prints [arXiv:1409.1225]

Sheth, R. K., Chan, K. C., & Scoccimarro, R. 2013, Phys. Rev. D, 87, 083002

Silvestri, A., Pogosian, L., & Buniy, R. V. 2013, Phys. Rev. D, 87, 104015

Simon, T., Zhang, P., Poulin, V., & Smith, T. L. 2023, Phys. Rev. D, 107, 123530

Simpson, F. 2010, Phys. Rev. D, 82, 083505

Song, Y.-S., Taruya, A., Linder, E., et al. 2015, Phys. Rev. D, 92, 043522

Springel, V. 2005, MNRAS, 364, 1105
Taruya, A. 2016, Phys. Rev. D, 94, 023504

Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522

Taruya, A., Nishimichi, T., Bernardeau, F., Hiramatsu, T., & Koyama, K. 2014, Phys. Rev. D, 90, 123515

Tsedrik, M., Moretti, C., Carrilho, P., Rizzo, F., & Pourtsidou, A. 2023, MNRAS, 520, 2611

Vainshtein, A. 1972, Phys. Lett. B, 39, 393

Valogiannis, G., & Bean, R. 2019, Phys. Rev. D, 99, 063526

Valogiannis, G., Bean, R., & Aviles, A. 2020, JCAP, 01, 055

Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019, JCAP, 12, 040

Verza, G., Carbone, C., & Renzi, A. 2022, ApJ, 940, L16

Verza, G., Carbone, C., Pisani, A., & Renzi, A. 2023, JCAP, 12, 044

Viel, M., Haehnelt, M. G., & Springel, V. 2010, JCAP, 10, 015

Vielzeuf, P., Calabrese, M., Carbone, C., Fabbian, G., & Baccigalupi, C. 2023, JCAP, 08, 010

Villaescusa-Navarro, F., Marulli, F., Viel, M., et al. 2014, JCAP, 03, 011

Villaescusa-Navarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53

Vlah, Z., White, M., & Aviles, A. 2015, JCAP, 09, 014

Wadekar, D., Ivanov, M. M., & Scoccimarro, R. 2020, Phys. Rev. D, 102, 123521

Will, C. M. 2018, Theory and Experiment in Gravitational Physics (Cambridge University Press)

Wright, B. S., Koyama, K., Winther, H. A., & Zhao, G.-B. 2019, JCAP, 06, 040

Zennaro, M., Bel, J., Villaescusa-Navarro, F., et al. 2017, MNRAS, 466, 3244

Zennaro, M., Bel, J., Dossett, J., Carbone, C., & Guzzo, L. 2018, MNRAS, 477, 491

Zhang, P., D'Amico, G., Senatore, L., Zhao, C., & Cai, Y. 2022, JCAP, 02, 036

Zucca, A., Pogosian, L., Silvestri, A., & Zhao, G.-B. 2019, JCAP, 05, 001

23 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland

24 Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France

25 School of Mathematics and Physics, University of Surrey, Guildford, Surrey GU2 7XH, UK

26 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

27 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

28 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy

29 INFN-Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

30 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

31 INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy

32 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

33 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy

34 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica, Edificio G. Marconi 00185, Roma, Italy

35 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, (Barcelona), Spain

36 Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra, (Barcelona), Spain

37 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

38 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

39 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

40 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain

41 University of Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, 69622 Villeurbanne, France

42 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

43 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

44 UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France

45 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande 1749-016, Lisboa, Portugal

46 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande 1749-016, Lisboa, Portugal

47 Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

48 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy

49 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

50 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

51 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

52 Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain

53 Institut de Ciencies de l'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s/n Cerdanyola del Vallés, 08193 Barcelona, Spain

54 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy

55 INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

1 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

2 Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, Viale delle Scienze 7/A, 43124 Parma, Italy

3 Technion Israel Institute of Technology, Haifa, Israel

4 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, TS, Italy

5 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy

6 IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy

7 INFN Gruppo Collegato di Parma, Viale delle Scienze 7/A, 43124 Parma, Italy

8 Dipartimento di Fisica "Aldo Pontremoli", Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy

9 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy

10 School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK

11 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

12 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany

13 Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126 Napoli, Italy

14 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

15 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy

16 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

17 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, (TO), Italy

18 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK

19 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

20 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

21 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

22 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irne-rio 46, 40126 Bologna, Italy
56 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany

57 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy

58 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern 0315 Oslo, Norway

59 von Hoerner & Sulger GmbH, SchloßPlatz 8 68723, Schwetzingen, Germany

60 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

61 Cosmic Dawn Center (DAWN), Copenhagen, Denmark

62 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

63 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

64 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland

65 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

66 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

67 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-sur-Yvette, France

68 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

69 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland

70 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

71 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands

72 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

73 Dipartimento di Fisica e Astronomia “Augusto Righi” – Alma Mater Studiorum Università di Bologna, Via Piero Gobetti 93/2, 40129 Bologna, Italy

74 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK

75 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

76 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, 8000 Aarhus C, Denmark

77 Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

78 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

79 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

80 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France

81 Space Science Data Center, Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy

82 Centre National d'Etudes Spatiales – Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

83 Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov 077125, Romania

84 Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, Via Marzolo 8, 35131 Padova, Italy

85 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany

86 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile

87 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain

88 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain

89 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

90 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

91 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain

92 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

93 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

94 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy

95 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

96 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

97 Institut de Astrofísica de Canarias, Calle Vía Láctea s/n, 38204 San Cristóbal de La Laguna, Tenerife, Spain

98 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany

99 Université St Joseph; Faculty of Sciences, Beirut, Lebanon

100 Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France

101 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France

102 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, TS, Italy

103 Institut de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

104 CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA

105 Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France

106 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

107 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

108 Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France

109 Dipartimento di Fisica – Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

110 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA

111 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France

112 Institute Lorentz, Leiden University, PO Box 9506, Leiden 2300 RA, The Netherlands

113 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

114 Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

115 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada

116 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain

117 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

118 Department of Computer Science, Aalto University, PO Box 15400, Espoo 00 076, Finland

119 Department of Physics and Astronomy, Vesilinnantie 5, 20041 University of Turku, Finland

120 Serco for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

121 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia

122 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Victoria 3122, Australia
123 W.M. Keck Observatory, 65-1120 Mamalahoa Hwy, Kamuela, HI, USA

124 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa

125 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm SE-106 91, Sweden

126 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

127 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

128 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy

129 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal

130 Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany

131 Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, Roma, Italy

132 INFN, Sezione di Roma 2, Via della Ricerca Scientifica 1, Roma, Italy

133 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

134 Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

135 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA

136 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
## 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 here. Following Bose & Koyama (2016) we have

$$A(k, \mu) = -(k\mu) \int d^3 q \left[ \frac{\mu_q}{q} B_{\sigma}(\mathbf{q}, \mathbf{k} - \mathbf{q}, -\mathbf{k}) + \frac{k\mu - q\mu_q}{|\mathbf{k} - \mathbf{q}|^2} B_{\sigma}(\mathbf{k} - \mathbf{q}, \mathbf{q}, -\mathbf{k}) \right], \quad (A.1)$$

where $\mu_q = \hat{q} \cdot \hat{z}$ and the cross bispectrum is given in terms of the general kernels and linear bias as

$$\begin{aligned} B_{\sigma}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = & 2 \\ & \times \left\{ \begin{aligned}[t] & [b_1 F_1(\mathbf{k}_2) - \mu_{\mathbf{k}_2}^2 G_1(\mathbf{k}_2)] [b_1 F_1(\mathbf{k}_3) - \mu_{\mathbf{k}_3}^2 G_1(\mathbf{k}_3)] G_2(\mathbf{k}_2, \mathbf{k}_3) \\ & + [b_1 F_1(\mathbf{k}_3) - \mu_{\mathbf{k}_3}^2 G_1(\mathbf{k}_3)] [b_1 F_2(\mathbf{k}_1, \mathbf{k}_3) - \mu_{\mathbf{k}_2}^2 G_2(\mathbf{k}_1, \mathbf{k}_3)] G_1(\mathbf{k}_1) \\ & \times P_{11,i}(\mathbf{k}_2) P_{11,i}(\mathbf{k}_3) \\ & + [b_1 F_1(\mathbf{k}_2) - \mu_{\mathbf{k}_2}^2 G_1(\mathbf{k}_2)] [b_1 F_2(\mathbf{k}_1, \mathbf{k}_2) - \mu_{\mathbf{k}_3}^2 G_2(\mathbf{k}_1, \mathbf{k}_2)] G_1(\mathbf{k}_1) \\ & \times P_{11,i}(\mathbf{k}_1) P_{11,i}(\mathbf{k}_2) \end{aligned} \right\}. \end{aligned} \quad (A.2)$$

The B term is given as

$$B(k, \mu) = (k\mu)^2 \int d^3 q F(q) F(k - q), \quad (A.3)$$

with

$$F(q) = \frac{\mu_q}{q} G_1(q) [b_1 F_1(q) - \mu_q^2 G_1(q)] P_{11,i}(q). \quad (A.4)$$

Finally, the C term is given as

$$C(k, \mu) = (k\mu)^2 \int d^3 p \, d^3 q \, \delta_D(\mathbf{k} - \mathbf{p} - \mathbf{q}) P_{11,i}(p) P_{11,i}(q) \\ \times \frac{\mu_p^2}{p^2} G_1^2(p) [b_1 F_1(q) - \mu_q^2 G_1(q)]^2. \quad (A.5)$$

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

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

$$\mu(k, a) = 1 + \left(\frac{k}{a}\right)^2 \frac{1}{3\Pi(k, a)}, \quad (B.1)$$

$$\gamma_2^{\text{app}}(\mathbf{k}-\mathbf{q}, \mathbf{q}, a) = \gamma_2(\mathbf{k}, 0, a) \frac{\Pi(0, a_f)}{\Pi(q, a_f)} \left\{ \frac{\Pi(0, a_f)}{\Pi(q, a_f)} + \frac{[\Pi(k, a_f) - \Pi(0, a_f)][\Pi(q, a_f) - \Pi(0, a_f)]}{\Pi^2(q, a_f)} \right\},$$

$$\gamma_2^{\text{app}}(\mathbf{k}, -\mathbf{q}, a) = \gamma_2(\mathbf{k}, 0, a) \frac{\Pi(0, a_f)}{\Pi(q, a_f)},$$

$$\gamma_2^{\text{app}}(\mathbf{k}, \mathbf{q} - \mathbf{k}, a) = \frac{q^2}{k^2} \gamma_2^{\text{app}}(\mathbf{k} - \mathbf{q}, \mathbf{q}, a), \quad (B.2)$$

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

$$\begin{aligned} \gamma_3^{\text{app}}(\mathbf{q}, \mathbf{k}, -\mathbf{q}, a) &= \gamma_3^{\text{app}}(\mathbf{q}, -\mathbf{q}, \mathbf{k}, a) = \gamma_3(0, \mathbf{k}, 0, a) \frac{\Pi^2(0, a_f)}{\Pi^2(q, a_f)}, \\ \gamma_3^{\text{app}}(\mathbf{k}, \mathbf{q}, -\mathbf{q}, a) &= \gamma_3(\mathbf{k}, 0, 0, a) \frac{\Pi^2(0, a_f)}{\Pi^2(q, a_f)}. \end{aligned} \quad (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

$$a_f^2 \Pi^2(q, a_f) = q^2 + m^2(a_f), \quad (B.4)$$

where, from Eq. (21),

$$m^2(a_f) := \frac{a_f^2 \gamma^3(a_f)}{2f_0 (3\Omega_{\text{m},0} - 4)^2}, \quad (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 $\chi^2$ 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_{\text{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 \le 0.5 \, h \, \text{Mpc}^{-1}$ 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^{-3} \, \text{Gpc}^3$ at $z=0.5$, estimated using Laureijs et al. (2011).<fnref n="11" /> In the ΛCDM case, we find sub-0.5% agreement between MG-Copter, PBJ and an additional code by Atsushi Taruya (AT).<fnref n="12" /> 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.

<footnote n="11">¹¹ Note the difference in volume taken here and in the main text for a Euclid-like setup comes from the different redshifts considered.</footnote>

<footnote n="12">¹² Download this code: http://www2.yukawa.kyoto-u.ac.jp/~atsushi.taruya/cpt_pack.html.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. C.1.** The ratio of the SPT 1-loop monopole (blue solid), quadrupole (orange dashed) and hexadecapole (green dotted) between different codes for ΛCDM (left) and DGP (right) at z = 0.5 without infrared-resummation applied. The **MG-Copter** prediction is taken as the reference in all plots. The top left shows the prediction by Atsushi Taruya's code, the middle left and top right plots show the PBJ predictions and the bottom left and bottom right show the **Pybird** predictions. Atsushi Taruya's code and PBJ predictions are computed using the EdS approximation and without screening in the DGP case. We also adopt these approximations in their respective reference spectra from **MG-Copter**. **Pybird** predictions are exact and so are the respective **MG-Copter** reference spectra. The DGP parameter is set to Ω<sub>rc</sub> = 0.25. Blue and beige bands indicate errors on the monopole and quadrupole assuming a ΛCDM Gaussian covariance with V = 3.8 h<sup>-3</sup> Gpc<sup>3</sup> and no shot noise contribution. We note that the hexadecapole error fills the plot and so we have omitted it.</figcaption></figure>

In the bottom left panel we compare the exact kernel predic-
tions between **MG-Copter** and **Pybird**. We find sub-1% agree-
ment in this case, with the monopole being sub-0.5%. Interest-
ingly, 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 imple-
mentation in PBJ. Again we observe sub-0.5% agreement. For
the comparison with **Pybird** we use the exact time and screen-
ing predictions for both codes and observe sub-1% agreement
between the predictions, with some oscillations caused by inter-
polation.

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} \le
0.2 h \, \text{Mpc}^{-1}$).

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) decompo-
sition approach (see Sect. 3.4), both within the EdS approxima-
tion, 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 dif-
ferences, 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 esti-
mated measurement errors within the predictions validity range
and for all multipoles. As discussed in Sect. 3.4.3 the pres-
ence of the baryon acoustic oscillations scale, $\ell_{osc}$, provides
a clear criterion to identify the leading contributions in $k^2$ 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 resumma-
tion, 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 Lagrangian resumma-
tion may be degenerate with counterterms, which is left to be
seen. A full comparison between PBJ and Pybird at the pos-
terior level will be conducted in a forthcoming Euclid paper
(D'Amico et al., in prep.).
<figure><img src="image_11.png" /><figcaption>**Fig. C.2.** The ratio of the SPT 1-loop monopole (blue solid), quadrupole (orange dashed) and hexadecapole (green dotted) between different codes for ΛCDM at z = 0.5. The MG-Copter prediction is taken as the reference in all plots. The left panels show the ratios using the WnW or OPTIRESUM resummation. The top left panel shows the PBJ prediction and the bottom panel the Pybird prediction, both predictions being computed using the EdS approximation. The right panel shows the same as the bottom left plot but now with Pybird applying the full Lagrangian resummation method. Blue and beige bands indicate errors on the monopole and quadrupole assuming a ΛCDM Gaussian covariance with V = 8.8 h⁻³ Gpc³ and no shot noise contribution. We note that the hexadecapole error fills the left plot and so we have omitted it but include it in the right hand plot as a green band.</figcaption></figure>