© ESO 2019

1. Introduction

Two decades have passed since the discovery of the present-day acceleration of the Universe (Riess et al. 1998; Perlmutter et al. 1999). Its physical origin, however, is still a mystery. The simplest choice is that this acceleration is the result of the vacuum energy, in the form of a cosmological constant, in Einstein’s equations. This model is in agreement with current observations, but is plagued by the fine-tuning problem. Another simple option is that dark energy is a dynamical field, as in quintessence models. This field may also interact with the dark matter sector, giving rise to interacting dark energy models (Amendola & Tsujikawa 2010).

Instead of modifying the particle content of the Universe, an alternative solution is extending Einstein gravity in such a way that the current acceleration of the Universe is accounted for. However, precise measurements on Earth, in the solar system, in binary pulsars, and of gravitational waves (LIGO Collaboration & Virgo Collaboration 2016; Bertotti et al. 2003; Will 2006; Williams et al. 2004) strongly constrain deviations); strongly constrain deviations from general relativity (GR) at these scales (Will 2014). This is a challenge to theories of gravity beyond Einstein GR, which claim to be an explanation of the present-day acceleration (Clifton et al. 2012; Joyce et al. 2015; Bull et al. 2016).

We note, however, that all these tests probe extensions to GR in astrophysical systems that reside in dense galactic environments. Conditions are therefore far from the cosmological background density, curvature, and even the gravitational potential differs from that of the background by several orders of magnitude (Baker et al. 2015).

A way of evading the high-density environmental constraints while still allowing deviations from GR on cosmological scales is to hide modifications to Einstein’s gravity in these environments. The idea is that the geometry is modified by introducing a new degree of freedom that drives the acceleration of the Universe. This degree of freedom, in the simplest case a scalar field, would be suppressed in high-density or curvature environments, however.

Several screening mechanisms have been proposed in the literature, and they can be classified in several ways (Joyce et al. 2015; Brax 2013; Brax & Davis 2015; Koyama 2016). We here investigate screening mechanisms that result from one of the following three properties: i) Weak coupling, in which the coupling to matter fields becomes weak in regions of high density, hence suppressing the fifth force. At large scales, where the density is low, the coupling becomes large and the fifth force acts strongly. Examples of theories of this type are symmetron models (Hinterbichler & Khoury 2010; Pietroni 2005; Olive & Pospelov 2008) and varying dilaton (Damour & Polyakov 1994; Brax et al. 2011). ii) Large mass, when the mass of the field becomes large in regions of high density, suppressing the range of the fifth force. At low densities, the scalar field mass is small and mediates a long range fifth force. Examples of this type are chameleons (Khoury & Weltman 2004a,b). Finally, iii) large inertia, when the scalar field kinetic function depends on the environment, being large in high dense regions (and then the Yukawa potential is suppressed) and small in low dense regions. Examples of the former are K-mouflage models (Babichev et al. 2009, 2011), and examples of the latter are Vainshtein models (Vainshtein 1972).

Screening mechanisms are a relatively generic prediction of viable modified gravity (MG) theories (Brax et al. 2012). Detecting them would therefore be a signature of physics beyond GR. Observational tests of screening focus on the transition between the fully screened and unscreened regimes, where deviations from GR are expected to be most pronounced. For viable cosmological models, this transition occurs in regions where the matter density and gravitational potential are nonlinear and start to approach their linear or background values: this can be observed in the outskirts of dark matter halos and their properties (e.g. Shirata et al. 2005; Davis et al. 2012; Oyaizu et al. 2008; Schmidt et al. 2009; Barreira et al. 2013; Wyman et al. 2013; Zhao 2014; Clifton et al. 2005; Lombriser et al. 2015, 2012; Martino & Sheth 2009; Clampitt et al. 2013; Llinares & Mota 2013; Hellwing et al. 2013; Gronke et al. 2014, 2015; Stark et al. 2016).

Screening mechanisms are in fact a nonlinear effect. Therefore, predictions of their signatures in astrophysical systems are computed by integrating the fully nonlinear equations of motion for the extra-degree of freedom of gravity and also using N-body simulations in order to simulate nonlinear structure formation. These computations are extremely time-consuming, however. It is therefore not viable to use them for a parameter estimation or even to place constraints on the parameter space of such theories using nonlinear codes of structure formation. It is then imperative to find other methods to describe and probe screening mechanisms in faster and accurate ways. Although screenings are more commonly studied in the highly nonlinear regime, they leave imprints on quasi-linear scales that can be captured by cosmological perturbation theory (PT) and should be considered in the low-density regions of large-scale structure formation; see for example (Koyama et al. 2009).

On the other hand, PT has experienced many developments in recent years (Matsubara 2008a; Baumann et al. 2012; Carlson et al. 2013), in part because it can be useful to analytically understand different effects in the power spectrum and correlation function for the dark matter clustering. These effects can be confirmed or refuted, and further explored with simulations to ultimately understand the outcomes of present and future galaxy surveys, such as eBOSS (Zhao 2016), DESI (Aghamousa et al. 2016), EUCLID (Amendola et al. 2013), and LSST (LSST Dark Energy Science Collaboration 2012). Two approaches have been used to study PT: the Eulerian standard PT (SPT) and Lagrangian PT (LPT), which both have advantages and drawbacks, but they are complementary in the end (Tassev 2014). The nonlinear PT for MG was developed initially in (Koyama et al. 2009), and has been further studied in several other works (Taruya et al. 2014a,b; Brax & Valageas 2013; Bellini & Zumalacarregui 2015; Taruya 2016; Bose & Koyama 2016, 2017; Barrow & Mota 2003; Akrami et al. 2013; Fasiello & Vlah 2017; Aviles & Cervantes-Cota 2017; Hirano et al. 2018; Bose et al. 2018; Bose & Taruya 2018; Aviles et al. 2018). The LPT for dark matter fluctuations in MG was developed in Aviles & Cervantes-Cota (2017), and further studies for biased tracers in Aviles et al. (2018). The PT for MG has the advantage that it allows us to understand the role that these physical parameters play in the screening features of dark matter statistics. We here study some of these effects through screening mechanisms by examining them at second- and third-order perturbation levels using PT for some MG models. To this end, we build on the LPT formalism developed in Aviles & Cervantes-Cota (2017), which was initially posited for MG theories in the Jordan frame, in order to apply it to theories in the Einstein frame. Because of the direct coupling of the scalar field and the dark matter in the Klein–Gordon equation, the equations that govern the screening can differ substantially from those in Jordan-frame MG theories. In general, screening effects depend on the type of nonlinearities that are introduced in the Lagrangian density. We present a detailed analysis of screening features and identify the theoretical roots of their origin. Our results show that screenings possess peculiar features that depend on the scalar field effective mass and couplings, and that may in particular cases cause anti-screening effects in the power spectrum, such as in the symmetron model. We perform this analysis by separating the growth functions into screening and non-screened parts. We note, however, that we do not compare the perturbative approach with a fully nonlinear simulation. We refer to (Koyama et al. 2009), for instance, for investigations like this at the level of the power spectrum.

This work is organized as follows: in Sect. 2 we set up the formalism for perturbation theory in both the Einstein and Jordan frames; in Sect. 3 we apply these methods to the specific gravity models investigated here; in Sect. 4 we show the matter power spectra and in Sect. 5 the screening growth function analysis. We conclude in Sect. 6 with a discussion of our results. Some formulae are displayed in Appendix A.

2. Perturbation theory in the Einstein frame

In this section we are interested in MG theories defined in the Einstein frame with action

