Spherically symmetric static solutions in a non-local infrared modification of General Relativity

We discuss static spherically symmetric solutions in a recently proposed non-local infrared modification of Einstein equations induced by a term $m^2g_{\mu\nu}\Box^{-1} R$, where $m$ is a mass scale. We find that, contrary to what happens in usual theories of massive gravity, in this non-local theory there is no vDVZ discontinuity and classical non-linearities do not become large below a Vainshtein radius parametrically larger than the Schwarzschild radius $r_S$. Rather on the contrary, in the regime $r\ll m^{-1}$ the corrections to the metric generated by a static body in GR are of the form $1+{\cal O}(m^2r^2)$ and become smaller and smaller toward smaller values of $r$. The modification to the GR solutions only show up at $r>m^{-1}$. For $m={\cal O}(H_0)$, as required for having interesting cosmological consequences, the non-local theory therefore recovers all successes of GR at the solar system and lab scales.


Introduction
Recent years have witnessed very intense activity on the study of infrared modifications of General Relativity (GR). This is motivated by the aim of explaining the present phase of accelerated expansion of the Universe through a modification of gravity at distances r ∼ H −1 0 (where H 0 is the present value of the Hubble parameter), and also turns out to be a rich and challenging theoretical subject. To modify GR in the far infrared it is natural to introduce a mass parameter m ∼ H 0 . At first sight, this could be achieved by giving a mass to the graviton. However, constructing a consistent theory of massive gravity turns out to be remarkably difficult, and the problem has a long history that goes back to classic papers by Fierz and Pauli in 1939 [1] and Boulware and Deser in 1973 [2]. In recent years there has been significant progress in this direction, in particular with the construction of the ghost-free dRGT theory [3,4] (see also [5][6][7][8], and [9,10] for reviews), although a number of difficulties and open problems persist; the dRGT theory (as well as galileon theories) very likely admits superluminal excitations over some backgrounds [11][12][13][14][15][16][17][18][19]. Furthermore, even if the sixth ghost-like degree of freedom is absent in any background, the fluctuations of the remaining five degrees of freedom can become ghost-like over non-trivial backgrounds [20][21][22]. Another open problem of dRGT theory is that it is not clear whether a satisfying cosmology emerges. Homogeneous and isotropic spatially flat Friedman-Robertson-Walker (FRW) solutions do not exist, and are in fact forbidden by the same constraint that removes the ghost [23]. There are open isotropic FRW solutions, which however suffer of strong coupling and ghost-like instabilities [24]. It is presently unclear whether there are stable and observationally viable inhomogeneous solutions in the full non-linear theory away from the decoupling limit (see the discussion in [10]).

JHEP08(2014)029
A peculiar aspect of massive gravity theories is that they require the introduction of an external reference metric. In the recent papers [25,26] it has been proposed a different approach, in which gravity is deformed by the introduction of a mass parameter m in such a way that no external reference metric is introduced, and general covariance is preserved. This can be achieved by adding nonlocal terms to the Einstein equations. The introduction of nonlocal terms for producing IR modifications of gravity has been suggested by various authors, following different lines of reasoning. In particular, nonlocal operators that modify GR in the IR appear in the degravitation proposal [27,28] (see also [29,30]). Non-local covariantizations of the Fierz-Pauli theory were discussed in [31]. A different nonlocal cosmological model has been proposed in [32], and has been further studied in a number of recent papers [33][34][35][36][37][38][39][40][41][42][43] (see [44] for a recent review). Another interesting nonlocal model has been studied in [45][46][47]. The model that has been proposed by one of us in [26] is defined by the classical equation of motion The inverse of the d'Alembertian is defined with the retarded Green's function, which ensures causality. The superscript T denotes the extraction of the transverse part of the tensor and ensures that the left-hand side of eq. (1.1) has zero divergence, and therefore T µν is automatically conserved. The extraction of the transverse part exploits the fact that, in a generic curved space-time, any symmetric tensor S µν can be decomposed as where ∇ µ S T µν = 0 [48,49]. The factor 1/3 in eq. (1.1) is a convenient normalization of the parameter m 2 in d = 3 spatial dimensions (and becomes (d − 1)/(2d) for generic d). Some conceptual aspects of this model have been discussed in [26,50]. Its cosmological consequences, at the level of background evolution, have been studied in [26,51], while a study of its cosmological perturbations will be presented in [52].
At the conceptual level, it is important to stress that the −1 operator in eq. (1.1) is defined with the retarded Green's function. This ensures causality, and also has the important consequence that eq. (1.1) cannot be the equation of motion of a fundamental nonlocal QFT. Indeed, the variation of an action involving −1 always gives rise to an equation of motion involving a symmetrized Green's function, rather than a retarded one [25,32,47]. Equation (1.1) should rather be understood as a classical effective equation of motion. Non-local effective equations involving a retarded Green's function govern for instance the dynamics of the in-in matrix elements of quantum fields, such as 0 in |φ|0 in or 0 in |ĝ µν |0 in , and encode quantum corrections to the classical dynamics [53,54]. Thus, issues of quantum vacuum decay induced by ghost instabilities, or non-linearities induced by quantum corrections, cannot be addressed directly from a study of eq. (1.1), but should rather be addressed in the fundamental underlying (local) QFT. This, however, is not conceptually very different from what happens in dRGT, where the UV completion is needed to address the causality issue.

