fkPT: constraining scale-dependent modified gravity with the full-shape galaxy power spectrum

Modified gravity models with scale-dependent linear growth typically exhibit an enhancement in the power spectrum beyond a certain scale. The conventional methods for extracting cosmological information usually involve inferring modified gravity effects via Redshift Space Distortions (RSD), particularly through the time evolution of fσ 8. However, classical galaxy RSD clustering analyses encounter difficulties in accurately capturing the spectrum's enhanced power, which is better obtained from the broad-band power spectrum. In this sense, full-shape analyses aim to consider survey data using comprehensive and precise models of the whole power spectrum. Yet, a major challenge in this approach is the slow computation of non-linear loop integrals for scale-dependent modified gravity, precluding the estimation of cosmological parameters using Markov Chain Monte Carlo methods. Based on recent studies, in this work we develop a perturbation theory tailored for Modified Gravity, or analogous scenarios introducing additional scales, such as in the presence of massive neutrinos. Our approach only needs the calculation of the scale-dependent growth rate f(k,t) and the limit of the perturbative kernels at large scales. We called this approximate technique as fk-Perturbation Theory and implemented it into the code fkpt, capable of computing the redshift space galaxy power spectrum in a fraction of a second. We validate our modeling and code with the f(R) theory MG-GLAM and General Relativity NSeries sets of simulations. The code is available at https://github.com/alejandroaviles/fkpt.


Introduction
During the last two decades galaxy surveys such as 2dF [1], WiggleZ [2], BOSS [3], eBOSS [4], and DES [5] have contributed to the inference of cosmological parameters and the test of models using techniques complementary to other probes, including the Cosmic Microwave Background (CMB), Supernovae type Ia, among others.These analyses are expected to be even more important in the near future with the upcoming stage IV experiments, such as the Dark Energy Spectroscopic Instrument (DESI) [6], Euclid [7] and the Legacy Survey of Space of Time (LSST) of the Vera C. Rubin Observatory [8], which are foreseen to achieve subpercentage accuracy in various cosmological parameters.Therefore, models for the analysis are demanded to be sufficiently precise.
The above-mentioned surveys can be employed to test gravity.As it is known, General Relativity (GR) is well tested in local, planetary scales using for instance Parametrized Post Newtonian (PPN) parameters [63], and several other approaches in astrophysical scales [64].This seems to be the case also in Cosmology, since the standard ΛCDM model is in general consistent with current cosmological data, and GR is believed to be correct at large scales.Most of the gravity theory tests so far have employed RSD, since this effect is intimately related to the growth rate of matter perturbations that depends upon the gravity model.RSD has been developed for many years, from the pioneering work of Kaiser [65] to modern nonlinear perturbation theory, see e.g.[55,[66][67][68].However, with the advent of more precise Large Scale Structure (LSS) data, new opportunities to test gravity at cosmological scales arise, yet the analysis demands methods of non-linear PT, such as the full-shape technique.This is because one of the main potential observables of MG theories is an enhancement of the power spectrum above a certain scale due to a larger strength of gravity, hence a major impact over the power spectrum is expected to lie in its broad-band shape.Let us shortly discuss about it.The sound horizon at decoupling is practically the same for ΛCDM and viable MG models, since they only affect late-time physics, and hence the BAO scale is unaltered, with the exception of a small extra degradation due to differences in the large scale displacement fields [69].On the other hand, RSD analysis measures the amplitude of fluctuations f σ 8 by breaking the degeneracy of the linear bias with the relative size between the different multipoles of the power spectrum.It is for this reason that the standard analysis (BAO + RSD) alone, we claim, is not efficient for detecting an MG signal, but one still needs the information carried by the broad-band of the power spectrum, that is precisely provided by the full-shape analysis.
The PT of MG has been explored extensively over the last decade [69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84].Recently, a few of us [85] proposed an accurate theoretical PT/EFT for MG models having a scale-dependent linear growth.This model was tested against the Elephant suite of N-body simulations for ΛCDM and the Hu-Sawicki (HS) f (R) model [86], finding very good agreement up to k ∼ 0.2 h Mpc −1 at redshift z = 0.5.Unfortunately, the computation of non-linear corrections to the power spectrum within this methodology is quite slow since the perturbative kernels are not known analytically, as for the ΛCDM case.Instead, they should be obtained from a set of differential equations that not only depend on the wave vectors configuration but also on the cosmological parameters.That is, to find the one-loop corrections to the power spectrum, having a form similar to one must solve differential equations to find the kernels K at each volume element of the integration.This complexity hinders the use of a Monte Carlo Markov Chain (MCMC) analysis for parameter estimation.A similar challenge arises in cosmological scenarios involving massive neutrinos, where scale-dependence is introduced through free-streaming and has been addressed in [37,38].In this work we build upon the ideas of these references and develop a method that is able to obtain the power spectrum in a fraction of a second.To do this, we first identify three effects that modify the standard Einstein-de Sitter (EdS) kernels: i) The failure of the velocity field θ = ∇ • v/(aHf ), to be equal to the overdensity δ at the linear perturbation level because the growth rate f at a given time is no longer constant, as in ΛCDM, but becomes scale-dependent; see fig. 1 below.The key relation is given by eq.(3.44), which also enters into the non-linear kernels.
ii) The abundance of matter Ω m cannot be approximated by the square of the growth rate, f 2 .This is particularly true for gravity models with scale-dependent linear growth, but may also be the case for some only time-dependent modified gravity models.
iii) The presence of non-linear screenings that should drive the MG theories at small cosmological scales.
The effects ii) and iii) give rise to the differential equations that considerably slow the computation of the one-loop power spectrum.Fortunately, we realize that the breakdown of the approximation Ω m ̸ = f 2 is not very harmful and can be healed by taking only timedependent corrections to the EdS kernels, as it is shown below in fig. 2. Furthermore, the scaledependent screenings operate on scales where the theoretical power spectrum is dominated by the EFT counterterms and become widely degenerate with them.We show that this is the case for HS-f (R) gravity model,5 but as long the screenings are small at mildly non-linear scales and below, this can be a good approximation for any MG theory.Although the validity of this assumption is expected to break above some wave-number, and should be preferably checked case by case.Finally, in §4 we notice that the dominant contribution is due to the effect described in point i).This allows us to construct a method that first evaluates f (k, t) using MG linear theory and then the perturbative kernels are evaluated ignoring the effects of ii) and iii).This method was proposed initially for massive neutrinos in [37] and further developed in [38] where the Python code Folps-nu was released.In the present work we test this methodology for the HS-f (R) theory, particularly we validate it against the state-ofthe-art MG-GLAM simulations [87,88].Once we ensure that our method is able to detect the MG signal, we test it against high precision NSeries simulations, which exists only for ΛCDM and were widely used in the past to test the pipeline of BOSS survey [10].
Together with this paper we release the C language code fkpt, 6 that computes the redshift space one-loop power spectrum of MG theories in a fraction of a second, including biasing terms, EFT counterterms, shot noise, as well as IR-resummations.Contrary to the PT codes enlisted above, fkpt does not use FFTLog methods since we want a flexible method that allows for a future incorporation of theories that do not reduce to ΛCDM at very large scales, and hence their kernels cannot be approximated as EdS when they are evaluated at small wave vectors.This, for example, is done for the scale-independent nDGP gravity model in [45].As explained in [37], the large scale kernels modifications almost double the number of matrix multiplications that should be done if using an FFTLog method.Hence, we adopt here a brute force approach, but equally accurate and almost as fast as FFTLog based codes, for solving the loop integrals.We expect our perturbative method and our code can be useful for testing gravity with the full-shape galaxy power spectrum of future spectroscopic surveys, particularly for DESI and Euclid.
The rest of this work is organized as follows: In §2 we discuss the basic equations of the f (R) model that are relevant to us.In §3 we delve into the general non-linear theory for models with scale-dependent linear growth, putting special emphasis on MG.In §4 we present the relevant approximations that lead to the fk-Perturbation Theory as well as our code fkpt.In §5 we validate our methodology using the MG-GLAM and NSeries simulations.Some of the formulas are displayed in a separate appendix.Final remarks and conclusions are put forward in §6.

