Scale-Dependent Gravitational Couplings in Parameterised Post-Newtonian Cosmology

Parameterised Post-Newtonian Cosmology (PPNC) is a theory-agnostic framework for testing gravity in cosmology, which connects gravitational physics on small and large scales in the Universe. It is a direct extension of the Parameterised Post-Newtonian (PPN) approach to testing gravity in isolated astrophysical systems, and therefore allows constraints on gravity from vastly different physical regimes to be compared and combined. We investigate the application of this framework to a class of example scalar-tensor theories of gravity in order to verify theoretical predictions, and to investigate for the first time the scale-dependence of the gravitational couplings that appear within its perturbation equations. In doing so, we evaluate the performance of some simple interpolating functions in the transition region between small and large cosmological scales, as well as the uncertainties that using such functions would introduce into the calculation of observables. We find that all theoretical predictions of the PPNC framework are verified to high accuracy in the relevant regimes, and that simple interpolating functions perform well (but not perfectly) between these regimes. This study is an important step towards being able to use the PPNC framework to analyse cosmological datasets, and to thereby test if/how the gravitational interaction has changed as the Universe has evolved.


Introduction
Tests of relativistic theories of gravity can be performed in physical environments that range from the Solar System [1,2], through to gravitational wave emission from binaries [3,4], and cosmological tests [5,6]. The motivation for these tests are varied, covering the basic scientific need to test physical theories in all areas in which they are to be applied, through to the desire to understand the constraints that observations can impose upon new propositions. On these grounds, one could argue that it is particularly important to test gravity in cosmology, due to the vast extrapolations that are required to apply gravitational theory on cosmological scales, and the relative abundance of new theories that appear in this field.
The wide array of spatial and temporal scales that need to be considered in cosmology, as well as the different physical environments that can exist, and the large number of theoretical proposals, also provide complications. These are factors that do not exist in the same way in other areas of gravitational physics, and which therefore present new challenges that need to be overcome. Historically, progress has been made in this field by either performing detailed predictions for a specific theory, before comparing to the data, or by constructing frameworks that cover multiple different theories or classes of theories. The latter of these approaches is more efficient, in the sense that fewer computations need to be done to constrain multiple different cases. It does, however, require us to be able to determine generic behaviour produced by different theories of gravity.
The gold standard in frameworks for testing gravity is the Parameterised Post Newtonian (PPN) formalism, which has seen great success in a variety of weak field systems [7]. However, it was not designed with cosmology in mind, and is not immediately applicable to observations made over cosmologically interesting scales. Work in producing frameworks for testing gravity in cosmology has therefore often tended to focus on creating new independent frameworks (see e.g. [5,8,9]). Such approaches usually rely on cosmological perturbation theory, and must therefore be extrapolated into the non-linear regime 1 . They also separate out the expansion of the "background" from other aspects of gravity, which overlooks a vital source of information, as well as obscuring the relationship with frameworks that are used in non-cosmological contexts (such as PPN).
Our approach to this problem is to take as our starting point the PPN framework, and assume this to be an appropriate description of gravitational physics on scales where Newtonian (and post-Newtonian) gravity are valid. We then construct the set of cosmological models that are consistent with this hypothesis, and which are statistically homogeneous and isotropic on large scales, by joining together many regions of space-time that are described in this way [13,14]. The result is a cosmology in which the background expansion and the weak gravitational field in the non-linear and mildly non-linear regime are consistent, and described by the same set of PPN parameters (suitably extended for the self-consistency of the cosmological model, and to include dark energy). As a corollary of this approach we also gain information about the behaviour of cosmological perturbations on the very largest scales [15,16]. We call the result 'Parameterised Post-Newtonian Cosmology' (PPNC).
The PPNC framework consistently parameterises gravity in the non-linear regime and at the level of the background, as well as on super-horizon scales. It is constructed from a set of time-dependent parameters that directly reduce to the PPN parameters in the appropriate limits, and that therefore allows the same aspects of the gravitational interaction to be constrained in both cosmological and non-cosmological systems. Furthermore, it provides information about both the small and large-scale limits of linear cosmological perturbations. It is therefore the natural extension of the PPN framework into the cosmological regime. In this paper we study the transition between the small and large-scale cosmological regimes of this approach using a pedagogical set of example theories, and how this transition should be expected to impact the growth of overdensities and influence observables.
This paper is laid out as follows. In section 2 we recap the mathematical formalism of the PPNC approach and the class of theories we will be considering. In section 3 we examine the behaviour of the perturbations and examine the quality of some simple interpolation functions. In section 4 we further examine the use of these interpolation functions to evolve the perturbations and calculate some simple observables. We conclude in section 5. We choose units such that c = 1, and use base 10 for the logarithms in all plots.

Mathematical formalism
Here we will present the salient details of the PPNC approach, which is based on an extended version of the PPN formalism suitably transformed for use in cosmology, before discussing a class of scalar-tensor theories of gravity as an example. For full details of the PPNC approach the reader is referred to Refs. [13][14][15][16]. For the scalar-tensor theories of gravity the reader may refer to, for example, Ref. [5].

