Testing nonlocal gravity with Lunar Laser Ranging

TL;DR

This paper evaluates various nonlocal gravity models against Lunar Laser Ranging constraints on the variability of G, confirming the exclusion of RR and DW models while supporting the viability of the RT model.

Contribution

The study provides a detailed analysis of nonlocal gravity models in light of Lunar Laser Ranging data, confirming the exclusion of RR and DW models and supporting RT model viability.

Findings

RR model is ruled out by LLR constraints.

DW model is ruled out due to ineffective screening.

RT model remains consistent with phenomenological tests.

Abstract

We study the impact of the limit on $|\dot{G}|/G$ from Lunar Laser Ranging on "nonlocal gravity", i.e. on models of the quantum effective action of gravity that include nonlocal terms relevant in the infrared, such as the "RR" and "RT" models proposed by our group, and the Deser-Woodard (DW) model. We elaborate on the analysis of Barreira et al. [1] and we confirm their findings that (under plausible assumptions such as the absence of strong backreaction from non-linear structures), the RR model is ruled out. We also show that the mechanism of "perfect screening for free" suggested for the DW model actually does not work and the DW model is also ruled out. In contrast, the RT model passes all phenomenological consistency tests and is still a viable candidate.

Figures

Equations

$$\Gamma =\frac{m_{\mathrm{Pl}}^2}{2}\int d^{4}x\sqrt{-g}\left[R-g^{2}R{\square }^{-1}R\right],$$

$${\Gamma }_{\mathrm{DW}}=\frac{m_{\mathrm{Pl}}^2}{2}\int d^{4}x\sqrt{-g}\left[R-Rf({\square }^{-1}R)\right],$$

$$f(X)=0.245\left[\tanh \left(0.0350(X+16.5)+0.032(X+16.5{)}^2+0.003(X+16.5{)}^3\right)-1\right],$$

$$\left(1-\frac{m^2}{\square }\right)G_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$S_{\mu\nu}=S_{\mu\nu}^{\rm T}+\frac{1}{2}({\mbox{\boldmath$\nabla$}}_{\mu}S_{\nu}+{\mbox{\boldmath$\nabla$}}_{\nu}S_{\mu})\,,$$

$$G_{\mu \nu }-m^2{\left({\square }^{-1}{G}_{\mu \nu }\right)}^{\mathrm{T}}=8\pi G{T}_{\mu \nu },$$

$$G_{\mu \nu }-\frac{m^2}{3}{\left(g_{\mu \nu }{\square }^{-1}R\right)}^{\mathrm{T}}=8\pi G{T}_{\mu \nu },$$

$${\Gamma }_{\mathrm{RR}}=\frac{m_{\mathrm{Pl}}^2}{2}\int d^{4}x\sqrt{-g}\left[R-\frac{m^2}{6}R\frac{1}{{\square }^2}R\right],$$

$$S_{\mathrm{Proca}}=\int d^{4}x\left[-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{1}{2}{m}_{\gamma }^2{A}_{\mu }{A}^{\mu }\right],$$

$$S_{\mathrm{nonlocal}}=-\frac{1}{4}\int d^{4}x\left[F_{\mu \nu }{F}^{\mu \nu }-m_{\gamma }^2{F}_{\mu \nu }\frac{1}{\square }{F}^{\mu \nu }\right],$$

$$\frac{m^2}{2}\int d^{4}x\mathrm{Tr}\left[F_{\mu \nu }\frac{1}{D^2}{F}^{\mu \nu }\right],$$

$$R\frac{1}{{\square }^2}R=36{\sigma }^2+\mathcal{O}({\sigma }^3).$$

$$\frac{\dot{G}}{G}=(7.1\pm 7.6)\times 10^{-14}{\mathrm{yr}}^{-1}.$$

$$\frac{\dot{G}}{G}=(0.99\pm 1.06)\times 10^{-3}\left(\frac{0.7}{h_0}\right)H_0.$$

$$G_{\mu \nu }+\Delta G_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$G_{\mu \nu }=8\pi G_{\mathrm{eff}}\left(T_{\mu \nu }-\frac{1}{8\pi G}{\mathcal{G}}_{\mu \nu }\right),$$

$$G_{\mathrm{eff}}(\varphi )=\frac{G}{1+\alpha (\varphi )}$$

$$H^2=\frac{8\pi }{3}{G}_{\mathrm{eff}}(\varphi )\left[\rho +\tilde{\rho }(\varphi )\right],$$

$$H^2=\frac{8\pi G}{3}\left[\rho +{\rho }_{\mathrm{DE}}(\varphi )\right].$$

$$d{s}^2=-(1+2\Psi )d{t}^2+(1+2\Phi )a^2(t)d{\mathbf{x}}^2,$$

$${\mbox{\boldmath$\nabla$}}^{2}\Phi=-4\pi G_{\rm eff}(\phi)a^{2}\,\left[\delta\rho+\delta\tilde{\rho}(\phi)\right]\,,$$

$$S_{\mu\nu}=S_{\mu\nu}^{\rm T}+\frac{1}{2}({\mbox{\boldmath$\nabla$}}_{\mu}S_{\nu}+{\mbox{\boldmath$\nabla$}}_{\nu}S_{\mu})\,.$$

$$(g_{\mu\nu}\Box^{-1}R)^{\rm T}=-Ug_{\mu\nu}-\frac{1}{2}({\mbox{\boldmath$\nabla$}}_{\mu}S_{\nu}+{\mbox{\boldmath$\nabla$}}_{\nu}S_{\mu})\,,$$

$$\displaystyle G_{\mu\nu}+\frac{m^{2}}{6}\,\left(2Ug_{\mu\nu}+{\mbox{\boldmath$\nabla$}}_{\mu}S_{\nu}+{\mbox{\boldmath$\nabla$}}_{\nu}S_{\mu}\right)$$

$$\square U$$

$$\displaystyle(\delta^{\mu}_{\nu}\Box+{\mbox{\boldmath$\nabla$}}^{\mu}{\mbox{\boldmath$\nabla$}}_{\nu})S_{\mu}$$

$$\frac{G_{\mathrm{eff}}}{G}=1+\mathcal{O}\left(\frac{1}{{\hat{k}}^2}\right),$$

$$G_{\mu \nu }$$

$$\square U$$

$$\square S$$

$$K_{\mu\nu}=2SG_{\mu\nu}-2{\mbox{\boldmath$\nabla$}}_{\mu}\partial_{\nu}S-2Ug_{\mu\nu}+g_{\mu\nu}\partial_{\rho}S\partial^{\rho}U-\frac{1}{2}g_{\mu\nu}U^{2}-(\partial_{\mu}S\partial_{\nu}U+\partial_{\nu}S\partial_{\mu}U)\,.$$

$$G_{\mathrm{eff}}=\frac{G}{1-(1/3)m^{2}S}.$$

Full text

**Testing nonlocal gravity with Lunar Laser Ranging

**

Enis Belgacem1, Andreas Finke1, Antonia Frassino1,2 and Michele Maggiore1

1Département de Physique Théorique and Center for Astroparticle Physics,

Université de Genève, 24 quai Ansermet, CH–1211 Genève 4, Switzerland

2Departament de Física Quàntica i Astrofísica, Institut de Ciències del Cosmos,

Universitat de Barcelona, Martí i Franquès 1, E-08028 Barcelona, Spain

We study the impact of the limit on $|\dot{G}|/G$ from Lunar Laser Ranging on “nonlocal gravity”, i.e. on models of the quantum effective action of gravity that include nonlocal terms relevant in the infrared, such as the “RR” and “RT” models proposed by our group, and the Deser-Woodard (DW) model. We elaborate on the analysis of Barreira et al. [1] and we confirm their findings that (under plausible assumptions such as the absence of strong backreaction from non-linear structures), the RR model is ruled out. We also show that the mechanism of “perfect screening for free” suggested for the DW model actually does not work and the DW model is also ruled out. In contrast, the RT model passes all phenomenological consistency tests and is still a viable candidate.

1 Introduction

The study of infrared (IR) modifications of General Relativity (GR) is largely motivated by the aim of obtaining a dynamical understanding of dark energy (DE). One direction of investigation, that has been very much explored, is based on the addition of new degrees of freedom to GR as, for instance, in massive gravity, bigravity, or scalar-tensor theories. A different route consists in asking whether, already in GR, quantum effects could modify the IR behavior and give rise to interesting cosmological consequences. The idea that quantum gravity at large distances could induce cosmological effects is quite old, see e.g. Taylor and Veneziano [2]. A first-principle understanding of the IR structure of quantum gravity is a highly non-trivial task, but hints of non-trivial physics in space-times of cosmological relevance emerge for instance from the study of their IR divergences [3], from the fact that strong IR effects appear for the conformal mode [4, 5], and from arguments against the quantum stability of de Sitter space [6, 7, 8, 9, 10] (see also recent discussions in the context of the “swampland” [11, 12, 13]).

When studying quantum effects, the relevant quantity is the quantum effective action, rather than the original fundamental action of the theory. Working in terms of the quantum effective action opens new phenomenological possibilities, related to the introduction of nonlocal terms. Indeed, according to the basic principles of QFT, the fundamental action of a theory must be local. However, when in the spectrum there are massless particles (such as the graviton in GR), or more generally particles light with respect to the relevant energy scales, the quantum effective action has also nonlocal terms. In particular, terms involving inverse powers of the d’Alembertian are relevant in the IR and can affect the behavior of the theory at cosmological scales. Given that a first-principle computation of the IR limit of the quantum effective action of gravity is quite difficult, especially if one is looking for non-perturbative effects, it makes sense to begin by adopting a phenomenological point of view and study what sort of nonlocal terms could give rise to a viable theory of dark energy.

In this spirit, a first nonlocal gravity model was proposed by Wetterich [14], who considered a quantum effective action of the form

[TABLE]

where $\Box$ is the covariant d’Alembertian.111More precisely, in [14] was considered a term of the form $R[\Box+\xi R+{\cal O}(R^{2})]^{-1}R$ in Euclidean space. Since $\Box^{-1}R$ is dimensionless, the associated constant $g^{2}$ is also dimensionless. The model however did not produce a viable cosmological evolution [14]. Deser and Woodard [15, 16] considered a more general nonlocal model of the form

[TABLE]

with $f(X)$ a dimensionless function (see [17] for review). This function is tuned so to obtain the desired background evolution. Requiring that the cosmological evolution closely mimicks that of $\Lambda$CDM leads to a function well-fitted by [18]

[TABLE]

for $X\equiv\Box^{-1}R<0$. To comply with solar system constraints, the proponents of the model set $f(X)=0$ for $X>0$, a crucial point to which we will come back later. At the level of background evolution the DW model is not predictive, given that $f(X)$ can be chosen to match any desired expansion. However, cosmological perturbations over the background determined by eq. (1.3) have been studied, and it has been found that the model can fit structure formation data [19, 20]. Similar nonlocal models, involving operators such as $R^{\mu\nu}\Box^{-1}G_{\mu\nu}$, have been discussed in [21, 22, 23].

Both the Wetterich and the Deser-Woodard proposals involve only dimensionless parameters associated to the nonlocal term. Another interesting possibility is that the nonlocal term involves a mass scale. Such a mass scale could be generated dynamically in quantum gravity within GR itself (or else be induced in the quantum effective action in extensions of GR where there is already an explicit mass or length scale in the fundamental action, such as massive gravity or theories with extra dimensions). Some indication for the appearance of an IR mass scale in GR at the non-perturbative level comes from lattice gravity [24, 25, 26, 27] and a recent simulation using causal dynamical triangulations [28]. Recently it has also been pointed out that a dynamical IR mass scale also appears using exact renormalization group equations, as a reflection of the instability of the conformal mode in Euclidean space [29]. A further hint comes from $N=1$ supergravity compactified on a circle, which is a model that allows one to perform non-perturbative computations that otherwise, in general, are not feasible. In this setting, an IR mass scale emerges because of instanton effects, from the running of the coefficient of the Gauss-Bonnet term [30]. While these results must be taken with a grain of salt, due to the difficulty of dealing with non-perturbative effects in gravity, still, they can at least be considered as a proof of principle of the fact that a mass scale could be dynamically generated in quantum gravity.

At a purely phenomenological level, a first nonlocal model associated to a mass scale was proposed by Arkani-Hamed, Dimopoulos, Dvali, and Gabadadze [31] to introduce the degravitation idea, and consisted in modifying the Einstein equations into

[TABLE]

where $m$ is the new mass scale.222Actually, to preserve causality, this must be understood as the equation of motion derived from a quantum effective action for the in-in vacuum expectation value of the metric, which automatically ensures that the Green’s function in the $\Box^{-1}$ operator is the retarded one; see [32, 33] for discussion. The same will hold for similar equations, such as eq. (1.7) that defines the RT model. However, a drawback of eq. (1.4) is that the energy-momentum tensor is not automatically conserved, since in curved space [{\mbox{\boldmath\nabla}}^{\mu},\Box^{-1}]\neq 0 and therefore the Bianchi identity {\mbox{\boldmath\nabla}}^{\mu}G_{\mu\nu}=0 no longer ensures {\mbox{\boldmath\nabla}}^{\mu}T_{\mu\nu}=0. In [34] it was however observed that it is possible to cure this problem by making use of the fact that any symmetric tensor $S_{\mu\nu}$ can be decomposed as

[TABLE]

where $S_{\mu\nu}^{\rm T}$ is the transverse part of $S_{\mu\nu}$, {\mbox{\boldmath\nabla}}^{\mu}S_{\mu\nu}^{\rm T}=0. Such a decomposition can be performed in a generic curved space-time [35, 36]. One can then modify eq. (1.4) into

[TABLE]

where the superscript “T” denotes the extraction of the transverse part from of the tensor $\left(\Box^{-1}G_{\mu\nu}\right)$. Energy–momentum conservation then becomes automatic. However, the cosmological evolution of this model turned out to be unstable, already at the background level [37, 38], so this model is not phenomenologically viable.333This instability is due to some growing modes in the homogeneous equation for the auxiliary fields used to put the model in a local form; see [32, 33] and section 3.1 for review. Then, in [37] a model was considered based on the nonlocal equation of motion

[TABLE]

where $\Box^{-1}$ is rather applied the Ricci scalar $R$ (the factor $1/3$ is a useful normalization for the mass scale $m$). We refer to it as the “RT” model, where R stands for the Ricci scalar and T for the extraction of the transverse part. A closed form for the quantum effective action from which eq. (1.7) could be derived is currently not known. This model is however related to another nonlocal model, proposed in [39], and defined by the quantum effective action

[TABLE]