Modified gravity
Although the perturbation theory method presented throughout this work is very general and encompasses a large variety of MG theories, we will concentrate on the case of f (R) gravity [86,89], which is perhaps the most studied scale-dependent MG model in Cosmology and for which we have the most precise simulations.It is defined through the action where R is the Ricci scalar and a function f (R) is added to the Einstein-Hilbert GR action.Further, S m the action of matter fields ψ i and depends also on the metric g.Variations with respect to the metric of this action lead to the field equations where f R = ∂f (R)/∂R, and with energy momentum tensor We use the Friedmann-Robertson-Walker metric with additional scalar perturbations in Newtonian gauge Equation (2.5) is the Poisson equation, modified by the term δf R , which then acts as a fifth force.We have defined such that eq.(2.6) governs the evolution of the scalar field, a.k.a. the Klein-Gordon equation, from which we can read that the mass of scalar field m = (M 1 /3) 1/2 , such that one recovers the ΛCDM model for scales larger than m −1 .Using eq.(2.6) to eliminate δf R from eq. (2.5), we arrive at with and δI. (2.12) As k ≫ m, the function µ goes to 4/3, and hence, in f (R) theories the strength of gravity is enhanced by a factor of one third at the smallest scales.This, of course may rule out any f (R) theory.However the nonlinear interactions, collected in the term δI provide a screening mechanism responsible to drive the theory to GR at specific limits.The screening mechanism behind f (R) and other scalar-tensor theories is the chameleon [90,91], that turns off the fifth force in regions with sufficiently deep potentials.In perturbation theory one expands the non-linear self-interaction δI as 7 δI(δf R ) = 1 2 where the functions M i are in general scale and time dependent.For the f (R) theories, these are only time dependent and given by the coefficients of the Taylor expansion Throughout this work we apply our results to the specific case of HS-f (R) model [86], with n = 1, defined by In order to have an effective ΛCDM model at background level, the energy scale is chosen to be M 2 = Ω m H 2 0 and the remaining constants comply with M 2 (a) = 9 4 M 3 (a) = 45 8 (Ω m a −3 + 4Ω Λ ) 7  (Ω m + 4Ω Λ ) 6 . (2.19) These functions depend only on the background evolution since they are the coefficients of the expansion of a scalar field potential about its background value, see eq. (2.15).The second order source entering eq.(3.33) becomes where Π(k) ≡ 1 3a 2 (3k 2 + M 1 a 2 ).The strength of the MG is given by the amplitude f R0 , the smaller is the value f R0 , the smaller the effect of MG.As a matter of notation the models with f R0 = −10 −6 , −10 −5 are called F6 and F5, respectively.While f R0 = 0 corresponds to GR.
Finally, we stress that results for other gravity models are straightforward to develop following the methods of the present and previous works; see e.g.Appendix B of [78].Moreover, theories defined in the Einstein frame can be put also within our framework by using field redefinitions [83]. 7We use the shorthand notation

Theoretical framework and perturvative kernels
In this section we present the general perturbative framework for theories with scale-dependent linear growth.To do this we first find the Lagrangian Perturbation Theory (LPT) kernels, afterwards we map them into the Eulerian frame, which is where we compute the multipoles of the power spectrum.
The LPT for MG was developed in [69].More recently, ref. [92] presented an LPT scheme for studying the clustering of dark matter particles in the presence of massive neutrinos.This framework is very general and can be applied for cosmologies with additional scales, and as such, under a few amendments it can be rewritten as a theory for MG, reducing to that of [69].In §3.1- §3.4 we present such scheme up to second order, since this is sufficient to understand the whole development idea, leaving the third order final results to Appendix A. After that, in §3.5 we show how to perform a mapping to obtain the kernels in Standard Perturbation Theory (SPT), that we use to obtain the power spectrum non-linear corrections.

Evolution in Lagrangian space
In the Lagrangian framework there is a particular interest in the map between the initial (Lagrangian) positions q of the cold dark matter particles and the final, or moment of observation, Eulerian positions x.This map between frames is given by the Lagrangian displacement vector field Ψ x(q, t) = q + Ψ(q, t).
The displacement evolves according to the Geodesic equation, where the first equality serves to define the differential operator T [93], and ∇ x = ∂/∂x are the derivatives with respect to the Eulerian coordinates.The gravitational potential Φ obeys the Poisson equation which here is modified by the term S(x, t).For example, in theories in the presence of massive neutrinos one identifies δ = f cb δ cb and S(x, t) = 4πGρ m f ν δ ν [92], where the label cb refers to the combined baryons and cold dark matter fluid and the label ν to the neutrino component.Furthermore, the relative abundances are defined as f cb = Ω cb /Ω m and f ν = Ω ν /Ω m .Also, in this case Ψ → f cb Ψ, so the Lagrangian displacement follows the cb particles.
In several MG theories one can write δ = δ m and S = 1 a 2 ∇ 2 x ϕ, where ϕ is an extra degree of freedom.We will follow this route below.
Taking the divergence of eq.(3.2) and moving to Fourier space we write where we omitted to write the time dependencies explicitly.Here, That is, for a function f (q) of the Lagrangian coordinates q, the notation means Further, for notational convenience, we write a tilde over Fourier transforms as δ(k) and S(k) to point out that they are "q-Fourier" transforms of functions defined over the Eulerian coordinates.That is, for f (x) we have Finally, we note that the "standard-Fourier" transform can be computed in any frame as Now, the equation of motion (3.4) contains functions that take Eulerian coordinates as arguments and functions that take Lagrangian coordinates.They are related by the Taylor expansion whose q−Fourier transform yields and the inverse relation [92] We will use the above equations to write eq.(3.4) in terms of Lagrangian coordinates only.But before that, in the next subsection we will revisit the standard approach used in ΛCDM, where the linear theory is scale independent.

