text:
# 1 Introduction

## 1 Introduction

Over the past few years there has been considerable discussion about theories with light massive gravitons in their spectra. These kinds of theories first arose in the context of brane-world models  and at the linearized approximation they predicted modifications of Newton’s constant or even of Newton’s law itself at cosmological scales. The evidence for a dark energy component in our universe and the associated cosmic acceleration  made the study of these models rather topical.

Indeed, generalizing these theories at the non-linear level , revealed that they naturally give rise to a period of late time acceleration of the universe . The non-linear analogue of a collection of light massive gravitons is a theory with many interacting metrics, one of which only couples to the matter fields of our universe . Depending on the particular form of the coupling of the several metrics, the resulting acceleration could have interesting testable differences from a (scalar field) quintessence model, as for example anisotropic features .

Despite the interest in these theories regarding cosmology, there are potentially dangerous issues associated with the presence of extra polarization states of the massive gravitons . In particular, the massive gravitons have a scalar-like polarization state whose coupling to matter does not depend on the mass of the graviton. This scalar coupling to matter is formally analogous (in the linearized approximation) to a Jordan-Fierz-Brans-Dicke coupling $`\omega =0`$. This coupling modifies (by a factor $`(2\omega +3)/(2\omega +4)=3/4`$) the usual general relativistic relation between interaction of matter and light. If General Relativity (GR) is modified by a mass term, or if GR is augmented by the addition of an extra massive graviton which dominates444In the case where the massive graviton component were sufficiently subdominant there would be no observational discrepancy, but also no phenomenological interest for considering such a massive component. the matter couplings, the discrepancy between the theory of massive gravitons and GR in e.g. the bending of light by the Sun, would be at the $`25\%`$ level, when the current observational accuracy is better that one part in $`10^4`$.

In addition to this (very serious) experimental difficulty, it has also been shown that massive gravity has serious theoretical defects . As soon as one goes beyond the linearized level, massive gravity has six degrees of freedom, instead of the expected five ($`2s+1`$, for a massive spin $`s=2`$). This sixth degree of freedom is problematic both because it represents a jump, with respect to the linearized theory, in the number of degrees of freedom, and because its energy has no lower bound. These arguments seem to exclude the existence of dominantly contributing massive gravitons, however small their mass might be. It seems, in that case, that gravity is a unique example of an “isolated” theory whose massive deformation can be excluded to infinite accuracy.

One could question several points in the arguments of  recalled above. The first one is that these arguments are based on the study of models containing an explicit breaking of the (linearized) gauge invariance $`\delta h_{\mu \nu }=_\mu \xi _\nu +_\nu \xi _\mu `$ by a Pauli-Fierz mass term for the light graviton(s), which reads:

$$\frac{m^2}{4}(h_{\mu \nu }h^{\mu \nu }h^2)$$
(1)
One could hope that the peculiar result of the discontinuity is linked to this explicit breaking and that continuity in predictions might be restored if the mass for the spin-2 field were generated spontaneously. For example, if light gravitons are generated by the compactification of a higher dimensional theory, one could imagine a higher dimensional mechanism which would restore the continuity in the theory555For an alternative way of generating spontaneously masses for gravitons see .. We will not explore such possibilities in this paper, but study instead theories with explicit mass terms with the aim of clarifying whether such theories are indeed physically sick.

A second potentially weak point of the above argument is that the effects of Pauli-Fierz mass terms were considered only around flat space backgrounds. One could imagine that the discontinuity that was found is a peculiarity of just the flat background and that if one considered a background with curvature, some of the difficulties might be evaded. Indeed, it was found that for instance in constant curvature backgrounds \[$`(A)dS`$ spaces\], the extra polarizations of the massive gravitons have a coupling $`m/H`$ where $`m`$ is the mass of the graviton and $`H`$ the “Hubble” constant of the $`(A)dS`$ space . In that case, the predictions of the massive theory were indistinguishable666Note that at the quantum level the discontinuity formally reappears , but as a quantum effect is suppressed and unobservable in any conceivable near future experiment. from the massless one as long as $`mH`$. This, however, prevents the massive gravitons from being cosmologically interesting since their Compton wavelength is much bigger than the “Hubble” radius of the “bare” cosmological constant.

In constant curvature backgrounds the parameter controlling the smoothness of the limit is $`m^2/`$ where $`H^2`$ is the curvature of the (A)dS space. One might conjecture that the same is true for any curved space and that the smoothness is controlled by a parameter $`m^2/`$, where $``$ is some curvature invariant. In the case of generic vacuum Einstein backgrounds (with vanishing scalar curvature and Ricci tensor) $``$ might be defined as some invariant of the Riemann tensor (such as the square root of the square of the Riemann tensor). In the case of a Schwarzschild background, $`R_S/R^3`$, where $`R_S`$ is the Schwarzschild radius of the source. Then a smooth limit might be expected to exist for distances $`RR_p`$, where the upper bound $`R_p`$ of the expected interval of existence of smooth solutions is:

$$R_p=\frac{R_S}{(mR_S)^{2/3}}=(\lambda _m^2R_S)^{1/3}$$
(2)
where $`\lambda _mm^1`$ denotes the Compton wavelength associated to the mass $`m`$.

A similar conjecture was made a long time ago by Vainshtein  based, however, on a different argument. That paper made two basic points: (i) it questioned the use of perturbation theory (and especially of the linearized approximation) in the derivation of observable consequences of massive gravity by showing that non-linear effects, proportional to a negative power of $`m^2`$, are important in a wide domain around the source, and (ii) it conjectured for the first time that the resummation of non-linear effects might actually restore continuity in a domain near the source. More precisely, this paper sketched the construction of solutions which can be expressed, at least within some intermediate range of distances $`R_SR_{}RR_V`$, as a series in positive powers of $`m^2`$ and whose leading term coincides with the one of GR. The interval where this expansion was constructed was bounded on the left by the radius of the star, $`R_{}`$, and on the right by the length scale:

$$R_V=\frac{R_S}{(mR_S)^{4/5}}=(\lambda _m^4R_S)^{1/5}$$
(3)
which differs from the one obtained by the conjecture made in the previous paragraph.

In the simple example of a massive graviton with Compton wavelength of the order of the universe’s horizon and of a source as massive as the Sun, the distance scale (3) is much bigger than the distance scales on which relativistic predictions of GR are tested. Thus, if the conjecture were correct, the massive graviton proposal would not be excluded. \[Actually, the same holds if continuity is restored on the scale (2), though $`R_pR_V`$.\] The arguments of  were repeated, with more intermediate details, in the recent paper of . This paper also provided a cosmological analogue of the continuity-discontinuity interplay in the context of the brane-induced gravity models (see also  for a cosmic string example for the same kind of brane models).

However, the arguments of  are questionable, and one of the main thrusts of the present paper will be to show that they actually are inconclusive. The questionable aspects concern both the local existence of continuous solutions, their global relevance, and even their global existence. For instance, we shall show that the action principle postulated in  does not admit , in the intermediate range $`R_SR_{}RR_V`$, local solutions of the type conjectured in  for restoring continuity. Moreover, it is not at all clear if such local solutions, when they do exist (for some other action principles), can be globally continued to the asymptotically flat ones of the linearized theory, a concern first raised in  and confirmed later in . Actually, it is not even clear if these solutions are indeed solutions, i.e. if they do not contain naked singularities at some finite distance.

Let us note here that there have been several papers  where the Vainshtein idea was investigated in brane-induced gravity for static spherically symmetric sources and novel phenomenology was proposed  777See  for a demonstration of the discontinuity at the linearized level for spherically symmetric sources in brane-induced gravity models. For the same models, see also  for special spherically symmetric solutions.. In our view, the approximations made in these papers, although plausible, do not address the global properties of the solutions. We will not explore the specific brane-induced gravity model in this paper.

In the present paper we will discuss, within the context of purely massive gravity, the claims of  as well as the more expanded version of them that appeared in . We will analyse in detail the procedure that was used in  to obtain the results, shedding light to some obscure points, and repeating the same calculation for several mass terms. We will show that the method has serious limitations and that there are contradictions between the statements made in  and in . Indeed, we shall point out that the field equations studied in  do not derive from the action principle postulated in . Actually, they cannot be derived from any action principle. We shall show that, if one starts from the simple action principle written in , there exist no expansions of the type postulated in  as a way to restore continuity. However, we shall show that such expansions do exist for other action principles (and do exist for the field equations written down in , after correction of some misprints). We will also see that the scale $`R_V`$ up to which the expansion of  makes sense, is not universal for different potential terms. We will then study (numerically) the global aspects of such solutions. Our numerical simulations strongly suggest that: (i) the Vainshtein-type local (approximate) solutions do not match, as they were supposed to, the asymptotically flat approximate solutions generated by normal perturbation theory (in agreement with ), and in fact that (ii) the asymptotically flat solutions cannot be extended inwards to globally regular solutions. Even if we forget about the issue of mass-continuity, we find the striking result that all asymptotically flat solutions run into naked singularities as the radial coordinate $`R`$ decreases.

However, our conclusions will not all be negative. Indeed, we shall explicitly construct some globally regular solutions of massive gravity which are continuously connected, when $`m0`$, to GR solutions, and which are phenomenologically consistent with experimental tests of relativistic gravity. The global solutions we shall (numerically) construct are generalizations (to the inclusion of a central source) of the black-hole-type solutions of massive gravity constructed long ago by Salam and Strathdee . \[The latter solutions were generalized by Isham and Storey  for a particular class of bigravity theories .\] The main difference between these solutions and the ones studied in Refs.  is the behaviour at infinity. The solutions are not required to be asymptotically flat, but instead to match to a cosmological solution of massive gravity (which is de Sitter, in the case at hand). In fact, a general argument of  (see section 4 there) has shown how to construct (at least for a limited time) general classes of solutions which are continuously connected, when $`m0`$, to GR solutions representing local gravitating systems (such as the solar system), embedded in some global cosmological background. Another big difference with the solutions studied in , and in the first part of this paper, is that these cosmologically-matched solutions are of a special “symmetry-breaking” type (see below).

As we shall discuss in our conclusions, our results leave open several important issues which must be tackled before a firm conclusion can be reached concerning the physical consistency (or inconsistency) of massive gravity theories.

## 2 Action and potentials for massive gravity

Our starting point is a generic action for (four dimensional) massive gravity. This can be obtained, as it was shown in , from a four dimensional bigravity action if we send the gravitational constant associated with the second metric (the one not coupled to “our world”) to zero. Then the corresponding metric gets (formally) frozen, i.e. it becomes a non-dynamical Einstein space background. We assume that this non-dynamical metric is flat, and we accordingly denote it by $`𝐟`$. \[We do not necessarily assume that the metric $`𝐟`$ is written in Lorentzian coordinates; i.e. the components $`f_{\mu \nu }(x)`$ are not necessarily assumed to be simply $`\eta _{\mu \nu }`$.\] The only remaining dynamical metric is denoted $`𝐠`$. Then the effective action for the purely massive gravity is:

$$𝒮=\frac{1}{16\pi G}d^4x\left(\sqrt{g}R[𝐠]\frac{m^2}{4}𝒱(𝐟^1𝐠)\right)+𝒮_{\mathrm{matt}}[𝐠]$$
(4)
where $`G`$ is the “massive” version of Newton’s constant and $`m`$ a mass parameter, which for the specific normalization of the potentials $`𝒱`$ which we will consider, is the canonically normalized Pauli-Fierz mass of the graviton. In the following, we explore the following possible potential terms:

$`𝒱^{(1)}=\sqrt{g}\left\{\mathrm{tr}[((𝐠𝐟)𝐟^1)^2](\mathrm{tr}[(𝐠𝐟)𝐟^1])^2\right\}`$ (5)
$`𝒱^{(2)}=\sqrt{f}\left\{\mathrm{tr}[((𝐠𝐟)𝐟^1)^2](\mathrm{tr}[(𝐠𝐟)𝐟^1])^2\right\}`$ (6)
$`𝒱^{(3)}=\sqrt{g}\left\{\mathrm{tr}[((𝐠^1𝐟^1)𝐟)^2](\mathrm{tr}[(𝐠^1𝐟^1)𝐟])^2\right\}`$ (7)
$`𝒱^{(4)}=\sqrt{f}\left\{\mathrm{tr}[((𝐠^1𝐟^1)𝐟)^2](\mathrm{tr}[(𝐠^1𝐟^1)𝐟])^2\right\}`$ (8)
Note that the above quantities are not scalars, but scalar densities. Comparing with the notation of , $`𝒱=(\mathrm{density})V`$ where $`V`$ is a scalar. All of the four above potentials are of the Pauli-Fierz type near $`𝐠𝐟`$. Note that the first one was written to be the starting point in .

Each contribution $`𝒮_x`$ to the action gives rise to a corresponding contribution to the energy-momentum tensor: $`T_{\mu \nu }^x(2/\sqrt{g})\delta 𝒮_x/\delta g^{\mu \nu }`$. The “gravitational” energy-momentum tensors arising from the mass terms will be denoted $`T_{\mu \nu }^{(g)}`$. For the specific mass terms above, they are respectively:

$`T_{\mu \nu }^{(g1)}={\displaystyle \frac{m^2}{16\pi G}}\left[{\displaystyle \frac{1}{4}}g_{\mu \nu }(h^{\alpha \beta }h_{\alpha \beta }h^2)(g_{\mu \kappa }g_{\nu \lambda }h^{\kappa \lambda }hg_{\mu \kappa }g_{\nu \lambda }f^{\kappa \lambda })\right]`$ (9)
$`T_{\mu \nu }^{(g2)}={\displaystyle \frac{m^2}{16\pi G}}{\displaystyle \frac{\sqrt{f}}{\sqrt{g}}}(g_{\mu \kappa }g_{\nu \lambda }h^{\kappa \lambda }hg_{\mu \kappa }g_{\nu \lambda }f^{\kappa \lambda })`$ (10)
$`T_{\mu \nu }^{(g3)}={\displaystyle \frac{m^2}{16\pi G}}\left[{\displaystyle \frac{1}{4}}g_{\mu \nu }(H^{\alpha \beta }H_{\alpha \beta }H^2)+(H_{\mu \nu }Hf_{\mu \nu })\right]`$ (11)
$`T_{\mu \nu }^{(g4)}={\displaystyle \frac{m^2}{16\pi G}}{\displaystyle \frac{\sqrt{f}}{\sqrt{g}}}(H_{\mu \nu }Hf_{\mu \nu })`$ (12)
where we have defined $`h_{\mu \nu }=g_{\mu \nu }f_{\mu \nu }`$ and $`H^{\mu \nu }=g^{\mu \nu }f^{\mu \nu }`$. The indices of $`h_{\mu \nu }`$ are raised by $`f^{\mu \nu }`$ and the ones of $`H^{\mu \nu }`$ are lowered by $`f_{\mu \nu }`$. Note that $`H^{\mu \nu }=h^{\mu \nu }+h^{\mu \kappa }h_\kappa ^\nu +\mathrm{}`$, which explains the sign differences in the above expressions.

In this paper we wish to consider spherically symmetric stationary (SSS) solutions for these mass terms, i.e. SSS solutions of:

$$R_{\mu \nu }\frac{1}{2}Rg_{\mu \nu }=8\pi G(T_{\mu \nu }^{(g)}+T_{\mu \nu }^{(\mathrm{matt})})$$
(13)
Note that the separate diffeomorphism invariance of the matter action $`𝒮_{\mathrm{matt}}`$ (which is assumed to couple only to $`𝐠`$) implies the separate conservation (on matter shell) of the material energy-momentum tensor:

$$^\mu T_{\mu \nu }^{(\mathrm{matt})}=0$$
(14)
A consequence of Eqs. (13), (14) and of the Bianchi identities is then the separate conservation of the gravitational energy tensor:

$$^\mu T_{\mu \nu }^{(g)}=0$$
(15)
## 3 Poincaré-covariant perturbation theory

Note that all the actions (4) admit as exact solution, in the absence of matter, the “trivial vacuum” $`𝐠=𝐟`$. \[See, however, section 8 below for a discussion of non-trivial vacua.\] When representing the non-dynamical flat background metric $`𝐟`$ in Lorentzian coordinates, i.e. $`f_{\mu \nu }=\eta _{\mu \nu }`$, one can then develop a Poincaré-covariant perturbation theory: $`h_{\mu \nu }=g_{\mu \nu }\eta _{\mu \nu }=h_{\mu \nu }^{(1)}+h_{\mu \nu }^{(2)}+\mathrm{}`$. Within such an approach it is implicitly required that the massive gravitational field $`h_{\mu \nu }`$ decay at infinity, so that:

$$g_{\mu \nu }\eta _{\mu \nu },\mathrm{as}r\mathrm{}$$
(16)
In the first part of this paper, we shall impose this standard requirement, the problem being to assess the existence and continuity of solutions of massive gravity matching the trivial vacuum at infinity. \[However, we shall relax this requirement in the second part of the paper.\]

Before entering the details of our investigation of non-linear SSS solutions, let us, as a warmup, recall the basic features of perturbation theory. At the linearized approximation, the equations of motion read:

$$\mathrm{}h_{\mu \nu }+h_{\mu ,\lambda \nu }^\lambda +h_{\nu ,\lambda \mu }^\lambda \eta _{\mu \nu }h_{,\kappa \lambda }^{\kappa \lambda }h_{,\mu \nu }+\eta _{\mu \nu }\mathrm{}h+m^2(h_{\mu \nu }\alpha \eta _{\mu \nu }h)=16\pi GT_{\mu \nu }$$
(17)
Here all index raisings and lowerings are made with $`\eta _{\mu \nu }`$, and the source term $`T_{\mu \nu }`$ is equal, at the linearized approximation, to the matter energy-momentum tensor. \[See below, for the generalization of Eq. (17) to the non-trivial vacuum cases.\] Following  we have introduced a parameter $`\alpha `$ to study the special rôle of the Pauli-Fierz mass term, which corresponds to $`\alpha _{\mathrm{PF}}=1`$.

Let us recall the consequences of Eq. 17). The divergence of this equation yields, when using the conservation of the source, $`^\nu T_{\mu \nu }=0`$, the constraint:

