{'fraction_non_alphanumeric': 0.04646231850656833, 'fraction_numerical': 0.0179918568026427, 'mean_word_length': 4.032941540364986, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 36, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymptotically flat hairy black holes in massive bigravity

We study asymptotically flat black holes with massive graviton hair within the ghost-free bigravity theory. There have been contradictory statements in the literature about their existence -- such solutions were reported some time ago, but later a different group claimed the Schwarzschild solution to be the only asymptotically flat black hole in the theory. As a result, the controversy emerged. We have analysed the issue ourselves and have been able to construct such solutions within a carefully designed numerical scheme. We find that for given parameter values there can be one or two asymptotically flat hairy black holes in addition to the Schwarzschild solution. We analyze their perturbative stability and find that they can be stable or unstable, depending on the parameter values. The stable hairy black holes that would be physically relevant can be neither very small nor very large and always have the mass and size close to those for the ordinary black holes of mass ~10^6 solar masses. One of their two geometries is very close to Schwarzschild. If the bigravity theory indeed describes physics, this would imply that the supermassive astrophysical black holes hide inside the"hairy features"that should manifest themselves in violent processes like black hole collisions.


I. Introduction
Theories with massive gravitons provide a natural modification of the General Relativity (GR) in the infrared regime and can be used to explain the current acceleration of our Universe [1,2]. Such theories have a long history pioneered by the work of by Fierz and Pauli [3] and marked by subsequent discoveries of many interesting features, such as the vDVZ discontinuity [4,5], the Vainshtein mechanism [6], the Boulware-Deser ghost [7], culminating in the discovery of the ghost-free massive gravity [8] and ghost-free bigravity [9] theories.
The ghost-free bigravity theory is the most interesting physically. It contains two dynamical metrics, usually called g µν and f µν , describing together two gravitons, one of which is massive and the other is massless. The theory admits self-accelerating cosmological solutions [10][11][12] whose properties agree with the observations [13][14][15][16][17], with the Λ-term mimicked by the graviton mass. The theory also admits solutions describing black holes [18], wormholes [19], and other interesting solutions (see [20] for a review). In what follows we shall be discussing black holes.
The bigravity black holes can be either "bald" of "hairy". The bald black holes are described by the known GR metrics. Such solutions were first discovered long ago [21][22][23] within the old bigravity theory inspired by physics of strong interactions [24]. In the simplest case, their two metrics are both Schwarzschild-(anti)-de Sitter and can be conveniently represented in the Eddington-Finkelstein coordinates as [25,26] g µν dx µ dx ν = −Σ g dv 2 + 2dvdr + r 2 dΩ 2 , f µν dx µ dx ν = C 2 −Σ f dv 2 + 2dvdr + r 2 dΩ 2 , (1.1) with Σ g = 1 − 2M g /r + Λ g /(3r 2 ) and Σ f = 1 − 2M f /r + Λ f /(3r 2 ), where values of constants C, Λ g , Λ f are fixed by the field equations. Passing to the Schwarzschild coordinates, one can diagonalize one of the two metrics, but not both of them simultaneously. Such solutions have been much studied [27], they exist also within the ghost-free bigravity [28], and they admit the charged [29] and spinning [30] generalizations. These solutions also admit the massive gravity limit where M f = Λ f = 0 hence the f-metric is flat while the g-metric remains non-trivial, and this yields all possible static black holes in the ghost-free massive gravity theory (it seems there can be also time-dependent black holes in this theory [31]).
Next, it was noticed in [18] that if the parameters of the potential are suitably adjusted, then the ghost-free bigravity reduces to the vacuum GR in the sector where the two metrics coincide, g µν = f µν . Therefore, all vacuum GR black holes, as for example the Schwarzschild solution, g µν dx µ dx ν = f µν dx µ dx ν = − 1 − 2M r dt 2 + dr 2 1 − 2M /r + r 2 dΩ 2 , (1.2) or its spinning generalization can be imbedded into the ghost-free bigravity. A Λ-term can be included by assuming the two metrics to be proportional to each other [18,20,26]. Such solutions are different from solutions of type (1.1), for example they do not admit the massive gravity limit. In addition, the solution (1.1) is linearly stable [32], whereas the solution (1.2) is unstable with respect to fluctuations which do not respect the condition g µν = f µν [33].
These facts essentially exhaust the available knowledge of the "bald" black holes in the bigravity theory. At the same time, more general "hairy" black holes not described by the classical GR metrics can exist as well. The first example of hairy black holes in physics was found long ago [34], followed by many other examples (see [35,36] for a review), so that nowadays hairy black holes are considered as something usual. One may therefore wonder if they exist in the ghost-free bigravity theory as well.
A systematic analysis of hairy black holes in the ghost-free massive bigravity has been carried out for the first time by one of the authors [18], but none of the solutions found were asymptotically flat. In that analysis both metrics were assumed to be static and spherically symmetric. If they are not simultaneously diagonal, then the most general solution is given by (1.1). If they are simultaneously diagonal, then one of the solutions is given by (1.2), but other more general black hole solutions exist as well.
Such solutions possess an event horizon -a hypersurface that is null simultaneously with respect to g µν and f µν . Therefore, both metrics share the horizon [37,38], but its radius r g H measured by g µν can be different from the radius r f H measured by f µν . One can set r H ≡ r g H to unit value via rescaling the system (rescaling at the same time the graviton mass), but the ratio u = r g H /r f H is scale-invariant. Choosing a value of u completely determines the boundary conditions at the horizon, which allows one to integrate the equations starting from the horizon toward large values of the radial coordinate r. As a result, the set of all black hole solutions can be labeled by just one parameter u, and integrating the equations for different values of u gives all possible black holes.
Choosing u = 1 yields the Schwarzschild solution (1.2). For u = 1 one finds more general black holes supporting a massive graviton "hair" outside the horizon, but in the asymptotic region their two geometries do not become flat [18]. The latter property is generic, and trying different values of u always gives either solutions with a curvature singularity somewhere outside the horizon, or solutions which exist for all values of r but show non-flat asymptotics.
At the same time, these facts do not completely exclude a possibility of some other asymp-totically flat black hole solutions different from (1.2) which would correspond to some special values of u different from u = 1. However, even if they exist, one does not find them by a brute force via trying many different values of u, and the reason is the following. The field equations reduce to three coupled first order ODE's [18], whose local at large r solution which approaches flat geometry as r → ∞ has schematically the following structure (A, B, C being integration constants): A r + Be −r + Ce +r .
(1. 3) Here r = mr is the dimensionless radial coordinate, with m and r being the graviton mass and dimensionful radial coordinate (we assume the graviton mass to have the dimension of inverse length, so that this is rather the inverse Compton wavelength mc/ ). The Newtonian mode A/r in (1.3) arises due to the massless graviton present in the theory, while the decaying mode Be −r and the growing mode Ce +r are due to the massive graviton. Now, when integrating from the horizon, the growing mode Ce +r will be inevitably present in the numerical solution at large r and will drive the solution away from flat space. This is why one does not find asymptotically flat solutions in this way.
To get them, one should suppress the growing mode by setting C = 0, hence the local solutions at large r will comprise a two-parameter set labeled by A and B. The next step is to numerically extend this local solution toward small r, extending at the same time the local solution at the horizon labeled by u toward large r, until the two solutions meet at some intermediate point. For these solutions to agree, three (the number of the ODE's) matching conditions should be satisfied via adjusting the three parameters u, A, B. In practice, this can be done within the numerical multiple-shooting method [39]. Once u, A, B are adjusted, this yields global asymptotically flat solutions.
The difficulty, however, is that the numerical scheme requires some input values of u, A, B which should be close to the "true values", otherwise the iterations do not converge. It was apriory unclear how to choose these input values, whereas choosing them randomly does not give the convergence. Some additional information was needed to properly choose these values, but at the time of writing the article [18] such an information was not available. As a result, the conclusion of that work was that asymptotically flat hairy black holes may exist, but they should be parametrically isolated form the Schwarzschild solution (1.2).
It is interesting that adding an extra matter source to obtain not a black hole but a regular object like a star, asymptotically flat solutions can be easily constructed, as was shown first in [18] and later in [40,41]. The black hole case is more difficult.
Fortunately, the additional information was later obtained within the analysis of pertur-bations of the Schwarzschild solution (1.2) [33,42]. Denoting g S µν the Schwarzschild metric, the two perturbed metrics are g µν = g S µν + δg µν and f µν = g S µν + δf µν . Linearizing the field equations with respect to δg µν and δf µν , one finds that perturbations grow in time and hence the background Schwarzschild black hole is unstable if r H ≡ mr H ≤ 0.86. On the other hand, for r H > 0.86 the perturbations are bounded in time so that the background is stable [33].
Curiously, the mathematical structure of the perturbation equations is identical [33] to that previously discovered by Gregory and Laflamm (GL) in their analysis of black strings in D=5 GR [43]. We shall therefore refer to the Schwarzschild solution with r H = 0.86 as GL point.
This change of stability at the GL point suggests that for r H close to 0.86 there could be two different asymptotically flat solutions: the Schwarzschild solution (1.2) and also some other solution which can be approximated by the zero perturbation mode that exists at the GL point. This new solution is different from Schwarzschild although close to it, hence it describes an asymptotically flat hairy black hole. To get this solution within the numerical scheme outlined above, one should choose the input parameters u, A, B to be close the GL point, u ≈ 1, r H ≈ 0.86, A ≈ −r H /2, B ≈ 0, and it is this essential piece of information that was missing when writing Ref. [18]. As soon as the solution is obtained, one can change the value of r H iteratively, thus obtaining "fully fledged" hairy black holes which may deviate considerably from the parent Schwarzschild solution.
Remarkably, this program was accomplished by the Portuguese group [44] via explicitly constructing asymptotically flat hairy black holes in the theory in the region below the GL point, for r H < 0.86. However, some time later spherically symmetric bigravity solutions were analyzed by the Swedish group [45], and it was claimed that the Schwarzschild solution (1.2) represents the unique asymptotically flat black hole in the theory. As a result, a controversy emerged and it was unclear if asymptotically flat hairy black holes exist or not.
We have therefore reconsidered the issue ourselves and below are our results. In brief, we were able to construct asymptotically flat hairy black holes in the theory, thereby confirming the finding of [44]. We apply a very carefully designed numerical scheme to exclude any ambiguities and to take into account the arguments of [45]. In fact, these arguments correctly point to some drawbacks of numerical analysis commonly present in many publications. From the mathematical viewpoint, one has to solve a non-linear boundary value problem where the boundaries are singular points of the differential equations (horizon and infinity). Since it is difficult to approach such points numerically, various approximations are used in practice, which may give reasonable results in some cases but inevitably increases the error and leads to a numerical instability. Only very rarely does one find in the literature a correct treatment of the problem (apart from the relaxation approach), as for example in [46][47][48]. We therefore pay a special attention to the details of our numerical scheme and describe them in a very explicit way. As a result, from the methodological viewpoint, our paper gives an example of how one should properly tackle a non-linear boundary value problem with singular endpoints.
We cross-check our results with two different numerical codes written independently by two of us. Our results strongly suggest that the hairy solutions exist and are indeed asymptotically flat and regular. We discover many of their new features, in particular we obtain hairy black holes also above the GL point for r H > 0.86, and we study for the first time the perturbative stability of the solutions. We were able to identify regions in the parameter space which correspond to stable solutions, and we determined subsets of these regions which agree with the constraints imposed by the cosmological observations. We find that the viable hairy black holes should be described by the g-metric that is very close to Schwarzschild, but their f-metric is different. The mass and size of these black holes should be close to those for the usual black holes of mass M ∼ 10 6 M ⊙ . Therefore, if the bigravity theory indeed describes physics, the supermassive astrophysical black holes should hide inside the "hairy features".
We have also attentively considered the arguments of Ref. [45]. In brief, this work seems to agree that the solutions exist but judges them physically unacceptable. We analyse the arguments and we think some of them are interesting and should be taken into consideration, but none of them is decisive, so that they should rather be viewed as a conjecture. To understand its origin, we notice that the numerical procedure adopted in that work is not suitable for suppressing the growing at infinity mode, which generates artificial numerical singularities.
This must be reason why the solutions were judged unacceptable in that work. However, no singularities appear within the properly chosen numerical scheme, and we specially adapt our scheme to be able to cope with the arguments of [45]. We shall postpone a more detailed discussion of Ref. [45] untill the end of this text to be able to make a comparison with our results.
The rest of the text is organized as follows. In Section II we introduce the massive bigravity theory of Hassan and Rosen [9]. The field equations, their reduction to the static and spherically symmetric sector, and the simplest solutions are described in Sections III-V. In Sections VI,VII we describe in detail our analysis of boundary conditions at the horizon and at infinity, and then summarize in Section VIII the structure of our numerical procedure. In Section IX we show our solutions for asymptotically flat hairy black holes and also describe the duality relation yielding the solutions above the GL point. After that, we discuss in Section X the perturbations of the hairy backgrounds and the analysis of the negative perturbation modes.
Our discussion culminates in Section XI where we describe various limits and identify regions in the parameter space where the solutions exist and where they are stable. In Section XII we give a brief summary of our results and discuss the arguments of Ref. [45]. The two Appendices contain the description of the desingularization of the equations at the horizon, as well as the complete set of the field equations in the time-dependent case.

II. The ghost-free bigravity
The theory is defined on a four-dimensional spacetime manifold endowed with two Lorentzian metrics g µν and f µν with the signature (−, +, +, +). The action is [9] S[g, f] = 1 where κ 1 and κ 2 are the gravitational couplings, κ is a parameter with the same dimension, and m is a mass parameter. The interaction between the two metrics is expressed by a scalar function of the tensor (the hat denotes matrices) Here the matrix square root is understood in the sense thatγ 2 =ĝ −1f , which can be written in components as If λ a (a = 1, 2, 3, 4) are the eigenvalues of γ µ ν then the interaction potential is where b k are dimensionless parameters while U k are defined by the relations Here [γ] = tr(γ) ≡ γ µ µ and [γ k ] = tr(γ k ) ≡ (γ k ) µ µ . The two metrics actually enter the action in a completely symmetric way, since the action is invariant under (2.5) The action is also invariant under rescalings κ → ±λ 2 κ, b k → ±b k , m → λ m, and this allows one to impose, without any loss of generality, the normalization condition κ = κ 1 + κ 2 .
Varying the action with respect to the two metrics gives two sets of Einstein equations, where κ 1 ≡ κ 1 /κ and κ 2 ≡ κ 2 /κ, and the normalization of κ implies that κ 1 + κ 2 = 1. The source terms in (2.6) are obtained by varying the interaction potential U, where f µα is the inverse of f µα and There is an identity relation following from the diffeomorphism-invariance of the interaction term in the action, where (g) ∇ρ and ∇ρ are the covariant derivatives with respect to g µν and f µν .
Equations (2.6) describe two interacting gravitons, one massive and one massless. This can be easily seen in the flat space limit. Setting g µν = f µν = η µν (the Minkowskoi metric), Eqs.(2.6) reduce to 0 = −m 2 κ 1 (P 0 + P 1 ) η µν , 0 = −m 2 κ 2 (P 1 + P 2 ) η µν , (2.10) with P m ≡ b m +2b m+1 +b m+2 . Therefore, the flat space will be a solution if only the parameters b k fulfil the conditions P 1 = −P 0 = −P 2 . Assuming this to be the case, let us set g µν = η µν + δg µν and f µν = η µν + δf µν where the deviations δg µν and δf µν are small. Linearizing the equations (2.6) with respect to the deviations yieldŝ (2.13) Therefore, one will have m FP = m if P 1 = 1. (2.14) This condition can be solved together with the conditions P 0 = P 2 = −1 implied by (2.10) to express the five b k in terms of two independent parameters, sometimes called c 3 and c 4 , At the same time, the theory has exactly 7 propagating degrees of freedom also away from the flat space limit and for arbitrary b k (see [49][50][51] for its Hamiltonian formulation).
Let us finally pass from the dimensionful spacetime coordinates x µ to the dimensionless ones, This is equivalent to the conformal rescaling of the metrics, after which the field equations (2.6) reduce to where T µ ν and T µ ν are still given by (2.7),(2.8) withγ = ĝ −1f . The Bianchi identities for these equations imply that which is consistent with (2.9). All fields and coordinates are now dimensionless and no trace of the mass parameter m is left in the equations. However, one has to remember that the unity of length corresponds to the dimensionful 1/m, which is the physical length scale.
In what follows we shall be analyzing equations (2.18) without making any assumptions about values of κ 1 , κ 2 and b k . However, when integrating the equations numerically, we shall assume that κ 1 + κ 2 = 1 and choose b k according to (2.15). Therefore, our solutions depend on three parameters of the theory, c 3 , c 4 and η, where κ 1 = cos 2 η, κ 2 = sin 2 η. (2.20)

III. Spherical symmetry
Let us introduce coordinates (x 0 , x 1 , x 2 , x 3 ) = (t, r, ϑ, ϕ) and choose both metrics to be static, spherically symmetric, and diagonal, where dΩ 2 = dϑ 2 + sin 2 ϑdϕ 2 while Q, ∆, R, q, W, U are functions of the radial coordinate r = mr. In fact, this is not the most general form of the spherically symmetric fields, since one could also include the off-diagonal metric element f 01 as shown by Eq.(B.1) in Appendix B. However, in the static case this would imply that (1.1) is the only possible solution [18] (the situation changes in the time-dependent case). Therefore, we choose the static metrics to be both diagonal, which leads to non-trivial solutions.
The tensor γ µ ν in (2.2) then reads and one obtains from (2.7) where Here u = U/R and The independent field equations are plus the conservation condition (g) ∇µ T µ ν = 0 , which has only one non-trivial component, where the prime denotes differentiation with respect to r. The conservation condition for the second energy-momentum tensor also has only one non-trivial component, but this follows from (3.7) due to the identity relation (2.9). As a result, there are 5 independent equations in (3.6), (3.7), which is enough to determine the 6 field amplitudes Q, ∆, R, q, W, U, because the freedom of reparametrizations of the radial coordinate r →r(r) allows one to fix one of the amplitudes.

IV. Field equations
Let us introduce new functions in terms of which the two metrics read The advantage of this parametrization is that the second derivatives disappear from the Einstein tensor and the four Einstein equations (3.6) become The conservation condition (3.7) reads and using Eqs.(4.5), (4.6), this reduces to with dP m = 2 (b m+1 + b m+2 u) (m = 0, 1). (4.10) The conservation condition (3.8) becomes This constraint can be resolved with respect to q to give where Σ(N, Y, R, U) is the (negative) ratio of the coefficients in front of Q and q in (4.9).
As a result, we obtain four differential equations (4.3)-(4.6) plus one algebraic constraint (4.12). The same equations can be obtained by inserting the metrics (4.2) directly to the action (2.1), which gives S = 4π m 2 κ L dtdr , (4.14) where, dropping a total derivative, Varying L with respect to N, Y, Q, q gives Eqs.(4.3)-(4.6), while varying it with respect to R, U reproduces conditions (4.8) and (4.11). The equations and the Lagrangian L are invariant under the interchange symmetry (2.5), which now reads Equation (4.3)-(4.6) contain R ′ and U ′ which are not yet known. One of these two amplitudes can be fixed by imposing a gauge condition, but the other one should be determined dynamically. We need therefore one more condition, and the only way to get it is to differentiate the constraint. Since the constraint should be stable, this gives the secondary constraint: Expressing here the derivatives N ′ , Y ′ , Q ′ , q ′ by Eqs.(4.3)-(4.6) and using the relation (4.13), this condition reduces to where the functions A(R, U, N, Y ) and B(R, U, N, Y ) are rather complicated and we do not show them explicitly. When the radial coordinate change, both R ′ and U ′ change, r →r(r), but the relation (4.18) between R ′ and U ′ remains the same. The secondary constraint can be resolved with respect to U ′ , We can now use the gauge symmetry (4.19) to impose the coordinate condition and then (4.20) reduces to Now, U ′ appears in the right-hand sides of Eqs.
The amplitudes Q, q are determined as follows. Injecting (4.13) to (4.5) yields the equation which determines Q, and when its solution is known, q is determined algebraically from (4.13).
In what follows we shall mainly focus on the three coupled equations (4.23) determining N, Y, U. As soon as their solution is obtained, the amplitudes Q, q are determined from (4.24),(4.13).

V. Analytical solutions
Some simple solutions of the equations can be obtained analytically [18], [52], for which it is convenient to use equation in the form (4.3)-(4.6).

A. Proportional backgrounds
Choosing the two metrics to be conformally related [18], [52], with a constant C, the solution is given by which describes two proportional Schwarzschild-(anti)-de Sitter geometries. The constant C and the cosmological constant Λ(C) are determined by where R Hub is the Hubble radius of our Universe. One way to fulfill this relation is to assume that the graviton mass is extremely small such that the Compton length is of the order of the Hubble radius, However, the relation can also be fulfilled by assuming that Λ is very small, which is possible if there is a hierarchy between the two couplings: κ 1 ≪ κ 2 = 1 − κ 1 ∼ 1. Eq.(5.3) implies then that Λ ∼ κ 1 and that C should be very close to a root of P 1 + CP 2 . The hierarchy between the two couplings is in fact necessary to reconcile with the observations the perturbation spectrum of the massive bigravity cosmology, for which one should assume that [13][14][15][16][17], where M ew ∼ 100 GeV is the electroweak energy scale and M Pl ∼ 10 19 GeV is the Planck mass. As a result, which is of the order of the solar size. However, in what follows we shall not be always assuming κ 1 to be small and shall present our results for arbitrary κ 1 ∈ [0, 1].

B. Deformed AdS background
Choosing U, q to be constant, The g-metric approaches the AdS metric in the leading O(r 2 ) order, but the subleading terms do not have the AdS structure.
It turns out that solutions of Eqs.(4.3)-(4.6) generically approach for r → ∞ either (5.2) or (5.9) (or they show a curvature singularity at a finite r), hence they are not asymptotically flat [18].

VI. Boundary conditions at the horizon
Let us require the g-metric to have a regular event horizon at some r = r H by demanding the metric components g 00 = Q 2 and g rr = N 2 to show simple zeroes at this point. Therefore, we demand that close to this point one has Q 2 ∼ N 2 ∼ r − r H and we consider only the exterior region r ≥ r H where Q 2 > 0 and N 2 > 0. Such a behaviour is compatible with the field equations if only the f-metric also shows a regular horizon at the same place, hence As a result, both metrics share a horizon at the same place r = r H , in agreement with [37,38]. However, the horizon radius measured by the g-metric, r H , can be different from the radius measured by the second metric, U(r H ). We therefore introduce the As a result, the local solutions close to the horizon are expected to have the form the two other amplitudes being The equations then allow one to recurrently determine the coefficients a n , b n , c n , d n , e n . It turns out they all can be expressed in terms of a 1 , which should fulfil a quadratic equation where A, B, C are functions of u, r H and of the theory parameters b k , κ 1 , κ 2 . It turns out that one should choose σ = +1, since choosing σ = −1 always yields singular solutions. Therefore, for a chosen a value of the horizon size r H , the local solutions (6.1),(6.2) comprise a set labeled by a continuous parameter u. These local solutions determine the boundary conditions at the horizon, and they can be numerically extended to the region r > r H .
The surface gravity for each metric is [18] and using the values of the expansion coefficients determined by the equations yields the relation κ g = κ f , hence the two surface gravities coincide, as coincide the Hawking temperatures, One has close to the horizon N(r) ∼ Y (r) ∼ √ r − r H hence the derivatives N ′ and Y ′ are not defined at the horizon. The usual practice would then be to start the numerical integration not at r = r H but at a nearby point r = r H + ǫ. However, although the dependence on ǫ is expected to be small, still its presence in the procedure may lead to numerical instabilities.
This point was emphasised in [45]. This difficulty can be resolved as follows. Setting where N 0 , Y 0 , U 0 are given by (5.2), while the deviations δN, δY, δU approach zero. In the linear approximation, the latters are described by where A, B 1 , B 2 are integration constants and real parts of λ 1 and λ 2 are negative. All of these three perturbation modes vanish for r → ∞, and since the number of equations (4.23) is also three, it follows that the AdS background is an attractor at large r.
b) Solutions extending up to arbitrarily large values of r and asymptotically approaching a deformed AdS background (5.9). The latter is also an attractor at large r. but for larger values of r they deviate from the Schwarzschild metric more and more [18] (this means the Schwarzschild solution is Lyapunov unstable [45]). All of this does not mean that the Schwarzschild is the only asymptotically flat black hole solution. There may be others, but they are not parametrically close to the Schwarzschild solution and should correspond to some discrete values of u which are difficult to detect by a "brute force" method.

VII. Boundary conditions at infinity
Let us suppose the solutions to approach flat space with g µν = f µν = η µν at large r and set In fact, a more general possibility would be to require the g-metric to approach the flat Minkowski metric diag(−1, 1, 1, 1) and the f-metric to approach just a flat metric, as for This would lead to solutions whose Lorentz invariance is broken in the asymptotic region [27,28]. However, we shall not analyze this option.
Inserting (7.1) to (4.23) yields where N N , N Y , N U are the non-linear in δN, δY, δU parts of the right-hand sides D N , D Y , D U in (4.23). Neglecting the non-linear terms, the solution of these equations is and describe the massive graviton, hence they contain the Yukawa exponents (remember that r = mr).
As one can see, among the three modes only two are stable for r → ∞ while the third one diverges in this limit, hence flat space is not an attractor. This is why one cannot get asymptotically flat solutions by simply integrating from the horizon -trying to approach flat space in this way, the unstable mode e +r rapidly wins and drives the solution away from flat space. The only way to proceed is to suppress the unstable mode from the very beginning by requiring the solution at large r to be where the dots denote non-linear corrections. The usual practice would be to neglect the dots and assume that the linear terms approximate the solution everywhere for r > r ⋆ , where r ⋆ is some large value. However, one can check that already the quadratic correction contains an additional factor of ln(r) and hence dominates the linear part for r → ∞. Therefore, nonlinear corrections are important, but if all of them are taken into account, it is not obvious that the solution will remain asymptotically flat.
Fortunately, problems of this kind have been studied -see, e.g., [47]. To take the non-linear corrections into account, the procedure is as follows. Let us express δN, δY, δU in terms of three functions Z 0 , Z + , Z − : Eqs. (7.2) then assume the form Terms on the left in these equations are linear in Z 0 , Z ± , while those on the right are nonlinear. Neglecting the non-linear terms, the solution is Z 0 = 1/r, Z + = e −r , Z − = e +r , and if we set this reproduces the linear part of (7.4). Now, to take the non-linear terms into account, one converts Eqs.(7.6) into the equivalent set of integral equations, where r ⋆ is some large value. These equations determine the solution for r > r ⋆ , and they are solved by iterations. To start the iterations, one neglects the non-linear terms, which gives the configuration (7.7). The next step is to inject this configuration to the integrals, which gives the corrected configuration, and so on. In practice, one introduces variables x = r ⋆ /r andx = r ⋆ /r assuming values in the interval [0, 1], and then one discretizes the interval to compute the integrals.
To see the convergence of the iterations, we compute for each Z and for each discretization shown on the right panel in Fig.1: the amplitudes Z 0 and Z ± against x = r ⋆ /r (for r ⋆ = 25 r H ).
One can see that the amplitude Z − is always small but non-vanishing, and that all the three amplitudes vanish for x = 0, hence the solution is indeed asymptotically flat.
This yields an asymptotically flat solution in the region r > r ⋆ . To extend this solution to the region r H < r < r ⋆ one only needs its values at r = r ⋆ , To recapitulate, the described above procedure yields the boundary values for the fields at a large r ⋆ and makes sure that the solution for r > r ⋆ exists and is indeed asymptotically flat.
It is worth noting that the parameter A determines the ADM mass, (which determines also Q(r), q(r)) with the following boundary conditions. At the horizon Far from the horizon, at r = r ⋆ ≫ r H , one has where Z 0 (r ⋆ ) and Z − (r ⋆ ) are functions of A, B determined by (7.9) via iterating the integral equations (7.8).
As a result, we have the boundary conditions at r = r H labeled by u and the boundary conditions at r = r ⋆ labeled by A, B. We use them to construct solutions in the region r H ≤ r ≤ r ⋆ . To this end, we choose some value of u and integrate numerically the equations starting from r = r H as far as some r = r 0 < r ⋆ and we obtain at this point some values which will depend on r H and u: Then we choose A and B and numerically extend the large r data (8.2) form r = r ⋆ down to r = r 0 , thereby obtaining (A, B). within the Newton-Raphson procedure [39]. However, unless the input configuration is close to the solution, the numerical iterations do not converge, hence some additional information is necessary to specify where to start the iterations.
As explained in the Introduction, the additional information is provided by the stability analysis of the Schwarzschild solution (1.2) [33,42]. In this analysis one considers the two metrics of the form (4.2) with where S = 1 − r H r while the perturbations δQ, δN, δq, δY, δU, δα are assumed to be small and depend on t, r. It turns out that at the GL point, for r H = 0.86, the perturbation

IX. Asymptotically flat hairy black holes
Applying the procedure outlined above, we were able to construct asymptotically flat hairy solutions. We confirm the results of Ref. [44] and obtain many new results.
First of all, we find that for r H approaching from below the GL value, of its Riemann tensor diverge where the zeros are located. An example of this is shown on the lower two panels in Fig. 2. As seen in Fig. 3 where results for different values of η are shown, small r H solutions become singular when η approaches π/2, but they are regular for small η. Fig. 2 are shown up to large but still finite values of the radial coordinate, r/r H ≤ 100 or r/r H ≤ 1000. What is shown is the combination of the solutions of differential equations (4.23) in the region r H ≤ r ≤ r ⋆ and of the integral equations (7.6) for r > r ⋆ where r ⋆ /r H = 25. At the same time, our procedure yields solutions in the whole region r ∈ [r H , ∞).

Solutions in
Introducing the compactified radial variable we plot in Fig. 3 the amplitudes N, Y, Q, q against ξ. As seen, the amplitudes approach unity as ξ → 1 (same is true for U ′ ) hence the solutions are indeed asymptotically flat. The solution for the f-metric is shown on the lower right panel in Fig. 3. Similarly, for η = 0 one has κ 2 = 0 and the f-metric is Schwarzschild, while the g-metric is a solution of the massive gravity on the Schwarzschild background shown on the upper left panel in Fig. 4. One should emphasise that the radii of the background Schwarzschild black holes for η = 0 and for η = π/2 are not the same. For example, for η = π/2 the Schwarzschild black hole has r H = 0.18 for the solution shown in Fig. 3, while for η = 0 the event horizon size is determined not by r H but by U H = ur H where u ≈ 7 (as seen in Fig.4) hence this time the Schwarzschild black hole is much larger. As a result, solutions on these different backgrounds look quite different -the solution for the f-metric on the lower right panel in Fig. 3 shows zeros hence is singular, while the solution for the g-metric on the upper left panel in Fig. 4 is regular.  In case I asymptotically flat hairy black hole exist only for 0 < r min H < r H < 0.86 hence they cannot be arbitrarily small. In case II they exist for any 0 < r H < 0.86, although their f-metric may be singular for small r H . We shall explain below in Section XI what happens when r H approaches the lower bound.

Duality relation
The results described above in this Section essentially reproduce those of Ref. [44], the only important difference is that we show solutions for different values of η, whereas Ref. [44] shows them only for η = π/4. However, starting from this moment and in the following two Sections we shall be describing new results.
Ref. [44] finds solutions only below the GL point, for r H ≤ 0.86. At the same time, the consistency of the procedure requires that there should be asymptotically flat hairy black holes also for r H > 0.86. This follows from the symmetry (4.16) of the equations, which now reads then forη = π/2 − η,c 3 = 3 − c 3 ,c 4 = 4c 3 + c 4 − 6 there should be the "dual" solution described bỹ Q(r) = q(w(r)),q(r) = Q(w(r)),Ñ(r) = Y (w(r)),Ỹ (r) = N(w(r)),Ũ(r) = w(r), where w(r) is the function inverse for U(r) such that U(w(r)) = r. This duality correspondence relates to each other black holes of difference size, since (4.24) has the horizon at r = r H while the horizon of (9.12) is located where w(r) = r H that is at r =r H = U(r H  It is unclear why solutions with r H > 0.86 were not found in [44].
The duality is in fact a powerful tool for studying the solutions, because sometimes their properties may look puzzling in one description but become completely obvious within the dual description.

X. Stability analysis
In this section we analyze the stability of the hairy solutions by studying their perturbations within the ansatz described in Appendix B,  Linearizing similarly the G 0 1 (f ) = κ 2 T 0 1 equation yields a linear relation between δα(r), δY (r), δU(r). Using these two algebraic relations one finds that the three equations G 0 0 (g) = κ 1 T 0 0 , G 0 0 (f ) = κ 2 T 0 0 and Taking all of this into account and linearizing similarly the remaining 3 equations G 1 1 (g) = κ 1 T 1 1 , G 1 1 (f ) = κ 2 T 1 1 and (g) ∇µ T µ 1 = 0, one finds that all 6 perturbation amplitudes δQ(r), δq(r), δN(r), δY (r), δU(r) and δα(r) can be expressed in terms of a single master amplitude Ψ(r) subject to the Schrödinger-type equation, Here the potential V (r) is a rather complicated function of the background amplitudes that we do not show explicitly, and the tortoise radial coordinate r * ∈ (−∞, +∞) is defined by the relation dr * = 1 a(r) dr, (10.7) where the function a(r) (also complicated) varies from 0 to 1 as r changes from r H to ∞. The potential V always tends to zero at the horizon, for r * → −∞, and it approaches unit value at infinity, for r * → +∞. One should remember that our dimensionless variables are related to the dimensionfull ones via r =mr, r H = mr H , V = m 2 V, ω = ω/m.
In the Schwarzschild case, when the static background amplitudes are Q = q = N = Y = 1 − r H /r and U = r, one has a(r) = Q 2 (r) and the potential reduces to (10.8) in agreement with Ref. [42]. In the flat space limit r H → 0 this reduces to V (r) = 1 + 6/r 2 , which is the potential of a massive particle of unit mass (in units of the graviton mass) with spin s = 2.
Eq.(10.6) defines the eigenvalue problem on the line r * ∈ (−∞, +∞). Solutions of this problem with ω 2 > 0 describe scattering states of gravitons. In addition, there can be bound states with purely imaginary frequency ω = iσ and hence with ω 2 = −σ 2 < 0. For such solutions the wavefunction Ψ is everywhere bounded and square-integrable, because one has e +σr * ← Ψ → e − √ 1+σ 2 r * as −∞ ← r * → +∞, respectively. Such bound state solutions grow in time as e iωt = e ±σt . Therefore, they correspond to unstable modes of the background black holes.

Computing the eigenfrequencies
Our aim is to investigate a potential existence of negative modes with ω 2 < 0 in the spectrum of the eigenvalue problem (10.6). If such modes exist, then the background black holes are unstable. If they do not exist, then the black holes are stable with respect to spherically symmetric perturbations, which strongly suggests that they should be stable with respect to all perturbations. Indeed, in most known cases the S-channel is usually the only place where the instability can reside (of course, this should be proven case to case).
The first thing to check is the shape of the potential V (r), because if it is everywhere positive, then there are no bound states. We therefore show in Fig.7   values in its vicinity, and then approaches unity as r → ∞. Therefore, since the potential is not positive definite, bound states may exist, but their existence is not yet guaranteed.
We know that a bound state certainly exists for the Schwarzschild background with r H < 0.86 [33,42]. When looking at the potentials for the hairy solution with r H = 0.78 in Fig.7, we notice that they are close to the Schwarzschild potential, hence a bound state could exist for these potentials as well.
In order to know whether bound states exist or not, one can use the well-known Jacobi criterion [54] and construct the solution of the Schrödinger equation ( solution (see [55,56] for a review on BHs perturbation theory and the tools that can be used to solve the perturbation equation).
The eigenfunctions Ψ against the ordinary radial coordinate r are shown in Fig.8. They vanish at the horizon then show a maximum, sometimes very close to the horizon, and then zero again for r → ∞.
As a result, we find the eigenfrequencies ω 2 < 0 and the corresponding eigenfunctions Ψ(r) for all hairy black holes obtained in [42]. Therefore, all these solutions are unstable.
It is worth emphasising that all of them correspond to the particular choice η = π/4, hence κ 1 = κ 2 = 1/2. In order to test our method, we have also computed the negative mode for the Schwarzschild solution [42].
As seen in Fig.9, the absolute value of the negative mode eigenvalue for the Schwarzschild solution is always larger than that for the hairy solutions. Therefore, the instability growth rate for the hairy black holes is not as large as for the Schwarzschild solution. In all cases, since one has ω = ω/m where ω is the dimensionful physical frequency, the instability growth time is 1/ω = 1/(ωm). If we assume the graviton mass m to be very small and given by (5.5), then the instability growth time will be cosmologically large, hence the instability will not play any role. However, as we shall see below, it is preferable to assume that 1/m ∼ 10 6 km according to (5.7), in which case the instability growth time will be of the order of 10 3 seconds, hence the instability is dangerous and should be avoided.
As seen in Fig.9, the eigenvalue ω 2 (r H ) < 0 approaches zero when r H → 0.86, therefore all hairy black holes become then stable. However, they are no longer hairy in this limit -they zero, as seen in the insertion in Fig.9.
The instability of hairy black holes is in fact somewhat puzzling, since it is unclear what they may decay into. Since the hairy solutions with r H < 0.86 are more energetic than the Schwarzschild solution, they probably decay via absorbing and/or radiating away their hair and approaching the "bald" Schwarzschild solution. However, the latter is also unstable for r H < 0.86 and should decay in its turn.
The perturbative instability of the Schwarzschild solution in the massive bigravity theory is mathematically equivalent [33] to the Gregory-Laflamme instability of the vacuum black string in D=5 [43]. It is known that the non-linear development of the latter leads to the formation of an infinite string of "black hole beads" in D=5, but the event horizon topology does not change [57]. This fact being established within the D=5 vacuum GR, a similar scenario is not possible in the D=4 bigravity theory, hence the fate of the bigravity black holes should be different. One possibility is that the black hole radiates aways all of its energy within the S-channel (some radiative solutions are known explicitly [58,59]), but it is unclear what happens to the horizon -if it disappears or not. In GR the horizon cannot disappear via a classical process [60], but in the bimetric theory the situation might be different.
Remarkably, we find that these puzzling issues are not omnipresent and the black holes can be stable for η = π/4. In Fig.10 we show ω 2 against η for several values of r H for solutions with c 3 = −c 4 = 2. One can see that ω 2 (η) < 0 approaches zero and the negative mode disappears when η approaches π/2. The background solutions then do not disappear and become stable.
However, the f-metric becomes then singular and shows oscillations of the amplitudes Y, q, U ′ , as seen in Fig.3, hence this case is not interesting. At the same time, ω 2 approaches zero also when η become small enough, if only r H is also small, as seen in Fig.10. The background solutions then become stable and remain regular.
Summarizing the above discussion, for some parameter values the hairy black holes are unstable, but for other parameter values they can be stable and regular. Below we shall describe a parameter choice leading to a large set of stable and regular solutions.  where all solutions bifurcate with the Schwarzschild solution but the f-metric for these solutions is not Schwarzschild, even though both metrics have the same mass. (The f-metric for η = π/2 shown in Fig. 3 (right-lower panel) is singular since q, Y oscillate, but the digram in Fig. 11 takes into account only values of r H corresponding to regular solutions.) For η = π/2 the mass depends non-linearly on r H .
Introducing the mass function M(r) via N 2 (r) = 1 − 2M(r)/r, the G 0 0 (g) = κ 1 T 0 0 Einstein equation in (2.18) assumes the form from where the ADM mass Here the "bare" mass M bare = r H /2 is determined only by the horizon radius and coincides with the mass of the η = π/2 solution, while the mass M hair expressed by the integral is the contribution of the massive hair distributed outside the horizon. As one can see in Fig. 11, one has M > r H /2 if r H < 0.86, hence the "hair mass" is positive and the hairy solutions are more energetic than the bare Schwarzschild black hole. However, the mass of the hair becomes of that close to the tachyon limit, with D ∼ 10 −6 (right). One has S 2 = 1 − r H /r. negative above the GL point, where r H > 0.86, and the hairy solutions are then less energetic than the bare one. Therefore the energy density T 0 0 can be negative. In fact, there are no reasons for which the standard energy conditions should be respected within the bigravity theory.
Each curve in Fig. 11 is defined only in a finite interval r H ∈ [r min H (η), r max H (η)]. The terminate points of the curves correspond to limits beyond which one cannot continue, and it is very instructive to understand what happens in these limits.

A.
The singular lower limit r H → r min H (η) It turns out that r min H (η) > 0 if η > 0.93 hence κ 1 ≤ 0.35 (these values would be different for different c 3 , c 4 ). In this case the f-metric of the limiting for r H → r min H (η) solutions becomes singular because the q, Y -amplitudes start to develop a second zero outside the horizon. This phenomenon has already been discussed above and is similar to what is shown in Figs. 2, 3.
Solutions exist also for r H < r min H (η) but they are all singular in the same sense and should be excluded. Therefore, one should require that r H > r min H (η). If η < 0.93 then r min H (η) = 0 and the solutions extend down to arbitrarily small values of r H . Remarkably, as seen in Fig. 11, the mass M does not vanish when r H → 0 but approaches a finite value, even though the bare mass M bare = r H /2 → 0. Therefore, all mass is contained in the "hair mass" in this limit, hence something remains when the horizon size r H shrinks to zero. A similar phenomenon is actually known, since in many non-linear field theories there are solutions describing a small black hole inside a soliton (for example, inside the magnetic monopole) [35]. Sending the horizon size to zero the black hole disappears, but its external non-linear matter fields remain and become the gravitating soliton containing in its center a regular origin instead of the horizon. Therefore, the r H → 0 limit of such a hairy black hole is a regular soliton.
One may think that the situation could be similar also in our case, hence there should be a limiting configuration to which the black hole solutions approach pointwise when r H → 0.
It seems that such a limiting solution indeed exists, however, it is singular and not of the regular soliton type. First, as seen in Fig. 11, the value of U H which determines the size of the f-horizon remains finite when r H → 0, hence the f-geometry remains black hole even in the limit. Secondly, as seen in Fig. 12, one has N 2 /S 2 ∼ r for r ≤ 0.5 for a solution with a very small r H . However, one has S = 1 − r H /r → 1 as r H → 0, hence one has in this limit N 2 ∼ r and the limiting form of the g-metric is something like a "zero size black hole". The numerical profiles shown in Fig. 12 suggest this limiting configuration to have the following structure at small r: The g-geometry is singular since its Ricci invariant R(g) = 2/r 2 + O(1/r) at small r, but the f-geometry remains of the regular black hole type because U does not vanish. Curiously, the temperature remains finite for r H → 0 and is aways the same for both metrics. The limiting g-temperature can be formally computed by assuming N 2 = αr, Q 2 = βr with α ≈ 0. decreases and vanishes for some r H = r tach H (η), then it becomes positive again, decreases again and vanishes for the second time for r H = r max H (η) > r tach H (η), after which it becomes negative and the procedure stops. Specifically, it turns out that the determinant of (6.3) factorizes, 8) in Fig. 11, the mass of all solutions, even of those with r H → 0, is actually bounded by the limiting values of the mass of the η = π/2 solutions, 1 2 r min H (η = π/2) ≈ 0.23 < M <  .7), which is consistent with the cosmological observations if κ 1 is parametrically small as expressed by (5.6), yields a physically acceptable result. The masses of the hairy black holes are then close to the mass of the Schwarzschild black hole of radius ∼ 10 6 km, that is ∼ 10 6 M ⊙ . This is typical for supermassive astrophysical black holes observed in the center of many galaxies.
Therefore, if the massive bigravity theory indeed describes physics, the "hairy features" can be present only in supermassive black holes and not in black holes of a smaller mass. Let us now collect all the facts together. The diagram in Fig. 13 shows the region in the (r H , η) plane within which there are regular hairy black hole solutions. The low boundary of this region at η = 0 corresponds to solutions whose f-metric is Schwarzschild, while the upper boundary at η = π/2 corresponds to solutions whose g-metric is Schwarzschild. The left boundary consistes of two pieces -the curve separating the left upper corner where the solutions are singular because their f-metric shows additional zeroes, and the portion of the η-axes corresponding to the "zero size black holes". Finally, the right boundary corresponds to the "tachyon limit" beyond which the solutions become complex-valued. where ω 2 > 0 and hence the solutions are stable, from sectors where ω 2 < 0 and the solutions are unstable. There are altogether two stable and two unstable sectors. It is worth noting that the stability region is now much larger than for solutions with c 3 = −c 4 = 2 considered in the previous Section. One also notices that the tachyonic solutions are in the unstable sector.
Finally, the diagram shows the "physical region" corresponding to physically acceptable solutions. As explained above, for such solutions the coupling κ 1 = cos 2 (η) should be very small for their mass not to be too large. Therefore, η should be very close to π/2. The solutions should be stable, hence they should correspond to the sector where ω 2 > 0. These conditions specify the physical region to be the thick (green online) line at the top of the diagram.
Physical solutions are therefore described by the g-metric which is extremely close to Schwarzschild, since The "hairy features" of the solutions hidden in the f-metric should be difficult to observe, unless in violent processes like black hole collisions producing large enough T µν (g, f ) to overcome the 10 −34 suppression. Summarizing, the static bigravity black holes should be extremely similar to the GR black holes, but their strong field dynamics is expected to be different.  Let us now see how the described above solutions look after the duality transformation (9.10). This transformation converts the parameter values (11.1) into (11.2), flips the sign of η − π/4 and swaps the Q, N, r with q, Y, U. Graphically, this amounts to relabelling the functions and plotting them against U instead of r. The ADM mass and temperature are invariant under duality. The stability property also does not change since, for example, if a solution is unstable and admits growing in time perturbations, then its dual version will contain the same growing modes and hence will be unstable as well. The solutions still exist in a finite interval r H ∈ [r min H (η), r max H (η)], but the limits are now different. The lower limit r min H (η) can no longer extend down to zero but always corresponds to the tachyon solutions with vanishingly small horizon determinant D. The upper limit r max H (η) corresponds for small η to solutions whose g-metric (which used to be the f-metric before the duality) starts being singular. For larger values of η the upper limit r max H (η) corresponds rather to the point where the two different solutions with the same r H but with different U H merge to each other. The zero size black holes still exist but now correspond to internal points of We notice that the solutions below the GL point, for r H < 0.86, are still more energetic than the solution with η = π/2, hence their "hair mass" M hair is positive, whereas above the GL point it becomes negative. Finally, Fig. 15 shows the existence diagram in the (r H , η) plane, together with the stability regions. The diagram looks quite different as compared to that in Fig. 11, although it corresponds to essentially the same solutions, up to the duality transformation. Although the duality does not change stability, it interchanges positions of the stability sectors. Therefore, the physical region corresponding to stable solutions with η close to π/2 is now above the GL point, where the hair mass is negative. The physical solutions are again characterised by the g-metric that is extremely close to Schwarzschild, but the novel feature is that now for each value of r H from the physical region there are two different solutions whose g-metrics are almost the same but the f-metrics are different.
Although the described above features correspond to just two particular sets of values of c 3 , c 4 , our numerics indicate that they are essentially generic. The solutions always exist only within a finite region of the (r H , η) plane whose borders correspond to the tachyon limit, to the singular solutions, etc. The physical solutions always correspond to values of η very close to π/2 and to values of r H close to the GL value r H = 0.86. As a result, the hairy solutions always have the g-geometry close to that for the supermassive astrophysical black holes.

XII. Concluding remarks
To recapitulate, we presented above a detailed analysis of static and asymptotically flat black holes in the ghost-free massive bigravity theory. Extending the earlier result of [44], we find that for given values of the theory parameters c 3 , c 4 , η and for a given even horizon The dimensionless ADM mass of the hairy solutions can be neither too large nor too small and is always not far from the GL value, M = 0.43. Therefore, to avoid the hairy black holes being unphysically heavy, one is bound to assume the massive graviton Compton length to be of the order of 10 6 km, in which case the agreement with the cosmological data is achieved by assuming that cos 2 η ∼ (M ew /M Pl ) 2 ∼ 10 −34 . The stable hairy black holes are then described by the g-metric which is extremely close to the Schwarzschild metric, although the f-metric is different, and they have the size and mass close to those of ordinary black holes of mass M ∼ 10 6 M ⊙ . As a result, if the bigravity theory indeed applies to describe physics, the supermassive astrophysical black holes should hide inside the "hairy features", which should become manifest in violent processes like black hole collisions.
Finally, we should discuss the paper [45] that also considers black holes in the ghost-free massive bigravity theory. This paper presents essentially the same classification of different types of black holes as the one previously given in [18] but in a more refined way, extending it and paying attention to some subtle points. The paper addresses in particular the issue of convergence of the solutions to the flat background in the asymptotic region. Among other things, it claims that the Schwarzschild solution is the only asymptotically flat black hole in the theory. At the same time, the paper does not contain a rigorous proof of this statement but gives just a number of plausibility arguments, so that the claim should rather be viewed as a conjecture, as actually explicitly stated in some places of [45]. These arguments are as follows.
First of all, it was emphasised in [45] that the usual practice of starting the numerical integration not at the horizon r = r H , which is a singular point of the differential equations, but at a regular nearby point r = r H +ǫ, as was done in [44], could in principle lead to numerical instabilities. We agree with this, and it is for this reason that we use the desingularization procedure (described in Appendix A below) which allows us to start the numerical integration exactly at r = r H (initial conditions exactly at r = r H were described also in [45]).
The paper [45] makes also another remark concerning the behaviour at the horizon. It is known that in order to be able to cross the horizon, for example when studying geodesics, one cannot use the Schwarzschild coordinates and one should introduce instead regular at the horizon coordinates. These can be, for example, Eddington-Finkelshtein (EF) coordinates in which g 00 = g 11 = 0, g 01 = g 10 = 0. It was noticed in [45] that the f-metric, when expressed in the same coordinates, generically does not have the same form, since it has f 11 = 0, hence the two metrics cannot be simultaneously EF. We understand this, but this does not invalidate the background solutions (Ref. [45] agrees on this). The horizon geometries are regular, and if one wishes, one can use the same boundary conditions at the horizon to integrate inside the horizon to recover the interior solutions. Within the parametrization described in Appendix A, this is achieved by simply changing the sign of the numerical integration step.
Next, small initial deviations from the Schwarzschild solution via setting at the horizon u = U H /r H = 1+ǫ were considered in [45]. Integrating the equations toward large r then yields metrics whose components diverge as r → ∞ instead of approaching finite values. In addition, they show curvature singularities at a finite r where the amplitudes q, Y, U ′ change sign. This observation, made already in [18], shows that there are no regular and asymptotically flat solutions in a small vicinity of the Schwarzschild solution. However, there can be regular solutions corresponding to u considerably deviating from unit value.
Finally, the paper [45] reproduces and analyzes (in Appendix A) one of the asymptotically flat solution found in [44]. It chooses for this a singular solution with oscillating q, Y, U ′ amplitudes -in order to illustrate the pathologies that may arise. The regular solutions are not specially considered in [45] because, within the numerical approach adopted in that work, they also look pathological. Specifically, Appendix D of [45] describes the numerical method -a straightforward integration starting from the horizon with the standard routine of MATHEMATICA. This adequately produces the solution with a given precision, but only within a finite range of the radial coordinate r. If one integrates farther trying to approach flat space, then the growing mode Ce +r generically present in the solution leads to a rapid accumulation of numerical errors triggering a numerical instability. Trying to suppress this mode by adjusting the horizon boundary conditions, one typically observes the derivatives of some functions in the solution growing without bounds at some finite r. Precisely this type of behaviour at the end of the integration interval is seen in Fig.11 in [45]. Therefore, one cannot get asymptotically flat solutions within the numerical scheme adopted in [45], since pathological features arise in this scheme inevitably. This must be the reason behind the conviction that such solutions do not exist. However, all pathologies can be eliminated within the more elaborate numerical scheme described above -via suppressing the growing Ce +r mode from the very beginning.

Appendix
In the following two appendices we describe the desingularization procedure and the timedependent case.
A. Desingularization at the horizon The horizon r = r H is a singular point of the differential equations -the derivatives N ′ and Y ′ expressed by Eqs. (4.23) are not defined at this point. The usual practice to handle this difficulty is to use the local power series expansions (6.1),(6.2) to start the numerical integration not exactly at r = r H but at a nearby point with r = r H + ǫ where ǫ is a small number. One may then hope that the results will not be very sensitive to the value of ǫ. However, in such an approach ǫ remains an arbitrary parameter not defined by any prescription. This inevitably affects the stability of the numerical procedure, which becomes evident when one studies the dependence of the solutions on the theory parameters. where C 1 = (r − r H ν 2 − κ 1 r 3 P 0 ) y − κ 1 r 3 P 1 U ′ ν , C 2 = ν r 2 (1 − κ 2 r 2 P 2 ) U ′ − κ 2 r 4 P 1 y − r H Uν y 2 . (A.3) Summarizing the above discussion, the equations in the desingularized form read ν ′ = − ν 2r + C 1 2νyr 2 S 2 ≡ F ν (r, U, ν, y), y ′ = − yU ′ 2U + C 2 2νyr 2 US 2 ≡ F y (r, U, ν, y), U ′ = D U (r, U, Sν, Sy) ≡ F U (r, U, ν, y), (A. 13) where C 1 and C 2 are defined by (A.3) while D U is the same as in (4.23). These equations apply for r > r H , while at r = r H they should be replaced by  6). This formulation allows one to start the integration exactly at the horizon r = r H and then continue to the r > r H region.