Parameterised Post-Newtonian Cosmology
The geometry of the constructed space-time in the PPNC approach can be written as where Φ 1 and Ψ 1 are both scalar functions of τ and x i , and a(τ ) is the scale factor as a function of conformal time. The reader will note that vector and tensor gravitational fields have been neglected at leading order, which is consistent with expectations from cosmological perturbation theory, as well as the order at which these objects appear in the post-Newtonian expansion that is used to generate this metric. We have also opted to consider geometries in which the background in spatially flat; this is motivated by CMB observations, but still allows the inclusion of small amounts of spatial curvature through the perturbation Ψ.
The metric presented in Eq. (2.1) is written in the longitudinal (conformal Newtonian) gauge, which has the benefit of being well defined in both the cosmological perturbation theory and post-Newtonian sectors of the expansions that we require [17]. It also corresponds to the leading-order part of the 'post-Newtonian gauge', in which the PPN approach is typically formulated. We note that the metric in Eq. (2.1) is not a starting ansatz in this approach, but is constructed from small regions of space-time described by the PPN "test metric", which is itself a post-Newtonian expansion about Minkowski space. In this sense, the behaviour of the scale factor a(τ ) and the potentials Φ and Ψ are emergent quantities, related directly to the gravitational potentials that occur in the test metric. We expect this geometry to be applicable to all astrophysical objects, including everything down to the scale of neutron stars and black holes, as all such bodies are expected to be well modelled by Newtonian (and post-Newtonian) gravity.
The gravitational parameters required for this construction, at the stated level of accuracy, are as follows: {α, γ, α c , γ c } , where each of these should be understood to be a function of conformal time τ , but not spatial coordinates x i . The first of these quantities appears in the Newtonian limit of the field and geodesic equations in exactly the same place as Newton's constant G. In the context of weak-field gravity in non-cosmological systems, it is therefore routinely set equal to one. In the cosmological context, we can do this at the present time τ 0 , but this does not necessarily mean that α will retain this value at different points in our cosmic history. It must therefore be a function of time, with the boundary condition α(τ 0 ) = 1 at the present time τ 0 . The parameter γ(τ ) is the 'curvature of space' parameter, which is tightly constrained by solar system observations to lie within one part in 10 5 of its GR value of one [18] at the present time (but not necessarily in the past). Finally, α c and γ c are cosmological parameters, which are required for consistency of the relevant equation when cosmological evolution is permitted, and which must obey the constraint whereρ is the average cosmological mass density, and the hats denote differentiation with respect to the number of e-foldings, ln a. These two parameters include information about the dark energy (which here is not taken to be a part ofρ, the mass density of baryons and dark matter).
The Friedmann equations that emerge from averaging the small-scale geometry take the following form: where H = a /a is the conformal Hubble rate, and primes denote differentiation with respect to τ . The parameter which determines the magnitude of the effective Newton's constant α can be seen to enter into the second Friedmann equation, multiplying the average mass density. Conversely, it is the curvature of space parameter γ which enters into this position in the first Friedmann equation. It can be seen that Eq. (2.3) acts as an integrability condition on these equations, and that α c and γ c take the place of dark energy terms. The equations governing the perturbations Φ and Ψ are given by [14] −H 2 Φ − HΨ + 1 3 and [16] where ∇ 2 denotes the spatial Laplacian constructed from partial derivatives, and where µ and ξ are generalisations of Newton's constant 2 . The density contrast is written in these equations as δ = δρ/ρ, and the 3-velocity it written as v i 3 . The G is a further function, which determines the form of the momentum constraint equation (2.8) [16].
The PPNC approach tells us the small-scale limit that the effective Newton's constant and momentum constraint parameters must take is given by while adiabatic perturbations on very large scales require (2.12) The PPNC approach therefore gives us a remarkable set of equations, parameterised by four functions of time (2.2) subject to one constraint (2.3), from which we can write both the background equations (2.4)-(2.5) and the small and large-scale limits of the perturbation equations (2.6)-(2.8), including the fully non-linear regime.
The equations presented above are expected to represent the cosmological behaviour of any theory of gravity that fits into the PPN framework, in terms of direct generalizations of the PPN parameters. What remains to be understood is how the functions µ, ζ and G interpolate between the large and small-scale limits presented above. In this study we will examine how some simple theories behave in this transition regime, and compare their behaviour to the simplest elementary functions. Knowledge of how these functions behave as they transition between the large and small-scale limits is important for understanding the formalism, and may be significant for the calculation of some cosmological observables, as well as direct integration of the parameterised equations.
The PPNC framework does not require density contrasts to be small. However, in this work we are interested in studying the transition regime between large and small scales, which is within the regime we expect to well approximated by linear perturbation theory. This is because the small-scale end of the linear perturbation regime overlaps with the regime in which post-Newtonian expansions are valid. As all non-linear cosmological structure formation happens within the post-Newtonian regime, the transition away from small-scales must occur within the perturbative regime. As a result, for the remainder of this work we can safely expand the mass density as ρ =ρ + δρ so that the gravitational physics is tractable without having to perform N-body simulations (see Refs. [10,11] for a discussion of this in the context of modified gravity). The small-scale end of the range of k-values that we will consider (with L 100 Mpc) is then assumed to be well-approximated by post-Newtonian expansions so that the results in Eq. (2.9) can be applied [21], while on the largest scales we can use the limits given in Eqs. (2.10)- (2.12). Linearising in all perturbations then allows us to write the Fourier space versions of Eqs. (2.6)-(2.8) as where we have followed the usual convention of using the same symbols for quantities defined in Fourier space and real space to avoid further cluttering our notation. It is the scale dependence of the Fourier space couplings µ, ξ, and G that will be investigated below 4 .

