{'fraction_non_alphanumeric': 0.03902568331762488, 'fraction_numerical': 0.012046418473138549, 'mean_word_length': 4.295227698066126, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Restoring General Relativity in massive bi-gravity theory

We study static spherically symmetric solutions of massive bi-gravity theory, free from the Boulware-Deser ghost. We show the recovery of General Relativity via the Vainshtein mechanism, in the weak limit of the physical metric. We find a single polynomial equation determining the behavior of the solution for distances smaller than the inverse graviton mass. This equation is generically of the seventh order, while for a specific choice of the parameters of the theory it can be reduced to lower orders. The solution is analytic in different regimes: for distances below the Vainshtein radius (where General Relativity is recovered), and in the opposite regime, beyond the Vainshtein radius, where the solution approaches the flat metric.


I. INTRODUCTION
Modification of General Relativity (GR) giving mass to the graviton started from the work of Fierz and Pauli [1]: they considered a linear theory of a single massive spin-2 field living in flat spacetime. The first nonlinear realization of the massive graviton was presented much later [2], although in a completely different context. To extend at the nonperturbative level the action of Fierz and Pauli, adding to the Einstein-Hilbert action a non-derivative self-coupling for the metric g, it is required the introduction of an additional metric f that may be a fixed external field, or a dynamical one. When f is non-dynamical we are in the framework of aether-like theories where diff. invariance can be restored by the introduction of a suitable set of Stuckelberg fields; on the other hand, if it is dynamical, we enter in the context of bigravity theories.
Unfortunately, for a generic potential the theory has a ghost propagating degree of freedom (d.o.f.) [3], the so-called Boulware-Deser ghost, associated with the Ostrogradski ghost in more general setup. Notably, the Fierz-Pauli theory was constructed so that it has five healthy propagating degrees of freedom, while the sixth mode is removed due to the specific choice of the coefficients in the mass term. When the theory is promoted to nonlinear level or considered around non-flat background, the sixth mode reappears leading to ghost instability. This problem was solved only recently by a careful choice of the massive gravity potential [4] such that, on fully non perturbative level, the theory propagates only five degrees of freedom [5] (see also [6][7][8]). We will refer to this theory as the de Rham-Gabadadze-Tolley (dRGT) model. Of course, one has to be cautious about the rest -the five propagating degrees of freedom -however, at least there is no a priori Ostrogradski instability associated with the sixth mode. The original dRGT construction of massive gravity takes the additional metric as a flat, non-dynamical field; then it was extended in the bi-gravity context supplying an extra Einstein-Hilbert term for the second metric [9]. The bi-metric approach to massive gravity is the main subject of our work. We would like to emphasize here that the bi-gravity formulation of massive gravity is not just a theoretical entertainment, but also cosmology calls for it. When the second metric is non-dynamical and Minkowski there is no homogeneous spatially flat FRW solution [10,11], on the contrary in the bigravity formulation flat FRW homogeneous solutions do exist [12][13][14]. Moreover, the cosmological perturbations are far less problematic [15]: all the d.o.f. propagate at the linear level without ghost instabilities [16]. For more recent works see [17].
Another problem which arises in massive gravity models is generic for theories with extra propagating degrees of freedom. Since the graviton mass turns on (at least) three extra degrees of freedom, it is expected that the extra interaction change the Newtonian limit or/and the light deflection. This can be easily seen in the so-called decoupling limit -the scalar part of the graviton is directly coupled to the matter with approximately the same coupling constant of the helicity-2 piece. The extra scalar behaves similar to the Brans-Dicke field ruling out the theory on observational ground. Moreover, a naive way to recover GR sending the mass of the graviton to zero, does not solve the problemthe so called vDVZ discontinuity [18] -since the theory with arbitrary small (but non-zero) graviton mass contains the extra propagating scalar, absent in the massless theory (GR). A way to overcome this difficulty was proposed by Vainshtein [19] in 1972. Vainshtein noticed that the linear approximation breaks down at some distance far from the source (now called the Vainshtein radius) and therefore one cannot approximate the solution by linearizing it close to the source. On the contrary, he showed that a solution can be found by expanding around the GR solution in powers of the graviton mass. This construction indicates the possibility that the extra propagating scalar mode can be hidden close to the source by non-linear effects. The question still remained if this close-to-GR solution could have been matched to an asymptotically flat solution [3]; it was in fact argued that it was not possible [20,21]. Only recently it has been realized that the Vainshtein solution close to the source matches the one obtained by linearization far from the source. Therefore GR is restored locally for asymptotically flat solutions, at least for some potentials in massive gravity [22,23] 1 (see also [26] for a more recent work). This matching was shown for potentials giving rise to the sixth dangerous mode. Later it was found, both analytically in the so-called decoupling limit [27][28][29] and numerically [30], that GR is also restored in the dRGT model. In the framework of the same model, the Vainshtein mechanism was studied in the decoupling limit for asymptotically non-flat spacetimes in [31], where these solutions were shown to be the only stable ones (see also a related work on Galileons [32] and on Horndeski theory [33]); and in the quasi-dilaton extension in [34]. For a recent review on the Vainshtein mechanism see [35]. 2 Other possible issues of massive gravity we are not going to discuss here include superluminality [37] (see, however, discussion in [38,39] on relation between causality and superluminality), strong coupling problem [40,41], instability of black holes [42].
For what concerns the bi-gravity formulation of the dRGT model, the numerical study of the Vainshtein mechanism was presented in [30], while the far-distance analytic expansion valid outside the Vainshtein radius was found in [43] and then studied up to the second order in [44]. Some estimates of the Vainshtein suppression have been put forward also in [45] in order to calculate the emission of gravitations waves. However, there is still a lack of an analytic analysis of the Vainshtein mechanism, that we fill in this work.
In this paper we analytically study the spherically symmetric solutions in the bi-gravity extension of the dRGT model and show that the Vainshtein mechanism indeed works also with a dynamical second metric. We make our analysis in the approximation of the weak gravitational field of the physical metric -the one coupled to the matter -which is always correct for (nonrelativistic) weak matter sources. This assumption allows us to treat the problem mostly analytically. The key equation we obtain is an algebraic polynomial equation, generically of the 7th order, for a function of the radius entering, in our ansatz, in the second metric. The other functions that parameterize the metrics are expressed in terms of this key function. In two regimes -inside and outside the Vainshtein radius -we can solve (approximately) the algebraic equation that gives different branches of the solution, then we identify the one which ensures the flat asymptotic behavior. We exhibit a solution featuring the GR behavior for the physical metric inside the Vainshtein radius and that matches the asymptotically flat (Yukawa decaying) solution of the linearized equations.
The paper is organized as follows. In Sec. II for completeness, we reanalyze the Vainshtein mechanism for the original dRGT model (with fixed reference metric) in a slightly different manner than in [27][28][29]. This approach will be generalized in the main part of the paper (Sec. III) for the bi-gravity extension, where we rigorously study the Vainshtein mechanism with the second metric dynamic. Our conclusions are formulated in Sec. IV.

II. DRGT MODEL
In this section we obtain static spherically symmetric solutions in the limit of weak gravitational field for the dRGT model. We reproduce the results already found in [27][28][29] for distances smaller than the Compton wavelength of the graviton, where the decoupling limit (DL) is a good approximation. In our approach however, originally introduced in [23] and called "weak-field approximation", the DL scheme is not used. The weak-field approximation allows to capture both the DL and the Yukawa part of the solution outside the Compton wavelength, where the DL ceases to operate. 3 In Sec. III, this scheme -with appropriate modifications -will be applied to the bi-gravity extension of the model.

A. Action and equations of motion
The action for the dRGT model can be written as follows [4] It is convenient to express the interaction potential U[g, f ] in terms of the matrix K, such that K µ ν = δ µ ν − γ µ ν , where the matrix γ µ ν is the square root of the product of the inverse physical metric g µα and the fiducial metric f αν , i.e. γ µ ν = √ g µα f αν , in the sense that γ 2 µ ν = γ µ α γ α ν = g µα f αν . As it is often assumed we will consider the fiducial metric to be flat 4 . The potential U consists of three pieces, each of them, in terms of K, reads where we introduced the notations [K] ≡ tr(K) =K ρ ρ and [K n ] ≡ tr(K n ) = (K n ) ρ ρ . Varying the action with respect to g µν one obtains where G µν is the Einstein tensor and on the right hand side there are the contributions from the energy-momentum tensor for the matter, T The last can be computed and gives

B. Static Spherically Symmetric Solutions
In this section we study spherically symmetric solutions for the case where the non-dynamical second metric f µν parametrizes a flat Minkowski space-time. The study of the Vainshtein mechanism in this case has been already done in a number of papers [28,29,46], in the decoupling limit. Here we reproduce these results, moreover we will give some additional new upshots outside the DL regime. The procedure is to consider the full equations of motion and then make reasonable approximations valid for the regimes in which we are interested. Following the weak-limit approximation scheme [23], we take the following ansatz: This ansatz is not the most general, indeed we do not consider the case with one of the two metrics off-diagonal, but (4) is where we find the Vainshtein mechanism at work. Since we are interested in a recovery of GR solutions, we require for weak matter sources to have weak gravity, i.e. that the functions ν and λ are small as well as their derivatives. So the first step is to consider: and to retain all the non-linearities in µ and µ . The tt, rr and θθ components of the Einstein equations in this approximation read, where we introduced Note that the pressure in the r.h.s. of (7) and (8) disappears as a consequence of the conservation of the matter energy momentum tensor in the weak field regime (5). Of course, like in GR, the three Einstein equations and the Bianchi one are not all independent, so we can consider (6), (7) and (9) as our independent set to be solved. From this set we are able to obtain one second order ODE on µ only. Indeed we can solve (7) for λ and then, substituting into (6) and (9), we end up with two equations, one for ν , ν, µ, and the other for ν , ν , ν, µ , µ. Taking the first equation and its first and second derivative, together with the second equation and its first derivative, we have a system of five equations in ν , ν , ν , ν, µ , µ , µ. We can then solve algebraically four of them for ν and all its three derivatives and, substituting in the last equation, it will be a second order differential equation only on µ of the form A, B, C, D, E are functions of µ whose form is not particularly illuminating, so we can omit it here. It is worth to stress the difference with respect to the equivalent equation found in the ghosty massive gravity theories [23]. Our second order equation needs two initial conditions in order to be solved, fixing to one the number of degree of freedom that it describes. For other kind of potentials [40] that exhibits the BD instability, the equation is of the fourth order [23], meaning four initial conditions and therefore two d.o.f. One of these modes is absent for the dRGT potential. Eq. (10) is hard to solve, even numerically. Without solving it, however, the equation clearly indicates that the weak field approximation (5) gives the relevant features of the fields for all ranges of distance for non-relativistic sources. In order to understand the behavior of the solutions, below we will consider various regimes.

Linear regime
Since we require asymptotically flat solutions, we expect that far away from the source also the field µ becomes small. Therefore, assuming µ 1 in (6), (7) and (9), as well as its derivative, we get the solutions: where C is an integration constant that in the following we will see equal to the Schwarzschild radius r S . Clearly the gravitational potentials (12) and (13) exhibit the vDVZ discontinuity in the limit of m → 0 , instead the field µ shows a singularity in the same limit. Actually the linear regime is nowhere allowed in the vanishing mass limit, indeed this regime exists only for values of radii for which µ 1. For non-zero small m, as can be easily seen from (11), the linear regime is valid for r r V , where r V is the Vainshtein radius, In the limit m → 0 we have r V → ∞ , making unreliable the condition to be outside r V . Hence, to look at the small mass limit, we need to consider a non-linear regime in µ. Finally, it is remarkable that the weak field approximation is able to retain the asymptotically Yukawa decay at large distances from the source, as well the Vainshtein crossover as we will see in the next paragraph.

Inside the Compton wavelength
In order to study the behavior of the solutions for which the m → 0 limit is well defined and the Vainshtein mechanism operates, we need to consider the distances inside the Compton wavelength, i.e. r 1/m. In this regime we can neglect terms ∼ m 2 λ and ∼ m 2 ν in the r.h.s. of equations (6)- (8). In this approximation, we can integrate equation (6) to obtain where R is the radius of the source, the Schwarzschild radius r S reads and the integration constant in (14) has been chosen to insure the continuity of the solution at the surface of the star. In the following for simplicity we will consider only a constant density source. Neglecting ∼ m 2 λ and ∼ m 2 ν in Eq. (7) one obtains while the integration of (8) gives (15) up to an integration constant. From (15) and (14) we find Finally, combining (9), (14) and (16) we get a single algebraic equation on µ, The last equation (17) corresponds to the one found in the DL of the model [28,29,46]. This confirms that the approximations we made here correspond to the DL in the full equations of motion. All the physics is hence enclosed in equation (17): once we have its solutions we can determine the gravitational potentials through equations (14) and (16). Eq. (17) is a fifth order algebraic equation and its solutions can not be presented in a closed form. Choosing the parameters of the theory α and β, we can lower the degree of such equation in order to get analytically solvable ones. E.g., for the minimal massive gravity potential with only U 2 -it corresponds to β = 0 and α = 1 -Eq. (17) becomes a third order algebraic equation. This special case was first studied in [27], where it was shown that the Vainshtein mechanism works properly reproducing GR inside the Vainshtein radius.
Since the solutions of the set of equations (17), (14) and (16) have been largely studied [28,29,46], we only report schematically the behavior of the solutions for different sub-regimes. To consider the most interesting case, for which the Vainshtein mechanism takes place, we set β > 0 so that Eq. (17) has two complex and three real solutions. Only one of the three real branches of the solution recovers GR and is asymptotically flat, so we give it in Table (18).
Inside the Vainshtein radius we find that µ = −1/ √ β at the leading order and the gravitational potentials are of the GR form plus small corrections. Clearly, the vanishing mass limit is well defined inside the Vainshtein radius and the corrections to GR smoothly vanish as m → 0. Outside the Vainshtein radius we find the asymptotically flat weak field solution where also µ 1: this solution match the one obtained in the linear regime (11)-(13) provided C = r S and r 1/m. To get some understanding of the solutions from Eq. (17), note the r.h.s of itself. The ratio (r V /r) 3 that appears there becomes large inside the Vainshtein radius, and small otherwise. For r r V we have either |µ| 1 retaining so the higher power of µ in the equation, or the leading order of µ cancels the r.h.s. itself, i.e. µ = ±1/ √ β. For r r V instead we can neglect at the first order the r.h.s. of (17), obtaining therefore three real constant values for µ: obviously only µ = 0 gives asymptotically the flat metric.
The branch of the solution presented in (18) is unphysical for β < 0 due to the square root in µ. A complete description of the solutions for the whole range of the free parameters α and β can be found in [46].

III. BI-GRAVITY
This section is devoted to the study of the Vainshtein mechanism in the Hassan-Rosen bi-gravity extension of the dRGT model [9]. In the bi-gravity approach the second metric f µν , which was fixed before, now becomes dynamical. The weak-field approximation, that we applied in the previous section, is also useful here. The trick is to deform the ansatz (4) to include the dynamics of the metric f µν . After writing down the action and the field equations in Sec. III A, we introduce the ansatz and identify the functions that can be treated linearly in the limit of weak sources, together with a fully non-perturbative function. In this way we are able to obtain again one algebraic polynomial equation, which captures the distances inside the Compton wavelength of the graviton.

A. Action and equations of motion
We consider the dRGT model where the additional metric f µν is dynamical thanks to its own Einstein-Hilbert term in the action [9], The interaction potential in (19) is given by the same expressions in (2) and (3) 5 . Note an extra parameter κ in the action (19), which accounts for a possible difference in Planck masses for the two gravitational sectors. One can realize that the limit κ → ∞ corresponds to the freezing of the metric f , therefore recovering the model with flat fiducial metric. Varying the action (19) with respect to g µν , we obtain 5 Note that in principle we can add terms U 0 = 1 and U 1 = [K]. These terms, however, account for cosmological terms for the g and f metric. Since we aim to find asymptotically flat solutions, we exclude those terms.
where G µν is the Einstein tensor associated with g µν and T µν reads, On the other hand, the variation of the action with respect to f µν gives, One can observe a useful relation when working with up-down indices, namely, for T µ ν ≡ g µα T αν and T µ ν ≡ f µα T αν , we have T µ ν = −T µ ν + Uδ µ ν .

B. Static Spherically Symmetric Solutions
Continuing the idea of the weak-field approximation scheme, we consider the following parametrization for the two metrics ds 2 = −e ν dt 2 + e λ dr 2 + r 2 dΩ 2 , df 2 = −e n dt 2 + e l (r + rµ) where ν, λ, n, l and µ are the r-dependent functions that describe the spherically symmetric foliation of the space-time in a common coordinate system. Again, this ansatz is not the most general one, since we are considering bi-diagonal metrics. Let us note that our ansatz (20) is compatible with the bi-flat solution obtained when one imposes T = T = 0, and that will be the reference point for the asymptotic behavior.

Static spherically symmetric ansatz and equations of motion.
Following the procedure used for the dRGT model, we consider a non-relativistic matter source with constant density. This allow us to assume ν, λ, n and l small, as well as their derivatives, {λ, ν, l, n} 1, {r λ , r ν , r l , r n } 1 , and to retain all the non-linearities in the field µ and µ . The tt, rr and θθ components of the Einstein equations in this approximation are for the first metric and for the second metric − l 2r + 1 + µ 2 (r + r µ) n + 1 + (r + r µ) (r + r µ) Of course, these equations are not all independent, indeed a combination of (25) and (28) can be obtained taking a suitable combination of (24) and (27) and its derivative. The Bianchi identities, ∇ Note that assuming l = n = 0 in (29) we get back (9), the Bianchi identity for the model with one dynamical metric. It is worth to mention that we were not able to obtain an analogue of equation (10), because, applying the similar approach that we described there, we find a system of linear equations which is not linearly independent. Therefore, in order to find analytical solutions to this set of equations, we need to do other approximations. This means to look at more specific regimes inside the weak field one (22).

Linear regime
Since asymptotically we want to find the bi-flat solutions (21), we expect that far away from the source also the field µ becomes small. Hence, assuming µ 1 and rµ 1 in (23), (24), (26), (27) and (29), we find the solutions: where C 1 and C 2 are two integration constants that will be determined in the next paragraph to be It is important to stress that (once (35) is taken into account) taking the limit κ → ∞, which freezes the dynamics of the second metric, we recover the solutions found in the same regime of the original dRGT model, i.e. (11)-(13) and l = n = 0 . For the m → 0 limit, the vDVZ discontinuity appears with the divergence in µ, the same arguments as in the previous section apply here: the linear regime is nowhere allowed in the vanishing mass limit. Again, to properly consider this limit, we need to rely on nonlinearities in the field µ .
Again, all the information is retained in an algebraic equation for µ only: for bi-gravity this equation is of the seventh order compared to the fifth order equation for the dRGT model. Generically equation (40) has seven solutions that can be real or complex.
For β > 0, it has three real and two complex conjugates solutions for all values of r; the remaining two solutions are real inside some radius r * and join together in a complex conjugates pair for r > r * , see Fig. 1. The value of r * depends on the choice of the free parameters α and β. The three everywhere real solutions for µ, shown in Fig. 1, have different asymptotic behavior, however all three recover GR inside the Vainshtein radius. One of these three is asymptotically flat (Vainshtein-Yukawa solution), and the others (dashed) have non-flat asymptotics (solutions three and four in Appendix B). One of these last two solutions, (B4), may be of interest in the context of cosmology. The asymptotically non-flat solution and a possible match to a cosmological one, however, deserves a separate study and will be discussed elsewhere.
For β < 0 two of the three solutions which are real for β > 0 become complex conjugates, therefore only one everywhere real solution is left. The other solutions show the same behavior as in the case β > 0. For β < 0 the  only everywhere real solution is given by the solution three of Appendix B and it does not show the asymptotically flat behavior nor the expected weak field solution inside the Vainshtein radius. Therefore, in the following, we will consider only the case β > 0. The special case β = 0, that for α = 1 corresponds to the minimal massive gravity potential with only U 2 , will be presented in Appendix A. However we do not guarantee that solutions exist for the whole range of parameters space; i.e. we do not exclude that for some range of the parameters α and β the asymptotic solution might not match the solution inside the Vainshtein radius. This issue would require a complete analysis for each value of α and β that is beyond the scope of this work. Indeed, contrary to Eq. (17) of the dRGT model, where a symmetry allows to study easily the whole range of parameters space [46], in Eq. (40) we were not able to find a similar strategy that facilitates the scan of solutions for general α and β.
Again, one can understand the behavior of the solutions analyzing two regimes: well inside and outside the Vainshtein radius. For r r V the ratio (r V /r) 3 in the r.h.s. of (40) becomes large leaving us with two possibilities: either |µ| 1 in order to compensate the large r.h.s. (so at the leading order the higher powers of µ dominate); or the leading order of µ cancels out the r.h.s., this happens for µ = −1 (double root) and µ = ±1/ √ β. For r r V , the ratio (r V /r) 3 is small and it suggests to neglect at first order the r.h.s of (40): this determines the three real constant asymptotical values of µ. One of these values is obviously zero and it is the one approached by the linear solution (30).
In Table ( 41) we present the branch of the solution which features the Vainshtein recovery of GR and the Yukawa decay. The other branches are schematically presented in Appendix B. The way how these regimes match together is understood from the numerical study that will be presented immediately after.
Inside the Vainshtein radius, for the metric g the GR solution is recovered (plus small corrections) and for the second metric the potentials l and n are of the same order of the corrections to GR of the metric g. Outside the Vainshtein radius we find the asymptotically flat solution with µ 1: this solution matches the one obtained in the linear regime (30)-(34) for r 1/m and as we already anticipated in (35). Note that inside the Vainshtein radius the vanishing mass limit, m → 0, is well defined and gives exactly GR for the first metric while zero for l and n potentials of the second metric. It is also important to stress that in the limit κ → ∞ for which f µν is frozen, this branch of the solution reproduces the Vainshtein-Yukawa one of the dRGT model, namely Eq. (18). In addition to the Fig. 1, in Fig. 2, 3 and 4 we depict correspondingly the potentials of the physical metric outside the source, the potentials of the second metric outside the source and the potentials for both metrics inside the source, only for the Vainshtein-Yukawa solution.  fig. (a) and, in fig. (b), is given their ratio for a better comparison. The analytic behavior for the regime well inside and outside the Vainshtein radius is given in Tab. (41).

IV. CONCLUSIONS
In this paper we studied the Vainshtein mechanism in the massive bi-gravity model with no Boulware-Deser ghost.
To attack the problem, we applied the "weak-field approximation scheme" where the metric coefficients are separated into two parts. One contains functions of the radial coordinate which remain small (i.e. the quadratic and higher order terms are negligible in comparison to the linear ones) for non-relativistic sources, even inside the Vainshtein radius; the other part is fully non-linear, with nonlinearities crucial for the existence of a GR like solution. This approach allows to capture all the important features of the solutions: the Vainshtein regime, the linear regime and the Yukawa decay.
In Sec. II we demonstrated how this scheme works for the original dRGT model -where the auxiliary metric is fixed to Minkowski. Inside the Compton wavelength, our results are in agreement with previous studies in the decoupling limit [27][28][29], see in particular (18). On the other hand, we can also describe the solution outside the Compton wavelength, which is beyond the validity of the decoupling limit. Moreover, for the function µ(r) which is introduced as the non-linear piece of the metric (4), we derived a single ordinary differential equation of the second order (10) valid for all radii.
For bi-metric massive gravity we modified the approach to incorporate the dynamics of the second metric. Notably, the ansatz we introduced (20), again separates the metric coefficients in the linearizable part and in the fully non-linear part µ(r). If we additionally assume that far from the source also µ is in the linear regime, then we readily obtain the linearized solution that shows the Yukawa decay. Other regimes can be obtained assuming radii smaller than the Compton wavelength, i.e. r 1/m. In this case, it is possible to derive one algebraic equation of the seventh order on µ, Eq. (40), while the other metric functions are given in terms of it. Using this master equation, we analyzed the behavior of the solutions in various sub-regimes and identified several branches of the solution. For β > 0 the only solution that has the desired behavior -asymptotic flatness -is presented in the Table (41): this solution shows the recovery of GR inside the Vainshtein radius.
It is worth to make a comment about the asymptotically non-flat solution (B4). Although we did not study (B4) in detail in the present work, since we concentrated on asymptotically flat solutions, this solution may have a physical meaning if matched to a cosmological solution at r → ∞. The same comment applies to the choice of the potential α 3 = α 4 = 0, which we discuss in Appendix A. While for the dRGT model this simplest potential gives an asymptotically flat solution recovering GR inside the Vainshtein radius, in the case of the bi-metric massive gravity, such an asymptotic solution does not exist and the solution featuring the Vainshtein behavior becomes asymptotically non-flat. The behavior and the physical meaning of these solutions, together with a possible match to cosmological ones, deserves a separate study.
To summarize, in the bi-gravity formulation of the dRGT massive gravity, with matter coupled to only one (physical) metric, we have found the recovery of GR for the physical metric inside the Vainshtein radius and the Yukawa decay outside. At the same time, the second metric is nontrivial because of the indirect coupling to matter via the interaction (mass) term; its deviation from flat spacetime is highly suppressed and it reaches non-negligible values only around the Vainshtein radius.   fig. (a) plot of the absolute value of the all five solutions of the function µ vs. r/rV , for R < r < 1/m. There is only one everywhere real solution. In fig. (b) plot of the functions µ, λ and l vs. r/rV , for R < r < 1/m. As shown analytically in Tab (A1), the asymptotic behavior is not flat. We take for both the plots the following values for the parameters: m · rV = 10 −2 and κ = 1.
Inside the Vainshtein radius we find that |µ| 1 and, in contrast to the general case, it is this behavior of µ that gives a well working Vainshtein mechanism with GR like solutions for the metric g and small potentials l and n for the other metric. For the dRGT model with the simplest potential a similar recovery of GR was found in [27], with flat asymptotic. For the bi-gravity model we find that, for r r V , µ asymptotically approaches a non-zero constant value c 0 , giving non-decaying tails for all the other gravitational potentials. This solution does not match the asymptotically flat linear solution (30)- (34). It does not mean, however, that the solution (A1) is non-physical, it may in fact match a nontrivial cosmological solution for large r. This possibility, however, deserves a separate study and goes beyond the scope of this paper. In Fig. 5 we show also the numerical study of this solution for all the range of distances inside the Compton wavelength.
The fact that the choice of parameters β = 0 does not allow for an asymptotically flat everywhere real solution, while with the non-dynamical second metric such a solution exists, might seem surprising. There is however a simple explanation for this effect: when κ → ∞ the outer part of the (nonphysical) asymptotically flat branch (dashed † cν is an integration constant fixed by the matching condition to be: cν = − 3 r S 2 R + m 4 r S  curve in Fig. 6) and the inner part of the everywhere real solution (thick curve in Fig. 6) join together to make one asymptotically flat solution recovering GR for small radii. This can be easily seen from Fig. 6, where the absolute value of the solutions of µ are plotted for bigger values of κ.
It is also instructive to compare our findings with the numerical work of Volkov [30], where he obtains an asymptotically flat solution for the simplest potential (α = 1, β = 0). It looks as our results contradict the ones in [30]. In fact, a possible explanation lies in some specific choice of m, r S and R : it seems that in [30] these parameters are chosen such that the size of the source R is larger than r * -the point below which the asymptotically flat solution (dashed curve in Fig. 6) becomes complex -therefore avoiding the problem since inside the source µ is constant. If the source is made more compact though still non-relativistic, we expect that this problem of complex-valued solution comes back, although without complete numerical analysis of the full equations of motion we cannot prove this statement.

Appendix B: Other branches in massive bigravity
In this appendix we report the analytical solutions, for the regime inside and outside the Vainshtein radius, of the four branches of the massive bigravity solution that we omitted in the main text.