JHEP08(2014)029
Observe that one can still in principle derive non-local causal equations from an action, using the formal trick of replacing by hand −1 → −1 ret after performing the variation. This is indeed the procedure used in [32,55], in the context of non-local gravity theories with a Lagrangian of the form Rf ( −1 R), see also the recent discussion in [56]. However, any direct connection to a fundamental quantum field theory is then lost. One can ask what is the classical non-local action that, in the above sense, reproduce eq. (1.1). If one linearizes eq. (1.1) over flat Minkowski space, writing g µν = η µν + h µν , one gets [26] where E µν,ρσ is the Lichnerowicz operator, and is now the flat-space d'Alembertian. This linearized equation can of course be derived from the action where S M is the matter action. This can be rewritten in a covariant form, as Thus, the action corresponding to eq. (1.1) contains a full series of cubic and higher-order terms in the curvature, and we do not have a compact closed-form expression. The model obtained truncating eq. (1.6) to order R 2 is interesting in its own right, and has been studied in [57]. Observe that the cubic and higher order terms are suppressed by powers of −1 R = O(h) (and not of R/M 2 Pl , where M Pl is the Planck mass). They are therefore on the same footing as the usual non-linearities of GR and cannot be neglected on non-trivial backgrounds, e.g. in cosmological backgrounds. So, the model defined by with no cubic and higher-order terms, and the model defined by eq. (1.1) are different. In a sense, eq. (1.1) provides the simplest non-local equation of motion in this class of model involving −1 R and a mass scale m, while eq. (1.7) provides the simplest action. At the phenomenological level, eq. (1.1) turns out to have rather interesting consequences. In particular, at the level of background evolution it admits flat FRW solutions, in which furthermore a dynamical dark energy emerges automatically. By fixing the free parameter m to a value m 0.67H 0 the model reproduces the observed value Ω DE 0.68. This leaves us with no free parameter and we then get a pure prediction for the EOS parameter of dark energy. Using the standard fit of the form w DE (a) = w 0 +(1−a)w a [58,59], the model predicts w 0 −1.04 and w a −0.02 [26], consistent with the Planck data [60], and on the phantom side. This should be compared with models such as that of ref. [32], JHEP08(2014)029 which involves an arbitrary function f ( −1 R), which can be chosen so to reproduce any expansion history. The model (1.1) is therefore very predictive, and passes remarkably the non-trivial test of giving an equation of state consistent with the existing limits. 1 The purpose of the present paper is to continue the investigation of the model defined by eq. (1.1), addressing in particular the issue of its classical non-linearities by studying the static spherically symmetric solutions. A typical issue of massive gravity theories is that they become non-linear when r is smaller than the Vainshtein radius, which is parametrically larger than the Schwarzschild radius r S of the source. For instance, in the theory defined by adding a Fierz-Pauli mass term to the Einstein-Hilbert action, one finds that the classical non-linearities become large below the Vainshtein radius r V = (GM/m 4 ) 1/5 [61,62]. In the dRGT theory [3,4] the strong-coupling energy scale is raised and the corresponding critical distance is lowered, to r V = (GM/m 2 ) 1/3 [63], which is still of order of 100 pc for m = O(H 0 ) and M = M . In order to recover the successes of GR at solar system and shorter scales, one must then show that a Vainshtein mechanism is at work, i.e. that the inclusion of classical non-linearities restore continuity with GR at r r V . Explicit examples of this type have indeed been found for the dRGT theory [64,65].
In the nonlocal model (1.1) the situation seems however different. Indeed, linearizing the theory over flat space, one finds that the matter-matter gravitational interaction mediated by this theory is given by [26] The first term is the usual GR result due to the exchange of a massless graviton, while the extra term vanishes for m → 0. Therefore this theory has no vDVZ discontinuity, and no Vainshtein mechanism is needed to restore continuity with GR. Of course, the fact that we do not need a Vainshtein mechanism does not necessarily mean that non-linearities will remain small down to the Schwarzschild radius r S , where also the classical non-linearities of GR become large. So, the purpose of this paper is to study static spherically symmetric solutions in this theory, and compare with the corresponding solutions of GR. In our problem we have two independent length-scales, the Schwarzschild radius r S of the source, and the length-scale m −1 . To have interesting and viable cosmological applications we must have m = O(H 0 ) (indeed, the analysis of [26] shows that the model generates a dynamical dark energy with the observed value of Ω DE if we choose m 0.67H 0 ). Thus, between these two scales there is a huge separation, r S m −1 . At scales r ∼ m −1 we expect that the nonlocal theory (1.1) will differ from GR. Indeed, the motivation for such 1 Similar interesting cosmological consequences follow from the model (1.7), and will be explored in [26,52]. Of course, as always in model building, there is always a freedom in the choice of the model itself, and in this sense the study of any single model in this class should be considered as an example of the typical consequences of nonlocal terms that can be associated to a mass parameter m.

JHEP08(2014)029
a model is just to produce a modification of GR in the far infrared, that could account for the observed acceleration of the Universe. The main motivation of this paper is to see if, in the region r S r m −1 , the theory remains linear and close to GR (while, of course, as r → r S , even GR becomes non-linear), and to compute explicitly the deviations from GR in the region r ∼ m −1 , where they could become relevant for comparison with structure formation.
The paper is organised as follows. In section 2 we write down the equations of motion in spherical symmetry. In section 3 we solve these equations analytically in the regime r m −1 using a low-m expansion. In section 4 we solve them in the region r r S , using the linearization over the Minkowski background, and we show that in the region r S r m −1 the solution overlaps with that found in section 3. The results are confirmed through a numerical analysis in section 5. Finally, in section 6 we give a discussion of the radiative and non-radiative degrees of freedom of the nonlocal theory, which is useful for a physical understanding of the results obtained. Section 7 contains our conclusions.