Conceptual comparison to other approaches
Before continuing, we note that the PPNC approach is different in both concept and execution to the Effective Field Theory (EFT) approaches to modified gravity in cosmology that have recently found prominence in the literature [22][23][24].
One key difference is that PPNC does not require us to specify the field content of the theory, but instead relies on the specification of the possible relationships between gravitational potentials and matter fields, which are often in the form of Coulomb potentials (or generalisations of Coulomb potentials 5 ). This approach is borrowed from the PPN formalism upon which PPNC is constructed, and which has proven so successful in theory-independent tests of gravity in the Solar System. Conversely, in EFT approaches one specifies the field content from the outset, which is often taken to be a single additional scalar field. A second key difference is that EFTs are constructed using cosmological perturbation theory, in which background and perturbation equations are treated separately. On the other hand, the PPNC approach has its foundation in post-Newtonian expansions, from which one can construct both cosmological background and perturbations [31].
This differences in starting point have some important consequences. For one, the structure of the EFT approach means that functions that are introduced in parameterisations will appear in either the background or the perturbation equations, but not both (as they are treated separately). This is in contrast to the PPNC approach where parameterizing functions are introduced into the post-Newtonian sector initially, and hence simultaneously appear in both the cosmological background and perturbation equations that emerge. Of course, for any specific theory that fits into both approaches the EFT background and perturbation functions would be found to depend on the underlying parameters of the theory in such a way that they could not be chosen independently. The PPNC approach makes this interdependence explicit, and as a consequence does not (in its current formulation) appear to be consistent with designer theories in which cosmological background and perturbations can be modified independently.
The fact that the PPNC approach is not built explicitly on either cosmological perturbation theory or an EFT expansion has some further advantages. The first is that it allows the time-dependence of the PPN parameters to be constrained with cosmology, and therefore for the constraints that are imposed in the Solar System to be directly extended into a much wider context [7]. The second is that the equations governing the growth and behaviour of inhomogeneities do not require the density contrast to be small, which means that they can provide a complete description of cosmological structure formation on all scales, such that the non-linear regime does not have to be handled in a different or separate manner. This could prove to be an important benefit for upcoming surveys such as Euclid, and we leave to future work the creation of N-body (or COLA [25]) simulations within this framework. Lastly, the fact that the PPNC formalism does not require the specification of the field content of the theory from the outset gives it the potential to have considerably more theory independence.

Scalar-Tensor Theories of Gravity
The example theories we will use to investigate the transition regime from small to large scales will be the Bergmann-Wagoner class of canonical scalar-tensor theories [26,27]. These have the action where ω = ω(φ) and Λ = Λ(φ) are functions that specify the coupling between the scalar and tensor degrees of freedom in the theory, and the scalar field potential, respectively. This class of theories contains within it the very well studied Brans-Dicke theory (ω =constant and Λ = 0) [28], as well as GR with a cosmological constant (ω → ∞ and Λ =constant). These are fully conservative theories, in that they do not have any preferred frame or preferred location effects, while they possess the full set of conservation laws (energy, momentum angular momentum and centre-of-mass motion). They are the sub-class of Horndeski theories with G 2 = G 2 (φ) and G 4 = G 4 (φ), and with all other G i = 0 [29].
In the perturbed Robertson-Walker geometry given by Eq. (2.1), the Friedmann equation in these theories can be written as and the Klein-Gordon equation for the scalar field as whereρ is again the mass density and radiation has been neglected. The second Friedmann equation follows from differentiating Eq. (2.17), and eliminating the second derivative of the scalar field using Eq. (2.18).
The perturbation equations around a Robertson-Walker background can be split up into constraint and evolution equations; we present these equations up to linear order in Fourier space. The Hamiltonian and momentum constraints can be written as and where v is the velocity potential. We have also abused notation to write φ and δφ as the background and perturbation to the scalar field. As is usual in longitudinal gauge, the shear evolution equation also reduces to a constraint, which in this case can be written in the particularly useful form The evolution equations for the matter variables can now be written as while the perturbed Klein-Gordon equation gives (2.24) and the Raychaudhuri equation gives This provides the complete set of linearized scalar equations in these theories.
The relevant PPN parameters α and γ are obtained by taking a post-Newtonian expansion about Minkowski space and comparing the resultant geometry to the test metric (see Ref. [1] for details). The result is Usual practise, within the literature on weak-field gravity, would be to set φ (by a choice of units) so that α = 1. Such a choice means that Newton's constant is recovered in the proper place in the Newton-Poisson equation and the geodesic equation. The appropriate choice in cosmology is to make this choice at τ = τ 0 , so that Newton's constant takes the value G at the present time, but is allowed to take different values at other points in cosmic history. This naturally incorporates the idea of a "varying Newton's constant" [30].
In addition to the usual PPN parameters, we are required to introduce the following additional two parameters in order to construct the cosmological equations [31]: One can verify that these parameters satisfy the constraint equation (2.3). In the limits ω → ∞ and Λ → constant, we can see that these equations reduce to α c = −2γ c = Λ and α = γ = 1, thereby recovering GR with a cosmological constant, as expected 6 .
In what follows we will take ω and Λ to be constants (unless otherwise specified), and integrate the background and perturbation equations of these theories directly in order to verify the theoretical predictions of the large and small-scale limits that are given by the PPNC formalism. We will further use our numerical solutions to investigate how the effective Newton's constants µ and ξ vary between these two limits within this class of theories. The results will then be used to determine how well a simple interpolating function would have performed in the PPNC equations, when they themselves are directly integrated. We will use the resulting density contrasts, and some simple observables, to determine the accuracy that the PPNC equations would have achieved in comparison to the direct integration of the example theories outlined above. Due to theories that introduce new scales not being contained in the set currently covered by the PPNC approach (see comments in Section 2.1.1), we expect this transition behaviour to be broadly representative of the general case.