$$\begin{array}{c}\hfill S={\displaystyle \int {\mathrm{d}}^{4}x\sqrt{-g}(\frac{M_{\mathrm{Pl}}^2}{2}R-{(\mathrm{\nabla }\phi )}^2-V(\phi ))+S_{\mathrm{m}}[{\stackrel{\sim }{g}}_{\mu \nu }]},\end{array}$$(1)

with the conformal metric

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_{\mu \nu }={(C(\phi ))}^2{g}_{\mu \nu },\end{array}$$(2)

where C(φ) is a conformal factor (in the literature it is more common to find A(φ); we use C instead because A is used below to characterize the strength of the gravitational force). By taking variations of action Eq. (1) with respect to the scalar field, we obtain the Klein–Gordon equation

$$\begin{array}{c}\hfill \square \phi =\frac{\mathrm{d}{V}_{\mathrm{eff}}}{\mathrm{d}\phi },V_{\mathrm{eff}}=V(\phi )-T(C(\phi )-1),\end{array}$$(3)

where T is the trace of the energy momentum tensor of matter. In PT we split the scalar field into background $\overline{\phi }$ and perturbed δφ parts,

$$\begin{array}{c}\hfill \phi (\mathbf{x},t)=\overline{\phi }(t)+\delta \phi (\mathbf{x},t).\end{array}$$(4)

Hereafter, a bar over a dynamical quantity means that we refer to its homogeneous and isotropic background value; we also assume a dark matter perfect fluid with T = − ρ. In the following we adopt the quasi-static limit for the perturbed piece, which relies on neglecting temporal derivatives in the Klein–Gordon equation, thus Eq. (3) becomes

$$\begin{array}{cc}\hfill \frac{1}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}^2\delta \phi & ={\displaystyle \sum _{n=1}^{\infty }\frac{1}{n!}[V^{(n+1)}(\overline{\phi })+\overline{\rho }{C}^{(n+1)}(\overline{\phi })]{(\delta \phi )}^n}\hfill \\ \hfill & +\overline{\rho }\delta {\displaystyle \sum _{n=0}^{\infty }\frac{1}{n!}{C}^{(n+1)}(\overline{\phi }){(\delta \phi )}^n}\hfill \\ \hfill & =\frac{\beta }{M_{\mathrm{Pl}}}\overline{\rho }\delta +m^2(\overline{\phi })\delta \phi +{\displaystyle \sum _{n=2}^{\infty }\frac{1}{n!}{M}_{\mathrm{Pl}}^{1-n}{\kappa }_{n+1}{(\delta \phi )}^n}\hfill \\ \hfill & +\overline{\rho }\delta {\displaystyle \sum _{n=1}^{\infty }}\frac{1}{n!}{M}_{\mathrm{Pl}}^{-1-n}{\beta }_{n+1}{(\delta \phi )}^n,\hfill \end{array}$$(5)

where we have subtracted the background evolution. Here $V^{(n)}(\overline{\phi })$ and $C^{(n)}(\overline{\phi })$ denote the n-ésime derivative of V and C functions evaluated at background values. In the above equation we also introduced the matter overdensity δ, defined through $\rho (\mathbf{x},t)=\overline{\rho }(t)(1+\delta (\mathbf{x},t))$. We also introduce following (Brax & Valageas 2013)

$$\begin{array}{c}\hfill {\kappa }_n(\overline{\phi })=M_{\mathrm{Pl}}^{n-2}(V^{(n)}(\overline{\phi })+\overline{\rho }{C}^{(n)}(\overline{\phi })),m^2={\kappa }_2.\end{array}$$(6)

$$\begin{array}{c}\hfill {\beta }_n(\overline{\phi })=M_{\mathrm{Pl}}^n{C}^{(n)}(\overline{\phi }),\beta ={\beta }_1.\end{array}$$(7)

We work in Lagrangian space, where the position x of a dark matter particle or fluid element with initial Lagrangian coordinate q is given by

$$\begin{array}{c}\hfill \mathbf{x}(\mathbf{q},t)=\mathbf{q}+\mathbf{\Psi }(\mathbf{q},t),\end{array}$$(8)

where Ψ is the Lagrangian displacement vector field. We further assume that Ψ is longitudinal and that it is a Gaussian-distributed variable at linear order. Dark matter particles follow geodesics of the conformal metric g̃,

$$\begin{array}{c}\hfill \ddot{\mathbf{\Psi }}+2H\dot{\mathbf{\Psi }}=-\frac{1}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}{\psi }_{\mathrm{N}}-\frac{1}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}logC(\phi ),\end{array}$$(9)

where ψ_N denotes the Newtonian potential, which obeys the Poisson equation¹

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathbf{x}}^2{\psi }_{\mathrm{N}}=4\pi G{a}^2\overline{\rho }\delta (\mathbf{x}).\end{array}$$(10)

A noticeable difference between MG theories defined in the Einstein and Jordan frames is that in the latter case, the new scalar degree-of-freedom is a source of the Poisson equation instead of the geodesic equation.

Equation (8) can be regarded as a coordinate transformation between Lagrangian and Eulerian coordinates, with Jacobian matrix J_ij = (∂xⁱ/∂q^j)=δ_ij + Ψ_i, j and Jacobian determinant J = det(J_ij). From mass conservation, a relation between matter overdensities and Lagrangian displacement can be obtained (Bouchet et al. 1995),

$$\begin{array}{c}\hfill \delta (\mathbf{x},t)=\frac{1-J(\mathbf{q},t)}{J(\mathbf{q},t)}·\end{array}$$(11)

The set of Eqs. (3), (9), and (10) is treated perturbatively in order to solve for the displacement field. Instead of working with the field δφ, however, we choose to define the rescaled field

$$\begin{array}{c}\hfill \chi (\mathbf{q},t)\equiv -\frac{2\beta }{C}\frac{\delta \phi (\mathbf{q},t)}{M_{\mathrm{Pl}}},\end{array}$$(12)

hereafter, we denote $C\equiv C(\overline{\phi })$ unless otherwise explicitly stated. We note that Eq. (12) is not a conformal transformation since the functions β and C are not free functions, but are evaluated at the background. By taking the divergence of Eq. (9), we have

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathbf{x}}·(\ddot{\mathbf{\Psi }}+2H\dot{\mathbf{\Psi }})=-4\pi G\overline{\rho }\delta (\mathbf{x})+\frac{1}{2a^2}{\mathrm{\nabla }}_{\mathbf{x}}^2\chi ,\end{array}$$(13)

which has the structure of an MG theory in the Jordan frame, and the PT formalism developed in Aviles & Cervantes-Cota (2017) applies directly. In Lagrangian space, we work in q-Fourier space, in which the transformation is taken with respect to q-coordinates, thus when transforming gradients with respect to x-coordinates, frame-lagging terms are introduced,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathbf{x}}^2\chi ={\mathrm{\nabla }}^2\chi +({\mathrm{\nabla }}_{\mathbf{x}}^2\chi -{\mathrm{\nabla }}^2\chi ).\end{array}$$(14)

These frame-lagging terms are necessary to obtain the correct limit of the theory at large scales, particularly for the theories in which the associated fifth force is short-ranged and the ΛCDM limit at large scales should be recovered. Since at linear order spatial derivatives with respect to Eulerian and Lagrangian coordinates coincide, the frame-lagging contributes nonlinearly to the theory. For comparison with other works, it is worthwhile to introduce the quantities

$$\begin{array}{cc}& A(k,a)=A_0(1+\frac{2{\beta }^2}{C}\frac{k^2/a^2}{k^2/a^2+m^2}),\hfill \end{array}$$(15)

$$\begin{array}{cc}& A_0(a)=4\pi G\overline{\rho },\hfill \end{array}$$(16)