Standard approach in vanilla ΛCDM
Before delving into more complex scenarios involving additional scales that arise when we incorporate the source term in eq.(3.3), let us revisit the case of the vanilla ΛCDM model.To achieve this, we take the divergence of the equation of motion (3.4) and utilize the Poisson equation (eq.(3.3) with S = 0), yielding Our objective is to write this equation in terms of pure Lagrangian coordinates, so that we can subsequently expand it perturbatively in terms of the displacement Ψ.Using the Jacobian transformation matrix, where a comma means partial derivative with respect to Lagrangian coordinates, we obtain the transformation rule between spatial derivatives with (J −1 ) ji the inverse of the Jacobian matrix.The LHS of (3.13) becomes up to cubic terms in the Lagrangian displacement.Now, for the RHS of eq. ( 3.13) we use the matter conservation where J = det J ij is the determinant of the Jacobian matrix, and the evolution is thought to be initiated sufficient early so we can neglect the Lagrangian fluctuation, 1 + δ(q) ≈ 1.
Then, one is able to relate the density field with the Lagrangian displacements which is again written in pure Lagrangian coordinates.Now, we can put all the ingredients together and write where we defined We notice eq.(3.19) is not the standard form to write the equation of motion, e.g.see [93].Indeed, it is notational simpler to omit the symbols [ • • • ](k).However, it is better to keep the above form for later comparison when including the source S(x) in the next subsection.
Because of the conservation equation (3.17), we had no need to use eq.(3.9) for the density field δ(k).This happens obviously for the vanilla ΛCDM model only, and hence the fact that one uses Fourier transforms in different frames is rarely noticed in the literature.However, below we will be forced to use eq.(3.9) when we include the source S.

Source
Now, let us consider the source term in the Poisson equation.We assume that S(k) has an expansion in the density fluctuations, that in Fourier space takes the form where A 0 α(k) is the factor of proportionality of the linear term and S NL (k, t) encapsulates the non-linearities.
Adding the source to the density field as appearing in the Poisson equation, where we performed the x-Fourier transform and in the last line we defined the function Now, based on eq.(3.9), the function α(k) in Eulerian space differs from its counterpart in Lagrangian space at first order in the relevant fields, in this case the Lagrangian displacement.Hence the combinations α(k)δ(k) and α(k) δ(k) differ only on non-linear terms in the density fluctuation.We refer to these non-linear terms as Frame-Lagging (FL). 8This property will allow us to obtain the function α(k) in Eulerian space where the physical phenomena are commonly described.For example, in massive neutrino cosmologies it is constructed with transfer functions coming out from a Boltzmann code [92].Thus, an important observation when transitioning to Lagrangian space is that we can retain the functional form of α(k) obtained in Eulerian space and add the corrections through FL terms.This process can be accounted for as follows.
For a Lagrangian space framework, we take the q-Fourier Transform (FT) We now look for an expression that schematically is written as Now, we omit for the moment the non-linear source SNL , which typically carries the screening effect, and develop up to second order,

S(k)
In the first equality we have expanded S(k) into S(k) using eq.(3.9).Then, in the second equality, we have written S/A 0 = αδ = α δ + α(δ − δ).In the last equality we have used eq.(3.12) and wrote δ , and finally grouped second order terms into the integral.Notice also that the term (δ − δ)Ψ i is third order and has been dropped out of this calculation.
The equation above still retains Eulerian coordinate dependencies through δ, which we can eliminate in favor of the Lagrangian displacement using eq.(3.18), leading to with frame-lagging kernels where we used δ 18).Equation (3.24) give the source in the Poisson equation (3.3) transformed to the Lagrangian frame, with the peculiarity of having the function α(k) the same functional form as in Eulerian coordinates.

Lagrangian displacement evolution equation
Finally, putting the above ingredients together, the geodesic equation becomes which is the equation of motion written in terms of pure Lagrangian coordinates.
To solve this equation perturbatively, we express the Lagrangian displacements as a series expansion Each term in the expansion can be further expanded as Here, δ (1) represents the first-order linear overdensity in Eulerian space.The integral involves a specific configuration of wavevectors is the kernel associated with the Lagrangian displacement at the n-th order.It is important to note that the kernel may have a time dependence, which is not considered in the standard approach that assumes GR and the relation Ω m = f 2 .

Linear order
After linearizing, eq.(3.27) can be expressed as T − A(k) [Ψ i,i ](k) = 0.This equation bears resemblance to the equation for the ΛCDM linear growth function D + , but modulated by the k dependence.Additionally, from eq. (3.18) one can confirm the relation between linear density and displacement fields is the same as in ΛCDM: where the linear growth function D + is the fastest growing solution to which we normalize to unity at large scales and present time: Therefore, the modification to the linear displacement arises solely from the extra kdependence introduced by the linear growth function, which is also present in the linear density field.As a result, the first-order Lagrangian Perturbation Theory (LPT) kernel can be expressed as as in the ΛCDM model.

2LPT
Keeping up to second order displacements in eq.(3.27), we obtain the second order kernel [69, 92] The functions A and B are obtained by first solving the equations subject to appropriate initial conditions to isolate the pure second order contributions, and ) is scale independent, both functions are equal A = B.In EdS background evolution we recover A(t) = B(t) = 1; whereas for the ΛCDM model, A(t) = B(t) are only weakly dependent on time and close to one.For standard cosmologies one finds that nowadays A ΛCDM (t 0 ) ≃ 1.01.
The term S 2 (k 1 , k 2 ) in Eq. (3.33) arises from the non-linear source term S NL and is responsible for the screening mechanism to second order in perturbation theory.For the HS-f (R) model this is given by with M 2 given by eq.(2.18).As will be discussed in sec.4.2, these terms are degenerate with EFT counterterms and will not be considered in this work.An alternative scenario that considers the screenings, and promises to be fast, was recently given in [46].But, we emphasize that one of the objectives of this paper is to show that screenings can be avoided in some cases because of their degeneracies with EFT parameters.Now, in one-loop integrals only the longitudinal piece of the displacement field shows up [94,95], so we can define the kernels Γ n through [82] where C n are constant coefficients that we fixed to the values C 1 = C 3 = 1, and C 2 = 3/7.Γ n and L (n) i kernels are related by A rapid computation yields the kernels Γ 1 (k) = 1 to first order, and to second order.The third order kernels can be found in Appendix A.

Velocity fields
For longitudinal flows, only the divergence of the velocity field is significant, hence it is useful to define where the peculiar velocity of particles relative to the Hubble flow is given in terms of the Lagrangian displacement time derivative by v = a Ψ.We also defined f 0 as the large scale limit of the growth rate, that is, where In the EdS case, we obtain Γ f n = Γ n , recovering the well-known relation Ψ(n) = nHf Ψ (n) [94].At first order, we have Γ f 1 (k) = f (k)/f 0 , which leads to the linear relationship between the velocity and density fields This relation will play a central role in the subsequent analysis.

Map to SPT Kernels
An alternative and more widely used path to cosmological Large Scale Structure Perturbation Theory is given by the Eulerian SPT formalism [15].Here, one expands the fields θ = θ (1) + θ (2) We can use eqs.(3.18), (3.37) and (3.43) to relate Lagrangian and Eulerian kernels.Following [84,85,96] we obtain to second order For third order kernels the relation is more cumbersome.But it simplifies for the particular configuration of wavevectors used in one-loop integrals