Basic equations
We look for static spherically symmetric solutions of eq. (1.1). As in [26], we define Then the original non-local equation (1.1) can be formally rewritten as a system of local equations for the variables g µν , S µ and U , and where the latter equation is obtained taking the divergence of eq. (1.2). Some subtleties involved in this localization procedure will be discussed below. We write the most general static spherically symmetric metric in the form ds 2 = −e 2α(r) dt 2 + e 2β(r) dr 2 + r 2 (dθ 2 + sin 2 θ dφ 2 ) . (2.7) Observe that the nonlocal equation (1.1) is generally covariant. Therefore, just as in GR, we can use the invariance under diffeomorphisms to set to one a function e 2µ(r) that otherwise, in the most general spherically symmetric solution, would multiply the term JHEP08(2014)029 r 2 (dθ 2 + sin 2 θ dφ 2 ), and that indeed must be kept in local massive gravity theories. 2 We use the labels (0, 1, 2, 3) for the indices (t, r, θ, φ) and we denote df /dr by f . The nonvanishing Christoffel symbols in the metric (2.7) are Γ 0 01 = α , Γ 1 00 = e 2(α−β) α , plus those related by the symmetry Γ ρ µν = Γ ρ νµ . Using these expressions we can compute B µν . In a spherically symmetric spacetime, for symmetry reasons it is clear that S 2 = S 3 = 0, and S 0 = S 0 (r), S 1 = S 1 (r). Then, the non-vanishing components of B µν are (2.9) To compute S 0 and S 1 we use eq. (2.6). The equations with ν = 2, 3 are automatically satisfied by S 2 = S 3 = 0. Setting ν = 0 gives an equation that only involves S 0 , which has the solution S 0 = 0. In principle it also admits other solutions. For instance, the non-vanishing solution of (∂ r − 2α )S 0 = 0 is S 0 = c 0 e 2α(r) , with c 0 a constant. However, the correct solution is uniquely specified by the condition that S 0 must be equal to zero when U = 0, since for U = 0 we have S µν = 0, which is already trivially transverse, and a non-vanishing S µ is in this case a spurious solution. Such spurious solutions typically arise when writing the original nonlocal equation as a system of local differential equations, introducing auxiliary fields such as, in our case, U and S µ [41,47,50,51,66,67] (see also the discussion in section 6). Thus, in our case we only retain the solution S 0 = 0, and the only non-vanishing auxiliary fields are U (r) and S 1 (r). Taking the ν = 1 component of eq. (2.6) we get The equation for U is obtained from eq. (2.1), written in the form U = −R. In general, on a scalar function U , However, on a function U (r) in a diagonal metric such as (2.7), U = (−g) −1/2 ∂ r [ √ −gg rr ∂ r U (r)] is just a covariant Laplacian. In the metric (2.7) this gives 2 In massive gravity models, where there is both a dynamical metric and a reference metric, the assumptions of staticity and of spherical symmetry are not sufficient to put both of them in this form, and in one of them remains a term 2D(r)dtdr [9]. As discussed in [64,65] in dRGT there are two possible branches of solutions: a branch with D(r) = 0, which is asymptotically flat, exhibits a vDVZ discontinuity at r rV , and recovers GR at r rV , thereby giving an explicit example of the Vainshtein mechanism; and a branch where D(r) is a non-vanishing function, which corresponds to a Schwarzschild-de Sitter solution. In our case, however, there is no reference metric, and for a static and spherically symmetric source we can set D(r) = 0.
Finally, the system of equations for the four functions {α, β, U, S 1 } is completed by taking any two independent components of the nonlocally modified Einstein equations (1.1).
We now study these equation in the region outside the source, where T µν = 0. Let us recall that in GR one typically takes the combinations e 2(β−α) R 00 + R 11 and R 22 (see e.g. [68]). In vacuum this gives e 2(β−α) R 00 + R 11 = 0 and R 22 = 0. The former equation gives α = −β and then the latter gives a differential equations for α. In our nonlocal theory, using eq. (2.14) we get instead Observe, from eq. (2.16), that now α = −β, unless S 1 is constant. Equations (2.11), (2.12), (2.16) and (2.17) provides four differential equations for the four functions α(r), β(r), U (r) and S 1 (r). Finally, it is convenient to trade S 1 for a field V (r) defined by S 1 (r) = e β rV (r), and work with the four dimensionless functions α(r), β(r), U (r) and V (r). In terms of V the final form of our equations is Observe that replacing S 1 by V eliminates the term β from eq. (2.11), and therefore from the whole system of equations.
3 Solution for r m −1 We now study eqs. (2.18)-(2.21) in the region mr 1, performing a low-m expansion. We assume that in the limit mr → 0 the terms on the right-hand side of eqs. (2.20) and (2.21)

JHEP08(2014)029
are negligible, and we will then check a posteriori the self-consistency of the assumption. Consider first the equations in the external region, where T µν = 0. In this case, to lowest order eqs. (2.20) and (2.21) reduce to their standard GR form, whose solution is given by Plugging these expressions into eq. (2.19) we get with u 0 , u 1 some constants, that parametrize the solution of the associated homogeneous equation. Observe that the inhomogeneous solution vanishes since the right-hand side of eq. (2.19) is zero on the unperturbed Schwarzschild solution (as it is also obvious from the fact that it is just r 2 R, and the Ricci scalar R vanishes on the Schwarzschild solution). The choice of homogeneous solution is a delicate point that requires some discussion. As discussed in detail in [26,50,51], this issue is related to the fact that, in order to complete the definition of the nonlocal model (1.1), we must specify what we actually mean by −1 .
In general, an equation such as U = −R is solved by where U hom (x) is any solution of U hom = 0 and G(x; x ) is any a Green's function of the operator. To define our nonlocal integro-differential equation we must specify what definition of −1 we use, i.e. we must specify the Green's function and the solution of the homogeneous equation. In our static setting the definition of G(x; x ) is irrelevant, since if R(x ) and −g(x ) are independent of t , all possible definitions the above equation reduce to is the Green's function of the covariant Laplacian. Still, the definition of the −1 operator is completed only once we specify U hom (x). In general, this will be fixed by the boundary conditions of the specific problem that we consider, which in our static case will therefore fix u 0 and u 1 . Similarly, we must specify the homogeneous solution associated to the definition of S µ , eq. (2.6), which is equivalent to completely define the nonlocal operation of taking the transverse part. In other words, u 0 and u 1 (and the similar constants that characterize the homogeneous solution of S µ ) are not parameters that can be varied and that classify all possible classical solutions of a given theory. Rather, they are fixed once and for all by the definition of the original nonlocal theory and the boundary conditions of the problem at hand. In particular, a non-vanishing value of u 0 corresponds to introducing a cosmological constant term in the theory. Indeed, denote by U old and U new two different definitions of U related by U new = U old + u 0 . Then the nonlocal theory using the definition U new , and is therefore the same as the old theory, in which we add to the right-hand side a cosmological constant Λ = −(1/3)m 2 u 0 . In this paper we consider the nonlocal model defined by setting u 0 = 0. For such a model, our aim is to show that the static spherically symmetric solutions of the theory reduce to the Schwarzschild solution of GR as m → 0.
Of course, if we switch on Λ, we should rather show that they approach the corresponding Schwarzschild-dS or (depending on the sign of Λ) Schwarzschild-AdS solutions. The appropriate choice of u 1 is more subtle. We will for the moment keep it generic, and we will later see how it can be uniquely determined. We therefore write We now plug these expressions for α, β and U into eq. (2.18), and we get plus the solutions of the corresponding homogeneous equation, that we set to zero on the ground that, when U = 0, we must have V = 0, since in this case U g µν = 0 and there is no transverse part to extract. Observe that, for r r S , this expression reduces to V (r) −u 1 r S /(2r). We have therefore obtained the solution for α, β U and V to zero-th order in m. To get the first correction to α and β we plug these expression for U, V back into eqs. (2.20) and (2.21) and solve them. In principle this can be done for generic r, as long as r m −1 , but the resulting solutions, involving polylog functions, are quite long and not very illuminating, so we we write down the correction term only in the limit r S r. Since we are treating mr perturbatively, we are then actually studying the solution in the region r S < ∼ r m −1 . The region r S r m −1 is particularly interesting in view of the fact that, in typical massive gravity theories, the classical theory becomes non-linear below a Vainshtein radius r V parametrically larger than r S . Studying the solution in this region allows us to investigate whether the same phenomenon happens in the nonlocal theory. In this regime we can use the zero-th order expressions for α, β, U and V and expand them to leading order in r S /r. Plugging these expressions on the right-hand side of eq. (2.20), to leading order in r S /r we still find 3 α and, to first non-trivial order, eq. (2.21) becomes (3.10)