Application in Example Theories
In this section we will solve the differential equations that describe the background cosmology and all relevant perturbed quantities in the example theories described above. We first solve for the background cosmology 7 to redshifts larger than 1 100, and verify that the evolution we obtain obeys the background PPNC equations (2.4)-(2.5). We then solve for the perturbations, using suitable initial conditions, by evolving Ψ and φ using the Raychaudhuri and Klein-Gordon equations, while calculating Φ, δρ and v at each time step using the anisotropic stress, Hamiltonian and momentum constraint equations. We work in code units where c = 1, and such that the conformal time today is unity, τ 0 = 1.

Background Cosmology and Parameter Evolution
The background scalar-tensor equations (2.17)-(2.18) can be integrated once values for ω and Λ have been provided, along with suitable boundary conditions for a(τ ) and φ(τ ). For illustrative purposes, the integrations we present in this section will be for several different combinations of constant values of ω and Λ. We will not concern ourselves with whether or not the values we choose are compatible with observations (some of them will not be), as at this stage we are concerned only with exemplifying how the PPNC formalism works in example theories that differ from GR, and not with modelling any particular physical phenomena. For these purposes it is useful to consider theories that are allowed to differ substantially from GR. The boundary conditions we will use are that a(τ 0 ) = 1 and that φ(τ 0 ) = (4 + 2ω)/(3 + 2ω) in the attractor solution for φ(τ ). This last condition is chosen so that α(τ 0 ) = 1 from Eq. (2.26). Both choices can be thought of as specifications of units.
The behaviour of a(τ ) and φ(τ ) have been well studied in the literature on these theories 8 , and we will not present them here beyond commenting that φ increases with time for positive values of ω. Once these solutions are known, however, they can be used to derive values for the PPNC parameters specified in Eqs. (2.26)-(2.28). In turn, these can be used to derive the values of the parameters specified in Eqs. (2.9)-(2.12). We plot the large and small-scale values of µ(τ ) and ξ(τ ) in Figures 1 and 2, for several different choices of ω and Ω Λ . The reader will note that the small-scale values of these two quantities (i.e. the solid lines) correspond directly to the values of the PPN parameters γ and α at different times in cosmic history. We consider how such time dependence could be fitted in Appendix A.
In both Figures 1 and 2 we can see that at early times the behaviour of µ and ξ is very similar on both large and small spatial scales. This is because the γ and α terms dominate the expressions for these quantities on large scales. In fact, it is only in the presence of a non-zero Λ that the large and small-scale limits are different. This appears to be a feature of the particular theories we are considering, and is not expected to be a generic property of theories of gravity in general. The curves that depict the small-scale limit of ξ in Figure 2 can be seen to converge at τ = τ 0 (today), which is due to our choice of boundary condition on φ. This also explains why the solid curves in the right-hand plot in Figure 1 converge today; as the value of γ depends only on ω, once φ 0 has been fixed.
Going back in time, α and γ increase for all of our theories. The differences between small and large scales, however, increase into the future. Increasing the value of Ω Λ can be seen to exacerbate this situation, while increasing ω has the opposite effect. In general, larger values of ω can be seen to result in a slower evolution in the values of µ and ξ, as expected (as this is the limit in which GR is recovered). While α is always greater than the GR value of unity, it can be seen that γ crosses the GR value of unity during its evolution (starting above it and finishing below it, in the cases we consider). It will be the behaviour between the large and small-scale limits depicted in these plots that will concern us for much of the rest of this study. For this we will also need to understand the evolution of perturbations.

Perturbations and Scale-Dependent Couplings
In this section we will evolve the perturbations in our example theories, and use them to determine the scale dependence of the couplings that appear in the PPNC parameterisation. This will help to provide guidance of how the interpolation between different spatial scales works within a mathematically well-defined alternative to GR.