$$\begin{array}{cc}& \mathrm{\Pi }(k,a)=\frac{C}{6a^2{\beta }^2}(k^2+m^2{a}^2),\hfill \end{array}$$(17)

$$\begin{array}{cc}& 3+2{\omega }_{\mathrm{BD}}(a)=\frac{C}{2{\beta }^2},\hfill \end{array}$$(18)

$$\begin{array}{cc}& M_1(a)=\frac{C}{2{\beta }^2}{m}^2.\hfill \end{array}$$(19)

These equations can be used as a translation table for different PT works in MG. Particularly in Aviles & Cervantes-Cota (2017), the M₁ and ω_BD functions are used extensively instead of β and m.

We return to the Klein–Gordon equation, which in q-Fourier space for a field χ is

$$\begin{array}{cc}\hfill -\frac{k^2}{2a^2}\chi (\mathbf{k})=& -(A(k)-A_0)\stackrel{\sim }{\delta }(\mathbf{k})+\frac{k^2/a^2}{6\mathrm{\Pi }(k)}\delta \mathcal{I}(\chi )-\hfill \\ \hfill & -\frac{C}{2{\beta }^2}\frac{k^2/a^2}{3\mathrm{\Pi }(k)}\frac{1}{2a^2}[({\mathrm{\nabla }}_{\mathbf{x}}^2\chi -{\mathrm{\nabla }}^2\chi )](\mathbf{k}),\hfill \end{array}$$(20)

where [(⋯)](k) means the q-Fourier transform of (⋯)(q), and we also note that δ̃(k) ≡ ∫ d³ qe^−ik·qδ(x). To avoid confusion with the q-Fourier transform of δ(q) or the x-Fourier transform of δ(x), we write a tilde over this overdensity. Equation (20) is derived directly using Eq. (12), the second term on the right-hand side (RHS) encodes all nonlinear terms of Eq. (5), while the last term arises when derivatives are transformed from Eulerian to Lagrangian coordinates. Specifically, the contribution from screenings is given by

$$\begin{array}{cc}\hfill \delta \mathcal{I}(\mathbf{k})& ={\displaystyle \sum _{n=2}^{\infty }\frac{{(-1)}^{n+1}}{2^{n}n!}\frac{C^n{\kappa }_{n+1}}{{\beta }^{n+1}}[{\chi }^n](\mathbf{k})}\hfill \\ \hfill & +{\displaystyle \sum _{n=1}^{\infty }}\frac{{(-1)}^{n+1}}{2^{n}n!}\frac{2A_0{C}^n{\beta }_{n+1}}{{\beta }^{n+1}}[{\chi }^n\delta ](\mathbf{k}).\hfill \end{array}$$(21)

We formally expand quantities as²

$$\begin{array}{cc}\hfill \delta \mathcal{I}(\mathbf{k})=& \frac{1}{2}{\displaystyle \underset{{\mathbf{k}}_{12}=\mathbf{k}}{\int }{M}_2({\mathbf{k}}_1,{\mathbf{k}}_2)\chi ({\mathbf{k}}_1)\chi ({\mathbf{k}}_2)}+\hfill \\ \hfill & +\frac{1}{6}{\displaystyle \underset{{\mathbf{k}}_{123}=\mathbf{k}}{\int }{M}_3({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)\chi ({\mathbf{k}}_1)\chi ({\mathbf{k}}_2)\chi ({\mathbf{k}}_3)+\cdots },\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill -\frac{1}{2a^2}& [({\mathrm{\nabla }}_{\mathbf{x}}^2\chi -{\mathrm{\nabla }}^2\chi )](\mathbf{k})=\frac{1}{2}{\displaystyle \underset{{\mathbf{k}}_{12}=\mathbf{k}}{\int }{\mathcal{K}}_{\mathrm{FL}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2){\delta }_L({\mathbf{k}}_1){\delta }_L({\mathbf{k}}_2)}\hfill \\ \hfill & +\frac{1}{6}{\displaystyle \underset{{\mathbf{k}}_{123}=\mathbf{k}}{\int }{\mathcal{K}}_{\mathrm{FL}}^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3){\delta }_L({\mathbf{k}}_1){\delta }_L({\mathbf{k}}_2){\delta }_L({\mathbf{k}}_3)}+\cdots .\hfill \end{array}$$(24)

Expressions for the frame-lagging kernels are derived in Aviles & Cervantes-Cota (2017). For example, to second order, we have

$$\begin{array}{cc}\hfill {\mathcal{K}}_{\mathrm{FL}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)& =2x^2(A(k_1)+A(k_2)-2A_0)+x\frac{k_2}{k_1}(A(k_1)-A_0)\hfill \\ \hfill & +x\frac{k_1}{k_2}(A(k_2)-A_0),\hfill \end{array}$$(25)

with x = k̂₁ · k̂₂.

There is a crucial difference between theories in the Jordan and Einstein frames. In the latter there is a direct coupling χδ between the scalar field and the matter density, as can be seen from Eq. (3) or from Eq. (21). This leads us to expand the overdensity in terms of the scalar field as

$$\begin{array}{c}\hfill \stackrel{\sim }{\delta }(\mathbf{k})={\displaystyle \sum _{n=1}^{\infty }\frac{1}{n!}\underset{{\mathbf{k}}_{1\cdots n}=\mathbf{k}}{\int }{\mathcal{K}}_{\chi \delta }^{(n)}({\mathbf{k}}_1,\cdots ,{\mathbf{k}}_n)\chi ({\mathbf{k}}_1)\cdots \chi ({\mathbf{k}}_n)}.\end{array}$$(26)

After some iterative manipulations of Eqs. (5), (20), (21), (23), and (26), we arrive at

