Convolution Lagrangian perturbation theory for biased tracers beyond general relativity

TL;DR

This paper evaluates and extends Lagrangian Perturbation Theory methods to accurately predict the clustering statistics of biased tracers in modified gravity models, incorporating environmental effects and screening mechanisms.

Contribution

It introduces a convolution Lagrangian perturbation theory approach that effectively models biased tracers in scalar-tensor modified gravity theories beyond general relativity.

Findings

CLPT and SPT provide accurate predictions for correlation functions and power spectra.

The methods characterize BAO, power-law, and small-scale regimes effectively.

Analytic bias parameters can incorporate environmental sensitivity in modified gravity.

Abstract

We compare analytic predictions for real and Fourier space two-point statistics for biased tracers from a variety of Lagrangian Perturbation Theory approaches against those from state of the art N-body simulations in $f(R)$ Hu-Sawicki and the nDGP braneworld modified gravity theories. We show that the novel physics of gravitational collapse in scalar tensor theories with the chameleon or the Vainshtein screening mechanism can be effectively factored in with bias parameters analytically predicted using the Peak-Background Split formalism when updated to include the environmental sensitivity of modified gravity theories as well as changes to the halo mass function. We demonstrate that Convolution Lagrangian Perturbation Theory (CLPT) and Standard Perturbation Theory (SPT) approaches provide accurate analytic methods to predict the correlation function and power spectra, respectively, for biased tracers in modified gravity models and are able to characterize both the BAO, power-law and small scale regimes needed for upcoming galaxy surveys such as DESI, Euclid, LSST and WFIRST.

Figures

Tables

Table 1. Table 1: The table presents the values for the best-fit Sheth-Tormen parameters ( q , p ) 𝑞 𝑝 (q,p) for the halo mass function ( 47 ), with respect to the simulations, through the criterion ( 82 ), as well as the bias factors b 1 subscript 𝑏 1 b_{1} and b 2 subscript 𝑏 2 b_{2} evaluated through ( LABEL:pbsbiasMGavg ) using the best-fit values. All the evaluations for the Group I simulations were performed at redshift z = 0.5 𝑧 0.5 z=0.5 and for the Group II simulations at z = 1 𝑧 1 z=1 . The labels low, mid. and high indicate reference to the three mass bins, specified in the text. For all f ​ ( R ) 𝑓 𝑅 f(R) models, the bias values shown are environmentally averaged as described in the text.

	Best-fit ST		Predicted Biases
Models	$q$	$p$	$b_1$	$b_2$
$GroupI:GR$	0.726	0.345	0.301	-0.501
$GroupI:F4$	0.671	0.361	0.120	-0.435
$GroupI:F5$	0.765	0.321	0.211	-0.470
$GroupI:F6$	0.670	0.362	0.230	-0.449
$GroupI:N1$	0.701	0.369	0.224	-0.661
$GroupI:N5$	0.702	0.357	0.268	-0.503
$GroupII:GRLow$	0.674	0.362	0.345	-0.183
$GroupII:GRMid.$	0.728	0.342	0.925	-0.05
$GroupII:GRHigh$	0.806	0.594	1.720	1.900
$GroupII:F5Low$	0.733	0.314	0.295	-0.170
$GroupII:F5Mid.$	0.788	0.282	0.909	-0.033
$GroupII:F5High$	0.746	0.305	1.491	0.416

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{R+f(R)}{16\pi G}+{\mathcal{L}}_m\right],$$

$$f(R)=-m^2\frac{c_1{\left(R/m^2\right)}^n}{c_2{\left(R/m^2\right)}^n+1}.$$

$$\bar{R}=3{\Omega }_m{H}_0^2\left(1+4\frac{{\Omega }_{\Lambda 0}}{{\Omega }_{m0}}\right).$$

$${\bar{f}}_{R_0}=-n\frac{c_1}{c_2^2}{\left(\frac{{\Omega }_{m0}}{3({\Omega }_{m0}+{\Omega }_{\Lambda 0})}\right)}^{n+1}.$$

$${\nabla }^2{\Phi }_N$$

$${\nabla }^2{f}_R$$

$$S=\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi G}+{\mathcal{L}}_m\right]+\int d^{5}x\sqrt{-g_5}\left(\frac{R_5}{16\pi G{r}_c}\right).$$

$$ϵ\frac{H}{r_c}=H^2-\frac{{\rho }_m}{24\pi G},$$

$${\nabla }^2{\Phi }_N$$

$${\nabla }^2\phi$$

$$\beta (a)=1+2H(a)r_c\left(1+\frac{\dot{H(a)}}{3H(a{)}^2}\right).$$

$$x(q,t)=q+\Psi (q,t).$$

$$1+{\delta }_m(x,t)=\int d^{3}q{\delta }_D\left[x-q-\Psi (q,t)\right]=\frac{1}{J(q,t)},$$

$$J_{ij}=\frac{\partial x^i}{\partial q^j}={\delta }_{ij}+\frac{\partial {\Psi }^i}{\partial q^j}.$$

$$\ddot{x}+2H\dot{x}$$

$$\frac{1}{a^2}{\nabla }_x^2\psi (x,t)$$

$$\Psi (q,t)=\sum _{n=1}^{\infty }{\Psi }^{(n)}(q,t)={\Psi }^{(1)}(q,t)+{\Psi }^{(2)}(q,t)+{\Psi }^{(3)}(q,t)...$$

$$(3+2{\omega }_{BD})\frac{1}{a^2}{k}_x^2\varphi ({\mathbb{k}}_x,t)=8\pi G{\bar{\rho }}_m\delta ({\mathbb{k}}_x,t)-\mathcal{I}(\varphi ),$$

$$\mathcal{I}(\varphi )$$

$${\nabla }_x\hat{\mathcal{T}}\Psi$$

$$\frac{1}{a^2}{\nabla }_x^2\psi (x,t)$$

$${\left(J^{-1}\right)}_{ij}\hat{\mathcal{T}}{\Psi }_{i,j}(\mathbb{k})=-A(k)\delta (\mathbb{k})+\frac{k^2}{a^{2}3\Pi (k)}\delta \mathcal{I}(\mathbb{k})$$

$$+\frac{M_1(k)}{3\Pi (k)}\frac{1}{2a^2}\left({\nabla }_x^2\varphi -{\nabla }^2\varphi \right)(\mathbb{k}),$$

$$A(k)$$

$$\Pi (k)$$

$$\mathbb{k}\cdot {\Psi }^{(1)}=i{D}^{(1)}(k,t){\delta }^{(1)}(\mathbb{k},t=0),$$

$${\Psi }^j(k,t)=\frac{i{k}^j}{k^2}{D}^{(1)}(k,t){\delta }^{(1)}(\mathbb{k},t=0).$$

$$\hat{\mathcal{T}}{D}^{(1)}(k,t)=A(k)D^{(1)}(k,t).$$

$$\mathbb{k}\cdot {\Psi }^{(2)}=\frac{i}{2}\int \frac{d^3{k}_1{d}^3{k}_2}{{\left(2\pi \right)}^3}{\delta }_D(\mathbb{k}-{\mathbb{k}}_{12})D^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2){\delta }_1^{(1)}{\delta }_2^{(1)},$$

$$D^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2)$$

$$\left(\hat{\mathcal{T}}-A(k)\right)D_a^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2)$$

$$\left(\hat{\mathcal{T}}-A(k)\right)D_b^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2)$$

$$\left(\hat{\mathcal{T}}-A(k)\right)D_{\delta \mathcal{I}}^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2)$$

$$\left(\hat{\mathcal{T}}-A(k)\right)D_{FL}^{(2)}({\mathbb{k}}_1,{\mathbb{k}}_2)$$

$$M_n=\frac{d^n\bar{R}(f_{R_0})}{d{f}_{R_0}^n},$$

Full text

We show that the novel physics of gravitational collapse in scalar tensor theories with the chameleon or the Vainshtein screening mechanism can be effectively factored in with bias parameters analytically predicted using the Peak-Background Split formalism when updated to include the environmental sensitivity of modified gravity theories as well as changes to the halo mass function.

We demonstrate that Convolution Lagrangian Perturbation Theory (CLPT) and Standard Perturbation Theory (SPT) approaches provide accurate analytic methods to predict the correlation function and power spectra, respectively, for biased tracers in modified gravity models and are able to characterize both the BAO, power-law and small scale regimes needed for upcoming galaxy surveys such as DESI, Euclid, LSST and WFIRST.

I Introduction

About 20 years after the discovery of cosmic acceleration (Perlmutter et al., 1999; Riess et al., 2004), the $\Lambda$CDM model of the universe is rightfully called the standard model of cosmology, after having successfully passed a broad range of tight observational tests (Eisenstein et al., 2005; Percival et al., 2007, 2009; Kazin et al., 2014; Spergel et al., 2013; Ade et al., 2013, 2016). This model augments the cosmic baryonic and leptonic matter components in the Standard Model, with two additional ingredients that dominate the cosmic energy budget, cold dark matter (CDM) and a cosmological constant, $\Lambda$, with the latter playing the role of vacuum energy that is responsible for the observed late-time acceleration. Central to the above paradigm is also the assumption that gravitational evolution is governed by Einstein’s General Theory of Relativity (GR), which led to the observed inhomogeneous pattern of the Large-Scale Structure (LSS) of the Universe through the process of gravitational instability. The value of the vacuum energy predicted by quantum field theory, however, is orders of magnitude larger than the best-fit one that is necessary to explain cosmic acceleration, so $\Lambda$ needs to be fine-tuned; the so-called cosmological constant problem (Weinberg, 1989). Such an unfortunate mismatch, together with the need to fully explore the space of all theoretical alternatives, has generated growing interest in testing gravity models that self-accelerate though a large-scale modification to GR, instead of dark energy, the so-called Modified Gravity (MG) models (Koyama, 2016; Ishak, 2018).

Modifying the Einstein-Hilbert action, however, introduces in principle an additional degree of freedom that is conformally coupled to matter and can produce significant deviations from the predictions of GR, which have passed a wide array of precise observational tests, especially in the Solar System (Will, 2006). Furthermore, the recent simultaneous detection of gravitational waves and EM counterparts by the LIGO/Virgo collaboration (Abbott et al., 2017a; Goldstein et al., 2017; Savchenko et al., 2017; Abbott et al., 2017b, c), has placed additional constraints (Sakstein and Jain, 2017; Ezquiaga and Zumalacárregui, 2017; Creminelli and Vernizzi, 2017; Baker et al., 2017) into the form of the most general expression of a scalar-tensor theory that produces second order equations of motion, described by the Horndeski action (Horndeski, 1974; Deffayet et al., 2011).

In order to be able to confront such tight constraints successfully, while at the same time provide a stable self-accelerative cosmic mechanism, viable candidates contain a restoring property, called “screening” (Khoury, 2010, 2013), which is a dynamical mechanism that weakens the additional fifth forces in high-density environments through the corresponding scalar field self-interactions. In the Vainshtein mechanism (Vainshtein, 1972; Babichev and Deffayet, 2013), GR is recovered thanks to the second derivative terms in the scalar field Lagrangian, that become large in high density environments and effectively weaken the coupling to the matter sources. The Vainshtein mechanism is very efficient in the vicinity of a massive source and contains a rich phenomenology, which makes it particularly attractive. Another popular class of screening consists of the chameleons (Khoury and Weltman, 2004a, b), where in regions of high potential the scalar fields become massive and cannot propagate, resulting thus in suppression of the fifth forces. Despite the fact that chameleons cannot produce self-acceleration (Wang et al., 2012), their very interesting phenomenology makes them serve as ideal testbeds for gravity and considerable efforts have been put into their study in the past decade (Burrage and Sakstein, 2018). Other screening classes include the symmetrons (Hinterbichler and Khoury, 2010; Olive and Pospelov, 2008), that employ spontaneous symmetry breaking and share qualitative similarities with the chameleons and the K-Mouflage (Dvali et al., 2011; Babichev and Deffayet, 2013), in which deviations are suppressed when scalar field gradients exceed a certain value.

The observed inhomogeneous LSS of the Universe, is the outcome of the subsequent nonlinear gravitational evolution of the primordial density fluctuations, partially modulated by the late-time acceleration. As a result, it provides us with an observational window into the fundamental physics that shaped this process, including sensitive tests of the underlying large-scale gravitational law, making it particularly valuable for constraining the various MG models. Indeed, as we are entering the era of “precision cosmology”, multiple spectroscopic and photometric surveys of the LSS, both already operating, such as the DES (Abbott et al., 2005) and also about to be commissioned in the next decade, like DESI (Levi et al., 2013), Euclid (Laureijs et al., 2011), the LSST (Abell et al., 2009) and WFIRST (Spergel et al., 2013), will provide us with particularly precise maps of the LSS that will shed light on the mysterious nature of the dark sector. Taking full advantage of this wealth of cosmological information poses a great challenge for experiment and theory alike.