JHEP08(2014)029
Since we are computing the m 2 correction only in the limit r r S , in the term proportional to m 2 we can approximate e 2β 1, so To the order at which we are working, the solution can be written in the form Thus, in the region r S r m −1 , The above result shows that our perturbative procedure is self-consistent, since the corrections to linearized theory are indeed small, as long as mr 1. It is clear that the procedure can be iterated, obtaining a systematic expansion in the small parameter (mr) 2 . This should be contrasted with what happens in massive gravity, when one considers the Einstein-Hilbert action plus a Fierz-Pauli mass term. Then the analogous computation gives, to first order in the non-linearities, [9,61] (3.14) The factor 4/3 in front of r S /r is due to the extra contribution coming from the exchange of the helicity-0 graviton, and gives rise to the vDVZ discontiuity. In contrast, no vDVZ discontinuity is present in eq. (3.13). Furthermore the correction terms are crucially different. In eq. (3.14) the correction explodes at low r, i.e. for r below the Vainshtein radius r V = (GM/m 4 ) 1/5 . In eq. (3.13), in contrast, the correction becomes smaller and smaller as r decreases, and perturbation theory is valid at all scales r m −1 , until we arrive at r ∼ r S where also GR becomes non-linear. In summary: • In the nonlocal theory defined by eq. (1.1) there is no vDVZ discontinuity. This confirms the result that was found in [26] by expanding over flat space.
• The linearized expansion is fully under control for all distances r S r m −1 . Contrary to (local) massive gravity theories, the classical theory stays linear for all distances down to r ∼ r S , where eventually also the usual GR non-linearities show up. For r S r m −1 the corrections to GR are actually smaller and smaller as r decreases, the classical theory never becomes strongly coupled, and recovers all successes of GR at the solar system and lab scales.
The m 2 expansion discussed in this section allowed us to obtain perturbatively the solution in the region r m −1 . As r approaches m −1 , the corrections become of order one and the small-m expansion breaks down. Furthermore we have written explicitly the correction terms only in the limit r r S . However, this is not due to an intrinsic limitation JHEP08(2014)029 of the perturbative expansion but is only done for simplicity, since the full expressions are somewhat long. In any case, we found that the corrections are proportional to m 2 r 2 , and becomes smaller and smaller as r decreases toward the horizon, so the terms that we have omitted are in fact negligible even close to the horizon (as we will also check in section 5 comparing the perturbative solutions (3.13) to the result of the numerical integration).
Thus the results of this section are a good approximation to the exact solution in the whole range r S < ∼ r m −1 .
4 Solution for r r S . The Newtonian limit The solution in the region r r S , with no limitation of the parameter mr, can be obtained with a different expansion, namely considering the effect of the source as a perturbation of Minkowski space, adapting the standard analysis performed in GR to recover the Newtonian limit. This will allow us to obtain analytically the solution in the region mr > ∼ 1, which is not accessible to the low-m expansion. Furthermore, in the region r S r m −1 both the low-mass and the Newtonian expansions are valid, and therefore we can match the solutions. We will then confirm the validity of these expansions with the numerical integration in section 5.
Thus, in this section we start from a backgroundḡ µν = η µν ,Ū = 0 andS µ = 0, and T µν = 0, and we perturb it adding the energy-momentum tensor of a localized source. We limit ourselves to a static non-relativistic source, in which case δT 00 = ρ, while δT 0i and δT ij vanish. We are interested in the scalar perturbations so, using the Newtonian gauge, the perturbed metric can be written as (4.1) We also expand the auxiliary fields as U =Ū + δU , S µ =S µ + δS µ . Since the background valuesŪ =S µ = 0, we simply write the perturbations as U and S µ , keeping however in mind that they are first-order quantities, just as Φ and Ψ. Furthermore, for a static source we necessarily have S 0 = 0, as before, both for the background and the perturbation. The vector S i can instead be decomposed as usual as S i = S T i + ∂ i S where S T i is a transverse vector, ∂ i S T i = 0, which only contributes to vector perturbations, while S is a scalar. Since we are studying the scalar sector we only retain S, and we write S i = ∂ i S. Thus, a static source induces scalar perturbations which are described by the four functions Ψ, Φ, U and S. Observe that we do not need to restrict to spherical symmetry. The vanishing of S 0 is just a consequence of the fact that the source is static, so in this problem nothing depend on time, and ∂ 0 U = 0, so eq. (2.6) with ν = 0 is a homogeneous equation that has the solution S 0 = 0.
Observe also that the radial coordinate used in this section is different from that used in section 3, since we are in a different gauge. In fact, if we linearize eq. (2.7) we get This expression is not in the Newtonian gauge. In the Newtonian gauge the factor (1 + 2Ψ) multiplies the whole term dx 2 , while in eq. (4.2) the factor (1 + 2β) only multiplies dr 2 .