that we denote as the RR model, after the two occurrences of the Ricci scalar in the nonlocal term. The RT and RR models are related by the fact that they coincide when linearized over Minkowski space. However beyond linear order, or when linearized over a different background such as FRW, they are different.

It turns out that, in both the RR and RT models, the auxiliary fields used to put the model in local form have a stable cosmological evolution at the background level [the condition that ruled out the model (1.6)]. In both models the nonlocal term acts effectively as dark energy, producing accelerated expansion in the recent epoch, and their cosmological perturbations are well-behaved in all phases of the cosmic evolution [40] (e.g. an early inflationary phase, radiation dominance, matter dominance and the recent dark-energy dominated epoch).444In a paper of our group [41] we incorrectly stated that the RT model is not viable if its evolution is started from an early inflationary phase, because of an instability at the level of scalar perturbations. This statement has been corrected in the v3 arXiv version of that paper, where it is shown that this apparent instability is numerically irrelevant because of the smallness of the scale associated to the nonlocal term compared to the inflationary scale. The RT model is then perfectly viable even if its evolution is started during an early inflationary phase. Detailed comparison with CMB, BAO, SNe and structure formation data [42, 43, 44, 33] show that both models fit the data at a level statistically indistinguishable from $\Lambda$CDM. Furthermore, in these models the tensor perturbations propagate at the speed of light, as they should to comply with the limits from the neutron star binary coalescence GW170817, but still with a different “friction term”, which leads to modified GW propagations and distinct predictions that could be tested with LISA or with the Einstein Telescope [45, 46]. Thus, these models have a very interesting phenomenology.

As for their physical meaning, it is interesting to observe that nonlocal terms of this type allow us to write down mass terms for gauge fields or for some components of the gravitational field, while preserving gauge or diffeomorphism invariance, respectively. This was observed in [47] for massive electrodynamics, where it was shown that the Proca action for a massive photon,

[TABLE]

in which the mass term breaks explicitly the $U(1)$ gauge invariance, is equivalent to the action

[TABLE]

which is nonlocal but gauge-invariant. Indeed, in  [48, 49, 50] it has been found that, in pure Yang-Mills theory, non-perturbative results obtained from lattice QCD for the gluon propagator and the running of the strong coupling constant can be reproduced, at the level of the quantum effective action, by the introduction of the non-abelian generalization of the above nonlocal term,

[TABLE]

(where $D_{\mu}^{ab}=\delta^{ab}\partial_{\mu}-gf^{abc}A_{\mu}^{c}$ is the covariant derivative). In this case the nonlocal term is generated dynamically in the IR by the strong interaction, and has the meaning of an effective gluon mass.

The nonlocal terms in eqs. (1.7) and (1.8) have a similar interpretation as a gauge-invariant mass for the conformal mode of the metric [51, 52]. Indeed, if we restrict for simplicity the dynamics to the conformal mode $\sigma(x)$, writing $g_{\mu\nu}(x)=e^{2\sigma(x)}\eta_{\mu\nu}(x)$, for the Ricci scalar we get $R=-6e^{-2\sigma}\left(\Box\sigma+\partial_{\mu}\sigma\partial^{\mu}\sigma\right)$. Therefore, to linear order in $\sigma$, $R=-6\Box\sigma+{\cal O}(\sigma^{2})$ and, upon integration by parts,

[TABLE]

Thus, the structure $R\Box^{-2}R$ in eq. (1.8) corresponds to a nonlocal but diffeomorphism-invariant mass term for the conformal mode, plus higher-order interaction terms (that are nonlocal even when written in terms of $\sigma$), that are required to reconstruct a diffeomorphism-invariant quantity.555Observe that, in contrast, since to linear order $R\Box^{-1}R\simeq 36\sigma\Box\sigma$, the nonlocal term in eq. (1.1) corresponds to a kinetic term for the conformal mode. The same conclusion holds for the RT model since it coincides with the RR model when linearized over Minkowski space (although the higher-order interaction terms differ). Therefore, as long as we consider the RT model linearized over Minkowski, we can consider eq. (1.7) as derived from the RR quantum effective action plus higher-order non-linear terms. Since nonlocal terms such as $R\Box^{-2}R$ cannot appear in the fundamental action, but only at the level of the quantum effective action, they describe mass terms generated dynamically by quantum effects. It is quite interesting to observe that numerical results from causal dynamical triangulations recently presented in [28] provide some evidence in favor of a mass term for the conformal mode, of the type described by eqs. (1.7) or (1.8).

Other variants of the RR model have been explored, such as terms in the effective action of the form $m^{2}R^{\mu\nu}\Box^{-2}R_{\mu\nu}$ or $m^{2}C_{\mu\nu\rho\sigma}\Box^{-2}C^{\mu\nu\rho\sigma}$, where $C_{\mu\nu\rho\sigma}$ is the Weyl tensor [53]. However, the former does not have a stable background cosmological evolution, just as the model (1.6); the latter does not affect the background evolution in FRW, and has well-behaved scalar perturbations, but it turns out to be unstable at the level of tensor perturbations. Observe that these structures also give a mass to the helicity-2 mode, so a physical mechanism that generates only a mass for the conformal mode would not produce them.666In this context, it is useful to clear up a misconception that occasionally appears in discussions. The quantum effective action should not be confused with a low-energy effective action. The quantum effective action is an object that includes the effects of quantum fluctuations (in principle exactly, if one were able to compute it), and is valid at all energy scales at which the original fundamental action is valid. To compute the quantum effective action one is not integrating out some fields from the theory, but rather one is integrating out their quantum fluctuations, remaining with a functional of the vacuum expectation values of all fields (see [32] for a pedagogical discussion). While a low-energy effective action in general has all terms consistent with the symmetries of the fundamental action, the quantum effective action will have a very specific structure. For instance, some mode could become massive while others could remain massless. Thus, it makes perfect sense to study a specific nonlocal structure; assuming that non-perturbative effects generate a term proportional to $R\Box^{-2}R$ does not imply that then also a term such as, e.g., $R^{\mu\nu}\Box^{-2}R_{\mu\nu}$ will automatically appear, even more considering that these terms have different physical meaning. Given that models that have the physical meaning of a mass for the conformal mode appear to be the viable ones, in [54] a model with a term $m^{2}R\Delta_{4}^{-1}R$ was also studied, where \Delta_{4}=\Box^{2}+2R^{\mu\nu}{\mbox{\boldmath\nabla}}_{\mu}{\mbox{\boldmath\nabla}}_{\nu}-\frac{2}{3}R\Box+\frac{1}{3}g^{\mu\nu}{\mbox{\boldmath\nabla}}_{\mu}R{\mbox{\boldmath\nabla}}_{\nu} is called the Paneitz operator. This operator enters the quantum effective action obtained by integrating the conformal anomaly in four dimensions [55, 56]. In flat space it reduces to $\Box^{2}$, so again this nonlocal term has the meaning of a mass for the conformal mode. The model turns out to have a viable evolution at the level of background, as well as in the scalar perturbations sector. However, the study of its tensor perturbations revealed that GWs do not propagate at the speed of light, and therefore the $\Delta_{4}$ model is ruled out by GW170817 [33]. A review of the conceptual aspects and phenomenological consequences of these nonlocal models associated to a mass scale is given in [32, 33] (see also [57, 58, 59, 60, 61, 62, 63, 64] for further works).

These considerations show that the construction of a phenomenologically viable nonlocal model is a highly non-trivial task. In this paper, elaborating on the discussion in [1], we examine another constraint on these models, coming from the bound on the time variation of Newton’s constant obtained from Lunar Laser Ranging (LLR). We will find that, within some caveats to be discussed below, this bound rules out both the DW model and the RR model, while the RT model survives.

The paper is organized as follows. In section 2 we review the current status of the limit on $\dot{G}/G$ from LLR, and we compare with a first ‘naive’ analysis for the RR model (while we will see that the RT model is not concerned by the bound, and passes the LRR test automatically). In section 3 we put the analysis for the RR model on firmer grounds, by examining different ways of matching a cosmological solution at large scales to a static solution at short scales, and we will conclude that the RR model is indeed ruled out. In section 4 we analyze a screening mechanism that was proposed for the Deser-Woodard model. We show that this mechanism actually does not work and that the Deser-Woodard model is also ruled out. We present our conclusions in section 5. Some technical material is relegated in three appendices. We use units $c=1$ and MTW conventions for the metric signature, Riemann tensor, etc.

2 Limits from $\dot{G}/G$ and modified gravity

Lunar Laser Ranging provides a stringent limit on the time variation of Newton’s constant. The most recent result is [65]

[TABLE]

Note that the LLR limit improves significantly with observation time and the above limit is one order of magnitude stronger than the previously published limit $(4\pm 9)\times 10^{-13}\,{\rm yr}^{-1}$ [66], which was also used in [1] to discuss constraints on nonlocal gravity. The cosmological implications of this result become more evident if we re-express this value in terms of the Hubble parameter today $H_{0}=h_{0}\times 100\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\simeq h_{0}\times(9.777752\,{\rm Gyr})^{-1}$. Then, eq. (2.1) reads

[TABLE]

Quite generally, in modified gravity models Newton’s constant becomes time dependent on cosmological scale. The scale for the time variation today is given by $H_{0}$, so, generically, on cosmological scales one finds $\dot{G}/G\simeq H_{0}$, with a coefficient typically of order one. If, in a given modified gravity model, such a result persists down to the scale of the solar system and of the Earth-Moon system, then the model is ruled out by eq. (2.2). A modified gravity models can however be immune to this bound, either because it simply does not induce a time-dependence in Newton’s constant at small scales (we will see that this is the case for the RT model), or because it exhibits a non-linear screening mechanism, such as the Vainshtein mechanism or the chameleon mechanism (see e.g. [67] for review).

2.1 Effective Newton’s constant in modified gravity

Before specializing to the non-local models that will be the main focus of our work, it is useful to discuss in more general terms why in modified gravity Newton’s constant is typically modified, and what is the relation between the effective Newton constant measured on cosmological scales and that measured from the motion of test masses on shorter scales, as in the case when we take the Moon as a test mass in the field of the Earth. Generically, in modified gravity the Einstein equations become

[TABLE]

where $\Delta G_{\mu\nu}$ can depend on extra fields of the theory, e.g. on some scalar fields in the case of scalar-tensor theories. As we will recall below, in non-local gravity $\Delta G_{\mu\nu}$ rather depends on non-local functional of the metric, and can be rewritten as a local expression in terms of some auxiliary fields. The set of equations of the theory is then completed by the equations of motion of the extra fields (or, in non-local gravity, by equations that follow from the definition of the auxiliary fields, see below). Suppose that, among the various terms in $\Delta G_{\mu\nu}$, there will be terms proportional to the Einstein tensor $G_{\mu\nu}$ itself, so $\Delta G_{\mu\nu}=\alpha(\phi)G_{\mu\nu}+{\cal G}_{\mu\nu}$, where $\alpha(\phi)$ is some function of the various fields or auxiliary fields of the theory, here denoted collectively as $\phi$ (and of the metric, that we do not write explicitly), and ${\cal G}_{\mu\nu}$ denotes the remainder, that is not proportional to $G_{\mu\nu}$. Of course, this assumes a choice of basis in the space of operators; in particular, it is convenient to take $G_{\mu\nu}$ and $g_{\mu\nu}$ as independent (rather than, say, $G_{\mu\nu}$ and $R_{\mu\nu}$) because the former, as we will see, contributes to the effective Newton constant, while terms proportional to $g_{\mu\nu}$ contributes to a cosmological constant or a potential term, so this separation has a more clear physical interpretation. Equation (2.3) can then be rewritten as

[TABLE]

where

[TABLE]

plays the role of an effective Newton’s constant. Whether this separation of $\Delta G_{\mu\nu}$ into $\alpha(\phi)G_{\mu\nu}$ and ${\cal G}_{\mu\nu}$ is useful or not depends on the context. For instance, at the level of background evolution, eq. (2.4) produces a modified Friedmann equation of the form

[TABLE]

where $\tilde{\rho}(\phi)$ is the contribution from ${\cal G}_{00}$. However, the time dependence of $\tilde{\rho}(\phi)$ will in general be as important as that in $G_{\rm eff}(\phi)$, so this separation is not very useful. It is more meaningful to put all the contribution from $\Delta G_{\mu\nu}$ together, defining an effective dark energy density $\rho_{\rm DE}(\phi)$ from $G_{\rm eff}(\phi)\left[\rho+\tilde{\rho}(\phi)\right]=G\left[\rho+\rho_{\rm DE}(\phi)\right]$, so that eq. (2.6) reads

[TABLE]

At the level of cosmological perturbations, this separation can however be more useful. Indeed, working with a perturbed FRW metric

[TABLE]

where $\Phi(t,{\bf x})$ and $\Psi(t,{\bf x})$ are the usual Bardeen potentials, the linearization of the $(00)$ component of eq. (2.4), for modes inside the horizon, gives a modified Poisson equation of the form

[TABLE]

where again $\delta\tilde{\rho}(\phi)$ is the contribution from $\delta{\cal G}_{00}$. However, in several modified gravity models it happens that the term $\delta\tilde{\rho}(\phi)$ is negligible for modes well inside the horizon. As we will recall below, this is in particular the case for the RR and RT nonlocal models. In that case, $G_{\rm eff}$ defined from eq. (2.4) captures the main effect on the growth of the perturbations [in general, together with a second indicator such as the gravitational slip $(\Phi+\Psi)/\Phi$].

However, the situation where this separation is most useful is when we study a modified gravity model at the scale of the solar system or the Earth-Moon system, as for LLR. Indeed, in this case, the term ${\cal G}_{\mu\nu}/(8\pi G)$ in eq. (2.4) is totally negligible with respect to $T_{\mu\nu}$. For instance, its $(00)$ component is of the order of the energy density of dark energy and, at the solar system scales, as an additive contribution its effect on the dynamics is totally irrelevant compared to $T_{00}$, which is the energy density $\rho$ of the bodies, such as the Earth, that dominate the dynamics at the scales relevant for LLR. In contrast, the effective Newton constant gives a multiplicative contribution to $\rho$. Even if numerically tiny, the time dependence that it induces can be tested through the very accurate LLR bound (2.2).

To sum up, in modified gravity the effective Newton constant relevant for comparison with LLR can be read directly from the modified Einstein equations, by selecting the terms proportional to $G_{\mu\nu}$ in $\Delta G_{\mu\nu}$, and its functional dependence on the fields $\phi$, given by eq. (2.5), is the same as that of the effective Newton constant that governs linear cosmological perturbation theory through the modified Poisson equation (2.9). In general, different screening mechanisms will be possible depending on the specific dynamics of the fields $\phi$, which could suppress the value of $\alpha(\phi)$ in the regime relevant for solar system tests. We next examine some concrete cases, namely the RT, RR and DW models.