From a theoretical standpoint, models of structure formation rely upon accurately tracing the nonlinear evolution of dark matter perturbations. In the linear regime and when gravity is governed by GR, different modes evolve independently, with the time evolution encapsulated in the scale-independent growth factor. At nonlinear scales, however, the dynamics of the self-gravitating dark matter system can only be tracked accurately by full-blown N-body simulations, which are highly computationally expensive . Additional complexity arises in accounting for the fact that the observed galaxies do not perfectly trace the underlying dark matter density field, but are biased tracers of it (Kaiser, 1984). While this effect can be easily captured by introducing a multiplicative bias factor in the large scales (Efstathiou et al., 1988), in the regime of nonlinear dynamics, empirical modeling needs to be combined with sophisticated simulations in order to predict the spatial distribution of galaxies inside gravitationally collapsed dark matter halos (Berlind and Weinberg, 2002; Kravtsov et al., 2004; Zheng et al., 2005). Furthermore, when MG configurations are considered, the fifth forces introduce an additional layer of complexity, scale-dependent growth is generally present even at the linear level and in the nonlinear scales, one needs to solve the scalar field Klein-Gordon (KG) equation that is highly nonlinear and adds to the computational costs significantly. The intermediate, quasi-linear, scales can fortunately be analytically accessed by higher order Perturbation Theory (PT) (Bernardeau et al., 2002; Carlson et al., 2009) approaches or hybrid methods (Tassev et al., 2013; Valogiannis and Bean, 2017).

The Lagrangian Perturbation Theory approach (Zeldovich, 1970; Buchert, 1989; Bouchet et al., 1995; Hivon et al., 1995; Taylor and Hamilton, 1996; Matsubara, 2008a, b; Carlson et al., 2013; Matsubara, 2015; Desjacques et al., 2018) to structure formation is one the oldest and most popular analytical frameworks in the literature, that has been been particularly successful at describing the Baryon Acoustic Oscillation (BAO) peak (Eisenstein et al., 2005), observed at a comoving scale of $\sim$110 Mpc/$h$, and the power-law correlation function, on comoving scales $\sim$20-90 Mpc/$h$, in $\Lambda$CDM (Vlah et al., 2015a). Combined with a model for halo bias (Kaiser, 1984; Mo et al., 1997; Sheth and Tormen, 1999; McDonald and Roy, 2009; Schmidt et al., 2013; Matsubara, 2008b), it can be used to predict the 2-point statistics for halos in the real and redshift space (Carlson et al., 2013; Matsubara, 2008b; Vlah et al., 2016), which serves as a crucial step to the theoretical description of galaxy clustering. Additional contributions from small-scale physics can also be included using techniques inspired by effective field theory theory (Vlah et al., 2015b, 2016). In the context of MG cosmologies, extensive studies have been performed in the framework of Eulerian Standard Perturbation Theory (SPT) (Koyama et al., 2009; Taruya et al., 2014a; Brax and Valageas, 2013; Taruya et al., 2014b; Bellini and Zumalacárregui, 2015; Fasiello and Vlah, 2017; Bose and Koyama, 2017; Bose and Taruya, 2018). LPT was first found to work very well within the COLA hybrid framework for chameleon and Vainshtein MG cosmologies (Valogiannis and Bean, 2017), while third order LPT for dark matter was recently developed in the case of scalar-tensor theories in (Aviles and Cervantes-Cota, 2017).

In this work, we perform a comprehensive study of how LPT can be used to make predictions for biased tracers in modified gravity theories, and how well the predictions compare with full numerical simulations for a variety of modified gravity models. We study chameleon and Vainsthein MG theories, focusing on the $f(R)$ Hu-Sawicki (Hu and Sawicki, 2007) and $n$DGP braneworld models (Dvali et al., 2000) as popular, representative examples for each category. The underlying dark matter clustering is described using the formalism in (Aviles and Cervantes-Cota, 2017). We then consider the Convolution Lagrangian Perturbation Theory (CLPT) (Carlson et al., 2013) for biased tracers, as well as the variants, using the particular resummation scheme in (Vlah et al., 2015b, 2016). We extend the peak-background split formalism (PBS) (Kaiser, 1984; Bardeen et al., 1986; Schmidt et al., 2013) in which the Lagrangian bias factors are calculated as responses, to account for modifications to the halo mass function in modified gravity theories and the environmental dependence of screening effects. We show how this formalism should be extended in the case of each screening mechanism. Finally, we cross-check and validate our results in terms of the 2-point statistics against state-of-the-art cosmological N-body simulations that allow us to assess the LPT predictions in the nonlinear, quasi non-linear and baryon acoustic oscillation peak regimes. Comparison with simulations is required to ensure that the predictions will be sufficiently robust for upcoming surveys such as DESI, Euclid, LSST and WFIRST.

We note that a recent paper, submitted to the arXiv while this work was being finalized, also studied biased tracers in MG using CLPT (Aviles et al., 2018) but used a different bias scheme and results were not compared against full simulations.

Our paper is structured as follows: in Section II.1, we present the MG models we studied and the details of the N-body simulations employed to test our LPT implementation. We summarize the formalism for developing biased tracer statistics for GR in Section II.2. In Section III.1, we discuss the modifications required to the perturbative schemes to predict the real and Fourier space 2-point statistics and the associated bias parameters for tracers of different masses for models beyond GR. The LPT predictions are compared to statistics derived from simulations in Section III.2. Finally, we conclude, and discuss implications for future work, in Section IV. More detail on the derivations in the paper are given in the Appendix.

II Formalism

II.1 Modified gravity models

II.1.1 The $f(R)$ model

Despite the latest surge in the field of MG, deformations to GR, together with associated experimental tests, have been around for almost as long as GR itself (Will, 2006). One of the oldest attempts consisted of adding a nonlinear function $f(R)$ of the Ricci scalar $R$ to the standard Einstein-Hilbert action, the so-called $f(R)$ theories (De Felice and Tsujikawa, 2010), with a resulting action $S$ of the form:

[TABLE]

with $\mathcal{L}_{m}$ denoting the matter sector Lagrangian and G the gravitational constant. Since such an action frees up an additional degree of freedom, the latest interest in this class of theories comes from the idea that a modification of this type is responsible for cosmic acceleration, rather than dark energy (Carroll et al., 2004).

Such a model is the Hu-Sawicki $f(R)$ model (Hu and Sawicki, 2007), which we study in this paper and the functional form of which is given by:

[TABLE]

In equation (2), $\Omega_{m0}$ denotes the matter fractional energy density and $H_{0}$ the Hubble Constant, both evaluated today, $m=H_{0}\sqrt{\Omega_{m0}}$, which has dimensions of mass and $n$, $c_{1}$ and $c_{2}$ are free parameters.

In order to match the $\Lambda$CDM expansion history and for sufficiently small values of $|f_{R_{0}}|$, the background value of the Ricci scalar, $\bar{R}$, becomes equal to:

[TABLE]

where $\Omega_{\Lambda 0}$ is the dark energy fractional density evaluated at the present time.

The derivative $f_{R}=\frac{df(R)}{dR}$ becomes a functional of the cosmological parameters when evaluated today, through the relationship

[TABLE]

The above mapping allows us to reduce the number of free parameters and the Hu-Sawicki model is commonly parametrized by quoting the values chosen for $n$ and $|f_{R_{0}}|$, with the latter being the background value of the fifth potential evaluated today.

The reason the Hu-Sawicki model is so popular is that it can be cast into the form of a scalar-tensor theory that realizes the chameleon screening mechanism (Khoury and Weltman, 2004a, b), through a conformal transformation (Brax et al., 2008), with the quantity $f_{R}$ identified as the scalar field that is coupled to matter. Through the interplay between matter and the self-interaction potential, the scalar chameleon field becomes massive near high over-densities and the associated fifth forces get exponentially suppressed due to the Yukawa effect. A stronger screening effect is manifested in lower values of the parameter $|f_{R_{0}}|$, with the limit of $|f_{R_{0}}|\rightarrow 0$ exactly recovering GR. Following the literature and also because of the available simulations (as we will discuss below in II.1.3), we choose $n=1$ and consider three different $f(R)$ models with $\mathinner{\!\left\lvert\bar{f}_{R_{0}}\right\rvert}=\{10^{-6},10^{-5},10^{-4}\}$, which we shall refer to, from now on, as F6, F5 and F4.

When considering perturbations around a homogeneous and isotropic FRW metric in the conformal Newtonian gauge, the resulting system of the Poisson and KG equations becomes (Hu and Sawicki, 2007):

[TABLE]

where $\Phi_{N}$ is the Newtonian potential and $\delta\rho_{m}$ the matter density perturbation. $\delta R$, the perturbation to the Ricci scalar, can be written as a function of the scalar field $f_{R}$, which will play a central role in our screening implementation in Section II.2.

II.1.2 The $n$DGP model

The second MG model under consideration comes from the realm of higher-dimensional braneworld cosmology. The simplest example of such a configuration is the so-called Dvali-Gabadadze-Porrati (DGP) model (Dvali et al., 2000), which is described by an action of the following form:

[TABLE]

In this model, the spacetime consists of 5 dimensions, rather than the usual 4, but the standard model fields are restricted to a 4-dimensional (4D) brane and the free parameter, $r_{c}$, denotes the length-scale below which gravity becomes 4D. $R_{5}$ and $g_{5}$ denote the corresponding 5D versions of the Ricci scalar and metric tensor determinant, respectively. The resulting Friedman equation from (6) is:

[TABLE]

where $\epsilon=\pm 1$. Each of the two values of $\epsilon$ represents a particular branch of the model, with $\epsilon=+1$ producing the self-accelerating solution, which, however, suffers from undesirable ”ghost” instabilities (Koyama, 2007) and is thus an unphysical model to consider. The value of $\epsilon=-1$ corresponds to the so-called normal branch, hereafter called $n$DGP, which is well behaved, but does not self-accelerate and can match a $\Lambda$CDM expansion history only in the presence of dark energy.

When focusing on the normal branch and in the quasi-static limit for sub-horizon scales, the perturbations in the conformal Newtonian gauge give the modified Poisson system of equations (Koyama et al., 2009):

[TABLE]

with the coupling $\beta$ given by

[TABLE]

The $n$DGP model is a typical example of a scalar-tensor theory that realizes the Vainshtein screening mechanism (Vainshtein, 1972; Babichev and Deffayet, 2013), in which the modifications to gravity are suppressed in the existence of large second derivatives of the scalar field. In the second equation of (8), in particular, it can be seen how, once the second derivatives of the scalar field become large in high densities, the second term on the right hand side (r.h.s.) becomes significant and effectively weakens the source strength, resulting in strong screening of the fifth forces.

We consider the $n$DGP model for two choices of the free parameter $r_{c}$, corresponding to $n\equiv H_{0}r_{c}=1$ and $n=5$, which we shall call, from now on, N1 and N5, respectively.

II.1.3 N-body Simulations

Accurate realizations of structure formation in the nonlinear regime can only be achieved by performing N-body simulations. In this paper, we test our results against two groups of state-of-the-art N-body simulations, that serve complementary purposes to each other, as explained below.

The first group of simulations, to which we shall refer as Group I from now on, are the ELEPHANT simulations, presented in (Cautun et al., 2018). These span the parameter space of both MG models we study: the F4, F5 and F6 $f(R)$ models and N1 and N5 $n$DGP models. The $f(R)$ simulations were performed using the ECOSMOG code (Li et al., 2012; Bose et al., 2017), while the $n$DGP ones (Hellwing et al., 2017) using the ECOSMOG-V version (Li et al., 2013; Barreira et al., 2015), both of which were based on the GR code RAMSES (Teyssier, 2002) and where suitably extended to integrate the scalar field KG equation for the corresponding models using adaptive-mesh-refinement techniques. The parameters describing the background $\Lambda$CDM cosmology are the best-fit ones given by the 9-year WMAP release (Hinshaw et al., 2013) and have the following values: $\Omega_{b}=0.046$, $\Omega_{cdm}=0.235$, $\Omega_{m}=0.281$, $\Omega_{L}=0.719$, $h=0.7$, $n_{s}=0.971$ and $\sigma_{8}=0.82$. $N_{p}=1024^{3}$ equal mass particles were placed in a simulation box with a side $L_{box}=1024$ Mpc/h and the density field was resolved in a $1024^{3}$ resolution grid. Furthermore, the simulations were initialized at redshift $z_{i}=49$ using the Zel’dovich approximation (Zeldovich, 1970) and evolved through $z_{f}=0$.