$$h_{\mu \nu }^{,\nu }=\alpha h_{,\mu }$$
(18)
Then the trace gives:

$$2(1\alpha )\mathrm{}h+m^2(14\alpha )h=16\pi GT$$
(19)
From the latter equation we see that if $`\alpha =1`$, i.e. when we sit on the Pauli-Fierz point, the trace of $`h_{\mu \nu }`$ is locally determined by a constraint:

$$h=\frac{16\pi GT}{3m^2}$$
(20)
and thus substituting in Eq. (17) one gets:

$$\frac{1}{16\pi G}h_{\mu \nu }=\frac{1}{\mathrm{}+m^2}\left(T_{\mu \nu }\frac{1}{3}T\eta _{\mu \nu }\right)+\frac{1}{3m^2}\frac{1}{\mathrm{}+m^2}T_{,\mu \nu }$$
(21)
This result for the massive gravity field in terms of the source exhibits the difficulties of massive gravity. First, there are the phenomenological difficulties associated with the $`1/3`$ instead of the Einsteinian $`1/2`$ factor in front of $`T\eta _{\mu \nu }`$. Second, there are the theoretical difficulties associated with the presence of a factor $`m^2`$ in some terms of the solution.

Before discussing further these difficulties beyond the linearized level, let us also briefly recall the results of  concerning the non-Pauli-Fierz mass terms, i.e. the case $`\alpha 1`$. In this case, we see from Eq. (19) that the trace of $`h_{\mu \nu }`$ becomes a dynamical field, i.e. a new degree of freedom. Then the corresponding result for the expression of $`h_{\mu \nu }`$ in terms of the source is:

$`{\displaystyle \frac{1}{16\pi G}}h_{\mu \nu }={\displaystyle \frac{1}{\mathrm{}+m^2}}\left(T_{\mu \nu }{\displaystyle \frac{1}{3}}T\eta _{\mu \nu }\right){\displaystyle \frac{1}{6}}\eta _{\mu \nu }{\displaystyle \frac{1}{\mathrm{}+m_0^2}}T`$
$`+{\displaystyle \frac{2\alpha 1}{2(1\alpha )}}{\displaystyle \frac{1}{\mathrm{}+m^2}}{\displaystyle \frac{1}{\mathrm{}+m_0^2}}T_{,\mu \nu }`$ (22)
where

$$m_0^2=m^2\frac{4\alpha 1}{2(1\alpha )}$$
(23)
is the mass of the extra scalar degree of freedom. \[It is non-tachyonic if $`\alpha `$ is in the interval $`1/4\alpha <1`$\]. Note that this field is always ghostlike (independently of the value of $`\alpha `$), as seen by considering the various contributions to the action, $`h_{\mu \nu }T^{\mu \nu }`$, where the term $`(1/6)\eta _{\mu \nu }(\mathrm{}+m_0^2)^1T`$ contributes with the opposite sign from the first term (note that the last, double gradient term, does not contribute, after integration, because of the conservation of the source, $`^\nu T_{\mu \nu }=0`$). As emphasized in  the case $`\alpha 1`$ does not exhibit the difficulties of the Pauli-Fierz case. On the one hand, the expression (3) contains no dangerous denominators, vanishing with $`m^2`$, and on the other hand, it smoothly merges (thanks to the identity $`1/3+1/61/2`$) into the GR result for the integral of $`h_{\mu \nu }T^{\mu \nu }`$ when $`m^20`$. Those nice continuity properties of the case $`\alpha 1`$ pointed out in  are particularly evident in the case $`\alpha =1/2`$, but hold in all cases $`\alpha 1`$ , independently of whether the extra scalar degree of freedom is tachyonic or not. \[These continuity properties have been recently further studied in \]. Nevertheless, the presence of a new ghostlike degree of freedom makes the theory for $`\alpha 1`$ pathological as a quantum theory. In view of this quantum ghost instability, we restrict ourselves to the ghost-free Pauli-Fierz case. We note in passing that Pauli-Fierz mass terms naturally arise from higher-dimensional gravity models. Indeed, the structure of the Pauli-Fierz mass term is already encoded in the structure of the gradient terms in Einstein’s action (see, e.g., section 3.1 of  ).

Let us now briefly indicate how the Poincaré-covariant perturbation theory would proceed beyond the linear approximation. As is done in the Poincaré-covariant perturbation theory in GR one can define an effective energy tensor $`T_{\mu \nu }^{\mathrm{eff}}`$, which combines the matter energy tensor with the non-linear parts $`hh+m^2hh+\mathrm{}`$ of the left-hand side of the exact field equations (13), such that the full field equations read as the linearized ones, Eq. (17), with the replacement $`T_{\mu \nu }T_{\mu \nu }^{\mathrm{eff}}`$. This effective energy tensor is not exactly conserved but rather satisfies a relation of the type $`^\nu T_{\mu \nu }^{\mathrm{eff}}=m^2(hh+\mathrm{})_\mu `$. This modifies both Eq. (18) and Eq. (19), and thereby Eq. (21), by non-linear terms. By iteration, one can then deduce the following structure for the non-linearity expansion of the solution:

$$h\left(1+\frac{1}{mR}+\frac{1}{(mR)^2}\right)U+\left(1+\frac{1}{mR}+\mathrm{}+\frac{1}{(mR)^6}\right)U^2+\mathrm{}$$
(24)
Here, we have denoted by $`U`$ a typical gravitational potential which is, outside the source and for distances $`R\lambda _m`$, of order $`UR_S/R`$. The expansion (24) contains increasingly high inverse powers of $`m`$. We can rewrite it as a sum of terms of the type $`(R_{p,q}/R)^n`$ with a sequence of length scales defined as:

$$R_{p,q}(\lambda _m^pR_S^q)^{1/(p+q)}$$
(25)
where, as above, $`\lambda _mm^1`$.

The expansion above can numerically make sense only if the radius $`R`$ is (much) larger than all the relevant length scales $`R_{p,q}`$. In view of the extremely large difference between $`\lambda _m`$ ($`10^{28}`$ cm, say) and $`R_S`$ ($`10^5`$ cm, say, for the Sun) the various length scales tend to be extremely large, which confirms the first point made in , namely the irrelevance of perturbation theory for describing the massive gravitational field near the source. For instance, the light deflection by the Sun probes the gravitational field just outside the radius of the Sun, while the expansion above represents $`h`$ there as a sum of terms, starting with a small first term $`UR_S/R10^6`$, but continuing with extremely large “corrections” such as $`U^2/(mR)^210^{+22}`$ !

Actually, the straightforward perturbation expansion (24) can be significantly improved by gauging away some of the worse terms. Indeed, the primitive source for the “bad” inverse powers of $`m^2`$ is the last term in Eq. (21). However, this term is mainly a “gauge term”, in the sense that a suitable coordinate transformation $`x_{\mathrm{old}}^\mu x_{\mathrm{new}}^\mu `$ can remove this term from the physical “massive” metric $`g_{\mu \nu }(x)`$, at the cost of introducing it in the non-dynamical flat metric $`f_{\mu \nu }(x)`$ (now written in the “curved” coordinate system $`x_{\mathrm{new}}^\mu `$ instead of the original Lorentzian coordinate system). However, the crucial point is that this “gauging away” cannot remove all inverse powers of $`m^2`$. Indeed, if we write $`x_{\mathrm{old}}^\mu =x_{\mathrm{new}}^\mu +\xi ^\mu (x_{\mathrm{new}})`$, we have $`g_{\mu \nu }^{\mathrm{new}}(x)=_\mu (x^\alpha +\xi ^\alpha )_\nu (x^\beta +\xi ^\beta )g_{\alpha \beta }^{\mathrm{old}}(x+\xi )`$, so that the expansion of the new field $`h^{\mathrm{new}}g^{\mathrm{new}}\eta `$ contains, besides the terms $`\xi `$ which can remove the last term in Eq. (21), the new terms $`\xi h+\xi \xi `$. Therefore, with $`\xi U/m^2U/(m^2R)`$ chosen to remove the worst term $`m^2`$ in Eq. (21) , we still end up with several terms that contain inverse powers of $`m^2`$, namely:

$$h^{\mathrm{new}}U+\left(1+\frac{1}{mR}+\mathrm{}+\frac{1}{(mR)^4}\right)U^2+\mathrm{}$$
(26)
One could continue to apply suitable gauge transformations to remove some of the worse terms arising from computing higher iterations of $`h`$, but it is clear that once some negative powers of $`m^2`$ have entered $`h`$ they will keep generating higher powers in higher iterations. In view of the vast difference in scale between $`\lambda _m`$ ( $`10^{28}`$ cm) and phenomenologically tested scales $`R`$ ($`10^{11}`$ cm for the Sun radius) the dangerous inverse denominators are of order of powers of $`(mR)^210^{+34}`$ and completely outgrow the small numerators which are powers of $`U10^6`$. Finally, the improved expansion (26) is again irrelevant for describing the massive gravitational field near the source, and can only be valid in some neighbourhood of infinity. On the other hand, note that if we consider, for simplicity, the case of stationary sources the basic “Newtonian” potential $`U`$ in the expansions above will always be a Yukawa-type potential (obtained by the action of $`(\mathrm{}+m^2)^1`$ on the material source $`T`$). This potential will decrease exponentially as $`R\mathrm{}`$. Therefore it is clear that perturbation theory will be formally well defined, to all orders, when $`R\mathrm{}`$, and will generate a numerically meaningful series (say in the whole domain $`R>m^1`$). This series can be thought of as uniquely defining a specific solution of massive gravity. The main question of concern being to know whether this solution, initially defined only in a neighbourhood of infinity, can be continued into a regular solution everywhere.

Before leaving the topic of Poincaré-covariant perturbation theory, it is important to notice the following. If we apply this perturbation theory to the particular case of static spherically symmetric sources (with no motion in the Lorentzian coordinate system), it is clear that each term of perturbation theory will be static and spherically symmetric. In particular, each term will be invariant under time-reversal $`tt`$. This implies that, in Lorentzian coordinates (i.e. with $`f_{\mu \nu }=\eta _{\mu \nu }`$,) there cannot be off-diagonal terms $`g_{0i}`$ in the massive metric. Therefore, the exact solution defined (as above) by perturbation theory can be written in a bi-diagonal form, say:

$$𝐟=dt^2+dr^2+r^2d\mathrm{\Omega }^2,𝐠=a^2dt^2+b^2dr^2+c^2d\mathrm{\Omega }^2,$$
(27)
with $`a,b,c`$ some functions of $`r`$, and $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$. We shall later reconsider this bi-diagonality.

## 4 The $`\lambda `$, $`\mu `$, $`\nu `$ gauge

When discussing SSS solutions, it is instructive to work in two different gauges for the metrics. We have just introduced the $`a,b,c`$ gauge (27) (which will be used below to introduce some convenient variables $`c`$, $`\overline{c}`$, $`\overline{b}`$). It is also convenient to work in the $`\lambda `$, $`\mu `$, $`\nu `$ gauge introduced now. The latter is useful for comparing our results with the literature, while the former is similar to the gauges used in cosmological studies, and has also the merit to give a convenient Lagrangian formulation to the system. Let us write the observable line element $`ds^2=g_{\mu \nu }dx^\mu dx^\nu `$ and the reference flat metric one $`ds_{\mathrm{fl}}^2=f_{\mu \nu }dx^\mu dx^\nu `$ in the following (“Schwarzschild”) gauge:

$`ds^2=e^{\nu (R)}dt^2+e^{\lambda (R)}dR^2+R^2d\mathrm{\Omega }^2`$ (28)
$`ds_{\mathrm{fl}}^2=dt^2+\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2e^{\mu (R)}dR^2+e^{\mu (R)}R^2d\mathrm{\Omega }^2`$ (29)
where $`{}_{}{}^{}\frac{d}{dR}`$. This gauge was used in part of the literature of spherically symmetric solutions in massive gravity. It has the advantage of separating the directly observable gravitational variables $`\nu (R),\lambda (R)`$ from the “gauge” function $`\mu (R)`$. The coordinate redefinition $`r=Re^{\mu /2}`$ shows trivially that the second metric is indeed flat, and that the form (28), (29) is actually equivalent to the form (27). \[See below for the explicit link between the two sets of variables.\]

The most general matter energy-momentum tensor which respects the spherical symmetry is of the form:

$$T_{(\mathrm{matt})\nu }^\mu =\mathrm{diag}(\rho ,P_r,P_t,P_t)$$
(30)
where we have left open the possibility of having different different radial $`P_r`$ and tangential $`P_t`$ pressures. In the case of a “fluid” source, we shall have $`P_r=P_t`$. Then the $`(t,t)`$ and $`(R,R)`$ components of the Einstein equations read:

$`e^{\nu \lambda }\left[{\displaystyle \frac{\lambda ^{}}{R}}+{\displaystyle \frac{1}{R^2}}(e^\lambda 1)\right]=m^2f_t(\lambda ,\mu ,\nu ,\mu ^{},R)+8\pi G\rho e^\nu `$ (31)
$`{\displaystyle \frac{\nu ^{}}{R}}+{\displaystyle \frac{1}{R^2}}(1e^\lambda )=m^2f_R(\lambda ,\mu ,\nu ,\mu ^{},R)+8\pi GP_re^\lambda `$ (32)
where the quantities $`f_t`$ and $`f_R`$ are proportional to the $`(t,t)`$ and $`(R,R)`$ components of the energy momentum tensor $`T_{\mu \nu }^{(g)}`$ generated by the potential term, and are defined as:

$$f_t\frac{1}{4}\frac{1}{\sqrt{g}}\frac{\delta 𝒱}{\delta g^{tt}}\frac{8\pi G}{m^2}T_{tt}^{(g)}\mathrm{and}f_R\frac{1}{4}\frac{1}{\sqrt{g}}\frac{\delta 𝒱}{\delta g^{RR}}\frac{8\pi G}{m^2}T_{RR}^{(g)}$$
(33)
From the conservation of the energy-momentum tensor generated by the mass term, i.e. from $`^\mu T_{\mu \kappa }^{(g)}`$=0, we obtain a third equation from the $`\kappa =r`$ component (the $`\kappa =t,\theta ,\varphi `$ components are identically zero). Defining:

$$f_g\frac{8\pi G}{m^2}^\mu T_{\mu r}^{(g)}$$
(34)
we have the following constraint equation for $`m0`$:

$$f_g(\lambda ,\mu ,\nu ,\lambda ^{},\mu ^{},\nu ^{},\mu ^{\prime \prime },R)=0$$
(35)
Note that this equation has no $`m`$ dependence. In the Appendix A we present the functions $`f_t`$, $`f_R`$ and $`f_g`$ for all the potentials that we are considering.

Finally, from the conservation of the matter energy-momentum tensor we have the final equation which closes the system:

$$2P_r^{}+\frac{1}{R}(3P_r\rho 4P_t)+\left(\nu ^{}+\frac{1}{R}\right)(\rho +P_r)=0$$
(36)
Note that in the “fluid” case, i.e. for an isotropic pressure, the latter equation simplifies to:

$$P^{}=\frac{\nu ^{}}{2}(\rho +P)$$
(37)
### 4.1 The General Relativity $`m=0`$ limit

In the GR limit, one must discard Eq. (35), and set $`m^2`$ to zero in Einstein’s equations (31), (32). It will be helpful to remind ourselves of the solution of the Einstein equations in GR for a star of constant density $`\rho =\rho _0`$ and isotropic pressure $`P=P_r=P_t`$. In this case (considered here for simplicity), the exterior (vacuum) solution , after absorbing a constant by a time rescaling, reads:

$$e^\nu =1\frac{R_S}{R},\lambda =\nu $$
(38)
where $`R_S2GM`$ is the Schwarzschild radius and appears here as an integration constant.

On the other hand, inside the star, eq. (31) can be integrated to:

$$e^\lambda =1\frac{2Gm(R)}{R}$$
(39)
where we have defined $`m(R)=\frac{4}{3}\pi \rho _0R^3`$. Note that we have imposed the constraint $`\lambda (R=0)=0`$, which is necessary to avoid a conical singularity at $`R0`$. If we denote the radius of the star by $`R_{}`$, matching to the external solution gives the links $`m(R_{})=M`$ and $`R_S=\frac{8}{3}\pi G\rho _0R_{}^3`$. We then can solve the remaining equations and obtain the well known result for the pressure:

$$P=\rho _0\frac{\sqrt{1\frac{R_S}{R_{}^3}R^2}\sqrt{1\frac{R_S}{R_{}}}}{3\sqrt{1\frac{R_S}{R_{}}}\sqrt{1\frac{R_S}{R_{}^3}R^2}}$$
(40)
and the redshift function $`\nu `$ (corresponding to the same rescaled time as the external solution above):

$$e^\nu =\left(\frac{3}{2}\sqrt{1\frac{R_S}{R_{}}}\frac{1}{2}\sqrt{1\frac{R_S}{R_{}^3}R^2}\right)^2$$
(41)
### 4.2 The perturbative limit for $`m0`$

Let us now see how the results change if we consider the massive theory. We will firstly deal with the linearized equations of motion with the aim of understanding analytically how one can match the unique decaying exterior solution into a regular interior one. In other words we will assume that $`\lambda ,\mu ,\nu 1`$. We will solve these linearized equations in the interior and the exterior of the star and then match them to determine completely the various integration constants.

#### 4.2.1 The exterior star solution

The equations (31), (32), (35) at linear order in $`\lambda ,\mu ,\nu `$ are the following (for all four mass terms):

$`{\displaystyle \frac{\lambda ^{}}{R}}+{\displaystyle \frac{\lambda }{R^2}}={\displaystyle \frac{m^2}{2}}(\lambda +3\mu +R\mu ^{})`$ (42)
$`{\displaystyle \frac{\nu ^{}}{R}}{\displaystyle \frac{\lambda }{R^2}}={\displaystyle \frac{m^2}{2}}(\nu +2\mu )`$ (43)
$`\nu ^{}={\displaystyle \frac{2\lambda }{R}}`$ (44)
Note that the linearization of the equations has removed entirely the function $`\mu `$ from the constraint equation (35). These equations, give a simple second-order differential equation for $`\nu `$:

