Cosmological gravity on all scales: simple equations, required conditions, and a framework for modified gravity
aa r X i v : . [ g r- q c ] A p r Cosmological gravity on all scales:simple equations, required conditions, and a framework for modified gravity
Daniel B. Thomas Jodrell Bank Centre for Astrophysics, School of Physics & Astronomy,The University of Manchester, Manchester M13 9PL, UK (Dated: April 29, 2020)The cosmological phenomenology of gravity is typically studied in two limits: relativistic per-turbation theory (on large scales) and Newtonian gravity (required for smaller, non-linear, scales).Traditional approaches to model-independent modified gravity are based on perturbation theory,so do not apply on non-linear scales. Future surveys such as Euclid will produce significant dataon both linear and non-linear scales, so a new approach is required to constrain model-independentmodified gravity by simultaneously using all of the data from these surveys. We use the higherorder equations from the post-Friedmann approach to derive a single set of “simple 1PF” (firstpost-Friedmann) equations that apply in both the small scale and large scale limits, and we ex-amine the required conditions for there to be no intermediate regime, meaning that these simpleequations are valid on all scales. We demonstrate how the simple 1PF equations derived here can beused as a model-independent framework for modified gravity that applies on all cosmological scales,and we present an algorithm for determining which modified gravity theories are subsumed underthis approach. This modified gravity framework provides a rigorous approach to phenomenologicalN-body simulations, and paves the way to consistently using all of the data from upcoming surveysto constrain modified gravity in a model-independent fashion.
I. INTRODUCTION
The current cosmological paradigm, ΛCDM, contains two hypothesised forms of matter: cold dark matter, anddark energy. The existing observational evidence for these is solely gravitational in nature, so it is natural to considerwhether modified gravity could be responsible for these observations. Furthermore, although General Relativity (GR)is well tested on smaller (Solar System) scales [1], it is less well tested on cosmological scales. Thus, there has beenmuch interest in using cosmology to constrain alternative gravity theories [2].One of the problems with constraining modified gravity theories is that the space of possible theories is vast. Onesolution to this is to use model-independent parameterisations [3–8], in order to constrain generic modifications toGeneral Relativity. These parameterisations have been used in many works to forecast results from future surveysand to derive constraints from existing data, see [2, 9] for reviews.Gravitational calculations are typically carried out in two regimes: on large scales using relativistic perturbationtheory (see e.g. [10]), and on small scales using the Newtonian limit and N-body simulations (see e.g [11]). Nearlyall of the existing parameterisations are built upon cosmological perturbation theory. Whilst this is accurate forpredictions on the largest scales in the universe, it breaks down on smaller scales where the density contrast is nolonger smaller than unity. The data from future surveys will include a significant contribution from non-linear scalessuch as Euclid [12], LSST [13] and the SKA [14], so to fully exploit these surveys for constraining modified gravitywe need to be able to parameterise modified gravity on these scales.The model-independent parameterisations can be used in several ways. The functional forms of these parametersin different theories can be used to evolve the perturbations in specific theories of gravity, and thus easily andquickly generate predictions for these theories, using common machinery (typically Boltzmann codes). In addition,it has been investigated how different functional forms or values of the parameters relate to different classes oftheories, in order to constrain the maximum number of theories with the minimal amount of calculation, such as theexpressions characterising the Horndeski classes of models (e.g. [15]). These parameterisations also allow for genericphenomenological modifications to gravity to be investigated, and thus for modified gravity to be investigated in amaximally model-independent way. In this way, these parameterisations act as a strong null test of ΛCDM. Oneparticular way to proceed along these lines is to bin the parameters into time and space dependent “pixels” that areallowed to have arbitrary values, and then perform a PCA analysis to determine which combinations are constrainedby current data (e.g. [16]) and future experiments (e.g. [17, 18]). This approach has several advantages, includingacting as a potential check on systematics [16], showing where the constraining power of specific data combinationsis, and requiring (almost) no assumptions as to the nature of the modifications to gravity.The forecasts [17, 18] restrict themselves to linear scales, as the parameterisations are not defined on non-linearscales and there is no theoretical justification for them on non-linear scales. This significantly affects how useful thesesurveys will be for constraining deviations from General Relativity. Conversely, works such as [16, 19] (and others) usethe linear theory parameterisations to justify the structure of modifications on non-linear scales; this is problematicas it isn’t clear a priori what these parameters mean and whether parameterising gravity in this same way across allscales is meaningful. We will resolve this issue in this work.The linear theory parameterisations have also been used to inspire phenomenologically modified N-body simulations[20–22] in an attempt to calculate model-independent modified gravity predictions on non-linear scales (see [23, 24] forsome earlier phenomenological modifications to N-body simulations). However, these works are lacking in theoreticaljustification for the modifications, in terms of when it is theoretically consistent and sufficient to parameterise thePoisson and slip equations on all scales, what the parameters mean, and whether they map to any theories of modifiedgravity. We will also address these issues in this work.One existing parameterisation that is valid on non-linear scales is [25], however it seems to only apply to theoriesthat fit into the Parameterised Post Newtonian (PPN) framework. It isn’t clear that many of the currently popularcosmological models fit into this approach, including theories with screening scales, such as f ( R ). Another recentdevelopment is [26], which parameterises the gravitational equations based on detailed spherical collapse calculationsfor known screening mechanisms. This has been implemented in N-body simulations, but it seems to primarilyinterpolate between known theories in particular regions of theory space, preventing it from achieving sufficientgenerality for comprehensive null tests and definitive statements about modified gravity from future surveys.Here we use the post-Friedmann formalism [27–32] to create a simple set of equations that apply on all cosmologicalscales and naturally incorporate both the Newtonian and perturbative limits. A key step in developing this set ofequations is the absence of an intermediate scale, i.e. a scale in the universe in which neither of these limits apply,which we discuss in detail. We will use this set of equations to create a theoretically justified set of parameterisedEinstein equations that are valid on all cosmological scales, including scales where the density contrast is arbitrarilylarge.We recap the Newtonian limit in cosmology and the post-Friedmann formalism in section II, before deriving ourequations in section III. This is used as the basis of a framework for modified gravity in IV. In section V we considersome of the practicalities of this framework, including presenting an algorithm for determining whether a givenmodified gravity theory fits into this framework. We conclude in section VI. Note that throughout, we use “allcosmological scales” to mean all scales where a perturbed FLRW metric with weak fields is a reasonable description ofthe spacetime, and by “non-linear scales” we mean scales of around 10Mpc and below, where the cosmological densitycontrast δ becomes greater than 1. We use “structurally non-linear” to refer to terms that have the form “variable × variable”, such as ρ~v or Φ . This is to distinguish from terms such as “ δ ”, which can be “non-linear” in the senseof being arbitrarily large, but is not of structurally non-linear form. II. THE NEWTONIAN LIMIT IN COSMOLOGY AND THE POST-FRIEDMANN FORMALISM
We will take the “Newtonian approximation” to be shorthand for the quasi-static (i.e. not evolving quickly in time,and therefore relatively down-weighting time derivatives and terms that are relevant on horizon scales), weak fieldand low velocity regime; these are the physical conditions typically associated with a c expansion. The Newtonianlimit is the leading order equations under these approximations.The Newtonian approximation is a good approximation for a ΛCDM universe (see e.g. [29, 33–37]), both onlinear scales and on smaller scales where the density contrast can become greater than unity. This is importantbecause Newtonian N-body simulations are a standard tool for cosmology. Recently, work has begun to investigatethe Newtonian limit in cosmology more thoroughly and build towards more complete general relativistic simulationsof our universe, see e.g. [29, 36–38].One of these approaches is the post-Friedmann formalism [27]. First proposed in [27, 28], the post-Friedmannformalism consists of a post-Newtonian type expansion of the Einstein equations in powers of the speed of light c ,altered compared to a “Solar-System” type expansion in order to apply to a FLRW cosmology. It is thus equivalentto taking a weak-field, low-velocity, quasi-static expansion of the Einstein field equations. This approximation workswell for ΛCDM cosmologies (see e.g. [29, 38] for calculations of the leading order corrections to this limit from N-bodysimulations, where these corrections are shown to be small).The starting point for the formalism is the perturbed FLRW metric, which is considered in Poisson gauge (see II C)and expanded up to order c − , g = − (cid:20) − U N c + 1 c (cid:0) U N − U P (cid:1)(cid:21) (2.1) g i = − aB Ni c − aB Pi c (2.2) g ij = a (cid:20)(cid:18) V N c + 1 c (cid:0) V N + 4 V P (cid:1)(cid:19) δ ij + h ij c (cid:21) . (2.3)Here, the two scalar potentials have each been split into their leading order (Newtonian) ( U N , V N ) and post-Friedmann( U P , V P ) components. The gauge freedom is chosen such that the vector potential appears in the 0 i part of the metric,and this has also been split up into B Ni and B Pi . The three-vectors B Ni and B Pi are both divergence-less, B Ni,i = 0and B Pi,i = 0. In addition, the tensor perturbation h ij is transverse and tracefree, h ii = h ,iij = 0.As is common for post-Newtonian type expansions, the terms of order c − and c − are considered to be leadingorder and the terms of order c − and c − are considered to be next-to-leading order. A factor of c − is attached totime derivatives. The matter content is taken to be pressure-less dust, the four-velocity of which is used to constructthe energy-momentum tensor, which is also expanded in powers of c (see [27] for details). The parameters describingthe pressure-less dust fluid are the background density ¯ ρ , the density contrast δ , and the peculiar velocity v i . Theassumption of pressure-less dust, applies to all of the results of this manuscript; in a universe dominated by radiation,the framework developed here will not apply . Note that non-linear behaviour in the velocity field manifests on larger scales than in the density field. However, an effective pressure caused by shell crossing (i.e. the breakdown of single streaming) is not a problem, as this would be higherorder and structurally non-linear ( ∼ v ) [39, 40], so would not contribute to the simple 1PF equations derived in section III. A. Leading order equations
The leading order ( c − ) equations from the post-Friedmann formalism correspond to the standard equations usedin a Newtonian cosmology, dδdt + v i,i a (1 + δ ) = 0 (2.4) dv i dt + ˙ aa v i = 1 a U N,i (2.5)1 c a ∇ V N = − πGc ¯ ρδ (2.6)2 c a ∇ ( V N − U N ) = 0 (2.7)Note that the time derivative in these equations is the convective derivative ( dA/dt = ∂A/∂t + v i A ,i /a , for anyquantity A ).A major advantage of the post-Friedmann formalism is that equations beyond those of standard Newtonian cos-mology can be derived. In particular, there is also a constraint equation for the vector potential that is present at the c − order, 1 c ∇ B Ni = − πG ¯ ρa c [(1 + δ ) v i ] | v , (2.8)where “ | v ” denotes the pure vector, i.e. divergence-free, component of the field. This vector potential is required inthe metric for the consistency of the Einstein equations in the Newtonian limit [27, 41], and it represents the leadingorder correction to Newtonian cosmology. It has been measured from N-body simulations [29, 31] and found to besmall, which provides a quantitative check of the Newtonian approximation, and thus the use of Newtonian N-bodysimulations, for GR+ΛCDM cosmology. The smallness of this vector potential will be important later on. B. Higher order equations
The post-Friedmann approach can be used to construct higher order equations, beyond the set of leading orderequations in the previous section. Keeping terms up to c order results in the first post-Friedmann order (“1PF”)equations, which describe structure formation on all scales (section VIII in [27]) and include terms that are higherorder in the c expansion, and terms that are structurally non-linear. In addition, no terms have been removed dueto any quasi-static type approximation. See [27] for full details of the relevant equations and notational conventions;here we present the equations with the background subtracted. The 00 Einstein equation is given by − c a ∇ V N + 1 c " ˙ aa ˙ V N + (cid:18) ˙ aa (cid:19) U N + 13 a ∇ V N − a V N,i V ,iN − a ∇ V P = 1 c πG ρδ + 1 c πG ρ (1 + δ ) v , (2.9)and the trace of the space-space Einstein equation is − c (cid:20) a ∇ ( V N − U N ) (cid:21) + 1 c ( − a ∇ ( V P − U P ) − a U N,k U ,kN − a V N,k V ,kN + 2 a U N,k V ,kN + 4 a V N ∇ ( V N − U N ) + 6 " ˙ aa ( ˙ U N + 3 ˙ V N ) + 2 ¨ aa U N + (cid:18) ˙ aa (cid:19) U N + ¨ V N = − πGc ¯ ρ (1 + δ ) v . (2.10)For reference, we present the conservation equations and 0 i Einstein equation in appendix A.Crucially, we note that the density contrast in these equations is not required to be small. These equations containboth the complete scalar and vector parts of linear cosmological perturbation theory, as well as containing a well-defined Newtonian limit that describes the growth of structure in a FLRW Universe on sub-horizon scales whether ornot the density contrast is large. Both of these limits are derived explicitly in [27].
C. A note on gauges
The Poisson gauge is a generalisation of the Newtonian gauge that gauge fixes the vector perturbations, so is anatural choice for the post-Friedmann approach. Moreover, as shown recently by [42], this is one of the few gaugesthat can be realised without the requirement of a small density contrast. For simplicity and clarity, we present all ofour results in this gauge, although we note that our modified Poisson equation (equation 3.16a; see later) containsthe gauge invariant density contrast. We leave to future work a full investigation of the impact of gauge choice (andconversion between gauges) on the equations presented here.
III. SIMPLE GRAVITATIONAL EQUATIONS ON ALL SCALES
Although the 1PF equations are quite complicated, there is a simpler version of them that still contains all of therequired information for the two limits (perturbation theory and Newtonian cosmology). In this section we will derivethese “simple 1PF” equations and comment on the conditions required for their validity. The resulting equations arenot just a curiosity: in section IV we will use them to create a simple parameterised framework for modified gravity,such that the same parameterisation applies on all cosmological scales.We start from the 1PF equations from [27], and we will use the “re-summed” potentials as defined in that work ψ P = − V N − c V P (3.1) φ P = − U N − c U P (3.2) ~ω P = ~B N + 2 c ~B P . (3.3)We rewrite the 1PF equations in terms of these resummed potentials (neglecting terms of order c − and higher), andcombine all of the “structurally non-linear” (see introduction) terms that are higher order in c into a single schematicterm, denoted [non-linear terms]. The time-time component of the 1PF field equations gives1 c a ∇ ψ P − c " ˙ aa ˙ ψ P + (cid:18) ˙ aa (cid:19) φ P + 1 c [non-linear terms] = 1 c πG ρδ , (3.4)where the “non-linear terms” includes both matter and metric variables; only their non-linear structure is importantfor our purposes here. The time-space component of the equations can be expressed in similar fashion1 c (cid:20) − a ∇ ω P i (cid:21) − c (cid:20) aa φ P,i + 2 ˙ ψ P,i (cid:21) + 1 c [non-linear terms] = 1 c πGaρv i − c πG ¯ ρaω P i . (3.5)We want to consider the scalar and pure vector parts of this equation separately, so we take the divergence to yield − c (cid:20) aa ∇ φ P + 2 ∇ ˙ ψ P (cid:21) + 1 c [non-linear terms] = 1 c πGa ( ρv i ) ,i , (3.6)for the scalar part of the equation. We write the vector part of this equation by removing the terms that are purelya divergence, leaving − c a ∇ ω P i + 1 c [non-linear terms] | v = (cid:20) c πGaρv i (cid:21) | v − c πG ¯ ρaω P i . (3.7)We will also require the space-space component of the 1PF field equations1 c (cid:20) a ∇ ∇ ( ψ P − φ P ) (cid:21) + 1 c [non-linear terms] = 0. (3.8)Finally, we have the Euler and continuity equations dv i dt + ˙ aa v i − φ P,i a + 1 c [non-linear terms] = 1 c a ddt ( aω P i ) (3.9) dδdt + v i,i a (1 + δ ) − c dψ P dt + 1 c [non-linear terms] = 0. (3.10)Focussing initially on the time-time and space-space equations, we now switch to Fourier space (a tilde “ ˜ A ” denotesquantities in Fourier space), in order to replace the derivatives and pull out the required structure of the scalarequations, − c a k ˜ ψ P − c " ˙ aa ˙˜ ψ P + (cid:18) ˙ aa (cid:19) ˜ φ P + 1 c h ^ non-linear terms i = 1 c πG ρ ˜ δ (3.11)1 c (cid:20) ˙ aa ˜ φ P + ˙˜ ψ P (cid:21) + 1 c h ^ non-linear terms i = 1 c πGk ¯ ρ ^ ([1 + δ ] v i ) ,i (3.12)1 c (cid:20) a k (cid:16) ˜ ψ P − ˜ φ P (cid:17)(cid:21) + 1 c h ^ non-linear terms i = 0. (3.13)We then substitute for the time derivatives in the time-time equation to arrive at two constraint equations, − c k ˜ ψ P + 1 c [non-linear terms] = 1 c πGa ¯ ρ ˜ δ + 1 c a ˙ aa πGk ¯ ρ ˜ θ (3.14)1 c (cid:20) a k (cid:16) ˜ ψ P − ˜ φ P (cid:17)(cid:21) + 1 c [non-linear terms] = 0. (3.15)We remind the reader that the density contrast is not required to be small in these equations, hence these equationsdescribe gravitational interactions on all scales. Note that when doing the substitution, the structurally non-linearparts of the ^ ([1 + δ ] v i ) ,i term have been moved into the “non-linear” terms, and we have defined θ = v i,i .To proceed further, we note that the terms denoted [non-linear terms] will vanish in both the small scale (Newtonian)limit, and the large scale (linear limit), because in the former case these terms are higher order in the c − expansion,and in the latter case because these terms are non-linear. This leaves us with the following equations1 c k ˜ ψ P = − c πGa ¯ ρ ˜ δ − c a ˙ aa πGk ¯ ρ ˜ θ (3.16a)1 c (cid:20) a k (cid:16) ˜ ψ P − ˜ φ P (cid:17)(cid:21) = 0, (3.16b)which are valid both on non-linear scales where the density contrast is large, and on larger scales where the Newtonianlimit (in particular the quasi-static approximation) is not valid. The structure of these equations is familiar: theyfollow the same structure as the standard linear perturbation equations, except that the potentials are the re-summedpotentials, and the density contrast is not required to be small. As is usual with such equations, the latter of theseequations can of course be replaced with the equation ˜ ψ P = ˜ φ P . We can combine these equations with the divergence-free part of the 1PF 0 i equation, and the 1PF Euler and continuity equations, from all of which the structurallynon-linear terms have also been dropped, yielding12 c a k ˜ ω P i = (cid:20) c πGa f ρv i (cid:21) | v − c πG ¯ ρa ˜ ω P i (3.17) dv i dt + ˙ aa v i − φ P,i a = 1 c a ddt ( aω P i ) (3.18) dδdt + v i,i a (1 + δ ) − c dψ P dt = 0. (3.19) In a slight abuse of notation, we will leave the continuity and Euler equations in real space, as this is the form they are typically seenin in the Newtonian limit and N-body simulations.
At linear order in perturbation theory, the scalars, vectors and tensors decouple from each other (e.g. [43]). Moreover,the vectors and tensors can be treated as higher order as long as (as is the case in ΛCDM) no matter sources thatactively generate linear vector and tensor perturbations (such as topological defects [44]) are present. Since the vectorpotential term in the Euler equation is also higher order in the c expansion, this suggests that this term can also bedropped in both of these two limits. This decouples the evolution of the scalar, vector and tensor components of themetric. It was shown in [29, 31] that the vector potential is negligibly small on all scales in ΛCDM, which justifiesthis choice. We now concentrate solely on the four equations containing only scalar variables, which we will refer toas the “simple 1PF equations”, 1 c k ˜ ψ P = − c πGa ¯ ρ ˜ δ − c a ˙ aa πGk ¯ ρ ˜ θ (3.20a)˜ ψ P = ˜ φ P (3.20b) dv i dt + ˙ aa v i − φ P,i a = 0 (3.20c) dδdt + v i,i a (1 + δ ) − c dψ P dt = 0. (3.20d)These four equations give a complete (leading order) description of structure formation in both the Newtonian andperturbative limits, for any cosmology where the Newtonian limit applies and there are no linear sources of the vectorand tensor in perturbation theory. If the vector potential is also small enough to not affect light bending significantlyat the level of our current observations [30] (which is a more stringent constraint than just being small enough for theNewtonian approximation to be good), then the scalar potential evolution described by these equations is sufficientto describe all large scale structure observables, including weak lensing, in these two limits.However, the real power of these simple equations is that they are of much more general validity than just the largeand small scale limits: In a realistic ΛCDM-like cosmology, these equations are actually valid on all cosmologicalscales. This is because, if the Newtonian limit is a good description of the Universe on all non-linear scales, and thegravitational dynamics are suitably “quasi-static” on linear scales well within the horizon (both of which known to betrue in a ΛCDM Universe) there is no “intermediate regime” in which the non-linear terms that we have discardedare important for structure formation. In other words, the two limits actually overlap, and thus the “simple 1PF”equations describe structure formation on all cosmological scales. We will now expand on this “intermediate regime”and its absence in ΛCDM. A. The intermediate regime
In this section we discuss the “intermediate regime” between the two limits usually used in cosmology, namelyperturbation theory (large scales) and the Newtonian limit (small scales). The intermediate regime refers to anyscales where both where non-linearity and higher order c effects are important for the gravitational equations, andthus the gravitational equations from neither limit apply. In ΛCDM, the lack of an intermediate regime is due to acombination of similarities between the equations in the two limits (including the decoupling of scalar, vector andtensor perturbations), the ability to ignore some terms in the perturbative equations on sub-horizon but linear scales,and the validity of the Newtonian approximation on all non-linear scales. This lack of an intermediate regime meansthat the simple 1PF equations apply on all cosmological scales in ΛCDM, which will also be important when usingthem as the basis of a parameterisation of modified gravity in section IV.Firstly, let is consider the similarities between the two limits. Both limits are weak field regimes, and thereforethe metric perturbations appear linearly in both limits. This manifests in the fact that the “re-summed” potentials,which appear linearly in the simple 1PF equations, describe all of the relevant metric degrees in the two limits. Ifthe Newtonian equations (for example) were a constraint equation for Φ and Ψ , then there would be no simplecorrespondence between the two limits. For there to be no intermediate regime on linear scales, the similarity betweenthe scalar equations needs to be even more precise. There need to be no linear terms in the Newtonian equationsthat aren’t present in the linear perturbation equations. In addition, the terms that are present in the linear limit but It is known that in a ΛCDM cosmology the perturbations to the metric are small [35, 37], except in the vicinity of black holes and othersuch objects, which is a negligible fraction of the spatial volume. not in the linearised Newtonian limit, need to be negligible below a certain scale ( k ∗ ) − where perturbation theoryis still valid, i.e. ( k ∗ ) − > ( k NL ) − , where k NL is the scale on which perturbations become non-linear. In GR, theseterms correspond to the time derivative of the metric potential in the continuity equation and the contribution of thevelocity to the gauge-invariant density contrast. It is well known that this is the case in a ΛCDM cosmology and thatthe linearised Newtonian limit equations match the scalar sub-horizon perturbation equations in a dust cosmology.Tensor-type potentials do not appear at leading order in either limit, but the vector potential is more complicated.In both limits it doesn’t influence the scalar dynamics, which is another key similarity between the two limits. Inprinciple, it is still present in the metric in the Newtonian limit and could influence observables. However, it has beenshown [30] that the vector potential in ΛCDM is too small to influence observables at the level of upcoming surveys;if this were not the case then the simple 1PF equations would need to be supplemented by the additional equationfor the vector potential in order to be a complete description of the quantities required for cosmological large scalestructure observables, such as light deflection.That there is no intermediate regime on linear scales is implicit in the use of cosmological N-body simulations, andthe known results that the output of these simulations matches that of perturbation theory on linear sub-horizonscales [33]. At this stage, there is just one check left to determine the lack of an intermediate regime, which is thatthe Newtonian approximation holds up to the non-linear scale, i.e. that none of the terms in the 1PF equations thatare higher order in the c expansion contribute on scales smaller than k NL . This has been tested using the smallnessof the induced vector potential, showing that the Newtonian approximation does indeed hold [29, 31, 38]. Thus wecan see that there is no intermediate regime where the gravitational dynamics require inclusion of terms that are bothstructurally non-linear, and higher order in the c expansion (a similar conclusion is reached differently in [37]).The importance of the smallness of the vector potential to the validity of the Newtonian approximation, andrestriction of the gravitational equations to dealing with scalars only, is why the vector potential will play an importantdiagnostic role when determining whether a particular modified gravity theory is contained within the approachpresented in section IV. B. The simple 1PF equations and N-body simulations
The simple 1PF equations are a generalisation of the usual Newtonian cosmological equations, so it is interestingto consider what the effect might be of implementing these equations in N-body simulations in lieu of the usualNewtonian equations. Comparing the simple 1PF equations to the standard Newtonian equations, there are two keydifferences: the time derivative of the potential in the continuity equation, and use of the gauge invariant densitycontrast in the Poisson equation. The requirements during the derivation that the Newtonian limit is good, and thatthere is no intermediate regime, mean that implementing these equations instead of the standard N-body equationswill make little difference except near the horizon scale where the c expansion breaks, and therefore won’t be relevantfor most of the applications of N-body simulations. Since these equations reduce to the standard perturbation theoryequations on large scales, then any differences between the simple 1PF and standard N-body equations on sub-horizonscales would already have been seen on linear scales in the simulation, in contrast to what is found (e.g. [33]). Thisresult will also apply to the modified gravity parameterisation (see later). However, we expect that on scales towardsthe horizon, the simple 1PF equations would affect the results from N-body simulations. For example, it might beexpected that the matter power spectrum in the simulation on the largest scales, which is known to be sensitive tothe gauge choice, might now match that calculated in the Poisson gauge. We leave a full investigation of these issuesto future work. IV. SIMPLE MODIFIED GRAVITY EQUATIONS ON ALL SCALES
The full 1PF equations govern gravitational structure formation on all scales: they don’t require either limit tobe specified, and they don’t require the density contrast to be small. As such, these equations can be used as thebasis of a framework to describe modified gravitational dynamics on all scales, for modified gravity theories where theNewtonian approximation is a good description of the small scale (non-linear) dynamics. A possible version of how Note that the vector potential can be small enough for the Newtonian approximation to be valid and still large enough that it must beconsidered for observables. Therefore the test described above for the vector potential to be ignored for all observables [30] is a stricterrequirement on the size of the vector potential than that required here for the Newtonian limit to be valid. these modified gravity equations might look like is1 c a ∇ V N + 1 c " ˙ aa ˙ V N + (cid:18) ˙ aa (cid:19) U N + 13 a ∇ V N − a V N,i V ,iN − a ∇ V P = 1 c πG A ( a, ~x ) ¯ ρδ + 1 c πG B ( a, ~x ) ¯ ρ (1 + δ ) v (4.1a) − c (cid:20) a ∇ ( V N − U N ) (cid:21) + 1 c ( − a ∇ ( V P − U P ) − a U N,k U ,kN − a V N,k V ,kN + 2 a U N,k V ,kN + 4 a V N ∇ ( V N − U N )+6 " ˙ aa ( ˙ U N + 3 ˙ V N ) + 2 ¨ aa U N + (cid:18) ˙ aa (cid:19) U N + ¨ V N = − πGc C ( a, ~x ) ¯ ρ (1 + δ ) v + D ( a, ~x ) (4.1b)1 c (cid:20) − a ∇ B Ni + 2 ˙ aa U N,i + 2 ˙ V N,i (cid:21) + 1 c (cid:20) − a ∇ B Pi + 4 ˙ aa U P,i + 4 ˙ V P,i + 2 ˙ V N U N,i + 4 ˙ aa U N U N,i +4 ˙ V N,i V N + 12 a B Ni ,k ( V N − U N ) ,k − a B Nk ,i ( U N + V N ) ,k + 1 a ∇ B Ni ( V N − U N ) + 12 a B Ni ∇ V N + 1 a B Nk V N,ki (cid:21) = 8 πGc E ( a, ~x ) ρav i + 8 πGc F ( a, ~x ) ρa (cid:8) v i (cid:2) v + 2( U N + V N ) (cid:3) − B Ni (cid:9) , (4.1c)where A, B, C, D, E, F are time and space dependent parameters that represent the effects of the modified gravity.We note that we are not presenting this as an optimal, rigorous or well-motivated parameterisation of these equations.It is merely intended to be a schematic representation of how such a set of equations might look. Considering thecomplexity of the equations (4.1), it is unclear how useful such a parameterisation would be. Moreover, motivatedby the derivation of the simple 1PF equations in the previous section, the fact that ΛCDM seems to be a reasonablemodel given our observations, and the lack of an intermediate regime in ΛCDM, we note that such a parameterisationis likely to be overkill.Instead, we will perform the same process as was performed in section III to write the parameterised gravitational1PF equations as1 c k ˜ ψ P + 1 c [parameterised non-linear terms] = − c πGa ¯ ρα (cid:16) a, ~k (cid:17) (cid:18) ˜ δ + ˙ aa c k e θ (cid:19) (4.2)˜ ψ P = β (cid:16) a, ~k (cid:17) ˜ φ P + 1 c [parameterised non-linear terms] (4.3)12 c a k ˜ ω P i + 1 c [parameterised non-linear terms] = (cid:20) c πGaγ (cid:16) a, ~k (cid:17) f ρv i (cid:21) | v − c πGǫ (cid:16) a, ~k (cid:17) ¯ ρa ˜ ω P i . (4.4)Here, we have switched to Fourier space and α (cid:16) a, ~k (cid:17) , β (cid:16) a, ~k (cid:17) , γ (cid:16) a, ~k (cid:17) and ǫ (cid:16) a, ~k (cid:17) are time- and space-dependentfunctions that we don’t expect to necessarily have simple functional forms in general. As before, the schematicterms representing the higher-order structurally non-linear terms are hiding the majority of the complexity in theseequations, except that they now contain additional complexity due to the extra modified gravity parameters thatmust necessarily be present. We now neglect these terms as in section III, leaving us with1 c k ˜ ψ P = − c πGa ¯ ρα (cid:16) a, ~k (cid:17) (cid:18) ˜ δ + ˙ aa c k e θ (cid:19) (4.5)˜ ψ P = β (cid:16) a, ~k (cid:17) ˜ φ P (4.6)12 c a k ˜ ω P i = (cid:20) c πGaγ (cid:16) a, ~k (cid:17) f ρv i (cid:21) | v − c πGǫ (cid:16) a, ~k (cid:17) ¯ ρa ˜ ω P i . (4.7)This is equivalent to a parameterisation of the gravitational equations in the simple 1PF equations derived in sectionIII, and should be used in conjunction with the Euler and continuity equations derived there (see equations (3.20)).Before continuing we note that, as before, the scalar, vector and tensor sectors are decoupled. In particular, the onlyeffect of the vector potential will be on photon trajectories, not on the evolution of structure in the universe. We canthus drop the vector equation as before, and work with just the scalar potentials. We comment further on this insection V A.0These parameterised equations describe the leading order gravitational dynamics in both the Newtonian and linearlyperturbed limits: the equations in each limit have the same structure and contain the same parameters, thereforeconsistent calculations can be carried out in the two limits. The only difference is that in the linear limit, the righthand side of the parameterised Poisson equation corresponds to the gauge invariant linear density contrast, as expectedin the Poisson gauge, and in the Newtonian limit the right hand side reduces to the standard Newtonian density.More importantly, for modified gravity cosmologies where there is no intermediate regime (as is the case for ΛCDM),we can go further and apply these equations on all cosmological scales, thereby having parameterised modified gravityequations where the parameters are meaningful on all scales. This means that the same parameters can be constrained using data that combines linear and non-linear scales ; this is the first time that such an approach has been shown tobe valid.For the sake of clarity, we recap the argument used to derive these equations. The starting point is the full 1PFequations, which describe structure formation on all scales and don’t require the density contrast to be small, thusthey are the natural choice for a parameterised set of equations to describe modified gravity on all scales. We rewritethese equations into a form that separates the terms that dominate in the large and small scale limits from theadditional terms that are required in an intermediate regime, and parameterise the resulting equations. We thenneglect the additional terms; crucially, the same pair of parameters describe the modified gravity effects in the twolimits. Then we note that in a modified gravity cosmology with no intermediate regime, these simple parameterisedequations describe the modified gravitational dynamics on all scales.
A. A convenient re-writing of the parameterised equations
Much investigation has been done into different choices of how to parameterise gravitational equations with thisstructure, and the relative merits of the different choices [2, 9]. As long as the system of equations has a constraintequation for one of the potentials, and an algebraic closure relation for the other, there are multiple choices thatsuffice. From here onwards we will work with the alternative parameters1 c k ˜ φ P = − c πGa ¯ ρG eff (cid:16) a, ~k (cid:17) (cid:18) ˜ δ + ˙ aa c k e θ (cid:19) (4.8a)˜ ψ P = η (cid:16) a, ~k (cid:17) ˜ φ P . (4.8b)The advantage of this choice is that it is φ P that governs the trajectories of massive particles, as can be seen fromthe 1PF Euler equation above. As a result, only G eff needs to be implemented into the evolution of the N-body code;the effects of η can all be included during post-processing steps. Therefore, these equations define the optimal wayto include modified gravity parameters in N-body simulations. Indeed, some theories have zero “slip” (e.g. [45]),meaning that η = 1, making these theories particularly easy to simulate in this approach. Another advantage ofmodified gravity simulations based on this framework is that there are no equations to solve for the additional fields,and thus the N-body simulations should take a similar amount of time as ΛCDM simulations. Note that, as discussedpreviously, only the leading order (in terms of c ) parts of these equations need to be implemented to model theoriesthat have a good Newtonian limit and no intermediate regime. These two equations replace the first two equationsin (3.20), and we refer to the resulting set of equations as the “Parameterised Simple 1PF” (PS1PF) equations.The G eff and η parameters will reduce to the functional forms already known for modified gravity theories on largescales, because the perturbative limit of the PS1PF equations reduces to the known linear theory parameterisation.Any complicated behaviour on small scales of additional fields that may be present is hidden in these functions,thus the time and space dependences are not expected to reduce to simple functional forms on small scales. Sinceour parameterisation reduces to that in [5] on super-horizon scales, this parameterisation satisfies the super-horizonconsistency condition [8, 46], as long as G eff and η tend towards finite, scale-independent functions on large scales.For convenience (because many cosmological observables are computed in Fourier space), and to connect to thelinear theory parameterisations, we have defined our modified gravity parameters in Fourier space. In principle it isalways possible to convert between parameters defined in Fourier space and in real space, and indeed when there isno scale dependence the conversion will be trivial. However, in the general scale-dependent case, converting betweenthe parameters becomes non-trivial because of the required convolutions.An important consequence of these equations is that they justify the use of phenomenological modified gravityN-body simulations, i.e. simulations based on a parameterised Poisson equation, rather than on the equations of1motion derived in a specific theory. The process used to arrive at these equations can be reversed to elucidate theconditions required for a parameterisation of the Poisson and slip equations to be a sensible and sufficient descriptionof the gravitational dynamics; see section V A. Under the phenomenological approach, the consequences of differentforms for G eff and η can be explored using N-body simulations for these values; we discuss a few possibilities for theseforms in section V C, including the possibility of binning G eff and η into time and space dependent “pixels” that canbe independently varied and constrained. Note that this does not justify the use of existing Halofit [47] and similarΛCDM semi-analytic fitting formulas for modified gravity cosmologies: N-body simulations still need to be employed. B. The Newtonian approximation and intermediate regime in modified gravity
The nature of the Newtonian approximation in modified gravity theories is important if the parameterisation ofthe simple 1PF equations is to apply to any “real” theories. Moreover, the equations are at their most powerful formodified gravity theories that (like ΛCDM) have no intermediate regime. Here we examine what is known aboutmodified gravity theories in these regards, and how it justifies the choice to use this approximation in building theparameterisation.The starting point of the majority of modified gravity cosmologies is to assume a weak field metric as in ΛCDM, andit is usual to consider no linear-order sources of vector and tensor perturbations. In addition, most of the currentlypopular theories are driven by the search for dark energy [2], and assume that cold dark matter is still the dominantmatter component of the universe. Cold dark matter is still non-relativistic in these modified gravity scenarios andthe velocities will still be small, suggesting that the Newtonian approximation will be good in these cosmologies.However, few detailed examinations of the Newtonian approximation in modified gravity have been performed onnon-linear scales. N-body simulations have been run for several specific modified gravity theories, including f(R) andDGP (see e.g. [48] and references therein) and Galileon theories [49, 50]. In the usual quasi-static and weak-fieldapproach used to derive the equations used in cosmological N-body simulations, several assumptions are made. Thedensity is taken to be the only quantity in the stress-energy tensor that contributes to the Einstein equations, andthe non-scalar elements of the metric are ignored. These assumptions amount to implicitly performing a c expansion,although it is usually not discussed in such terms; without this expansion, there is no reason to neglect the non-scalarparts of the metric when the density contrast is allowed to become arbitrarily large.As a result of the implicit c expansion, we note a couple of subtleties that are not normally considered, and whichare relevant to the framework developed here. Firstly, the expansion needs to be performed consistently for theorieswith dimensionful coupling constants, which means understanding the power of c that should be applied to thecoupling. We will discuss this in more detail in the context of the cubic Galileon in section V B 2. Secondly, for themajority of theories for which N-body simulations have been run, there remain two untested assumptions: that theNewtonian limit is valid on all non-linear scales, and that only the scalar potentials contribute to the bending oflight. The exception to this is Hu-Sawicki f ( R ) gravity [51], where these assumptions have been explicitly checked[32, 52]. In the former of these, the post-Friedmann formalism was applied to f ( R ) gravity, providing an alternativejustification for the equations used in N-body simulations and additionally deriving the equation for the leading ordercorrection, the vector potential. This potential was measured from f ( R ) N-body simulations and found to be similarlynegligible as in ΛCDM, showing that the Newtonian approximation is good on all non-linear scales for this theory,and that even on non-linear scales only scalar gravitational potentials are important.The work [26], which parameterises the Poisson equation in a way that generalises the behaviour seen in N-bodysimulations for specific theories, is implicitly built upon the same assumptions such as the absence of relevant vectorand tensor perturbations, and the lack of an intermediate regime. The derivation of the simple 1PF equations andthe lack of intermediate regime described above for ΛCDM show that the choice to parameterise the Poisson equationin [26] is valid for theories that fit into the framework considered here. In this sense, the parameterisation in [26] isa subset of the framework considered here, and can be used to guide functional forms of the PS1PF equations (seesection V C 2) for theories with screening mechanisms.One important area that has been investigated in detail is the Quasi-Static Approximation (QSA) [53–60], whichcorresponds to a sub-horizon approximation and treating time derivatives of the metric and extra field perturbationsas small compared to their spatial derivatives [55, 61], | ˙ X | ≤ HX (4.9a) |∇ X | ≫ H X or k ≫ H , (4.9b)where X represents either the metric or field perturbations. Typically it is found to be valid on scales that are smallerthan the sound horizon of the dark energy perturbations [58], but more interestingly it is found in [54, 55] that the2QSA is best in theories that are close to ΛCDM behaviour at the background level. There have also been severalinvestigations of the QSA in N-body simulations, see e.g. [48, 62–64]. This is an important consistency test of thesimulations, and it is a partial test of the Newtonian approximation in N-body simulations because it tests some ofthe approximations of the c expansion (corresponding to the relative down-weighting of time derivatives).The QSA is important as it is generally equivalent to neglecting the terms in the linear perturbation equations thatare not present in the Newtonian limit, therefore showing that the QSA is valid is a key part of showing that thereis no intermediate regime. This connection can be understood as follows. Newtonian theory is acausal, so it can onlyresemble GR in a “local inertial frame” or “causal region” where the system can be considered to be slowly varying(i.e. close-to-static) and in casual contact. This is part of the physical meaning of the c expansion. In GR, in amatter dominated universe, the “local inertial frame” is set by the Hubble horizon (and the speed of light/gravity c ), so it is expected that in this regime the Newtonian approximation could be a good description. The Newtonianapproximation is therefore naturally associated with the quasi-static conditions (equations (4.9)) applied to the metricperturbations. Put differently, the terms that are removed by the quasi-static conditions are ones that would not beexpected to be present in Newtonian theory. In modified gravity theories, there is an analogous “local inertial frame”or “causal region” set by the sound speed of the dark energy perturbations, within which the QSA is valid becausethe system can be considered to be slowly varying on these scales [58]. Therefore, within this region, the terms thatcan be neglected are those that would not be present in a “Newtonian-like” approximation to the modified gravitytheory. As a result, the existence of an intermediate regime on linear scales is unlikely for theories where the QSA isvalid on scales larger than k NL .We can reverse the argument from [55] for theories that are not close to ΛCDM behaviour. For these theories, thequasi-static approximation doesn’t hold on linear sub-horizon scales, so we expect that the Newtonian approximationis less likely to be valid. In particular, this means that there is more likely to be an intermediate regime and that thevector and tensor components of the metric are likely to be sourced.The combination of all these results suggest that the Newtonian approximation is a sensible choice for examiningmodified gravities that appear observationally similar to ΛCDM, and that we should expect the majority of popularmodified gravity theories to fit into this framework. In section V A we make this more concrete and present analgorithm for determining whether a particular theory fits into this framework. The results discussed here for Hu-Sawicki f ( R ) show that this theory does fit into the framework of the PS1PF equations. C. Background cosmology
The framework developed here focusses purely on the inhomogeneities and their evolution, and assumes a ΛCDMexpansion history. This choice has been made for both observational and theoretical reasons.Although it is possible to rule out some modified gravity theories using the background expansion history, manymodified gravity theories can be designed to mimic the ΛCDM expansion history (e.g. [51]). In addition, it hasbeen found that the quasi-static approximation (which is an indicator of whether a theory fits into this framework) istypically better for theories with a ΛCDM-like expansion history [55].Moreover, the background expansion is well understood theoretically and, as measurements of the expansion historybecome more precise, modified gravity theories that make meaningful modifications to the background will be ruled out(or ΛCDM will be). Conversely, the inhomogeneities in the Universe contain a much greater amount of information,and permit a much greater spectrum of possible behaviours that are currently much less understood. In particular,future surveys such as Euclid will deliver large data volumes on the inhomogeneities in the Universe on non-linearscales.
V. USING AND APPLYING THE MODIFIED GRAVITY FRAMEWORK
The discussion so far has focussed on the derivation of the simplified 1PF equations in standard gravity, and usingthem as the basis of a framework to describe modified gravity. In this section we delve further into the practicalitiesof using the framework presented here. We present an algorithm for determining whether a given theory fits into theframework proposed here, and we apply the post-Friedmann formalism to some example theories in order to illustratethe algorithm with respect to these theories. Finally, we consider some possible functional forms for the parameters.3
A. Algorithm for modified gravity theories
We begin by recapping the properties of ΛCDM elucidated earlier that result in there being no intermediate regime,and thus that the simple 1PF equations apply on all scales, • A weak field metric is appropriate on all cosmological scales; • There are no sources of vector and tensor perturbations in linear perturbation theory; • The linearised Newtonian scalar equations have no terms that don’t appear in linear perturbation theory; • There is a scale k ∗ > k NL , below which the terms that are present in linear perturbation theory, but not presentin the linearised Newtonian scalar equations, are negligible; • The Newtonian approximation is valid on all non-linear scales.For any theory where all of these criteria are met, the two limits will be valid, only scalar variables need to be evolvedto compute cosmological structure formation, and there will be no intermediate regime. This means that the simple1PF equations for the scalar sector can be utilised on all scales. If in addition the following extra criterion is satisfied, • The effect of the vector potential (e.g. on photon trajectories) is negligible for upcoming large scale structureobservations,then the additional constraint equation for the vector potential is not required, and the equations (3.20) are a completedescription of the variables required for computing large scale structure observables.The post-Friedmann vector potential is a key diagnostic here. Not only does it check the validity of the Newtonianapproximation on all non-linear scales, but it is required to be small for the scalar gravitational potentials to be theonly relevant gravitational degrees of freedom for structure formation. Therefore its smallness directly relates to thenon-existence of an intermediate regime. It is also the quantity that needs to be calculated for the extra criterion.We now use these criteria, the post-Friedmann vector potential, and the work carried out for f ( R ) in [32] tobuild an algorithm to determine whether a given theory fits under the assumptions and approximations behind theparameterisation of the simple 1PF equations in section IV.Algorithm1. Check that a weak field metric is appropriate on all cosmological scales.2. Derive the linear perturbation equations and check that there are no sources of vector and tensor perturbations.3. Derive the equations in the Newtonian limit using the post-Friedmann approach.4. Check for the existence of a scale k ∗ > k NL , below which the only terms that contribute significantly in thelinear perturbation equations are the same terms that are present in the linearised Newtonian equations (formost theories this will be equivalent to testing the QSA).5. Derive the equation for the post-Friedmann vector potential in the modified gravity theory, use it to calculate thevector potential from N-body simulations, and check that this is small enough for the Newtonian approximationto be valid on all non-linear scales. The extra criterion above for the effect of the vector potential can be tested easily from step 5,6. Determine if the effect of the vector potential on photon trajectories is negligible for upcoming large scalestructure observations; a reasonable rule of thumb for this is to compare the E-mode of cosmic shear from thevector potential to that generated by higher order deflections by the scalar [30, 65] To some extent, “how small is small enough?” is an arbitrary question. To be in line with claims about 1% precision cosmology, wepropose the conservative requirement P ω ( a, k ) < − P Φ ( a, k ) for all a, k . G eff and η functions for thattheory as a function of time and scale, by extracting the two sides of equations 4.8 at different snapshots, and takingthe ratio. Note that this requires outputting the gravitational potentials at each snapshot, in addition to the particleposition and velocity data that is the typical output of the simulations.In comparison to other approaches such as [25], the full procedure to determine whether a particular theory iscovered by this approach is somewhat more laborious. However, the trade-off is that the approach here does not a priori exclude any areas of theory space (known or otherwise), or rely on extrapolating from known theories andscreening mechanisms. Instead, the restrictions in theory space are based on the phenomenology of the theories; weexpect the approach here to apply to any theories that can describe a ΛCDM-like universe and are thus not alreadyruled observationally. As a result, a null detection of the parameters in this approach would be a stronger and moredefinitive constraint on modified gravity effects.From the arguments considered here, we can envisage modified gravity theories that wouldn’t be described by theparameterised simple 1PF equations. For example, a theory that generates vector or tensor perturbations at firstorder, perhaps due to having an additional vector or tensor field. Alternatively, a theory such as the massive gravitytheory in [67] might not fit into this framework, as it was pointed out in [68] that the small sound speed is likely tocause problems for the QSA. B. Post-Friedmann scalar-tensor theories and applying the algorithm
A substantial number of the modified gravity theories under investigation are scalar-tensor theories (see e.g. [2]),where an additional scalar field φ is included in the action. This includes f ( R ) gravity, which has already beenshown to fit into this framework. In this section we apply the post-Friedmann formalism to several more scalar-tensortheories, in preparation for examining them in terms of our algorithm. We will particularly be interested in thepost-Friedmann equation for the vector potential. Note that throughout this section we work in the Jordan frame.Before examining any theories in detail, we first consider what result we might expect. In the f ( R ) case, theconstraint equation for the vector potential is the same as in GR; the (inhomogeneous) additional scalar field does notenter. The vector potential still has a different value in f ( R ) because the scalar field causes the density and velocityfields to evolve differently. This result is not surprising: if the scalar field only enters the leading order Einsteinequations linearly, then no matter what derivatives are applied to it, a divergence-free quantity cannot be created,and thus the divergence-free part of the G i Einstein equation will not contain the scalar field. The scalar fields inthese theories typically behave like the scalar metric potentials in the presence of non-relativistic matter, so it is usualto expand them in powers of c as (see e.g. [1, 69]) φ = φ ( t ) (cid:18) c φ ( ~x, t ) + 1 c φ ( ~x, t ) + ... (cid:19) , (5.1)where φ is the background cosmological value of the field. Note that due to the weak field approximation, the leadingorder inhomogeneous part of the scalar field is assigned order c − , as was done in the post-Friedmann expansion for f ( R ) [32]. Assuming an expansion of this form, the inhomogeneous scalar field can only appear linearly in the leadingorder Einstein equations, and thus we expect the constraint equation for the vector potential to have the same formas in GR, except for a possible re-scaling of G by the background field.There is a possible loophole in this argument, which is that the scalar field could contain a dimensionful couplingconstant that is assigned a positive power of c in the expansion. A theory with such a coupling can have leading order Although this probably needs to be a truly divergenceless vector field, rather than a theory like Generalised Einstein-Aether [66], wheredue to the timelike constraint the vector field primarily contributes a scalar mode. i equation despite only having scalars present in the theory. We will discuss this case in more detail below whenapplying the post-Friedmann approach to the cubic Galileon.
1. Post-Friedmann Brans-Dicke
The prototypical scalar-tensor theory is Brans-Dicke theory [70], where the modified action is S = 116 πG Z d x √− g (cid:18) φR − wφ ( ∇ φ ) (cid:19) + S m , (5.2)where S m denotes the matter action and ( ∇ φ ) = g µν ∇ µ φ ∇ ν φ . For simplicity (and because it is true for the Galileoncase), we will take w to be a constant, and not a function of φ . This choice does not change the equation for thevector potential that we will derive. The Einstein equations and scalar field equation are given by φG µν − ∇ µ ∇ ν φ + g µν (cid:3) φ − wφ ∇ µ φ ∇ ν φ + w φ g µν ( ∇ φ ) = 8 πGc T µν (5.3)(3 + 2 w ) (cid:3) φ = 8 πGc T . (5.4)We now expand these equations according to the post-Friedmann description, see [27] for the components of theEinstein and stress-energy tensors, and we expand the scalar field as in equation (5.1). Taking terms up to order c − ,and subtracting off the background (homogeneous) terms, the equations become1 c ∇ V N = 4 πGφ c a ¯ ρδ + 12 c ∇ φ (5.5)1 c ∇ ( V N − U N ) = 12 c ∇ φ (5.6)1 c ∇ B Ni = − πG ¯ ρa φ c [(1 + δ ) v i ] | v (5.7)1 c (3 + 2 w ) ∇ φ = 8 πGφ c ¯ ρδ . (5.8)The structure of these equations match those obtained in PPN. As anticipated above, the constraint equation for thevector potential is the same as in GR, except for a re-scaling of G by the background field φ .
2. Post-Friedmann cubic Galileon
We now consider the cubic Galileon as an example of a more complicated scalar-tensor theory. Following [69], wedefine Y = − ( ∇ φ ) and write the action as S = Z d x √− g πG (cid:18) Rφ − wφ ( ∇ φ ) − α Y (cid:3) φφ (cid:19) . (5.9)The Einstein and scalar field equations are give by φG µν −∇ µ ∇ ν φ + g µν (cid:3) φ − wφ ∇ µ φ ∇ ν φ + w φ g µν ( ∇ φ ) + α (cid:18) g µν Y (cid:3) φ φ − Y ∇ ( µ ∇ ν ) φφ + (cid:3) φ ∇ µ φ ∇ ν φ φ (cid:19) = 8 πGc T µν (5.10)(2 w + 3) (cid:3) φ + α φ (cid:18) Y φ + 5 ∇ µ φ ∇ µ Yφ − ∇ µ φ ∇ µ (cid:3) φ − (cid:3) Y − ( (cid:3) φ ) − Y (cid:3) φφ (cid:19) = 8 πGc T , (5.11)where as expected these equations differ from the Brans-Dicke equations due to the new terms involving the dimen-sionful coupling α . We will proceed by applying the post-Friedmann approach as before, but first let us consider6what the lowest order contributions to each equation are from the new terms. The leading order contributions to thedifferent components of the Einstein equations are(00) : α ˙ φ ∇ φ a c φ (5.12a)(0 i ) : α ˙ φ c φ (cid:18) ˙ aa ˙ φ φ ,i − ˙ φ φ , i − ¨ φ φ ,i + φ φ ,i ∇ φ a (cid:19) (5.12b)( ij ) : α ˙ φ c φ (cid:18) δ ij ∇ φ − φ ,ij (cid:19) . (5.12c)The scalar field equation has many more terms at leading order; here we just note that they have order c , andschematically include terms such as (cid:0) ∇ φ (cid:1) . See appendix B for the full set of leading order post-Friedmannequations.We now consider the issue of which power of c should be assigned to α . In [69], detailed consideration was givento assigning orders to α ; the majority of the insight developed there will apply here, so we will just note a couple ofkey points in order to proceed: • The leading order behaviour of the metric is Newtonian, even deep inside the Vainshtein region. • The c order of α essentially comes from specifying that the leading order terms contribute to the leading orderscalar field equation. • The c order of α is independent of whether one is dealing with the “standard” or “dual” formulation, andwhich Vainshtein limit is being considered. As a result, the leading order gravitational and scalar field equationsin the c expansion are the same in the different Vainshtein limits. Following this logic for the post-Friedmann case considered here, we assign an order c to α . This agrees with theorder assigned in the PPNV case, as might be expected: it would be strange if this parameter required a different orderin the two expansions, because the PPNV expansion should result from taking the “no cosmological evolution” limitof the post-Friedmann equations (and then expanding with respect to the Vainshtein radius). Since α = M p Λ (see[69]), this is the same order that would be anticipated using M p ∼ c , (as appears in front of the stress-energy tensorin the Einstein field equations), which makes sense as we are not expanding in the Λ scale. The same considerationof the power of M p comprising α yields the correct c order for both the quartic and quintic Galileons as well. Thusthe assignment of c to α appears to be a consistent and robust choice.Using this order for α , the c terms in the scalar equation all contribute at the same order as the Brans-Dicketerms, by design, matching the result from [69]. In addition, the leading order additions to the gravitational equationselucidated in equations (5.12) will also contribute to the leading order gravitational equations. This last result isdifferent to the PPNV case, and arises due to the cosmological evolution of the background value of the field. Ofcourse, if this evolution is slow, we would expect the standard Brans-Dicke term to dominate.This might seem like an un-necessarily laborious process, just to select the leading order terms. However, we notethat in the PPNV case, the dimension of the coupling is important because it shows that at leading order there isno contribution to the gravitational equations beyond those of Brans-Dicke theory [69, 71]. In that case, just takingthe leading order Galileon terms (as decided by applying the weak field and quasi-static approximations with noconsideration for the order of α ) would result in the wrong equations. Deriving the equations in the consistentmanner described here allows us to be sure of whether these terms should be present at leading order. In the FLRWsituation considered here, unlike the Minkowski situation that is the basis of the PPNV approach, the backgroundscalar field is time-dependent, which is why we obtain a different result to the PPNV case. If it were constant, therewould be no Galileon terms contributing to the Einstein equations beyond the usual Brans-Dicke terms. Note thatthere is no such issue for the case of f ( R ) gravity, because there is no new dimensionful coupling constant in theaction, only the usual πGc term that is present in the GR case.We can obtain our desired result for the vector potential by taking the pure vector part of the cubic Galileoncontribution to the 0 i equation, which will contribute to the constraint equation for the vector potential. Knowing Taking the different Vainshtein limits is a simplification that allows a solution to be obtained for the scalar field. We do not require thisin the present work as we wish the equations to be valid at all locations, irrespective of proximity to the Vainshtein radius. α , we define α = c α ∗ in order to make the order of the extra terms explicit. From equation (5.12b),the additional terms are α ∗ ˙ φ c φ (cid:18) ˙ aa ˙ φ φ ,i − ˙ φ φ , i − ¨ φ φ ,i + φ φ ,i ∇ φ a (cid:19) . (5.13)We can see that these additional terms are of order c , so they will contribute to the leading order 0 i equation.The firstthree terms are exactly as anticipated earlier, and have no divergence-free part. However, consider the fourth term: ∇ × ( A ∇ B ) contains ∇ A × ∇ B , so the curl will not vanish. Instead we will have a term of the form ∇∇ φ × ∇ φ .The full vector constraint equation thus deviates from the Brans-Dicke form, c ∇ B Ni = − πG ¯ ρa φ c [(1 + δ ) v i ] | v − α ∗ ˙ φ c a φ (cid:2) φ ,i ∇ φ (cid:3) | v . (5.14)This equation can be used to extract the vector potential from cubic Galileon N-body simulations. This quantity hasseveral uses: firstly it is a consistency check of the approximations used in this approach and is part of the algorithmdescribed above for determining whether the cubic Galileon fits into the modified gravity framework described here,which we will use in the next section. In addition, this vector potential should be taken into account when ray-tracingthrough N-body simulations [30] and can generate unique lensing signals [30, 72]. We leave it to future work toperform this extraction.We remind the reader that the vector potential equation is not required to have the Brans-Dicke (or GR) form inorder for the theory to fit into the modified gravity approach described here. However, as we will describe in thenext section, for theories where the equation matches that in Brans-Dicke, it is possible to construct a plausibilityargument for determining whether these theories fit into the framework described here that does not require runningN-body simulations. For theories where the vector equation does not match Brans-Dicke (or GR), no such shortcutcan be made.
3. Applying the algorithm
Having applied the post-Friedmann formalism to these theories, we can now start examining them using the algo-rithm in section V A. We note that points 1 and 2 of the algorithm naturally come out of the extensive work on linearperturbation theory in scalar tensor theories, see e.g. [2, 57, 73].Steps 3-4 require us to compare the equations in the perturbative limit with the linearised equations from theNewtonian limit. From the equations in the literature, or applying perturbation theory to the equations above, wecan see that there are additional terms present in the perturbation equations. Even for Brans-Dicke theory, there are asubstantial number of extra terms [2, 73], and there are even more for the cubic Galileon [57]. However, the weak-fielddescription that applies in both limits strongly shapes the form of these extra terms: since the metric potentials andthe inhomogeneities in the scalar field are all order one, there can be no structurally non-linear terms, and all of theadditional terms will involve time derivatives, either because of pre-factors such as ˙ aa or ˙ φ , or because they containtime derivatives of the perturbations themselves. For example, some of the additional Brans-Dicke contributions tothe 00 equation have forms such as [2] ˙ aa ˙ φ ; ωφ ˙ φ φ ! . (5.15)These additional terms are of course exactly the terms that are neglected in the QSA, as discussed in section IV B.This simplifies step 4 of the algorithm considerably, due to the substantial effort that has gone into examining the QSA(see section IV B above, or e.g. [74] for a recent summary): as long as the terms that are missing in the linearisedNewtonian equations are those associated with the QSA (as they are in this case), then step 4 simply reduces tochecking the QSA, which is known to be fine for the theories considered here (see e.g. [57]). So these theories satisfysteps 1-4 of the algorithm.We now move on to considering steps 5-6 of the algorithm. These tests involve checking the smallness of thevector potential on non-linear scales, and comparing it to the scalar quantities in order to determine the validity This difference to Brans-Dicke and f ( R ) gravity can be understood in terms of the phenomenology of the theory. The dimensionfulcoupling constant results in a theory that clusters much more strongly. C. Different functional forms for the parameters
Whilst on linear scales the parameters G eff and η can be analytically calculated for specific theories or classesof models, on non-linear scales they are best thought of as emergent prescriptions for the complicated non-lineardynamics within particular theories. Here we mention a few sensible choices for the functional forms of the parametersin the simple 1PF equations, including some known cases and some possible functional dependencies that are worthinvestigating.
1. Linear scales
In the linear perturbation limit, the parameters in the simple 1PF equations will reduce to the known functionalforms in linear theory.
2. Inclusion of specific scales
In [26], building on earlier work [75], the authors propose specific functional forms for the parameterised Poissonequation, and show that specific choices for some of the parameters can recreate the phenomenology of particulartheories that have been studied. These particular forms are based around detailed studies of spherical collapse intheories with screening mechanisms. The modifications to the Poisson equation take the form b kk a (cid:20) (cid:18) k k (cid:19) a (cid:21) /b − ! , (5.16)where a and b are time varying constants, and k is a particular scale associated to the modifications (which is alsoexpected to be time dependent).An alternative approach was suggested in [3], where the authors use k π P ( a, k ) = ∆( a, k ) to measure the degree ofnon-linearity at a particular scale. This defines the scale of non-linearity k NL ( a ) through ∆( a, k NL ) = 1. The validity9of this form for predicting the power spectrum has been tested (e.g. [76]; see [75] for a nice summary) and foundto be good for weakly non-linear scales. It would be interesting to investigate functional forms for the modificationsto the Poisson equation as a function of k NL ( a ): these functions would naturally incorporate time dependence andreturn to the linear forms on the appropriate scales.
3. Environmental dependence
Rather than explicitly including specific scales, theories with screening mechanisms can be understood by allowingthe modifications to the Poisson equation to have additional environmental dependence. This was used in [77] andtested in [63] in the form ∇ Ψ = 4 πGa ¯ ρδ (1 + ∆( δ )) , (5.17)where the modification ∆( δ ) is explicitly dependent on the local density, which naturally occurs in the Vainshteinscreening mechanism. This form of ∆( δ ) was constructed to represent DGP, and included the “crossover scale” r c explicitly, and an approximate form was found to be given by ∆( δ ) = β ( t ) (cid:0) √ . δ − (cid:1) (where β ( t ) is a timedependent function of the background cosmology).Another form of environmental dependence that would be interesting to explore is a dependence on the localacceleration, which links to MONDian [78] theories, although an examination of how the post-Friedmann approachapplies to such theories is beyond the scope of this paper.
4. Maximally phenomological pixels
As described in the introduction, a standard approach to exploring constraints on parameters whose functionalforms are not known, in a maximally model independent fashion, is to express the parameters as piecewise constantfunctions in a set of bins, or pixels, in time and space. This technique has been successfully used for modified gravity,[16–18] and has also been used to constrain the dark matter equation of state in a model-independent fashion [79].As well as being the most model-independent way to analyse observational data, pixelised functions can be usefulto understand the phenomenology of different theories, by isolating the ranges in scale and time where changesto the gravitational law result in specific effects. Similarly, when used in combination with a specific dataset (orforecasted dataset), a pixelised approach can determine at which times and scales that dataset is most sensitive tothe modifications.
5. Hybrid pixels
Since the functional forms of the parameters are well known on linear scales, a hybrid scheme can be defined, wherethe pixel boundaries are defined with respect to k NL ( a ) rather than k . This reduces the number of pixels required bymaking the time dependence implicit, and naturally returning to known forms on linear scales. Schematically, suchhybrid pixels would have the form G eff = G eff,linear ; k < A k NL ( a ) (5.18) G eff = G eff,1 ; A k NL ( a ) < k < A k NL ( a ) (5.19) G eff = G eff,2 ; A k NL ( a ) < k < A k NL ( a ) (5.20) G eff = G eff,3 ; A k NL ( a ) < k , (5.21)for a set of constants { A , A , A , G eff,1 , G eff,2 , G eff,3 } with no time dependence. It would be interesting to see whetherthis approach more easily captures the phenomenology of existing theories. VI. DISCUSSION
We have presented a simplified set of equations that describe structure formation on all cosmological scales andthat do not require the density contrast to be small, equations 3.16. This set of equations depends on the validity of0the Newtonian approximation in a given cosmology, and requires that there is no “intermediate regime” (where boththe Newtonian limit and perturbation theory fail). This is the case for a matter dominated ΛCDM cosmology.We have used these equations to create a simple framework for modified gravity (equations 4.8) that is valid on all cosmological scales, and which is known to include currently popular models such as f ( R ). This set of equations isimportant for upcoming cosmological surveys such as Euclid [12], as it allows them to constrain model-independentmodified gravity by simultaneously using the data from all scales, not just linear scales. This maximises the constraintson model-independent modified gravity that are achievable by these surveys.In section V A we present an algorithm for investigating the validity of our approximations for a given modifiedgravity theory, in particular whether there is an intermediate regime. This algorithm determines whether a theoryof gravity is subsumed under this approach, and therefore whether constraints determined using this parameterisedapproach are applicable to that theory of gravity. We have illustrated this algorithm by applying the post-Friedmannapproach to the cubic Galileon, and calculating the equation for the vector potential that is required for the algorithm.We note that throughout this work, we are primarily focussed on modified gravity theories that replace dark energy,and thus on cosmologies that still contain dark matter. We will investigate the application of the post-Friedmannformalism to modified gravity theories without dark matter in future work.An important corollary of the equations and framework developed here is that it is theoretically consistent to performphenomenological N-body simulations based on a modified Poisson equation, including work that was performed (butnot justified) in [20–22]. In particular, the derivation outlined here (and used to construct the algorithm in sectionV A) elucidates the conditions required for a parameterisation of the Poisson and slip equations to be a sensible andsufficient description of the gravitational dynamics.We have presented some functional forms for the parameters that could be used in simulations, including the possi-bility of exploring the modified gravity parameter space by binning these parameters into time and space dependent“pixels”, and examining the resulting phenomenology. This pixel approach enables forecasts (and data analysis) tounderstand which scales and redshifts different surveys (and types of surveys) are sensitive to.The set of equations in this work is the starting point for bringing the full constraining power of future surveys tobear on modified gravity theories in a model-independent way. This will provide a significant null test of the ΛCDMparadigm and enable general conclusions about gravity to be drawn from the data, rather than statements aboutspecific models or areas of parameter space. ACKNOWLEDGEMENTS
The author thanks Ali Mozaffari, Johannes Noller and Michael Kopp for useful comments on the manuscript,and Richard Battye, Francesco Pace, Sankarshana Srinivasan and Marco Bruni for useful discussions. The authoracknowledges support from Science and Technology Facilities Council (STFC) grant ST/P000649/1.
Appendix A: 1PF equations
For brevity in section II we present a subset of the 1PF equations. Here we present the additional equations thatare required for section III.1PF Continuity equation, dδdt + v i,i a ( δ + 1) + 1 c (cid:20) ( δ + 1) (cid:18) a v j U N,j − ˙ aa v + 3 dV N dt (cid:19)(cid:21) = 0 . (A1)1PF Euler equation, dv i dt + ˙ aa v i − U N,i a + 1 c (cid:20) v i ddt ( U N + 2 V N ) + 2 a U N,i ( U N + V N ) − a v V N,i − a U P,i + 1 a v i v j U N,j − ˙ aa v v i − a ddt ( aB Ni ) + B Nj,i v j a (cid:21) = 0 . (A2)11PF 0 i Einstein equation,1 c (cid:20) − a ∇ B Ni + 2 ˙ aa U N,i + 2 ˙ V N,i (cid:21) + 1 c (cid:20) − a ∇ B Pi + 4 ˙ aa U P,i + 4 ˙ V P,i + 2 ˙ V N U N,i + 4 ˙ aa U N U N,i +4 ˙ V N,i V N + 12 a B Ni ,k ( V N − U N ) ,k − a B Nk ,i ( U N + V N ) ,k + 1 a ∇ B Ni ( V N − U N ) + 12 a B Ni ∇ V N + 1 a B Nk V N,ki (cid:21) = 8 πGc ρav i + 8 πGc ρa (cid:8) v i (cid:2) v + 2( U N + V N )] − B Ni (cid:9) . (A3) Appendix B: The leading order post-Friedmann cubic Galileon equations
For posterity, we present here the complete leading order gravitational and scalar field equations for the cubicGalileon. According to the post-Friedmann approach, these are the equations that should be implemented in N-bodysimulations in order for the non-relativistic truncation of the equations to be consistent.1 c ∇ V N + α ∗ ˙ φ ∇ φ c φ = 4 πGφ c a ¯ ρδ + 12 c ∇ φ (B1)1 c ∇ ( V N − U N ) + α ∗ ˙ φ c φ ∇ φ = 12 c ∇ φ (B2)1 c ∇ B Ni = − πG ¯ ρa φ c [(1 + δ ) v i ] | v − α ∗ ˙ φ c a φ (cid:2) φ ,i ∇ φ (cid:3) | v (B3)1 c (3 + 2 w ) ∇ φ + α ∗ φ " ˙ φ a c ∇ φ + φ c a ( φ ,ki φ ,ik ) − c a φ ( ∇ φ ) + 2 c a φ ¨ φ ∇ φ = 8 πGφ c ¯ ρδ . (B4)We note these equations may not be consistent with those used in [49]. We leave a detailed comparison between thetwo approximation schemes, and resulting sets of equations, to future work. [1] C. M. Will, Living Rev. Rel. , 3 (2006), arXiv:gr-qc/0510072 [gr-qc].[2] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, phys. rep. , 1 (2012), arXiv:1106.2476.[3] W. Hu and I. Sawicki, Phys. Rev. D , 104043 (2007), arXiv:0708.1190.[4] C. Skordis, Phys. Rev. D , 123527 (2009), arXiv:0806.1238 [gr-qc].[5] G.-B. Zhao, L. Pogosian, A. Silvestri, and J. Zylberberg, Phys. Rev. D , 083513 (2009), arXiv:0809.3791.[6] S. F. Daniel, E. V. Linder, T. L. Smith, R. R. Caldwell, A. Cooray, A. Leauthaud, and L. Lombriser,Phys. Rev. D , 123508 (2010), arXiv:1002.1962.[7] M. A. Amin, R. V. Wagoner, and R. D. Blandford, MNRAS , 131 (2008), arXiv:0708.1793.[8] E. Bertschinger and P. Zukin, Phys. Rev. D , 024015 (2008), arXiv:0801.2431.[9] M. Ishak, Living Rev. Rel. , 1 (2019), arXiv:1806.10122 [astro-ph.CO].[10] K. A. Malik and D. Wands, Phys. Rept. , 1 (2009), arXiv:0809.4944 [astro-ph].[11] E. Bertschinger, Ann. Rev. Astron. Astrophys. , 599 (1998).[12] A. Refregier, A. Amara, T. D. Kitching, A. Rassat, R. Scaramella, J. Weller, and f. t. Euclid Imaging Consortium, ArXive-prints (2010), arXiv:1001.0061 [astro-ph.IM].[13] Z. Ivezic and for the LSST Collaboration, ArXiv e-prints (2008), arXiv:0805.2366.[14] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. W. Lazio, IEEE Proceedings , 1482 (2009).[15] E. Bellini and I. Sawicki, JCAP , 050 (2014), arXiv:1404.3713 [astro-ph.CO].[16] G.-B. Zhao, T. Giannantonio, L. Pogosian, A. Silvestri, D. J. Bacon, K. Koyama, R. C. Nichol, and Y.-S. Song,Phys. Rev. D , 103510 (2010), arXiv:1003.0001 [astro-ph.CO].[17] A. Hojjati, G.-B. Zhao, L. Pogosian, A. Silvestri, R. Crittenden, and K. Koyama, Phys. Rev. D , 043508 (2012),arXiv:1111.3960 [astro-ph.CO].[18] G. Zhao, D. Bacon, R. Maartens, M. Santos, and A. Raccanelli, Advancing Astrophysics with the Square Kilometre Array(AASKA14) , 165 (2015), arXiv:1501.03840.[19] F. Simpson, C. Heymans, D. Parkinson, C. Blake, M. Kilbinger, J. Benjamin, T. Erben, H. Hildebrandt, H. Hoekstra,T. D. Kitching, Y. Mellier, L. Miller, L. Van Waerbeke, J. Coupon, L. Fu, J. Harnois-D´eraps, M. J. Hudson, K. Kuijken,B. Rowe, T. Schrabback, E. Semboloni, S. Vafaei, and M. Velander, MNRAS , 2249 (2013), arXiv:1212.3339. [20] W. Cui, P. Zhang, and X. Yang, Phys. Rev. D , 103528 (2010), arXiv:1001.5184 [astro-ph.CO].[21] D. B. Thomas and C. R. Contaldi, ArXiv e-prints (2011), arXiv:1112.6378 [astro-ph.CO].[22] Y. Zhang, P. Zhang, X. Yang, and W. Cui, Phys. Rev. D , 023521 (2013), arXiv:1301.3255 [astro-ph.CO].[23] H. F. Stabenau and B. Jain, Phys. Rev. D , 084007 (2006), astro-ph/0604038.[24] I. Laszlo and R. Bean, Phys. Rev. D , 024048 (2008), arXiv:0709.0307.[25] T. Clifton and V. A. Sanghai, Phys. Rev. Lett. , 011301 (2019), arXiv:1803.01157 [gr-qc].[26] F. Hassani and L. Lombriser, arXiv e-prints , arXiv:2003.05927 (2020), arXiv:2003.05927 [astro-ph.CO].[27] I. Milillo, D. Bertacca, M. Bruni, and A. Maselli, Phys. Rev. D , 023519 (2015), arXiv:1502.02985 [gr-qc].[28] I. Milillo, Linear and non-linear effects in structure formation , Ph.D. thesis, University of Portsmouth (2010).[29] M. Bruni, D. B. Thomas, and D. Wands, Phys.Rev.
D89 , 044010 (2014), arXiv:1306.1562.[30] D. B. Thomas, M. Bruni, and D. Wands, JCAP , 021 (2015), arXiv:1403.4947 [astro-ph.CO].[31] D. B. Thomas, M. Bruni, and D. Wands, Mon. Not. Roy. Astron. Soc. , 1727 (2015), arXiv:1501.00799 [astro-ph.CO].[32] D. B. Thomas, M. Bruni, K. Koyama, B. Li, and G.-B. Zhao, JCAP , 051 (2015), arXiv:1503.07204 [gr-qc].[33] N. E. Chisari and M. Zaldarriaga, Phys. Rev. D , 123505 (2011), arXiv:1101.3555 [astro-ph.CO].[34] S. R. Green and R. M. Wald, Phys. Rev. D , 063512 (2012), arXiv:1111.2997 [gr-qc].[35] S. R. Green and R. M. Wald, Classical and Quantum Gravity , 234003 (2014), arXiv:1407.8084 [gr-qc].[36] J. Adamek, D. Daverio, R. Durrer, and M. Kunz, JCAP , 053 (2016), arXiv:1604.06065 [astro-ph.CO].[37] W. E. East, R. Wojtak, and T. Abel, Phys. Rev. D , 043509 (2018), arXiv:1711.06681.[38] J. Adamek, R. Durrer, and M. Kunz, Classical and Quantum Gravity , 234006 (2014), arXiv:1408.3352.[39] D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, JCAP , 051 (2012), arXiv:1004.2488 [astro-ph.CO].[40] A. Erschfeld and S. Floerchinger, JCAP , 039 (2019), arXiv:1812.06891 [astro-ph.CO].[41] M. Kopp, C. Uhlemann, and T. Haugg, JCAP , 018 (2014), arXiv:1312.3638 [astro-ph.CO].[42] T. Clifton, C. S. Gallagher, S. Goldberg, and K. A. Malik, Phys. Rev. D , 063530 (2020), arXiv:2001.00394 [gr-qc].[43] H. Kodama and M. Sasaki, Progress of Theoretical Physics Supplement , 1 (1984).[44] R. Durrer, M. Kunz, and A. Melchiorri, Phys. Rept. , 1 (2002), arXiv:astro-ph/0110348 [astro-ph].[45] E. V. Linder, JCAP , 005 (2018), arXiv:1801.01503 [astro-ph.CO].[46] E. Bertschinger, Astrophys. J. , 797 (2006), arXiv:astro-ph/0604485 [astro-ph].[47] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas,G. Efstathiou, and H. M. P. Couchmann (VIRGO Consortium), Mon. Not. Roy. Astron. Soc. , 1311 (2003),arXiv:astro-ph/0207664 [astro-ph].[48] H. A. Winther, F. Schmidt, A. Barreira, C. Arnold, S. Bose, C. Llinares, M. Baldi, B. Falck, W. A. Hellwing, K. Koyama,B. Li, D. F. Mota, E. Puchwein, R. E. Smith, and G.-B. Zhao, MNRAS , 4208 (2015), arXiv:1506.06384.[49] A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh, and S. Pascoli, JCAP , 027 (2013), arXiv:1306.3219.[50] B. Li, A. Barreira, C. M. Baugh, W. A. Hellwing, K. Koyama, S. Pascoli, and G.-B. Zhao, JCAP , 012 (2013),arXiv:1308.3491 [astro-ph.CO].[51] W. Hu and I. Sawicki, Phys. Rev. D76 , 064004 (2007), arXiv:0705.1158 [astro-ph].[52] L. Reverberi and D. Daverio, JCAP , 035 (2019), arXiv:1905.07345 [astro-ph.CO].[53] A. Cardoso, K. Koyama, S. S. Seahra, and F. P. Silva, Phys. Rev. D , 083512 (2008), arXiv:0711.2563 [astro-ph].[54] A. de La Cruz-Dombriz, A. Dobado, and A. L. Maroto, Phys. Rev. D , 123515 (2008), arXiv:0802.2999.[55] J. Noller, F. von Braun-Bates, and P. G. Ferreira, Phys. Rev. D , 023521 (2014), arXiv:1310.3266 [astro-ph.CO].[56] H. Oyaizu, Phys. Rev. D , 123523 (2008), arXiv:0807.2449.[57] A. Barreira, B. Li, C. M. Baugh, and S. Pascoli, Phys. Rev. D86 , 124016 (2012), arXiv:1208.0600 [astro-ph.CO].[58] I. Sawicki and E. Bellini, Phys. Rev. D , 084061 (2015), arXiv:1503.06831.[59] S. Peirone, K. Koyama, L. Pogosian, M. Raveri, and A. Silvestri, Phys. Rev. D97 , 043519 (2018),arXiv:1712.00444 [astro-ph.CO].[60] D. Traykova, E. Bellini, and P. G. Ferreira, JCAP , 035 (2019), arXiv:1902.10687 [astro-ph.CO].[61] A. Silvestri, L. Pogosian, and R. V. Buniy, Phys. Rev.
D87 , 104015 (2013), arXiv:1302.1193 [astro-ph.CO].[62] S. Bose, W. A. Hellwing, and B. Li, JCAP , 034 (2015), arXiv:1411.6128 [astro-ph.CO].[63] F. Schmidt, Phys. Rev. D80 , 043001 (2009), arXiv:0905.0858 [astro-ph.CO].[64] H. A. Winther and P. G. Ferreira, Phys. Rev. D , 064005 (2015), arXiv:1505.03539 [gr-qc].[65] E. Krause and C. M. Hirata, Astronomy and Astrophysics , A28 (2010), arXiv:0910.3786 [astro-ph.CO].[66] T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys. Rev. D75 , 044017 (2007), arXiv:astro-ph/0607411 [astro-ph].[67] D. Comelli, M. Crisostomi, and L. Pilo, JHEP , 085 (2012), arXiv:1202.1986 [hep-th].[68] M. Lagos, E. Bellini, J. Noller, P. G. Ferreira, and T. Baker, JCAP , 021 (2018), arXiv:1711.09893 [gr-qc].[69] A. Avilez-Lopez, A. Padilla, P. M. Saffin, and C. Skordis, JCAP , 044 (2015), arXiv:1501.01985 [gr-qc].[70] C. Brans and R. H. Dicke, Phys. Rev. , 925 (1961).[71] N. Bolis, C. Skordis, D. B. Thomas, and T. Zlosnik, Phys. Rev. D99 , 084009 (2019), arXiv:1810.02725 [gr-qc].[72] D. B. Thomas, L. Whittaker, S. Camera, and M. L. Brown, Mon. Not. Roy. Astron. Soc. , 3131 (2017),arXiv:1612.01533 [astro-ph.CO].[73] R. Nagata, T. Chiba, and N. Sugiyama, Phys. Rev.
D66 , 103510 (2002), arXiv:astro-ph/0209140 [astro-ph].[74] T. Baker et al. , (2019), arXiv:1908.03430 [astro-ph.CO].[75] L. Lombriser, JCAP , 039 (2016), arXiv:1608.00522 [astro-ph.CO].[76] K. Koyama, A. Taruya, and T. Hiramatsu, Phys. Rev. D , 123512 (2009), arXiv:0902.0618 [astro-ph.CO].[77] J. Khoury and M. Wyman, Phys. Rev. D80 , 064023 (2009), arXiv:0903.1292 [astro-ph.CO]. [78] M. Milgrom, Astrophys. J. , 365 (1983).[79] M. Kopp, C. Skordis, D. B. Thomas, and S. Ilic, Phys. Rev. Lett.120