JHEP08(2014)029
For clarity, we will continue to denote by r the radial coordinate of the metric (2.7), while we denote by r N the radial coordinate in the Newtonian gauge, so eq. (4.1) reads As discussed above, in the region r S r m −1 both the low-mass expansion of the previous section and the Newtonian expansion of this section hold, and we can therefore match the results. In order to perform the matching, we will however need the relation between the two coordinates, as well as between α, β and Ψ, Φ in this regime. This is easily found observing that, when r S r m −1 , the full Schwarzschild-like metric (2.7) can be linearized and written as in eq. (4.2). We then rewrite the metric (4.3) as Comparing with eq. (4.2) we see that r = (1 + Φ)r N and (1 + Φ)dr N /dr = 1 + β. Inserting here r N = (1 + Φ) −1 r we get β = −rΦ /(1 + Φ) which, to the linearized order at which we are working, is equivalent to β = −rΦ (and, since rΦ is already a first-order quantity, we do not need to distinguish r from r N here). In summary, in the overlapping region r S r m −1 we can compare the results of the two approaches, using the relations After these preliminary remarks, we perform the actual linearization. We write eq. (1.1) in the form Linearizing the (00) components of eq. (4.6) and setting to zero all time dependences, as appropriate for a static source, we get Observe that the Laplacian is in principle with respect to the coordinate r N , but since all quantities in eq. (4.8) are first-order in the perturbations, we can equivalently use r. The same will be true for all other equations below. The (0i) component of eq. (4.6) vanishes identically on time-independent perturbations. The linearization of the (ij) equation gives, setting again to zero all time derivatives, where, since we are working to linearized order over Minkowski space, we are free to write all spatial indices as lower indices, and we have used the fact that, for a Newtonian source,

JHEP08(2014)029
T 0i can be neglected. Applying to this equation the projector (∇ −2 ∂ i ∂ j − 1 3 δ ij ) to obtain the traceless part, we get One might be tempted to rewrite this equations as Φ + Ψ − (m 2 /3)S = 0, but this would not be correct. In general, if a function f satisfies ∇ 2 f = 0 over all of space, and we further impose the boundary condition that f vanishes at infinity, then f = 0. However, the equations that we are writing in this section are only valid for r r S . Of course, from the fact that a function f satisfies ∇ 2 f = 0 at large r we cannot conclude that f itself is identically zero at large r. Indeed, any function that, at large r, approaches the form f (r) = c 0 + c 1 r S /r satisfies ∇ 2 f = 0 at large r. In our problem, the constant term c 0 is eliminated requiring that the functions Φ, Ψ and S vanish at infinity. We remain however with the possibility of a 1/r term. Thus, eq. (4.10) only implies that, at r r S , for some constant c 1 , which can be determined by matching the solution with those found in section 3 for r m −1 , as we will do below. Taking the trace of eq. (4.9) and combining it with eqs. (4.10) and (4.8) we get   This is an inhomogeneous Helmholtz equation, and we can solve it writing where

JHEP08(2014)029
The Green's function of the inhomogeneous Helmholtz equation is well known. Writing G(r) = −[1/(4πr)]f (r) one finds with β arbitrary. The corresponding solution for U is In the r r S limit this becomes More generally, even if ρ(x) is not a Dirac delta, at distances r much larger than the source size we can write |x − x | r − x ·n, wheren = x/r, so cos(m|x − x |) cos(mr − mx ·n). For m = O(H 0 ), all over the source m|x | is negligibly small with respect to one (and not just with respect to mr) and we can replace cos(m|x − x |) by cos(mr). Therefore, at large distances the coefficient of the 1/r term for a generic ρ(x) is the same as for ρ(x) = M δ (3) (x), just as in GR. The appropriate Green's function, and therefore the value of β, is fixed by the boundary conditions. In most problems in which the inhomogeneous Helmholtz equation appears, the Green's function is fixed imposing a no-incoming wave boundary condition at infinity, which selects G(r) = −e imr /(4πr), i.e. β = i. However, such a boundary condition is not appropriate to our problem, since U (r) is real. In our case, for an extended source, β must rather be fixed by matching this large distance solution to the solution in the inner source region, as we will discuss below.
We can now plug this solution for U into eq. (4.8). Using for simplicity ρ(x) = M δ (3) (x), we get where c Φ = −(c 2 + 2c Ψ ) is a second independent constant. To compare with the functions A(r) and B(r) of the previous section we use the fact that, in the region r S r, A(r) = 1 + 2α(r) and B(r) = 1 + 2β(r). Using eq. (4.5) we therefore have A(r) = 1 + 2Ψ and B(r) = 1 − 2rΦ , so Matching these expression with the solution (3.13), which is valid for r m −1 , we get β = 0, c Ψ = −4/3 and c Φ = 2/3. On the other side, comparing the terms m 2 r 2 , allows us to fix u 1 in the solution (3.13), and we get u 1 = 1. The latter result could have also been derived more directly matching the small mr limit of eq. (4.22) to the large r/r S limit of eq. (3.7).
In conclusion, plugging the value of these constant into the full solutions (4.26) and (4.27) we find that, in the Newtonian limit, ds 2 = −A(r)dt 2 + B(r)dr 2 + r 2 (dθ 2 + sin 2 θ dφ 2 ) , (4.30) where while the auxiliary field U = − −1 R is given by U (r) = r S r cos mr .