Gravitationally bound dark matter halos were identified using the ROCKSTAR halo finder (Behroozi et al., 2013). Finally, so as to get an estimate of the variance, each model was simulated for 5 random realizations, corresponding to different random phases in the initial density field.

The $1024^{3}$ $(Mpc/h)^{3}$ volume simulations’ results in Group I become noisy at scales $r>100$ Mpc/h. To probe the BAO scales, where LPT has been previously found to perform very well for GR (Vlah et al., 2015a), we also test our results against the largest volume $f(R)$ simulations performed to date for the modified gravity lightcone simulation project (Arnold et al., 2018). In these simulations, which we will call Group II from now on, the box side is $L_{box}=1536$ Mpc/h with $2048^{3}$ equal mass particles used, for GR and the $\mathinner{\!\left\lvert\bar{f}_{R_{0}}\right\rvert}$=$10^{-5}$ model. The parameters describing the background $\Lambda$CDM cosmology are the best-fit ones given by the Planck collaboration (Planck Collaboration et al., 2016) and have the following values: $\Omega_{b}=0.0486$, $\Omega_{m}=0.3089$, $\Omega_{L}=0.6911$, $h=0.6774$, $n_{s}=0.9667$ and $\sigma_{8}=0.8159$. The simulations were performed using the MG code MG-GADGET (Puchwein et al., 2013), which is a MG extension to the code P-GADGET3, an improved version of the code GADGET-2 (Springel, 2005), created for GR simulations. Dark matter halo catalogues were produced using the SUBFIND algorithm (Springel et al., 2001). Each model has been simulated for one random realization.

For more detailed discussions on the N-body implementations, we refer interested readers to the corresponding publications.

II.2 Convolution Lagrangian Perturbation Theory for biased tracers in MG

II.2.1 LPT for dark matter

The Lagrangian Perturbation Theory approach to structure formation has been extensively studied (Zeldovich, 1970; Buchert, 1989; Bouchet et al., 1995; Hivon et al., 1995; Taylor and Hamilton, 1996; Matsubara, 2008a, b; Carlson et al., 2013; Matsubara, 2015; Desjacques et al., 2018) in the context of $\Lambda$CDM scenarios. Opposite to the Eulerian picture, in which one monitors the evolution of the desired quantities at a given, fixed, position, in LPT one instead tracks down the evolution of a given fluid element over time. Starting from an initial, Eulerian, comoving position $\mathbb{q}$ at a desired early time $t_{0}$, each mass element is mapped to its comoving Lagrangian position $\mathbb{x}(\mathbb{q},t)$ at time $t$, through the relationship

[TABLE]

The Lagrangian displacement $\mathbb{\Psi}(\mathbb{q},t)$, taken to be 0 at the initial time $t_{0}$, is the fundamental quantity of interest in LPT. Furthermore, enforcing mass conservation, through the continuity equation, between the initial and final infinitesimal volume elements centered around $\mathbb{q}$ and $\mathbb{x}$, respectively, gives $\rho_{m}(\mathbb{x},t)d^{3}x=\rho_{m}(\mathbb{q},t_{0})d^{3}q$. Assuming $t_{0}$ refers to an epoch early enough that the density perturbations around the background density $\bar{\rho}$ are negligible, meaning $\rho_{m}(\mathbb{q},t_{0})=\bar{\rho}_{m}$, allows us to obtain the dark matter fractional overdensity, $\delta_{m}$, in the Lagrangian picture:

[TABLE]

with $\delta_{D}$ being the Dirac delta function and $J(\mathbb{q},t)$ the determinant of the deformation matrix

[TABLE]

For an irrotational flow, which is a good approximation for cold dark matter and assuming that the gravitational evolution is governed by GR, perturbations around a flat FRW metric give the geodesic and Poisson equations, in the quasi-static approximation and for sub-horizon scales, as:

[TABLE]

We should point out that in (13) $\psi(\mathbb{x},t)$ denotes the metric perturbation, which should not be confused with the Lagrangian displacement field $\mathbb{\Psi}(\mathbb{q},t)$.

In the LPT picture, we perturbatively expand $\mathbb{\Psi}$ as

[TABLE]

and equations (10)-(14) form a closed system, that is recursively solved for the various orders of $\mathbb{\Psi}$. The first order solution is the so-called Zel’dovich approximation (Zeldovich, 1970).

In MG theories, as also explained in Section II.1, an additional degree of freedom is present, that directly couples to matter and causes particles to deviate from the nominal geodesics of GR. Consequently, equations (13) are, in principle, modified for a scalar-tensor theory and so is the LPT framework presented above. In (Valogiannis and Bean, 2017), the LPT approach was expanded for chameleons and symmetrons, including the first order contribution to the Klein-Gordon equation and was shown to perform very well in the context of the COLA hybrid framework. The LPT approach for scalar-tensor theories up to third order was presented in (Aviles and Cervantes-Cota, 2017), the main results of which we summarize in this section.

In (Bellini and Zumalacárregui, 2015), which developed an SPT framework for studying MG theories with a screening mechanism in the nonlinear regime, based on the closure theory approximation in (Taruya and Hiramatsu, 2008), the scalar field KG equation for a Brans-Dicke-like (BD) theory with interactions of a scalar field, $\phi$, was written as:

[TABLE]

where $\mathcal{I}(\phi)$ denotes the perturbative form of the field self-interaction term:

[TABLE]

where we adopted the standard notation $\mathbb{k}_{ijk}=\mathbb{k}_{i}+\mathbb{k}_{j}+\mathbb{k}_{k}$ and $M_{1}(\mathbb{k},t)$, $M_{2}(\mathbb{k}_{1},\mathbb{k}_{2})$ and $M_{3}(\mathbb{k}_{1},\mathbb{k}_{2},\mathbb{k}_{3})$ are mass terms. The higher order piece in the Fourier space representation of the interaction term in equation (16), incorporates the screening effect, up to third order, that is responsible for recovering GR at small scales. The mapping between a given scalar-tensor theory and the BD form above can be easily performed through assigning appropriate values to the mass terms and the BD coupling $\omega_{BD}$ above, as we will later show for our two models of study.

The perturbed modified Einstein equations have the form (Aviles and Cervantes-Cota, 2017)

[TABLE]

where we use (10) and introduce the time derivative operator $\mathcal{\hat{T}}=\frac{d^{2}}{dt^{2}}+2H\frac{d}{dt}$, as in (Matsubara, 2008a). The last term in the second line of (17), called $\it{frame}$-$\it{lagging}$ in (Aviles and Cervantes-Cota, 2017), is a geometrical term that occurs due to the fact that, in LPT, the KG equation should be expressed in Lagrangian $\left(\nabla^{2}\right)$, rather than Eulerian coordinates $\left(\nabla^{2}_{\mathbb{x}}\right)$. Taking this into account, equations (15)-(17) are combined to give

[TABLE]

where, following the definitions in (Bellini and Zumalacárregui, 2015), we have

[TABLE]

and all the quantities are Fourier transforms in the Lagrangian $\mathbb{q}$-space. The inverse Jacobean in (LABEL:eq:finpsi) reflects the derivative transformation to the $\mathbb{q}$-space, where the Einstein notation is adopted. Furthermore, $\delta(\mathbb{k})$, $\delta\mathcal{I}(\mathbb{k})$ and $\left(\nabla^{2}_{\mathbb{x}}\phi-\nabla^{2}\phi\right)(\mathbb{k})$ are the Lagragian Fourier transformations of the Lagrangian-transformed overdensity (11), the higher order interaction Kernels in (16) and the frame-lagging Kernel, correspondingly. The expression for the latter is given in (Aviles and Cervantes-Cota, 2017). Equation (LABEL:eq:finpsi) forms a closed system with (11), (12), (15) and (16) that is solved, perturbatively, to obtain the MG solution up to various orders in $\mathbb{\Psi}$, as in (14).

Solving for the first order solution, one gets (Valogiannis and Bean, 2017):

[TABLE]

which can be easily solved for the displacement field, as:

[TABLE]

We see that the r.h.s of (21) can be conveniently decomposed into a product of the first order density mode at very early times, $\delta^{(1)}(\mathbb{k},t=0)$, early enough to be gaussian and a space-time dependent growth factor $D^{(1)}(k,t)$, given by:

[TABLE]

In the GR limit, $A(k)=A(k=0)=4\pi G\bar{\rho}_{m}$, $D^{(1)}$ becomes scale independent and is nothing else than the first order growing mode for GR; the Zel’dovich approximation.

Moving on to the second order piece, we have:

[TABLE]

where, for compactness, we adopted the notation $\delta^{(1)}_{1}=(\mathbb{k}_{1},0)$ and the second order growth factor, $D^{(2)}(\mathbb{k}_{1},\mathbb{k}_{2})$, is given by (Aviles and Cervantes-Cota, 2017):

[TABLE]

The four individual components are given by:

[TABLE]

The two last terms represent the second order contributions to the growth factor, given by the screening and frame-lagging effects, correspondingly, while the expression for $K^{(2)}_{FL}$ is given in (Aviles and Cervantes-Cota, 2017). Despite its lengthier expression, when taking the GR limit we get $D^{(2)}_{FL}=D^{(2)}_{\delta\mathcal{I}}=0$ and $D^{(2)}_{a}=D^{(2)}_{b}$, allowing $D^{(2)}$ to become scale-independent, reducing to the known GR result, which can be well approximated by $D^{(2)}(t)=-\frac{3}{7}\left(D^{(1)}(t)\right)^{2}$ for $\Lambda$CDM cosmologies (Bouchet et al., 1995).

Solving for the third order piece in (LABEL:eq:finpsi) results in a lengthy differential equation for the third order MG growth factor, $D^{(3)}(\mathbb{k}_{1},\mathbb{k}_{2},\mathbb{k}_{3})$, that also needs to be symmetrized. The result is given by equation (LABEL:D3MG) in the appendix A.1. It should be also noted that equations (22), (25) and f(LABEL:D3MG) can be either solved by inverting the linear operator $\left(\mathcal{\hat{T}}-A(k)\right)$ using its Green function, as done in (Matsubara, 2015; Aviles and Cervantes-Cota, 2017), or by numerically solving the corresponding differential equations. Even though both methods give results that agree with each other well, we chose to proceed with the latter because it is computationally faster. The differential equations were solved using a $5^{th}$ order Runge-Kutta scheme, implemented in $\it{Mathematica}$ (Inc., ).

We finish this section by showing the particular expressions for the mass terms in (16) and the sources in (19) for the $f(R)$ and $n$DGP models we study. For the $f(R)$ Hu-Sawicki model, the Brans-Dicke scalar is simply $\omega_{BD}=0$, while the mass terms are given by the expansion (Koyama et al., 2009):

[TABLE]

which, using (2), gives (Aviles and Cervantes-Cota, 2017):

[TABLE]

In the case of the $n$DGP braneworld model, a similar procedure, informed by (8), gives the relevant expressions (Koyama et al., 2009; Aviles et al., 2018):

[TABLE]

with $\beta$ defined in (9). It is interesting to notice that, even though the interaction term in (8) contains only second order derivatives, in (LABEL:eq:nDGPmasses) a $3^{rd}$ order mass contribution is now present, that arises when transforming the Eulerian derivatives to the Lagrangian space through (12).

II.2.2 2-point statistics for biased tracers in GR

The perturbative theory of galaxy clustering (Desjacques et al., 2018), which aims to describe the statistics of biased tracers in the quasi-linear regime, consists of a perturbative description for the evolution of the underlying dark matter density field, combined with an analytical description for the bias parameters at each given order. In the case of cold dark matter, the calculation of the 2-point statistics, even in LPT, is a more straightforward process, since one just needs to plug the ${\bf q}$-space Lagrangian overdensity, mapped to the Eulerian frame through (11), into the common expressions for the autocorrelation function:

[TABLE]

and its Fourier space counterpart, the matter power spectrum

[TABLE]

When studying biased tracers, like for example dark matter halos, we need an analytical model to describe their statistical prevalance with respect to the underlying density field. Following (Matsubara, 2008b; Carlson et al., 2013), we employ a model of a local in matter density Lagrangian bias in which the positions of biased tracers are purely specified by a distribution of the underlying CDM density field $\delta(\mathbb{q},t=0)\equiv\delta(\mathbb{q})$, encoded through a function $F\left[\delta_{R}(\mathbb{q})\right]$, as