$$\begin{array}{c}\hfill M_2({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{2C{\beta }_2{A}_0}{{\beta }^2}{\mathcal{K}}_{\chi \delta }^{(1)}({\mathbf{k}}_1)-\frac{C^2{\kappa }_3}{4{\beta }^3},\end{array}$$(27)

$$\begin{array}{cc}\hfill M_3({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)& =\frac{3C{\beta }_2{A}_0}{{\beta }^2}{\mathcal{K}}_{\chi \delta }^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)\hfill \\ \hfill & -\frac{3C^2{\beta }_3{A}_0}{2{\beta }^3}{\mathcal{K}}_{\chi \delta }^{(1)}({\mathbf{k}}_1)+\frac{C^3{\kappa }_4}{8{\beta }^4},\hfill \end{array}$$(28)

with the kernels given by

$$\begin{array}{cc}& {\mathcal{K}}_{\chi \delta }^{(1)}({\mathbf{k}}_1)=\frac{3}{2A_0}\mathrm{\Pi }(k_1),\hfill \end{array}$$(29)

$$\begin{array}{cc}& {\mathcal{K}}_{\chi \delta }^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{1}{2A_0}{M}_2^{\mathrm{FL}}({\mathbf{k}}_1,{\mathbf{k}}_2),\hfill \end{array}$$(30)

where the frame-lagged M₂ function is defined in Eq. (A.3). It is convenient to symmetrize these M functions over their arguments, as we do below.

In theories of gravity defined in the Jordan frame, M₂ and M₃ are k-dependent if non-canonical kinetic terms or higher derivatives of the scalar field are present in the Lagrangian; a known case with such scale dependencies is the Dvali Gabadadze Porrati (DGP) braneworld model, whereas no scale dependencies in the Ms is found in f(R) Hu-Sawicki gravity; see Koyama et al. (2009). For theories in the Einstein frame, the k-dependence arises through the couplings χδ in the Klein–Gordon equation, even if no derivatives other than the standard kinetic term appear in their defining action.

It is worth mentioning that functions M₂ and M₃ encode the physics of particular theories, and they determine the screening properties as well; these are the coefficients of Taylor-expanding the nonlinearities of the Klein–Gordon in Fourier space. As we note, these functions can be positive or negative, which will be responsible for the screening properties of a model. Moreover, if β₂ and β₃ are zero, as happens for theories with a conformal factor that is linear in the scalar field, both M₂ and M₃ become scale independent. This is, for example, the case of the first proposed chameleon model (Khoury & Weltman 2004b).

We define the linear differential operator (Matsubara 2015)

$$\begin{array}{c}\hfill \widehat{\mathcal{T}}=\frac{{\partial }^2}{\partial t^2}+2H\frac{\partial }{\partial t},\end{array}$$(31)

and the equation of motion for the displacement field divergence (Eq. (13)) becomes (Aviles & Cervantes-Cota 2017)³

$$\begin{array}{cc}\hfill [{(J^{-1})}_{\mathit{ji}}\widehat{\mathcal{T}}{\mathrm{\Psi }}_{i,j}](\mathbf{k})& =-A(k)\stackrel{\sim }{\delta }(\mathbf{k})+\frac{k^2/a^2}{6\mathrm{\Pi }(\mathbf{k})}\delta I(\mathbf{k})\hfill \\ \hfill & +\frac{M_1}{3\mathrm{\Pi }(k)}\frac{1}{2a^2}[({\mathrm{\nabla }}_{\mathbf{x}}^2\chi -{\mathrm{\nabla }}^2\chi )](\mathbf{k}).\hfill \end{array}$$(32)

We perturb the displacement field as Ψ = λΨ⁽¹⁾ + λ²Ψ⁽²⁾ + λ³Ψ⁽³⁾ + 𝒪(λ⁴), and solve the above equation order by order. Stopping at third order allows us to calculate the first corrections to the linear power spectrum. Hereafter, we absorb the control parameter λ into the definition of Ψ. To first order, Eq. (32) yields

$$\begin{array}{c}\hfill (\widehat{\mathcal{T}}-A(k))(i{k}_i{\mathrm{\Psi }}_i^{(1)}(\mathbf{k}))=0.\end{array}$$(33)

This equation has the same form as the linear equation for the matter overdensity δ(k, t). Therefore, we obtain

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_i^{(1)}(\mathbf{k},t)=i\frac{k_i}{k^2}{D}_{+}(\mathbf{k},t){\stackrel{\sim }{\delta }}_L(\mathbf{k},t_0),\end{array}$$(34)

with D₊(k) the fastest growing solution of equation $(\widehat{\mathcal{T}}-A(k))D(k)=0$ normalized to unity as D₊(k = 0, t₀) = 1. This normalization is useful for theories that reduce to ΛCDM at very large scales, which is the case when the fifth-force range is finite. The initial condition ${\stackrel{\sim }{\delta }}_L(\mathbf{k},t_0)$ is fixed by noting that by linearizing the RHS of Eq. (11), we have ${\delta }^{(1)}(\mathbf{x})=-{\mathrm{\Psi }}_{i,i}^{(1)}(\mathbf{q})$. Because we study linear fields, we can safely drop the tilde over the overdensity in Eq. (34). To second order, Eq. (32) leads to the solution

$$\begin{array}{c}\hfill {\mathrm{\Psi }}^{i(2)}(\mathbf{k})=\frac{i{k}^i}{2k^2}{\displaystyle \underset{{\mathbf{k}}_{12}=\mathbf{k}}{\int }\frac{3}{7}({\overline{D}}_{\mathrm{NS}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)-{\overline{D}}_{\mathrm{S}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)){\delta }_1{\delta }_2},\end{array}$$(35)

where we denote δ_1, 2 ≡ δ_L(k_1, 2). Momentum conservation implies k = k₁ + k₂, as is explicit in the Dirac delta function, cf. Eq. (22). We split the second-order growth into non-screened (NS) and screening (S) parts. These growth functions D⁽²⁾ are solutions, with the appropriate initial conditions, to the equations

$$\begin{array}{cc}& (\widehat{\mathcal{T}}-A(k))D_{\mathrm{NS}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)=(A(k)-(A(k_1)+A(k_2)-A(k))\frac{{({\mathbf{k}}_1·{\mathbf{k}}_2)}^2}{k_1^2{k}_2^2}\hfill \\ \hfill & +(A(k)-A(k_1))\frac{{\mathbf{k}}_1·{\mathbf{k}}_2}{k_2^2}+(A(k)-A(k_2))\frac{{\mathbf{k}}_1·{\mathbf{k}}_2}{k_1^2})D_{+}(k_1)D_{+}(k_2)\hfill \end{array}$$(36)

$$\begin{array}{cc}& (\widehat{\mathcal{T}}-A(k))D_{\mathrm{S}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)={\left(\frac{2A_0}{3}\right)}^2\frac{k^2}{a^2}\frac{M_2({\mathbf{k}}_1,{\mathbf{k}}_2)D_{+}(k_1)D_{+}(k_2)}{6\mathrm{\Pi }(k)\mathrm{\Pi }(k_1)\mathrm{\Pi }(k_2)}\hfill \\ \hfill & \equiv {\mathcal{S}}_{\mathrm{S}}^{(2)},\hfill \end{array}$$(37)

and the normalized growth functions are defined as

$$\begin{array}{cc}\hfill {\overline{D}}_{\mathrm{S},\mathrm{NS}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2,t)& =\frac{7}{3}\frac{D_{\mathrm{S},\mathrm{NS}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2,t)}{D_{+}(k_1)D_{+}(k_2)}·\hfill \end{array}$$(38)

In an EdS universe, we obtain the well-known result ${\overline{D}}_{\mathrm{NS}}^{(2)}=1-{({\widehat{\mathbf{k}}}_1·{\widehat{\mathbf{k}}}_2)}^2$, while in ΛCDM, we obtain the same result multiplied by a function that varies slowly with time, such that today, ${\overline{D}}_{\mathrm{NS}}^{(2)\mathrm{\Lambda }\mathrm{CDM}}\simeq 1.01{\overline{D}}_{\mathrm{NS}}^{(2)\mathrm{EdS}}$. The screening second-order growth, ${\overline{D}}_{\mathrm{S}}^{(2)}$, is zero in both EdS and ΛCDM models.

The function ${\overline{D}}_{\mathrm{S}}^{(2)}$ is important for our discussion. It encodes the nonlinearities of the “potential” of the scalar field, and it yields the second-order screening effects that drive the theory to GR at small scales. The total second-order growth function, as can be read from Eq. (35), is given by $D^{(2)}=D_{\text{NS}}^{(2)}-D_{\text{S}}^{(2)}$, such that negative values of $D_{\text{S}}^{(2)}$ enhance the growth of perturbations (anti-screening effects), while positive values of it yield the standard suppression of the fifth force.

Analogously, each higher perturbative order carries its own screening and is efficient over a certain k interval. The third-order Lagrangian displacement field can be computed to give

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_i^{(3)}(\mathbf{k})=\frac{i{k}_i}{6k^2}{\displaystyle \underset{{\mathbf{k}}_{123}=\mathbf{k}}{\int }{\overline{D}}^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2{\mathbf{k}}_3){\delta }_1{\delta }_2{\delta }_3},\end{array}$$(39)

with normalized growth

$$\begin{array}{c}\hfill {\overline{D}}^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2{\mathbf{k}}_3)=\frac{D^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2{\mathbf{k}}_3)}{D_{+}(k_1)D_{+}(k_2)D_{+}(k_3)}·\end{array}$$(40)