2.2 RT model

Let us first recall how the RT model can be recast in local form [37]. We begin by introducing an auxiliary field $U=-\Box^{-1}R$ and we define $S_{\mu\nu}=-Ug_{\mu\nu}=g_{\mu\nu}\Box^{-1}R$.777The notation that has been traditionally used for the RR and RT model is different from the one that is standard in the DW model. In particular, the auxiliary field called $U$ in the RT and RR models correspond to $-X$ in the notation used for the DW model. We then split this tensor into a transverse part $S_{\mu\nu}^{\rm T}$, which by definition satisfies {\mbox{\boldmath\nabla}}^{\mu}S_{\mu\nu}^{\rm T}=0, and a longitudinal part,

[TABLE]

Thus, the nonlocal quantity $S_{\mu\nu}^{\rm T}=(g_{\mu\nu}\Box^{-1}R)^{\rm T}$ that appears in eq. (1.7) can be written as

[TABLE]

in terms of an auxiliary scalar field $U$ and an auxiliary four-vector field $S_{\mu}$. Then the equation of motion (1.7) can be rewritten as the coupled system of equations

[TABLE]

where the last equation is obtained by taking the divergence of eq. (1.5). We now observe that eq. (2.12) has the form (2.3), but with a $\Delta G_{\mu\nu}$ in which there no term proportional to $G_{\mu\nu}$. Therefore, according to the discussion in section 2.1, the effective Newton constant relevant for LLR is equal to $G$.

This can be confirmed by expanding the equations to linear order over FRW, and reading the effective Newton constant from the Poisson equation. For the RT model this computation has been performed in [58, 40]. The result is that, for sub-horizon modes,

[TABLE]

where $\hat{k}=k/(aH)$ is the physical momentum $k/a$ divided by $H$. This confirms again that, for sub-horizon modes, $G_{\rm eff}$ reduces to $G$ and the correction, for the modes relevant to LLR, is totally negligible. Thus the RT model passes the test (2.2).

2.3 RR model

We next consider the RR model. Let us recall that, for this model, it is convenient to introduce two auxiliary fields $U$ and $S$ defined from $U=-\Box^{-1}R$ and $S=-\Box^{-1}U$ [39]. Then, the equation of motion of the theory become equivalent to the following coupled system

[TABLE]

where

[TABLE]

Now in the $K_{\mu\nu}$ there is a term proportional to $G_{\mu\nu}$. As a consequence, the effective Newton constant relevant for LLR is given by

[TABLE]

Again, we can compare the the result obtained by expanding eqs. (2.16)–(2.18) to linear order over FRW. The result, for sub-horizon modes, is [40, 1]

[TABLE]

where $\bar{S}(t)$ is the background cosmological solution for the auxiliary field $S$. Equation (2.21) shows that the effective Newton constant that enters the Poisson equation for the cosmological perturbations has a dependence on time through the background value of the auxiliary field $\bar{S}(t)$. The effective Newton constant relevant for LLR, given by eq. (2.20), has the same functional dependence on $S$, as expected from the discussion in section 2.1, so in this case there is a potential problem with LLR, as observed in [1]. Of course, the crucial issue, also recognized in [1], is that it is not evident whether the cosmological time-dependent solution for $S$ is still appropriate at solar system or LLR scales, where one could rather imagine that the cosmological solution matches to a static solution $S_{\rm static}(r)$, i.e. it is not obvious whether, when we apply eq. (2.20) to LLR, we have the right to replace $S$ by the cosmological solution $\bar{S}(t)$.

We will examine this issue in detail in the following sections. For the moment, if one just simply assumes that eq. (2.21) is applicable at LLR scales and uses the background solution for the $S$ field obtained from the numerical integration of the equations of motion (2.16)–(2.18), one finds $\dot{G}/G\simeq 0.12\,H_{0}$. As observed in [1], this already exceed by a factor ${\cal O}(10)$ the bound $(4\pm 9)\times 10^{-13}\,{\rm yr}^{-1}$ [66], that in terms of $H_{0}$ reads $\dot{G}/G\simeq(5\pm 12)\times 10^{-3}H_{0}(0.7/h_{0})$. The more recent bound (2.2) is then exceeded by a factor ${\cal O}(100)$.

Actually, as we will recall in section 3.1, the evolution in the RR model is also characterized by one free parameter, which is the initial value $u_{0}$ of the auxiliary field $U$, at an initial time set deep in radiation dominance, and represents a marginal direction in parameter space (while, in the space of initial conditions $\{U(t_{i}),\dot{U}(t_{i}),S(t_{i}),\dot{S}(t_{i})\}$, the remaining three directions correspond to irrelevant directions in parameter space). The above result for $\dot{G}/G$ was obtained for the “minimal model”, defined by $u_{0}=0$. Larger values of $u_{0}$ bring the evolution of the model closer and closer to that of $\Lambda$CDM, and could be naturally generated by a phase of primordial inflation. Indeed, $U$ grows during inflation in such a way that, starting from a value of order one at the beginning of inflation, at the end of inflation it reaches a value $U\simeq 4\Delta N$, where $\Delta N$ is the number of efolds. This value can then be taken as the initial value $u_{0}$ for the subsequent evolution during radiation dominance. For $\Delta N=60$ this gives a large value $u_{0}\simeq 240$. It is therefore interesting to see if the bound (2.2) can be avoided with large values of $u_{0}$. We have therefore computed $\dot{G}/G$ from the numerical evolution for several values of $u_{0}$. Using our standard notations (see section 3.1) $\gamma\equiv m^{2}/(9H_{0}^{2})$ and $V\equiv H_{0}^{2}S$, parametrizing the time evolution in terms of $x\equiv\ln a(t)$, and denoting by a prime the derivative with respect to $x$, for modes well inside the horizon we have

[TABLE]

and therefore

[TABLE]

Thus, at the present time, $\left(\dot{G}_{\rm eff}/{G_{\rm eff}}\right)_{0}=\alpha H_{0}$, where

[TABLE]

and the subscript [math] denotes the present time. This quantity is easily computed from the solution for $\bar{V}(x)$ obtained from the numerical integration of the equations of motion of the model, specialized to a FRW background [39, 32, 33], that we will recall in section 3.1. In Fig. 1 we compare the numerical results for $\alpha$, obtained for different values of the initial condition $u_{0}$, with the observational limit (2.2). We see that, under the assumption that eq. (2.21) is indeed still valid at the scale of the Earth-Moon system, the RR model is ruled out even for very large values of $u_{0}$, say $u_{0}\simeq 500$. Of course, with a sufficiently large value of $u_{0}$, the model would eventually comply with the limit. According to the fit $\alpha\simeq 2.11/u_{0}$ (red line), this requires u_{0}\,\raise 1.72218pt\hbox{>}\kern-6.99997pt\lower 2.6694pt\hbox{\sim}\,2000. However, such values becomes unnaturally large even from the point of view of a previous inflationary phase. Furthermore, the model would become so close to $\Lambda$CDM that it would become basically indistinguishable from $\Lambda$CDM in its cosmological predictions.

The question is therefore whether the extrapolation of eq. (2.21) from cosmological scales down to the scale of the Earth-Moon system is justified. We address this issue in section 3.

2.4 DW model

We next consider the DW model. The variation of eq. (1.2) gives a modified Einstein equation of the form (2.3), with [15, 16]

[TABLE]

where $f^{\prime}(X)=\partial f/\partial X$. The model can be written in local form by introducing two auxiliary fields $X=\Box^{-1}R$ and $Y=\Box^{-1}\left[Rf^{\prime}(X)\right]$, so that the equations for the auxiliary fields are $\Box X=R$ and $\Box Y=Rf^{\prime}(X)$.888In the literature on the DW model the auxiliary field $Y$ is denoted by $U$. We change the notation to avoid confusion with the auxiliary field denoted by $U$ in the RR and RT models, that rather corresponds to $-X$. Picking the term proportional to $G_{\mu\nu}$ one finds the effective Newton constant [17]

[TABLE]

The same result can be found by writing the Poisson equation for the perturbations around FRW in the large-$k$ limit, which gives [19, 20]

[TABLE]

where as usual the overbar denotes the field evaluated on the background cosmological solutions. Note that in this case, contrary to the RR and RT models, besides an effective Newton constant there is also an extra contribution proportional to $[\delta U+f^{\prime}(\bar{X})\delta X]$, that does not vanish in the sub-horizon limit, i.e. the term $\delta\tilde{\rho}(\phi)$ in eq. (2.9) is not suppressed for large $k$.

The time dependence of the auxiliary fields $X$ and $Y$ therefore potentially induces a time dependence in the Newton constant that govern the dynamics at LLR scales. In the attempt at introducing a screening mechanism, the proponents of the model set $f(X)=0$ when $X>0$, a point that will be discussed in detail in section 4.

3 Matching of cosmological and static solutions

A first natural reaction to the result found in section 2.3 for the RR model is that, after all, the FRW metric has no relevance for the Earth-Moon system, which is rather governed by a static Schwarzschild metric. However, as correctly observed in [1], the issue is more subtle. The complete metric will, of course, be one that, at the largest cosmological scales, reduces to FRW and, as we go down in scales, is gradually affected by the various structures and inhomogeneities of the Universe (superclusters, clusters, galaxies, etc.), until at the solar system scale is equal to the Schwarzschild metric generated by the Sun (plus the corrections from the other relevant bodies), and finally at the Earth-Moon scale it reduces to the static metric generated by the Earth. The auxiliary fields $U$ and $S$ of the RR model satisfy eqs. (2.17) and (2.18), where the $\Box$ operator and the Ricci scalar are computed in this metric. At the largest scales, where the metric is given by FRW, we get the corresponding cosmological solutions $U_{\rm cosmo}(t)$ and $S_{\rm cosmo}(t)$, that depend on time but not on ${\bf x}$, and whose explicit form we will recall in section 3.1. At solar system or Earth-Moon scales, we can study the solution of eqs. (2.17) and (2.18) in a static Schwarzschild metric. In this background we can look for static solutions, that we denote here by $U_{\rm static}(r)$ and $S_{\rm static}(r)$, and that will be discussed in section 3.2 and app. A, extending the discussion of refs. [39, 68]. These solutions are static by construction since they have been found by looking for solutions with no time dependence. However, it is not obvious whether, in the full metric that interpolates between the FRW metric at the largest cosmological scales and the Schwarzschild metric at the Earth-Moon scale, the complete solution for $S$ will correspondingly interpolate between $S_{\rm cosmo}(t)$ and $S_{\rm static}(r)$. It is in principle possible that at the Earth-Moon scale the solution for $S$ could still retain a time dependence, inherited from the time evolution at large scales (and which is lost by definition when searching for a purely static solution, matched at infinity to flat space). If this is the case, then Newton’s constant will indeed have a time dependence at the Earth-Moon scale relevant for the LLR observations. Indeed, this is precisely what happens in models, such as galileons or k-essence, in which a field $\varphi$ has a shift symmetry $\varphi\rightarrow\varphi+{\rm const}$. In this case, thanks to the shift symmetry, near the present epoch $t_{0}$ the equation of motion of the field admits a solution with separation of variables of the form $\varphi(t,r)=\varphi_{\rm static}(r)+\varphi_{\rm cosmo}(t_{0})+(t-t_{0})\dot{\varphi}_{\rm cosmo}(t_{0})$ [69]. The RR model does not have such a shift symmetry, so the applicability of this result is not obvious, but, still, this shows that the issue is non-trivial.

An equivalent way of posing the problem is to follow the evolution of the fields $S$ and $U$ starting from radiation dominance (RD), before structure formation. At that time the FRW metric holds everywhere. The results that we will review in section 3.1 reveal that during RD the auxiliary fields still show no significant time evolution. However, as soon as we enter in matter dominance (MD), the $U$ field begins to evolve with time. It is not clear what would stop the time evolution of these fields in the regions where the metric is dominated by the contribution of massive bodies, but still the overall gravitational field is weak, which is the case everywhere except in the vicinity of compact objects such as black holes or neutron stars.

To understand the problem, we need to study the equations for $U$ and $S$ in a metric that interpolates between the cosmological and the static solution. We begin however by first recalling, in the next two subsections, the form of the cosmological and static solutions for the auxiliary fields in the RR model.

3.1 Cosmological solutions for the RR model

Let us recall the form of the cosmological solutions for the RR model [39, 40, 33]. We consider a flat FRW metric $ds^{2}=-dt^{2}+a^{2}(t)d{\bf x}^{2}$, and we look for cosmological background solutions for the auxiliary field that depend only on time, $U_{\rm cosmo}(t)$ and $S_{\rm cosmo}(t)$ in the notation of the above discussion. For ease of notation, however, in this subsection we write $U_{\rm cosmo}(t)=\bar{U}(t)$ and $S_{\rm cosmo}(t)=\bar{S}(t)$, and in general we use an overbar for quantities corresponding to the background cosmological evolution. We define $h(t)=H(t)/H_{0}$, where $H(t)=\dot{a}/a$ and $H_{0}$ is the present value of the Hubble parameter, we use $x=\ln a(t)$ to parametrize the temporal evolution, and we denote $df/dx\equiv f^{\prime}$. Observe that $U$ is dimensionless while $S=-\Box^{-1}U$ has dimension of $({\rm length})^{2}$. To show that the nonlocal term gives rise to an effective dark energy density it is convenient to introduce the dimensionless field $W(t)=H^{2}(t)S(t)$. The equations of motion (2.16)–(2.18) then become

[TABLE]

where $\Omega_{M},\Omega_{R}$ are the present values of $\rho_{M}/\rho_{\rm tot}$ and $\rho_{R}/\rho_{\rm tot}$, respectively, $\gamma=m^{2}/(9H_{0}^{2})$, $\zeta=h^{\prime}/h$ and $\bar{Y}$ is a known functional of $\bar{U},\bar{V}$, whose explicit expression is $\bar{Y}\equiv(1/2)\bar{W}^{\prime}(6-\bar{U}^{\prime})+\bar{W}(3-6\zeta+\zeta\bar{U}^{\prime})+(1/4)\bar{U}^{2}$. In this form, one sees that there is an effective dark energy density $\rho_{\rm DE}=\rho_{0}\gamma\bar{Y}$ where, as usual, $\rho_{0}=3H_{0}^{2}/(8\pi G)$. Actually, to perform the numerical integration of these equations, and also to study the cosmological perturbations, it can be more convenient to use a variable $V(t)=H_{0}^{2}S(t)$ instead of $W(t)=H^{2}(t)S(t)$. Then eqs. (3.1)–(3.3) become