Numerical integration
We now study the equations numerically, in order to confirm the above analytic results. In the numerical analysis it is convenient to trade U for a field W defined by U = W + 2α.
The advantage of this transformation is that now α disappeared from eq. (5.2). To integrate the equations we need to assign the initial conditions. To this purpose, we take advantage of the fact that we know the zero-th order solution is quite close to the exact solution in the region r S r m −1 . We therefore choose a value r in in this region, and we assign α(r in ), β(r in ), U (r in ), U (r in ) V (r in ) and V (r in ) using the zero-th solution given by eqs.  result) while the Newtonian result becomes less accurate. These plots confirm that the theory never becomes non-linear in the region r S r m −1 . The exact numerical solution follows the analytic solution obtained in an expansion in powers of mr, until mr becomes of order one. There is no Vainshtein radius r V r S below which the low-mass expansion fails.
It is also interesting to study the stability under perturbations of the solution that we have found. Equation (4.20) and the discussion below eq. (4.22) show that, if we perturb the source replacing ρ → ρ + δρ, while still preserving the fact that the source has compact support, no instability develops, and the only effect of δρ is to replace the source mass M by the corresponding value M + δM , just as in GR. Concerning the near-horizon region of a BH solution, our analytic and numerical results indicate that the corrections to the Schwarzschild solution near the horizon are O(m 2 r 2 S ), which for m ∼ H 0 is negligibly small, so we do not expect any instability to develop in the BH quasi-normal modes.

Degrees of freedom of the nonlocal theory
To better understand the meaning of the results obtained, it is useful to discuss what are the radiative and non-radiative degrees of freedom of the theory. This issue has already been examined in detail in refs. [26,50] (see also [57] for a related analysis in a similar nonlocal model). However, we find useful to summarize here the main results discussed in the above papers, and compare them with what we have learned above.
Counting the propagating degrees of freedom in a nonlocal theory (or even in a local theory when we transform to nonlocal variables) involves some subtleties, of which one must be aware in order to get correct results. The simplest example of what can go wrong with nonlocal transformations is provided by a theory in which, among other fields, also appears a scalar field φ that satisfies a Poisson equation ∇ 2 φ = ρ [50]. This field is clearly non-radiative. If, in the classical equation, we set the source ρ = 0, we simply have φ = 0. There are no associated plane waves freely propagating in empty space associated to this field, and at the quantum level there are no creation and annihilation operators associated to it. However, if we define a new fieldφ fromφ = −1 φ, the original Poisson equation

JHEP08(2014)029
can be rewritten as φ = ∇ −2 ρ ≡ρ , (6.1) so nowφ looks like a propagating degree of freedom. However, for ρ = 0 our original equation ∇ 2 φ = ρ only has the solution φ = 0. If we want to rewrite it in terms ofφ without introducing spurious degrees of freedom we must therefore supplement eq. (6.1) with the condition that, when ρ = 0,φ = 0. In other words, the solutions of the associated homogeneous equation φ = 0 must be discarded (or, more generally, is uniquely fixed by the boundary conditions of the problem). Correspondingly, the coefficients a k , a * k of the general plane-wave solution of the equations φ = 0 cannot be considered as free parameters that, upon quantization, give rise to the creation and annihilation operators of the quantum theory, and there are no propagating quanta associated toφ.
A similar situation also appears in general relativity. Consider GR linearized over flat space, g µν = η µν + h µν , and decompose as usual h µν as where i and β i are transverse vectors, ∂ i β i = ∂ i i = 0, and h TT ij is transverse and traceless, ∂ j h TT ij = 0 and δ ij h TT ij = 0. With these variables we can form the gauge-invariant combinations Φ = −φ − (1/6)∇ 2 λ and Ψ = −ψ +γ − (1/2)λ, which describe two degrees of freedom in the scalar sector, and the gauge-invariant transverse vector Ξ i = β i − (1/2)˙ i , which describes two degrees of freedom in the vector sector. These gauge-invariant quantities are the usual Bardeen's variables specialized to flat space. The remaining degrees of freedom (two in the scalar sector and two in the vector sector) can be set to zero with a gauge transformation, while in the helicity-2 sector h TT ij is gauge invariant. We have therefore split the 10 components of h µν into four pure gauge modes and six physical (i.e. gauge-invariant) degrees of freedom. Of course, not all these six degrees of freedom are radiative. To see this, we decompose similarly the energy-momentum tensor as where ∂ i σ i = 0, ∂ i Σ i = 0, ∂ i σ ij = 0 and δ ij σ ij = 0. The linearized equations of motion can then be written as [69] ∇ 2 Φ = −4πGρ , (6.6) This shows that only the tensor perturbations obey a wave equation and are therefore radiative. The gauge-invariant scalar and vector perturbations obey a Poisson equation,

JHEP08(2014)029
and therefore represent physical but non-radiative degrees of freedom. Now, compare this result with that obtained decomposing h µν as where h TT µν is transverse and traceless with respect to the Lorentz indices, ∂ µ h TT µν = 0, η µν h TT µν = 0, and therefore has five independent components. The 10 components of the metric perturbations are therefore split into the five components of h TT µν , the four components of µ , plus the scalar s. Under a linearized diffeomorphism h µν → h µν − (∂ µ ξ ν + ∂ ν ξ µ ) we have µ → µ − ξ µ while the tensor h TT µν and the scalar s are gauge invariant. We now plug this decomposition into the quadratic Einstein-Hilbert action, where E µν,ρσ is the Lichnerowicz operator (and we rescaled h µν → κh µν , where κ = (32πG) 1/2 , in order to have a canonically normalized kinetic term). The result is Performing the same decomposition in the energy-momentum tensor, the interaction term can be written as At first sight this result is surprising, because it seems to suggest that the five components of the transverse-traceless tensor h TT µν and the scalar s are all radiative fields. Note that these degrees of freedom are gauge invariant, so they cannot be gauged away. 4 Furthermore, according to eq. (6.12) the scalar s should be a ghost! Of course these conclusions are wrong, and the correct conclusion is the one drawn from eqs. (6.6)-(6.9), namely that in GR there are two radiative and four non-radiative degrees of freedom. What went wrong is the following. Eqs. (6.6)-(6.9) can be inverted, to give Φ, Ψ, etc. in terms of the original field h µν . The inversion involves the inverse Laplacian, and is therefore nonlocal in space, 4 Observe also that we have not used linearized gauge invariance to set µ = 0 in the action (in which case, one should have been worried that we might have lost the equations obtained performing the variation with respect to µ ). Rather, µ disappears automatically from the action, so there is no equation of motion associated to it. The fact that µ disappears automatically from the action is just a property of E µν,ρσ when applied to a tensor of the form ∂ρ σ . Of course, this property of the E µν,ρσ is just what guarantees the invariance of the quadratic action under linearized gauge transformations hµν → hµν − (∂µξν + ∂ν ξµ).