[TABLE]

where, consistent with the literature, we use the subscript $X$ to indicate biased tracers. $\delta_{R}(\mathbb{q})$ denotes the primordial density field smoothed over some scale $R$, while $\bar{\rho}_{X}$ is the mean density of tracers. Density conservation provides the equivalent of equation (11) for tracers,

[TABLE]

This model of local Lagrangian bias, which corresponds to a non-local bias in the Eulerian space, can be extended to include a biasing scheme that is non-local in the Lagrangian space (Matsubara, 2011). Combining (32) and (29) and after some transformations one gets the general expression for the 2-point correlation function for biased tracers in LPT,

[TABLE]

with

[TABLE]

where $\tilde{F}_{1}$, $\tilde{F}_{2}$ are the Fourier space representations of $F$, with corresponding wavemodes $\lambda_{1}$ and $\lambda_{2}$ and $\mathbb{\Delta}=\mathbb{\Psi}_{2}-\mathbb{\Psi}_{1}$. The notation $\delta(\mathbb{q}_{1})\equiv\delta_{1}$ has been adopted for all quantities. The ensemble average K in (34) can be cast into an exponent of a power series in cumulants, through the cumulant expansion theorem, $\langle e^{iX}\rangle=\exp\left[\sum_{N=1}^{\infty}\frac{i^{N}}{N!}\langle X^{N}\rangle_{c}\right]$, which, combined with a multinomial expansion, gives:

[TABLE]

in terms of a series of correlators

[TABLE]

where we adopted the commonly used notation for the Lagrangian cumulants in (36).

Keeping all terms in (35) that contain cumulants up to third order, which is the equivalent of the one-loop correction to the linear power spectrum, results in a highly oscillatory integrand that presents challenges when ensuring the integral is fully converged. (Matsubara, 2008a) proposed expanding all contributions to the exponent but the scale-independent “zero-lag” piece of $A_{ij}^{00}$, which results in a non-perturbative resummation scheme that is simpler to handle analytically. Building upon this result, (Carlson et al., 2013) proposed keeping all the terms of $A_{ij}^{00}$ in the exponent, in their Convolution Lagrangian Perturbation Theory (CLPT) scheme. Keeping only the linear component of $A_{ij}^{00}$ exponentiated, as done in (Vlah et al., 2015b, 2016), and performing the $\lambda$ and $k$ integrations in (33) gives a CLPT expression for the 2-point real space correlation function,

[TABLE]

with

[TABLE]

Furthermore, in (37) we define $U^{10}_{i}\equiv U_{i}$, $W^{000}_{ijk}\equiv W_{ijk}$ and use superscript numbers in brackets to indicate the various orders of contribution. The $1^{st}$ and $2^{nd}$ order Lagrangian bias factors, $b_{1}$ and $b_{2}$, are the expectation values of the $1^{st}$ and $2^{nd}$ order derivatives of the Lagrangian bias function $F$, respectively, through the identity (Matsubara, 2008b, 2011),

[TABLE]

In the case of dark matter, we have $F=1$ and $\tilde{F}(\lambda)=2\pi\delta_{D}(\lambda)$ (Carlson et al., 2013), and we recover $b_{1}=b_{2}=0$ for the unbiased, dark matter distribution.

The Fourier transform gives the CLPT power spectrum for biased tracers (Vlah et al., 2015b, 2016):

[TABLE]

In addition to the one-loop expressions for the two-point statistics, we also calculate the Zel’dovich ($1^{st}$ order LPT) approximation (White, 2014) for biased tracers in the configuration and Fourier space, which can be identified as the subset of terms in (37) and (40) that are linear in $P_{L}(k)$. These are the terms that depend on combinations of $\xi_{L}$, $U^{(1)}$ and $A^{L}_{ij}$.

While, for GR, CLPT does a very good job at modeling the configuration space $\xi(r)$, it is known to perform less well in reconstructing clustering in the Fourier space (Vlah et al., 2015a). Expanding the resummed exponent in (40) and performing the resulting integrals gives the Eulerian one-loop Standard Perturbation Theory (SPT) power spectrum for biased tracers in GR (Matsubara, 2008b),

[TABLE]

where

[TABLE]

is the 1D variance of the Lagrangian displacement and the functions $Q_{n}$ and $R_{n}$ were defined in (Matsubara, 2008b) for GR. The SPT power spectrum has been shown to follow the power spectrum much better than the CLPT prediction in GR (Vlah et al., 2015a).

The correlation function obtained from Fourier transforming (41) is, unfortunately, known to be ill-behaved (Matsubara, 2008a). However, if one performs an alternative resummation proposed in (Matsubara, 2008a, b), known as Lagrangian Resummation Theory (LRT), the resulting power spectrum,

[TABLE]

which differs from (41) only by the exponential prefactor, can be Fourier transformed to the configuration space and is found to characterize the BAO scales well for both dark matter (Matsubara, 2008a; Vlah et al., 2015a) and biased tracers (Matsubara, 2008b; Carlson et al., 2013). It decays sharply for large values of $k$ however.

II.2.3 Calculation of bias parameters in GR

In this section, we present an analytical model for the calculation of the bias parameters (39) in GR, which will be extended to include MG in section III.1.2. It should be noted though, that even in the complete absence of an analytical model, one or both of the bias parameters in CLPT can be treated as free parameters, to be fitted over simulations, as for example done in (Carlson et al., 2013; Reid and White, 2011). With regards to analytical models for bias, arguably the most popular one is the halo approach (Mo and White, 1996; Mo et al., 1997; Sheth and Lemson, 1999; Sheth and Tormen, 1999; Ma and Fry, 2000; Peacock and Smith, 2000; Scoccimarro et al., 2001; Cooray and Sheth, 2002), that is based on the extended Press-Schechter (PS) formalism (Press and Schechter, 1974; Bond et al., 1991), in combination with the Peak-Background Split (PBS) approach (Kaiser, 1984). A discussion of the accuracy of such approaches can be found in (Manera et al., 2010; Hoffmann et al., 2015). In what follows, we briefly summarize the main ingredients of this prescription in GR.

Let $M_{0}$ be the mass of a collapsed region at a redshift of interest $z$, that is enclosed in a Lagrangian region of radius $R_{0}$, which, given the mean matter density $\rho_{m0}$, will be given by

[TABLE]

The variance of matter density fluctuations in this region is,

[TABLE]

with $P_{L}(k,z=0)$ the linear matter power spectrum evaluated today and $W\left(kR_{0}\right)$ the top-hat smoothing Kernel,

[TABLE]

For GR, density perturbations are evolved in time, relative to present time using the GR linear growth factor $D(z)$.

Based on the PS theory and its variants (Press and Schechter, 1974; Bond et al., 1991), the comoving mean number density of halos per logarithmic mass bin $d\ln M$, $\bar{n}_{h}$, can be analytically modeled as:

[TABLE]

where $\bar{N}_{h}$ is the mean number of halos with mass $M$, in a bin of width $dM$, enclosed in a comoving volume $V$.

The quantity $\nu_{c}(M)$, the peak significance, is given by

[TABLE]

where $D(z)$ is the linear growth factor at the time of collapse $z$, normalized so that $D(z=0)=1$. In (48), $\sigma(M,z)=D(z)\sigma(M)$ is the variance at redshift $z$, with $\sigma(M)$ the variance (45) evaluated today and $\delta_{cr}$ is the critical overdensity for collapse at redshift $z$. For an Einstein De-Sitter (EDS) cosmology, the latter is always $\delta_{cr}=1.686$, which turns out to be a very good approximation for $\Lambda$CDM cosmologies and will be adopted here.

$f\left[\nu_{c}(M)\right]$ is the multiplicity function, that in the original PS theory is given by:

[TABLE]

The prescription (47), often referred to as the universal mass function, is exact in an EDS universe with a power law power spectrum. While (49) has been used to describe the halo mass function for a broad range of cosmologies, it lacks the necessary accuracy for precision predictions. For this reason, Sheth and Tormen (ST) (Sheth and Tormen, 1999), introduced an alternative function:

[TABLE]

where $A(p)=\left[1+\pi^{-\frac{1}{2}}2^{-p}\Gamma(0.5-p)\right]^{-1}$ and $q,p$ are free parameters that can be fitted over N-body simulations. The best fit pair was initially proposed to be $(q,p)=(0.707,0.3)$ which was later updated to $(q,p)=(0.75,0.3)$. These are considered to be the “standard” ST parameters (Sheth and Tormen, 1999; Sheth et al., 2001). For the PS function, $q=1,p=0$.

Based on the PBS argument (Kaiser, 1984), a large-wavelength density perturbation $\Delta$ (that is effectively constant on small scales) has the same effect on the formation of biased tracers as a modification to the mean background density by this offset. If by $\bar{n}_{h}(M,\Delta)$ we denote the halo mass function’s response to such a perturbation, also sometimes called the conditional mass function, then the fractional overdensity of halos will be given by (Mo and White, 1996; Mo et al., 1997; Sheth and Lemson, 1999):

[TABLE]

where $\bar{n}_{h}(M,0)=\bar{n}_{h}(M)$, is the standard, unconditional halo mass function. It also worth noticing that equation (51) also defines $F\left[\delta_{R}(\mathbb{q})\right]$, through (31). The Lagrangian bias of order $n$, is then given by

[TABLE]

where the time argument in the above is assumed and omitted for simplicity. Equation (52) is the rigorous definition of the bias parameters, that is exact even in the absence of an analytical description for the halo mass function and can be calculated numerically, for example employing separate universe simulations (Schmidt et al., 2013). In the presence of a universal mass function, the conditional mass function is given by the same expression (47), but with the modified peak significance $\nu^{\prime}_{c}(M)=\frac{\delta_{cr}-\Delta}{D(z)\sigma(M)}$, in which case the the bias factors are easily calculated (Sheth and Lemson, 1999) by (51), as:

[TABLE]

which are the so-called PBS Lagrangian biases. Applied on the ST function (50), the PBS biases are:

[TABLE]

These PBS biases are identical to the Lagrangian bias factors defined within the context of CLPT through (39) (Matsubara, 2008b, 2011).

III Results

III.1 Lagrangian Perturbation Theory for Biased Tracers in MG

III.1.1 2-point statistics for Biased Tracers in MG

In section II.2.2, we showed the expressions for the calculation of 2-point statistics in CLPT and its variants, under the assumption that gravitational evolution is governed by GR. Here we explain how each of these relationships have to be modified in the case of MG theories. We note that these results are consistent with those recently presented in (Aviles et al., 2018).

The two-point statistics for biased tracers in MG are given, as in GR, by the definitions 29 and 30, in the configuration and Fourier space, respectively. Considering biased tracers in the Lagrangian space, through (31), the overdensity of biased tracers in LPT is given by

[TABLE]

where we used density conservation. The above equation is similar to (32) for GR, but differs in that the Lagrangian field $\mathbb{\Psi}(\mathbb{q})$ follows the MG evolution presented in Section II.2.1. In particular, if we work in terms of the Fourier transform of $\mathbb{\Psi}(\mathbb{q})$, labeled as $\tilde{\mathbb{\Psi}}(\mathbb{p})$, the $n^{th}$ order LPT solutions in MG will be given by,

[TABLE]

where $\tilde{\delta}_{L}(\mathbb{p}_{n})$ are the linear-density Fourier transformed fields at the time of evaluation and the Kernels $L_{j}^{(n)}(\mathbb{p}_{1},..,\mathbb{p}_{n})$ are given by (Aviles and Cervantes-Cota, 2017):

[TABLE]

The MG growth factors in (57), are the ones given by (22), (24) and (LABEL:D3MG). Plugging the overdensity (55) into (29) and (30) and working as in (33)-(36), we arrive at equations (37) and (40), that give the 2-point statistics for biased tracers in MG and depend on the MG Lagrangian correlators

[TABLE]

For the MG correlators (58), we use the same definition and index structure as in GR, but these functions differ from their GR counterparts, because the quantities in the cumulants follow the MG LPT solutions (57). In particular, plugging (56) and (57) into (58), we get the Lagrangian correlators in MG:

[TABLE]

where we additionally performed the decompositions

[TABLE]

and defined $a(k)=P_{L}(k)+\frac{9}{98}Q_{1}(k)+\frac{10}{21}R_{1}(k)$. The functions $Q_{n}(k)$ and $R_{n}(k)$ that appear in the r.h.s of (III.1.1), are defined, as in GR, to be:

[TABLE]