Initial Conditions
To begin, let us specify the initital conditions for our numerical integration. We wish to specify these as adiabatic fluctuations that satisfy [33] δφ φ = δρ ρ . (3.1) We will set initial conditions for our perturbations by starting the evolution at an early enough time that all Fourier modes of interest start off outside the horizon. Isocurvature modes are neglected, as they are generally expected to be small on these scales [34]. Under the adiabatic condition in Eq. (3.1), it can be seen that the large-scale limit of Eqs. (2.19)-(2.25) admit the solution with the momentum constraint giving v = H −1 ζ and where is the comoving curvature perturbation on uniform density hypersurfaces. The c 1 and c 3 in these expressions are constants, and taking ζ(τ i ) = c 2 to be another constant at the initial time τ i allows us to write initial conditions for δφ and δφ as and the initial conditions for Ψ and Ψ as Together with the constraint equations, these conditions are sufficient to integrate the Raychaudhuri and Klein-Gordon equations, and therefore serve as initial conditions for superhorizon adiabatic perturbations. The two constants c 1 and c 2 in the equations above can be thought of as specifying the initial amplitude of the growing and decaying modes within this setup. To ensure that the decaying mode is under control, we impose on these constants the supplementary condition 9 δφ | τ i = 0, which leaves us with a single constant that specifies the amplitude of our perturbations at the initial time τ i .

Present-Day Values of Gravitational Couplings
In this section we examine the gravitational couplings µ, ξ and G at the present time τ = τ 0 , by numerically evolving our perturbation equations from the initial conditions specified above. This is achieved by taking the values of the Φ, Ψ and δρ from our integration, and substituting them into the following re-arranged versions of Eqs. (2.6)-(2.8), obtained after linearisation and Fourier transformation: The results, plotted in Figures 3 to 5, show the values of the coupling functions µ, ξ and G at z = 0 in these theories. In each of these figures the wavenumber k of the perturbative modes is plotted on the xaxis in units of the Hubble rate, H. We can clearly see the scale-independent limits predicted by Eqs. (2.9)-(2.12) are being approached as k H on the right of each plot, and k H on the left. Each line in each plot represents a particular combination of constants in the theory, but in each case we can see a smooth transition between these two limits. The vertical black lines in each plots correspond to the Hubble scale k = H, which can be seen to be the scale about which the transition takes place. In each figure the left-hand plot shows the effect of changing the value of the coupling parameter ω, while the right-hand plot shows the effect of changing the amount of dark energy Ω Λ .    In Figures 3 to 5 we have included simple interpolation functions, to explore how well they fit through the transition scale. We considered several possibilities for interpolating between the large and small-scale limits, but found the simplest to be a tanh function 10 , as suggested in [15]: wheref = (S + L)/2 is the mean of the small (S) and large (L) scale limits of the coupling, and A f = (S − L)/2 is its amplitude. Using the definition of the tanh function, in terms of exponentials, allows us to write this as which is similar to previous attempts to parameterise the scale dependence of gravitational couplings using Padé approximants [35,36]. This interpolation function is fixed to reach a value halfway between the two endpoints at the Hubble scale, but otherwise to require no additional parameters (i.e. they are a 'zero-parameter' fit). While they give a reasonable first approximation, these simple functions clearly miss the detailed (oscillating) features. We will examine these features in more detail in later sections, but note that they appear to be similar to the Gravity Acoustic Oscillations (GAO) discussed in Ref. [37]. In all cases, we can see that deviations from GR reduce to zero in the limit ω → ∞ (as expected), and that Ω Λ acts to increase the difference between the values of couplings at small and large scales (with increasing values of Ω Λ causing a bigger difference). For the values of ξ plotted in Figure 4 we can see that all curves approach unity as k → ∞, as required from the present-day value of α and Eq. (2.9). Correspondingly, the values of µ plotted in Figure 3 approach the values of γ given in Eq. (2.26), again as anticipated from Eq. (2.9). Of course, when using this formalism to constrain gravity using cosmological observations one may wish to use observations performed in nearby astrophysical systems to constrain the present day value of γ, which requires it to be very close to unity [18].
The final quantity we examine over this transition region is the gravitational 'slip', which is the ratio of the two gravitational potentials η = Φ/Ψ, and is shown in Figure 6. The value this quantity should take on small scales can be readily calculated from the PPN approach, from which we obtain η = α/γ. This value can be seen to be obtained from the large-k limit of our numerical evolution of the perturbations. On large scales, however, it is so far not clear how to obtain an analytic expression for η (such an expression was missing from section 2.1). At our present level of understanding, we therefore take the value of the slip on large scales to be an additional parameter of the framework, which can be discarded if and when an analytic expression in terms of other PPNC parameters is obtained. We can see from Figure 6 that the small-scale value of η varies with ω, but not Ω Λ (as expected). On the other hand, the large-scale limit varies with both parameters, and is less than unity for all the parameter values we have examined within our class of example theories. Interestingly, the slip is the only one of our coupling functions that has quantitatively different large and small-scale limits even in the matter-dominated era 11 . As expected, it can be seen that varying Ω Λ only makes a difference at late times, where it acts to increase (decrease) the large-scale value of µ (ξ) compared to the small scales. The effect of changing ω is a little more complicated. During the matter-dominated era the large-scale behaviour of µ and ξ is dominated by parameters γ and α, so that there is no difference between the behaviour on large and small scales (note that interesting features around the horizon scale still exist). During Λ domination, however, we see that decreasing ω results in increasing deviations from GR 12 .