JHEP08(2014)029
but is local in time. In particular this means that there is a one-to-one correspondence between the initial condition assigned on h µν on a given time slice, and those assigned on Φ, Ψ, etc. on the same time slice. In contrast, when we invert eq. (6.10), we find that the inversion involves the inverse d'Alembertian, and is therefore non-local even in time. Thus, this one-to-one correspondence on the initial conditions is lost, and the counting of degrees of freedom goes wrong. In particular, the inversion of eq. (6.10) gives 16) and the relation between s and the Bardeen variables Φ and Ψ is We see that the situation is exactly the same as that illustrated by eq. (6.1): the nonradiative field Φ − Ψ is transformed into an apparently radiative field s by the nonlocal relation (6.17), which involves −1 . The bottom-line is that, if in GR we wish to use the decomposition (6.10) and the fields s and h TT µν we can do it, provided that we supplement eq. (6.15) with the condition that, when ρ = σ = 0, we must have s = 0, i.e. we must discard the solution of the homogeneous equation s = 0, since when ρ = σ = 0 we see from eqs. (6.6) and (6.7) that Φ = Ψ = 0 (and similarly for the component h TT µν with helicities 0 and ±1). So, again, there are no creation and annihilation operators associated to these fields in the quantum theory, and they cannot appear in external lines nor in loops. 5 Having understood these simple but important points in the familiar context of GR, we are now well-armed for understanding the situation in the nonlocal theory. As we discussed in section 2, using the auxiliary field U = − −1 R as well as the auxiliary four-vector field S µ that enters in the extraction of the transverse part in eq. (1.2), eq. (1.1) can be rewritten as together with U = − −1 R and 19) which is obtained taking the divergence of eq. (1.2). Naively, one would say that U = − −1 R is equivalent to U = −R, and therefore the original nonlocal model can be written as a set of local equations for g µν , U and S µ , given by eqs. (6.18) and (6.19) together with U = −R. However, it is just in this "localization" step then one is introducing spurious solutions. The point is that, as we already discussed in section 3, an equation such as U = −R is solved by eq. G(x; x ) is any a Green's function of the operator, and the definition of the nonlocal model is completed only once we have fully specified the definition of −1 , by specifying U hom (x). The specific choice of U hom (x) can depend on the specific problem that we are studying, e.g. a different choice will appropriate for the study of a static solution, or for the study of FRW solutions, simply because the boundary conditions of these problems are different. In any case, in any given problem a choice must be performed, and this fixes the homogeneous solution. The solutions of the equations in the local formulation are also solutions of the original integro-differential equations only for this specific choice of U hom (x), and all other solutions are spurious. Thus the solution of the homogeneous equation U = 0 cannot be interpreted as a free field, which is expanded in plane wave, and whose coefficients a k , a * k are then interpreted as creation and annihilation operators in the corresponding quantum theory. Indeed, we have seen in section 4 that in the specific problem of a static spherically symmetric solution, the homogeneous solutions is fixed to the expression given in eq. (4.33). This is similar to what happens for the scalar metric perturbations Φ and Ψ, which are also fixed by the boundary conditions. In contrast, for the tensor perturbations, on any given background we always have the freedom to add gravitational waves freely propagating to infinity, i.e. to add to the solution of eq. (6.9) an arbitrary solution of the homogeneous equation h TT ij = 0. A similar choice must be made when we define what it means exactly to extract the transverse part in eq. (1.1). Indeed, the solution of the homogeneous equation is not a free radiative degree of freedom, but it is part of the definition of extraction of the transverse part. In this case the most obvious definition is to set this homogeneous solution to zero, corresponding to the fact that, if S µν = 0, there is no transverse part to extract, and S T µν = 0, so we define the extraction of the transverse part so that, when S µν = 0, also S µ = 0.
Having realized that the fields U and S µ do not carry radiative degrees of freedom, it becomes clear that the content of eq. (6.18), as far as radiative degrees of freedom are concerned, is the same as in GR, namely two massless graviton states with helicities ±2, which are now coupled also to extra non-radiative fields. This can be checked explicitly by looking at the linearized version of the theory. Linearizing eq. (1.1) over Minkowski space we get where and is now the flat-space d'Alembertian. Let us examine first the scalar sector. In the linearized limit we can put the theory in a local form in an even simpler way, namely introducing two auxiliary scalar fields U = − −1 R and S = − −1 U . Then, writing again the metric in terms of the Bardeen variables and the energy-momentum tensor as in eqs. (6.4) and (6.4), the equations of the linearized theory in the scalar sector can be rewritten as [57]  The equations for U is just the linearization of U = −R and, as explained above, we must discard its homogeneous solution. Hence, it does not describe radiative degrees of freedom. The same holds for S, which plays the role that in the full nonlinear theory is played by the four-vector S µ , and is defined so that S = 0 when U = 0. The Bardeen variables Φ and Ψ still satisfy Poisson equations, as in GR, and therefore remain non-radiative. This should be contrasted with what happens in massive gravity with a Fierz-Pauli mass term, where Φ becomes radiative and satisfies a massive Klein-Gordon equation ( − m 2 )Φ = 0 [69][70][71]. Thus, in the nonlocal theory there is no radiative degree of freedom in the scalar sector, while the vector and tensor sector are obviously not affected by the presence of U and S. In particular, the graviton remains massless, as we see from eqs. (1.8) and (1.9). Indeed, the first term in eq. (1.9) describes the matter-matter interaction induced by a massless helicity-2 field, while the second term contributes to the matter-matter interactioñ These two terms are induced by the massless field S and by the massive field U , respectively. If U were a radiative field, its contribution would be ghost-like, and one should worry about vacuum decay in the quantum theory. However, we have seen that U is not radiative, and at the quantum level there are no creation and annihilation operators associated to it, and no quantum vacuum instability. 6 Another way to understand this result, again discussed in [50], is to observe that, linearizing around flat space, g µν = η µν + h µν , we have R = R (1) + O(h 2 ), where R (1) = ∂ µ ∂ ν (h µν − η µν h). Therefore in the linearized theory Comparing with eq. (6.16) we see that, in the linearized theory, U is the same as the non-radiative field s. The quadratic Lagrangian corresponding to the linearized equa- 6 Of course, such a field can induce instabilities at the classical level. However, while the quantum vacuum decay would be a disaster for the consistency of the theory, classical instabilities must be examined on a case-by-case basis, and can in fact even be welcome. This is particularly true in a cosmological setting, where the emergence of a phase of accelerated expansion is in a sense a classical instability. Indeed, the result of [26] show that, at the background level, the cosmological evolution of this theory is perfectly viable.
In [52] we will examine the cosmological perturbations in these nonlocal models, and we will see again that they are perfectly viable, and in agreement with the observation. Similarly, the results of the present paper showed that no dangerous instability develops in static solutions, and that GR is smoothly recovered at r m −1 .