The scale and redshift dependency of the growth factors alters the evaluation of these expressions relative to GR:

[TABLE]

and

[TABLE]

The functions $Q_{1}$-$Q_{13}$ in (LABEL:Qsmain) and $R_{1}$, $R_{2}$ in (LABEL:Rsmain) differ from GR as they depend on the MG growth factors (22), (24) and (LABEL:D3MG). In addition, the functions $\left[Q_{1}\right]_{MG}$, $\left[R_{1}\right]_{MG}$ and $\left[R_{1}(k)+R_{2}(k)\right]_{MG}$ are new ones that arise in MG. In the GR limit, that corresponds to $\bar{D}_{a}^{(2)}=\bar{D}_{b}^{(2)}=1$ and $\bar{D}_{FL}^{(2)}=\bar{D}_{\delta\mathcal{I}}^{(2)}=0$, $\left[Q_{1}\right]_{MG}=Q_{1}$, $\left[R_{1}\right]_{MG}=R_{1}$, $\left[R_{1}(k)+R_{2}(k)\right]_{MG}=R_{1}(k)+R_{2}(k)$ and the functions $Q_{n}$ and $R_{n}$ reduce to their GR expressions in (Matsubara, 2008a). In that limit, furthermore, the correlators (III.1.1) recover their GR forms presented in (Carlson et al., 2013). The derivations of the above, along with a more detailed discussion, are presented in Appendices A.1 and A.2.

In Figure 1, we show the contributions of the different terms in (37) as a function of $r$ for the F4 model (which predicts the largest deviations from GR), evaluated at $z=0.5$.

As in the GR case, proceeding to expand the resummed exponent in (40) and performing the resulting integrals, as shown in Appendix A.3, gives the equivalent of the Eulerian one-loop SPT power spectrum for biased tracers in MG:

[TABLE]

Equation (LABEL:PkXSPT), that depends on the functions (61), is the MG version of (41). Appendix A.1 provides a more thorough discussion. While we refer to this as the “SPT” expression, we note that, unlike in GR, where equation (41) has been shown to be identical to the SPT expression, in MG, additional terms that appear when transforming the Klein-Gordon equation from Eulerian to Lagrangian coordinates need to be considered to show the equivalence Aviles et al. (2018).

The LRT power spectrum for MG theories is obtained, just like in GR, by keeping the zero-lag term exponentiated:

[TABLE]

The derivation is discussed in Appendix A.3. The configuration space counterpart, $\xi^{LRT}_{X}(r)$, is obtained by Fourier transforming (65),

[TABLE]

with $j_{0}(kr)$ the zeroth-order Bessel function.

To evaluate the expressions (37), (40), (LABEL:PkXSPT) and (65) we modified a publicly available code released by (Vlah et al., 2016) in 111https://github.com/martinjameswhite/CLEFT_GSM, that efficiently performs the 2D integrals in (37), (40) using Haskel transformations, as well those in (III.1.1) and (61). On top of the functions $Q_{1}$-$Q_{13}$ and $R_{1}$, $R_{2}$, the code was extended to evaluate the new functions $\left[Q_{1}\right]_{MG}$, $\left[R_{1}\right]_{MG}$ and $\left[R_{1}(k)+R_{2}(k)\right]_{MG}$, as well as the modified correlators (III.1.1). We make this code publicly available in 222https://github.com/CornellCosmology/bias_MG_LPT_products. Our modified version accepts the modified gravity model growth factors, $D_{MG}(k,z)$ as input along with the linear power spectrum given by:

[TABLE]

The linear power spectrum for the background $\Lambda$CDM cosmology is calculated using the publicly available code CAMB (Lewis et al., 2000).

After calculating the necessary MG growth factors using our $\sc{Mathematica}$ notebook, we feed our modified version of the code with tabulated values of the growth factors for the various values of $k$, $r$ and $x$ needed at a given cosmological redshift $z$. The PYTHON module computes the various $Q_{n}(k)$ and $R_{n}(k)$ functions through equations (LABEL:Qsmain) and (LABEL:Rsmain), which are then used to calculate the various components of the CLPT power spectrum $P_{X}(k)$. The k functions can then be simply combined to give the SPT and LRT power spectra, by equations (LABEL:PkXSPT) and (65), respectively. Finally, the modified C$++$ counterpart follows a similar procedure to compute the configuration space two-point correlation function given by CLPT, through (37). The procedure is explained in more detail in the Appendix A.2.

We finish this section by noting that, even though we restrict our model to the case of a local Lagrangian bias, our framework can be extended to include a curvature bias and/or corrections from EFT, as in (Vlah et al., 2016). For modified gravity theories with scale-dependent growth, a general expansion bias is not possible in principle, though for some theories, such as the $f(R)$, the effects of the fifth force can be perturbatively absorbed in terms of higher-order derivatives (Desjacques et al., 2018; Aviles et al., 2018).

III.1.2 Calculation of bias parameters in MG

In this section we turn to the final necessary ingredient to describe biased tracers, an analytical model for the calculation of the bias parameters in MG.

Central to the GR derivation in Section II.2.3, is the assumption that spherical collapse is independent of the exterior spacetime, which, in the case of GR evolution, is given by Birkhoff’s theorem. In MG theories, however, the additional degree of freedom violates Birkhoff’s theorem, which will have important consequences for the modeling of the halo mass function, as well as on its response to an external density perturbation.

The Press-Schechter formalism (47) relies upon the assumption that the linearly evolved critical overdensity, $\delta_{cr}$, is always constant at the time of collapse, for example $\delta_{cr}=1.686$ for an EDS evolution. This is not the case for MG due to the presence of the additional scalar field that generates the fifth forces. Following (Li and Efstathiou, 2012; Lombriser et al., 2013), if we define $1+\delta_{h}=y_{h}^{-3}$, the evolution of a spherically symmetric halo density perturbation, will be given by:

[TABLE]

where ′ denotes derivatives with respect to $\ln(a)$ and $G_{eff}$, the modified Newton’s constant, is given by

[TABLE]

For the $n$DGP braneworld model,(Schmidt et al., 2010)

[TABLE]

where $\beta(a)$ was defined in (LABEL:eq:nDGPmasses) and

[TABLE]

with the $n$DGP parameter, $n=H_{0}r_{c}$. Note that the fifth force modification does not depend on either the mass or the environment, a property of the Vainshtein mechanism, which means that the collapse barrier for this model is redshift and scale independent (Schmidt et al., 2010; Barreira et al., 2013). Thanks to this property, the halo biases for this model can be easily calculated by the corresponding GR expressions (54), with a different value for the $\it{constant}$ threshold $\delta_{cr}$. For the background cosmology of the Group I simulations, at $z=0.5$, we integrate equation (68) and find $\delta_{cr}$ to have values (linearly extrapolated at $z=0.5$) of $\delta_{cr}=1.571$ and $\delta_{cr}=1.657$ for the N1 ($n=1$) and N5 ($n=5$) cases, respectively.

For $f(R)$ models, for a collapsing sphere of mass $M$ and radius $R_{th}$, $E$ is given by

[TABLE]

with $\beta=\frac{1}{\sqrt{6}}$. Finally, $\frac{\Delta R}{R_{th}}$ is given by

[TABLE]

Here the overdensity related to $y_{h}$ is embedded on a longer wavelength environment with over density, $1+\delta_{env}=y_{env}^{-3}$. In models that realize the chameleon screening mechanism, like the $f(R)$ Hu-Sawicki, only a thin shell of a massive sphere contributes to the fifth force, with a fractional thickness $\frac{\Delta R}{R_{th}}$, that causes an enhancement given by $E$. It is this factor $E$ that causes the environmental dependence of spherical collapse in MG, through (68). When set equal to zero, we recover the standard GR solution, which is the one that describes the evolution of the environment, since on such a long wavelength perturbations $G_{eff}\approx G$ and we have:

[TABLE]

In light of the coupling between a collapsing halo and the background on which it evolves, as seen through (68), it is clear that in MG the collapse barrier, $\delta_{cr}$, that is constant in GR, should now be promoted to a function of both the mass $M$ and the environment overdensity $\delta_{env}$. For each choice of $M$ and $\delta_{env_{i}}$, equations (68) and (74) form a system of coupled differential equations that we solve simultaneously, as in (Li and Efstathiou, 2012; Lombriser et al., 2013). The critical overdensity is identified as the smallest value of the boundary condition $\delta_{h}(a_{i})$ at the initial scale factor $a_{i}=0.002$, that causes a singularity i.e. signifying the onset of nonlinear collapse, at the scale factor of interest, which is then evolved to that scale factor $a$ through (Li and Efstathiou, 2012; Lombriser et al., 2013):

[TABLE]

with $D^{(1)}$ the linear GR growth factor

In Figure 2, we show the critical density, $\delta_{cr}(z,M,\delta_{env})$ as a function of mass, $M$, for the three $f(R)$ models at $z=0.5$, as it was obtained for a variety of environmental densities, $\delta_{env}$. Because the gravitational strength, $G_{eff}$, is greater in the modified models this allows haloes to form more easily, translating into a lower value of $\delta_{cr}$ than the value in $\Lambda$CDM. $G_{eff}$ is scale dependent and tends towards the GR value on large scales (which would collapse into large mass haloes). GR is also recovered for highly screened models and high density environments. For this reason the critical threshold tends towards the GR value $\delta_{cr}=1.686$ for smaller values of $f_{R0}$, increasing values of $M$ and more positive values of $\delta_{env}$ (regions with larger screening) as found in (Li and Efstathiou, 2012; Lombriser et al., 2013). The deviations from GR become progressively less pronounced as we move from weaker (F4, top panel) to stronger (F6, lower panel) screening, as one would expect.

Having calculated the function $\delta_{cr}(z,M,\delta_{env})$ for the three $f(R)$ models, the MG halo mass function can be again given by the universal prescription (47), but now with a modified peak significance

[TABLE]

The growth factors in MG are scale-dependent, meaning that in equation (76) one should in principle use $D^{(1)}_{MG}(k,z)$, however, it has been shown (Li and Efstathiou, 2012; Lombriser et al., 2013), that it is sufficient to use the GR growth factor, $D^{(1)}(z)$ to define the $\nu_{c}$ parameters. Modifications beyond this assumption are accounted for later in the free parameters of the halo mass function.

The dependence on the mass and the environment alters, and makes environment dependent, not only the unconditional halo mass function in MG, but also the conditional one. In the presence of a long-wavelength density perturbation, $\Delta$, the conditional halo mass function will be again described by (47), but now with a modified $\nu^{\prime}_{cMG}$ given by:

[TABLE]

Application of the bias definition (52) in that case, gives the first and second order bias parameters

[TABLE]

In the case of a sample of halos in a mass range of width $dM$ around a single value $M$, the conditional and unconditional mass functions need to be first averaged over the mass range,

[TABLE]

with

[TABLE]

The details of the peak-background split derivations, for (LABEL:pbsbiasMG), (LABEL:pbsbiasMGavg), are shown in more detail in Appendix B.

We note that another popular method, as an alternative to PBS, for calculating halo biases is the excursion set approach (Bond et al., 1991). Here the universal halo mass function is associated with the Brownian-walk, first crossing distribution of a collapse threshold. In the GR case, the redshift and scale independence of the collapse barrier leads to an analytical solution that recovers the common PBS biases. Given the potential scale, environment and redshift dependence that $\delta_{cr}$ has in MG models however, there is no analytical solution for the excursion set approach in MG and one would need to perform numerical simulations, rather than analytic prediction available for PBS, to determine the predicted biases (Li and Efstathiou, 2012; Lam and Li, 2012; Li and Lam, 2012).

In order to make predictions to compare against simulations, given that the correlation statistics sample a distribution of environments, rather than a specific value of $\delta_{env}$, we average all environment-dependent quantities over a probability distribution for environments in which halos form and reside, defining these on a fixed scale $\zeta$ which we set to 8 Mpch/$h$ (Lam and Sheth, 2008; Li and Lam, 2012; Lam and Li, 2012),

[TABLE]

where $\beta\equiv{\left(\frac{\zeta}{8}\right)^{3}}/(\delta_{cr}\sigma_{8}^{(2/\omega)})$ with $\omega\equiv\delta_{cr}\frac{n_{s}+3}{3}$.

The list of environments over which we average the other dependent quantities for the rest of this work, is $\delta_{env}=\left[-1,-0.72,-0.43,-0.15,0.13,0.42,0.7,0.98,1.27,1.55\right]$.