[TABLE]

Given the initial conditions, that we will set for definiteness deep in RD, the numerical integration of these equations is straightforward. The initial conditions on the auxiliary fields $U$ and $V$ (or, equivalently, $U$ and $W$), are in one-to-one correspondence with the solutions of the homogeneous equations associated to eqs. (3.2) and (3.3), i.e.

[TABLE]

In any given cosmological epoch $\zeta$ has an approximately constant value $\zeta_{0}$, with $\zeta_{0}=\{0,-2,-3/2\}$ in de Sitter (dS), RD and MD, respectively. Taking $\zeta$ constant the homogeneous solutions are obtained analytically,

[TABLE]

During the whole cosmological evolution $-2\leq\zeta_{0}\leq 0$, so $u_{1},w_{0}$ and $w_{1}$ parametrize irrelevant directions in parameter space. Note also that there is no exponentially growing solution, so no instability at this level of the analysis.999This is the condition that was not fulfilled by the model (1.6). In contrast, $u_{0}$ parametrizes a marginal direction in parameter space and can be taken as a free parameter of the model. To have an idea of a natural value for this parameter we observe, from eq. (3.5) with $\zeta_{0}=0$, that, during a pre-existing phase of de Sitter inflation, $\bar{U}(x)=4x$. Then, if $\bar{U}(x)={\cal O}(1)$ at the beginning of inflation, by the end of inflation $\bar{U}(x)\simeq 4\Delta N$, where $\Delta N$ is the number of inflationary efolds. Thus $u_{0}$, which is the initial value for the subsequent evolution starting from RD, could naturally be expected to be of order $4\Delta N$. For typical inflationary models $\Delta N={\cal O}(50-60)$, so inflation could naturally produce a value of $u_{0}={\cal O}(200-250)$. We can then study the evolution for the simplest choice $u_{0}=0$ [which is representative of any value $u_{0}={\cal O}(1)$] and for a large value of $u_{0}$, say $u_{0}=200$.

The result of the numerical integration is shown in Fig. 2 for $u_{0}=0$ and in Fig. 3 for $u_{0}=200$. For illustration, in both cases the parameter $\gamma$ is adjusted so to get a value $\rho_{\rm DE}/\rho_{0}=0.7$ today, which is obtained by setting $\gamma\simeq 9.126\times 10^{-3}$ for $u_{0}=0$ and $\gamma\simeq 5.907\times 10^{-5}$ for $u_{0}=200$.101010Of course, for each value of $u_{0}$, the exact value of $\rho_{\rm DE}/\rho_{0}$, and therefore $\gamma$, must eventually be fixed by comparing the model with a set of cosmological data and performing Bayesian parameter estimation, see the analysis in [42, 43, 44]. In terms of the variable $x=\ln a$, radiation-matter equilibrium corresponds to $x\simeq-8.1$, while the present epoch corresponds to $x=0$. From the left panel in Fig. 2 we see that, if we set $u_{0}=0$, the auxiliary fields, and therefore $\rho_{\rm DE}$, start from zero during RD, and then grow rather steeply during MD and the present DE-dominated era. This is due to the fact that $R=0$ in RD, so with the initial condition $u_{0}=0$ we have $\bar{U}(x)=0$. As we enter in MD, $R$ becomes non-zero and drives the growth of $\bar{U}(x)$. In turn, $\bar{U}(x)$ will eventually drive the growth of $\bar{V}(x)$, through eq. (3.6). We see indeed from the right panel in Fig. 2 that the growth of $\bar{V}(x)$ is delayed compared to the growth of $\bar{U}(x)$. For large $u_{0}$, we see from Fig. 3 that $\bar{U}(x)$ [and therefore $\rho_{\rm DE}(x)$] is non-vanishing already in RD, and it evolves much more slowly. The evolution of the DE density, and therefore the whole cosmological evolution, is then much closer to that of $\Lambda$CDM. This can also be quantified in terms of the DE equation of state $w_{\rm DE}(x)$, defined as usual from the conservation equation ${\rho}^{\prime}_{\rm DE}+3(1+w_{\rm DE})\rho_{\rm DE}=0$. Observe that, since both $\rho_{\rm DE}(x)$ and its derivative ${\rho}^{\prime}_{\rm DE}=d\rho_{\rm DE}/dx$ are positive, necessarily $[1+w_{\rm DE}(x)]$ is negative, i.e. the DE equation of state is on the phantom side. Fig. 4 shows $w_{\rm DE}$ as a function of redshift $z$ [related to $x=\ln a$ by $x=-\log(1+z)$] for $u_{0}=0$ and $u_{0}=200$. We see that, for $u_{0}=200$, $w_{\rm DE}(z)$ differs from the $\Lambda$CDM value $-1$ only by about $1\%$, with $w_{\rm DE}(z=0)\simeq-1.01$, while for $u_{0}=0$ the deviation from $\Lambda$CDM is much larger, and $w_{\rm DE}(z=0)\simeq-1.15$.

3.2 Static solution for the RR model

We next study the static Schwarzschild solution for the RR model, expanding the analysis performed in [68, 39]. In this subsection, for notational simplicity, we denote simply by $U(r)$ and $S(r)$ the solutions denoted above as $U_{\rm static}(r)$ and $S_{\rm static}(r)$. We consider a static spherically-symmetric metric of the form

[TABLE]

Using the notation $A(r)=e^{2\alpha(r)}$ and $B(r)=e^{2\beta(r)}$ and proceeding as in GR, we get two independent equations for $\alpha(r)$ and $\beta(r)$ by taking the combinations $e^{2(\beta-\alpha)}R_{00}+R_{11}$ and $R_{22}$. In the RR model, using eq. (2.16), this gives

[TABLE]

which reduce to their GR counterparts for $m=0$. In the metric (3.11) the equations for the auxiliary fields, eqs. (2.17) and (2.18), become

[TABLE]

The solution of these equations has been briefly discussed in [39], following the detailed analysis for the RT model performed in [68]. In app. A we give a more detailed derivation specific to the RR model, while here we just give the results. We are only interested in the solution at $r\gg r_{S}\equiv 2GM$, which is the regime relevant for studying solar system tests or LLR. For the metric, the solution for $r\gg r_{S}$ (and $mr$ generic) is the same as for the RT model,

[TABLE]

In particular, for $r_{S}\ll r\ll m^{-1}$ we have

[TABLE]

and $B(r)\simeq 1/A(r)$. Observe that this reduces smoothly to the Schwarzschild solution of GR in the limit $m\rightarrow 0$, contrary to what happens in massive gravity where, at the linearized level, there is the vDVZ discontinuity. The solution for the auxiliary fields $U(r)$ and $S(r)$, again in the limit $r\gg r_{S}$, is

[TABLE]

In the region $r_{S}\ll r\ll m^{-1}$ these solutions reduce to

[TABLE]

Note that the field $S$ has dimensions of inverse mass squared, so the combination $m^{2}S$ in eq. (3.20) is dimensionless, and is the one that enters the effective Newton constant (2.20). Observe that, through eq. (2.20), a spatial dependence of $S$ induces a spatial dependence of Newton’s constant. However, the resulting modification of Newton’s constant must be contained in the motion of a test particle in the Schwarzschild metric (3.11), (3.16), (3.17). The latter shows that any deviation from GR vanishes in the limit $mr\rightarrow 0$. This is indeed in agreement with eq. (3.20), that shows that, in the static context that we are considering here (i.e. a Schwarzschild solution that at infinity reduces to flat Minkowski space, rather than being matched to FRW), in the region $r_{S}\ll r\ll m^{-1}$ we have

[TABLE]

so the deviations are suppressed not only by $r_{S}/r\ll 1$ but also by a further factor $(mr)^{2}$ which, for $m\sim H_{0}$ (as determined by the cosmological solution) and $r$ of order, say, of the solar system scale, is totally negligible.

3.3 A simple model for the matching: McVittie’s metric

In the two previous subsections we have studied separately the cosmological solution of the RR model, in which the source is just given by the cosmic fluid, and the static solution, in which the source is just a static central mass and the metric at infinity approaches the flat Minkowski metric. We now wish to put together these solutions, to understand what is the short-distance solution for the auxiliary field $S$ that determines the effective Newton constant, when at large distances the metric approaches FRW. We begin with a somewhat idealized but still useful setting, which is provided by the McVittie metric [70] (see e.g. [71] for review). This metric can be taken as a model of a spacetime sourced by a central mass $M$ plus a cosmic fluid, at distances much greater than the Schwarzschild radius of the central mass. To introduce it let us recall first of all that, if we denote by ${\bf x}$ the comoving coordinates in FRW and by ${\bf r}=a(t){\bf x}$ the physical coordinates, the flat FRW metric can be rewritten as

[TABLE]

where $r=|{\bf r}|$. Physical coordinates are particularly useful for matching this solution to the Schwarzschild metric generated by a central mass $M$. Indeed, in physical coordinates, the McVittie metric takes the very simple form [72]

[TABLE]

where $r_{S}=2GM$. This metric reduces to the FRW solution (3.23) as $M\rightarrow 0$, and to the Schwarzschild solution as $H\rightarrow 0$. Observe also that, in $g_{00}$, the relative weight of the two terms $r_{S}/r$ and $H^{2}r^{2}$ is determined by the combination $r_{c}(t)$, defined by $r_{S}/r_{c}(t)=r_{c}^{2}(t)H^{2}(t)$, i.e.

[TABLE]

For $r\gg r_{c}$ we are in the cosmological regime, where the FRW term dominates, while for $r\ll r_{c}$ the metric is dominated by the central mass.

Note however that, unlike the Schwarzschild metric, which in $r=r_{S}$ only has a coordinate singularity, the McVittie metric has a genuine physical singularity at $r=r_{S}$. Indeed, if one computes the energy-momentum tensor that sources this metric, one finds a fluid with a pressure $p(t,r)$ that, at $r\gg r_{S}$, correctly reduces to the cosmological fluid corresponding to the chosen time dependence of $H(t)$, but that diverges as $r\rightarrow r_{S}$. For instance, when matching to a matter-dominated universe at cosmological scales, one finds that the relation between pressure and energy density is [72]

[TABLE]

This correctly reduces to the result $p=0$ for non-relativistic matter at $r\gg r_{S}$, but diverges as $r\rightarrow r_{S}$. Thus, the McVittie metric becomes pathological, and is no longer a useful model for a mass embedded in an expanding universe, when $r$ is close to $r_{S}$. However, this is only relevant if we want to describe black holes or compact objects in an expanding universe. For a central object of size much larger than its Schwarzschild radius, the metric in the exterior region can be very well approximated by eq. (3.24), and the pressure $p(t,r\gg r_{S})$ is an excellent approximation to that of the desired cosmic fluid.

Thus, we can use the McVittie metric as a first simple toy model to study the effect of the cosmic expansion on the Earth-Moon system (with the Earth playing the role of the central mass $M$ and the Moon of a test mass), and in particular to study whether, in this metric, the field $S$ of the RR model, that determines the effective Newton constant, retains a time dependence at small distances. Of course, the McVittie metric is a solution of the Einstein equations of GR, and not of those of the RR model (as it is clear from the fact that it does not even depend on the parameter $m$). This means that, when $mr$ is no longer negligible, in the RR model we should expect deviations from it. In this section we will then consider a slightly more general metric

[TABLE]

with a more general function $\Phi(r)$. The McVittie metric corresponds to $2\Phi(r)=r_{S}/r$. The simplest possibility is that, in the RR model, the function $\Phi(r)$ will receive corrections $1+{\cal O}(m^{2}r^{2})$. More generally, the various occurrences of $\Phi$ in eq. (3.27) could be modified in different ways as $mr$ becomes non-negligible. However, as a first toy model, eq. (3.27) will already be sufficient to understand the behavior of the solutions of the equations $\Box U=-R$ and $\Box S=-U$ on a background that interpolates between a Schwarzschild metric at small scales and FRW at large scales. Furthermore, from the analysis in the similar setting performed in app. A, we understand that $\Phi(r)$ will be suppressed by a factor $r_{S}/r$ even at cosmological scales, i.e. when $mr$ is not small, see eq. (A.23) and is therefore extremely small at cosmological scales. This also means that we do not need knowledge of the exact way in which the solution of the RR model differs from the McVittie solution of GR as $mr$ becomes of order one.

Actually, eq. (3.27) must be understood as a first step in an iterative procedure in which $\Phi(r)$ is more generally replaced by a function that also has an adiabatic dependence on time, $\Phi=\Phi(t,r)$. Indeed, if the analysis in the metric (3.27) should lead to the conclusion that $S$ depends on time even at short scales, as indeed we will show below, it would follows that the dynamics at short distances is governed by $G_{\rm eff}(t)$, rather than by $G$. In turn, this means that the occurrencies of $G$ that are implicitly present in $\Phi(r)$, e.g. the factor $G$ in $\Phi(r)=r_{S}/(2r)=GM/r$, must be replaced by $G_{\rm eff}(t)$, and correspondingly $\Phi(r)$ becomes a function $\Phi(t,r)$ with a dependence on time. Since in any case $G_{\rm eff}(t)/G$ is very close to one, an iterative approach of this kind, in which we start at first with a function $\Phi(r)$ in eq. (3.27), is appropriate. We will see explicitly how this iterative procedure works in the discussion below eq. (3.37).

One might argue that this setting is however still oversimplified because, to move from the Earth-Moon system scale to cosmological scales we have to go through a hierarchy of embeddings, from Solar System to the Galaxy, Local Group, Clusters and Superclusters, and one might wonder whether each such embedding could provide some form of screening, so that the matching of the local Schwarzschild metric to the cosmological FRW solution will be more complicated [71]. We will address this question in section 3.4 and, for the moment, we just proceed with the exercise in the McVittie-like metric (3.27). The use of a generic function $\Phi(r)$, rather than the explicit expression $r_{S}/(2r)$ corresponding to the McVittie solution of GR, will also allow us to have a better understanding of the structure of the equations, and to more easily make contact with the analysis of section 3.4.

We will assume that, in the region of interest (such as $r\gg r_{S}$ for a single central mass with Schwarzschild radius $r_{S}$), we have

[TABLE]

and also

[TABLE]

This is the case in particular for the McVittie metric when $r\gg r_{S}$, since in that case $\Phi(r)=r_{S}/(2r)$ and $r\partial_{r}\Phi(r)=-r_{S}/(2r)$. In the metric (3.27), neglecting terms ${\cal O}(\Phi^{2})$, the equations $\Box U=-R$ and $\Box S=-U$ become