Power spectrum
The apparent position s of an object is distorted from its true position x because of the Doppler effect induced by its peculiar velocity, such that we observe it at a redshift space coordinate Since the map from real to redshift coordinates conserves the number of tracers, 1+δ s (s) and the redshift-space PS becomes [66,97] (2π with the velocity moments generating function where ∆u = u(x 2 ) − u(x 1 ) and x = x 2 − x 1 .Function M (or its Fourier transform) plays a central role in RSD.Different expansion procedures of eq.(3.54) yield different approaches to RSD modeling [66].Our approach follows the moment expansion (ME) approach of [97], that uses a Taylor expansion of the generating function.That is, the m-th density weighted velocity field moment of the generating function is an m-rank tensor defined as [66,97] with δ 1 = δ(x 1 ) and δ 2 = δ(x 2 ).Hence, from eq. (3.53), the power spectrum in the moment expansion approach becomes where the Ξm i 1 •••im (k) are the Fourier moments of the generating function -the Fourier transforms of their configuration space counterparts, Ξ m i 1 •••in (x).Finally, one can write [37,85] for some functions I m n (k).Hence, the momentum expansion redshift space power spectrum can be written as up to a Dirac delta function localized at k = 0.
When accounting for EFT corrections and shot-noise, the expression for the one-loop power spectrum is given by which is composed by the following elements: 1.The momentum expansion perturbation theory power spectrum where the one-loop real space power spectra P δδ , P δθ , and P θθ are presented below in §3.7.The function A TNS is defined in [98] as with the bispectrum B σ given through Meanwhile, the function D(k, µ) is given by with µ k−p the cosine of the angle between the wave-vector k − p and the line-of-sight direction n, and function F is given by [98] is the linear Kaiser power spectrum [99], but with the additional k-dependence in the growth rate.
We notice that to linear order, we can use the relation between velocity and densities given in eq.(3.44) to write and further, since functions A TNS and D are pure non-linear, one recovers the Kaiser power spectrum.In the following, we refer to the density-density linear power spectrum P L δδ simply as P L . 2. The EFT counterterms with n the average number density of galaxies, such that for a Poissonian distribution Despite the success of SPT-EFT in modeling the broadband power spectrum, the theory yet gives poor results in modeling the BAO since long-wavelength displacement fields, though being essentially linear, stream largely contributing to damp features in the power spectrum in a manner that is non-perturbative under an SPT scheme.Then, in order to model the spread and degradation of the BAO oscillations due to large scale bulk flows, we employ IR-resummations [25] as implemented by Ivanov et al. [28,29].Here, we split the linear power spectrum in a piece with the wiggles removed, P nw , and a wiggles piece, P w , such that the spectrum can be written as P L = P nw + P w .Defining P s (k, µ) as the non-resummed full-spectrum computed through eq.(3.59) using the complete linear power spectrum P L , and analogously P s,nw (k, µ), but using only the non-wiggle piece P nw , the IR-resummed power spectrum is [28] and ) The BAO peak scale is roughly given by ℓ BAO ≃ 105 h −1 Mpc, with j n the spherical Bessel function of degree n.The choice of the scale k s that splits between long and short modes is arbitrary, but the results are very robust for k s ≳ 0.1 h Mpc −1 .Our code fkpt uses the value k s = 0.4 h Mpc −1 .
Finally, to fit the data, we use the multipoles from the equation where L ℓ is the Legendre polynomial of degree ℓ.

Biasing
It is well known that there is no complete biasing theory for general theories with linear scaledependent growth (e.g.[24]), and one must add higher order derivative bias operators of the form ∇ 2 δ, . . . .However, these terms become degenerate with EFT counterterms.Hence, to describe the galaxy-matter connection, it suffices to use the EFT theory of bias of [16,17], with some tweaks studied in [85].We have the biased spectra where the function σ 3 (k) is given by eq.(3.26) of [85] and the biased spectra of the form P XY are given by eqs.(3.43)-(3.49) of the same reference.For the biased functions A TNS and D(k, µ; f ) we use

fk-perturbation theory
In the following, we describe the fk-Perturbation Theory (fkPT), which approximates the full theory described in the previous section.We further present its implementation in a fast C code, fkpt, that allows us to sample a large space of parameters with standard MCMC recipes.

fk-kernels
The power spectrum can be written as a sum of k-functions multiplied by powers of the linear growth and the cosine angle µ and the functions I(k, p) are invariant under spatial rotations.That is they only depend on the magnitudes p = |p| and k = |k| and on the angle between p and q, and therefore the functions I are indeed only functions of the wave-vector magnitude k.There are two kinds of these functions When using EdS evolution, or more precisely when Ω 2 m = f , one has explicit analytical expressions for the kernels K.However, for theories that introduce new scales one has to solve differential equations for obtaining its precise form.These equations depend on the cosmological parameters, but also on the wave-vectors.For example, one of the I functions is with (see eq. (3.48)) This kernel will serve us to illustrate our approach: We notice, there are two new types of contributions in the G 22 kernel that are not present in the EdS case: 1.The first type comes from the linear growth rates inherited from linear theory.In particular, the dipole term is determined by the advection of density fields.As discussed in [37], the second-order velocity field has a term given by and in virtue of eq.(3.44), 2. The second type of contribution arises from the presence of the functions A(k 1 , k 2 , t), B(k 1 , k 2 , t), and their derivatives.These functions appear when the equation f 2 = Ω m fails to hold true, and, in the case of MG, they carry the screening effects.
It turns out that the dominant corrections to the EdS kernels are of the first type.In the case of f (R) theories, the corrections to G 22 induced by the functions A and B are about 1-2% [69].Indeed, the lack of precise numerical values for these terms is not much more harmful than the use of EdS kernels when fitting the ΛCDM, which is a method that has proven to be quite accurate when fitting simulations and real data.On the other hand, the corrections provided by the growth rates can exhibit significantly larger magnitudes.To illustrate this, in figure 1 we plot the function f (k)/f 0 for the scenarios F4, F5, and F6 at redshifts z = 0.5 and z = 1.5.This shows that the corrections to G 22 can surpass the 15%.
The above arguments point out that the use of EdS kernels to model the LSS in scaledependent MG may not be correct.How bad is this method, clearly depends on the particular MG model.On the other hand, the use of the precise exact kernels is not computationally viable for an MCMC analysis.Hence we adopt the use of fk-kernels, which were introduced for massive neutrino cosmologies in [37] and further developed in [38].The fk-kernels, which consider the exact f (k, t) functions, can be defined as and similar for n > 2. That is, in this approach one fixes functions A and B to their large scale counterparts A LS (t) and B LS (t), obtained by evaluating all momenta at zero value.By eliminating the k-dependence from the function A(k, µ) in eqs.(3.33) and (3.34) we observe that the large scale values of A and B are equal when the screening vanishes.This behavior is expected in theories that converge to GR for large scales, such as scale-dependent theories like f (R), but not in theories with a massless scalar field mediator, such as DGP.
One can set the value of A LS and B LS to unity as in EdS, which could be a good idea for the cases in which the additional gravitational scalar degree of freedom is massive and hence the associated fifth force has a finite range, as in f (R).However, preserving their exact largescale values does not significantly impact computational time, as the differential equations only need to be solved once in the limit where all momenta tend to zero.This approach may be particularly beneficial, if not necessary, in theories that do not converge to the standard ΛCDM model at large scales, such as DGP and cubic Galileons.For the specific case of the scale-independent normal branch of DGP, ref. [45] uses the correct large scale values of the kernels and find good agreement when confronting to simulations.
Throughout this work we will use the exact large scales values of the functions A LS , B LS , ..., instead of the ones in EdS where these functions are unity.The latter scheme assumes the approximation Ω m = f 2 to be valid at large scales, and we call it EdS-fk.Figure 2 shows the power spectrum ratios for EdS-fk and fk kernels for multipoles ℓ = 0 and 2 (black and red dotted lines) for the model F5 at redshift z = 0.38.It shows that the difference is smaller than around 1% in the range of interest for full-shape analyses.