In Figure 3 we demonstrate the impact of the reduced $\delta_{cr}$ values in modified gravity models, separately from modifications to the halo mass function itself, by plotting the environmentally averaged bias parameters $b_{MG}^{1}$ and $b_{MG}^{2}$ for GR and three $f(R)$ models we study calculated from relations (54) and (LABEL:pbsbiasMG), while assuming the same underlying halo mass model, with $(q,p)=(0.707,0.3)$ for $\it{all}$ models. As described earlier, the increased gravitational strength in modified gravity models, parameterized by $G_{eff}$ in (69), allows haloes of a given mass to collapse more readily, yielding a lower critical threshold values for $\delta_{cr}$ and a lower bias relative to the background (reduced $b_{1}$ relative to GR). The deviations from GR are most pronounced in the models with least screening, like the $F4$ model and are suppressed in the presence of stronger screening.

To accurately predict the biases in the various modified gravity models, we also need to accurately characterize the halo mass function. The values $(q=0.707,p=0.3)$ used to characterize the ST halo mass function in GR were fixed by fitting to $\Lambda$CDM simulations (Sheth and Tormen, 1999) and would not be expected to predict the mass function for modified gravity theories, given the different physics involved in the growth rate and collapse of nonlinear structures. Since the form of the halo mass function is critical to evaluating the biases for the LPT correlation function and power spectra predictions, we determine the best fit values for $(q,p)$ from the simulated halo mass functions for each model in the mass ranges considered. These then uniquely determine the predictions for the biased tracer correlation and power spectra. This approach minimizes the errors that would be introduced at the outset of the LPT modeling from an inaccurate halo mass function. To do this, we evaluate the environment averaged Sheth-Tormen halo mass function (47) for various pairs values of $(q,p)$, each using the MG prescription (76), over the distribution of environments (LABEL:probenv) and identify the pair of values that best fits the corresponding halo mass functions from the simulations, through the simple criterion that minimizes the quantity

[TABLE]

previously used to fit the halo mass function in Galileons (Barreira et al., 2013, 2014). In (82), $i$ is the number of mass bins over which the sum is performed, which, can be tuned to model a narrow mass range rather than fitting the whole range with a single set of parameters.

In Figure 4, we plot the halo mass functions for the GR and three $f(R)$ MG models we consider at $z=0.5$, together with the best fit ST halo mass functions obtained through (82). In the case of $f(R)$ gravity, we plot both the mean ST halo mass function, as well as the halo mass function for each value of $\delta_{env}$, to give a sense of the variation among various different environments.

In Table 1, we show the best fit values of $(q,p)$ for the different models, mass ranges and simulations considered. For Group I simulations, we study predictions at $z=0.5$ for halos in a mass range $(2-3.5)\ \times 10^{12}M_{\odot}/h$, for all models. For the Group II simulations we analyze a $z=1$ snapshot, where we consider halos in three separate mass bins: a lower mass bin of $9\times 10^{11}-2\times 10^{12}\ M_{\odot}/h$, an intermediate bin of $5\times 10^{12}-1\times 10^{13}\ M_{\odot}/h$ and a higher mass bin, $1.1\times 10^{13}-9\times 10^{13}\ M_{\odot}/h$.

In Figure 5 we present the predicted correlation function, $\xi_{env}$ for different environments for the F6 model, along with the correlation function for the environment averaged bias values. Changing the environment can have a notable effect on the predicted correlation function, with variations of $\sim 10\%$ for the F6 model presented in the figure. A decrease in the background environmental density corresponds to a reduction in the correlation function because lower values of $\delta_{env}$ give rise to a lower $\delta_{cr}$ and consequently a lower value for the linear bias $b_{1}$.

III.2 Comparison with simulations

In this Section, we compare the performance of the various LPT resummation schemes under consideration, combined with our bias model, against the corresponding results from the Group I and Group II simulations, discussed in Section II.1.3, with respect to the correlation function and the power spectrum. All correlation functions from the simulations have been calculated employing the publicly available code $\sc{CUTE}$ (Alonso, 2012), using 30 linearly spaced bins in the range $0-140$ Mpc/$h$. The power spectra, on the other hand, have been extracted using a Cloud-In-Cell (CIC) mass assignment scheme, on a grid with resolution of $N_{grid}=1024^{3}$ and $N_{grid}=1200^{3}$ cells, for the Group I and Group II simulations, respectively. The power was binned in 30 linearly spaced points in the $k$ range of $0.008-0.3$ h/Mpc.

As discussed in Section III.1.2, for Group I simulations, we study predictions at $z=0.5$ for halos in the mass range $(2-3.5)\ \times 10^{12}M_{\odot}/h$, for all models. For the Group II simulations we analyze a $z=1$ snapshot, where we consider halos in three separate mass bins: a lower mass bin of $9\times 10^{11}-2\times 10^{12}\ M_{\odot}/h$, an intermediate bin of $5\times 10^{12}-1\times 10^{13}\ M_{\odot}/h$ and a higher mass bin, $1.1\times 10^{13}-9\times 10^{13}\ M_{\odot}/h$.

III.2.1 Correlation function

To assess the accuracy and performance of the LPT predictions, we first utilize the Group II simulations as a comparison dataset. Given that we only have 1 realization available for the F5 and GR models, the correlation functions exhibit noise due to the random initial phase, but still facilitate valuable conclusions about the performance of the methods tested given the appropriate combination of large volume and high resolution. These allow us to evaluate the performance of CLPT and its variations simultaneously across a wide range of scales, including both the BAO scales and the region roughly following a power-law scaling relation, down to $r\sim 5$ Mpc/$h$.

The LPT predictions use the PBS biases evaluated from the best fit halo mass function, when fitted over the specific mass ranges, for each bin as summarized inTable 1). Both the LPT and simulation results are compared to the Zel’dovich prediction for the correlation function.

To benchmark our findings we consider both the GR simulations as well as those for the modified gravity model. Figure 6 shows these comparisons for the three mass bin ranges, and Figure 7 shows the fractional variations with respect to the Zel’dovich component of CLPT.

For both the GR and modified gravity cases it is important to carefully understand the form of the halo mass function to get accurate LPT predictions. We find that simply adopting the standard ST pair of values, $(0.75,0.3)$, gives a poor approximation to the first order bias $b_{1}$ (for the various values of halo mass), consistent with the findings of the simulation creators (Arnold et al., 2018) when they extracted the bias estimate from the simulations and compared it to a standard ST prediction. For the results with the best fit halo mass parameters, we find that the full CLPT results for both the GR and the F5 model, incorporating the bias parameters evaluated using our PBS model (Table 1) and the environment averaging where necessary, does a very good job, in describing the power-law correlation function, $20-80$ Mpc/h, for all three mass bins and significantly improves upon linear theory at the BAO peak. For all three mass ranges, shown in Figure 7, the simulated correlation function falls below the Zel’dovich approximation and we find that the CLPT predictions are reflecting this better than both the linear and LRT predictions. The three approaches, CLPT, LRT and Zel’dovich, all perform well in characterizing the BAO peak for the low mass and intermediate mass bins, for the largest differences being in the highest mass bin for F5, where the LRT approach performs slightly better. The LRT performs poorly at the smaller scales, under $20Mpc/h$, significantly overshooting the observed correlation function, consistent with the results reported in previous studies performed on $\Lambda$CDM cosmologies (Vlah et al., 2015a; Matsubara, 2008b; Carlson et al., 2013; Reid and White, 2011). The Zel’dovich approximation provides the best agreement at scales below $10Mpch/h$.

It is also interesting to notice that, in Figure 6, while the correlation functions have similar values for GR and the F5 model in the lowest mass bins, the F5 result is noticeably lower than the GR one for the highest mass bin. The amplitude of the correlation function depends on the interplay between the dark matter component (which has higher values in MG) and mostly the linear bias factor, $b_{1}$, which is lower for MG. In the lower mass bins, the combination of the above two is such that the difference between the GR and the F5 curves is almost neutralized, while in the highest mass bins the linear bias factor $b_{1}$ is, relatively, even lower, causing the F5 two-point function to have clearly lower values than in the GR case.

To expand on the results from the large volume simulations, we now look into the comparison with the Group I simulations, that, while spanning a smaller volume, allow us to test our schemes for a wider range of models. For each model, five different random realizations are available and the error bars represent the standard deviations over these realizations. Starting with the $f(R)$ family, which is plotted in Figure 8, we see that the picture painted for the $F4$ an $F5$ models here is similar with the one for the $F5$ model (at $z=1$), with CLPT performing the best at scales $r>20$ Mpc/h and the Zel’dovich result being superior at capturing the smaller scales, while the trend is more pronounced in the F4 case. For the F6 model however, not only does CLPT perform better at these larger scales, but it seems to trace the simulation results more accurately compared to Zel’dovich down to r$\sim 7$ Mpc/h.

To explore this small scale performance sensitivity to screening in more detail note we consider what makes the linear order LPT and its one-loop extensions perform differently. In (Vlah et al., 2015a) (and also in (Tassev, 2014)), it was argued, in dark matter-only studies, that LPT does a poor job at estimating the higher order corrections to the linear displacement dispersion, given by (42) and the one-loop correction piece in LPT given by

[TABLE]

Comparison with simulations in (Tassev, 2014; Vlah et al., 2015a), found CLPT to overestimate the size of this correction at small scales, through (83), with the true value of the total $\sigma_{L}^{2}+\sigma_{1loop}^{2}$ being closer to $\sigma_{L}^{2}$, which is what the Zel’dovich result uses. Because one-loop CLPT strongly depends on these corrections, through its zero-lag terms (as can be seen in Appendix A.2), it performs worse at the smaller scales compared to its Zel’dovich counterpart. Calculated from our theory prediction for GR at $z=0.5$, the ratio $\frac{\sigma_{1loop}^{2}}{\sigma_{L}^{2}}=0.126$, close to the value in Figure 5 of (Vlah et al., 2015a). In comparison, for the $f(R)$ models at $z=0.5$, the ratio $\frac{\sigma_{1loop}^{2}}{\sigma_{L}^{2}}$ is $(F4,F5,F6)=(0.212,0.180,0.152)$. The higher values for F4 and F5 lead to an overestimation in these cases that is responsible for the worse behavior at smaller $r$. It Is worth noting here that if we do not include screening, the ratio $\frac{\sigma_{1loop}^{2}}{\sigma_{L}^{2}}=0.17$ in the $F6$ model, as opposed to the full value of $0.152$. For the $z=1$ Group II simulations, the ratios are $(GR,F5)=(0.08,0.103)$, which explains the lower discrepancy and better performance of CLPT for F5. This also is consistent with considering that this is an earlier reference in which clustering differences between the theories will be less pronounced. From a physical standpoint, the overestimation reflects an inability in LPT (including the Zel’dovich result), to trap dark matter particles within halos (Vlah et al., 2015a), which seems to be more apparent in the LPT approach for stronger MG chameleons. Fortunately, as we said earlier, these models are the ones that violate the observational constraints and are thus less interesting from an astrophysical standpoint.

Finally, we test our LPT approaches applied on the $n$DGP models, that represent the Vainshtein screening mechanism, and the correlation functions of which are presented in Figure 9, for all LPT schemes and the Group II simulations. Just like in the $f(R)$ models, CLPT does a very good job at describing the correlation function for large scales and beyond that, it even seems to perform at least equally well as the Zel’dovich curve down to scales $r\sim 10$ Mpc/h, similar to the F6 and GR cases in Figure 8 discussed earlier. The measurement of the 1-loop statistic discussed in the previous paragraph is consistent with this; for the $n$DGP models, the ratio$\frac{\sigma_{1loop}^{2}}{\sigma_{L}^{2}}=(0.129,0.122)$ for (N1, N5) respectively, very consistent with the GR value = 0.126. This is the case even in the weaker screening case, the N5 model, and is a very promising sign, given that the Vainshtein mechanism is highly efficient at screening modifications to gravity at smaller scales and contains viable candidates that self-accelerate (even though this particular model does not). The relative performance among the different resummation schemes is very similar to the one observed in the GR and $f(R)$ cases, with all LPT models improving the accuracy at the BAO peak upon linear theory, with the LRT scheme giving more power that CLPT and then the Zel’dovich result. The characterization of the BAO peak on scales $r>100$ Mpc/h is limited in the simulation box with side 1,024 Mpc/h; larger-volume simulations for the $n$DGP model, comparable to the Group II simulations for F5 or GR, will allow us to more clearly trace the region between $100-140$ Mpc/h and draw stronger conclusions about how our models perform at the BAO scales.