$$\nu ^{\prime \prime }+\frac{2}{R}\nu ^{}m^2\nu =0$$
(45)
Then the other variables are obtained from $`\nu `$ via:

$$\lambda =\frac{R}{2}\nu ^{},\mu =\frac{\nu ^{}}{2m^2R}\frac{\nu }{2}$$
(46)
The generic solution of the $`\nu `$ equation above contains two integration constants, $`C_1`$ and $`C_1^{}`$. The solution proportional to $`C_1`$ is exponentially decaying at infinity, while the one proportional to $`C_1^{}`$ is exponentially growing. In keeping with our boundary condition (16) we shall reject the exponentially growing solution and consider only the other one, namely:

$`\nu ={\displaystyle \frac{C_1}{R}}e^{mR}`$ (47)
$`\lambda ={\displaystyle \frac{mC_1}{2}}\left(1+{\displaystyle \frac{1}{mR}}\right)e^{mR}`$ (48)
$`\mu ={\displaystyle \frac{C_1}{2R}}\left(1+{\displaystyle \frac{1}{mR}}+{\displaystyle \frac{1}{(mR)^2}}\right)e^{mR}`$ (49)
Note that, if needed, the exponentially growing solution is simply obtained from the latter solution by changing $`mm`$. The integration constant $`C_1`$ is proportional to the Schwarzschild radius $`R_S2GM`$ of the source (the precise link will be seen later when we match with the interior). In the region $`Rm^1`$ we obtain the limits of the above solution for $`\lambda ,\mu ,\nu `$:

$$\nu \frac{C_1}{R},\lambda \frac{\nu }{2}\frac{C_1}{2R},\mu \frac{C_1}{2m^2R^3}$$
(50)
Note that the second relation differs from the corresponding GR relation $`\lambda \nu `$. These show the well known fact that there is a finite discontinuity between this metric and the one obtained in GR, independent of the mass of the graviton as was shown in . Additionally, the last relation shows that, in the $`\lambda ,\mu ,\nu `$ gauge, and at the linearized approximation, there is a dangerous denominator $`m^2`$ only in the “gauge” function $`\mu `$. Actually, one can see that the gauge function $`\mu `$ is similar to the gradient $`\xi `$ of the coordinate transformation that we introduced above to remove the leading $`𝒪(m^2)`$ terms in the (linearized) perturbative field $`h`$. We had above $`\xi U/(mR)^2`$ which indeed corresponds to $`\mu R_S/(m^2R^3)`$. Therefore, as above, the choice of gauge has allowed us to “regularize” the linearized physical metric $`\nu ,\lambda `$, but we expect that non-linear effects will re-introduce some negative powers of $`m^2`$ in $`\nu ,\lambda `$.

#### 4.2.2 The interior star solution

In order to be able to apply our perturbative calculation to the problem of source matching, we need to formally consider a star with sufficiently small density, and sufficiently large radius, so that our assumptions $`\lambda ,\mu ,\nu 1`$ remain satisfied everywhere. We shall give later (after having taken into account non-linear effects) the precise conditions that the star characteristics must satisfy. Evidently, such a low density “star” is not physically relevant to the discussion of real stars, but the point of this subsection is to show in detail how, in principle, one can match the well-defined exterior perturbative solution to a unique interior one.

The interior linearized equations of motion (again for all four mass terms) are:

$`{\displaystyle \frac{\lambda ^{}}{R}}+{\displaystyle \frac{\lambda }{R^2}}={\displaystyle \frac{m^2}{2}}(\lambda +3\mu +R\mu ^{})+8\pi G\rho `$ (51)
$`{\displaystyle \frac{\nu ^{}}{R}}{\displaystyle \frac{\lambda }{R^2}}={\displaystyle \frac{m^2}{2}}(\nu +2\mu )+8\pi GP`$ (52)
$`\nu ^{}={\displaystyle \frac{2\lambda }{R}}`$ (53)
and additionally the continuity equation is:

$$P^{}=\frac{\nu ^{}}{2}(\rho +P)$$
(54)
These equations, give a differential equation for $`\nu `$:

$$\nu ^{\prime \prime }+\left(\frac{2}{R}+\frac{8}{3}\pi G(\rho +P)R\right)\nu ^{}m^2\nu =8\pi G\left(2P+\frac{4}{3}\rho \right)$$
(55)
Let us assume, for the sake of simplicity, that $`\rho =const.\rho _0`$ and $`P\rho _0`$, i.e. a constant density star of low pressure, so that we can neglect the pressure $`P`$ in the above equation. This assumption means that we are considering a non-compact start, with $`R_{}R_S`$, i.e. $`R_{}(8\pi G\rho _0)^{1/2}`$. We can then safely neglect the second addendum in the $`\nu ^{}`$ parentheses. Thus, we have the simple differential equation:

$$\nu ^{\prime \prime }+\frac{2}{R}\nu ^{}m^2\nu =8\pi G\frac{4}{3}\rho _0$$
(56)
which differs from the exterior equation (45) above by the source term on the right-hand side (RHS). Under the same approximations, the other metric functions are determined in terms of a solution of the interior $`\nu `$ equation by means of:

$$\lambda =\frac{1}{2}R\nu ^{},\mu =\frac{\nu ^{}}{2m^2R}\frac{\nu }{2}$$
(57)
The general solution for the $`\nu `$ function that is regular at $`R=0`$ depends on only one integration constant, say $`C_2`$ , and reads:

$$\nu (R)=\frac{4}{3}\frac{8\pi G\rho _0}{m^2}+\frac{C_2}{R}\mathrm{sinh}(mR)$$
(58)
The continuity of $`\nu `$ relates the two integration constants, $`C_1`$ of the exterior solution and $`C_2`$ of the interior one, as:

$$C_2=\frac{R_{}}{\mathrm{sinh}(mR_{})}\left[\frac{4}{3}\frac{8\pi G\rho _0}{m^2}\frac{C_1}{R_{}}e^{mR_{}}\right]$$
(59)
The continuity of the first derivative of $`\nu `$, determines the constant $`C_1`$ in terms of the mass and the radius of the star:

$$C_1=\frac{4R_S}{(mR_{})^3}[mR_{}\mathrm{cosh}(mR_{})\mathrm{sinh}(mR_{})]$$
(60)
Consistently with our problem of studying the $`m0`$ limit, we shall assume that the radius of the star is such that $`R_{}m^1`$. We then get that:

$$C_1\frac{4}{3}R_S$$
(61)
\[In the opposite limit, $`m^1R_{}`$ we would have obtained that $`C_1\frac{4R_S}{(mR_{})^2}e^{mR_{}}`$\]. Note that the factor $`\frac{4}{3}`$ here is linked to the well-known fact that, if one assumes the validity of linearized theory, one must renormalize the “bare Newton constant” $`G`$ appearing in the massive gravity action by $`\frac{4}{3}GG_N`$ to recover the usual Newtonian constant $`G_N`$.

Finally, the full solution is obtained by inserting the matched solution for $`\nu `$ into the expressions (57) above for the functions $`\lambda `$ and $`\mu `$.

Let us finally sketch the structure of the (unique) solution of perturbation theory, in the $`\lambda ,\mu ,\nu `$ gauge, when taking into account the next order in perturbation theory. If we introduce $`Ue^{mR}R_S/RR_S/R`$ (in the region $`R_{}<Rm^1`$), we can write:

$`\nu \lambda U+\left(1+{\displaystyle \frac{1}{mR}}+\mathrm{}+{\displaystyle \frac{1}{(mR)^4}}\right)U^2+\mathrm{},`$ (62)
$`\mu {\displaystyle \frac{U}{(mR)^2}}+\left({\displaystyle \frac{1}{(mR)^2}}+\mathrm{}+{\displaystyle \frac{1}{(mR)^6}}\right)U^2+\mathrm{},`$ (63)
This expansion is related to the general perturbation expansions (24), (26) above. In fact, roughly speaking the expansion for the “physical variables” $`\nu ,\lambda `$ corresponds to the gauge-improved expansion (26), while the expansion for the gauge variable $`\mu `$ corresponds to the expansion of the gradient $`\xi `$ of the “improving” gauge transformation $`\xi `$ introduced in section 2 above.

One would need to study more carefully the structure of higher order terms in the expansions above to delineate what are the most relevant length scales $`R_{p,q}`$ , see Eq. (25), determining the range of validity of the $`\lambda ,\mu ,\nu `$ perturbation theory, i.e. for determining the constraints on $`R_{}`$ and $`\rho _0`$ ensuring that the series above (considered both in the interior and the exterior) make sense. As it is clear that the crucial powers of $`mR`$ entering the successive monomials $`U^n/(mR)^p`$ will increase linearly with the perturbation order $`n`$, there will be a finite limit, when $`n\mathrm{}`$, to the sequence of relevant scales $`R_{p,q}`$, and therefore this sequence will have a finite, global least upper bound, or supremum, $`R_{\mathrm{sup}}`$. Perturbation theory is then valid (everywhere) if $`R_{\mathrm{sup}}R_{}`$.

It is evident that since we need (for being able to use perturbation theory) to constrain ourselves to stars of extremely small density, the solution that we have obtained is unrealistic. However, we have gone through this explicit derivation to show how, in principle, our various physical requirements (decay at infinity, matching at the star radius, and regularity at the origin) determine a unique solution in terms of the equation of state ($`\rho _0`$) and the mass of the star, and the mass of the graviton. Note that this uniqueness of the solution is a non-trivial consequence of the treatment based on perturbation theory. Indeed, perturbation theory allowed not only to provide boundary conditions at infinity that selected one solution from the second-order differential equation for $`\nu `$, but it also allowed one to kill the “sixth degree of freedom” linked (in the $`\lambda ,\mu ,\nu `$ gauge) to the fact that the equation $`f_g=0`$ is a second-order differential equation in $`\mu `$. The perturbative treatment allowed us to recursively determine the higher derivatives of $`\mu `$ from the algebraic calculation of $`\mu `$ in terms of $`\nu `$ and $`\nu ^{}`$ (see Eq. (46)).

## 5 On a claim by Vainshtein

In the old paper of Vainshtein , it was claimed that continuity in the massless limit could be obtained if one, instead of treating the problem perturbatively, used a different “non-perturbative” expansion, proceeding in positive powers of $`m^2`$, whose first terms for the physical variables $`\lambda `$ and $`\nu `$ were their usual GR values, and whose first term for the gauge variable $`\mu `$ were (in the exterior region) of order:

$$\mu \sqrt{\nu }\sqrt{\lambda }\sqrt{R_S/R}$$
(64)
The aim of this section is to critically reexamine this claim, to explain in detail some points which have remained unclear, and to discuss its serious difficulties.

A first point to clarify concerns the theoretical framework chosen to “define” massive gravity. Indeed, the starting point of the original reference  was not an action, but the Einstein equation with a postulated energy-momentum tensor for a massive graviton. However, it is easy to see (by checking the dissymmetry in $`(x,y)`$ of the functional derivative of the postulated $`\sqrt{g(x)}T_{\mu \nu }^{(V)}(x)`$ with respect to $`g^{\alpha \beta }(y)`$) that the field equations postulated in  cannot be derived from any action (and, in particular, they correspond to none of the four models written down at the beginning of this paper). This makes them quantum mechanically inconsistent. The energy-momentum tensor postulated in  reads (after correcting the sign of the $`m^2`$ terms as they appear in Eq. (1) there, which actually corresponded to a tachyonic mass term):

$$T_{\mu \nu }^{(V)}=\frac{m^2}{16\pi G}(h_{\mu \nu }h\eta _{\mu \nu })$$
(65)
Note also that, in the same paper, the explicit (t,t) and (R,R) components of the Einstein equations appear with the correct sign, consistent with the $`T_{\mu \nu }^{(V)}`$ given above. For this energy momentum tensor the functions $`f_t`$, $`f_R`$ and $`f_g`$ are given in Appendix A. Note also that in , $`f_g`$ was written down incorrectly, nevertheless the final result of  starting from Eq. (65) is mathematically correct.

The basic idea of  was the following. The root of the discontinuous behaviour between massive gravity and GR is the “new” constraint $`f_g`$, whose non-linearity expansion starts with (when keeping terms relevant when Eq. (64) holds):

$$\frac{1}{R}f_g=\frac{\lambda }{R^2}\frac{\nu ^{}}{2R}+\frac{4\mu \mu ^{}}{R}+\frac{7\mu ^2}{4}+\mu \mu ^{\prime \prime }+\mathrm{}=0$$
(66)
The ellipsis in (66) contain quadratic terms $`𝒪(\mu \nu +\mu \lambda +\nu ^2+\mathrm{})`$ as well as higher non-linearities.

When this equation is solved in the usual perturbation theory, i.e. by solving for the linear terms $`𝒪(\lambda )`$ and $`𝒪(\nu ^{})`$ and then adding the quadratic terms $`𝒪(\mu ^2+\mu \nu +\mathrm{})`$ as corrections, one gets the “discontinuous” result (50) above. The observation then was that one might cure the discontinuity by assuming that the $`𝒪(\mu ^2)`$ terms in Eq. (66) are comparable to the linear ones $`𝒪(\lambda )`$ and $`𝒪(\nu ^{})`$. In more detail, one assumes that the three functions $`\lambda `$, $`\mu `$, $`\nu `$ admit an expansion in positive powers of $`m^2`$, say:

$$f(R)=\underset{n=0}{\overset{\mathrm{}}{}}m^{2n}f_n(R)$$
(67)
where the leading terms, $`\nu _0`$, $`\lambda _0`$, $`\mu _0`$, are such that the first two are solutions of Eqs. (31), (32) with $`m^2`$ set to zero, i.e. (in the exterior):

$$\nu _0(R)=\lambda _0(R)=\mathrm{ln}\left(1\frac{R_S}{R}\right)=\frac{R_S}{R}+\frac{1}{2}\left(\frac{R_S}{R}\right)^2+\mathrm{}$$
(68)
while $`\mu _0`$ is a solution of the constraint equation (66) with the above given $`\nu _0`$, $`\lambda _0`$ substituted, i.e.

$$\frac{\lambda _0}{R^2}\frac{\nu _0^{}}{2R}+\frac{4\mu _0\mu _0^{}}{R}+\frac{7\mu _0^2}{4}+\mu _0\mu _0^{\prime \prime }+\mathrm{}=0$$
(69)
One is looking for a solution $`\mu _0`$ of Eq. (69) which satisfies (in the exterior region) Eq. (64) at lowest order. Note that the integration constant associated with $`\nu _0`$ has been set to zero so that one is consistent with the “trivial-vacuum” (see section 3). Note also that the quantity $`R_S`$ entering Eq. (68) is, at this stage, a GR-like integration constant whose exact link with the mass of the source is not important.

One then finds that the above requirements lead to a unique (exterior) solution for $`\mu _0`$, which further admits an expansion in powers of $`R_S/R`$ starting as:

$$\mu _0=\sqrt{\frac{8}{13}}\sqrt{\frac{R_S}{R}}\left(1+\frac{61}{12\sqrt{26}}\sqrt{\frac{R_S}{R}}+\frac{39281}{13104}\frac{R_S}{R}+\mathrm{}\right)$$
(70)
where the higher order terms are obtained by considering the constraint equation (69) augmented with higher-order terms in $`\lambda _0`$, $`\nu _0`$ and $`\mu _0`$.

Starting from this result for $`\mu _0`$ (where all higher-order coefficients are, in principle, uniquely determined), we then go back to Eq. (31) (in the exterior region), which yields a linear first-order differential equation for the $`𝒪(m^2)`$ term $`\lambda _1`$ in $`\lambda `$:

$$\left(\frac{1}{R^2}+\mathrm{}\right)\lambda _1+\left(\frac{1}{R}+\mathrm{}\right)\lambda _1^{}=\frac{5}{\sqrt{26}}\sqrt{\frac{R_S}{R}}\frac{257}{156}\frac{R_S}{R}+\mathrm{}$$
(71)
The general solution of this inhomogeneous linear ordinary differential equation (ODE) can be written as the sum of a particular inhomogeneous solution:

$$m^2\lambda _1=(mR)^2\left(\sqrt{\frac{2}{13}}\sqrt{\frac{R_S}{R}}\frac{257}{312}\frac{R_S}{R}+\mathrm{}\right)$$
(72)
plus the general solution of the homogeneous equation, which is of the form $`c_{\lambda 1}/R`$. The latter general homogeneous solution can be absorbed in a $`𝒪(m^2)`$ change of the integration constant $`R_S`$, and can therefore be physically ignored.

We can then insert the known values of $`\nu _0`$, $`\lambda _0`$, $`\mu _0`$ and $`\lambda _1`$ in Eq. (32) to get a linear first-order differential equation for the $`𝒪(m^2)`$ term $`\nu _1`$ in $`\nu `$:

$$\left(\frac{1}{R}+\mathrm{}\right)\nu _1^{}=\sqrt{\frac{2}{13}}\sqrt{\frac{R_S}{R}}\frac{13}{24}\frac{R_S}{R}+\mathrm{}$$
(73)
Again, the general solution of this inhomogeneous linear ODE can be written as the sum of a particular inhomogeneous solution:

$$m^2\nu _1=(mR)^2\left(\frac{2}{3}\sqrt{\frac{2}{13}}\sqrt{\frac{R_S}{R}}\frac{13}{24}\frac{R_S}{R}+\mathrm{}\right)$$
(74)
plus the general solution of the homogeneous equation, which is simply a constant $`c_{\nu 1}`$. The latter constant can be absorbed in a rescaling of the time variable, and can therefore be physically ignored.

The next step is more tricky. Indeed, we get a second-order linear ODE for the $`𝒪(m^2)`$ term $`\mu _1`$ in $`m`$ of the type:

$$(\mu _0+\mathrm{})\mu _1^{\prime \prime }+\left(2a\mu _0^{}+\frac{b}{R}\mu _0+\mathrm{}\right)\mu _1^{}+\left(\mu _0^{\prime \prime }+\frac{b}{R}\mu _0^{}+\mathrm{}\right)\mu _1=\frac{3}{\sqrt{26}}\sqrt{\frac{R_S}{R}}+\frac{813}{624}\frac{R_S}{R}+\mathrm{}$$
(75)
Here $`a=7/4`$ and $`b=4`$ are the coefficients appearing in the $`\mu `$-quadratic terms in the constraint (66). The problem now is that the general homogeneous solution of Eq. (75) a priori introduces two new integration constants $`c_{m1+}`$, $`c_{m1}`$ in a solution of the type $`\mu _1^{\mathrm{hom}}=c_{\mu 1+}R^{s_+}+c_{\mu 1}R^s_{}`$ where the exponents $`s_+`$, $`s_{}`$ are the two roots of the quadratic indicial equation $`s^2+(ba1)s+3/4b/2=0s_\pm =\frac{\sqrt{5}}{8}(\sqrt{5}\pm \sqrt{21})`$. Contrary to what happened above for the integration constants entering the homogeneous solutions in $`\lambda _1`$ and $`\nu _1`$, the $`\mu _1`$-integration constants cannot a priori be physically ignored. However, if we follow the spirit of  and of the leading requirement (64), i.e. if we try to construct the simplest vacuum solution which is entirely determined by the GR-like integration constant $`R_S`$, we can discard the above homogeneous solution and continue the iteration by selecting for $`\mu _1`$ the particular solution determined by the inhomogeneous “source terms” in the $`\mu _1`$-equation. This solution has the following solution:

$$m^2\mu _1=(mR)^2\left(\frac{1}{7}+\frac{1305}{644\sqrt{26}}\sqrt{\frac{R_S}{R}}+\mathrm{}\right)$$
(76)
Note that the exponents $`s_+`$, $`s_{}`$ entering the $`\mu _1`$-homogeneous solution are irrational. Taking into account the possibility of adding this homogeneous solution would introduce a new sequence of characteristic length scales (25) which would mix in a complicated way with the other (rational) scales.

Continuing in the same way (always discarding the homogeneous solutions entering the $`\mu `$ equations) we end up with a double series for the three functions $`f=\nu ,\lambda ,\mu `$ of the form:

$$f=\left(\frac{R_S}{R}\right)^a\underset{n=0}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{\mathrm{}}{}}f_{nk}\left(\frac{m^2R^{5/2}}{R_{S}^{}{}_{}{}^{1/2}}\right)^n\left(\frac{R_{S}^{}{}_{}{}^{1/2}}{R^{1/2}}\right)^k$$
(77)
where $`a=1`$ for $`\lambda `$, $`\nu `$ and $`a=1/2`$ for $`\mu `$. For notational simplicity, we have not indicated in Eq. (77) the presence of the logarithmic terms $`\mathrm{log}(R/S)`$ which appear each time the “source terms” for the perturbations in $`\lambda `$ and $`\nu `$ contain a power of $`R`$ which match (modulo a factor $`R^2`$) the power of $`R`$ of a homogeneous solution of the left-hand-side. The arbitrary scale $`S`$ that one can introduce in these logs is not physically important because a rescaling of $`S`$ corresponds to the addition of an (ignorable) homogeneous solution in $`\lambda `$ and $`\nu `$. \[Because of the irrational character of $`s_+`$, $`s_{}`$, no logs primitively enter from the perturbation equations for the $`\mu `$ variable.\] Without loss of generality, we can then choose $`S=R_S`$ as scale in all the logs.

The final result (77) thereby involves only two length scales $`R_S`$ and the length scale entering the monomial $`\frac{m^2R^{5/2}}{R_{S}^{}{}_{}{}^{1/2}}`$. For this expansion to make sense, we should make sure that both $`\frac{R_S}{R}`$ and $`\frac{m^2R^{5/2}}{R_{S}^{}{}_{}{}^{1/2}}`$ be (much) smaller than 1. This requirement is satisfied if we look at distances $`R_SRR_V`$, where $`R_V`$ is the length scale introduced in Eq. (3), which, using the general definition (25) can be written as:

$$R_V=R_{4,1}$$
(78)
To summarize so far: we have verified part of the claim of , namely the existence of a (particular) double expansion, a priori valid in the intermediate range $`R_SRR_V`$, that (i) represents a formal (vacuum) local solution of the Einstein equations modified by the mass term (65), (ii) is continuous as $`m^20`$, and (iii) satisfies Eq. (64).

The weak points of this argument are several. First, the specific local solution chosen in Eq. (77) corresponds to a postulated energy momentum tensor which does not stem from an action, and we are going to see that a similar solution does not exist for the action used as a starting point in . But the weakest point concerning the expansion (77) is the lack of rationale indicating that this local solution does approximate the unique global one satisfying the standard boundary conditions at infinity (16) and being regular at the origin. The detailed argument above, showing that two arbitrary integration constants have been set to zero to obtain Eq. (77), does not make it a priori probable that Eq. (77) happens to satisfy all the required boundary conditions.

Before tackling the more complicated “global” issues, let us first study the issue of the existence of formal Vainshtein-type expansions for action-based field equations. We start by considering the potential $`𝒱^{(1)}`$ that appeared in the recent paper of , which re-discussed the result of . Trying to repeat the same procedure as above we face the following problem. The equation for $`\mu _0`$ now reads:

$$\frac{\lambda _0}{R^2}\frac{\nu _0^{}}{2R}\frac{6\mu _0\mu _0^{}}{R}\frac{3\mu _0^2}{4}\frac{3\mu _0\mu _0^{\prime \prime }}{2}+\mathrm{}=0$$
(79)
The $`\mu `$-quadratic terms in this equation are different from those in Eq. (69). This difference is crucial because if we look again for a solution satisfying Eq. (64), i.e. admitting an expansion of the form $`\mu _0=\mu _{00}\sqrt{R_S/R}+\mu _{01}(R_S/R)+\mathrm{}`$ we find that there exists no such solution888This contradicts the statements of , which were supposed to concern the mass term $`𝒱^{(1)}`$. They, instead, seem to correspond to the postulated $`T_{\mu \nu }^{(V)}`$ of .. \[Formally, the only such solution starts with an imaginary coefficient: $`\mu _{00}=\pm i\sqrt{\frac{8}{27}}`$.\]

The same problem persists if we consider the potential $`𝒱^{(2)}`$. Then the same equation reads:

$$\frac{\lambda _0}{R^2}\frac{\nu _0^{}}{2R}\frac{2\mu _0\mu _0^{}}{R}\frac{\mu _0^2}{4}\frac{\mu _0\mu _0^{\prime \prime }}{2}+\mathrm{}=0$$
(80)
and would again formally give an imaginary solution: $`\mu _{00}=\pm i\sqrt{\frac{8}{9}}`$.

The situation is better, however, for the potentials $`𝒱^{(3)}`$ and $`𝒱^{(4)}`$. The approximated equation for $`𝒱^{(3)}`$ is then:

$$\frac{\lambda _0}{R^2}\frac{\nu _0^{}}{2R}+\frac{2\mu _0\mu _0^{}}{R}+\frac{\mu _0^2}{4}+\frac{\mu _0\mu _0^{\prime \prime }}{2}+\mathrm{}=0$$
(81)
which gives $`\mu _{00}=\sqrt{\frac{8}{9}}`$ and for $`𝒱^{(4)}`$ we have:

$$\frac{\lambda _0}{R^2}\frac{\nu _0^{}}{2R}+\frac{6\mu _0\mu _0^{}}{R}+\frac{3\mu _0^2}{4}+\frac{3\mu _0\mu _0^{\prime \prime }}{2}+\mathrm{}=0$$
(82)
which gives $`\mu _{00}=\sqrt{\frac{8}{27}}`$. Thus, at least in these cases, a Vainshtein-type expansion similar to (77) exists. Again, for these expansions to make sense one should be at the intermediate range of distances $`R_SRR_V`$, with $`R_V`$ given by Eq. (3).

Let us note in passing that the new scale selected by a Vainshtein-type construction is not at all universal for different mass terms. For example, let us consider the potential :

$$𝒱^{(\sigma )}=(fg)^{1/4}(\sigma _2\sigma _1^2)$$
(83)
where $`\sigma _n\mathrm{\Sigma }_a(\mathrm{ln}\lambda _a)^n`$, with $`\lambda _a`$ denoting the eigenvalues of $`𝐟^1𝐠`$ (see below). Then, the non-linearity expansion of the constraint (35) gives:

$$\frac{1}{R}f_g=\frac{\lambda }{R^2}\frac{\nu ^{}}{2R}+\frac{3\mu ^2\mu ^{}}{8R}+\frac{\mu \mu ^2}{4}+\frac{R\mu ^3}{48}+\frac{3\mu ^2\mu ^{\prime \prime }}{32}+\frac{R\mu \mu ^{}\mu ^{\prime \prime }}{16}+\mathrm{}=0$$
(84)
If we introduce the short-hand notations $`ϵmR_S`$ and $`xR/R_S`$, we end up with instead of Eq. (77), expansions of the type:

$$f=x^a\underset{n=0}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{\mathrm{}}{}}f_{nk}\left(ϵ^2x^{8/3}\right)^n\left(x^{1/3}\right)^k$$
(85)
with $`a=1`$ for $`\lambda `$, $`\nu `$ and $`a=1/3`$ for $`\mu `$. This expansion has the same form even if we add the stabilizing $`\lambda \sigma _2^2`$ term (see ). By the same argument as before, this gives the following range of distances for which the above expansion makes sense:

$$R_SR\frac{R_S}{(mR_S)^{3/4}}R_{3,1}$$
(86)
A final example would be to consider the following potential:

$$𝒱^{(0)}=(fg)^{1/4}\sigma _2^2$$
(87)
for which the mass of the graviton is zero, but still we have a non-trivial potential which breaks general covariance. Then the expansion at lowest order looks like:

$$f=x^1\left[f_{00}+f_{10}ϵ^2\mathrm{log}x+\mathrm{}\right]$$
(88)
for all three functions $`\mu `$, $`\nu `$ and $`\lambda `$. The formal range of validity of this expansion is now:

$$R_SRR_Se^{1/(mR_S)^2}$$
(89)
At this point, it is important to remark that, contrary to what one might think (in view of the controversy over the “discontinuity” issue) there exist many local solutions  of the massive-gravity field equations which are continuous as $`m^20`$ and which, therefore, might suggest that massive gravity is compatible with local tests of GR. The specific Vainshtein construction is (when it works) one way of exhibiting such solutions (constrained to the obtention of strictly stationary solutions). However, we have seen above that this specific way fails in important cases, such as the simplest mass term $`𝒱^{(1)}`$. However, even when there does not exist a $`\mu _0`$ satisfying the simple requirement (64), there might still exist other solutions. Indeed, the equation that $`\mu _0`$ has to satisfy is a second-order quasi-linear ODE and one can always construct local solutions (depending on two arbitrary parameters) for such an equation. \[The problem then is to show that the domain of existence of the solution is large enough to cover the phenomenological tests.\] In fact, if we forget for a moment about the restriction to strictly stationary solutions, Ref.  has indicated how to construct general classes of solutions, evolving on the time scale $`m^1`$ (thought of as the Hubble time scale), which are $`𝒪(m^2)`$ close to any desired GR solution everywhere in space. All these constructions leave open, however, two basic issues: (i) can these local solutions be extended to globally asymptotically flat solutions, or, if they cannot, (ii) can they be extended to a cosmological solution which is a natural “attractor” of the cosmological dynamics, so that it is indeed natural to use such a solution to describe the universe around us?

Before tackling the issue of whether the formal, particular Vainshtein-type expansions are part of a globally regular solution (both at the origin and at infinity), we study in the next section the structure of the field equations in a different gauge.

## 6 The $`c`$, $`\overline{c}`$, $`\overline{b}`$ variables

Let us consider the $`a,b,c`$ gauge (27), and replace the three basic variables $`a,b,c`$ by the equivalent combinations $`c`$, $`\overline{c}`$, $`\overline{b}`$ where we define:

$$\overline{b}=ab,\overline{c}=ca^2$$
(90)
Let us indicate the link between the $`\lambda `$, $`\mu `$, $`\nu `$ variables and the $`a,b,c`$ ones (or the equivalent $`c`$, $`\overline{c}`$, $`\overline{b}`$ ones). We first compute from the three functions $`a(r)`$, $`b(r)`$, $`c(r)`$ the functions $`\lambda (r)`$, $`\mu (r)`$, $`\nu (r)`$ as follows:

$$a^2(r)=e^{\nu (r)},b^2(r)=e^{\lambda (r)+\mu (r)}\left(1+\frac{r}{2}\frac{d\mu }{dr}\right)^2,c^2(r)=r^2e^{\mu (r)}$$
(91)
Then we make the coordinate transformation $`rR`$ by inverting $`R=c(r)=re^{\mu (r)/2}`$.

The advantage of introducing these variables is that they simplify the Lagrangian formulation of the dynamics. Indeed, the action can easily be checked to be:

$$𝒮=\frac{1}{2G}𝑑t𝑑r\left(\frac{\dot{c}\dot{\overline{c}}}{\overline{b}}+\overline{b}\right)\frac{m^2}{16G}𝑑t𝑑r\overline{𝒱}(c,\overline{c},\overline{b},r)+𝒮_{\mathrm{matt}}$$
(92)
where we have introduced $`\overline{𝒱}𝒱/\mathrm{sin}\theta `$. Note that $`𝒱`$ has the scaling behaviour $`𝒱=r^2f(c/r,\overline{c}/r,\overline{b})`$. We have denoted $`\dot{}d/dr`$. The pair $`(c,\overline{c})`$ appears here as “light-cone” coordinates for a “relativistic particle”. We refer to the cosmological studies of massive gravity for a discussion of similar “relativistic particle Lagrangians”. Note that the variable $`\overline{b}`$ is a “gauge-like” variable (radial lapse) which does not have a kinetic term, just as the (relative) time lapse variable $`e^\gamma `$ did in cosmological studies. Therefore the equation of motion for $`\overline{b}`$ is an algebraic equation, instead of the second-order ODE’s one gets for the dynamical variables $`(c,\overline{c})`$.

The equations of motion for this gauge choice are:

$`{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{dc}{dr}}\right)={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{\overline{𝒱}}{\overline{c}}}4\pi G{\displaystyle \frac{c^2}{\overline{c}}}(\rho +P_r)`$ (93)
$`{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{d\overline{c}}{dr}}\right)={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{\overline{𝒱}}{c}}+4\pi Gc(\rho +P_r+4P_t)`$ (94)
$`{\displaystyle \frac{1}{\overline{b}^2}}{\displaystyle \frac{dc}{dr}}{\displaystyle \frac{d\overline{c}}{dr}}1={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{\overline{𝒱}}{\overline{b}}}+8\pi Gc^2P_r`$ (95)
A consequence of these equations are the conservation equation of the matter energy-momentum tensor:

$$2\frac{dP_r}{dr}+\frac{d\mathrm{log}c}{dr}(3P_r\rho 4P_t)+\frac{d\mathrm{log}\overline{c}}{dr}(\rho +P_r)=0$$
(96)
Note that the structure of the equations of motion in vacuum ($`\rho =P_r=P_t=0`$) is relatively simple: Eq. (95) is an algebraic equation for $`\overline{b}`$ whose solution yields $`\overline{b}=B(dc/dr,d\overline{c}/dr,c,\overline{c},r)`$. Inserting this expression in the first two equations then yields two second-order, non-linear ODE’s for the radial evolution of $`c,\overline{c}`$.

We note that it is sometimes convenient to reformulate the field equations by changing the radial variable from $`r`$ to the radial analog of the “proper time”, namely $`\stackrel{~}{r}`$ defined by $`d\stackrel{~}{r}\overline{b}dr`$. Then the set of field equations read:

$`{\displaystyle \frac{d^2c}{d\stackrel{~}{r}^2}}={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{\overline{𝒱}}{\overline{c}}}4\pi G{\displaystyle \frac{c^2}{\overline{c}}}(\rho +P_r)`$ (97)
$`{\displaystyle \frac{d^2\overline{c}}{d\stackrel{~}{r}^2}}={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{1}{\overline{b}}}{\displaystyle \frac{\overline{𝒱}}{c}}+4\pi Gc(\rho +P_r+4P_t)`$ (98)
$`{\displaystyle \frac{dc}{d\stackrel{~}{r}}}{\displaystyle \frac{d\overline{c}}{d\stackrel{~}{r}}}1={\displaystyle \frac{m^2}{8}}{\displaystyle \frac{\overline{𝒱}}{\overline{b}}}+8\pi Gc^2P_r`$ (99)
$`{\displaystyle \frac{dr}{d\stackrel{~}{r}}}={\displaystyle \frac{1}{\overline{b}}}`$ (100)
Note the appearance of a fourth equation in the system linked to the definition of $`\stackrel{~}{r}`$ . The pressure-balance equation takes the same form as Eq. (96) with the replacement $`drd\stackrel{~}{r}`$.

The use of the variable $`\stackrel{~}{r}`$ corresponds to writing the two metrics as:

$`ds^2=a^2(\stackrel{~}{r})dt^2+{\displaystyle \frac{d\stackrel{~}{r}^2}{a^2(\stackrel{~}{r})}}+c^2(\stackrel{~}{r})d\mathrm{\Omega }^2`$ (101)
$`ds_{\mathrm{fl}}^2=dt^2+{\displaystyle \frac{d\stackrel{~}{r}^2}{\overline{b}^2(\stackrel{~}{r})}}+r^2(\stackrel{~}{r})d\mathrm{\Omega }^2`$ (102)
### 6.1 The General Relativity $`m=0`$ limit

To get some familiarity with the $`c`$, $`\overline{c}`$, $`\overline{b}`$ variables let us see how one can derive the Schwarzschild solution in this new language. In GR $`\overline{b}`$ is a gauge variable, and we can set $`\overline{b}=1`$ for the exterior solution which implies that, after fixing an integration constant, we can write $`\stackrel{~}{r}=r`$. Alternatively, we can start from the $`\stackrel{~}{r}`$-formulation and set to zero the terms $`𝒪(m^2)`$. This yields the very simple evolution equations $`d^2c/d\stackrel{~}{r}^2=0=d^2\overline{c}/d\stackrel{~}{r}^2`$, submitted to the constraint $`(dc/d\stackrel{~}{r})(d\overline{c}/d\stackrel{~}{r})=1`$. The general solution is $`c=\alpha \stackrel{~}{r}+\beta ,\overline{c}=\gamma \stackrel{~}{r}+\delta `$, with the constraint $`\alpha \gamma =1`$. The solution depends on three arbitrary constants. However, we can use the remaining rigid gauge symmetries of the problem ($`\stackrel{~}{r}^{}=a\stackrel{~}{r}+b,t^{}=ct`$) to set $`\alpha `$ to 1 and $`\beta `$ to zero. This leaves only one physical integration constant, $`\delta R_S`$, in terms of which the solution finally reads:

$$c=\stackrel{~}{r}=r,\overline{c}=\stackrel{~}{r}R_S=rR_S$$
(103)
For the interior solution, and for a constant density star, we choose the gauge variable to have the form:

$$\overline{b}=\frac{3}{2}\sqrt{\frac{1\frac{R_S}{R_{}}}{1\frac{R_S}{R_{}^3}r^2}}\frac{1}{2}$$
(104)
so that we continue having $`c=r`$. Then the $`\overline{c}`$ function reads:

$$\overline{c}=r\left(\frac{3}{2}\sqrt{1\frac{R_S}{R_{}}}\frac{1}{2}\sqrt{1\frac{R_S}{R_{}^3}r^2}\right)^2$$
(105)
Note that the interior solutions can be obtained from the ones of the $`\lambda `$, $`\mu `$, $`\nu `$ gauge \[Eqs. (39), (41)\], using the connecting relations (91) with $`\mu =0`$.

### 6.2 The perturbative limit for $`m0`$

We can look again how the results change in the linearized approximation for the massive case in terms of the $`c`$, $`\overline{c}`$, $`\overline{b}`$. In this case we linearize as:

$$c=r+\delta c,\overline{c}=r+\delta \overline{c},\overline{b}=1+\delta \overline{b}$$
(106)
We will only present here the exterior star solution to verify the correspondence with the $`\lambda `$, $`\mu `$, $`\nu `$ gauge. The vacuum equations (97), (98), (99), (100), using the $`r`$ variable become:

$`\delta c^{\prime \prime }\delta \overline{b}^{}={\displaystyle \frac{m^2}{2}}(r\delta \overline{b}+\delta c\delta \overline{c})`$ (107)
$`\delta \overline{c}^{\prime \prime }\delta \overline{b}^{}={\displaystyle \frac{m^2}{2}}(3r\delta \overline{b}+3\delta c+\delta \overline{c})`$ (108)
$`\delta c^{}+\delta \overline{c}^{}2\delta \overline{b}={\displaystyle \frac{m^2}{2}}r(3\delta c+\delta \overline{c})`$ (109)
for all four initial potentials (5), (6), (7), (8). From Eq. (109) we can solve for $`\delta \overline{b}`$ and substitute back to Eqs. (107), (108). Adding the resulting equations we obtain:

$$r(\delta c^{}\delta \overline{c}^{})+(1+m^2r^2)\delta \overline{c}(13m^2r^2)\delta c=0$$
(110)
Differentiating the above equation and substituting the difference $`\delta c^{\prime \prime }\delta \overline{c}^{\prime \prime }`$ from Eq. (107) (with $`\delta \overline{b}`$ substituted), we obtain a second independent first order equation:

$$r(2\delta c^{}+\delta \overline{c}^{})+\frac{1}{4}(2m^2r^2)\delta \overline{c}+\frac{1}{4}(223m^2r^2)\delta c=0$$
(111)
From Eqs. (110), (111) we can solve for $`\delta \overline{c}`$ and substitute to any of them. In this way we get a second order differential equation for $`\delta c`$:

$$\delta c^{\prime \prime }+\frac{4}{r(2+m^2r^2)}\delta c^{}\frac{2+m^2r^2}{r^2}\delta c=0$$
(112)
while $`\delta \overline{c}`$ and $`\delta \overline{b}`$ are then given as:

$$\delta \overline{c}=3\delta c\frac{4r}{2+m^2r^2}\delta c^{},\delta \overline{b}=\frac{1}{2}\delta c^{}+\delta \overline{c}^{}\frac{m^2}{4}r(3\delta c+\delta \overline{c})$$
(113)
Thus, although the system (107), (108), (109) at first sight seems to involve four integration constants, in fact only two of them are independent. Finally, we can integrate Eq. (112) and write the general solution decaying at infinity as:

$`\delta c={\displaystyle \frac{C_1}{4}}\left(1+{\displaystyle \frac{1}{mr}}+{\displaystyle \frac{1}{(mr)^2}}\right)e^{mr}`$ (114)
$`\delta \overline{c}={\displaystyle \frac{C_1}{4}}\left(3+{\displaystyle \frac{1}{mr}}+{\displaystyle \frac{1}{(mr)^2}}\right)e^{mr}`$ (115)
$`\delta \overline{b}={\displaystyle \frac{C_1}{2r}}\left(1+{\displaystyle \frac{1}{mr}}+{\displaystyle \frac{1}{(mr)^2}}\right)e^{mr}`$ (116)
where we have set the second integration constant corresponding to the exponentially growing solution to zero. It is easy to see, from the connecting relations (91) and after defining $`R=re^{\mu /2}`$, that these solutions are identical with the ones obtained in the $`\lambda `$, $`\mu `$, $`\nu `$ gauge \[Eqs. (47), (48), (49)\], with the same integration constant $`C_1`$. The latter, as we saw in section 4.2, for $`R_{}m^1`$ is $`C_1\frac{4}{3}R_S`$.

## 7 Singularity in the general asymptotically flat solution

Let us now discuss the central question of whether non-linear effects do cure or not the “discontinuity” (as $`m^20`$) with respect to GR exhibited by the linearized approximation. We recall that the specific conjecture made by Vainshtein in this respect was that the non-linear “dressing” of the unique perturbative solution (defined by resumming Poincaré-covariant diagrams, as sketched in ) leads to a solution of the type constructed in section 5. There are several ways one can take to investigate this conjecture. Either one can start from the center of the star with the requirement that it is locally flat and integrate outwards, or start from “infinity” (practically from some distance $`m^1`$) with the requirement of asymptotic flatness and integrate inwards. In addition, we can work in several different gauges.

Our work had to be necessarily numerical because it is a difficult dynamical question to control the global structure of a non-linear theory such as massive gravity. Though we have done simulations both ways (inwards and outwards), and with several gauges ($`\lambda `$, $`\mu `$, $`\nu `$, and $`a,b,c`$), we shall mainly report on the results obtained by integrating inwards, in the $`a,b,c`$ gauge. This way of proceeding has indeed two advantages: (i) the solution is essentially unique (once the equation of state of the matter is fixed) in the sense that its “initial state” at infinity is physically determined by giving oneself only one dimensionless parameter, $`mR_S`$, and (ii) the (inwards) radial evolution of the convenient variables $`c`$, $`\overline{c}`$, $`\overline{b}`$ is similar to the cosmological evolution studied in , which allows us to draw some intuition about its qualitative behaviour. Note that, by contrast, starting the evolution from the center obliges one to consider a two parameter family of solutions. One then needs to “shoot” from the center with a one-parameter family of initial data until one eventually matches the unique exponentially decaying solution at infinity.

Let us start by discussing the exterior region. \[As we shall see this suffices to conclude negatively.\] The dynamics of the $`c`$, $`\overline{c}`$, $`\overline{b}`$ variables is given by the Lagrangian:

$$=\frac{1}{2G}\left(\frac{\dot{c}\dot{\overline{c}}}{\overline{b}}+\overline{b}\right)\frac{m^2}{16G}\overline{𝒱}(\overline{b},c,\overline{c},r)$$
(117)
We start at infinity with the solution (114), (115), (116). For this solution, the “kinetic term” $`\dot{c}\dot{\overline{c}}`$ is positive (and near + 1), which means that the analogous “relativistic particle” initially moves on a spacelike worldline. Let us recall that, as $`\overline{b}`$ has no kinetic term, the value of $`\overline{b}`$ at any “moment” (i.e. at any radius $`r`$) is determined by extremizing the Lagrangian (117) considered as a function of $`\overline{b}`$. As long as the sign $`\dot{c}\dot{\overline{c}}>0`$ does not change, the first two terms in the Lagrangian (117) would suffice to “confine” $`\overline{b}`$ around a unique extremum $`\sqrt{\dot{c}\dot{\overline{c}}}`$. If ever the dynamics tends to decrease the value of $`\dot{c}\dot{\overline{c}}`$, this will tend to drive $`\overline{b}`$ towards smaller values. In fact the dynamics is due to the interplay between the “kinetic terms” in (117) and the potential terms and we must also take into account the $`\overline{b}`$dependence of $`\overline{𝒱}(\overline{b},c,\overline{c},r)`$. This dependence is somewhat complicated, and has a form which depends on the specific mass term considered. For our four fiducial mass terms $`𝒱^{(i)}`$, $`i=1,2,3,4`$ we have:

$`\overline{𝒱}^{(1)}=\overline{b}^3F^{(1)}+\overline{b}G^{(1)}`$ (118)
$`\overline{𝒱}^{(2)}=\overline{b}^2F^{(2)}+G^{(2)}`$ (119)
$`\overline{𝒱}^{(3)}={\displaystyle \frac{1}{\overline{b}}}F^{(3)}+\overline{b}G^{(3)}`$ (120)
$`\overline{𝒱}^{(4)}={\displaystyle \frac{1}{\overline{b}^2}}F^{(4)}+G^{(4)}`$ (121)
where, $`F^{(i)}=r^2F^{(i)}(\frac{c}{r},\frac{\overline{c}}{r})`$, $`G^{(i)}=r^2G^{(i)}(\frac{c}{r},\frac{\overline{c}}{r})`$. However, this structure is not enough for concluding about the “confining properties” of the potential terms, because the coefficients of the various powers of $`\overline{b}`$ above do not have a definite sign. Indeed, the explicit form of these functions is:

$`F^{(1)}={\displaystyle \frac{c^2}{r^2}}F^{(2)}=c^2\left(12+12{\displaystyle \frac{c^2}{r^2}}2{\displaystyle \frac{c^4}{r^4}}+6{\displaystyle \frac{\overline{c}}{c}}4{\displaystyle \frac{c\overline{c}}{r^2}}\right)`$ (122)
$`G^{(1)}={\displaystyle \frac{c^2}{r^2}}G^{(2)}=c^2\left(2+6{\displaystyle \frac{c}{\overline{c}}}4{\displaystyle \frac{c^3}{r^2\overline{c}}}\right)`$ (123)
$`F^{(3)}={\displaystyle \frac{c^2}{r^2}}F^{(4)}=c^2\left(12+12{\displaystyle \frac{r^2}{c^2}}2{\displaystyle \frac{r^4}{c^4}}+6{\displaystyle \frac{c}{\overline{c}}}4{\displaystyle \frac{r^2}{c\overline{c}}}\right)`$ (124)
$`G^{(3)}={\displaystyle \frac{c^2}{r^2}}G^{(4)}=c^2\left(2+6{\displaystyle \frac{\overline{c}}{c}}4{\displaystyle \frac{r^2\overline{c}}{c^3}}\right)`$ (125)
We can, however, have an analytical idea of the natural tendency of the dynamics by looking at the linearized solution (114), (115), (116) which results from the combined effect of kinetic and potential terms. If we look at Eq. (116) (remembering that $`C_1>0`$) we see that all the terms in the variation of $`\overline{b}`$ are negative, i.e. that $`\overline{b}`$ tends to decrease from its initial value $`\overline{b}=1`$. In addition one can deduce from Eqs. (114), (115) that the crucial quantity $`\dot{c}\dot{\overline{c}}`$, after increasing, starts decreasing when $`r\stackrel{<}{}m^1`$. As we said above, such a decrease of $`\dot{c}\dot{\overline{c}}`$ tends to further drive $`\overline{b}`$ down.

If we first consider the mass terms $`i=1,2`$, we see that they contain only positive powers of $`\overline{b}`$. Whatever be the sign of the coefficients of these terms, such functions are quite inefficient in preventing $`\overline{b}`$ from decreasing all the way towards zero, if the tendency of the “confining” term $`\dot{c}\dot{\overline{c}}/\overline{b}`$ in (117) is to drive $`\overline{b}`$ there because $`\dot{c}\dot{\overline{c}}`$ happens to decrease below zero. Therefore, in the cases $`i=1,2`$ it might a priori happen that the dynamics drives $`\overline{b}`$ toward zero at some finite radius, similarly to what we found in the cosmological study  where mass terms with weakly confining features tended to drive the gauge function $`e^\gamma `$ to zero in a finite time. The numerical simulations we performed for the potentials $`i=1,2`$, with various values of $`mR_S`$ (in the range $`10^5mR_S10^1`$) have shown that indeed such a behaviour is generic: we systematically find that the inward radial evolution ends up in a singularity999 Within the context of massive gravity the value $`\overline{b}=0`$ corresponds to an invariant singularity of the bi-gravity configuration $`(𝐠,𝐟)`$ because one of the eigenvalues of $`𝐟^1𝐠`$ is $`1/\overline{b}`$. where the variable $`\overline{b}`$ tends to zero, while at the same time $`\dot{c}\dot{\overline{c}}0`$.

We have also run similar simulations with the other mass terms $`i=3,4`$, and we found the same singular behaviour: $`\overline{b}0`$ at a finite radius, while, at the same time, $`\dot{c}\dot{\overline{c}}0`$. In the latter cases, the apparent “potential barriers” $`F^{(3)}/\overline{b}`$ and $`F^{(4)}/\overline{b}^2`$ that might have helped to prevent $`\overline{b}`$ to tend to zero turn out to be ineffective, because the (non-positive definite) coefficients are found to be driven by the dynamics towards zero at the same time that $`\dot{c}\dot{\overline{c}}0`$.

Our final conclusion (based on our numerical simulations) is therefore that for all values $`10^5mR_S10^1`$ (and probably also for all smaller values $`mR_S<10^5`$) the unique, asymptotically flat SSS solution (which is well defined near infinity) develops a singularity at some finite radius. The latter radius depends on the mass potential and is roughly between $`R_p`$ and $`m^1`$.

To complement our numerical study, we also performed simulations (still in the $`a,b,c`$ gauge) that start at the center of a star. This case is more delicate because generic (regular) initial data at $`r=0`$ depend now on two arbitrary parameters, and we expect that only a “line” of data can evolve into an asymptotically flat solution. If we consider a weakly self-gravitating star, it is relatively easy to choose initial data at the center (with $`\overline{b},\dot{c}`$ and $`\dot{\overline{c}}`$ all near one) that satisfy the constraint (95). Then, when considering the evolution system with respect to the “proper radius” $`\stackrel{~}{r}`$, general theorems on the continuity (with respect to small parameters appearing in the coefficients) of solutions of ODE systems guarantee that, if $`m^2R_{}^2`$ is small enough, any regular interior solution in the GR limit will be smoothly deformed into some massive gravity interior one. This provides a way to exhibit classes of local solutions of massive gravity which are continuously connected to GR solutions. However, a hard dynamical question (which is not covered by the local continuity theorems for ODE’s) is to know on which length scale these solutions stay close to GR solutions, and what are their asymptotic behaviour at infinity. We have addressed this dynamical question by numerical simulations (in the $`a,b,c`$ gauge). We found that all our simulations evolved a singularity at a finite radius (where $`\overline{b}`$ run away toward zero). This singularity occurred at distances between $`R_p`$ and $`m^1`$ depending on the mass term. In view of this singularity, we could not address the issue of “shooting” the initial data with the aim of matching an exponentially decaying solution at infinity.

We have also performed simulations in the $`\lambda `$, $`\mu `$, $`\nu `$ gauge. Their results have again all be negative: all solutions starting from the center developed some strange behaviour for large radii, and never showed any tendency to match the perturbative solution at infinity. Let us describe some of our results in the $`\lambda `$, $`\mu `$, $`\nu `$ gauge.

In view of the nature of the ODE system in this gauge we a priori need to give initial data, at $`R=0`$, for $`\lambda (0)`$, $`\nu (0)`$, $`\mu (0)`$ and $`\mu ^{}(0)`$. First of all, we know that in order not to have a conical singularity we should have $`\lambda (0)=0`$. Second, if we remember that the function $`\mu `$ defines the link between the “flat” and the “curved” radial coordinates $`r=Re^{\mu (R)/2}`$ , the regularity (in Cartesian-like coordinates) of this link implies that we must also require $`\mu ^{}(0)=0`$. Finally, we have only two initial data: say $`C\nu (0)`$ and $`D\mu (0)`$.

The leading terms for the three functions are:

$`\lambda =AR^2+\mathrm{}`$ (126)
$`\nu =C+BR^2+\mathrm{}`$ (127)
$`\mu =D+ER^2+\mathrm{}`$ (128)
In GR the quantities $`A`$ and $`B`$ are related to the density $`\rho _0`$ (for a constant density star) and the central pressure $`P_c`$ as:

$`3A_{GR}=8\pi G\rho _0`$ (129)
$`2B_{GR}A_{GR}=8\pi GP_c`$ (130)
while the constant $`C`$ is arbitrary. In the massive gravity case, we obtain the following relations:

$`3A=8\pi G\rho _0+m^2\mathrm{\Sigma }_1(C,D)`$ (131)
$`2BA=8\pi GP_c+m^2\mathrm{\Sigma }_2(C,D)`$ (132)
for all four potentials $`𝒱^{(i)}`$, $`i=1,2,3,4`$. The functions $`\mathrm{\Sigma }_1(C,D)`$ and $`\mathrm{\Sigma }_2(C,D)`$ are zero when we sit on the point $`C=D=0`$, but they may have non-trivial roots (we will see more in the next section). In order to study the conjecture of Vainshtein (which assumes that we are near GR around the source) we explored the region of the $`(C,D)`$ parameter space that is in the neighbourhood of $`(0,0)`$, namely $`0.1C0.1`$ and $`0.1D0.1`$.