Evolution of Perturbations and Couplings
As the universe transitions to a Λ-dominated phase, we can see that the value of µ becomes close (but not identical) to its GR value on large scales regardless of the value of ω, while on small scales decreasing ω still results in substantial deviations from the GR value of unity. This behaviour is reversed for ξ: the large-scale behaviour can be seen to deviate substantially from the GR value, while on small scales it approaches one (exactly, by construction, as described above). We note for the reader that the bottom two panels of Figures 7 and 8 correspond precisely to the left-hand plots of Figures 3 and 4, reproduced here for ease of comparison.
In the interests of brevity, we have not plotted the time dependence of the slip η or momentum constraint parameter G here. Both are relatively simple. The large-scale value of G is zero at all times, and the small-scale value varies monotonically from being close to zero in the matter-dominated era, to around −0.05 today (for ω = 10 and Ω Λ = 0.7). The oscillations close to the horizon scale shown in Figure 5 persist at a similar amplitude over time, but are shifted vertically as the small-scale value of G evolves. The small-scale value of η depends only on ω, with no time dependence (because the time dependence in α and γ is the same, therefore cancelling-out when the ratio is taken). The large-scale value of η decreases monotonically from about 0.95 deep in the matter-dominated era to around 0.9 today (for ω = 10 and Ω Λ = 0.7). Again, the oscillations around the Hubble scale do not change significantly over time.

Accuracy of Interpolating Functions
Let us now examine the performance of the interpolation function from Eq. (3.9), in reproducing the coupling functions from our set of example theories. For this, we will use as a specific example theory given by ω = 10 and Ω Λ = 0.7 (i.e. the common case in Figs. 7 and 8). We specify the fractional difference between the simple interpolating function and the true value for this theory using the notation where f here is being used to denote any one of the gravitational couplings (i.e. µ, ξ or G). Figure 9 shows the accuracy of the interpolation as a function of time and scale for µ and ξ, as quantified by this fractional error. The range of times being shown in this figure extends from deep in the matter-dominated era through to late times in the Λ-dominated era. In creating these plots we have interpolated between gridpoints, avoiding points where our solver produced large numerical errors. It can be seen from Figure 9 that the simple interpolation works well for both µ and ξ on most spatial scales (with errors below ∼ 1%), but with the oscillations around the Hubble scale leading to larger errors (up to ∼ 5%). The detail of these plots shows some interesting structure, if we look at how the error evolves with time in each of the two cases. For the coupling µ the situation is relatively simple: there is primarily a single peak that does not change substantially with time, and which tracks the position of the Hubble scale. The coupling ξ is more complicated: in this case there are several peaks that move position relative to the Hubble scale. In a cosmology with Λ = 0 this more complicated behaviour does not occur, and one recovers a plot that is qualitatively similar to the case involving µ.

Evolution using Parameterised Equations
Part of the purpose of having a parameterised framework for gravity in cosmology is to be able to model gravitational physics without having to specify a theory or class of theories a priori. Practical implementation of such a framework therefore requires us to be able to evolve perturbations using only the parameterised equations, which is what we will consider in this section. We will perform this evolution using the zero-parameter interpolating function from Eq. (3.9), which we will call the 'simple parameterised' model, and evaluate the accuracy of the results by comparing to our class of example theories only after the evolution has been calculated. The quantities we will consider when using this approach will be the density contrast, as well as some example observables.
With regard to implementation, we will evolve the density perturbation and the velocity using the conservation equations for these quantities, which do not depend on the theory of gravity for their functional form (assuming the underlying theory admits mass and momentum conservation at leading order). The parameterised Hamiltonian constraint equation (2.6), together with the slip η, will then be used to evolve the leading-order scalar metric potentials. For the functional form of the parameters {α, γ, α c , γ c } we will use the evolution determined by the background, as discussed in section 3.1. For the slip, the small-scale limit will be calculated from the parameters α and γ, while the large-scale limit will be extracted directly from the evolution of our example theory 13 . We start the evolution of our perturbations at redshift z = 1 100, and end at the present day τ 0 , so that we cover both the matter and Λ-dominated eras of cosmic history. We use the same initial conditions for evolving the parameterised equations as in our example theories, for ease of direct comparison.