The fact that CLPT performs well for all modified gravity models considered in the power-law and BAO scales is very encouraging. On the smaller scales, its robustness for highly screened models is also a positive result. If an MG model was ever detected, it would be a highly screened case, given the tight constraints placed on GR; models F4 and F5 are actually ruled out by observations (Burrage and Sakstein, 2018), but we include them in our analysis to fully investigate the chameleon phenomenology with LPT.

III.2.2 Power spectrum

Complementary to the correlation function, we also perform tests on its Fourier space counterpart, the halo power spectrum. The mass bins and bias values used in the power spectra calculations are exactly the same as the ones presented in the correlation function case, but with the additional step that all power spectra are shot noise corrected (Seljak et al., 2009): $\tilde{P}(k)=P(k)-\frac{1}{n_{h}}$, where $P(k)$ is the uncorrected power spectrum, and $n_{h}$ is the number density of halos in each bin. Especially for the higher mass bins that contain less halos, this effect is not negligible, especially at higher $k$. We also identify the scale at which perturbation theory starts to fail, $k_{NL}$, with the vertical dashed-dot blue line, using the definition (Matsubara, 2008a), $k_{NL}=(2\sigma_{L}^{2})^{-1}$, with $\sigma_{L}^{2}$ the linear power spectrum dispersion defined in (42).

In Figure 10, we present the F5, Group II snapshot at $z=1$. We find the expanded, SPT power spectrum (LABEL:PkXSPT) to perform very well at capturing the small $k$ and to follow the power spectrum until $k\sim 0.25h/$Mpc, where it starts to overestimate power compared to the simulations, for all three mass bins. This behavior is consistent with what was found in the GR case in earlier works on (Vlah et al., 2015a; Matsubara, 2008b) and also for dark matter in MG (Aviles and Cervantes-Cota, 2017). The linear theory result is only accurate at very large scales and quickly underestimates the power at $k>0.05h$/Mpc. Unlike in the correlation function comparison, where LPT was found to perform very well at a wide range of scales, here we see that the CLPT power spectrum (40) decays quickly and performs considerably worse than the SPT expansion. This is not unexpected, since the power spectrum in LPT has been found to receive, unlike in the configuration space, significant contributions from large, nonlinear $k$ modes, where LPT performs poorly and fails to trap particles inside dark matter halos (Vlah et al., 2015a). Our results show that this to also the case in our MG models. We find that this effect is even more pronounced in the LRT power spectrum (65), which decays sharply in $k$-space, for this reason we do not include the result in our plots.

In Figure 11 summarizing the GR and $f(R)$ cases from the Group I simulations, at $z=0.5$, while the CLPT consistently underestimates the power spectrum for all models, the SPT result tracing the simulation points well for F5 and F6 until $k\sim 0.2h$/Mpc, at which it starts to overestimate the power spectrum. The performance for the the F4 model is slightly worse on small scales. This earlier deviation is not surprising given that the model has the smallest $k_{NL}$ prediction, resulting from a comparatively higher 1D linear dispersion.

For the two $n$DGP models, also shown in Figure 11, we find that the SPT predictions perform well at scales $k<0.15$ Mpc/h but overestimate the power on small scales. The CLPT predictions consistently underestimate the power, and are broadly comparable to the linear prediction.

IV Conclusions

In this work, we modeled the two-point statistics of biased tracers in modified gravity (MG) up to one-loop order in the linear power spectrum, using the Convolution Lagrangian Perturbation Theory (CLPT) framework and its variants. Following standard methods in the literature, the linear piece of the two-point Lagrangian correlator for dark matter is kept exponentiated in the expression for the two-point correlation function, but everything else is expanded, leading to a series of convolution integrals, the expressions of which we derive for scalar-tensor theories.

The evolution of the underlying dark matter density field is described by the LPT framework for dark matter, suitably extended to study scalar-tensor theories, along with an analytical model for the calculation of the first and second order bias parameters in MG. To perform the bias calculations, we employ the Peak-Background split (PBS) approach, in which biases are modeled rigorously as responses of the universal Sheth-Tormen halo mass function in the presence of a long-wavelength density perturbation. This is extended in MG theories, to account for the dependence of the gravitational collapse on the environment and screening. Our PBS implementation, provides a quantitative prediction for the increased production, and the related lower biases, for haloes of a given mass. We apply this scheme to the $f(R)$ Hu-Sawicki and the $n$DGP braneworld models, that are representatives of the chameleon and Vainshtein screening mechanisms, respectively. We make the code used for the analytic predictions publicly available in 333https://github.com/CornellCosmology/bias_MG_LPT_products and evaluate their performance against state-of-the-art cosmological N-body simulations, for a variety of MG models at $z=0.5$ and $z=1$, with respect to the correlation function and the power spectrum in a variety of mass regimes and scales.

The CLPT implementation, in combination with the analytical bias model, gives good agreement with the simulations, with the only free parameters necessary being those to best-fit the Sheth-Tormen universal halo mass function at the given mass range. The CLPT predicts the correlation function across scales $20-80$ Mpc/$h$, tracing the simulation results at an accuracy of $2-3\%$ and better. At the BAO scales, that provide a valuable probe of fundamental physics, CLPT was found to improve significantly upon the linear theory and Zel’dovich predictions for the F5 models, just like in GR. The Lagrangian Resummation Theory (LRT) approach improved the accuracy a little further at BAO scales for the highest mass range considered. At scales of $r<20$ Mpc/$h$, the CLPT performed well for the highly screened model $F6$ and for the $n$DGP models, while the Zel’dovich predictions performed better for the weakly screened F5 and F4 models. The reason for this behavior was identified, being an overestimation in these low-screening chameleon models of the one-loop contributions to the zero-lag terms at small scales.

In Fourier space, consistent with findings for GR, the CLPT power spectrum was found to underestimate power quickly, compared to the simulations for all MG models. This is due to the power spectrum receiving significant contributions from large $k$, where LPT performs poorly. The Standard Perturbation Theory (SPT) approach, though, which is the low-$k$ expanded version of this power spectrum, performs very well and remains consistent with the simulation results down to $k\sim 0.2h$/Mpc for the $f(R)$ models and down to $k\sim 0.15h$/Mpc for the two $n$DGP models. Beyond these scales, the SPT curve overestimates the power spectrum, as has been found for GR previously.

While we have focused our analysis on LPT predictions for real space, our model can be expanded to capture the redshift-space distortions required for upcoming LSS surveys. Furthermore, even though we focused on a local in matter density bias scheme in the Lagrangian space, in which the bias is purely a function of the local density, one can extend this to include other factors determining bias into the formalism, such as curvature bias, and model them successfully by this PBS scheme. The same applies for potential extensions to include EFT corrections to our LPT model, as in (Vlah et al., 2016), which could also be used to calculate the components of the Gaussian Streaming Model for MG theories. Finally, our CLPT MG framework can be used to analytically predict marked statistics in MG and assess their ability to boost the MG signals carried in cosmic density fields, as in (White, 2016; Valogiannis and Bean, 2018). We leave these natural extensions to future work.

In the coming decade, a wide array of cosmological surveys will span a large part of the observable universe, searching for hints of new physics beyond $\Lambda$CDM. In this work we demonstrate that semi-analytical approaches, extensively employed in the context of standard GR, can serve as invaluable tools to predict structure formation in cosmologies with an extra degree of freedom in the gravitational sector. A next step for these approaches are to confront them in comparison to realistic simulations of galaxies and clusters that will be observed with surveys coming online in the coming year or two and assess survey ability to identify and constrain potential deviations from GR.

Acknowledgments

We wish to thank Wojciech Hellwing for kindly making available the $n$DGP simulations, on behalf of (Hellwing et al., 2017) and Baojiu Li for kindly providing the $\sc{ELEPHANT}$ simulations, on behalf of (Cautun et al., 2018), for numerous discussions and for making us awaref of the Lightcone simulations. We are also grateful to Christian Arnold, as well as to all the other authors of (Arnold et al., 2018), for kindly providing their Lightcone simulations and for numerous discussions. We also wish to thank Martin White for useful discussions on available simulations of $\Lambda$CDM cosmologies. The work of Georgios Valogiannis and Rachel Bean is supported by DoE grant DE-SC0011838, NASA ATP grants NNX14AH53G and 80NSSC18K0695, and NASA ROSES grant 12-EUCLID12-0004.

Appendix

Appendix A One loop corrections for biased statistics in MG

In section III.1.1 of the main text, it was stated that the two point statistics for biased tracers, up to one loop in CLPT, are given by equations (37) and (40), in the configuration space and the Fourier space, respectively. The expressions are convolutional integrals over a sum of individual terms that depend on the Lagrangian correlators (58), which are essentially the fundamental blocks of CLPT. These expressions (58), however, are functions of the Lagrangian coordinates $\mathbb{q}$, while the LPT solutions for the displacement fields up to various orders are found in the Fourier space (through the growth factors (22), (25) and (LABEL:D3MG)). As was also stated in III.1.1, substituting the LPT solutions (56) into (58) gives the integral expressions for the MG Lagrangian correlators (III.1.1), that depend on the functions (61), which are also the building blocks of the SPT and LRT power spectra (LABEL:PkXSPT) and (65). We start with the innermost layer of integration, deriving the expressions for the k-dependent functions (61) in section A.1, before showing how to get to the correlators (III.1.1) in section A.2. Finally, in section A.3, we show how the SPT and LRT expressions for the two-point statistics are derived. The notation and index structure is the one adopted in (Matsubara, 2008a, b; Carlson et al., 2013). These results are consistent with those recently presented in (Aviles et al., 2018).

A.1 Polyspectra and k-functions in MG

In LPT, we solve for the displacement field $\mathbb{\Psi}(\mathbb{q})$ across various orders in perturbation theory, as

[TABLE]

In the Fourier space representation, $\tilde{\mathbb{\Psi}}(\mathbb{p})$, the solutions are expanded as

[TABLE]

where $\tilde{\delta}_{L}(\mathbb{p}_{n})$ are the linear density fields in the Fourier space at the time of evaluation. For gravitational evolution governed by GR, the growth factors are only functions of time, and under the additional assumption of an Einstein-De Sitter evolution, the Kernels $L_{j}^{(n)}(\mathbb{p}_{1},..,\mathbb{p}_{n})$ admit simple scale-independent expressions and (LABEL:psifour) can be simply evolved in time by powers of the linear growth factor (Matsubara, 2008a). This assumption gives results accurate at the sub-percent level for $\Lambda$CDM cosmologies (Bouchet et al., 1995).

In MG, however, the simple description presented above does not hold, in principle, because the growth factors depend on both space and time, as we saw in Section II.2.1. Following (Aviles and Cervantes-Cota, 2017), we define

[TABLE]

where the MG growth factors are calculated through the prescription described in Section II.2.1 and their time arguments have been dropped for notational simplicity. Furthermore, the subscript in the third order Kernel is meant to emphasize on the fact that the configuration that enters the 2-point statistics should be symmetrized (Matsubara, 2008a; Aviles and Cervantes-Cota, 2017). The symmetrized third order growth factor that enters the two-point statistics is given by (Aviles and Cervantes-Cota, 2017):

[TABLE]

with $D^{(3)}(\mathbb{k},-\mathbb{p},\mathbb{p})$ given by:

[TABLE]

In (LABEL:D3MG), we additionally defined

[TABLE]

and also used the kernels $K^{(2)}_{FL}$ and $K^{(3)}_{\delta\mathcal{I}sym}$, the forms of which were shown in (Aviles and Cervantes-Cota, 2017). The second and third order LPT kernels (LABEL:LPTKernels) need to be evaluated numerically, now, for each value of $\mathbb{p}_{1},\mathbb{p}_{2}$ and $z$ and so should (LABEL:psifour), which is the main point of divergence between the MG implementation and the corresponding one in GR.

Now, following (Matsubara, 2008b), we define the mixed polyspectra $C_{i_{1}...i_{N}}$, as

[TABLE]

It is useful to decompose the various polyspectra into the various constituents, order by order in perturbation theory, as e.g.

[TABLE]

where the additional notation has been adopted

[TABLE]

and as previously, the bracketed numbers in the superscripts indicate the orders of contribution in each $\tilde{\Psi}_{i}(\mathbb{p})$. The various polyspectra can be expressed as functions of the Lagrangian kernels (LABEL:LPTKernels), by identifying the different contributions across each order in LPT, as in (LABEL:polyorders) and plugging the solutions (LABEL:psifour) into equation (LABEL:mixedployorder). The ones relevant for the two-point statistics of biased tracers are (Matsubara, 2008b):

[TABLE]

where by $P_{L}(p)$ we denote the MG linear power spectrum.