JHEP08(2014)029
tion (6.20) is h µν E µν,ρσ h ρσ − 1 3 m 2 (P µν h µν ) 2 , (6. 28) and we see that the term proportional to m 2 is just a mass term for s = P µν h µν . Writing the metric as in eq. (6.10), instead of eq. (6.12) we now obtain [50] S (2) = 1 2 We see that the effect of the non-local term at the linearized level can be described as follows. In usual massless GR we have six physical (i.e. gauge-invariant) degrees of freedom, that can be described by the five degrees of freedom of the transverse-traceless tensor h TT µν , plus the scalar s. Out of them, only two are radiative, i.e. the helicity ±2 components of h TT µν . The helicity-(±1) components of h TT µν , the helicity-0 component of h TT µν and the scalar s are all non-radiative, and all these six fields are massless. In the matter-matter interaction, the contribution from s has the opposite sign, compared to the one coming from the helicity-0 component of h TT µν , and the two cancel, while the helicity ±1 components of h TT µν , are coupled to ∂ µ T µν , which vanish, and therefore do not contribute. Thus, we only remain with the contribution from the helicity-(±2) components. In the nonlocal theory, the field s remains non-radiative, but becomes massive. Therefore, the cancelation with the helicity-0 component of h TT µν is only approximate, and only holds for |k 2 | m 2 . Thus, for m ∼ H 0 , well inside the horizon we recover GR and there is no vDVZ discontinuity, as indeed the computation of the previous sections showed explicitly. In contrast, at cosmological scales, there are departures from GR, but no new propagating degree of freedom.
It is also interesting to compare the above discussion with the results of refs. [72,73], where the author performed a very general analysis of ghost-free modified gravitational actions, linearized over Minkowski space, by including the most general form factors depending on the operator, i.e. terms such as h µν a( )h µν , h σ µ b( )∂ σ ∂ ν h µν , etc. The propagator can then be found in full generality, and one can impose conditions on the form factors a( ), b( ), etc. such that no ghost-like pole appears, and furthermore the UV behavior is improved. The analysis in [72,73] was mostly tuned toward the UV behavior, and therefore one is mostly concerned with positive powers of the operator. What we learn from our discussion in this section is that, when we apply this analysis to the IR, where non-local operators such as −1 become relevant, an apparent ghost-like pole in the propagator is not yet necessarily a sign of a trouble, since it could simply correspond to a non-propagating degree of freedom.
Observe also that, with our (−, +, +, +) signature, the operator ( + m 2 ) that appears in the linearized equation of motion for s (which is just eq. (6.24), given that at the linearized level s = U ) corresponds to a dispersion relation k 2 0 = −m 2 + k 2 . Therefore, static solution do not decay at large distances with a Yukawa suppression r −1 exp{−mr}, but are instead oscillatory, r −1 cos(mr), as indeed we found in section 4.

JHEP08(2014)029 7 Conclusions
The analytic and numerical results discussed in this paper show that, in the nonlocal theory defined by eq. (1.1), the linearized expansion is valid for all distances r in the range r S r m −1 , and in this region the corrections to GR are of the form 1 + O(m 2 r 2 ). This is in sharp contrast with what typically happens in local theories of massive gravity, where the linear expansion breaks down below a Vainshtein radius r V which is parametrically larger than r S , and which diverges as m → 0. In local massive gravity theories (whether in a Fierz-Pauli or dRGT form) this breakdown of linearity is necessary for their observational viability, since these theories have a vDVZ discontinuity at large distance. Without such a breakdown of linearity, this discontinuity would persist down to the solar system scale, and then the theory would be ruled out. In contrast, the nonlocal theory (1.1) has no vDVZ discontinuity, and it remains linear down to the near-horizon region. Therefore, all successes of GR at the solar system and lab scales are automatically recovered. This is an important consistency check of the nonlocal theory which, together with its cosmological properties discussed in [26,51], makes it a interesting candidate for a dynamical explanation of dark energy.
Furthermore, we have determined the behavior of the solution in the region (r r S and mr generic) using a Newtonian expansion. Equations (4.31) and (4.32) provide an analytic expression for the modifications of the static Newtonian forces at distances of order m −1 in the nonlocal model that we have studied, and could have potential applications in the study of structure formation at large scales in such a model.