Density Contrasts
Calculating the density contrast at τ 0 using the parameterised equations, and the method described above, leads to results that are in good agreement with the results we obtain directly from our example theories across a range of spatial scales. There are, however, errors introduced due to the inadequacy of the simple interpolating function (3.9) to model all of the detail around the horizon scale (as was found in section 3.3). The accuracy of using the simple interpolation is displayed in the plots in Figures 10, where we show the ratio of the density contrast computed by evolving the parameterised equations with the simple interpolation to the value obtained from our example theories directly. The inaccuracy of the interpolation can be seen to create an error of order ∼ 1% around the horizon scale.
Decreasing ω (and therefore increasing the deviation from GR) increases this error, whereas it appears to be mostly insensitive to the value of Ω Λ . These errors are smaller than those in µ and ξ, as displayed in Figure 9.
As well as considering the fractional error introduced from evolving the parameterised equations with the simple interpolation, it is also of interest to consider the relative size of this error compared to the difference in density contrast that occurs due to changing the underlying theory of gravity (which is, after all, the signal that we would ultimately hope to detect). These two things are displayed at τ 0 , and over a range of spatial scales, in the left-hand plot of Figure 11. Here we have chosen to plot the fractional error in δρ for an example theory with ω = 10 and Ω Λ = 0.7 (the blue curve in each of the panels in Figure  10), with the change in density contrast that one gets from increasing the value of ω to 100 and 10 6 (keeping Ω Λ = 0.7 for all curves). It can be seen that the effects from changing the value of ω is much larger than the error introduced from the imperfect interpolating function, suggesting that the bias on constraints on the parameters will be small compared to any detected deviation from General Relativity. Of course this bias should still be understood or modelled, which would presumably be at the expense of adding further parameters to the framework.  : Left: Accuracy of the density contrast at z = 0 from evolving the simple parameterised equations (blue), compared to the difference induced by changing ω to 100 (orange) and 10 6 (green). Right: Accuracy of the density contrast today from evolving the simple parameterised equations (blue), compared to the error introduced by assuming a ΛCDM background (red for ω = 100; purple for ω = 10) or ΛCDM perturbation equations (orange for ω = 10; green for ω = 100). In each plot we take Ω Λ = 0.7 and ω = 10, unless otherwise stated.
In the right-hand plot of Figure 11, we compare the inaccuracies introduced from evolving the simple parameterised equations with the inaccuracies induced by applying the modified gravity parameterisation to either the background only or the perturbations only (instead of both). The inaccuracy introduced by only modifying the background or the perturbations is much larger than that introduced by evolving the simple parameterised equations directly (this is true even for theories that are much closer to GR). Interestingly, the inaccuracy introduced by not including the correct background expansion appears to be even larger than that from not evolving the perturbations with the correct equations. This result is significant for the PPNC approach because the PPNC approach provides a way to consistently parameterise both the background and the perturbations in terms of the same underlying parameters. This is not true of any other attempts at creating a theory-independent framework for testing gravity in cosmology of which we are aware. These results suggest that when constraining modified gravity in a model independent or phenomenological way, the background should generally be allowed to deviate from ΛCDM (for example using a dark energy equation of state w(a)) in order to prevent biasing the constraints.

Lensing of Light
We now wish to consider the effect of evolving the perturbations using the simple parameterised equations on some cosmological observables. This is not intended to be a detailed observability study, nor is it carried out in detail for specific experimental setups. Instead we wish to make a preliminary examination of how well the parameterised equations fare (in comparison to the "true" values), when calculating observables. For the lensing of light we simply have to integrate (Φ + Ψ) 2 over the comoving distance χ, from χ = 0 out to some source χ * . Doing this, we see that the parameterised equations reproduce the numerical values that one would get using the specific equations of any of our example theories with errors that are entirely negligible. This is expected, as lensing is primarily operating well within the horizon, and within this regime the parameterised equations give very good results. We also note that we reproduce the same results if we include the lensing function (1 − χ/χ * ) /χ 4 in the integral over χ. This simple approach to weak-lensing calculations will of course break down on the largest angular scales; we leave a detailed examination of the effect of using the parameterised equations in this regime for future work.

Integrated Sachs-Wolfe
As the Integrated Sachs-Wolfe (ISW) effect primarily operates over larger scales, we expect that this observable will be more sensitive to the inaccuracy of our simple interpolation between small and large scales. To investigate this, for a given , we look at the integral giving the contribution of each k-mode to that mode: where j l is the Bessel function. We calculated this quantity with both the simple parameterised equations as well as with the correct evolution for specific example theories at = 10 and = 100. Our results are shown in the left-hand panel of Figure 12. As expected, there is a noticeable difference between the curves for some k values. We can also integrate our results over k, with the usual factor of k −1 to account for a scale-invariant primordial spectrum, in order to examine the accuracy of the ISW C s. This is shown by the blue curve in the right-hand plot of Figure 12. The error peaks around a couple of percent, similar to (but a little worse than) the inaccuracy we found in the density contrast. We note that these errors are never greater than 20% of cosmic variance over this range of . Finally, we again consider the consequence of changing ω, and show the results in Figure 12. Again, each curve is plotted as a ratio of the results obtained in a theory with ω = 10, and again we see that the errors introduced by integrating the simple parameterised equations directly is substantially smaller than the change in the observable due to changing the value of ω. This indicates that the errors introduced when integrating the simple parameterised equations directly is not likely to strongly bias any parameter constraints based on this observable. Figure 12: Left: The contribution of different k modes to the ISW effect, for = 10 and 100. For = 10, the value obtained from the theory with ω = 10 is blue, and the value from evolving with the simple parameterised equations is orange. For = 100, the theory with ω = 10 is green, and the parameterised version is red. Each pair of curves has been normalised to the highest point, to allow them to be shown on the same axes. In all cases Ω Λ = 0.7. Right: Errors induced in the ISW observable at τ 0 from using the simple parameterised equations with ω = 10, compared to different theories. Each curve is the ISW observable normalized by its value in the fiducial ω = 10 theory. Blue is obtained by evolving the parameterised equations, orange is ω = 100, and green is ω = 10 6 .