The scalar functions $Q_{n}(k)$ and $R_{n}(k)$ that contribute to the SPT power spectrum (and as we shall see in the next section, to the Lagrangian correlators (58)), are expressed as functions of the polyspectra (LABEL:mixedployorder) in GR (Matsubara, 2008b). Fortunately, since in MG the above picture is only modified through the modified Kernels in (LABEL:LPTKernels), the relationships that give the various scalar functions are of the same form as the ones presented in (Matsubara, 2008b) (in particular, equations (A50)-(A67) in Appendix A). However, one should be cautious at this point, because certain symmetries that are present in the GR solutions, are not preserved anymore. In particular, the integral in the l.h.s of eq. (A59) in (Matsubara, 2008b) will not be equal to $R_{1}(k)+R_{2}(k)$ anymore, because of the scale and redshift dependence of the MG growth factors. In a similar manner, the functions resulting from eq. (A57) and (A61), that used to be equal to $R_{1}(k)$ and $Q_{1}(k)$, respectively, in GR, will differ for our MG models and should be additionally calculated. We denote these by $\left[R_{1}(k)+R_{2}(k)\right]_{MG}$, $\left[R_{1}(k)\right]_{MG}$ and $\left[Q_{1}(k)\right]_{MG}$ to emphasize on their GR limit. In order to illustrate how these calculations are performed, we show the derivations for $\left[R_{1}(k)+R_{2}(k)\right]_{MG}$ and $\left[Q_{1}(k)\right]_{MG}$, that are both new in MG and serve as a representative example of each category. For the former, we have:

[TABLE]

where we defined the quantities $x=\hat{\mathbb{k}}\cdot\hat{\mathbb{p}}$, $p=kr$ and

[TABLE]

as was done in (Aviles and Cervantes-Cota, 2017). Similarly, for $\left[Q_{1}(k)\right]_{MG}$ we will have:

[TABLE]

where $|k-p|=k\sqrt{1+r^{2}-2rx}$ is used and similarly as before, we defined

[TABLE]

The above two Kernels, as well as all of the rest that we need, can be compactly written in the form (Matsubara, 2008b):

[TABLE]

and

[TABLE]

Using similar methods as the one presented above we get that the various $Q_{n}(k)$, and after some algebra, are given by:

[TABLE]

Similarly, in accordance with equation (99), we will have the $R_{n}(k)$ functions:

[TABLE]

The functions $Q_{1}$-$Q_{3}$, $R_{1}$ and $R_{2}$ are the only ones that are necessary to calculate LPT statistics for pure dark matter considerations in MG, with the rest of them that we present here, being the additional functions needed for statistics of biased tracers in MG (in the configuration space). It should be emphasized at this point, that even the functions that have the same integral structure as in GR (for example, $Q_{9}$-$Q_{13}$), do differ from their GR values, but this difference is manifested in the MG linear power spectra that appear in the integral relations (98) and (99).

Let us finish this section, by noting that in the GR limit, the above functions can be rather easily shown to recover their standard GR forms given in (Matsubara, 2008b), if one keeps in mind that $\bar{D}_{\delta\mathcal{I}}^{(2)}=\bar{D}_{FL}^{(2)}=0$ and $\bar{D}_{a}^{(2)}=\bar{D}_{b}^{(2)}=1$ (for Einstein-De Sitter) in this limit.

A.2 Lagrangian correlators and q-functions in MG

Having derived the expressions for the scalar functions $Q_{n}(k)$ and $R_{n}(k)$ in MG cosmologies, we will now derive the integral formulas for the Lagrangian correlators (58), that are the building blocks of the 2-point statistics in CLPT. We will adopt the notation of (Carlson et al., 2013) in this section and will show how the functions in their Appendix B will change for our MG models. To illustrate how the connection between the functions (58) and the ones presented in the previous section is drawn and also to show how these calculations are performed, we pick a reprsentative example of one these functions, $U_{11}^{(2)}$ and show the derivation below. Starting with the definition:

[TABLE]

we plug in the Fourier space representation of the field $\Delta_{i}=\Psi_{i}(\mathbb{q_{2}})-\Psi_{i}(\mathbb{q_{1}})$, as well of the linear overdensities and get:

[TABLE]

This expression contains two terms, which turn out to be equal. For this reason, we focus on the first one and notice that the cumulant can be expressed as a polyspectrum, through (LABEL:mixedployorder). The substitution gives:

[TABLE]

In the last line, we defined the integral $I_{C}$, that can be calculated by using the definitions (LABEL:mixedployorderspef):

[TABLE]

where we made use of the previous result (LABEL:R12plus). Plugging the result (LABEL:Urelpart2) into (LABEL:Urelpart) and relabelling $p_{1}$ as $k$, we have

[TABLE]

where we made use of the Bessel function identity $\frac{1}{2}\int_{-1}^{1}dxxe^{ixkq}=ij_{1}(kq)$. In exactly the same way, the second term in (LABEL:Urelpart) is equal to the first, which finally gives:

[TABLE]

This is the MG equivalent of equation (B$28$) in Appendix B of (Carlson et al., 2013), which is obviously recovered in the GR limit. In a similar manner, but after lengthy calculations, we get the expressions for all the correlators:

[TABLE]

where we defined $a(k)=P_{L}(k)+\frac{9}{98}Q_{1}(k)+\frac{10}{21}R_{1}(k)$ and, following (Carlson et al., 2013), we decomposed the matter terms as

[TABLE]

Now that the basic framework has been laid out, let us finish this section by briefly summarizing the steps followed to calculate the two point statistics for a given model: after calculating the necessary MG growth factors, (22), (25) and (LABEL:D3MG), using our $\sc{Mathematica}$ notebook, we feed our modified version of the code in 444https://github.com/martinjameswhite/CLEFT_GSM with tabulated values of the growth factors for the various values of $k$, $r$ and $x$ needed at a given cosmological redshift $z$, as well as with the MG linear power spectrum given by:

[TABLE]

The linear power spectrum for the background $\Lambda$CDM cosmology is calculated using the publicly available code CAMB (Lewis et al., 2000). The PYTHON module computes the various $Q_{n}(k)$ and $R_{n}(k)$ functions through equations (LABEL:Qs) and (LABEL:Rs), which are then used to calculate the various components of the CLPT power spectrum $P_{X}(k)$, through the integrations (LABEL:qfuncs) and (40). To calculate (LABEL:qfuncs), the $q-$function module is modified accordingly. The k functions can then be simply combined to give the SPT and LRT power spectra, by equations (LABEL:PkXSPT) and (65), respectively. Finally, the modified C$++$ counterpart follows a similar procedure to compute the configuration space two-point correlation function given by CLPT, through (37).

A.3 SPT and LRT Power spectra

In the main text, it was stated that the SPT and LRT power spectra, given by (LABEL:PkXSPT) and (65), correspondingly, are produced when one expands the resummed terms in the exponent of relation (40). Here we show how this derivation takes place, a process that once again shares a lot of similarities with the corresponding one in GR. When fully expanding the resummed Lagrangian correlator $A^{L}_{ij}$ in (40), one gets:

[TABLE]

As done previously, we pick one of terms that is modified, that is the term $b_{1}^{2}ik_{i}U^{11}_{i}$, and perform the integration as an example. Plugging in $U^{11}_{i}$ from (LABEL:qfuncs):

[TABLE]

where we also used the identity $\int_{0}^{\infty}dqq^{2}j_{1}(pq)j_{1}(kq)=\frac{\pi}{2pk}\delta_{D}(p-k)$. Similar computations for the rest of the terms give (LABEL:PkXSPT). In the LRT case, we expand everything but the $q$-independent, “zero-lag” piece of $A^{L}_{ij}$, which is equal to $2\sigma^{2}_{L}\delta_{ij}$, with $\sigma_{L}^{2}=\frac{1}{6\pi^{2}}\int_{0}^{\infty}dkP_{L}(k)$. Since this term is scale-independent, it can be pulled out of the q integral and all the other integrations can be performed in the same manner as above, resulting in (65), that differs from (LABEL:PkXSPT) only in terms of the resummed exponential factor.

Appendix B PBS biases in MG

In this section, we will explain the derivation of the Lagrangian PBS biases (LABEL:pbsbiasMG) in MG models with environmentally dependent gravitational collapse. The PBS argument is commonly employed in conjunction with the halo approach (Mo and White, 1996; Mo et al., 1997; Sheth and Lemson, 1999), where one states that the conditional halo mass function $\bar{n}_{h}(M,\Delta)$, modulated by a long-wavelength density perturbation $\Delta$, is modeled by the universal mass function prescription (47), with the collapse threshold shifted as $\delta_{cr}\rightarrow\delta_{cr}-\Delta$. The same result was also derived in (Schmidt et al., 2013), based on the rigorous definition (52), following a separate universe approach: with regards to galaxy clustering, a large-scale density perturbation, $\Delta$, can be viewed as a modification of the mean physical density $\bar{\varrho}_{m}$ by an offset, that is (Desjacques et al., 2018)

[TABLE]

where $\bar{\varrho}_{m}$ should not be confused with the mean comoving density $\bar{\rho}_{m}$ and in a similar manner, the fractional overdensity at a point $\mathbb{x}$, $\delta_{m}(\mathbb{x})$, is shifted as:

[TABLE]

This reasoning can be employed to calculate the conditional halo mass function $\bar{n}_{h}(\Delta)$, by noticing that (47) depends solely on comoving quantities (that will not change), with the only exception of the density threshold $\delta_{cr}$, through the peak significance

[TABLE]

that quantifies how rare a fluctuation above the density barrier $\varrho_{cr}=(1+\delta_{cr})\bar{\varrho}_{m}$ is, given an RMS amplitute at the time of collapse z, $\delta\varrho_{RMS}=D^{(1)}(z)\sigma(M)\bar{\varrho}_{m}$. Now following Birkhoff’s theorem, the critical density for collapse, $\varrho_{cr}$, will be unaffected in the presence a density perturbation $\Delta$, since a collapsing overdensity will not depend on the external spacetime. Combining this fact with equation (113), relation (115) gives (Schmidt et al., 2013)

[TABLE]

The conditional halo mass function $\bar{n}_{h}(M,\Delta)$ is now given by the universal prescription (47), but with the peak significance $\nu^{\prime}_{c}(M)$ in (116). Combining this fact with the rigorous definition of the bias (52), gives the known PBS biases (53) for GR.

Let us now turn to the MG case. Following the discussion in Section III.1.2, it is now clear how the density threshold will not be scale and redshift independent anymore, but will depend on the comoving halo mass $M$, as well as the environmental density $\delta_{env}$, that is, a function $\delta_{cr}(M,\delta_{env})$. As a consequence, the peak significance in MG will now become:

[TABLE]

In the presence of a long-wavelength density perturbation $\Delta$, $\bar{\varrho}_{m}$ will again change according to (113), but also $\delta_{env}$ will now change as dictated by (114), and so will $\delta_{cr}$ that depends on it, which will become $\delta^{\prime}_{cr}(M,\delta^{\prime}_{env})=\delta_{cr}(M,\delta_{env}+\Delta)$. From the equivalent of (116), the MG peak significance will now become:

[TABLE]

As in the GR case, the conditional halo mass function $\bar{n}_{h}(M,\Delta)$ for MG is now given by the universal prescription (47), but with the peak significance $\nu^{\prime}_{cMG}(M)$ given by (118). Calculating, now, the first and second order derivatives that we need in (52), we will have, starting with the linear order:

[TABLE]

where we used the chain rule a few times and also the fact that $\nu^{\prime}_{cMG}(M)=\nu_{cMG}(M)$ at $\Delta=0$. The definition (52), combined with the result (LABEL:deriv1) and the universal prescription (47), gives:

[TABLE]

Similarly, the second order derivative will be:

[TABLE]

which will give the expression for the second order bias factor:

[TABLE]

Equations (120) and (122) give the Lagrangian PBS biases of first and second order in MG, for any universal mass function $f\left[\nu_{cMG}\right]$. Applying these to the particular ST form (50), we arrive at the relationships (LABEL:pbsbiasMG) and (LABEL:pbsbiasMGavg) in the main text. The derivatives of the form $\frac{d\delta_{cr}(M,\delta_{env})}{\delta_{env}}$ can be easily calculated numerically, as soon as the function $\delta_{cr}(M,\delta_{env})$ is known from integrating equation (68). In MG models that still possess a scale and redshift $\it{independent}$ barrier, like the $n$DGP model (and possibly other models in the Vainshtein family), these derivatives will vanish and we recover the standard GR expressions for the PBS biases (53), but with a different $\delta_{cr}$.