[TABLE]

(where $\dot{H}=\partial_{t}H$), and

[TABLE]

In the RR model the Ricci scalar $R$, that appears on the right-hand side of eq. (3.30), can be determined by taking the trace of eq. (2.16), which gives

[TABLE]

For a single fluid $T=-\rho+3p$, so for a cosmological fluid with energy density and pressures $\rho_{\rm cosmo}(t)$ and $p_{\rm cosmo}(t)$, respectively, plus a central mass $M$,

[TABLE]

We now study whether eqs. (3.30) and (3.31) have an approximate solution

[TABLE]

where $U_{\rm cosmo}(t)$ is the cosmological solution discussed in section 3.1 and $U_{\rm static}(r)$ is the static Schwarzschild solution of section 3.2. The details of the analysis are given in App. B, while here we give the main results.

On the left-hand side of eq. (3.32) we can make use of the fact that $m^{2}S_{\rm static}(r)$ is always very small compared to one, as long as $r\gg r_{S}$, as we see from eq. (3.20), so we can replace $[1-(1/3)m^{2}S]$ by $[1-(1/3)m^{2}S_{\rm cosmo}(t)]$. We then introduce the notation

[TABLE]

Then, the analysis in App. B shows that, on the ansatz (3.34), the equation $\Box U=-R$ takes the form

[TABLE]

where the question mark stresses that we are just checking if the ansatz is correct, and

[TABLE]

By definition, $U_{\rm cosmo}(t)$ and $S_{\rm cosmo}(t)$ solve the equation

[TABLE]

while $U_{\rm static}(r)$ and $S_{\rm static}(r)$ solve

[TABLE]

Comparing with eq. (3.37) we see that the ansatz (3.34) would indeed be a solution of the equations, were it not for the term $G_{\rm eff}(t)/G$ on the right-hand side of eq. (3.37). Using, as in section 3.1, $\gamma=m^{2}/(9H_{0}^{2})$ and $V(t)=H_{0}^{2}S_{\rm cosmo}(t)$, we have

[TABLE]

Parametrically, this is not suppressed by any large factor. However, from the numerical solution of the cosmological equations discussed in section 3.1, we know that $G_{\rm eff}(t)/G$ is always very close to one, in the past and present cosmological evolution. For instance, for $u_{0}=0$ we have $\gamma\simeq 9.126\times 10^{-3}$ and $V(t_{0})\simeq 2.14$, so $3\gamma V(t_{0})\simeq 0.06$. For large values of $u_{0}$ this quantity decreases further. For instance, for $u_{0}=200$ we have $\gamma\simeq 5.907\times 10^{-5}$, $V(t_{0})\simeq 30.45$ and $3\gamma V(t_{0})\simeq 0.005$. Thus, we can devise a perturbative strategy in the small quantity $3\gamma V(t)$. To zero-th order we just neglect it, and then $G_{\rm eff}(t)=G$. Of course, to this order we have lost the effect that we are looking for, namely whether the effective Newton constant is time dependent. We must then go to the next order. To this purpose we start from the expression

[TABLE]

and then, to compute $S$ on the right-hand side, we set $G_{\rm eff}=G$ in eq. (3.37). In this case, we see that the ansatz (3.34) indeed solves (within the approximations discussed) $\Box U=-R$, and precisely the same analysis applied to $\Box S=-U$ shows that indeed $U=U_{\rm cosmo}(t)+U_{\rm static}(r)$ and $S=S_{\rm cosmo}(t)+S_{\rm static}(r)$ is an approximate solution of both $\Box U=-R$ and $\Box S=-U$. Plugging this solution for $S$ in eq. (3.42) we get $G_{\rm eff}/G$ to first order in $3\gamma V(t)$,

[TABLE]

which is indeed the expression used in [1] and in section 2.3. We now see more formally how this emerges, at least within the McVittie model, and under what approximations.

One could in principle iterate the process. It is clear that, at next order, the result will basically amount to replacing the static solutions $U_{\rm static}(r)$ and $S_{\rm static}(r)$ with a quasi-static solution that is obtained by replacing $G$ with the $G_{\rm eff}(t)$ computed at first order. The physics behind this iterative procedure is quite clear. If Newton’s constant has a slow time dependence on cosmological time-scales, on the time-scale ${\cal O}(10)$ yr probed by LLR the change in the Newton constant has the effect of an adiabatic perturbation, and the corresponding Schwarzschild solutions at $r\gg r_{S}$ is simply obtained by replacing $G$ with $G_{\rm eff}(t)$.

The conclusion of this section is that, within the modelisation of the problem provided by the McVittie solution, it is indeed true that, even at the Earth-Moon scale probed by LLR, Newton’s constant is replaced by an effective time-dependent Newton’s constant, just as on cosmological scales. To a precision of order of a few percent or better, the corresponding expression is given by eq. (3.43).111111Indeed, the accuracy is probably much better since, at next iterative order, we expect that the only change would be given by the replacement of $U_{\rm static}(r)$ and $S_{\rm static}(r)$ with the corresponding quasi-static solutions. Since, in any case, in $S=S_{\rm cosmic}(t)+S_{\rm static}(r)$ we have seen that $S_{\rm cosmic}(t)$ is numerically much larger than $S_{\rm static}(r)$, the result of successive iterations will have very little effect on $G_{\rm eff}(t)/G$. As we found in section 2.3, in this case the LLR bound is violated by a factor that, for $u_{0}=0$, is ${\cal O}(100)$, and even for $u_{0}=200$ is still violated by about one order of magnitude, so the level of accuracy of the first-order solution (3.43) is perfectly adequate and shows that, within this McVittie modelisation, the RR model is ruled out by LLR, unless we go to very large values of $u_{0}$, u_{0}\,\raise 1.72218pt\hbox{>}\kern-6.99997pt\lower 2.6694pt\hbox{\sim}\,{\cal O}(2000).

3.4 A more general analysis: perturbed FRW

The next question is whether the result found above is specific to our modelization of the problem. After all the McVittie metric, describing a single central mass in an expanding FRW universe, might look like a very simplified setting, compared to the actual setting in which the mass $M$ (the Earth for LLR) is embedded into a hierarchy of structures such as the solar system, the Galaxy, the Local Group, clusters and superclusters, before reaching the scales where the Universe can be approximated by a uniform FRW background. Before drawing our conclusions from a computation performed in the McVitties metric, we should investigate whether the presence of these non-linear structures can somehow provide a screening mechanism. This requires to study the equations $\Box U=-R$ and $\Box S=-U$ on a metric that gives a realistic description of the Universe from the Earth-Moon scale to the cosmological scales. The most natural choice is to consider a perturbed FRW spacetime (3.44). In the absence of anisotropic stress121212In the RR model the nonlocal term generates a small anisotropic stress, which however is completely negligible numerically, see Fig. 8 of ref. [40]. we have $\Psi=-\Phi$, so

[TABLE]

To make contact with the analysis of the previous subsection one could also re-express this metric using physical coordinates ${\bf r}=a(t){\bf x}$, which leads to

[TABLE]

For McVittie metric, physical coordinates were particularly transparent for showing how the Schwarzschild metric merges with the FRW metric. However, in this section we will use comoving coordinates, since the corresponding equations turn out to be somewhat simpler.

In recent years there has been much discussion on whether nonlinear structures can imprint a significant backreaction on the metric, in particular in the attempt of explaining the accelerated expansion of the Universe (see e.g. [73] for review). The issue is controversial, but there is by now a broad consensus on the fact that no large backreaction appears if one uses the Newtonian gauge (2.8), both from analytic arguments (see e.g. [74, 75]) and from $N$-body simulations [76, 77]. The metric (2.8) [or, neglecting anisotropic stresses, the metric (3.44)] gives an excellent approximation to the actual metric of the Universe (in the scalar-perturbation sector), and $|\Phi|$ never exceeds a value of order $10^{-4}$, except in strong-field regions, such as the vicinity of black holes or neutron stars. In particular, for a virialized system with characteristic velocity $v$, the virial theorem gives $\Phi\sim v^{2}$. However, the spatial derivatives are enhanced. For the first derivative, from the Euler equation we get $\partial_{i}\Phi\sim{\cal H}v$ [where ${\cal H}=(1/a)\partial_{\eta}a$, $\eta$ is conformal time, and $\partial_{i}$ is the derivative with respect to comoving coordinates] so, in the natural units provided by the Hubble parameter, $\partial_{i}\Phi$ is suppressed only by one power of $v$. Furthermore, from Poisson equation, we have {\mbox{\boldmath\nabla}}^{2}\Phi\sim{\cal H}^{2}(\delta\rho/\rho). In non-linear structures the overdensity $\delta\rho/\rho$ can be extremely large, say ${\cal O}(10^{30})$. Still, even in these cases $|\Phi|$ is small and eq. (3.44) provides a good description of the metric. We will then take the perturbed FRW metric (3.44) as a good description from cosmological scales down to the scale of the Earth-Moon system probed by LLR.

Of course, on the length-scales and time-scales relevant for LLR, the effect of the explicit occurrence of $a(t)$ in eq. (3.44) [or, equivalently, the effect of the $Hr$ and $H^{2}r^{2}$ terms in eq. (3.45)] is in any case completely negligible compared to the effect of the static or quasi-static potential $\Phi$ generated by the local mass distribution and, to excellent accuracy, on these scales we might as well use a perturbed Minkowski metric, setting $a(t)=1$. Furthermore, as we recalled in section 3.2, in the RR model there is no vDVZ discontinuity and the effects associated to the mass $m$ go smoothly to zero on scales $r$ such that $mr\ll 1$ (in contrast, e.g. to massive gravity, where, at the linearized level, there is the vDVZ discontinuity). Thus, one might be tempted to conclude that, in the RR model, the motion of a test mass (such as the Moon) in the field of the Earth is in practice indistinguishable from that in GR. As pointed out in [1], this is not necessarily correct if Newton’s constant (that implicitly enters the metric through $\Phi$, e.g. through the expression $\Phi(r)=GM/r$ in the case of a single massive body) retains a time dependence down to such scales, because the extraordinary accuracy of the LLR measurement allows us to detect a rate of change of Newton’s constant comparable to $H_{0}$ [and in fact even three order of magnitude smaller, see eq. (2.2)]. In other words, the potentially dangerous time dependence of the metric at short scales does not come from the explicit occurrences of $a(t)$ in eq. (3.44), or of $H(t)$ in eq. (3.45), but rather from a possible implicit time dependence that could be present in $\Phi$ through a time-dependent Newton constant. In turn, whether Newton’s constant evolves with time depends on whether the auxiliary fields $U$ and $S$, solutions of $\Box U=-R$ and $\Box S=-U$, retain a time dependence on these small scales, due to the fact that they are time dependent on cosmological scales. As we have seen, this is precisely what happens in the McVittie metric.

We then use a perturbed FRW metric to study the solution of the equations $\Box U=-R$ and $\Box S=-U$. As in the previous section, we give here the main results, deferring to App. C some technicalities. To first order in $\Phi$, for a generic function $f$, in the metric (3.44) we have

[TABLE]

where $\nabla$ is the gradient with respect to the comoving coordinates. We could now deal with the equation $\Box U=-R$ just as we did for the McVittie metric, looking for an approximate solution displaying a separation of variables.131313In section  3.3 we have shown the existence of a solution with an approximate separation of variables between cosmic time $t$ and physical coordinates ${\bf r}$. We could as well have used cosmic time and comoving coordinates ${\bf x}$. Indeed, for a differentiable function $f(t,{\bf x})$ we have $(\partial f/\partial x^{i})_{t}=a(\partial f/\partial r^{i})_{t}$ and

$\left(\frac{\partial f}{\partial t}\right)_{x}=\left(\frac{\partial f}{\partial t}\right)_{r}+Hr^{i}\left(\frac{\partial f}{\partial r^{i}}\right)_{t}\,.$

(3.47) Taking $f(t,r)=U_{\rm cosmo}(t)+U_{\rm static}(r)$, we have $(\partial f/\partial t)_{r}\sim HU_{\rm cosmo}$, and $U_{\rm cosmo}={\cal O}(1-10)$ in the recent cosmological epoch. In contrast, in a McVittie-like metric, $(\partial f/\partial r)_{t}\sim(1/r)U_{\rm static}$ and $U_{\rm static}\sim r_{S}/r\ll 1$ (one can also express this by observing that spatial derivatives vanish on the FRW background, so they pick an extra factor of order $\Phi$ with respect to time derivatives, and $\Phi\sim r_{S}/r\ll 1$). Then, the second term on the right-hand side of eq. (3.47) is negligible with respect to the first (and similarly for $S$). However, in the spirit of the use of a perturbed FRW metric, it is more natural to look for a perturbative solution of the form

[TABLE]

where $\delta U$ is a perturbation with respect to $U_{c}$ and similarly for $S$. [In this section, to keep the equations compact, we denote the cosmological background solutions $U_{\rm cosmo}(t)$ and $S_{\rm cosmo}(t)$ by $U_{c}(t)$ and $S_{c}(t)$, respectively]. Of course, we already know from the analysis in [58, 40] that a perturbative approach is consistent at cosmological scales, where the perturbations of the matter density are small. Here however we wish to see if such a perturbative solution is still consistent everywhere, including the ‘short distance’ scales where there can be large density perturbations, and therefore without linearizing with respect to the matter density contrast $\delta_{M}\equiv\delta\rho_{M}/\rho_{M}$, as in [58, 40]. The key to this more general treatment is that, as we have discussed, $|\Phi|$ still remains small even when its Laplacian, which is related to the density perturbations, is large.

We begin by discussing the equation $\square U=-R$. We drop the ${\cal O}(\Phi)$ corrections (as we did in the analysis for the McVittie metric), and we use the perturbative ansatz (3.48), (3.49). In addition, we assume a quasi-static approximation

[TABLE]

where $\delta\dot{U}\equiv\partial_{t}(\delta U)$, $\delta\dot{S}\equiv\partial_{t}(\delta S)$, and we also assume

[TABLE]

that will be checked a posteriori. Using these approximations, we show in App. C that the equation $\Box U=-R$ reduces to

[TABLE]

where $\rho_{b}$ is the background energy density of the cosmological fluid. We now make use of the fact that, in the sub-horizon limit, $\Phi$ satisfies the Poisson equation of the RR model,141414The Poisson equation is a combination of the $00$ and of the divergence of $0i$ components of the modified Einstein equation of the RR model, in the quasi-static and sub-horizon limit, and does not rely on the use of linear cosmological perturbation theory. Therefore, just as in GR, also for the RR model it holds even for large density contrasts.