For such a choice of initial data, one finds that, for small enough $`m`$, the solution stays very close to the GR result until a radius which is very roughly of order $`R_V`$ and then the solution starts to deviate considerably. Note that these numerical simulations provide still another way of constructing local solutions of massive gravity that are close to GR ones. Note also that, when using a matter density model where $`\rho `$ undergoes a finite jump at $`R=R_{}`$, the variables $`\lambda ^{},\nu ^{\prime \prime }`$ and $`\mu ^{\prime \prime }`$ suffer corresponding, correlated jumps.

However, even when changing in a continuous manner the initial data, we never found any solution which, at large radii, tend to match the perturbative one. In fact, all solutions tend to run into a numerically unstable behaviour, where the simulation stops after a finite radius. We saw no tendency that the functions at distances much bigger than $`R_p`$ match the perturbative solution ones. We could not run the simulation far enough and could not understand if this denotes that there is a singularity at some finite distance from the star. However, by numerically transforming the $`\lambda `$, $`\mu `$, $`\nu `$ variables into the $`c`$, $`\overline{c}`$, $`\overline{b}`$ ones we have found that at the point of numerical difficulties of the $`\lambda `$, $`\mu `$, $`\nu `$ gauge, there was a tendency for the $`c`$, $`\overline{c}`$, $`\overline{b}`$ variables to exhibit their usual singular behaviour ($`\overline{b}0`$, $`\dot{c}\dot{\overline{c}}0`$).

Summarizing: Our numerical simulations have falsified the conjecture put forward in  and . Far from helping to regularize the discontinuities and attendant $`𝒪(m^n)`$ terms exhibited by the linearized approximation, the non-linear effects in massive gravity aggravate the situation by developing singularities at a finite radius when starting from the unique, asymptotically flat solution.

This negative result closes, in our opinion, one possibility. However, it still leaves open several possibilities for eventually reconciling massive gravity with phenomenology. Indeed, our study so far has focused on only one type of possible asymptotic behaviour at infinity: namely, (16), i.e. the case where the physical metric asymptotes the “trivial” vacuum solution. In the next section, we shall point out that even the simple-minded mass terms (5)-(8) might admit other (translationally invariant) vacuum solutions. Finally, in the last section we shall consider the possibility that continuity be restored for solutions which asymptote some non-trivial cosmology at infinity.

## 8 Non-trivial translationally-invariant, spherically symmetric vacua

Returning to our original potentials, let us examine all their possible translationally invariant, spherically symmetric vacua. A translationally invariant vacuum is a solution (in absence of matter) for which, in some gauge, both $`g_{\mu \nu }`$ and $`f_{\mu \nu }`$ are constant (i.e. independent of the spacetime coordinates $`x^\mu `$). Requiring spherical symmetry will impose a further constraint on the symmetry group fixing the bi-metric configuration ($`g_{\mu \nu }`$, $`f_{\mu \nu }`$).

The first constraint translates in the statement that the matrix $`g_\nu ^\mu f^{\mu \sigma }g_{\nu \sigma }`$ can be globally diagonalized as diag$`(\lambda _0,\lambda _1,\lambda _2,\lambda _3)`$, where the eigenvalues $`\lambda _i`$ are constant. In other words, there exist coordinate systems where the two metrics read:

$$𝐟=dt^2+dx^2+dy^2+dz^2,𝐠=\lambda _0dt^2+\lambda _1dx^2+\lambda _2dy^2+\lambda _3dz^2$$
(133)
Requiring spherical symmetry then imposes the constraint $`\lambda _1=\lambda _2=\lambda _3`$. This leaves only two constants.

Finally, rewording the result in terms of the $`\lambda `$, $`\mu `$, $`\nu `$ gauge, one easily finds that the above requirements are equivalent to the constraints:

$$\lambda =0,e^\mu =\lambda _1=const.,e^\nu =\lambda _0=const.^{}$$
(134)
Technically, we are then left with looking for extrema of the mass terms with constant values of $`\lambda `$, $`\mu `$, $`\nu `$ of the type (134).

For the last three mass terms we found that they only admit the trivial vacuum $`(\lambda _1,\lambda _0)=(e^\mu ,e^\nu )=(1,1)`$.

On the other hand, we found that the first mass term $`𝒱^{(1)}`$ admits three possible vacua. This can be seen by looking at the $`\mathrm{\Sigma }_1(C,D)`$, $`\mathrm{\Sigma }_2(C,D)`$ functions in Eqs. (131), (132):

$`\mathrm{\Sigma }_1(C,D)={\displaystyle \frac{3}{4}}(e^D1)(e^D+3e^C2)`$ (135)
$`\mathrm{\Sigma }_2(C,D)={\displaystyle \frac{1}{4}}(5e^D(e^C3)+7e^{2D}3e^C+6)`$ (136)
remembering that $`\nu =C`$, $`\mu =D`$ in vacuum. Then the conditions $`\mathrm{\Sigma }_1(C,D)=0`$, $`\mathrm{\Sigma }_2(C,D)=0`$ have three solutions: the trivial vacuum $`(\lambda _1,\lambda _0)=(e^\mu ,e^\nu )=(1,1)`$, plus two new ones:

$$(\lambda _1,\lambda _0)=(e^\mu ,e^\nu )=(\frac{1}{2},\frac{1}{2}),(\lambda _1,\lambda _0)=(e^\mu ,e^\nu )=(\frac{3}{2},\frac{1}{6})$$
(137)
Note that the first non-trivial vacuum corresponds to $`g_{\mu \nu }=\frac{1}{2}f_{\mu \nu }`$. Such a conformal situation respects not only spherical symmetry, but in fact also the full Poincaré symmetry. By contrast, the last vacuum breaks the Poincaré symmetry to a mere Euclidean one.

Let us examine now the perturbations around these vacua for the first potential. For the trivial vacuum, we already know that the linearized perturbations feature a Pauli-Fierz mass term, i.e. a ratio $`\alpha =1`$ between the two terms, with a spin-2 mass, say $`m_2`$, identical to the mass parameter entering the action, i.e. $`m_2^2=m^2`$. For the $`(\frac{1}{2},\frac{1}{2})`$ vacuum, one finds that $`\alpha =\frac{1}{8}`$ and $`m_2^2=16m^2`$. The fact that we do not have a Pauli-Fierz mass term means, as we recalled above, that we have a ghost-like scalar field ($`h`$). In addition, the fact that the value of $`\alpha `$ is not between $`1/4`$ and $`1`$ means that this ghost degree of freedom is also tachyonic, since $`m_0^2=\frac{2}{3}m_2^2`$. As we said above, such excitations are pathological, and we do not wish to consider them (though we know that this is a cheap way of ensuring continuity as $`m^20`$).

Finally, the third vacuum $`(\frac{3}{2},\frac{1}{6})`$, does not give an equation of the type (17), but rather instead of the $`h_{\mu \nu }\alpha \eta _{\mu \nu }h`$ mass term, one finds the following contribution:

$$D_{\mu \nu }^{\kappa \lambda }h_{\kappa \lambda }+E_{\mu \nu }h$$
(138)
where the tensors $`D_{\mu \nu }^{\kappa \lambda }`$ and $`E_{\mu \nu }`$ are constructed from the two background metrics $`\overline{\eta }_{\mu \nu }`$ and $`\eta _{\mu \nu }`$, if $`g_{\mu \nu }=\overline{\eta }_{\mu \nu }+h_{\mu \nu }`$.

The important point is that these tensors for the background we are considering are not proportional to $`\delta _\kappa ^\mu \delta _\lambda ^\nu +(\mu \nu )`$ and $`\overline{\eta }_{\mu \nu }`$ respectively. This happens because the “metric condensate” is not proportional to the background of the observable metric, and therefore, as we said, breaks the underlying Poincaré invariance. Thus, as we see from Eq. (138), the equation of motion for the perturbation $`h_{\mu \nu }`$ will have different structure from Eq. (17). We leave to future work the task of redoing, in this new setting, the analysis of .

In the next section, we shall finally study other types of vacua, which admit a cosmological interpretation.

## 9 Continuous “non-co-diagonal” solutions with de Sitter asymptotics in massive gravity

In this section we explore another possibility for eventually reconciling massive gravity with phenomenology. This possibility is based on an interesting type of spherically symmetric stationary solutions of massive gravity found long ago by Salam and Strathdee . \[These types of solutions were then generalized to a full bigravity setting .\] These solutions tend, at spatial infinity, towards a de Sitter solution, instead of one of the translationally-invariant vacua considered above.

The aim of this section is two-fold: (i) to describe the “universality class”  of mass terms which admit such solutions, and (ii) to show how the vacuum solution of  can be extended to a globally regular spacetime representing a “star” embedded in a de Sitter background.

### 9.1 Spontaneous symmetry breaking in massive gravity

Up to this section, we have only considered spherically symmetric solutions which could be represented in a co-diagonal way, i.e. such that both $`𝐠`$ and $`𝐟`$ are diagonal in the same coordinate system. As we indicated, perturbation theory (around the trivial vacuum) necessarily select such a type of solution, if the matter sources are supposed to be “at rest”. However, non-perturbative solutions might introduce an interesting twist in the link between the two metrics $`𝐠`$ and $`𝐟`$. A convenient formalism for describing this possibility consists in introducing a “link field” $`Y(x)`$ (also called “Stueckelberg field”) relating the two metrics $`𝐠`$ and $`𝐟`$. This formalism first appeared in  and, in more detail in , and has been recently developed in . More precisely, instead of writing the theory of massive gravity in a coordinate system which is common to both metrics, and of writing the mass term in terms of the difference $`g_{\mu \nu }(x)f_{\mu \nu }(x)`$ between the two metric tensors (supposed to live on the same manifold, and expressed in the same coordinate system), one can explicitly introduce the map $`Y`$ between the $`𝐠`$-Riemannian manifold, and the $`𝐟`$-one. If we use independent coordinate systems on both manifolds (say $`x^\mu `$ on the $`𝐠`$ manifold $`_g`$, and $`X^\alpha `$ on the $`𝐟`$-one $`_f`$) the abstract map $`Y`$ takes the explicit form $`X^\alpha =Y^\alpha (x^\mu )`$. This map (from $`_g`$ towards $`_f`$) “pulls back” the metric $`𝐟`$ into an image metric $`Y_{}𝐟`$ which lives on $`_g`$. This leads to expressing the mass term in terms of:

$$H_{\mu \nu }[g,f,Y]=g_{\mu \nu }(x)_\mu Y^\alpha _\nu Y^\beta f_{\alpha \beta }(x)$$
(139)
Then the mass term (and therefore the action) becomes a functional of three fields: $`g,f`$ and $`Y`$. In the massive gravity context, $`𝐟`$ is a non-dynamical flat metric, and we can always choose to express it in Lorentzian coordinates, i.e. with $`f_{\alpha \beta }(x)=\eta _{\alpha \beta }`$. Then the action depends only on two dynamical fields: $`g`$ and $`Y`$. The usual formalism (used above) consists of using the gauge invariance of the theory to fix $`Y=`$ id (which destroys the explicit gauge invariance). Ref.  advocates the usefulness of not using this “unitary gauge”, but instead of working with an explicitly gauge invariant theory with two fields: $`g`$ and $`Y`$. The immediate problem with this proposal is then that the kinetic terms for the $`Y`$ field are very non-linear because they result from replacing (139) in the mass terms. \[For instance, the existence of a sixth degree of freedom is hidden in the highly non-linear nature of this kinetic term.\]

One can view $`Y`$ as a gravitational analogue of the Higgs field. Ref.  has indicated how to deal with perturbation theory around the “trivial Higgs configuration” $`Y_0^\alpha (x^\mu )=x^\alpha `$, by expanding the theory in powers of the “Goldstone” field $`\pi ^\alpha (x)`$ such that $`Y^\alpha (x^\mu )=x^\alpha +\pi ^\alpha (x)`$. The analogue, in this language, of the perturbative co-diagonality of two SSS metrics is the following. If we consider some static matter source (in $`_g`$), perturbation theory for the Goldstone field will never generate a time-component $`\pi ^0(x)`$ (which induces a non-zero mixed component $`g_{0i}(x)`$ in the physical metric).

In this language, one needs to be in a non-perturbative situation of “spontaneous symmetry breaking” to expect such a $`\pi ^0(x)`$ to develop. This situation was in fact recently considered, mutatis mutandis, in . There it was noticed that besides the stable family of bigravity cosmological solutions, where the two metrics $`𝐠_1`$ , $`𝐠_2`$ can be written in a time-orthogonal gauge ($`g_{0i\mathrm{\hspace{0.17em}1}}=0=g_{0i\mathrm{\hspace{0.17em}2}}`$) and where the remaining metric coefficients ($`g_{00},g_{ij}`$) depend only on time, there might also exist, under some conditions, families of “symmetry breaking” solutions where the relative “shift vector” $`b^i`$ (proportional to $`g_2^{0i}`$ in the gauge where $`g_{0i\mathrm{\hspace{0.17em}1}}=0`$) is not equal to zero. The possible non-zero values of $`b^i`$ were the extrema of the mass term $`V(𝐛)V_0+a𝐛^2+b𝐛^4+𝒪(𝐛^6)`$. Similarly to the Higgs mechanism, for certain shapes of the function $`V(𝐛)`$ (e.g. when $`a<0`$ and $`b>0`$) there might exist, besides the trivial extremum $`𝐛=0`$, non-trivial extrema ($`𝐛0`$) of $`V(𝐛)`$. The precise condition for this to happen was studied in the Appendix B of . We are going to recover the same condition for the possibility of a similar symmetry breaking in the context of SSS solutions.

In this context, the co-diagonal form (28), (29) is not the most general form compatible with spherical symmetry and stationarity. The most general form depends on 4 functions, and can be written in the following $`\lambda `$, $`\mu `$, $`\nu `$, $`\alpha `$ gauge:

$`ds^2=e^{\nu (R)}dT^2+e^{\lambda (R)}dR^2+R^2d\mathrm{\Omega }^2`$ (140)
$`ds_{\mathrm{fl}}^2=(dT+\alpha dR)^2+\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2e^{\mu (R)}dR^2+e^{\mu (R)}R^2d\mathrm{\Omega }^2`$ (141)
Here all the functions $`\lambda `$, $`\mu `$, $`\nu `$, $`\alpha `$ depend only on $`R`$. The above gauge has the advantage of separating the variables into two groups: the “physical” variables $`\lambda `$, $`\nu `$ (which are directly observable in gravitational experiments), and the “gauge” variables $`\mu `$, $`\alpha `$ which enter only the unobservable background flat metric. The gauge variables enter the action (4) only through the mass term $`(m^2/4)𝒱(𝐟^1𝐠)`$. Therefore one can immediately obtain the equations of motion of $`\mu `$, $`\alpha `$ by varying only the mass term:

$`{\displaystyle \frac{\delta }{\delta \mu }}𝒱(𝐟^1𝐠)=0`$ (142)
$`{\displaystyle \frac{\delta }{\delta \alpha }}𝒱(𝐟^1𝐠)=0`$ (143)
Equation (142) above generalizes the constraint (35) to the case where $`\alpha 0`$. Let us consider here the new equation (143). A crucial point is that $`\alpha `$ enters $`𝒱`$ only algebraically. Therefore, similarly to what we recalled above concerning the dynamics of the cosmological shift vector $`𝐛`$, non-trivial SSS solutions will exist only if $`𝒱(\alpha )`$ admits non-trivial extrema with respect to $`\alpha `$. To understand what this implies for the “universality class” of such mass terms, let us relate the $`\alpha `$ dependence on the eigenvalues of the gravitational energy-momentum tensor $`T_{\mu \nu }^{(g)}`$. Let us recall (from ) the following bi-geometrical facts: The bi-gravity configuration $`(𝐠,𝐟)`$ defines, at each spacetime point, a preferred moving frame with respect to which $`𝐠=\mathrm{diag}(\lambda _0,\lambda _1,\lambda _2,\lambda _3)`$ and $`𝐟=\mathrm{diag}(1,1,1,1)`$. \[Note that this preferred frame is not proportional to the coordinate frame used in the normal gauges such as the $`\lambda `$, $`\mu `$, $`\nu `$, $`\alpha `$ gauge.\] The eigenvalues $`\lambda _a`$ of $`𝐟^1𝐠`$ depend on $`R`$. They can be easily computed from the metric coefficients $`\lambda `$, $`\mu `$, $`\nu `$, $`\alpha `$. For instance:

$$\lambda _2=\lambda _3=e^{\mu (R)}$$
(144)
Let us focus on the $`\alpha `$ dependence of the eigenvalues. Note that $`\alpha `$ disappears when calculating the determinant of $`𝐟`$. Therefore the product:

$$\mathrm{\Delta }\lambda _0\lambda _1$$
(145)
which is the ratio of the two $`2\times 2`$ determinants of the $`(T,R)`$ sub-blocks of the two metrics, does not depend on $`\alpha `$. Indeed,

$$\mathrm{\Delta }=\frac{e^{\nu +\lambda +\mu }}{(1R\mu ^{}/2)^2}$$
(146)
On the other hand, $`\alpha `$ enters the trace of the $`(T,R)`$ sub-block of $`𝐠^1𝐟`$. Calculating this trace yields:

$$\frac{1}{\lambda _0}+\frac{1}{\lambda _1}=e^\nu +e^\lambda [(1\frac{R\mu ^{}}{2})^2\alpha ^2]=\frac{\lambda _0+\lambda _1}{\lambda _0\lambda _1}$$
(147)
When varying only $`\alpha `$, one therefore finds that $`\lambda _2=\lambda _3`$ and $`\lambda _0\lambda _1`$ do not change, but that $`\lambda _0+\lambda _1`$ varies proportionally to $`\alpha \delta \alpha `$. On the other hand, the variation of the mass term $`S^{(g)}𝒱`$ in the action is related, by definition, to the logarithmic changes of the eigenvalues by $`\delta S^{(g)}=(1/2)\sqrt{g}T_{}^{(g)}{}_{a}{}^{a}\delta \mathrm{log}\lambda _a`$, where $`T_{}^{(g)}{}_{a}{}^{a}`$ denote the eigenvalues (in the preferred frame which is easily seen to diagonalize also $`T_{\mu \nu }^{(g)}`$) of the gravitational energy tensor. The final result is that the variational derivative (143) is found to be proportional to the product:

$$\alpha \frac{T_{}^{(g)}{}_{0}{}^{0}T_{}^{(g)}{}_{1}{}^{1}}{\lambda _0\lambda _1}$$
(148)
As we said above, a universal solution of this $`\alpha `$-constraint is the trivial case $`\alpha =0`$. The spontaneous symmetry breaking (SSB) situation corresponds to the case where it is the second factor above which vanishes. Using the formulas given in  one can easily see that the condition for SSB does not depend on the choice of density prefactor in front of $`𝒱wV`$, where $`w`$ can be either $`\sqrt{f}`$, $`\sqrt{g}`$, $`(fg)^{1/4}`$ etc… Finally, we can write the condition for SSB in terms of the scalar potential $`V(\lambda _0,\lambda _1,\lambda _2,\lambda _3)𝒱/w`$ as the existence of solutions for the equation:

$$\rho (\lambda _0,\lambda _1,\lambda _2,\lambda _3)\frac{1}{\lambda _0\lambda _1}(\frac{V}{\mathrm{ln}\lambda _0}\frac{V}{\mathrm{ln}\lambda _1})=0$$
(149)
For instance, if we consider the class of potentials $`V`$ which depend only on $`\sigma _1`$ and $`\sigma _2`$, where $`\sigma _n\mathrm{\Sigma }_a(\mathrm{ln}\lambda _a)^n`$, we recover the condition $`_{\sigma _2}V=0`$ derived in .

It happens that all the mass terms written in (5)-(8) can satisfy the SSB condition. First, we note that they contain only two independent “scalar” potentials $`V`$. Moreover, the two types of scalar potentials can be mapped into each other by the exchange $`𝐠𝐟`$, which is equivalent to the exchange $`\lambda _a1/\lambda _a`$. Therefore, it suffices to consider the potential $`V^{(1)}`$ defined, say, by:

$$V^{(1)}\frac{1}{2}\left\{\mathrm{tr}[((𝐠𝐟)𝐟^1)^2](\mathrm{tr}[(𝐠𝐟)𝐟^1])^2\right\}=\frac{1}{2}[\underset{a}{}(\lambda _a1)^2(\underset{a}{}(\lambda _a1))^2]$$
(150)
The computation of the quantity $`\rho `$ defined above then yields:

$$\rho =1\underset{c0,1}{}(\lambda _c1)=d\underset{c0,1}{}\lambda _c$$
(151)
where $`d`$ denotes the space dimension. In our case, it is $`d=3`$, and we can also use $`\lambda _2=\lambda _3`$ to write $`\rho =32\lambda _2`$.

Finally, we conclude that for the mass terms $`𝒱^{(1)}`$ and $`𝒱^{(2)}`$:

$$\rho [𝒱^{(1)}]=\rho [𝒱^{(2)}]=32\lambda _2$$
(152)
while the corresponding results for the other ones (obtained by exchanging $`\lambda _a1/\lambda _a`$) read:

$$\rho [𝒱^{(3)}]=\rho [𝒱^{(4)}]=3+2/\lambda _2$$
(153)
Therefore the SSB condition (i.e. the equation of motion for $`\alpha 0`$) can be satisfied for all four potentials. In the first two cases, it fixes the value of $`\lambda _2`$ , i.e. of $`e^\mu `$, to be simply $`3/2`$, while in the last two cases one must have $`\lambda _2=e^\mu =2/3`$. \[The corresponding values in space dimension $`d`$ would be $`d/(d1)`$ and $`(d1)/d`$, respectively.\] It is a happy technical accident (which simplifies the subsequent derivations) that these values are constant. This is, however, not crucial. The essential point is that the SSB condition $`\rho (\lambda _a)=0`$ admits solutions. More general types of SSB potentials might admit solutions of $`\rho (\lambda _0,\lambda _1,\lambda _2)=0`$ along some curved hypersurface in $`\lambda _a`$-space, say $`\lambda _2=f(\lambda _0,\lambda _1).`$

### 9.2 Non-codiagonal solutions in massive gravity

Let us now study the global structure of SSB solutions (both in the exterior and the interior). There are two ways of proceeding. Either we continue using a “geometrical” approach based on the systematic use of the invariant eigenvalues, or we use a more pedestrian approach based on writing down explicitly all the field equations in some gauge. Before using a more pedestrian approach let us briefly indicate the results of the “geometrical” one in the simple exterior case.

Having obtained the value of $`\lambda _2=e^\mu `$ we can now get a simple constraint on the product of eigenvalues $`\mathrm{\Delta }\lambda _0\lambda _1`$. Indeed, we found above (as a consequence of the SSB condition) that the $`2\times 2`$ $`(T,R)`$ block of the gravitational energy tensor satisfies (in its diagonalizing frame) $`T_{}^{(g)}{}_{0}{}^{0}=T_{}^{(g)}{}_{1}{}^{1}`$. This condition actually means that, in any frame, the $`2\times 2`$ $`(T,R)`$ block of $`T_{\mu \nu }^{(g)}`$ is proportional the $`2\times 2`$ unit matrix. When going back to the physical gauge $`\lambda `$, $`\mu `$, $`\nu `$, $`\alpha `$ this means that the equation of state of the gravitational energy satisfy $`P_r^{(g)}=\rho ^{(g)}`$. Using this information in the usual Einstein equations (in the exterior) (31), (32) gives us the simple constraint that the product $`e^\nu e^\lambda `$, i.e., modulo a constant, the quantity $`\mathrm{\Delta }\lambda _0\lambda _1`$, is constant. At this stage we know all the eigenvalues modulo the knowledge of, say, $`\lambda _1(R)`$. To get the radial variation of $`\lambda _1(R)`$ it suffices to write down the conservation equation (36) for the gravitational part of the energy-momentum tensor. This is simplified by the result $`P_r^{(g)}=\rho ^{(g)}`$ and now reads

$$_RP_r^{(g)}+\frac{2}{R}(P_r^{(g)}P_t^{(g)})=0$$
(154)
As the components $`(P_r^{(g)},P_r^{(g)},P_t^{(g)},P_t^{(g)})`$ of the gravitational energy-momentum tensor are checked to depend only on the eigenvalues, the above equation gives an evolution equation for $`\lambda _1(R)`$. Actually, one can see immediately a particular solution: namely $`\lambda _1=\lambda _2`$. \[One then checks that this is the only possible solution.\] Indeed, the isotropy condition $`\lambda _1=\lambda _2=\lambda _3`$ guarantees that $`P_r^{(g)}=P_t^{(g)}`$. Then we get a solution of the above equation because all eigenvalues are constant (indeed, we know already that $`\lambda _2=\lambda _3`$ is constant, and that the product $`\lambda _0\lambda _1`$ is also constant).

Finally, we end up (in the exterior region) with 4 constant eigenvalues, three of which are determined to be $`\lambda _1=\lambda _2=\lambda _3=3/2`$ , or $`2/3`$ (depending on the choice of mass term), and the fourth $`\lambda _0`$ being an arbitrary constant (with $`\lambda _0\lambda _1=\mathrm{\Delta }=const.`$). Though the set of eigenvalues introduces a dissymmetry between the (eigen-)time direction and the space directions, there is a further simplification which comes (again because of the SSB constraint) from the result above $`T_{}^{(g)}{}_{0}{}^{0}=T_{}^{(g)}{}_{1}{}^{1}`$. This condition, together with the equality of the three “spatial” eigenvalues, shows that, finally $`T_{}^{(g)}{}_{0}{}^{0}=T_{}^{(g)}{}_{1}{}^{1}=T_{}^{(g)}{}_{2}{}^{2}=T_{}^{(g)}{}_{3}{}^{3}`$. Therefore the gravitational energy tensor is fully isotropic in spacetime: $`T_{}^{(g)}{}_{\mu }{}^{\nu }=\mathrm{\Lambda }\delta _\mu ^\nu `$, where the constant $`\mathrm{\Lambda }`$ is some function of the eigenvalues. Such a gravitational energy tensor is equivalent to a cosmological constant $`\mathrm{\Lambda }`$. Therefore we conclude, without calculations, that the only SSB spherically symmetric stationary solutions of massive gravity for the potentials (5)-(8) must be Schwarzschild-(Anti-)de Sitter solutions. This was, indeed, the result of .

Let us now generalize this exterior solution to a global solution, including the interior of a star. For concreteness, we shall consider the specific mass term $`𝒱^{(4)}`$, (8). The main message of this section would not be altered if we used any of the (5), (6) or (7) instead. We know already, from the previous reasoning, the general structure of the exterior solution. We only need to compute now the deviation from Schwarzschild-de Sitter (SdS) in the interior of the star. We could do the explicit calculations of the interior structure by working, as above, with the geometrical eigenvalues, and the corresponding gravitational energy tensor. We shall however use a more pedestrian approach, writing down the equations of motion in the gauge defined by the line elements (140), (141). \[In Appendix B we also present the results following the notation of the original derivation .\] For simplicity we will consider a star of isotropic pressure $`P`$. The $`(T,T)`$ and $`(R,R)`$ equations of motion, after setting $`e^\mu =2/3`$, then read:

$`e^\lambda \left[{\displaystyle \frac{\lambda ^{}}{R}}+{\displaystyle \frac{1}{R^2}}(e^\lambda 1)\right]=8\pi G\left({\displaystyle \frac{9m^2}{64\pi G}}\sqrt{{\displaystyle \frac{3}{2}}}e^{3(\lambda +\nu )/2}+\rho \right)`$ (155)
$`e^\lambda \left[{\displaystyle \frac{\nu ^{}}{R}}+{\displaystyle \frac{1}{R^2}}(1e^\lambda )\right]=8\pi G\left({\displaystyle \frac{9m^2}{64\pi G}}\sqrt{{\displaystyle \frac{3}{2}}}e^{3(\lambda +\nu )/2}+P\right)`$ (156)
The conservation of the matter energy-momentum tensor gives the familiar equation:

$$P^{}=\frac{\nu ^{}}{2}(\rho +P)$$
(157)
while from Eq. (142) (which replaces the previous constraint $`f_g=0`$ derived from the conservation of the energy-momentum tensor generated by the mass term) one gets an algebraic equation for $`\alpha `$ with solution:

$$\alpha ^2(R)=\frac{1}{4}e^\nu \left[3R(\lambda ^{}+\nu ^{})2(e^\lambda 1)(3e^\nu 2)\right]$$
(158)
Note that $`\alpha ^2(R)`$ undergoes a finite jump (downwards) as $`R`$ crosses the radius of the star (from the interior).

In the exterior, one trivially finds by adding Eqs. (155) and (156) that $`\mathrm{\Delta }^{}=\frac{2}{3}(e^{\lambda +\nu })^{}=0`$ which implies, as we said earlier, that $`\mathrm{\Delta }=\frac{2}{3}(e^{\lambda +\nu })=const.\mathrm{\Delta }_0>0`$. Then from Eq. (155) for example one finds the usual SdS solution:

$$e^\lambda =\frac{2}{3\mathrm{\Delta }_0}e^\nu =1\frac{R_S}{R}\frac{\mathrm{\Lambda }}{3}R^2,\mathrm{with}\mathrm{\Lambda }=\frac{3}{4}\frac{m^2}{\mathrm{\Delta }_0^{3/2}}$$
(159)
Then from Eq. (158) one can directly find $`\alpha `$:

$$\alpha ^2=\frac{2}{3\mathrm{\Delta }_0}(e^\lambda 1)\left(e^\lambda \frac{9\mathrm{\Delta }_0}{4}\right)$$
(160)
The requirement of positivity of $`\alpha ^2`$ then restricts $`0<\mathrm{\Delta }_0<4/9`$. Note also that at the position of the dS horizon, the gauge function $`\alpha `$ diverges since $`e^\lambda \mathrm{}`$. However, this is a mere coordinate singularity of the bigravity configuration $`𝐟`$, $`𝐠`$ as is expected from the fact that the eigenvalues of $`𝐟^1𝐠`$ are regular there, being equal, in the exterior, to $`(\lambda _0,\lambda _1,\lambda _2,\lambda _3)=(3\mathrm{\Delta }_0/2,2/3,2/3,2/3)`$. A coordinate transformation which regularizes both $`𝐟`$ and $`𝐠`$ is $`dTdT+\alpha (R)dR`$. \[Note that the black-hole solution of  is also regular at the event horizon by the same transformation, something not noticed in .\]

The exterior massive gravity solution above involves, as separate scales, only $`R_S`$ and $`\mathrm{\Delta }_0^{3/4}\lambda _m`$. Note, however, that (when $`\mathrm{\Delta }_01`$) the scale $`R_p`$ (2) enters indirectly the solution as the scale where the two separate contributions to $`e^\lambda =\frac{2}{3\mathrm{\Delta }_0}e^\nu `$ become comparable to each other.

In the interior of the star, one has to solve the closed system of differential equations (155), (156) and (157). Note that the contribution of the massive graviton in the interior equations of motion differs from a pure cosmological constant. A crucial feature of this system is that the graviton mass $`m^2`$ appears only with a positive power, and in front of a lower-derivative term. General theorems on the continuity (with respect to parameters) of solutions of ODE systems then guarantees that, if $`m^2R_{}^2`$ is small enough, any regular interior solution in the GR limit will be smoothly deformed into some massive gravity interior one (satisfying the regularity condition $`\lambda (R)R^2`$ at the center). The only requirement then that the solution is acceptable is that $`\alpha ^2`$ in the interior be positive. We have numerically verified that this criterion is not difficult to satisfy.

Finally, let us note that the only effect of the mass of the graviton in this kind of SSB solutions is the introduction of a cosmological tail with an associated scale $`\mathrm{\Delta }_0^{3/4}\lambda _m`$. Since $`\mathrm{\Delta }_0`$ is an arbitrary dimensionless integration constant (lying in the interval $`0<\mathrm{\Delta }_0<4/9`$), one obtains the peculiar result that the cosmological scale of the theory is undetermined. A priori, however, one should expect that $`\mathrm{\Delta }_0𝒪(1)`$ and that its precise value should finally be determined when trying to match this SSS solution with a cosmological one. Note also, that since the regular SSB solutions that we constructed differ from the GR ones in the presence of a cosmological constant $`\mathrm{\Lambda }`$ only inside the star, one does not have observable consequences of the type claimed in  that might test the model experimentally.

## 10 Conclusions

We have explored the issue of continuity in the purely massive gravity theory for several mass terms, which all reduce to a Pauli-Fierz one in the linearized limit. We have discussed in detail the claim of  that the resummation of non-linear effects make the massless limit smooth. We showed that the type of local expansion (considered in some bounded domain outside the source) proposed in  does not exist for the simplest, action-based, mass term considered in the subsequent, more detailed publication . However, it exists for other mass terms. We showed the double series nature of the expansion that was implicitly performed in  and noted that different mass terms give different kinds of expansions.

Then we discussed the global nature of massive gravity solutions. When considering the inward radial evolution of general asymptotically flat solutions, we find that they all end up in a singularity at some finite radius well outside the source. When, on the other hand, we consider the asymptotics of a two-parameter class of solutions which are close to GR in and near the source, we also find that they end up in a singularity at some large radius. These results were obtained in two different gauges, one of which allows one to qualitatively relate the singular runaway of the solution to a similar runaway found in recent cosmological studies of massive gravity .

These results falsify the continuity claim (as it was formulated) of , and vindicate the statements of  that all asymptotically flat solutions of massive gravity are discontinuous as $`m0`$, and phenomenologically excluded. Let us note, for completeness, that to fully vindicate this negative result one should still explore the properties of the new types of (asymptotically flat) vacua that we have discovered in section 8 above.

In addition, one should also explore the effect of using “more confining” potentials, such as the potential $`V\sigma _2\sigma _{1}^{}{}_{}{}^{2}+\lambda \sigma _{2}^{}{}_{}{}^{2}`$ considered in . There it was shown that such a mass potential can cure the runaway of the gauge variable $`e^\gamma `$ (when $`\lambda >0`$). One could hope that, within our context, it might similarly cure the $`\overline{b}`$ runaway by confining $`\overline{b}`$ within a finite interval separated from zero. However, we do not think that such a modification of the mass term is sufficient for curing the discontinuity problem. First, we anticipate that any modification of the mass term by additional “confining” terms containing only parameters of order unity (such as $`\lambda `$ in the example above) will have great difficulty in confining $`\overline{b}`$ because it will have to fight against the extremely large dimensionless number $`hGT/m^2U/(mR)^2`$ which, according to perturbation theory, will tend to enter the evolution of the gauge variables $`\overline{b}\mu h`$. \[And, indeed, our preliminary numerical simulations found a runaway of $`\overline{b}`$ toward zero with potentials $`V\sigma _2\sigma _{1}^{}{}_{}{}^{2}+\lambda \sigma _{2}^{}{}_{}{}^{2}`$.\] Second, even if one “confines by brute force” the evolution of $`\overline{b}`$ by adding sufficiently efficient terms in $`V`$, we see no reason why the resulting global solution should get close to a GR one. We anticipate that, even if the runaway singularity is cured, the resulting (non spontaneously symmetry breaking) solutions will be very different from GR solutions, and therefore phenomenologically unacceptable.

Our negative result still leaves open the possibility, proposed in , that continuity and phenomenological compatibility be restored by considering solutions that asymptote some non-trivial cosmological spacetime of curvature $`=𝒪(m^2)`$. \[Note that this proposal differs from what is suggested by the continuity results in (A)dS backgrounds  which need $`𝒪(m^2)`$ to restore continuity.\] We have provided evidence for this possibility by exhibiting a particular class of “symmetry breaking” solutions with matter sources, which generalize black-hole-type solutions found long ago in . As $`m0`$, these solutions smoothly tend to an asymptotically flat general relativistic spherical star model. The “price” for continuity is, however, twofold: (i) asymptotically the solution tends, when $`m0`$, to a de Sitter solution, and (ii) the mass term must belong to one of the special “universality classes” delineated in . The “price” (i) might actually be seen to be a virtue, as current cosmological data favor such an asymptotic state for our universe. On the other hand, the requirement (ii) is not innocent because it was shown in the latter reference that the mass terms coming from brane models do not belong to the needed “symmetry breaking” universality class.