B. Field equations with time dependance
Let us allow both metrics to depend on time, assuming that they are still spherically symmetric. The gauge freedom of reparametrizations of the t, r coordinates can be used to make the g-metric diagonal, but the f-metric will in general contain an off-diagonal term. The two metrics can be written as [10] ds 2 g = −Q 2 dt 2 + dr 2 ∆ 2 + R 2 dΩ 2 , where dΩ 2 = dθ 2 + sin 2 θ dφ 2 and Q, q, ∆, W , α, U, R are functions of r and t.
One can check that the tensor has the property γ µ σ γ σ ν = g µ σ f σ ν . This tensor is used to compute the energy-momentum tensors T µ ν and T µ ν in (2.7).
One can redefine the two amplitudes similarly to (4.1) where the prime denotes the derivative with respect to r, and one can impose the gauge condition R = r. (B.4) As a result, the independent field equations (2.18) become G 0 0 (g) = κ 1 T 0 0 , G 1 1 (g) = κ 1 T 1 1 , G 0 1 (g) = κ 1 T 0 1 , G 0 0 (f ) = κ 2 T 0 0 , G 1 1 (f ) = κ 2 T 1 1 , G 0 1 (f ) = κ 2 T 0 1 , (B.5) plus two non-trivial components of the the conservation condition Here one has explicitly where the dot denotes the partial derivative with respect to t, while where P m are defined in (3.5). The components of the second stress-energy tensor are T 0 0 = − r 2 NU 2 A P 1 qY + P 2 α 2 N 2 QY + qNU ′ , T 1 1 = − r 2 U 2 A P 1 QU ′ + P 2 α 2 NQY + qU ′ , where A = NQY α 2 + qU ′ . The components of the Einstein tensor for f µν , are complicated: (B.10) Finally, there are two non-trivial components of the conservation law, where dP m are defined in (4.10). Eqs.(B.5), (B.6) comprise a system of 8 equations for 6 functions Q, q, ∆, W , α, U. For this system not to be overdetermined, only 6 equations out of 8 should be independent. As shown in Section X, this indeed happens at least for small α, when the perturbative analysis of the equations shows that some of them coincide.