[TABLE]

Thus, eq. (3.53) can be rewritten as

[TABLE]

Next observe that, for modes with wavenumber $k$, the left-hand side is of order

[TABLE]

while, using the quasi-static approximation,

[TABLE]

and $U_{c}$ is ${\cal O}(1-10)$, so in momentum space

[TABLE]

Inside the horizon $H\ll k/a$, so the right-hand side of eq. (3.55) is parametrically smaller than the left-hand side, and we must have

[TABLE]

Thus, the smallness of $\Phi$ ensures the smallness of $\delta U$ with respect to $U_{c}$, which in contrast is ${\cal O}(1-10)$, and validates the use of the perturbative ansatz (3.48).

To conclude our argument we need to check the consistency of the assumption (3.52). This is obtained for a study of the equation $\Box S=-U$, and again the details are given in App. C. Note that eq. (3.52) also implies that $|\delta S|\ll|S_{c}|$, since $|\delta U|\ll U_{c}\sim H_{0}^{2}S_{c}\sim m^{2}S_{c}$ today, and therefore confirms the validity of the perturbative approach also for $S$.

Therefore, the perturbative ansatz gives a consistent solution to the equations for the auxiliary fields in terms of a generic metric perturbation $\Phi$. In turn, the latter can be determined from the modified Einstein equation (2.16), providing a fully consistent perturbative solution of the equations of the RR model. Furthermore, the perturbative ansatz obeys sensible initial conditions, matching to the linear perturbation theory solution at early times. Indeed, before the epoch of non-linear structure formation $\rho$ is very close to $\rho_{b}$ and, from eq. (3.53), $\delta U$ is basically zero. In other words what we have found is that, before structure formation, when the homogeneous FRW metric holds to great accuracy everywhere, the fields $U$ and $S$ evolve with time according to the cosmological solution. More precisely, as we recalled in section 3.1, deep in RD the cosmological solutions for the auxiliary fields is basically constant in time, but the evolution of $U$ becomes non-negligible as we enter in MD, even before structures become non-linear. As structures form and become non-linear, the solutions for $U$ and $S$ remains a small perturbations of the cosmological solution, of the form (3.48), (3.49), as a reflection of the fact that the metric perturbations never become large (except, of course, in the very limited regions around compact bodies). Thus, structure formation does not stop the time evolution of the auxiliary fields.

It is interesting to compare the results of this section with those obtained in section 3.3 in the McVittie metric. The result $\delta U=2\Phi$ that we have found inside the horizon in the very general setting of a perturbed FRW spacetime is confirmed in the special case of the McVittie-like solution, where $\delta U$ corresponds to $U_{\rm static}$. Indeed inside the horizon, i.e. for $mr\ll 1$, we found in eq. (3.21) that $U_{\rm static}=r_{S}/r$, while, again for $mr\ll 1$, the potential $\Phi$ is given by $r_{S}/(2r)$, see eq. (A.29). Thus, the results of the two approaches are fully consistent.

The result $\delta U=2\Phi$ is also consistent with that found using linear cosmological perturbation theory for the RR model. Indeed, in that case, for modes well inside the horizon, in the linearization approximation one finds $\delta U=2(\Psi+2\Phi)$ [see eq. (4.10) of [40], neglecting the time derivatives as appropriate to the quasi-static approximation and retaining only the leading term in the the limit $\hat{k}\equiv k/(aH)\gg 1$]. Together with $\Psi=-\Phi$ (which is due to the fact that in the RR model anisotropic stresses are negligible, see again [40]), this gives again $\delta U=2\Phi$. Once again, the advantage of the present discussion, with respect to the linearized treatment of [40], is that here we have not assumed that the density perturbations are small, but only that $|\Phi|$ is small, a condition that, in the absence of strong backreaction effects, includes the situation in which there is a full hierarchy of structures between the solar system scale and the cosmological scales.

To sum up, we have found a consistent solution of the form (3.48), (3.49), where the perturbations $\delta U$ and $\delta S$ induced by the perturbation of FRW (and therefore by non-linear structures, as long as one accepts that there is no large backreaction and non-linear structures are well described by a perturbed FRW metric) remain small. As a consequence, the solutions for $U$ and $S$ retain their time dependence down to short scales, and even on LLR scales we have $G_{\rm eff}\simeq G/[1-(1/3)m^{2}S_{c}(t)]$.

Finally, one might ask what happens if one tries to force an exactly static solution at the Earth-Moon scale. In GR this can be done, and the classic example is the Einstein-Straus construction, in which one has again a central mass $M$ embedded in FRW, as in the McVittie case, but, rather than finding a metric that interpolates smoothly between the two, one just matches the two solutions on a boundary surface $\Sigma$ (see e.g. [71, 78] for review). In the interior the metric is taken to be exactly the static Schwarzschild metric generated by the mass $M$, while in the exterior it is given by a FRW solution with energy density $\rho$. The two metrics are then matched by requiring that the induced metric on $\Sigma$ agrees on the two sides. This fixes the matching radius $r_{0}$, that, with respect to the Schwarzschild coordinates of the interior (that correspond to the physical FRW coordinate of the exterior region) turns out to be given by $M=(4/3)\pi r_{0}^{3}\rho$, where $\rho$ is the energy density in FRW. In this case the solution for the metric is exactly static in the interior region. Nevertheless, if one studies the propagation of a scalar field obeying the equation $\Box\phi=0$ in this metric, one finds that the solution for the field in the inner region is time dependent [78]. Technically, this comes from the fact that we must impose a matching condition for the scalar field at the surface $\Sigma$, and in this way the field inherits a time dependence even in the inner region. Similarly, in the RR model, because of the matching at $\Sigma$, the fields $U$ and $S$ will depend on time in the interior, thereby inducing a time-dependence in the Newton constant. This means that, in the RR model, there will be no exact Einstein-Straus solution, since even in the interior there will be a time-dependence in the metric induced by $G_{\rm eff}(t)$.

4 No perfect screening for free

We next examine how the Deser-Woodard model behaves with respect to LLR. In ref. [16] a mechanism (referred as ‘perfect screening for free’ in [17]) was proposed, based on the observation that the scalar d’Alembertian $\Box$ has different signs when acting on functions of time and on functions of space, as it is clear from the flat-space limit \Box\rightarrow-\partial_{t}^{2}+{\mbox{\boldmath\nabla}}^{2}. As we saw in eq. (1.2), the DW model can be written in terms of a function $f(X)$, where

[TABLE]

(Note that, in the notations that we use for the RR model, $X=-U$; in the following we will use both $X$ and $U$). We first observe that $R$ is positive both in a cosmological setting, where the main contribution to the energy-momentum tensor comes from the cosmological fluid, and in the local regions (such as the solar system, or the Earth-Moon system), where the main contribution comes from the mass of the non-linear structures, because in any case it is dominated by the positive energy density $\rho$ [see e.g. eq. (B.7) for the full expression in the RR model]. Then, in ref. [16, 17] it is argued that in cosmology (where the time derivative dominates) $X\simeq\frac{1}{-\partial_{t}^{2}}R$ is negative because of the minus sign in $-\partial_{t}^{2}$, while it is positive in the regime dominated by structure formation, where the spatial derivatives dominate and X\simeq\frac{1}{{\mbox{\boldmath\nabla}}^{2}}R. Deser and Woodard then complete the definition of their model by stating that the function $f(X)$ is given by eq. (1.3) for $X<0$, and by $f(X)=0$ for $X>0$. According to the authors, this would automatically set to zero any effect of the non-local term in the solar system, and should allow the model to pass the solar system (and LLR) tests. Ref. [17] then discusses why, in the author’s view, such a choice of $f(X)$ is not an unnatural fine tuning.

First of all we observe that, if this mechanism actually worked, we could extend it trivially to the RR model. Indeed, upon integration by parts, $m^{2}R\Box^{-2}R$ can be written as $m^{2}(\Box^{-1}R)^{2}$. We could then define the RR model writing $m^{2}f(X)$, with $f(X)=X^{2}$ for $X<0$ and $f(X)=0$ for $X>0$.

However, this mechanism does not work, neither for the RR model nor for the DW model. The first reason is that $1/\Box$ is not an algebraic quantity but an integral operator, and one cannot determine the sign of $X$ as done above. Equation (4.1) is equivalent to the differential equation $\Box X=R$. In a regime where the time derivative dominates this is equivalent to $-\partial_{t}^{2}X=R$, i.e., in our notation, $\partial_{t}^{2}U=R$ [in flat space; the corresponding equation in FRW is eq. (3.2)]. Note that the equation for $U$ is the same in the DW and RR models, since in both cases it simply follows from the definition $\Box U=-R$. So, first of all, the fact that $R$ is positive only provides a statement about the sign of the second derivative of $U$ (in flat space; or of the corresponding differential operator in FRW), and not directly about $U$ itself. To integrate $\partial_{t}^{2}U=R$, or eq. (3.2), we also need the initial conditions on $U$. In the DW model one sets $U(t_{i})=\dot{U}(t_{i})=0$ at an initial time $t_{i}$, chosen deep in RD (see e.g. eq. (37) of [17]). In this case, the integration of the differential equation indeed confirms that $U>0$, i.e. $X<0$.

Consider now the regime dominated by spatial derivatives. In that case, in the flat-space limit, {\mbox{\boldmath\nabla}}^{2}X\simeq R, i.e. {\mbox{\boldmath\nabla}}^{2}U\simeq-R. The argument of ref. [16, 17] amounts to saying that, since $R$ is positive, the solution $X$ of {\mbox{\boldmath\nabla}}^{2}X=R is also positive. In general, this is false. For instance,

[TABLE]

which shows that the positive function $1/r$ has a Laplacian given by a Dirac delta with negative coefficient. In other words, the Green’s function of the Laplacian is

[TABLE]

with a crucial minus sign. At short distances, both in the RR and DW models, the dominant contribution to the Ricci scalar $R$ comes from the energy density $\rho$ and can be determined by taking the trace of the Einstein equations, which gives $R\simeq 8\pi G\rho$.151515The fact that $R$ at short scales reduces to the GR expression $8\pi G\rho$ is necessary in any model that aims at modifying GR at the horizon scale without spoiling solar system and laboratory tests, unless there is a vDVZ discontinuity and a corresponding Vainshtein mechanism. This source term is therefore localized in space (e.g. it can be approximated with a Dirac delta for a single mass such as the Earth, at $r\gg r_{S}$). In this case, the solution of {\mbox{\boldmath\nabla}}^{2}X=8\pi G\rho is

[TABLE]

A positive $\rho$ therefore generates a negative $X$, i.e. a positive $U$, just as in the regime dominated by the time derivative.161616Observe that the requirement that $\rho$ has compact support is important because it ensures the convergence of the convolution with the Green’s function. The Ricci scalar also contains non-compact terms that depend on the auxiliary fields themselves (and that are specific to the theory considered, DW or RR; see eq. (B.7) for the expression in the RR model). Their contribution can be included iteratively, using first the zero-th order result for the auxiliary fields obtained including in $R$ only the dominant contribution coming from $\rho$, and using them to compute the source term $R$ to the next iterative order. Already to the first iterative order the source term $R$ no longer has compact support, and the solution cannot be obtained with the simple Green’s function (4.3). There are however other techniques to deal with the problem. Technically, the issue becomes very similar to the one encountered when one computes the gravitational-wave production from a localized system, such as a compact binary system, using a post-Newtonian/post-Minkowskian expansion: in that case, the gravitational waves produced at leading order act themselves as non-compact sources for gravitational-wave production at the next order (see e.g. section 5.3 of [79] for a pedagogical review). In any case, the point is that these corrections are calculable and finite. At short scales, where the dynamics is dominated by the matter density $\rho$, they are just a small correction to the leading term and they cannot affect the sign of the solution for $U$ computed to leading order. On cosmological scales the Ricci scalar $R$ is no longer dominated by $\rho$, but rather by the contribution of the auxiliary fields that are responsible for the effective dark energy, and $U$ is no longer necessarily positive. Indeed, in the RR model it oscillates as $\cos(mr)$, see eq. (3.19).

In terms of $U=-X$,

[TABLE]

Taking for instance $\rho({\bf x})=M\delta^{(3)}({\bf x})$, we get $U(r)=r_{S}/r$, as we already found in eq. (3.21). In conclusion, the static solution for the auxiliary field has the same sign as the cosmological solution, so that $U$ is always positive (i.e. $X$ is always negative), contrary to what is needed in the screening mechanism proposed in [16, 17].

Actually there is a second, more subtle, reason why this ‘perfect screening for free’ mechanism does not work. As we have seen in sections 3.3 and 3.4, once we look for a solution valid across the whole range of scales, from the short scales relevant for LLR and solar system experiments, up to the cosmological scales, we find that the solution at short scales is actually well approximated by the sum of the static solution and the cosmological solution, $U(t,r)\simeq U_{\rm cosmo}(t)+U_{\rm static}(r)$. Note, in particular, that the analysis of section 3.4 based on the perturbed FRW metric can be repeated with basically no modification for the DW model. Indeed the equation to be studied, $\Box U=-R$, is the same for the RR and DW model, and in both cases one just needs to assume the validity of a perturbed FRW metric. What changes between the two models is the explicit expression of the Ricci scalar in terms of the auxiliary fields of the theory. However, this just means that $U_{\rm cosmo}(t)$ will be the respective cosmological solution in each model.171717Furthermore, in the DW model, on cosmological scales the explicit expression of $R$ as a function of time is the same as in $\Lambda$CDM because, by construction, the DW model has the same background evolution as $\Lambda$CDM, while in the RR model the evolution turns out to differ from that in $\Lambda$CDM only at the level of about $10\%$, so in practice, on cosmological scales, we are integrating the equation $\Box U=-R$ with a very similar source $R(t)$ while, at short scales, in both the DW and RR model $R$ reduces to the GR expression. Now, numerically, for the Earth-Moon system $U_{\rm static}(r)\simeq r_{S}/r$ where $r_{S}$ is the Schwarzschild radius of the Earth, about $1~{}{\rm cm}$, while $r$ is the Earth-Moon distance. Thus, $r_{S}/r\sim 10^{-11}$ is a minuscule quantity, while $U_{\rm cosmo}(t)={\cal O}(1)$ so, numerically, the solution for $U$ even on the scale of the Earth-Moon system is completely dominated by the cosmological solution.