Discussion
We have investigated the scale-dependence of gravitational couplings in the PPNC approach to cosmological modelling, in which the coupling functions that appear in the cosmological perturbation equations can be directly linked to the parameters of the PPN framework (the workhorse of astrophysical tests of gravity). As well as just describing perturbations, however, the PPNC approach also specifies the evolution of the cosmological background in terms of the same set of underlying parameters. This link between gravitational physics on small and large cosmological scales has allowed the limits of the relevant gravitational couplings to be deduced, but (before the present work) without any information about how to interpolate between these different regimes.
Our study has included an investigation of the application of the PPNC framework within the class of canonical scalar-tensor theories of gravity given in Eq. (2.16). Within these example theories we have numerically determined the evolution of the PPNC parameters {α, γ, α c , γ c }, and used their values to determine the large and small-scale limits of the coupling parameters that govern the evolution of perturbations. These results agree perfectly with existing theoretical predictions in every case we considered, and therefore verify the ability of this framework to be applied to the cosmological background, the small-scale quasistatic linear and quasi-non-linear regime, as well as the large-scale super-horizon regime.
We have investigated the way in which the gravitational coupling parameters from the perturbation equations transition from small to large cosmological scales. This is the first time that this transition region has been investigated in the PPNC approach. We have found that, within our class of example theories, there is an underlying smooth transition between the small and large scales, which exists together with oscillations in these couplings around the Hubble scale (due to the existence of the scalar field). The smooth transition is reasonably well modelled by a simple tanh function, as proposed in Ref. [15], and the oscillations would appear to be linked to the Gravity Acoustic Oscillations (GAOs) from Ref. [37].
We have considered the inaccuracy introduced into some basic observational predictions by using the simple tanh function. For flat-sky lensing calculations we find this inaccuracy to be small, as expected since the modes that contribute most are well below the horizon scale thus outside of the transition regime. This is not the case for the ISW effect, where we have shown (see Section 4.3) that using the simple tanh function can introduce errors at the level of about a percent. However, this inaccuracy is significantly below the level of uncertainty in the ISW effect caused by cosmic variance, as well as the effect of modifying the theory of gravity itself. Nonetheless this inaccuracy represents a potential source of uncertainty in applying this formalism to interpret real observations, and may well also be larger for other observables we have not considered here. It should therefore be recognised and understood, or itself modelled (which would presumably be at the expense of adding further parameters to the framework). We find that these conclusions remain qualitatively unchanged when ω and Λ are made smooth functions of the scalar field φ.
Our study lays the groundwork for sensibly implementing the PPNC approach in Boltzmann codes and N-body simulations, and thus constraining gravity using all cosmological scales, and in a way that can be combined with gravitational constraints from other regimes and scales (due to its explicit links with PPN). Of course, much work remains to be done in working through the practicalities of such an application, including how to model unspecified functions of time, and whether methods such as binning or reconstruction from observations would be possible or appropriate. One might also consider applying this approach to constraining vector gravitational fields in cosmology, which could be achieved using results and methods in Refs [16] and [38][39][40][41], or higher-order gravitational fields [42][43][44]. Investigations of these questions, as well as remaining theoretical issues, will be the subjects of future work.

Appendices
A Time-Dependence of Parameters Figures 1 and 2 showed the behaviour of the PPN parameters α and γ as functions of conformal time τ , within some of our example theories of gravity. In general, however, the time dependence of these parameters is not known, and indeed is likely to vary considerably between different theories of gravity. It is therefore useful to pause and consider how well we might have reproduced the calculated evolutions of α and γ if we had instead posited them as simple functions, with degrees of freedom to be fixed. This situation is illustrated in Figures 13 and 14, where we again reproduce the behaviour of γ and α from our example theories with different values of ω and Ω Λ , but this time also with a simple fitting function f (t) = At −n +B where A, B and n are constants. Fixing the end points, as would be achieved by (for example) using astrophysical data of gravitational experiments at the present time, and adjusting the remaining freedom we can see it is possible to achieve reasonable fits. In reality, of course, we would not fit functions to pre-specified curves (as these would not be known), but instead would use data to constrain constants.  For completeness, let us also present the behaviour of the parameters γ c , α c and the large-scale limit of η = Φ/Ψ as functions of τ in Figure 15, for several different choices of ω and Ω Λ . The reader may note that we have chosen to plot the ratio of γ c and α c to the background matter densityρ(τ ), in order to make them dimensionless. This choice also more accurately reflects the importance of these parameters at different stages in the cosmological evolution, and is thus more comparable with Figures 1 and 2. The signs of α c and γ c are determined by the signs of the Λ term in Eq. (2.28) and (2.27), as this term dominates when these parameters are cosmologically relevant. Increasing Ω Λ results in the amplitude of these parameters increasing, while decreased ω (making the theory further from GR) decreases their amplitude. The behaviour of the large-scale limit of η can be seen to have a different behaviour to the other quantities. In all cases, the large-scale slip η is below the GR value, and becomes further from GR with decreasing ω (with this effect increasing at late times).
Increasing Ω Λ reduces the value of η compared to its value in GR at late times.

B Other Interpolating Functions
One possible generalization of the zero-parameter interpolation function introduced in Eq. (3.9) is to expand it to include a single additional constant A, as The fit of the interpolation can be improved by tailoring the value of A, as shown in Figure  16. However, the optimal A will certainly vary between different theories of gravity, so if we wish to remain as model-independent as possible then the zero-parameter interpolation seems most appealing. Should a more complicated interpolation function be required then the parameter A (or its equivalent in other possible generalizations) would need to be included in the data analysis as a nuisance parameter, and marginalised over. Whether or not this would be beneficial would require a study of relevant statistics, such as Bayesian information.