Screenings and EFT parameters
In cosmological large scale structure modeling, the screenings have played an important role and are necessary to properly match the simulated data [69,70].However, these effects can be quite degenerate with EFT counterterms.Using the F5 model at redshift z = 0.38, in fig. 2 we show with solid lines the ratio of the power spectrum multipoles ℓ = 0 and 2 without screenings and considering full kernels.Here we have chosen typical values for the cosmological, EFT, bias and noise parameters.We notice that up to k = 0.2 hMpc −1 the effect of the screenings is about 5% in the monopole and almost a 10% in the quadrupole.However, if we add a correction to the non-screened monopole of the form P 0 (k) → P 0 (k) + c 0 k 2 P L (k) and one to the non-screened quadrupole P 2 (k) → P 2 (k) + c 2 k 2 P L (k), the effect of the screenings is effectively counteracted.This is shown with dot-dashed lines, where we can see that the difference with the full kernels up to k = 0.2 hMpc −1 is smaller than 0.5%, and for the quadrupole than the 2%.The oscillations in the dot-dashed lines appear because the added linear power spectrum is not IR-resummed, and they should disappear if this is properly done.This analysis suggests there is a high degeneracy between EFT counterterms and the chameleon screening in f (R) that allows us to do not consider the latter in the analysis, which are very slow to compute.While we were finishing this work, an approximate method to treat the screenings accurately was proposed in Ref. [46].However, with the use of simulations we are capable of detecting the F5 signal with our simple prescription of ignoring the screenings.Figure 2. We show the screenings-EFT counterterms degeneracy by taking the ratio of the power spectrum multipoles ℓ = 0 and 2 without screenings and considering full kernels (solid lines).In dot-dashed lines we show the results of adding counterterms of the form αk 2 P L (k) to P 0 and to P 2 produce similar effects than the screenings.We also show the ratio of the power spectrum multipoles when using EdS-fk and fk-kernels (dotted lines), as explained in §4.1.We use the model F5 at redshift z = 0.38.
Typically, as larger are the MG strengths, the screening effects appear at lower k values, which may be a signal of the breakdown of the fkPT approximation.For f (R), this scale is given by k M 1 = a √ M 1 [83], where M 1 is given by eq.(2.17) and related to the effective mass of the new gravitational scalar degree of freedom by m = M 1 /3.Although it may be expected that the screenings operate mainly around this scale, limiting the extension of our formalism beyond it, we will not adhere to this rule-of-thumb here.Instead, we will validate fkPT against simulations for the HS-f (R) model.
Related to this topic, it's the concern of how much one can trust the fkPT method for a specific theoretical model.Ideally, one should validate the theory against simulations, if available.A more economical option is to perform a comparison against the full perturbation theory results, as we have done in this section.However, perhaps the most efficient way to work is by using linear parameterizations, and expect that EFT counterterms absorb the effect of the screening non-linearities, as a practical renormalization of the EFT coefficients.This is plausible for power spectra of the form P (k) = t(k)P MG, non-screenings (k) + (1 − t(k))P GR (k), where t(k) is a transition function admitting an expansion t(k) = 1 + ak 2 + • • • , such that the non-screened power spectrum, correctly modeled by fkPT, is obtained for large scales, GR is recovered for small scales, and intermediate scales having terms degenerate with EFT contributions.

fkpt code
fkpt is a C language code, public available at https://github.com/alejandroaviles/fkpt, that computes the one-loop redshift space tracers power spectrum using th fk-Perturbation Theory.
It receives as input the ΛCDM linear power spectrum and the cosmological parameters at the desired output redshift z.It solves eq.(3.30) to obtain linear growth function D + (k, t) and growth factor f (k, z).Then, it computes the MG power spectrum using This is an excellent approximation for several MG models, as f (R).However, the code can also receive directly the MG power spectrum obtained from another code.The code then splits the power spectrum in wiggle and non-wiggle pieces using the fast sine transform technique described in [100], then it computes the IR-resummed power spectrum given by eq.(3.69).fkpt does not use an FFTLog method since we want flexibility that allows future addition of theories that do not reduce to ΛCDM at very large scales, and hence their kernels cannot be approximated as EdS when they are evaluated at small wave vectors.That is, our code treats the large scales exact, and also serves for computing the GR power spectrum using the exact ΛCDM kernels.Despite we use a brute force approach, our code takes about 0.5 seconds in a standard personal computer to compute a single power spectrum, hence being capable of explore the parameter space with MCMC in reasonable time.
Other desired capabilities, as the Alcock-Paczyński effect and analytical marginalization, can be computed from the outside using a python interface that we provide together with the code in the github repository.