In conclusion, first of all it is not true that, in the static solution, $U$ is negative; the static solution has the same sign as the cosmological solution. And, furthermore, the sign of the static solution is irrelevant, because even at short scales $U$ is dominated by the cosmological solution. This means that the screening mechanism suggested in [16, 17], based on setting $f(X)=0$ for $X>0$, does not work: $X=-U$ is always negative, both at cosmological scales and at short scales. As a consequence, the time-dependence of the effective Newton constant obtained on cosmological scales holds down to the LLR scale. In the DW model the latter is given by eq. (2.26). This quantity has been computed numerically in [19, 20] and one finds that, today,

[TABLE]

where $z$ is the redshift, see Fig. 4 of [20]. In terms of the derivative with respect to cosmic time this means that, in the present cosmological epoch,

[TABLE]

with a proportionality coefficient of order one. Thus, the DW model violates the LLR bound (2.2) by about three orders of magnitude and, just as the RR model, is ruled out.

5 Conclusions

In this paper we have elaborated on the suggestion in ref. [1], that the limit on the time variation of Newton’s constant from Lunar Laser Ranging could be in conflict with the RR model. The main question was whether one could extrapolate the result for $G_{\rm eff}(t)$ found on cosmological scales down to the scales relevant for Lunar Laser Ranging. The analysis that we have presented shows that this is indeed the case, unless some strong backreaction effect from nonlinear structures comes into play and provides a screening mechanism, a hypothesis that seems unlikely in view of much recent work on the backreaction problem. We have found that this indeed rules out the RR model [unless the parameter $u_{0}$ is very large u_{0}\,\raise 1.72218pt\hbox{>}\kern-6.99997pt\lower 2.6694pt\hbox{\sim}\,{\cal O}(2000), which seems quite unnatural even from the point of view of an earlier inflationary phase]. However the RT model, defined by eq. (1.7), passes this test because in this case $G_{\rm eff}$ reduces to $G$ on small scales, see eq. (2.15). We have further examined the situation in the Deser-Woodard model and we found that it suffers from the same problem as the RR model. The screening mechanism that was proposed for this model, based on the idea that $\Box^{-1}R$ changes sign from the cosmological regime to the short-scale regime, turns out to be technically incorrect, and even the Deser-Woodard model, with the form of the distortion function chosen so to mimic $\Lambda$CDM, is ruled out by Lunar Laser Ranging (at least, unless it can be supplemented by a different screening mechanism).

As we briefly reviewed in the Introduction (see [33] for more details), in the last decades have been proposed several variants of the idea of modifying phenomenologically the quantum effective action of gravity, through the addition of non-local terms relevant in the infrared. At present, the RT model appears to be the only viable model that has passed all the tests performed, namely: (1) At the level of background cosmological evolution it provides a viable FRW solution, with the nonlocal term acting as an effective dark energy, thereby providing accelerated expansion without a cosmological constant. (2) Cosmological perturbations are well behaved, both in the scalar and in the tensor sector. (3) Implementing the model in a Boltzmann code and comparing with observations, one finds that the model fits CMB, BAO, SNe, structure formation data and local $H_{0}$ measurements at a level statistically equivalent to $\Lambda$CDM. (4) Gravitational waves propagate at the speed of light, complying with the limit from GW170817/GRB 170817A. (5) Solar system tests, and in particular the limit on the time variation of Newton’s constant from Lunar Laser Ranging, are satisfied.181818Note also that there is no problem with ghosts nor extra degrees of freedom, once one properly interprets these nonlocal models as derived from a quantum effective action. See in particular section 2.4 of [33], or sections 5-6 of [32], for detailed discussions.

The fact that these requirements are so stringent, to the extent that they appear to select a single model, is quite interesting, and gives potentially important hints for attempts at deriving these nonlocal terms from a fundamental theory, a step that is eventually essential to put these ideas on firmer grounds.

Acknowledgments. We thank Alex Barreira and Baojiu Li for very useful discussions, already several years ago, that eventually led us to fully understand the importance of their observation on LLR. We also thank Giulia Cusin, Yves Dirian and Stefano Foffa for useful discussions and comments on the manuscript. The work of E.B., A.Fi. and M.M. is supported by the Fonds National Suisse and by the SwissMap National Center for Competence in Research. The work of A.Fr. was supported by a Swiss Government Excellence Scholarship and is currently supported from ERC Advanced Grant GravBHs-692951 and MEC grant FPA2016-76005-C2-2-P.

Appendix A Static solution for the RR model

In this appendix we give a more detailed derivation of the results presented in section 3.2. To find the solution we follow the strategy used in [68] for the RT model: we first compute the solution in the limit $mr\ll 1$, obtaining the leading term (and, for the metric, the next-to-leading term) of the solution in the regime $r\ll m^{-1}$; we then compute it in a Newtonian expansion around flat space, which gives a solution valid for $r\gg r_{S}$, and we finally match the two solutions in the overlap region $r_{S}\ll r\ll m^{-1}$, which will allow us to fix all coefficients of the solutions.

A.1 Solution for $mr\ll 1$

Let us begin with the expansion in powers of $mr$. To zero-th order, eqs. (3.12) and (3.13) reduce to their standard GR form, whose solution is given by

[TABLE]

Plugging these expressions into eq. (3.14) we get

[TABLE]

so, at $r\gg r_{S}$,

[TABLE]

Until now, everything is just the same as for the RT model. Plugging these expressions in eq. (3.15) we get an equation for $S(r)$ whose solution is

[TABLE]

where $\operatorname{Li_{\scriptscriptstyle 2}}(z)$ denotes the polylogarithm of order 2. For $|z|<1$ (i.e. $r>r_{S}/2$ since, in our case, $z=1-\frac{r_{S}}{r}$) it is given by

[TABLE]

Expanding for $r\gg r_{S}$, eq. (A.4) reduces to

[TABLE]

We can now plug these expressions for $U$ and $S$ into eqs. (3.12) and (3.13) to get the corrections ${\cal O}(m^{2}r^{2})$ to the metric. Again, we limit ourselves to $r\gg r_{S}$. The result is

[TABLE]

and

[TABLE]

where

[TABLE]

In terms of $A(r)=e^{2\alpha(r)}$ and $B(r)=e^{2\beta(r)}$, again to leading order in the limit $r\gg r_{S}$, we get

[TABLE]

A.2 Solution in the Newtonian limit

As in [68], we next compute the static solution in a Newtonian expansion over flat space. The resulting solution will be valid for $r\gg r_{S}$, but with no restriction on $mr$. For scalar perturbations, in the Newtonian gauge, the perturbed metric can be written as

[TABLE]

where we denote by $r_{N}$ the radial coordinate in the Newtonian gauge, while we reserve the notation $r$ for the radial coordinate in Schwarzschild coordinates used in the previous subsection. Comparing with the linearization of eq. (3.11), which is

[TABLE]

one obtains  [68] $r=(1+\Phi)r_{N}$ and

[TABLE]

We also expand the auxiliary fields as $U=\bar{U}+\delta U(r)$, $S=\bar{S}+\delta S(r)$. However, since, when expanding over Minkowski space, the background values $\bar{U}=\bar{S}$ are equal to zero, we simply write the perturbations as $U(r)$ and $S(r)$, keeping in mind that they are first-order quantities, just as $\Phi(r)$ and $\Psi(r)$. The linearization of eq. (2.17) gives

[TABLE]

while the linearization of $\Box S=-U$ gives

[TABLE]

Note that, in these equations, the Laplacian is with respect to the Newtonian coordinate $r_{N}$. However, since all quantities on the left and right-hand sides are already of first-order in the perturbations, we can equivalently use $r$. The linearization of the $(00)$ component of eq. (2.16) gives

[TABLE]

that, combined with eq. (A.16), becomes

[TABLE]

Finally, from the linearization of the $(ij)$ component of eq. (2.16) we get, again upon use of eq. (A.16),

[TABLE]

In particular, upon contraction with $\delta_{ij}$, we get

[TABLE]

We now observe that eqs. (A.15), (A.16), (A.18) and (A.20) are identical to the equations found in [68] in the context of the analysis for the RT model, with the field $S$, that in the RR model is defined by $S=-\Box^{-1}U$, playing the role of the field that was called $S$ in the context of the RT model, and defined by the fact that the spatial component $S_{i}$ of the auxiliary field $S_{\mu}$ can be decomposed as $S_{i}=S_{i}^{\rm T}+\partial_{i}S$, where $\partial_{i}S_{i}^{\rm T}=0$. This is not surprising; indeed, it was already observed in [39], using directly the nonlocal formulation, that the RR and RT model are identical when linearized over flat space. Here we have found the same result using the formulation in terms of auxiliary fields.

We can therefore now read the solution of these equations from the corresponding results in [68]. Note that one might be tempted to rewrite eq. (A.20) as $\Phi+\Psi-(m^{2}/3)S=0$, or eq. (A.15) as $U-2\Psi-4\Phi=0$. However, these equations have been derived in the limit $r\gg r_{S}$, where the linearization over flat space is appropriate, and are not valid for generic $r$. If a function $f$ satisfies {\mbox{\boldmath\nabla}}^{2}f=0 over all of space (and we further impose the boundary condition that $f$ vanishes at infinity), then $f=0$. However, from the fact that a function $f$ satisfies {\mbox{\boldmath\nabla}}^{2}f=0 at large $r$ we cannot conclude that $f$ itself is zero, even if we impose that it vanishes at infinity: we still have a freedom in the solution due to the fact that any function that, at large $r$, approaches the form $f(r)\simeq c_{1}r_{S}/r$ satisfies {\mbox{\boldmath\nabla}}^{2}f=0 at large $r$. The correct procedure is to fix constants such as $c_{1}$ by performing the matching with the solution found with the expansion in powers of $mr$. We quickly review below the correct procedure for finding the solution, referring the reader to [68] for more details.

First of all, to derive the solution for the auxiliary fields it is convenient to combine eqs. (A.16), (A.18) and (A.20) into an equation involving $U$ only,

[TABLE]

Note that this equation implies that, for a localized source, as in the Schwarzschild case that we are studying here, $U$ must go to zero at large distances. From this we can already anticipate that, when we perform the matching with the solution (A.3) in the overlap region $r_{S}\ll r\ll m^{-1}$, we will get $u_{0}=0$. For $\rho({\bf x})=M\delta^{(3)}({\bf x})$ (or, more generally, at distances $r$ much larger than the size or the localized source) the most general solution of eq. (A.21) is

[TABLE]

with $\beta$ an arbitrary constant. Plugging this into eqs. (A.15) and (A.17) we get the solution for $\Phi$ and $\Psi$ (again, valid only in the linearization region $r\gg r_{S}$),

[TABLE]

where $c_{\Phi}$ and $c_{\Psi}$ are constants, for the moment arbitrary, that reflect the freedom of adding to the solution at large distances a term $\propto 1/r$, as discussed above. Note that we still require that $\Phi$ and $\Psi$ go to zero at infinity, as appropriate for the metric perturbations generated by a localized source. Using eq. (A.14) we can read the corresponding expression for $\alpha(r),\beta(r)$, and therefore for $A(r)\simeq 1+2\alpha(r)$ and $B(r)\simeq 1+2\beta(r)$. The result is

[TABLE]

Thus, the solution for the metric found with a Newtonian expansion around flat space, which is valid for $r\gg r_{S}$ but with no restriction on $mr$, has three undetermined parameters $\beta,c_{\Psi},c_{\Phi}$. In contrast, we have seen that the solution found with an expansion in powers of $mr$, again assuming $r\gg r_{S}$, given by eqs. (A.7) and (A.8) [or, equivalently, by eq. (A.11) for $B(r)$ and the corresponding expression for $A(r)$] has two free parameters $u_{0},u_{1}$. Quite nicely, all these parameters can be fixed by comparing the two solutions in the overlap region $r_{S}\ll r\ll m^{-1}$, where they are both valid. Indeed, taking $mr\ll 1$ in eqs. (A.25) and (A.26) we get

[TABLE]

Comparing eqs. (A.10) and (A.11) to eqs. (A.27) and (A.28) we get $c_{\Phi}=2/3$, $c_{\Psi}=-4/3$, $\beta=0$, $u_{0}=0$ and $u_{1}=1$. The fact that $u_{0}$ is forced by the matching to vanish has important consequences. Indeed, if the coefficient of the term $m^{2}r^{2}(r/r_{S})$ in eqs. (A.10) and (A.11) would have been of order one, this term would have become of order one for r\,\raise 1.72218pt\hbox{>}\kern-6.99997pt\lower 2.6694pt\hbox{\sim}\,r_{c}=(r_{S}m^{-2})^{1/3}. This scale is much smaller that $m^{-1}$ since, from the comparison with cosmological data, the mass scale $m$ of the nonlocal model is of order $H_{0}$ and therefore $r_{S}\ll m^{-1}$. Thus, the solution given by eqs. (A.10) and (A.11) would have entered a ‘strong coupling’ regime, where the expansion in powers of $mr$ breaks down, already at a critical value of $r_{c}$ such that $mr_{c}\ll 1$. However, since $u_{0}$ is fixed by the matching procedure to the value $u_{0}=0$, the coefficient of this term vanishes, both in $A(r)$ and in $B(r)$, and the solution (A.10), (A.11) is valid until $mr$ approaches a value of order one.

For $A(r)$ and $B(r)$ we therefore end up with the solution already announced in eqs. (3.16) and (3.17) while, for $\Phi$ and $\Psi$, eqs. (A.23) and (A.24) can be rewritten as

[TABLE]

while for $U(r)$ we get

[TABLE]

We can now compute the corresponding solution for $S(r)$ (that was not explicitly studied in [68, 39], but that we need here because we are eventually interested in matching the static solution for $S(r)$ to its cosmological solution, for obtaining the effective Newton constant). Using eq. (A.31), eq. (A.16) becomes

[TABLE]

which, as usual, is valid only for $r\gg r_{S}$. This can be rewritten as

[TABLE]

whose solution is

[TABLE]

with $s_{0},s_{1}$ dimensionless integration constants. In the limit $mr\ll 1$ this becomes

[TABLE]

From eq. (A.6), with $u_{0}$ and $u_{1}$ now fixed to $u_{0}=0$ and $u_{1}=1$, we know however that the asymptotic behavior of $S(r)$ when $mr\ll 1$ must be $S(r)\simeq-(1/2)rr_{S}$. Comparison with eq. (A.35) then shows that the term $s_{0}+(s_{1}+1)r_{S}/r$, which is a zero-th order term with respect to the small parameter $mr$, is actually absent and this fixes $s_{0}=0$ and $s_{1}=-1$. Thus, in the end

[TABLE]