The complete expression for the third-order growth function D⁽³⁾, equivalent to Eq. (35), is large and can be found in Aviles & Cervantes-Cota (2017), although in Eq. (69) below we show this function for a particular configuration of wave vectors. In Sect. 5 we split the growth into non-screened and screening parts. Unlike the second-order case, at third order, the decomposition cannot be performed directly through the linear differential equations that govern the growth. In this case, ${\overline{D}}_{\mathrm{NS}}^{(3)}$ is therefore obtained by setting M₂ and M₃ equal to zero,

$$\begin{array}{c}\hfill {\overline{D}}_{\mathrm{NS}}^{(3)}\equiv {\overline{D}}^{(3)}{|}_{M_2=M_3=0},\end{array}$$(41)

while ${\overline{D}}_{\mathrm{S}}^{(3)}$ is obtained through the relation

$$\begin{array}{c}\hfill {\overline{D}}^{(3)}={\overline{D}}_{\mathrm{NS}}^{(3)}-{\overline{D}}_{\mathrm{S}}^{(3)}.\end{array}$$(42)

In this way, third-order screenings are realized by having ${\overline{D}}_{\mathrm{S}}^{(3)}>0$, while anti-screening is realized by ${\overline{D}}_{\mathrm{S}}^{(3)}<0$.

3. Modified gravity theories with different screening mechanisms

As shown above, expressions M₂(k₁, k₂) and M₃(k₁, k₂) depend upon the explicit form of the Klein–Gordon equation, which in turn depends on the frame that is posed. The formalism developed here applies to the Einstein frame in which the coupling function C(φ) is nontrivial, but applies also to the Jordan frame by setting $C(\overline{\phi })\equiv 1$, as explicitly done in Aviles & Cervantes-Cota (2017).

In the following we treat examples of models with different screening properties. We start with symmetron models that are posed in the Einstein frame and follow with f(R) Hu-Sawicki and DGP models posed in the Jordan frame.

3.1. Symmetron model

The symmetron model can be introduced with the action of Eq. (1) with a self-interacting potential

$$\begin{array}{cc}\hfill V(\phi )& =V_0-\frac{1}{2}{\mu }^2{\phi }^2+\frac{1}{4}\lambda {\phi }^4,\hfill \end{array}$$(43)

and the conformal factor

$$\begin{array}{cc}\hfill C(\phi )& =1+\frac{1}{4}{\left(\frac{\phi }{M}\right)}^2.\hfill \end{array}$$(44)

Assuming that the background piece of the scalar field always sits in the minimum of the effective potential $V_{\mathrm{eff}}=V+\overline{\rho }C$, we obtain

$$\begin{array}{c}\hfill \overline{\phi }={\overline{\phi }}_0\sqrt{1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3},\end{array}$$(45)

where a_ssb is the scale factor at which the ℤ₂ symmetry is broken. The scalar field effective mass and the strength β of the fifth force are

$$\begin{array}{cc}\hfill m(a)& =m_0\sqrt{1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3},\beta (a)={\beta }_0\sqrt{1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3}.\hfill \end{array}$$(46)

These functions are commonly generalized to

$$\begin{array}{cc}\hfill m(a)& =m_0{[1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3]}^{\widehat{m}},\beta (a)={\beta }_0{[1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3]}^{\widehat{n}}.\hfill \end{array}$$(47)

In this way, a symmetron model can be characterized by the set of parameters $(a_{\mathrm{ssb}},m_0,{\beta }_0,\widehat{n},\widehat{m})$, although other equivalent parameters are also used in the literature; for example, in order to contain the parameters μ, λ, and M, instead. Since variations of fermion masses cannot vary too much over the Universe lifetime, we simply set $C(\overline{\phi })=1$.

The function M₁ plays an important role in the upcoming discussion. It is given by

$$\begin{array}{c}\hfill M_1(a)=\frac{Cm{(a)}^2}{2\beta {(a)}^2}=\frac{C{m}_0^2}{2{\beta }_0^2}{[1-{\left(\frac{a_{\mathrm{ssb}}}{a}\right)}^3]}^{2(\widehat{m}-\widehat{n})}.\end{array}$$(48)

We soon realize that if $\widehat{m}=\widehat{n}$, M₁ becomes a constant. Expressions for M₂ and M₃, given by Eqs. (27) and (28), depend on the conformal coupling, β_n, κ_n, and ${\mathcal{K}}_{\chi \delta }^{(n)}$. In the present case, these formulae are cumbersome and we do not show them here, but we solved them numerically when we integrated the differential equations for the growth functions and to construct power spectra. M₂ and M₃ are indeed important to determine the fate of nonlinearities. Cosmological screenings are encoded in these functions, and they therefore serve to distinguish among different screening types. In the present case, M₂ and M₃ are negative for a certain cosmological epoch and specific wave numbers that are reflected in an anti-screening effect in the power spectra shown in the next section.

3.2. f(R) Hu-Sawicki model: chameleon mechanism

Here we consider a Lagrangian density given by $\mathcal{L}=\frac{1}{2}{M_{\text{Pl}}}^2\sqrt{-\mathit{g}}(R+f(R))$, in contrast to Eq. (1). As is known, f(R) models can be brought to a Jordan-frame description, and our perturbation formalism (Aviles & Cervantes-Cota 2017) can then be applied. We analyze the Hu-Sawicki model with parameter n = 1. We can define the scalar degree of freedom to be χ = δf_R, with f_R = df/dR, for which we find a Klein–Gordon equation that in the quasi-static limit is

$$\begin{array}{c}\hfill \frac{3}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}^2\chi =-2A_0\delta +\delta R,\end{array}$$(49)

where $\delta R=R(f_R)-R({\overline{f}}_R)=M_1\chi +\frac{1}{2}{M}_2{\chi }^2+\frac{1}{6}{M}_3{\chi }^3+\cdots$. By developing Eq. (49) in terms of the scalar field χ, and setting C = 1, we arrive at Eq. (20) with 2β² = 1/3. The M functions are obtained by using $R(f_R)\simeq \overline{R}{(f_{R0}/f_R)}^{1/2}$ and $\overline{R}=3H_0^2({\mathrm{\Omega }}_{m0}{a}^{-3}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})$:

$$\begin{array}{c}\hfill M_1(a)=\frac{3}{2}\frac{H_0^2}{|f_{R0}|}\frac{{({\mathrm{\Omega }}_{m0}{a}^{-3}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^3}{{({\mathrm{\Omega }}_{m0}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^2},\end{array}$$(50)

$$\begin{array}{c}\hfill M_2(a)=\frac{9}{4}\frac{H_0^2}{|f_{R0}{|}^2}\frac{{({\mathrm{\Omega }}_{m0}{a}^{-3}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^5}{{({\mathrm{\Omega }}_{m0}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^4},\end{array}$$(51)

$$\begin{array}{c}\hfill M_3(a)=\frac{45}{8}\frac{H_0^2}{|f_{R0}{|}^3}\frac{{({\mathrm{\Omega }}_{m0}{a}^{-3}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^7}{{({\mathrm{\Omega }}_{m0}+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^6},\end{array}$$(52)

while the scalar field mass is given by $m=\sqrt{M_1/3}$. We use values f_R0 = −10⁻⁴, −10⁻⁸, and − 10⁻¹² that correspond to models F4, F8, and F12, respectively. For these values of f_R0, the expansion history is indistinguishable from that in ΛCDM, as we have assumed. This has been studied in Hu & Sawicki (2007).

We note that M₂ and M₃, which determine the screening properties, are independent of k and positive. This is a non-trivial feature since Eqs. (27) and (28) may also be negative, as we saw for the symmetron model. In the f(R) case the functions Ms are positive, implying a normal screening (in contrast to what we found in the symmetron case: anti-screening).

3.3. Cubic Galileons and DGP models: Vainshtein mechanism

Cubic Galileons (Nicolis et al. 2009) and DGP (Dvali et al. 2000) models stem from different physical motivations, with different background dynamics, but they share a similar structure. Both are theories defined in the Jordan frame with the Klein–Gordon equation in the quasi-static limit given by

$$\begin{array}{c}\hfill \frac{1}{a^2}{\mathrm{\nabla }}_{\mathbf{x}}^2\chi =-Z_1(a)({({\mathrm{\nabla }}_{\mathbf{x}}^2\chi )}^2-{({\mathrm{\nabla }}_{\mathbf{x}i}{\mathrm{\nabla }}_{\mathbf{x}j}\chi )}^2)-Z_2(a)\rho ,\end{array}$$(53)

where Z₁ and Z₂ are model-dependent functions of time. We note that the scalar field becomes mass-less, such that M₁ = 0, and therefore the linear growth D₊ is scale independent. The screening on these models is provided by the Vainshtein mechanism, which arises from the nonlinear derivative terms in the Klein–Gordon equation. The functions M₂ and M₃ become (Aviles et al. 2018)

$$\begin{array}{cc}& M_2({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{Z_1(a)}{{\beta }^2}[k_1^2{k}_2^2-{({\mathbf{k}}_1·{\mathbf{k}}_2)}^2]\hfill \end{array}$$(54)

$$\begin{array}{cc}& M_3({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)=\frac{3Z_1(a)}{2a^2{\beta }^4{A}_0}[2{({\mathbf{k}}_1·{\mathbf{k}}_2)}^2{k}_3^2+({\mathbf{k}}_1·{\mathbf{k}}_2)k_1^2{k}_3^2\hfill \\ \hfill & -({\mathbf{k}}_1·{\mathbf{k}}_2){({\mathbf{k}}_1·{\mathbf{k}}_3)}^2-2({\mathbf{k}}_1·{\mathbf{k}}_2)({\mathbf{k}}_2·{\mathbf{k}}_3)({\mathbf{k}}_3·{\mathbf{k}}_1)].\hfill \end{array}$$(55)

The function M₃ appears when we transform from Eulerian to Lagrangian derivatives in Eq. (53). Given this structure, the screenings in both cubic Galileon and DGP are similar. Functions β, Z₁, and Z₂ depend on the specific model. For definiteness, we consider here DGP:

$$\begin{array}{cc}\hfill {\beta }^2& =\frac{1}{6}{[1-2ϵH{r}_{\mathrm{c}}(1+\frac{\dot{H}}{3H^2})]}^{-1},\hfill \end{array}$$(56)

$$\begin{array}{cc}\hfill Z_1& =\frac{2{\beta }^2{r}_c^2}{a^4},\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill Z_2& =4A_0{\beta }^2,\hfill \end{array}$$(58)

where ϵ = +1 for the self-acceleration branch and −1 for the normal branch. r_c is the crossover scale below which the theory behaves as a scalar tensor theory. For a side-by-side comparison of both models and for the expression of the functions β, Z₁, and Z₂ in cubic Galileons, see Barreira et al. (2013).

4. Matter power spectrum

For MG models with an early EdS phase, as we posit here, the linear matter power spectrum is given by

$$\begin{array}{c}\hfill P_{\mathrm{L}}(k,t)={(D_{+}(k,t))}^2{P}_{\mathrm{L}}^{\mathrm{\Lambda }\mathrm{CDM}}(k,t_0).\end{array}$$(59)

The building blocks for loop matter statistics are the functions

$$\begin{array}{cc}\hfill Q_1(k)& ={\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}{({\overline{D}}_{\mathrm{NS}}^{(2)}-{\overline{D}}_{\mathrm{S}}^{(2)})}^2{P}_{\mathrm{L}}(p)P_{\mathrm{L}}(|\mathbf{k}-\mathbf{p}|),}\hfill \end{array}$$(60)

$$\begin{array}{cc}\hfill Q_2(k)& ={\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}\frac{\mathbf{k}·\mathbf{p}}{p^2}\frac{\mathbf{k}·(\mathbf{k}-\mathbf{p})}{{|\mathbf{k}-\mathbf{p}|}^2}({\overline{D}}_{\mathrm{NS}}^{(2)}-{\overline{D}}_{\mathrm{S}}^{(2)})P_{\mathrm{L}}(p)P_{\mathrm{L}}(|\mathbf{k}-\mathbf{p}|),}\hfill \end{array}$$(61)

$$\begin{array}{cc}\hfill Q_3(k)& ={\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}\frac{{(\mathbf{k}·\mathbf{p})}^2}{p^4}\frac{{(\mathbf{k}·(\mathbf{k}-\mathbf{p}))}^2}{{|\mathbf{k}-\mathbf{p}|}^4}{P}_{\mathrm{L}}(p)P_{\mathrm{L}}(|\mathbf{k}-\mathbf{p}|),}\hfill \end{array}$$(62)

$$\begin{array}{cc}\hfill R_1(k)& ={\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}\frac{21}{10}\frac{D^{(3)s}(\mathbf{k},-\mathbf{p},\mathbf{p})}{D_{+}(k)D_{+}^2(p)}{P}_{\mathrm{L}}(p)P_{\mathrm{L}}(k)},\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill R_2(k)& ={\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}\frac{\mathbf{k}·\mathbf{p}}{p^2}\frac{\mathbf{k}·(\mathbf{k}-\mathbf{p})}{{|\mathbf{k}-\mathbf{p}|}^2}({\overline{D}}_{\mathrm{NS}}^{(2)}-{\overline{D}}_{\mathrm{S}}^{(2)})P_{\mathrm{L}}(p)P_{\mathrm{L}}(k)},\hfill \end{array}$$(64)

where the normalized growth functions in Eqs. (60) and (61) are evaluated as D⁽²⁾(p, k − p), while in Eq. (64), they are evaluated as D⁽²⁾(k, −p). These functions were introduced in (Matsubara 2008b) for EdS evolution and extended in (Aviles & Cervantes-Cota 2017) for MG. From them, power spectra and correlation functions in different resummation schemes can be obtained; see, for example (Carlson et al. 2013; Matsubara 2008c; Vlah et al. 2015b). In particular, the SPT power spectrum is defined as

$$\begin{array}{c}\hfill P_{\mathrm{SPT}}(k)=P_{\mathrm{L}}+P_{22}+P_{13},\end{array}$$(65)

and it can be shown that for ΛCDM (Matsubara 2008b; Vlah et al. 2015a) and for MG (Aviles et al. 2018), the following expressions hold:

$$\begin{array}{cc}\hfill P_{22}(k)& =\frac{9}{98}{Q}_1(k)+\frac{3}{7}{Q}_2(k)+\frac{1}{2}{Q}_3(k)\hfill \end{array}$$(66)

$$\begin{array}{cc}\hfill P_{13}(k)& =\frac{10}{21}{R}_1(k)+\frac{6}{7}{R}_2(k)-{\sigma }_L^2{k}^2{P}_{\mathrm{L}}(k),\hfill \end{array}$$(67)

where the one-dimensional variance of the linear displacement fields is

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{L}}^2=\frac{1}{3}{\delta }_{\mathit{ij}}⟨{\mathrm{\Psi }}_i(0){\mathrm{\Psi }}_j(0)⟩=\frac{1}{6{\pi }^2}{\displaystyle {\int }_0^{\infty }}\mathrm{d}p{P}_{\mathrm{L}}(p).\end{array}$$(68)

In Eq. (63) the third-order growth function D^(3)s is the solution to

$$\begin{array}{c}\hfill D^{(3)s}(\mathbf{k},-\mathbf{p},\mathbf{p})={(\widehat{\mathcal{T}}-A(k))}^{-1}({\mathcal{S}}_1+{\mathcal{S}}_2+{\mathcal{S}}_3),\end{array}$$(69)

where the label“s” means that D⁽³⁾(k₁, k₂, k₃) is symmetrized over its arguments; afterward, it is evaluated in the double-squeezed configuration given by k₃ = −k₂ = p and k₁ = k. We note from Eq. (63) that these are the only quadrilateral configurations –subject to momentum conservation: k = k₁ + k₂ + k₃– that survive in the one-loop computations. Expressions for the sources 𝒮_i are written in Appendix A. We have split the sources since in ΛCDM only 𝒮₁ is present, and in general, it is the dominant contribution. Source 𝒮₂ is a mix of frame-lagging and terms coming from the scale-dependent gravitational strength, which yields a small contribution. Meanwhile, 𝒮₃ is only composed by screenings; indeed, this is the only source containing the function M₃.

In the left panel of Fig. 1 we plot power spectra for a symmetron model with a_ssb = 0.33, m₀ = 1 h Mpc⁻¹, and β₀ = 1, and $\widehat{n}=\widehat{m}=0.5$. We plot this for the following cases: full-loop SPT (solid red curves); without screenings (dashed black curves); considering only source 𝒮₁ in Eq. (69) (dot-dashed gray curves); and the linear power spectrum (dotted blue curves). The lower panel in this figure shows the ratios of the different power spectra to the one-loop ΛCDM power spectra, for which we assumed the reference cosmology as given by WMAP 2009 best fits (Ω_m = 0.281, Ω_b = 0.046, h = 0.697, and σ₈ = 0.82). An interesting observation from these plots is that the perturbative screenings do not always act in the “screening direction”; that is, although the power spectrum with only source 𝒮₁ has a lower screening contributions than the full loop curves, it is actually closer to the GR power spectrum. Only 𝒮₁ is equivalent (up to frame-lagging terms) to keep only the γ₂ term in the expression of δI of Bose & Koyama (2016), who Fourier-expanded it in powers of matter overdensities, instead of in powers of the scalar field, as we do in Eq. (23). We interpret this fact to mean that the nonlinearities of the Klein–Gordon equation that are due to the couplings χδ also provide anti-screening effects. This behavior can be observed in the figures for symmetron models in Brax & Valageas (2013), but unfortunately, it is not discussed in that paper. An analogous plot for the F4 chameleon model is shown in the right panel of Fig. 1. We note that the different screening contributions always act in the same “direction”, driving the theory to GR; we furthermore note that in this case, the contributions due to sources 𝒮₂ and 𝒮₃ are almost negligible. In the following section we discuss these effects in more detail by considering the functions ${\overline{D}}_{\mathrm{S}}^{(2)}$ and ${\overline{D}}_{\mathrm{S}}^{(3)}$ for the symmetron model, defined in the Einstein frame, and for the DGP and f(R) models that are defined in the Jordan frame.

Fig. 1. SPT power spectra for the symmetron model (left panels) with assb = 0.33, m0 = 1 h/Mpc and β0 = 1, and for the F4 model (right panels), evaluated at redshift z = 0. We plot the linear theory, the full SPT, the SPT without screenings, and the SPT considering only S1 (source1) screenings (see text for details). In the lower panels we take their ratios to the one-loop SPT power spectrum for the ΛCDM model.

5. Screening growth functions

In this section we study the main features of the normalized second- and third-order growth functions ${\overline{D}}_{\mathrm{S}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)$ and ${\overline{D}}_{\mathrm{S}}^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$. Since a Dirac delta function accompanies the second-order growth functions and ensures that k = k₁ + k₂, these three wave vectors form a triangle. Therefore, by assuming statistical homogeneity and isotropy, the growth functions depend only on three positive numbers, for example, the lengths of the sides of the triangles. We can therefore write ${\overline{D}}_{\mathrm{S}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)={\overline{D}}_{\mathrm{S}}^{(2)}(k,k_1,k_2)$. Three triangle configurations are considered: equilateral, k = k₁ = k₂; orthogonal, $k_1=k_2=\sqrt{2}k$; and squeezed, k ≃ k₁, k₂ ≃ 0. There is a second squeezed configuration with k ≃ 0 and k₁ ≃ k₂ corresponding to very large scales, where the fifth force vanishes for mass-less theories and the screenings are zero. Analogously, we consider third-order screening growths given by Eq. (42) in the doubled-squeezed configuration explained above.

We stress that ${\overline{D}}_{\mathrm{S}}^{(2)}$ and ${\overline{D}}_{\mathrm{S}}^{(3)}$ with positive values will screen the fifth force, while negative values will act anti-screening it instead. Given this, the M functions of each model will determine their screening properties.

5.1. Chameleon screening in the Jordan frame: the f(R) case

First, we consider the f(R) Hu-Sawicki model for different f_R0 = −10⁻⁴, −10⁻⁸, −10⁻¹² corresponding to the F4, F8, and F12 models. In the upper panel of Fig. 2 we show plots for the second-order screening growth functions, ${\overline{D}}_{\mathrm{S}}^{(2)}$, evaluated at redshift z = 0 and for different triangular configurations. The vertical lines correspond to the screening wave number

$$\begin{array}{c}\hfill k_{M_1}\equiv a\sqrt{M_1(a)},\end{array}$$(70)

Fig. 2. ${\overline{D}}_{\mathrm{S}}^{(2)}$: contributions to the second-order normalized growth functions due to screenings for the F4, F8, and F12 models (upper panel), considering different triangle configurations. The vertical lines are located at $k=a\sqrt{M_1}$, showing the characteristic scale at which the screenings are present. The lower panel shows the third-order growth ${\overline{D}}_{\mathrm{S}}^{(3)}$ for different angles $x=\widehat{\mathbf{k}}·\widehat{\mathbf{p}}$ and with k = p.

which provides us with a rough estimate to characterize the scale at which the screening is present; it is close to the maximum screening growth of the largest triangular contribution (squeezed modes). We note that in f(R), the screening scale has a simple dependence $k_{M_1}\propto \sqrt{1/|f_{R0}|}$, as can be seen from Eq. (50). The lower panel of Fig. 2 shows the growth $D_S^{(3)}(\mathbf{k},-\mathbf{p},\mathbf{p})$ for the double-squeezed configuration with additionally |k|=|p| and for different values of the cosine angle $x=\widehat{\mathbf{k}}·\widehat{\mathbf{p}}$.

k_M1 should be understood as a phenomenological scale that does not show where the screening effects start to be present, but its usefulness is that it characterizes the screening for the models considered in this work. To show this at different redshifts, we study the effects of the background evolution on the screening growth. In Fig. 3 we show this for the F4 model at different redshifts z = 0, 3, and 10. The upper panel uses a ΛCDM background cosmology and the lower panel an EdS background evolution. In ΛCDM the screening curves are narrower and reach smaller maxima. This is expected because the cosmic acceleration attenuates the clustering of dark matter, and hence the nonlinear effects. Instead, in an EdS background, the growth pattern is preserved and only shifted by the scale k_M1 (denoted again by the vertical lines). Although this scale is useful, we note that it does not provide the range over which the screening growth is present. This will be more evident in the following section, where we study the symmetron model and show that the behavior of the growth is very different below and above k_M1.

Fig. 3. ${\overline{D}}_{\mathrm{S}}^{(2)}$ functions for the F4 model for redshifts z = 0, 3, and10, and considering different triangle configurations. The vertical lines are located at $k=a\sqrt{M_1}$, showing the characteristic scale at which the screenings are present. The upper panel shows a ΛCDM background evolution, and the bottom panel presents an EdS background evolution.

It is worthwhile to observe that in f(R), the nonlinearities of the Klein–Gordon equation lead to screening for any configuration. This is manifested in Figs. 2 and 3 because the screening growth functions always take positive values. To show this behavior, we show in Fig. 4 the second-order normalized growth functions ${\overline{D}}^{(2)}$ and ${\overline{D}}_{\mathrm{S}}^{(2)}$ (upper panel) and their ratios (bottom panel) for the equilateral and orthogonal configurations. Because the ratios always take values lower than unity, the standard screening behavior always occurs. This a consequence of Eq. (51), which shows that M₂ is always positive. Below we show that this is not necessarily the case for models defined in the Einstein frame.

Fig. 4. Normalized second-order growth functions ${\overline{D}}^{(2)}$ (solid curves) and ${\overline{D}}_{\mathrm{NS}}^{(2)}$ (dashed curves) for the F4 model at redshift z = 0, and considering different equilateral (brown) and orthogonal (red) configurations. The bottom panel shows the ratio between them. Since these take values lower than 1, the growth is always suppressed in the screened ${\overline{D}}^{(2)}$ case.

5.2. Symmetron screening

We consider the case of symmetron models. In Fig. 5 we show the second-order screening growth functions for models with $\widehat{n}=\widehat{m}=1/2$, a_ssb = 0.33, and different values of m₀ and β₀. The vertical lines correspond again to the characteristic scale k_M1. We are interested in also studying the growth functions at different redshifts, but since for $\widehat{n}=\widehat{m}$ the function M₁ becomes a constant, the corresponding k_M1 values are only rescaled by the scale factor. For this reason, we do it with a model with $\widehat{n}=0.25$ and $\widehat{m}=0.5$. The other parameters are fixed to a_ssb = 0.33, m₀ = 1 h Mpc⁻¹, and β₀ = 1. For this case, we show plots in Fig. 7.

Fig. 5. ${\overline{D}}_{\mathrm{S}}^{(2)}$ functions for symmetron models with assb = 0.33, and considering different triangle configurations. The vertical lines are located at $k=a\sqrt{M_1}$, showing the characteristic scale at which the screenings are present. The upper panel shows the model with β0 = 1 and different values of m0, while the bottom panel shows the specific model with m0 = 1 h Mpc−1 and different values of β0.

Fig. 7. ${\overline{D}}_{\mathrm{S}}^{(2)}$ functions for symmetron models with assb = 0.33, m0 = 1 h Mpc−1, β0 = 1, $\widehat{n}=0.25$ and $\widehat{m}=0.5$, and considering equilateral and orthogonal triangle configurations. The vertical lines are located at $k=a\sqrt{M_1}$, showing the characteristic scale at which the screenings are present.

Unlike in f(R) theories, the nonlinear terms in the Klein–Gordon equation for symmetron models lead to anti-screening effects (negative regions of the plots, cf. Figs. 5 and 7), or more precisely, there are configurations of interacting wave modes that instead of driving the theory toward GR drive the theory away from it. From Eq. (37), all triangular configurations with

$$\begin{array}{c}\hfill k_1^2+k_2^2<k_{\mathrm{AS}}^2(a),\end{array}$$(71)

where the anti-screening wave number is defined as

$$\begin{array}{c}\hfill k_{\mathrm{AS}}^2(a)\equiv \frac{{\kappa }_3\beta a^2}{{\beta }_2}-2m^2{a}^2,\end{array}$$(72)

will contribute with a negative source to the second-order screening as long as the RHS of the above equation is positive. We note that for the symmetron model, the anti-screening effects appear at scales below k_M1, although there is no a priori evident reason for this to happen, and it might be the case that other Einstein-frame posited theories show anti-screening effects above k_M1 as well.

In Fig. 6 we show an equivalent plot to Fig. 4, where in the upper panel we show the growth functions ${\overline{D}}^{(2)}$ and ${\overline{D}}_{\mathrm{S}}^{(2)}$, and in the lower panel we present their ratio. Regions where the ratios are greater than 1 correspond to wave-number configurations that act as anti-screening, that is, that enhance the MG fifth force.

Fig. 6. Normalized second-order growth functions ${\overline{D}}^{(2)}$ (solid curves) and ${\overline{D}}_{\mathit{NS}}^{(2)}$ (dashed curves) for a symmetron model with assb = 0.33, m0 = 1 h Mpc−1 and β = 1 at redshift z = 0. We consider equilateral (brown) and orthogonal (red) triangle configurations. The bottom panel shows their ratio. Regions with values greater than 1 correspond to wave numbers for which the growth is enhanced by the screening contributions.

The situation is similar when we consider third-order screening growth functions $D_S^{(3)}(\mathbf{k},-\mathbf{p},\mathbf{p})$, as shown in Fig. 8. We again note that while f(R) models always lead to positive screenings, symmetron models pose configurations that enhance the MG fifth force.

Fig. 8. Third-order growth function ${\overline{D}}_{\mathrm{S}}^{(3)}(\mathbf{k},-\mathbf{p},\mathbf{p})$ for symmetron models with β = 1 and assb = 0.33 and two different masses m0 = 0.1,  10 h Mpc−1. The vertical lines are located at kM1 for each model, showing the characteristic scale at which the screenings are present.

We may think of the case of the (original) chameleons (Khoury & Weltman 2004b) defined in the Einstein frame. Here we have a conformal factor C(φ) = e^φ/M ≃ 1 + φ/M and a potential V(φ) that decays with the scalar field. The linear conformal coupling implies that functions β_n vanish for n ≥ 2 (or at least, if we consider the full coupling e^φ/M, they are quite small compared to β) and by virtue of Eqs. (27) and (28), the functions M₂ and M₃ become scale independent. Moreover, the effective mass in these models decays with time, implying that function κ₃(a) is negative, while function κ₄(a) is positive. Then M₂ and M₃ are always positive. Hence, the growth functions will have the same kind of behavior as we observe in f(R). Equivalently, this can be observed at second order from the anti-screening wavenumber in Eq. (72), which is always negative, and hence no scale contributes to anti-screening configurations. In contrast to this, in symmetron models the effective mass grows with time and therefore κ₃ is positive, but β₂ is different from zero because the conformal coupling is quadratic in the scalar field. Henceforth, we observe this mix of screening and anti-screenings effects in PT.

5.3. Vainshtein screening: the DGP case

The source $S_{\text{S}}^{(2)}$ in Eq. (37) for the DGP and cubic Galileons depends only on the angle $x=\widehat{\mathbf{k}}·\widehat{\mathbf{p}}$ of the triangle configuration, and not on the scales k and p. This is because the mass of the scalar field is zero in these cases. By combining Eqs. (37) and (54), we obtain

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{S}}^{(2)}=\frac{Z_1(a)D_{+}^2(a)}{{\beta }^2(a)}(1-x^2),\end{array}$$(73)

where we note that the linear growth D₊(a) is also scale independent.

In Fig. 9 we show ${\overline{D}}_S^{(2)}$ for the normal branch DGP with a cross-over scale equal to the Hubble size, $r_c=H_0^{-1}$. For both cubic Galileons and DGP, it is a function depending only on x and time.

Fig. 9. ${\overline{D}}_{\mathrm{S}}^{(2)}$ and ${\overline{D}}_{\mathrm{S}}^{(3)}$ functions as a function of cosine angle $x=\widehat{\mathbf{p}}·\widehat{\mathbf{k}}$ for a DGP model with a cross-over scale rc = 1/H0. We state that the behavior is qualitatively the same as for cubic Galileons.