Let us finally note that our positive results (on the re-establishment of continuity for symmetry-breaking asymptotically dS solutions) are still far from a satisfactory proof of the physical admissibility of massive gravity theories. Many subtle issues need to be checked, such as: (i) can the simple spherically symmetric solutions explored here be deformed into inhomogeneous spacetimes able to describe the solar system and the observed cosmological universe ? (ii) are such solutions attractors of the cosmological evolution ?, (iii) are they stable under quantum fluctuations101010See  for recent studies of massive spin fields in (A)dS backgrounds. ? In particular, we worry about the (classical and quantum) effect of the “sixth degree of freedom” present in massive gravity  whose rôle in the simple solutions explored so far might have been unduly downplayed.

Let us finally mention that a better setting for investigating these questions (especially the “cosmological attractor” issue) would be a full bi-gravity theory , which, contrary to the purely massive gravity explored here, does not contain anya priori fixed, non-dynamical background. In that case, there exist bi-de Sitter “locked” solutions  which are cosmological attractors (and which are more general than the above considered SSB ones). They can be trivially generalized to bi-Schwarzschild-deSitter solutions with correlated source masses. An interesting question then is whether these kinds of solutions are “spacelike” attractors as well, if we detune the special relation between the source masses.

Acknowledgments: We would like to thank B. Tekin and W. Siegel for bringing to our attention Refs.  and . I.K. is supported in part by PPARC rolling grant PPA/G/O/1998/00567 and the EC TMR grants HPRN-CT-2000-00152 and HRRN-CT-2000-00148. A.P. is supported in part by the European Community’s Human Potential Programme under contracts HPRN–CT–2000–00131 Quantum Spacetime, HPRN–CT–2000–00148 Physics Across the Present Energy Frontier and HPRN–CT–2000–00152 Supersymmetry and the Early Universe. I.K. and A.P. wish to thank IHES for hospitality during the inception of this work. T.D. dedicates this work to the memory of F. Damour.

Appendix A: The functions $`f_t`$, $`f_R`$ and $`f_g`$

In this Appendix we will present the functions $`f_t`$, $`f_R`$ and $`f_g`$ for the various massive gravity theories that we considered in the text. First, these quantities for the postulated energy-momentum $`T_{\mu \nu }^{(V)}`$ (see (65)) of  read:

$`f_t={\displaystyle \frac{1}{2}}\left[e^{\mu +\lambda }\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2+2e^\mu 3\right]`$ (A.1)
$`f_R={\displaystyle \frac{1}{2}}e^\mu \left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2(32e^\mu e^\nu )`$ (A.2)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^{\lambda \mu }}{2R}}\left[1e^\lambda \left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2\right]\left[e^{\lambda +\mu }+(e^\nu +2e^\mu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2\right]`$ (A.3)
$`+{\displaystyle \frac{1}{4}}\nu ^{}e^\nu \left[e^{\lambda \mu +\nu }(e^\nu +2e^\mu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2e^{\lambda +\mu }\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^22e^\mu +3\right]`$
$`+{\displaystyle \frac{1}{4}}e^{\lambda \mu }[R\mu ^2(e^\nu 3)2\lambda ^{}(1{\displaystyle \frac{R\mu ^{}}{2}})(e^\nu +2e^\mu 3)`$
$`\mu ^{}\left(4(e^\mu 3)+e^\nu (4+R\nu ^{})\right)`$
$`+6R\mu ^{\prime \prime }4e^\mu R\mu ^{\prime \prime }+2e^\nu (\nu ^{}R\mu ^{\prime \prime })\left]\right(1{\displaystyle \frac{R\mu ^{}}{2}})`$
For the theories which we have been studying, and for which the equations of motion stem from an action principle we obtain the following expressions.

For the potential $`𝒱^{(1)}`$ (see (5)) we have:

$`f_t={\displaystyle \frac{1}{4}}e^\nu \left[e^{\lambda +\mu }(3e^\nu +2e^\mu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2+3e^\nu (2e^\mu 3)+e^\mu (e^\mu 6)+6\right]`$ (A.4)
$`f_R={\displaystyle \frac{1}{4}}e^\lambda \left[3e^{\lambda +\mu }(2e^\mu +e^\nu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2+e^\nu (2e^\mu 3)+e^\mu (e^\mu 6)+6\right]`$ (A.5)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^\mu }{8R}}[8(e^\lambda 1)(e^\mu +e^\nu 3)+3R^2\mu ^2(4e^{\lambda +\mu }+e^\nu (e^\lambda +4)+4e^\mu 3e^\lambda 12)`$ (A.6)
$`+{\displaystyle \frac{1}{2}}R^3\mu ^3(e^\mu +e^\nu 3)(R\mu ^{}8)6R\lambda ^{}e^\lambda (2e^\mu +e^\nu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)`$
$`+R^2\mu ^{}\nu ^{}e^\lambda (2e^\mu +3e^\nu 3)+8R\mu ^{}(6+3e^\lambda 2e^\mu e^\nu (e^\lambda +2)4e^{\lambda +\mu })`$
$`6Re^{\lambda +\nu }(\nu ^{}+R\mu ^{\prime \prime })2Re^\lambda (2e^\mu 3)(\nu ^{}+3R\mu ^{\prime \prime })\left]\right(1{\displaystyle \frac{R\mu ^{}}{2}})^3`$
For the potential $`𝒱^{(2)}`$ (see (6)) we have:

$`f_t={\displaystyle \frac{1}{2}}e^{(3\nu 3\mu \lambda )/2}\left[e^{\lambda +\mu }\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^1+(2e^\mu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)\right]`$ (A.7)
$`f_R={\displaystyle \frac{1}{2}}e^{(3\lambda \mu \nu )/2}(2e^\mu +e^\nu 3)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^1`$ (A.8)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^{(\lambda +3\mu +\nu )/2}}{8}}[2\lambda ^{}e^{\lambda +\mu }(2e^\mu +e^\nu 3)(1{\displaystyle \frac{R\mu ^{}}{2}})+{\displaystyle \frac{2e^\nu }{R}}(4e^\mu +R\nu ^{}(2e^\mu 3))(1{\displaystyle \frac{R\mu ^{}}{2}})^3`$ (A.9)
$`+{\displaystyle \frac{e^{\mu +\lambda }}{R}}(8(e^\mu 3)4R\mu ^{}(e^\mu +3)+R^2\mu ^2(2e^\mu +3)2R^2\mu ^{}(2e^\mu 3))`$
$`{\displaystyle \frac{e^{\lambda +\mu +\nu }}{R}}(8+2R\nu ^{}+R\mu ^{}(R\mu ^{}R\nu ^{}4)+2R^2\mu ^{\prime \prime })`$
$`+{\displaystyle \frac{8e^\mu }{R}}(e^\mu 3)(1{\displaystyle \frac{R\mu ^{}}{2}})^3\left]\right(1{\displaystyle \frac{R\mu ^{}}{2}})^2`$
For the potential $`𝒱^{(3)}`$ (see (7)) we have:

$`f_t={\displaystyle \frac{1}{4}}e^{2\mu }\left[e^\lambda (e^\mu +e^\nu (3e^\mu 2))\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2+2e^\mu (3e^\nu +1)3e^{2\mu }(2e^\nu +1)e^\nu \right]`$ (A.10)
$`f_R={\displaystyle \frac{1}{4}}e^{\nu 2\mu +\lambda }\left[e^\nu (1+6e^\mu (e^\mu 1))e^\mu (3e^\mu 2)+e^\lambda (e^\nu (3e^\mu 2)e^\mu )\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2\right]`$ (A.11)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^{\lambda \nu 2\mu }}{8R}}\{e^\mu [8(e^\lambda 1)2(1{\displaystyle \frac{R\mu ^{}}{2}})R\lambda ^{}2R\nu ^{}+R\mu ^{}(R\nu ^{}+R\mu ^{}8)2R^2\mu ^{\prime \prime }]`$ (A.12)
$`+e^\nu [8(e^\lambda 1)(3e^\mu 1)R^2\mu ^2(3e^\mu 4)+2R\lambda ^{}(3e^\mu 2)(1{\displaystyle \frac{R\mu ^{}}{2}})`$
$`+R\mu ^{}(3e^\mu 2)(8+R\nu ^{})2R(3e^\mu 2)(\nu ^{}R\mu ^{\prime \prime })]\left\}\right(1{\displaystyle \frac{R\mu ^{}}{2}})`$
For the potential $`𝒱^{(4)}`$ (see (8)) we have:

$`f_t={\displaystyle \frac{1}{2}}e^{(3\lambda +5\mu +\nu )/2}\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)\left[e^\lambda (3e^\mu 2)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^2\right]`$ (A.13)
$`f_R={\displaystyle \frac{1}{2}}e^{(\lambda +7\mu +3\nu )/2}\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)^3\left[e^\nu (3e^\mu 2)e^\mu \right]`$ (A.14)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^{(3\nu +3\lambda +7\mu )/2}}{16R}}\{e^\mu [2(1{\displaystyle \frac{R\mu ^{}}{2}})(86R\lambda ^{}(1{\displaystyle \frac{R\mu ^{}}{2}})6R\nu ^{}6R^2\mu ^{\prime \prime }`$ (A.15)
$`+R\mu ^{}(3R\nu ^{}+5R\mu ^{}20))+4e^\lambda (R\nu ^{}(3e^\mu 2)4)]`$
$`+e^\nu [2(1{\displaystyle \frac{R\mu ^{}}{2}})(8+24e^\mu R\mu ^{}(60e^\mu 44)`$
$`+R^2\mu ^2(15e^\mu 14)6R\lambda ^{}(3e^\mu 2)\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)`$
$`6R^2\mu ^{\prime \prime }(3e^\mu 2))+16e^\lambda (3e^\mu 1)]\left\}\right(1{\displaystyle \frac{R\mu ^{}}{2}})`$
For the potential $`𝒱^{(\sigma )}`$ (see (83)) we have:

$`f_t={\displaystyle \frac{1}{8}}e^{(3\nu 3\mu \lambda )/4}[3\mu (\mu +\nu +4)+\lambda (3\mu +\nu +4)`$
$`2(2\mu +\nu +4)\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}})]\sqrt{1{\displaystyle \frac{R\mu ^{}}{2}}}`$ (A.16)
$`f_R={\displaystyle \frac{1}{8}}e^{(3\lambda 3\mu \nu )/4}[\nu (3\mu +\lambda +4)+\mu (3\mu +3\lambda +8)`$
$`2(2\mu +\nu )\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}})]\sqrt{1{\displaystyle \frac{R\mu ^{}}{2}}}`$ (A.17)
$`f_g`$ $`=`$ $`{\displaystyle \frac{e^{(\lambda +\nu +3\mu )/4}}{32\sqrt{1\frac{R\mu ^{}}{2}}}}\{(\mu ^{}+R\mu ^{\prime \prime })[\mu (8+3(\lambda +\mu )4\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}}))4(\nu +2\mu )`$ (A.18)
$`+\nu (4+\lambda +3\mu 2\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}}))]`$
$`+(1{\displaystyle \frac{R\mu ^{}}{2}})[(3\mu ^{}+\nu ^{}+\lambda ^{})(\mu (8+3(\lambda +\mu )4\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}}))`$
$`+\nu (4+\lambda +3\mu 2\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}})))`$
$`+\left(8\nu ^{}+{\displaystyle \frac{32}{R}}\right)\left(\lambda 2\mathrm{log}\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)\right)`$
$`4(3\mu ^{}(\mu +\nu )+2\nu ^{}(\nu \mu )+\lambda ^{}(\nu +3\mu ))`$
$`4\mu ^{}\left(8+3(\lambda +\mu )4\mathrm{log}\left(1{\displaystyle \frac{R\mu ^{}}{2}}\right)\right)`$
$`4\nu ^{}(4+\lambda +3\mu 2\mathrm{log}(1{\displaystyle \frac{R\mu ^{}}{2}}))]\}`$
Appendix B: The SSB solution in a different gauge

In this Appendix we present the SSB solution of section 9 in the notation of the original derivation . We now use the following gauge for the two metrics: The background flat metric is written as:

$$ds_{\mathrm{fl}}^2=dt^2+dr^2+r^2d\mathrm{\Omega }^2$$
(B.1)
while the observable metric is written in the following form:

$$ds^2=C(r)dt^2+A(r)dr^2+2D(r)dtdr+B(r)d\mathrm{\Omega }^2$$
(B.2)
The off diagonal term $`D/C`$ is equivalent (after going from the ‘flat” time $`dt`$ to the “curved” one $`dT=dt(D/C)dr`$, and changing $`dr=(dr/dR)dR`$ ) to the (gauge) variable $`\alpha `$ above.

Note the values of the eigenvalues of $`𝐟^1𝐠`$ in this new notation:

$`\mathrm{\Delta }\lambda _0\lambda _1=AC+D^2,`$ (B.3)
$`\lambda _0+\lambda _1=A+C,`$ (B.4)
$`\lambda _2=\lambda _3=B/r^2`$ (B.5)
For simplicity, we will consider a star with isotropic pressure $`P`$. The energy momentum tensor is then:

$$T_{\mu \nu }^{(\mathrm{matt})}=(\rho +P)u_\mu u_\nu +Pg_{\mu \nu }$$
(B.6)
where $`u_\mu =(\sqrt{C},\frac{D}{\sqrt{C}},0,0)`$ is the unit velocity along the timelike static Killing vector $`\xi _\mu =(C,D,0,0)`$. Explicitly, it reads:

$$T_{\mu \nu }^{(\mathrm{matt})}=\left(\begin{array}{cccc}\rho C& \rho D& & \\ \rho D& (\rho +P)\frac{\mathrm{\Delta }}{C}\rho A& & \\ & & PB& \\ & & & PB\mathrm{sin}^2\theta \end{array}\right)$$
(B.7)
The first step for finding the solution is to impose the SSB condition (149) (valid when $`\alpha 0`$, i.e. $`D0`$). For the mass term $`𝒱^{(4)}`$, (8) that we consider for concreteness, this yields:

$$\lambda _2=\lambda _3=\frac{B}{r^2}=\frac{2}{3}$$
(B.8)
Then we can follow the steps indicated above to derive the exterior solution. The exterior constancy of $`\mathrm{\Delta }`$ came from considering the effect of the source combination $`\rho ^{(g)}+P_r^{(g)}`$. In the present notation this corresponds to using $`AT_{tt}^{(g4)}+CT_{rr}^{(g4)}=0`$, to get $`AT_{tt}^{(\mathrm{matt})}+CT_{rr}^{(\mathrm{matt})}=\mathrm{\Delta }(\rho +P)`$, which implies that $`AR_{tt}+CR_{rr}=8\pi G\mathrm{\Delta }(\rho +P)`$. The latter equation explicitly reads:

$$\mathrm{\Delta }^{}=8\pi G\frac{\mathrm{\Delta }^2}{C}r(\rho +P)$$
(B.9)
The conservation of the gravitational energy-momentum tensor $`T_{\mu \nu }^{(g4)}`$ then gives:

$$3r\mathrm{\Delta }^{}9\mathrm{\Delta }^2+(6(A+C)4)\mathrm{\Delta }=0$$
(B.10)
while the conservation of matter energy-momentum tensor $`T_{\mu \nu }^{(\mathrm{matt})}`$ gives:

$$P^{}=\frac{1}{2}\frac{C^{}}{C}(\rho +P)$$
(B.11)
Combining Eqs. (B.9) and (B.10) we obtain the following expression for the combination (B.4):

$$A+C=\frac{3}{2}\mathrm{\Delta }+\frac{2}{3}4\pi G\frac{\mathrm{\Delta }}{C}r^2(\rho +P)$$
(B.12)
One finally needs only one non-trivial field equation to close the system, e.g. the (t,t) component of the Einstein equations:

$$\mathrm{\Delta }^2(632\pi Gr^2\rho )4\mathrm{\Delta }(C+rC^{})+4rC\mathrm{\Delta }^{}=3m^2r^2\sqrt{\mathrm{\Delta }}$$
(B.13)
The system of radial ODE’s (B.9), (B.11), (B.13) gives, for any given equation of state $`P=f(\rho )`$, a closed evolution system for the variables $`\mathrm{\Delta }`$, $`P`$ and $`C`$. A crucial feature of this system is that the graviton mass $`m^2`$ appears only with a positive power, and in front of a lower-derivative term (see (B.13)). General theorems on the continuity (with respect to parameters) of solutions of ODE systems then guarantees that, if $`m^2R_{}^2`$ is small enough, any regular interior solution in the GR limit will be smoothly deformed into some regular massive gravity interior one.

To complete this Appendix, let us give the explicit form, in the $`A,B,C,D`$ gauge of the exterior solution. It is given by $`\mathrm{\Delta }=const.\mathrm{\Delta }_0`$, $`A+C=\frac{3}{2}\mathrm{\Delta }_0+\frac{2}{3}`$, and by the nontrivial equation for $`C`$ whose general solution is easily found to be:

$$C=\frac{3}{2}\mathrm{\Delta }_0(1p(r))$$
(B.14)
where

$$p(r)=\frac{\stackrel{~}{R}_S}{r}+\frac{m^2}{6\mathrm{\Delta }_0^{3/2}}r^2$$
(B.15)
One then recognizes the Schwarzschild-de Sitter (rather than AdS) solution by changing the time variable to:

$$d\stackrel{~}{T}=\frac{1}{\sqrt{\beta +1}}\left(dtdr\frac{\sqrt{p(p+\beta )}}{1p}\right)$$
(B.16)
where $`\beta =(4/9\mathrm{\Delta }_0)1`$. The observable metric gets transformed to an exact Schwarzschild-de Sitter metric (here represented in a conformally scaled form):

$$ds^2=\frac{2}{3}\left[(1p)d\stackrel{~}{T}^2+\frac{1}{1p}dr^2+r^2d\mathrm{\Omega }^2\right]$$
(B.17)
The connection of these coordinates with the ones used in section 9 is simply $`R^2=\frac{2}{3}r^2`$ and $`dT=\frac{4}{9\mathrm{\Delta }_0}d\stackrel{~}{T}`$. Note also that $`\stackrel{~}{R}_S=\sqrt{\frac{3}{2}}R_S`$.