This result can be checked observing that, through eq. (2.16), the spatial dependence of $S(r)$ induces a spatial dependence of Newton’s constant, because of the term $m^{2}SG_{\mu\nu}$ in $m^{2}K_{\mu\nu}$. However, the resulting modification of Newton’s constant must be contained in the motion of a test particle in the Schwarzschild metric (3.11), (3.16), (3.17). The latter shows that any effect vanishes in the limit $mr\rightarrow 0$. This means that, in eq. (A.34), we must have $s_{0}=0$ and $s_{1}=-1$ so that, in the limit $mr\rightarrow 0$, $m^{2}S(r)$ vanishes.

Observe that we have found that the solution $S$ of {\mbox{\boldmath\nabla}}^{2}S=-U is everywhere negative, including at $mr\ll 1$, even if the solution for $U$, at $mr\ll 1$, is positive. This is not in contradiction with the analysis performed in section 4, where we found that, with $R$ positive, the solution for $U$ of {\mbox{\boldmath\nabla}}^{2}U=-R is positive, a result that invalidated the screening mechanism that had been proposed for the DW model. The point, as discussed in footnote 16 on page 16, is that for $R$ the leading contribution is given by the source density $\rho$, which has compact support. The solution, at leading order, can then be found with the Green’s function (4.3). In contrast, in {\mbox{\boldmath\nabla}}^{2}S=-U the source term $U$ does not have compact support and the solution can no longer be obtained with a convolution with this Green’s function. Rather, we have obtained it from the direct integration of eq. (A.32) that, thanks to the $1/r$ behavior of the source term, takes the simple form (A.33).

Appendix B Details of the computation in the McVittie metric

In this appendix we give the details of the derivation of eq. (3.37). We start from eq. (3.30) and we use the conditions (3.28) and (3.29). Then, in eq. (3.30), to lowest order all the occurrences of $\Phi$ can be dropped. Indeed, in $(1-2\Phi-H^{2}r^{2})$, the term $2\Phi$ is always subleading with respect to one, both in the region $r\ll r_{c}$ where the McVittie metric is dominated by the central mass, and in the region $r\gg r_{c}$ where the metric approaches FRW. This is different to what happens to the $H^{2}r^{2}$ term, that is totally negligible with respect to one in the mass-dominated region $r\ll r_{c}$, but is not parametrically small with respect to one at the horizon scale. Similarly, we can neglect the $\Phi$ and $r\partial_{r}\Phi$ terms in the coefficients of $\partial_{t}^{2}U$ and $\partial_{t}U$ in eq. (3.30) [or $\partial_{t}^{2}S$ and $\partial_{t}S$ in eq. (3.31)]. A similar analysis can be done for the dimensionless term $[2r\partial_{r}\Phi+r^{2}(\dot{H}+2H^{2})]$ in eqs. (3.30) and (3.31): even in the regime $r_{S}\ll r\ll r_{c}$, which is dominated by the central mass, the term $2r\partial_{r}\Phi$ can be neglected since it is much smaller than one and therefore, as of $(1/r)\partial_{r}U$ in eq. (3.30), is negligible compared to the term $(1-2\Phi-H^{2}r^{2})(2/r)\partial_{r}U\simeq(2/r)\partial_{r}U$ [and similarly for the term $(1/r)\partial_{r}S$ in eq. (3.31)].

Of course, this does not mean that we have lost all the dependence on the central mass $M$. The mass enters through the source term $R$, e.g. through a term proportional to $\rho({\bf x})=M\delta^{(3)}({\bf x})$ if, for the purpose of studying the solution at $r\gg r_{S}$, we approximate the mass as point-like. The situation is precisely the same as in the Newtonian approximation discussed in App. A.2 where we found that, to leading order in $r_{S}/r$, $U$ satisfies eq. (A.21), ({\mbox{\boldmath\nabla}}^{2}+m^{2})U=-8\pi G\rho, where on the left-hand side there is no occurrence of $M$ nor of the metric perturbation $\Phi$, but the mass enters through $\rho$ on the right-hand side. Dropping all occurrences of $\Phi$, eq. (3.30) becomes

[TABLE]

In the RR model the Ricci scalar $R$ can be determined by taking the trace of eq. (2.16), which gives (3.32), where the trace of the energy-momentum tensor $T$ is given by eq. (3.33).

We now wish to test whether (at least with these approximations) eq. (B.1) is solved by the ansatz

[TABLE]

where $U_{\rm cosmo}(t)$ is the cosmological solution discussed in section 3.1 and $U_{\rm static}(r)$ is the static Schwarzschild solution of section 3.2 (a completely analogous analysis can be performed for the equation obtained from $\Box S=-U$).

The terms that would forbid a separation of variables of this sort are those in which there is a dependence both on $r$ and on $t$. Several of these terms have already disappeared when we dropped $\Phi(r)$. Among the remaining terms, consider for instance the term

[TABLE]

that comes out when we insert the ansatz (B.2) in the first term in eq. (B.1). A priori it depends both on $t$ and $r$. However, in the regime $r_{S}\ll r\ll m^{-1}$ we have seen that $U_{\rm static}\simeq r_{S}/r$, and then the right-hand side of eq. (B.3) vanishes. More generally, the right-hand side of eq. (B.3) will be suppressed by a factor $r_{S}/r$. For instance, using the solution (3.19) valid in the actual metric for the RR model, we get

[TABLE]

so indeed in the cosmological regime this term is suppressed by a factor $r_{S}/r$. At cosmological distances this suppression factor is of the order of the ratio of the Schwarzschild radius of the central mass (about 1 cm for the Earth) to the size of the horizon of the Universe, so it is totally negligible. Therefore the term (B.3) is only relevant in the mass-dominated region $r\ll r_{c}$. However, in this regime the term $H^{2}(t)r^{2}$ is completely negligible compared to one, so we can replace everywhere the term (B.3) by

[TABLE]

that now depends only on $r$. The same analysis holds for the term $r^{2}(\dot{H}+2H^{2})(1/r)\partial_{r}U_{\rm static}$ in eq. (B.1), that is negligible at cosmological scales because of $(1/r)\partial_{r}U_{\rm static}$, and also in the mass-dominated region because there it is a correction ${\cal O}(H^{2}r^{2})$ to the term $(2/r)\partial_{r}U_{\rm static}$ in eq. (B.5). Thus, the term $r(\dot{H}+2H^{2})\partial_{r}U_{\rm static}$ can be dropped everywhere. Finally, the term $\partial_{t}\partial_{r}U$ vanishes on the ansatz (3.34).

Thus, on the ansatz (B.2), the left-hand side of eq. (B.1) can be replaced by

[TABLE]

which leads to the left-hand side of eq. (3.37). We see that the left-hand side of the equation nicely separates into the $r$-dependent part of $\Box U$ that determines the static solution and the $t$-dependent part relevant for the cosmological solution.

We next study the source term $-R$ that, according to eqs. (3.32) and (3.33), on the ansatz (B.2), becomes

[TABLE]

where, in the last line, in $3U+U^{2}-\partial_{\mu}S\partial^{\mu}U$, of course $U=U_{\rm cosmo}(t)+U_{\rm static}(r)$ and $S=S_{\rm cosmo}(t)+S_{\rm static}(r)$.

We now observe, from eq. (3.20), that $m^{2}S_{\rm static}(r)$ is always very small compared to one, as long as $r\gg r_{S}$. Indeed, at cosmological scales, where $mr$ is not small, it is anyhow suppressed by a factor $r_{S}/r$ while, if $r_{S}\ll r\ll m^{-1}$, it is suppressed both by a factor $r_{S}/r$ and by a factor $(mr)^{2}$. Thus, we always have $|m^{2}S_{\rm static}(r)|\ll 1$ (as usual, as long as $r\gg r_{S}$), and we can neglect it in eq. (B.7). Thus, in the denominator of the first and second term in eq. (B.7), we can replace $\{1-(m^{2}/3)[S_{\rm cosmo}(t)+S_{\rm static}(r)]\}$ by $[1-(m^{2}/3)S_{\rm cosmo}(t)]$.

This means the the first term in eq. (B.7) is exactly the same that one get in the cosmological solution, when one sets $U=U_{\rm cosmo}(t)$, while the second term is the one would get when one studies the static solution, setting $U=U_{\rm static}(r)$, except that now $G$ is replaced by $G/[1-(1/3)m^{2}S_{\rm cosmo}(t)]$, something which is missed when considering only the static solution. The combination $G/[1-(1/3)m^{2}S_{\rm cosmo}(t)]$ is nothing but the effective Newton constant found in eq. (2.21), so the second term in eq. (B.7) can be rewritten as $-8\pi G_{\rm eff}(t)M\delta^{(3)}({\bf r})$.

Finally, consider the last term in eq. (B.7). When evaluated on $U_{\rm cosmo}(t)$, $S_{\rm cosmo}(t)$ it gives the term that, in the cosmological solution, combines with the energy-momentum tensor of the cosmological fluid, and provides a dynamical dark energy, as studied in section 3.1. When evaluated on the purely static solution it rather gives the non-linear terms that contribute to the static Schwarzschild solution studied in section 3.2. However, when we evaluate it on the ansatz $U=U_{\rm cosmo}(t)+U_{\rm static}(r)$, $S=S_{\rm cosmo}(t)+S_{\rm static}(r)$, $G$ is replaced by $G_{\rm eff}(t)$. In principle, there is also a mixed term coming from $U^{2}-\partial_{\mu}S\partial^{\mu}U$, given by

[TABLE]

This term, however, can be neglected. In fact, at $mr\ll 1$, it is suppressed because of the $m^{2}$ factor in front of it, while at cosmological scales it is suppressed by the factors $r_{S}/r$ coming from the static solutions. As a result, $-R$ reduces to the expression given in the right-hand of eq. (3.37).

Appendix C Details of the computation in the perturbed FRW metric

In this appendix we present the details of the computation of section 3.4. We begin by discussing the equation $\square U=-R$. Using the ansatz (3.48), (3.49) and making use of the condition $|\Phi|\ll 1$, from eq. (3.46) we get

[TABLE]

As in section 3.3, to get the Ricci scalar we take the trace of eq. (2.16) in the perturbed FLRW metric, which gives191919We have neglected the pressure term and we have written $T=-\rho$, as appropriate for our problem where the sources are non-relativistic structures. The pressure term could be easily reinstated, but plays no role in the following.

[TABLE]

Dropping again the ${\cal O}(\Phi)$ corrections and using the perturbative ansatz (3.48), (3.49), we get

[TABLE]

where $\delta\dot{U}\equiv\partial_{t}(\delta U)$, $\delta\dot{S}\equiv\partial_{t}(\delta S)$. Note that the term proportional to $\vec{\nabla}\delta U\cdot\vec{\nabla}\delta S$ is formally of second order in the perturbations, but for the moment we keep it because it could be enhanced by the spatial derivatives (even if we will show below that it is actually negligible). Now we further assume a quasi-static approximation

[TABLE]

and we also assume

[TABLE]

whose validity will be checked a posteriori. Thus, the equation $\Box U=-R$ becomes

[TABLE]

We now observe that, using eq. (3.52),

[TABLE]

On the other hand $|\vec{\nabla}\delta U|^{2}/(3a^{2})$ is much smaller than the term $a^{-2}\nabla^{2}\delta U$ on the left-hand side of eq. (C.6), since it contains one more power of the small perturbation $\delta U$ and has the same number of spatial derivatives. So, the term $m^{2}|\vec{\nabla}\delta U\cdot\vec{\nabla}\delta S|/(3a^{2})$ on the right-hand side of eq. (C.6) can be dropped. Similarly, using first the quasi-static approximation (C.4) and then eq. (C.5), we can show that, for modes inside the horizon also the terms $\dot{S}_{c}\delta\dot{U}$ and $\dot{U}_{c}\delta\dot{S}$ can be dropped. As for the term $-\delta U-\frac{2}{3}U_{c}\delta U$ on the right-hand side of eq. (C.6), for modes well inside the horizon today, this is negligible compared to the term $a^{-2}\nabla^{2}\delta U$ on the left-hand side, since $m$ is of order $H_{0}$.

Thus, we finally get

[TABLE]

For vanishing $\delta U$ and $\Phi$ this reduces to the background equation for $U_{c}$ sourced by the background density $\rho_{b}$, which we subtract to get

[TABLE]

as claimed in eq. (3.53).

We finally need to check the consistency of the assumption (C.5). To this purpose we consider the equation $\Box S=-U$, that, in the same approximations become

[TABLE]

From this we subtract the background equation $-\ddot{S}_{c}-3H\dot{S}_{c}=-U_{c}$ to get

[TABLE]

The term $\delta U$ here must be dropped, since we found that $\delta U\simeq 2\Phi$ and we have already dropped the terms ${\cal O}(\Phi)$ from the wave operator.202020In any case, this term would source a value of $\delta S$ such that $(k/a)^{2}\delta S_{k}\sim\delta U_{k}$, i.e $m^{2}\delta S_{k}\sim(m^{2}/k_{\rm phys})^{2}\delta U_{k}$, where $k_{\rm phys}=k/a$. Since $m\sim H_{0}$, for modes well inside the horizon this would give $m^{2}\delta S_{k}\sim(m^{2}/k_{\rm phys})^{2}\delta U_{k}\ll\delta U_{k}$. As for the term $4\dot{\Phi}\dot{S}_{c}$, we can estimate its typical value as follows. The time evolution of $\Phi$ is determined by the typical time-scale $T$ of the evolving matter structures, $\dot{\Phi}\sim\Phi/T$, while the cosmological background solution evolves with a time-scale given by $H$, $\dot{S}_{c}\sim HS_{c}$. In the recent cosmological epoch, when $H\sim H_{0}$, we have

[TABLE]

where, as usual, $V(t)=H_{0}^{2}S(t)$, and today $V$ is of order one. The contribution of structures that evolve on a time-scale $T$ shorter than the Hubble time is therefore such that the term $4\dot{\Phi}\dot{S}_{c}$ in eq. (C.11) dominates over $\delta U$. Nevertheless, even this source term induces a value of $\delta S$ that satisfies eq. (3.52). Indeed, eq. (C.11) gives

[TABLE]

Writing $T\sim\lambda/v=1/(k_{\rm phys}v)$, where $\lambda=1/k_{\rm phys}$ is the typical size of a structure and $v$ its velocity, and recalling that $\Phi\sim\delta U$ and $m\sim H_{0}$, we get

[TABLE]

since both $H_{0}/k_{\rm phys}\ll 1$ and $v\ll 1$. At early times (approaching matter-radiation equality) $\delta S$ is even more suppressed because its sources, which are perturbations around the background, were even smaller.