Model validation
Following the previous sections where we introduced the fk-Perturbation Theory framework and its implementation code fkpt, we are now ready to assess their performance by fitting simulated halo power spectra obtained from state-of-the-art N -body simulations.One of the main objectives of this study is to determine the capability of our method to successfully recover the MG signal in specific scenarios.We anticipate that our method is able to doing so when the signal is sufficiently strong, being effective only for F5 power spectra at low redshifts, or when we perform a joint analysis with tracers at different redshifts.Unfortunately, we were unable to detect the signal from F6, which is very close to GR.On the other hand, our study also serves to gain insight on how much we can test GR with the full shape power spectrum, so we will devote some time to the analysis of ΛCDM simulations.
Before discussing the results, we lay down a brief overview of the adopted N -body simulations.For these we use a suite of simulations performed with the code MG-GLAM [87,88], which is an extension of the parallel Particle-Mesh GLAM (GaLAxy Mocks) pipeline for fast generation of synthetic galaxy catalogs [101].The set of simulations contains runs with the following ΛCDM reference cosmology {Ω b , Ω m , h, n s , ln(10 10 A s )} = {0.0486,0.3089, 0.6774, 0.9667, 3.01887} corresponding to the best fit values of Planck 2015 (last column in table 4 of [102]).Apart from the GR results, two instances of the HS-f (R) model were simulated (among other MG models not considered in this work), corresponding to |f R0 | = 10 −6 and 10 −5 , that we refer as to F6 and F5, respectively.The simulations are performed over a cubic box of a comoving volume of (1024 h −1 Mpc) 3 , with grid size of 4096 3 , and including 2048 3 dark matter particles in each realization.All realizations were initialized with Zeldovich Approximation initial conditions at redshift z = 100 using the z = 0 power spectrum of the ΛCDM reference cosmology extrapolated to the initial redshift.This implies that the σ 8 value is different for the different models, but the primordial amplitude remains the same.These values are σ GR 8 = 0.8161, σ F6 8 = 0.8292 and σ F5 8 = 0.8694.Gravitationally bound halos were identified using the bound density maxima (BDM) spherical overdensity halo finder [103], selecting halos within the mass range 10 12.5 < M h < 10 13 h −1 M ⊙ .
For our discussions, we have selected five specific redshifts widely used in the literature.These include z = 0 and z = 0.48, which is utilized as a proxy of z = 0.5.Furthermore, we considered redshifts z = 0.38 and 0.61 as they coincide with the distinct, non-overlapping bins known as z 1 and z 3 in the BOSS DR12 dataset [10].Finally, z = 1.5 corresponds to the redshift of QSOs in the BOSS DR16.These particular redshift values have been extensively employed in joint analyses of 2-point statistics.
In fig. 3 we show plots for the halo power spectrum at different redshifts with the error bars arising from the covariance matrix rescaled by a factor of 1/5 as explained below.We notice there is not a clear pattern in the halo power spectra of different gravity theories.For the matter power spectrum, the MG models present more power than GR at all scales because the strength of gravity is larger in MG.However, for the halo spectrum the situation is very different because of the halo large scale bias.It is known from previous works that bias evolution differs among MG models [84,[104][105][106][107], with a tendency for smaller linear bias b 1 as the MG strength increases, see for example fig. 4 of ref. [108].This phenomenon further complicates the differentiation between MG models.Figure 3 also emphasizes the significance of bias evolution in MG.Notably, for redshifts z = 0.38, 0.5 and 0.61, models F5 and F6 are nearly indistinguishable, as they overlap within the error bars of our simulations, even after rescaling the covariance matrix.Conversely, at redshift z = 1.5, models GR and F6 overlap, and model F5 apparently could be distinguished, however the observed differences are almost entirely due to a different large scale bias, since the clustering effects of MG are very small at such a large redshift.This can be confirmed by examining the ratio between the same multipole for two different models, e.g.GR to F5, which remains nearly constant with the wavenumber k.We conclude, that for our chosen halos the only redshift at which the three models could be clearly distinguished is z = 0.This shows the importance of employing joint analysis of different redshifts or tracers when testing gravity.
For each of the examined data sets, we fit the mean of the 100 realizations.The power spectrum multipoles covariance of a single realization is determined by where P (k i ) is the value of the power spectrum of a single realization at bin k i and P (k i ) = ⟨ P (k i )⟩ is the mean over the ensemble of realizations at bin k i .However, since the volume of each simulation is (1024 Table 1.Cosmological and nuisance parameters and Gaussian (G) and uniform (U) priors used for fitting the MG-GLAM simulated data.We fix the primordial spectral index to n s = 0.9667.Tidal bias b s 2 and b 3nl are obtained from coevolution using eqs.(5.4).
data as cov25 and cov100.With these we construct a Gaussian Likelihood L ∝ exp(−χ 2 /2), with χ 2 = (P theory − P sim ) T C −1 N (P theory − P sim ) ( with P sim = P sim ℓ (k i ) the data vector, the P theory = P theory ℓ (k i ) the model vector, and C N is the covariance matrix in eq.(5.1), rescaled either by the factor 1/N , with N = 25 or 100.
Our fitting set consists of 11 parameters: four cosmological {h, ω c , ω b , ln(10 10 A s )}, one MG parameter f R0 , two local biases {b 1 , b 2 }, two EFT counterterms {α 0 , α 2 } and the two shot noise {α shot 0 , α shot 2 }.The spectral index n s is fixed to the simulation value.Meanwhile, tidal bias b s 2 and third-order non-local bias b 3nl are fixed by co-evolution theory [109][110][111] to These expressions are valid within GR when assuming local Lagrangian bias.However, previous works have found that these same relations yield good results for MG simulations [85], as well as in massive neutrino cosmologies [38].
In our fittings we adopt uniform flat priors in all parameters with the exception of ω b , for which we use a Gaussian prior centered at the value of the simulations but a width given by Big Bang Nucleosynthesis (BBN) observations [112,113].None of the posterior distributions saturate the flat priors, indicating that they can be considered as uninformative.A list of Figure 5. F5 detection: We present the posterior distribution for fitting the F5 MG-GLAM simulations at redshift z = 0 in the cov25 and cov100 cases.The parameter f R0 is varied with a flat prior over the interval (−0.01, 0.01).The confidence intervals for the absolute values are given in table 2.
In our pipeline, the fkpt code treats only the absolute value of f R0 , ignoring its sign.This means that a signal detection would cause the MCMC chains to randomly sample the posterior around either the values f R0 = −10 −5 or 10 −5 (shown as dashed vertical gray lines).For the sake of transparency, we do not smooth the MCMC in this plot.Our analysis provides unequivocal evidence that our pipeline successfully detects the MG F5 signal.
all priors is provided in table 1.We have chosen a flat prior on |f R0 |, instead of the perhaps more natural option of a flat prior over log |f R0 |.We do this because our preliminary results have shown a better performance, getting a closer |f R0 | best fit value when fitting the F5 simulations for z = 0, as shown in fig. 4. In this plot we utilize a uniform prior U(−10, −1) When reporting the parameters in figures and the main tables, we do it only for the cosmological parameters and the linear bias b 1 .We omit ω b for which the posteriors are entirely dominated by the prior, which is very tight.Furthermore, for the matter density we combine the baryons and cold dark matter abundances and use Ω m instead of ω m = ω b + ω c to avoid showing the trivial degeneracy with h.
The linear ΛCDM power spectrum at redshift z 0 = 0 is obtained from the Einstein-Boltzmann code CLASS9 [114] which serves as an input of fkpt, which first obtains the MG power spectrum using eq.(4.9).Subsequently, the code computes the loop corrections in eqs.(3.69) and finally the halo power spectrum multipoles ℓ = 0 and 2 through eq.(3.73) that we compare against the simulated data.To sample the parameter space, we perform MCMC runs with the code emcee 10 [115] which is based in the affine-invariant ensemble sampler method [116].The contour and 1-dimensional posterior plots as well as the confidence intervals are computed using the GetDist Python package [117].For presentation purposes, in all figures with the exception of fig. 5 we use a GetDist smoothing scale of 0.7.However, the confidence intervals we present in all tables are computed without any previous smoothing.
In fig. 5 we show the posterior distribution when fitting the F5, MG-GLAM simulations at redshift z = 0, considering both the cov25 case (blue lines) and cov100 case (red lines), and utilizing a maximum value k max = 0.17 h Mpc −1 in the power spectrum.As detailed in Table 1, we vary the MG parameter f R0 with a flat prior ranging from −0.01 to 0.01.However, our fkpt code only considers the absolute value of f R0 , being insensitive to its sign.We opt for this symmetric prior instead of a simpler range (0, 0.01) to avoid encountering edge effects, which can arise when the MCMCs approach the boundaries of the interval.These boundary regions are precisely where we expect to find the signal, hence we aim to prevent any bias from such effects.We anticipate that a signal detection would cause the chains to randomly sample the posterior around either f R0 = −10 −5 or 10 −5 .This behavior is clearly observed in fig.5, particularly in the cov100 fitting case.Additionally, we omit to smooth the MCMC in this plot to avoid any visual ambiguity.Instead, we opt to exhibit a histogram.Overall, this analysis provides strong evidence that our theoretical model and numerical implementation have successfully detected the MG F5 signal.
In table 2 we display the one-dimensional posterior confidence intervals when performing the analyses at redshift z = 0 for the case of F5, as described above, as well as for F6 and GR.We notice that we recover all the cosmological parameters within the 1-or 2-σ intervals, with the exception of A s , for which we underestimate the true value.by fig.6, that shows triangular plots including the 2-dimensional contours and 1-dimensional distributions for the for cov25 case (left panel) and cov100 (right panel).We notice we obtain a 1-σ detection of the MG-F5 signal for the case cov25, and almost 2-σ for cov100.However, no detection is found in either case for the weaker gravity model F6.As explained above, in fig. 3 we displayed plots for the data fitted in this analysis, which suggest that the F6 and GR models yield very similar results on these halo mass range, while F5 can be clearly distinguished.
To compare our numerical fittings against a theoretical model, we calculated the linear large-scale halo bias using the Peak-Background Split (PBS) theory as described in [84].In scale-dependent MG, the critical threshold for collapse δ c -the minimum density fluctuation necessary for a space region to collapse and create a halo-varies as a function of the enclosed mass within that region.This differs from General Relativity (GR), where this threshold remains constant.We compute δ c (M ) using the formulae of [118].Further, we use a Sheth-Tormen halo mass function, which has proven to be universal also for HS-f (R) [84] when seen as a function of the peak significance ν = δ c /σ R , instead of the variance σ R , as is done in some works, e.g.[119].For a redshift of z = 0, we obtained the following values for the linear halo bias b F5, PBS These values are in good agreement with our results as can be observed from table 2.
To continue with the analysis, we performed fits to the F5 simulations at redshift 0.38   for covariance rescaled by a factor of 1/25 (cov25).Unlike the redshift z = 0, in this case we did not recover the MG signal, as can be seen in the left panel of fig.7 and in table 3.
For this reason, we performed a joint analysis, named z 1,3 , that includes halos at redshifts 0.38 and 0.61.In this case, we were able to obtain the value |f R0 | = 10 −5 within the 0.68 confidence intervals.We have seen in Figure 2 and discussed above that it is difficult, if not impossible, to discern between F6 and F5 models for these two redshifts.Therefore, we also ran chains for F6 and show our results in the right panel of fig.7 and in table 3. We notice that this analysis also detects a MG signal, but it is also located close to |f R0 | = 10 −5 , as we found for F5.However, the estimated linear bias is indeed different for cases F5 and F6, and these are in agreement with our theoretical results obtained using PBS, which are b F5, PBS Our approach can effectively estimate parameters in a standard ΛCDM cosmology, similar to other codes such as Class-PT, PyBird, FOLPSν, and Velocileptors.However, our method offers a potential advantage for future surveys by employing the exact ΛCDM kernels within our fkpt code, in contrast to the commonly used EdS kernels when fitting data.To demonstrate this, we conducted an analysis fitting General Relativity (GR) simulations while deactivating the effects of MG.This was achieved by setting f R0 to a very small value, such as 10 −10 or any other very small number in fkpt.We performed this analysis for two scenarios: the redshift z 0 = 0 and a joint analysis of two redshifts, z  7.  4. GR: One-dimensional constraints on the parameter f R0 in fits using GR simulations.The results are presented for redshift z 0 = 0 and for the joint analysis z 1,3 , the later corresponding to redshifts z = 0.38 and z = 0.61.The covariance is rescaled as cov25.The error bars indicate the 0.68 confidence intervals.The reported b 1 for the z 1,3 case represents the linear bias of the halos at z = 0.38.
the cosmological parameters within the 2-σ confidence intervals, with the exception of the underestimated A s .To explore the effects of activating MG, we repeated the fittings with a free parameter f R0 .The results, shown also in fig.8 and Table 4, are somewhat surprising.We observe no significant differences in the posterior distributions, except for a slight offset in the best fit of Ω m and broadening in its distribution.This indicates a lack of substantial degeneracy between the cosmological parameters and f R0 at large scales, with the possible exception of Ω m .
Finally, we want to compare the dependence of our fittings on the maximal wave-number.The results are presented in fig.9, where we show only the means of the full-shape analyses with k max = 0.15, 0.16, 0.17, 0.18, 0.19, 0.20 and 0.21 h Mpc −1 an their 1-σ error bars.For these analyses we have used the F5 model simulations at z = 0 with the maximum available volume given by the cov100 rescaling.This plot further justify our baseline choice of k max = 0.17 h Mpc −1 .

Comparing to NSeries simulations
Finally, since we are worried about the lack of precision in our fittings to MG-GLAM simulations, particularly in A s and to a lower extent into h, we opt to use a larger set of simulations, although they exist only for GR.Specifically, in this subsection we utilize the cubic boxes of the NSeries galaxy mocks, comprising 7 realizations, each one with a volume of V 1 = (2.7 h −1 Gpc) 3 [10]. 12These simulations were initially generated to investigate systematic effects within the BOSS data pipeline.The expanded volume offered by this dataset grants us greater confidence in testing our perturbative theoretical model, which extracts the cosmological information from the larger scales where simulations with smaller sizes, as MG-GLAM, are not optimal.The cosmological parameters are Ω M = 0.286, n s = 0.97, h = 0.7, ln(10 10 A s ) = 3.06619 and Ω b = 0.047.The covariance matrix is constructed using the NGC 1,000 EZmocks catalogues [120] available at the same URL.
We fit to the mean galaxy power spectra of the 7 realizations, with the covariance rescaled to the maximum allowed volume (N = 7), such that the effective volume of the simulations is V 7 = 137.8(h −1 Gpc) 3 .This analysis besides permitting us to have better confidence at the large scale fittings, allow us to test our modeling in a limiting case of unrealistic small errors.In fig. 10 we show triangular plots for our fittings when letting free |f R0 | and when keep it fixed to zero, these are accompanied by table 5. We notice that our fits are quite precise, matching the parameter of the simulations with high accuracy and precision even when the effects of MG are allowed.Table 5. NSeries: One-dimensional constraints on the parameter f R0 in fits using NSeries GR simulations.The results are presented for redshift z = 0.5.We fit the mean power spectrum of the 7 cubic boxes each with volume V 1 = (2.7 h −1 Gpc) 3 .

Conclusions
In recent years, full-shape methods, also known as full-modeling or direct-fit, have become a standard approach for extracting cosmological information from galaxy surveys.This shift from the fixed-template classical-analysis resulted from advancements in the EFT of Structure Formation.This framework extends PT by incorporating biases, counterterms and shot noise parameters, as well as the use of IR-resummations to model the smearing of the BAO.Much of this progress has assumed a ΛCDM model where linear growth is scale-independent.While this approach adequately describes scenarios involving additional scales which produce small effects in clustering, such as those with massive standard model neutrinos, it may fall short for more generalized MG theories.In these theories, the two scalar gravitational potentials (in Newtonian Gauge) differ even in the absence of anisotropic stresses in the matter content, introducing scale-dependent terms into the Poisson equation.The pioneering work by [70] introduced perturbative methods to accommodate these new scales within the kernels, further developed in various scenarios in the literature.However, a major challenge of these methods lies in the absence of algebraic expressions for the perturbative kernels.Instead, they must be derived by solving differential equations for each wave-vector configuration and each set of cosmological parameters.This significantly slows down the computation of statistic one-loop corrections, making parameter exploration, such as through MCMC, exceedingly time-consuming.The present work addresses this issue building-upon the formalism introduced in [37,38] for massive neutrinos.We identify two types of contributions within the PT kernels that are absent in the EdS model.Firstly, the introduction of scale and time-dependent functions A and B, and their third order counterparts, due to the deviation of the growth factor f from being equal to Ω 1/2 (a) and because of the screening contributions that should drive the theory to GR at high k.While these functions exist in the ΛCDM model, they only depend on time in that context.Secondly, contributions arise from the scale-dependence of the growth factors f = f (k, t), resulting in differences between velocity and density fields at linear order, encapsulated in θ (1) (k, t) = (f (k, t)/f (k = 0, t))δ (1) (k, t).Due to the advection of large scale density fields, this property is inherited to higher orders in the perturbative kernels.
In scenarios such as HS-f (R) models, the first type of contribution is predominantly influenced by non-linear screening effects, often degenerate with EFT counterterms.Thus, we approximate these functions with their largest scale values, computed where all wave-vector arguments go to zero.Notably, for scale-dependent theories converging to the ΛCDM model at large scales, the approximated values correspond to the genuine ΛCDM values, unlike the EdS approximation where these functions equal unity.Hence, this method also serves us to test GR models using exact kernels.
Our approach, named f (k)-Perturbation Theory (fkPT), keeps the factors f (k)/f (k = 0) that represent the dominant contribution to the kernels, while maintaining the rest of the features as in their large scale limit.We have developed and released the fkpt code, enabling the computation of the redshift space power spectrum for scale-dependent MG.We validate our method and code using the MG-GLAM simulations for HS-f (R) models F6 and F5, alongside GR.While we successfully recovered the MG signal for F5 at redshift z = 0, indicating its significant impact on the power spectrum, we were not able of doing so at z = 0.5, since the MG signal is weaker in this case.However, we obtain the signal at higher redshifts through a joint analysis of z 1 = 0.38 and z 3 = 0.61.Additional examinations using DESI-like simulations and estimation of parameters using real data will be addressed in a future work.
Additionally, the evolution of large-scale bias differs across various gravity theories, suggesting different values of b 1 for the halos used in our analysis.We utilized the Peak-Background-Split formalism with a Sheth-Tormen mass function to theoretically derive large scale bias by using non-constant threshold density δ c (M ).Our analytical results coincide reasonably well with values obtained through a MCMC analysis, further demonstrating the efficacy of our method in recovering the MG signal from simulated data.
However, we acknowledge a limitation in our ability to recover the amplitude of the primordial perturbations.This shortcoming may be attributed to the modest size (L = 1024 h −1 Gpc) of the MG-GLAM simulations, insufficient for thoroughly testing the large scales that are crucial for extracting cosmological information from full-shape analysis.Further, to find the MG signal we were forced to rescale the covariance by factors of 1/25 and 1/100, reaching effective volumes up to V = 100 (h −1 Gpc) 3 , for which the simulations may not be sufficiently accurate.Therefore, we additionally use the cubic boxes of the NSeries simulations, primarily used to validate the BOSS pipelines.Despite being available only for GR, we successfully utilized these simulations, allowing the MG parameter |f R0 | to vary freely and achieving excellent results using our methodology.
The advent of the Dark Energy Spectroscopic Instrument (DESI) and Euclid spectroscopic surveys opens new opportunities to test models beyond the standard ΛCDM with the clustering of galaxies.In particular, MG models that influence the late times clustering can elude probes based on CMB observations, or even kinematics tests in cosmological distances imposed by, e.g, Supernovae type Ia.Hence, new non-linear methodologies that account for additional scales, but at the same time being sufficiently fast to be implemented in MCMC samplers for parameter estimation, would prove valuable for investigating gravitational effects in the forthcoming years.Future endeavors in this line of research will steer towards this objective, together with the joint utilization of diverse cosmological probes.

. 4 )
considering the fluid perturbation ∆ρ = ρδ and the MG associated scalar field perturbation δf R = f R − fR , R = R + δR, where the bar indicates background quantities, and R ≡ R( fR ).The perturbative field equations in Fourier space are[70]

Figure 1 .
Figure 1.Linear growth function f as a function of k for different redshifts (z = 0 and z = 1.5) and different gravitational strengths (F4, F5 and F6).All cases are normalized to one at large scales by showing f (k, z)/f 0 (z).

Figure 3 .
Figure 3. MG-GLAM simulated halo power spectrum multipoles ℓ = 0 and 2. The central points are the mean of the 100 simulations and the error bars are their RMS which is further divided by 1/5 to show the errors for the rescaled covariance cov25.

Figure 4 .
Figure 4. Effect of priors on f R0 : We show fittings to F5 model power spectrum when opting for a flat prior on f R0 over the interval (-0.01,0.01)(blue solid line) and a flat prior on log 10 (|f R0 |) over the interval (-10,-1) (green dashed line).
we rescale the covariance C by factors of 1/25 and 1/100, considering effective volumes of ∼ 25 h −3 Gpc 3 and ∼ 100 h −3 Gpc 3 .We refer to these sets of

Figure 6 .
Figure 6.Triangle plots for fitting GR, F6 and F5 simulations at redshift z = 0, for the cov25 case (left panel) and cov100 (right panel).The shadows show the 0.68 and 0.95 confidence intervals.

Figure 7 .
Figure 7. Triangle contour plots from fits using F5 (left panel) and F6 (right panel) power spectra from simulated halo catalogs.The results are presented for redshift z = 0.38 and for the joint analysis of redshifts z = 0.38 and z = 0.61.The covariance is rescaled as cov25.The shadows indicate the 0.68 and 0.95 confidence intervals.This figure is accompained by table 3.

b 1 Figure 9 .
Figure 9. Impact of k max in the fkPT full-shape analysis, when fitting the MG-GLAM simulations, with F5 and redshift z = 0.The dots indicate the means and the bars are the 1-σ errors.The shadowed region corresponds to our baseline choice of k max = 0.17 h Mpc −1 throughout this work.