{'fraction_non_alphanumeric': 0.04373711637197035, 'fraction_numerical': 0.010959886303983202, 'mean_word_length': 4.321886072499319, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Classification and asymptotic structure of black holes in bimetric theory

We study general properties of static and spherically symmetric bidiagonal black holes in Hassan-Rosen bimetric theory. In particular, we explore the behaviour of the black hole solutions both at the common Killing horizon and at the large radii. The former study leads to a new classification for black holes within the bidiagonal ansatz. The latter study shows that, among the great variety of the black hole solutions, the only solutions converging to Minkowski, Anti-de Sitter and de Sitter spacetimes at large radii are those of General Relativity, i.e., the Schwarzschild, Schwarzschild-Anti-de Sitter and Schwarzschild-de Sitter solutions.


I. INTRODUCTION
General Relativity (GR) is a remarkably successful theory of gravitational interactions. It has passed many experimental tests so far both in the weak-and the strong-field regime [1,2]. Yet, GR cosmology cannot give a satisfactory explanation for the physical origin of the accelerated expansion of the Universe [3]. This issue contributes to the socalled "cosmological constant problem" [4], one of the most intriguing open problems in modern physics. It provides a strong motivation for studying alternative theories of gravity, which could shed light on the physical cause of the acceleration.
Moreover, GR cosmology needs the introduction of an unknown matter component, called the "dark matter", which is required for explaining many astrophysical and cosmological observations and has well established theoretical motivations [5][6][7]. Despite the effort that has been dedicated to revealing the dark matter, its origin is still unknown. So far, only gravitational interactions between the dark matter and the Standard Model sector have been detected; this suggests that the physical explanation for the origin of the dark matter can still be found within the realm of modified gravity. Some examples of this research front can be found in [8][9][10][11].
It is well-known that GR is the unique, nonlinear theory describing the self-interaction of a massless spin-2 tensor field (the "graviton", if quantised) and its interaction with non-gravitational fields [12,13]. Therefore, a natural way of generalising it is to add a massive spin-2 tensor field to the theory. The first attempt to describe massive spin-2 fields was in 1939 when Fierz and Pauli formulated a linear theory describing the dynamics of a massive spin-2 field [14,15]. In 1972, Boulware and Deser claimed that any nonlinear generalisation of the Fierz-Pauli theory has to suffer from a ghost mode, the Boulware-Deser (BD) ghost [16,17]. However, more recently a ghost-free nonlinear theory of a massive spin-2 field -de Rham-Gabadadze-Tolley (dRGT) massive gravity -was proposed [18,19]. The dRGT theory involves two rank-2 tensors (usually called metrics, for the sake of familiarity), one dynamical and the other a non-dynamical reference metric. This model has five propagating modes associated with the massive spin-2 field. A As already pointed out in [33], GR solutions are not unique since, given a set of values for the parameters in the action, we can have different conformal factors between the two metrics. Each conformal factor corresponds to one value of the cosmological constant for the GR solutions. If two conformal factors lead to Λ = 0, then we will have two Schwarzschild solutions. Similarly for SAdS when we have different Λ < 0 and for SdS when we have different Λ > 0. Note that this cosmological constant is different from the one defined in the action of GR. In this sense, a GR solution will not be completely specified by the global parameters of the theory and its mass, since, for fixed values of these, we can have a discrete and finite number of observationally distinguishable solutions, specified by the the different conformal factors. In addition, turning to non-GR solutions, to completely specify the BH, we need to specify possible deviations from the GR solutions regarding the initial conditions at the Killing horizon, even though such deviations are diverging from the GR solutions at large radii. Whether the differences between these solutions should be denoted hairs, depends on the exact definition of hair employed. We note, however, that for fixed values of the global action parameters and the BH mass, all found BH solutions always diverge from each other at large radii.
Given that: (i) the only bidiagonal BH solutions in HR bimetric theory asymptotically converging to Minkowski, dS and AdS spacetimes are GR solutions with conformal metrics, (ii) bidiagonal BH solutions display dynamical linear instabilities, (iii) spherically symmetric solutions including matter sources, e.g. star solutions, have non-conformal metrics, the outstanding remaining question is: What is the endpoint of gravitational collapse of matter in HR bimetric theory?

II. THE HASSAN-ROSEN BIMETRIC THEORY
The action of the HR bimetric theory in vacuum is [22], where M g , M f are the reduced Planck masses and R g , R f are the Ricci scalars respectively for g µν and f µν . The kinetic terms for g µν and f µν have the standard Einstein-Hilbert form. The interaction potential is given in terms of the square root matrix S := g −1 f (so that S µ ρ S ρ ν = g µρ f ρν ) through the elementary symmetric polynomials e n (S). These polynomials are the scalar invariants of S and can be expressed in the following way, e 0 (S) := 1, e n (S) := S [µ1 µ1 · · · S µn] µn , n ≥ 1, which can be expanded (noting that e n (S) ≡ 0 for n > 4), The β n parameters are arbitrary real numbers, whereas m is a parameter having the dimension of mass, setting the energy scale of the interaction. Since we will not consider any source of matter in this work, we do not include any matter Lagrangian in the action. Equations of motion in vacuum are obtained by varying the action (1) with respect to g µν and f µν separately, yielding, where R g µν and R f µν are the Ricci tensors of g µν and f µν , respectively, and, Finally, the Bianchi constraints read,

III. BIDIAGONAL ANSATZ AND FIELD EQUATIONS
This section is organised as follows. In section III A, we define the systems of our interest and clarify the assumptions. We then present a general result, stated as a proposition, about the relation between the curvature singularities and the regularity of the square root. We remark that this result is valid generically in HR bimetric theory, not only for BH solutions. In section III B, we study for the first time the physics of bidiagonal BHs at their common Killing horizon and propose a new classification of the solutions. In section III C, we expand the field equations as nonlinear ODEs and discuss their properties. Finally, in section III D, we consider the exact initial conditions at the Killing horizon (not determined in previous works).

A. Static and spherically symmetric bidiagonal BHs
Our goal is to find all possible bidiagonal static and spherically symmetric BH solutions with a common Killing horizon for the two metrics. The reason to investigate only the bidiagonal solutions is motivated by the fact that all non-bidiagonal solutions are known to correspond to GR solutions (Schwarzschild, SdS and SAdS) as shown, for example, in [27]. Whether this is the case also for bidiagonal solutions is the main question we want to address in this paper.
As proved in [44,45], if two static and spherically symmetric metrics are diagonal in a common coordinate system and they describe smooth geometry, they must share their Killing horizons. The statement is very general as it uses the invariant scalars composed from different contractions of two metrics, which are independent from the chosen coordinate system. This result will be later used to construct the ansatz for the metric f . We first consider the g-sector.
The g-sector is endowed with the static timelike Killing vector field K. Since the geometry of g is static and spherically symmetric, assuming that the vector field K is unique, we can adopt coordinates as prescribed in [49]: We begin by selecting a spacelike surface S orthogonal to the orbits of the one-parameter isometry generated by K. We use t to denote this parameter and refer to it as the 'time coordinate'. The image of this isometry are the surfaces S t also orthogonal to K, which we label by t so that all the points at S t have the same t. Now, we select a two-sphere on the spacelike surface S t and choose angular coordinates (θ, φ) where the metric on each 2-sphere takes the form r 2 dθ 2 + sin 2 θ dφ 2 . Although the parameter r is proportional to the square root of the total area of the sphere, it will be referred to as the 'radial coordinate' of the sphere. We carry these spherical coordinates to the other spheres of S t along the radial geodesics orthogonal to the selected sphere. Assuming a non-vanishing ∇ µ r = 0, we select (r, θ, φ) as the coordinates on S t , resulting in the spacetime coordinates (t, r, θ, φ). Allowing further for the reparametrisation of r in terms of a new coordinate ξ (with the unitary gauge r = ξ), we arrive at the coordinate chart x µ = (t, ξ, θ, φ). The introduction of ξ will become obvious later, when we obtain a system of non-linear ordinary differential equations (ODEs) as the field equations. In short, the system will become autonomous with ξ as the independent variable. In these coordinates, the metric g takes the form, with g tt (ξ) = −K 2 . Note that r(ξ) is necessarily monotonic [50]. Observe also that this chart breaks down at points where K or ∇ µ r vanish, or when they become collinear. Turning our attention to the f -sector, we begin by noting that the quadratic spaces of g and f (which are locally built on a common tangent space), are mapped to each other by the square root S (a linear transformation). For any two non-vanishing tangent space vectors U and V , due to fact that S is a square root f (U, V ) = g (U, S(S(V )), we have g (S(U ), V ) = g (U, S(V )) and f (U, V ) = g (S(U ), S(V )), that is, S is symmetric (self-adjoint) relative to g and it acts as a local pull-back of g to f . The self-adjointness of S follows from Corollary 1.34 in [51].
The equivalent expressions of the above statements in matrix notation are: f = gS 2 (S is the square root), (gS) T = S T g = gS (S is self-adjoint), and f = S T gS (S is a congruence). If S and g are both diagonal, then f is diagonal as well, and the dynamics of f can be easily transferred to S by defining, where τ (ξ), Σ(ξ) and R(ξ) are real fields. The set of independent ± n in (8) indicates the square root branch (giving 16 branches in total for n = 1, ..., 4). The branch will be selected later in section III C.
In the above, we required a non-vanishing τ (ξ) component of S. As we shall see, a reasonable requirement is that all the fields of S in (8) are non-vanishing as they contribute to the determinant of S. There are two reasons for that.
First, the square root is a multi-valued function and the branch of the square root is defined by the 'sign' variables in (8). Since different branches are different functions, changing the selected branch requires manual 'gluing' of the solutions. Note that the selection of the branch is a discrete variable which cannot come out as the solution of a differential equation.
Second, to assert if the obtained solutions are healthy, we can track the elementary symmetric polynomials of S, all of which are scalar invariants and independent of the coordinate system used. In particular, the determinant of S vanishes if and only if S is singular, and a non-invertible S will introduce curvature singularities. This is demonstrated by the following proposition. Proposition 1. Let g be non-singular with finite principal scalar invariants at some point in a region of spacetime, i.e., without curvature singularities. Let further S be a linear transformation which defines f (·, ·) = g(S(·), S(·)). Then a necessary condition for the principal scalar invariants of f to be finite is that S is invertible.
Proof. The proof is coordinate independent. A non-invertible S induces a curvature singularity in f since the index of f is reduced as the Lorentzian signature of g is spoiled by S (i.e., the singular S will create a non-empty kernel of f so that f cannot be inverted). The principal invariants of f (e.g., the Kretschmann scalar) cannot be constructed, since an invertible f is required for contractions. More precise, the connection is not unique for such f and we get a degenerate manifold with singular semi-Riemannian geometry [52,53]. Note that f and g are interchangeable.
As a corollary, this proposition obviously holds for the square root in the HR action (1).
One could think that in vielbein formulation, a continuous transition of the vielbein at the point where det S = 0 would not be a problem. However, the vielbeins are symmetrised by a fixed Lorentz transformation, and the square root branches correspond to reflections of such transformation (where the principal branch of the square root corresponds to the proper orthochronous Lorentz transformation). This can be seen, e.g., from eqs. (2.4), (2.7) and (2.16) in [33], after supplementing (2.7) with the symmetrising Lorentz transformation.
Note that the existence of a healthy connection for f and an invertible S is mandatory for writing the Bianchi constraints in the f -sector (6b). When the field equations are studied as a Cauchy problem, the Bianchi constraints will become the energy and the momentum constraint equations on a spacelike hypersurface, which will be propagated by the evolution equations [54]. Then during development, without the Bianchi constraint, one cannot cross the hypersurface where det S vanishes. Therefore, we postulate an invertible S by definition. However, even if det S is non-vanishing, some other elementary symmetric polynomial of S could be infinite. For instance, τ (ξ)Σ(ξ) could be finite with τ (ξ) → ∞ and Σ(ξ) → 0. In this case, the trace of S will diverge. Hence, we require that all the elementary symmetric polynomials of S, which are in the action (1) and the equations of motion (4), are finite. This imposes, τ (ξ) = 0, Σ(ξ) = 0 and R(ξ) = 0 for all ξ.
B. Choice of the metric fields and behaviour at the Killing horizon The fields g tt (ξ) and g ξξ (ξ) in (7) can be denoted as q(ξ) and F (ξ), so that F (ξ) = 0 corresponds to the Killing horizon, retaining the finiteness of the metric determinant at the Killing horizon (this will be suitable when going to Eddington-Finkelstein coordinates later). The ansatz becomes, with the induced, From now on, we employ the following naming conventions: • The fields at an arbitrary point ξ are denoted just by their name. Example: r := r(ξ), q := q(ξ), R := R(ξ) etc. • The derivatives of the fields (denoted with primes) are always with respect to ξ.
• The fields at some designated point, ξ = ξ label , inherit the same 'label'.
As earlier pointed out, the coordinate system (t, ξ, θ, φ) is not suitable to go through the Killing horizon where K 2 vanishes. For this purpose we introduce the ingoing Eddington-Finkelstein coordinates x µ = (v, ξ, θ, φ) adapted for g µν . In these coordinates, the metrics g and f are (see section B for more details), with the square root, The common Killing horizon is at F = 0. Clearly, g is regular at the Killing horizon in this coordinate system (which then becomes a null frame for g). We have, however, a possible problem with f for which the component f ξξ has F in denominator. We thus make the following proposition: Proposition 2 (Crossing condition).
Consider the metrics g, f and the square root S given by eqs. (13) and (14), wherein g is regular. The necessary and sufficient condition for f to smoothly cross the common Killing horizon F (ξ H ) = 0 is that the following limit is finite, Proof. The proof simply follows from (13). 2 We refer to the condition (15) as the crossing condition. 3 The meaning of 'smoothly cross' in the proposition will become clear later when discussing the relationship between the null cones of the two metrics in section IV C. From now on, we assume that the crossing condition is satisfied. In the limit ξ → ξ H when F (ξ) → 0, we have, 2 The form of S does not influence the result. Consider the 2 × 2 blocks of g and S in the form, The symmetrisation f = f T implies a = (Σ − τ )/F without any restriction on b. Then (with b = 0 in our case), Now, denoting non-zero components by •, two real symmetric matrices A and B cannot be simultaneously diagonalised if they have the following structure, This follows from the theorem on canonical pair forms (see [55,56] and the references therein). Hence, as a corollary, the following proposition holds: Proposition 3 (Proper bidiagonality condition).
Let g and f be the metrics of Proposition 2. Then g and f are simultaneously diagonalisable at the common Killing horizon F (ξ H ) = 0, if and only if, We refer to the condition (18) as the proper bidiagonality condition. A similar statement can be made by extending Proposition 1 in [44] (see section C for more details).
The other way to state Proposition 3 would be that the two metrics are simultaneously diagonalisable at the Killing horizon, if and only if, the Eddington-Finkelstein coordinates for g µν are the Eddington-Finkelstein coordinates for f µν too, and they are a null frame for both metrics. This is obvious from (16) after setting f ξξ = 0. In such case, having also g vv = f vv = 0, the (v, ξ) coordinates become null for both metrics.
Proposition 3 contains two pieces of information: The first is that, if (18) and consequently Σ 2 H = τ 2 H holds, then the two metrics are conformally related separately in their time-radial and their angular part at the Killing horizon. The conformal factor for the time-radial part is τ 2 H and for the angular part, R 2 H . We refer to the time-radial 2-block of the metric with g TR and to the angular 2-block with g ANG . At the Killing horizon, hence, the metrics are block-conformal with conformal factors τ 2 H and R 2 H . The second piece is that, even if we impose bidiagonality outside and/or inside the horizon, this is not enough to guarantee bidiagonality at the horizon. If the proper bidiagonality condition (18) is not satisfied whereas the crossing condition (15) is satisfied, the two metrics are bidiagonal outside and inside the horizon, but not exactly at the horizon. This is due to the fact that, at the horizon, the metrics written in Schwarzschild coordinates are not well-defined, and one cannot apply the inverse coordinate transformation from Eddington-Finkelstein to Schwarzschild to recover bidiagonality. 4 An example of coordinates that diagonalises the metrics, if the proper bidiagonality condition is true, is, T := (v + r)/ √ 2, X := (−v + r)/ √ 2 (which is the π/4 clockwise rotation of the two basis vectors ∂ v and ∂ ξ ). In these coordinates, we obtain the following two block-conformal metrics at horizon, If we use the same coordinates without imposing the proper bidiagonality condition, f will be manifestly non-diagonal with a non-zero dT dX component. This different behaviour at the horizon suggests a new classification for the solutions: (i) proper bidiagonal solutions, for which the proper bidiagonality condition holds, and (ii) improper bidiagonal solutions, for which only the crossing condition holds. The improper bidiagonal solutions are bidiagonal everywhere except at the Killing horizon.

C. Field equations
To expand the field equations (4), we need to specify the branch of the square root. We select the principal branch where ± n in (13) are all positive, and drop the absolute values in (4) assuming all the fields of S to be positive. The other branches will be treated elsewhere. With this choice of the branch, we further define an auxiliary field called the horizon crossing function, According to eq. (15), Φ(ξ) must be finite at the horizon. For increased readability, we also introduce the auxiliary field σ(ξ) by, Note that σ is not to be considered an independent field since its dynamics is governed by Σ. Also, given that rR must be monotonic [50] and Σ must be non-zero and non-diverging due to the regularity of the elementary symmetric polynomials of S, σ must also be non-zero and finite.
With the fields given in (10)- (12), gauge-fixing r(ξ) = ξ, and restricting ourselves to the principal square root branch, the field equations are, Equation (23e) is the Bianchi constraint for g. The Bianchi constraint for f has the same structure as of g, but with the S-fields in the denominator (i.e., divided by det S). Here, we have introduced the utility functions b n [R] and α n [R] defined as, which come out from the interaction between the two metrics; they contain the β-parameters and are always multiplied by the energy scale m of the interaction potential. 5 Since all the fields in (23) are only functions of ξ and there is no explicit dependence on ξ in the equations, the system is an autonomous system of nonlinear ODEs with ξ being the independent variable. Note that the field derivatives in (23) are grouped with r and that r is always given together with m. This shows that the system exhibits a scaling invariance with respect to r, provided that mr is kept unchanged. 6 By simply replacing r → 1 and m → m er, the system will be given in terms ofr = log r. This scaling invariance provides an intriguing possibility to trade r for m and consider m as the dynamical variable.
The system (23) is linear in the field derivatives and in an almost normal form, i.e., the form where first derivatives of each field depend only on the fields and not on their derivatives, with the dependence given as explicit functions. However, solving (23) for derivatives, we get R = 0 which implies that the system is dynamically overdetermined. In other words, there is an algebraic constraint rendering one of the fields non-dynamical. In principle, any of the fields can be rendered algebraic from any two equations, but the inspection of (23) shows that only τ and F have a linear dependence. From eqs. (23d) and (23e), we get, Taking the derivative of either τ or F from (25), and plugging it back into (23), enables us to solve for R (or Σ). 7 Although τ and F could be used on equal footing, the remnant equations are simpler for τ . Thus, we select F to be dynamical and τ to be determined from the other fields by (25a). Note that if we selected F to be determined algebraically, all dynamics would be determined by the square root, S. Once we know the algebraic expression for R , we also know Σ and Φ algebraically (again noting that σ took over dynamics from Σ). Hence, we got the normal form of (23). We also note that equations (23b)-(23f) do not depend on q. This reflects an invariance under the rescaling of t coordinate, e q dt → dt. As a consequence, we can integrate eq. (23a) and determine q after the other fields are known (up to an integration constant).
To summarise: τ , Σ and q can be determined once the fields r, F , σ and R are known, and for the latter we have an autonomous system of four nonlinear ODEs, where three of them are coupled, Here, we introduced two vector fields, V and F, to simplify expressions (these vector fields will be extensively used later on in the paper). The fields F F and F σ are given by (23b) and (23c), respectively. The expression for R = F R will not be quoted here as its total number of indivisible sub-expressions is around 2600.

D. Exact initial conditions at the Killing horizon
The system (26) is an initial value problem. Thus, to solve it, either analytically or numerically, we must specify four initial values (conditions). 8 With the choice of fields (10)-(12), we can impose the exact initial conditions at some common Killing horizon and integrate both outwards and inwards. Therefore, we can obtain the complete solutions in the radial interval r ∈ (0, +∞).
Denoting one of the Killing horizons ξ H , we can start integrating from some r H = r(ξ H ). 9 This is our first of four initial conditions. The second exact initial condition at the horizon is, of course, F (ξ H ) = 0.
We use the crossing condition in Proposition 2 to determine one of the two remaining initial conditions. In order to do this, we study the behaviour of the crossing function, at the horizon. The algebraic expression for τ is known and, as pointed out earlier, it is possible to derive an equation for the combination R + rR , which is inside Σ. It turns out that both τ and Σ depend on F as rational functions, where P indicates the polynomial in F at the numerator and Q the polynomial in F at the denominator. Therefore, we can write, 7 Following a procedure similar to that described in [33].
The coefficients p i and q j are functions of R, σ and r. In order to have a finite limit when F → 0, the term p 0 must vanish at the horizon, After fixing r H , eq. (30) is fourth-order in R H and third-order in σ H and we can solve it for either. We consider R H as a free parameter, and solve eq. (30) for σ H given a value of R H . It turns out that there are always only two independent non-zero solutions for σ H , corresponding to two possible disconnected branches of solutions.
The last thing to consider is the fact that the last term in eq. (29), should have a finite limit, Therefore, after having solved for σ H , we must check that this quantity is finite. If it is, then the crossing condition is satisfied and we can have a BH solution. Otherwise, we exclude these initial conditions. Note also that, since To summarise, besides r H , we only have one free parameter, R H . Once specified, and after choosing the branch for σ H , the system is completely determined and can, in principle, be integrated. However, the equations cannot be integrated analytically, in general. Nonetheless, in the next sections, we will describe a strategy allowing us to do analytical considerations.

IV. BLACK HOLE SOLUTIONS
In this section, we discuss the BH solutions. Subsection IV A is about the choice of global parameters, i.e., parameters appearing in the action (1). In section IV B, we describe the detailed analysis of found BH solutions. In particular, in subsection IV B 4, we apply the analysis to the model defined by the chosen global parameters and show that the only solutions converging to Minkowski, AdS and dS at large radii are the GR solutions. In section IV C, we study the causal structure of improper bidiagonal solutions both inside and outside the Killing horizon.

A. Global parameters
From now on, we use geometrical units, G = c = 1, together with the Planck constant, h = 1. We measure lengths in units for which the Fierz-Pauli mass [14,15] is m g = 10 −1 . To facilitate comparison with the results of [33], we use the default set of global parameters given in Table I. We have considered also many other parameter values, but the qualitative behaviours of the found solutions do not depend on their specific values.
Considering the coordinate chart x µ = (t, r, θ, φ), we demand that the bi-Minkowski solution, is a solution of eq. (26), and obtain the so-called asymptotic flatness conditions, 10 These conditions can also be deduced from cosmological solutions [57]. In addition, we use the normalisation condition, The normalisation (34) can always be assumed without loss of generality, since the action (1) is invariant under the imposing the normalisation condition together with the rescaling m 2 → m 2 δ. By using the normalisation condition, we effectively trade one of the free parameters β 1 , β 2 , β 3 in favour of m g . Usually, the asymptotic flatness conditions and the normalisation condition contain the conformal factor, often denoted by c, between the two metrics for the bi-Minkowski solution. However, any such c can be absorbed into the β parameters by rescaling β n → β n c n . By definition, in the 4-dimensional phase space (r, F, σ, R) of eq. (26), exact initial conditions at a Killing horizon must belong to the 3-dimensional hyperplane F = 0. Moreover, eq. (30) defines a 2-dimensional surface A embedded in the hyperplane F = 0. To be more precise, eq. (30) defines an algebraic variety [58], i.e., the set of solutions of a system of polynomial equations, in our case eq. (30). The surface A is shown in Figure 1a. The red and blue regions represent the two independent branches for σ H (which are the solutions of eq. (30)), as functions of R H and r H . Setting the Killing horizon to some r H = r * , defines a 1-dimensional curve C ⊂ A, as the intersection of the plane r H = r * and A (see Figures 1a and 1b, where r * = 1). 11 In Figure  an improper bidiagonal solution close to SAdS. The SdS solution is not considered in Figure 1 since it has R H < 0 for this parameter set. 12 In section IV B 2, we will explain in more detail why GR solutions correspond to these initial conditions.
indicates a proper bidiagonal solution which is non-GR (we come back to this solution in section IV B 3).
Summarising, for a BH to have a Killing horizon at some r * , the metric fields must assume values lying on C at r = r * . Starting from initial conditions belonging to A, we can obtain all possible static and spherically symmetric proper and improper bidiagonal BH solutions having a common Killing horizon. Note that only values of R H belonging to the domain of C are allowed (that is, for which eq. (30) has real solutions).
We choose the Killing horizon to be at r H = 1. Given that m g = 10 −1 , the Compton wavelength of the massive mode is λ = 10, and r H /λ = 10 −1 . The Fierz-Pauli mass is observationally constrained to be very small, m g 10 −30 − 10 −33 eV [59]. Therefore, the ratio between the Killing horizon and the Compton wavelength is expected to be much smaller than 10 −1 . Such small numbers are demanding to handle numerically. However, since the algebraic variety is only shrunk when r H becomes smaller (see Figure 1a), there are no qualitative changes in the structure of initial conditions. For this reason, we can reasonably extend the following results to more realistic BHs.

Lyapunov stability and the conjugate system of ODEs
Lyapunov stability is an important property of a solution to an autonomous system of ODEs. It is defined as follows [60]: Consider a generic autonomous system of ODEs with independent variable t and let S denote the set of all its solutions. A solution x(t) ∈ S is called stable if, for every > 0, there exist a δ > 0 and a t 0 such that the following holds, If, in addition, there exist a δ > 0 and at such that the following holds, then x(t) is called asymptotically stable.
One can consider both linear and nonlinear Lyapunov stabilities. In the former case, the equations of motion are linearised and stability properties refer to the linearised system, whereas in the latter case, equations are left in the original nonlinear form and stability properties refer to the full system. If not stated explicitly, we will be referring to nonlinear Lyapunov stability of generic solutions of eq. (26) simply as stability. Note that the independent variable t does not need to be time. Thus, this concept of stability does not necessarily concern the time evolution of the system.
When dealing with ODEs and partial differential equations in GR, the study of the stability of the solutions usually refers to Lyapunov stability [40, sections 10.2-10.5]. Indeed, this concept is very useful physically; for example, when considering linear perturbation around a fixed background, one typically expands the metric as the sum of the background metric and the linear perturbation. This is equivalent to study the absolute difference between the unperturbed and the (linearly) perturbed solution, as one has to do in studying (linear) Lyapunov stability. Our approach is quite similar, although it concerns nonlinear Lyapunov stability.
In general, when studying a system of ODEs, it is possible to change variables in order to simplify the equations or to get additional insights. In doing this, one applies a C k -diffeomorphism to the phase-space, with k at least 1. The resulting system is said to be conjugate to the original one and shares with the original system most of its qualitative properties, such as the number of fixed points [61]. We remind that a fixed point of a system of ODEs is a solution with all derivatives vanishing identically; in our case, this is when V (ξ) ≡ 0 in eq. (26). On the other hand, in order to study the stability of a generic solution x 0 (t), one can use the following is another generic solution of the system of ODEs [60]. Since the new system is conjugate to the old one, one can study the former in order to obtain information about the latter without loss of information. Clearly, if |x(t) − x 0 (t)| satisfies the requirements in the definition of Lyapunov stability, then x 0 (t) is a stable or asymptotically stable solution.
We now apply these concepts to our case. Consider the C ∞ -diffeomorphism given by, defined in r ∈ (0, +∞) and where Λ ∈ R is constant. Here, {ρ,F ,σ,R} are the variables of the conjugate system of ODEs. In doing this, we are: (i) shifting the fields with respect to GR solutions, and (ii) introducing the inverse of the radial coordinate in order to map the radial infinity to ρ = 0. 13 Plugging the shifted functions into eqs. (23b), (23c), (23f) and into the equation for R in eq. (26), we obtain the conjugate ODEs system, where all the fields are functions of ξ. Here, Σ is a function ofFR(ρ,F ,σ,R) rather thanR . Also, α n and b n now stand for α n [R + R 0 ] and b n [R + R 0 ], respectively. To find the necessary conditions that must be imposed in order to have fixed points of the conjugate system, we assume that F ,σ,R, ρ ≡ F 0 ,σ 0 ,R 0 , ρ 0 , with the subscript 0 indicating a constant value. IfF 0 = 0, then F is not a solution to the original system. Given that solutions of conjugate systems are uniquely mapped into each other through the diffeomorphism [61], we conclude that we can have a fixed point of the conjugate system only ifF 0 = 0. Besides, ρ ≡ 0 if and only if ρ ≡ 0, corresponding to radial infinity. Our fixed points are thus located at radial infinity and their study give information about the asymptotic structure of the solutions. After inserting F ,σ,R, ρ ≡ 0,σ 0 ,R 0 , 0 into eq. (38), the derivatives of the fields are identically zero only if the following conditions hold (otherwise the right-hand sides of the equations diverge, having ρ in the denominators), where R 0 +R 0 = R H from (37). In particular, condition (a) makesF vanish identically, whereas conditions (b) and (c) combined makeσ vanish identically, and condition (d) is equivalent toR = 0 identically. By imposing these necessary conditions, we are considering all the fixed points of the conjugate system of the form, Keeping ρ arbitrary and substituting ρ,F ,σ,R = ρ, 0, 0,R 0 in eq. (37), together with the conditions in eq. (39), we obtain (coming back to the original system), 13 Note that ρ = 0 is outside of the original domain, but it can be included by compactification.
Finally, using eqs. (25a) and (41), we obtain the following quartic equation in R which yields all possible R H , with the resulting (see eq. (39a)), If we require that R H = 1 is a flat solution, i.e., that 0 = Λ(1) , we get the asymptotic flatness conditions (33). Conversely, imposing (33) ensures that R H = 1 is always an element in the set (42). Now, the equations in (41) tells us that the fixed points correspond to the points at infinity reached by GR solutions when r → ∞. Second, this shows that GR solutions are exact solutions of eq. (26). Note that the cosmological constant of the GR solutions depend on R H , the free parameter in our initial conditions. Therefore, depending on R H , we will have Schwarzschild, SdS or SAdS solutions (see Figure 1b). The asymptotic structure is thus determined by the initial conditions at the Killing horizon. For example, for R H = 1, Φ ≡ 0 and Λ(R H ) = 0, i.e., the Schwarzschild solution. Note again that the asymptotic flatness conditions (33) guarantee that R H = 1 is always a possible initial value for R(ξ). Substituting the fields in eq. (41) into the metrics, the metrics become conformal, with conformal factor R 2 H . These results are all compatible with previous results (see, for example, [33]). After the considerations made in this subsection, we can better understand Figure 1, in which GR solutions have σ H = 1 and Φ H = 0.
We note that, Λ being determined by R H through eq. (43), all values of R H for which Λ(R H ) = 0 give us Schwarzschild solutions. Therefore, it is possible to find specific sets of global parameters for which more than one value of R H has real solutions for σ H and Λ(R H ) = 0. In such cases, the Schwarzschild solutions will differ in the value of the conformal factor between the two metrics. Of course, this can also be the case for SAdS and SdS solutions, with each solution having a different value of the cosmological constant.
Solutions are thus not totally determined by the size of their horizon (i.e., by their mass), but also by R H . In [33] such cases were already found and denoted "special" solutions. We show that the values of the global parameters, determine whether multiple Schwarzschild, SAdS or SdS solutions exist or not. In all cases, every possible solution can be found by integrating from the algebraic variety at the Killing horizon.
The constraint Φ H = 0 defines a finite number of points on C, hence a finite number of proper bidiagonal solutions (we saw three of them in Figure 1b). Proper bidiagonality does not imply GR solutions, as we can see in Figure 1b. The solution indicated by has Φ H = 0, σ H = 1, R H = 3/2 but R H = 0 at the Killing horizon; therefore, it is not a fixed point (i.e., a GR solution). However, we can prove the following, Assuming Φ ≡ 0 implies Σ ≡ τ and Substituting σ and τ into Φ by using the equations of motion, and substituting τ ≡ Σ = σ (rR + R), we obtain, where all the metric fields are functions of ξ. By hypothesis, Φ is identically zero. The second term depends on β parameters which are inside the utility functions α 1 , b 1 and b 2 , whereas the first term does not depend on them. Therefore, in order for Φ to be identically zero, independently of the β parameters, R must be constant (i.e., R ≡ R H ) and σ ≡ 1. This implies Σ ≡ τ ≡ R ≡ R H , i.e., g µν and f µν are conformal (as we can see in (13)), corresponding to GR solutions. The classification of bidiagonal static and spherically symmetric BH solutions is shown in Figure 2.

Lyapunov instability of GR solutions
In this section, we show that GR solutions are nonlinearly Lyapunov unstable. We will do this using both numerical investigations and analytical considerations.
We start with the numerical solutions (see section D for more details about the accuracy of these solutions). Using the parameter values given in Table I, the set of the fixed points (42) gets the values, This is consistent with [33]. When used as the initial values, these correspond to GR solutions. First, note that the value R H = 1 is exact, guaranteed by the asymptotic flatness conditions in eq. (33). Second, there are two negative values for R H . For GR solutions, R ≡ R H ; For non-GR solutions, instead, given that R cannot change sign (since this implies det S = 0), if R H is negative, then R will be negative for all ξ (i.e., all r). The same holds for σ which is positive (because σ H = 1). Then, from eq. (14), we can deduce that different signs of R H mean different disconnected square root branches, belonging to different model and different phase-spaces. Note that these square root branches are different from the two branches of solutions obtained by choosing the value of σ H . A negative R is equivalent to retaining the absolute value of R in eq. (14) together with ± 3 = ± 4 = −1. We emphasize that the equations of Although the solution shows a behaviour ∝ r 2 , it is different from SAdS, which is analysed below in the main text. Note that F pert exhibits a singularity in the σ field (not shown here) at finite r, therefore this is not a physical solution. The vertical solid black line indicates the Killing horizon; the vertical dashed black line indicates the Compton wavelength of the massive mode. The y axis is in sign(y) log 10 (1 + |y|) scale, which is ≈ y for small y and ≈ sign(y) log 10 (|y|) for large y.
motions (23) are valid only for R > 0. However, the transformation, is equivalent to changing the square root branch, but for different β parameters and different asymptotic flatness conditions (i.e., in a different phase space).
Choosing R H = R Schw H ≡ 1, we obtain the Schwarzschild solution. Perturbations around it can be obtained by considering R H = R Schw H (1 + ) and selecting the same branch for σ H as the Schwarzschild solution ( Figure 1b). Hence, these perturbations are improper bidiagonal (like the perturbations around SAdS and SdS, treated later). All found perturbations around Schwarzschild solutions have a diverging σ(r) at some finite radial coordinate r =r; the bigger , the smallerr. A diverging σ(r) means that det(S) is diverging, i.e., det(S −1 ) = 0. Therefore, these perturbations are not physical solutions. Nonetheless, we can introduce the change of variable σ(r) → s(r) := σ(r) −1 and solve for s(r), which is 0 when σ(r) diverges. This strategy allows us to study the solution asymptotically using numerical integration, despite the fact that perturbations are not physical. Figure 3 shows the difference between the analytical Schwarzschild solution, with R H = 1 and the numerical perturbation around it, F pert , with R H = 1.03. The fact that the difference is increasing with growing r, is true for all found perturbations around Schwarzschild solutions, regardless of the values of the global parameters. This shows that the Schwarzschild solution is Lyapunov unstable. Contrary to [34] where asymptotically flat solutions other than the Schwarzschild were found, our result shows that Schwarzschild solution is the unique asymptotically flat solution (see also section A). The numerical integration in Figure 3 extends to r = 10 16 . However, we will show later that there is strong evidence that the instability of Schwarzschild (and, more generally, of GR solutions) persists in the limit r → ∞, independent of the global parameter values. The y axis is in sign(y) log 10 (1 + |y|) scale, which is ≈ y for small y and ≈ sign(y) log 10 (|y|) for large y. However, the study of the difference F pert − F SAdS at large radii, which matters for Lyapunov stability, shows that it is unstable. Here, F SAdS indicates the analytic SAdS solution, which is 0 when r = r H . The difference between the improper bidiagonal perturbation around SAdS with R H = 3 and the analytic SAdS solution with R H = R SAdS H is plotted in Figure 4. F pert − F SAdS oscillates with increasing amplitude at radii larger than the Compton wavelength of the massive mode (λ = 10 in our case). A similar asymptotic behaviour around SAdS was found in [62].
Such a result shows that SAdS is an unstable solution. We emphasise that all the found perturbations around SAdS have a non-diverging σ(r) and a monotonically increasing R(r)r (which is a necessary condition to have a spherically symmetric spacetime [50]). Therefore, they do not suffer from any pathological behaviour and are physically acceptable solutions.
We now turn to the analysis of SdS solution, with R H = R SdS 0 ≈ −8.557 (as pointed out, belonging to another model). Since the SdS solution has both a Killing and a cosmological horizon, we can perturb it at either horizon. When perturbing the solution at the Killing horizon, one could expect it to have a cosmological horizon close to the corresponding SdS cosmological horizon. However, all the found perturbations around the Killing horizon displayed a diverging σ(r), before reaching a cosmological horizon. This means that they are not physically acceptable. These perturbations, defined on a compact domain, were found already in [33], where it was stated that there are no other perturbations to SdS solutions. However, solutions obtained when perturbing SdS at the cosmological horizon do not have a compact domain, since they are defined when r → ∞. On the other hand, integrating inwards from the cosmological horizon, σ(r) diverges before reaching a Killing horizon for all the found perturbations, showing that also these solutions are non-physical. Nevertheless, in order to study the asymptotic behaviour of perturbations (1 + ), with = 10 −3 . F pert − F SdS oscillates with a growing amplitude, showing that the SdS solution is Lyapunov unstable. The envelope of the oscillations is proportional to r 1/2 . Note that F pert exhibits a singularity in the σ field (not shown here) at finite r. Therefore, this is not a physical solution. The vertical dashed black line indicates the Compton wavelength of the massive mode. The y axis is in sign(y) log 10 (1 + |y|) scale, which is ≈ y for small y and ≈ sign(y) log 10 (|y|) for large y.
around SdS solution, we integrate numerically from the cosmological horizon, imposing initial conditions lying on the algebraic variety. Note that R C = R H for the SdS solution (because R is constant), of course corresponding to different r H and r C . The SdS solution with Killing horizon r H = 1 has the cosmological horizon at r = r C ≈ 3.128. In Figure 5 Figure 4, we conclude that also the SdS solution is Lyapunov unstable. Apart from the fact that amplitudes always grow with increasing radii, numerical integration shows that the amplitude of the oscillations is proportional to the size of the perturbation around R H .
We will now show how the numerically found Lyapunov instabilities of GR solutions can be motivated by looking at the equations of motion in eq. (38). Thanks to the autonomy of the system, we can study the 3-dimensional phase space section defined by ρ = 0, corresponding to infinite radii. When ρ = 0, the right-hand sides of the first two equations in eq. (38) are finite if the conditions in eq. (39) hold. That is, if we have a fixed point, corresponding to a GR solution.
For differentiability of the solutions, a non-GR solution will thus necessarily have a diverging gradient for infinitesimally small ρ, i.e., infinitely large r. Therefore it will not tend to a fixed point, since for this to be possible, the modulus of its gradient in the phase space should tend to 0. Hence, the non-GR solution will depart more and more from the GR solutions when r → ∞, strongly indicating that GR solutions are Lyapunov unstable, independently of the global parameter values.
Since we have not found a Lyapunov instability function, we do not have a formal proof of the instability of GR solutions. Nevertheless, the reasoning above combined with the empirical fact that all numerical solutions we have found show the instability, and the phase space analysis that we are going to pursue, are strong hints that GR solutions are Lyapunov unstable.
In the following, we show this instability by directly plotting their phase space portraits. Before doing that, we introduce our method by showing GR phase space portraits. In GR, for a static and spherically symmetric metric, Einstein equations reduce to, which can be formally rewritten as an autonomous system of two ODEs, Following the procedure described for example in [63], we introduce variables which compactify the phase space, Since ≥ 0 and 2 + ϕ 2 ≤ 1, with the last equality valid only in each of the limits r → ∞ and F → ∞, the whole phase space is compactified onto the unitary right half disk by stereographic projection. Depending on the value of Λ, we will have different phase spaces, plotted in Figure 6. In each of the three phase spaces, the family of trajectories is parametrised by r H (or, equivalently, by the mass of the BH), which is the only parameter we can vary to obtain a new solution in GR. All trajectories will eventually reach asymptotically stable fixed points, indicated by black dots in Figure 6, displaying the Lyapunov stability of the solutions. In HR bimetric theory, varying only r H keeping R H fixed and equal to one of the values in eq. (48), we will obtain a family of Schwarzschild solutions, SAdS solutions or SdS solutions (depending on the constant value of R H ). The solutions belonging to the same family will asymptotically approach each other, exactly as in GR.
On the other hand, varying the ratio between the lengths of the common Killing horizon with respect to f µν and g µν , i.e., R H , does change the asymptotic behaviour of solutions belonging to the same phase space (i.e., to the same model), as we saw explicitly in numerical solutions. This fact can be appreciated from the bimetric phase spaces shown in Figure 7, where the vector field in eq. (38) is plotted. That is, we can infer the behaviour of the solutions without the need for numerical integration. Note that the considered Schwarzschild and SAdS solutions belong to the same phase space, whereas the SdS solutions belongs to a separate phase space, as we already highlighted.
The ODE phase space is 4-dimensional, with variables ρ,F ,σ,R . Using a color code for one dimension, we can make 3-dimensional plots. Here, the ρ axis is represented by the color code. In Figure 7, we take 2-dimensional sections of the 3-dimensional plots in order to emphasise the features we are interested in. Each plot shows the phase space around one GR solution, indicated by a red sphere located at (0, 0, 0).
In the first row, we show the phase space around the Schwarzschild solution, in the second row we show the phase space around the SAdS solution, in the third row we show the phase space around the SdS solution. The first column contains the R ,σ projections in which we are looking at the phase space "from the top", the second column contains the R ,F projections in which we are looking the phase space "from the front".
The value of ρ = 1/r is indicated by the colour of the trajectories: the bluer the trajectory, the larger the ρ (smaller r), the redder the trajectory, the smaller the ρ (larger r). The asymptotic behaviour of the solutions is thus described by the red regions. In each plot, one arbitrary reference solution is highlighted in black.
From Figures 7a and 7b, we see that the Schwarzschild solution is Lyapunov unstable. In Figure 7a we clearly see the saddle point structure of the phase space. There is one stable direction (left-bottom to right-top diagonal) and one unstable direction (left-top to right-bottom diagonal). Looking at the black reference solution in Figure 7b, we see that the phase space flow is departing from the Schwarzschild solution also along theF axis when r → ∞. Note also that the plots represent a very small region around the Schwarzschild solution (as evident from the frame scales). Very small perturbations in R H thus lead to completely different asymptotic behaviours.
We comment on the last four plots together, because they share the same qualitative features in pairs. In Figures 7c  and 7e, we see the top views of the phase spaces around SAdS and SdS solutions, respectively. They both look asymptotically stable, since the flow seems to be approaching the point (0, 0, 0) when r → ∞ (i.e., when the flow becomes redder). This is radically different from the Schwarzschild case, where the flow was showing a saddle point structure. However, by looking at the black reference solutions in the front views of Figures 7d and 7f, it is evident that SAdS and SdS solutions are unstable. Although trajectories get closer to (0, 0, 0) in the R ,σ plane, they are oscillating around the SAdS and SdS solutions with a growing amplitude along theF axis. These results are consistent with what we found numerically, i.e., perturbations around Schwarzschild solutions diverge from it and perturbations around SAdS and SdS solutions diverge from them with growing amplitude oscillations.
As we already pointed out, perturbations around the Schwarzschild and SdS solutions are not physically acceptable because det(S −1 ) = 0 for some finite r. Phase space portraits tell us that, even if well-behaved perturbations were found, they would not approach GR solutions; they would follow the unstable flow in the phase space.
Of course, plots in Figure 7 describe trajectories only up to some finite r. However, the system is autonomous and therefore the qualitative behaviour of the flow is not changed when r grows. Although not a formal proof, this is a strong indication of instability.
Finally, we have also studied Lyapunov stability of GR solutions in the dRGT massive gravity limit, κ → ∞, and we found analogous results. However, for some values of the global parameters, the Schwarzschild solution happens to be Lyapunov stable, but not asymptotically stable, which is compatible with the results of [33]. This means that perturbations oscillate around it with constant amplitude, not being asymptotically flat.

C. The causal structure of the improper bidiagonal solutions
In the following, we investigate the causal structure of the black hole solutions by looking at the relationship between the null cones of the two metrics. In particular, we will construct a set of null frame vielbeins and plot the resulting null radial geodesics aligned with the associated generating null cones. This analysis will highlight the difference between the proper and the improper bidiagonal solutions in more detail.
Consider the metric fields g and f at some point covered by the ingoing Eddington-Finkelstein coordinates x µ = (v, r, θ, φ) adapted for g. We can construct a complex null frame field (vielbein) of the metric g which suits the ingoing null radial geodesics as, µ g = e −q/2 , It is easy to verify that g µν = −2 g is the inverse of the metric g from eq. (13). In general, the geodesics of g are obtained by the variational method from the action, where the re-parametrisation invariance is ensured by the einbein e(λ). Moreover, e(λ) acts as a Lagrangian multiplier imposing the mass-shell constraint so that µ 2 = 1, 0, −1 corresponds to time-, null-or space-like geodesics, respectively. For the null radial geodesics, besides µ 2 = 0, we have dθ/dλ = dφ/dλ = 0. Note also that the vielbein (56) solely depends on the radial coordinate x 1 ≡ r. Under these conditions, the variation of (57) with respect to x µ (λ) and e(λ) gives the following two geodetic equations, up to a suitably gauged e(λ), Hence, the above null vector fields are aligned with the tangent vector fields of null radial congruence matching the null radial rays. In other words, the outgoing and the ingoing null radial geodesics of g are the integral curves of the null vectors g and n g , correspondingly.
In similar fashion, a complex null frame vielbein of the metric f that is adapted to the ingoing null radial geodesics can be constructed as, f yields the inverse of the metric f from eq. (13). Then, the resulting outgoing and the ingoing null radial geodesics of f are the integral curves of the null vectors f and n f , respectively, with the corresponding geodetic equations, Given the solution to the equations of motion (23), the null radial geodesics of g and f can be determined by integrating eqs. (58) and (60), displaying the causal structure of the spacetime. Note that the g and f null cone fields are simply plotted as given by eqs. (56) and (59). In order for the null cones to be Minkowski-like in the weak field limit, however, it is customary to use the coordinates (t , r) = (v − r, r). In these coordinates, the null radial geodesics and the null cone fields for the improper bidiagonal solutions described in section IV B, are plotted in Figures 8-10.
Several observations can be made from these plots. First, the static symmetry is obvious; as the fields are dependent only on the radial coordinate r, the geodetic plots are invariant under translations along the t axis.
Second, far from the Killing horizon, but inside the Compton wavelength of the massive mode, the null cones are almost coinciding. While approaching the horizon, they differ more and more. At the horizon, the null cones share at least one null direction (Figures 8b, 9b, 10b, 10c), since the translational Killing vector becomes null for both the metrics. Inside the horizon, both null cones shrink while r → 0, showing that the Killing horizon is actually an event horizon. The ingoing null radial geodesics entering the Killing horizon do not come back; they approach the curvature singularity as in GR. However, starting from the same point, the ingoing null radial geodesics of g and f reach the curvature singularity at different finite t , that is, at different advanced times v.
Finally, the causal structure from the plots displays the difference between the proper and the improper bidiagonal solutions. As pointed out in [56], the algebraic type of the improper bidiagonal solutions correspond to the null cones having an odd number of common null directions with a common tangent plane (Figures 8b, 9b, 10b, 10c). On the other hand, the algebraic type of the proper bidiagonal solutions correspond to the null cones having an even number of common null directions with a common tangent plane (the GR solutions have coinciding null cones, therefore they share all their null directions).
While crossing the Killing horizon, the null cones of the improper bidiagonal solutions deform continuously (as can be seen directly from eqs. (56)-(60)), therefore the fact that they cannot be simultaneously diagonalised only at the horizon is not a singular behaviour. This clarifies the meaning of the expression "smoothly cross" used in Proposition 2 (section III B).

V. CONCLUSIONS
In this paper, we studied static and spherically symmetric bidiagonal black hole solutions in HR bimetric theory. Although the subject in part has been treated previously, this study is motivated by the fact that a nonlinear phase space stability analysis and a causal classification of these solutions was lacking. The former study shows that all the found non-GR solutions either diverge from the GR solutions at large radii, or have a singular square root S = g −1 f at finite r. These results are complementary to those in [33] and extend them. There, the author considers the asymptotic properties of BH solutions in the leading order by studying the relative differences between GR solutions and their perturbations; here, we study the absolute differences between solutions and focus on their convergence properties in the phase space of our system of ODEs, i.e., we study the Lyapunov stability of the solutions, which was not analysed before.
We emphasize that the phase space analysis allowed us to look at all solutions of the fields equations in a qualitative way, for a given set of the global parameters appearing in the action of the theory. Exploring the phase spaces for several different sets of global parameters gives equivalent results. Hence, we conjecture that GR solution are Lyapunov unstable independently on the values of the global parameters of the theory. This conjecture is also motivated by analytical considerations about the behaviour of the fields equations at large radii. In particular, this conjecture implies that there cannot be asymptotically flat bidiagonal solutions other than Schwarzschild.
We pointed out that choosing a bidiagonal ansatz, is equivalent to imposing that the two metrics share their Killing horizons. We also presented an original proposition stating that, provided that one metric is regular, a necessary condition for the other metric not to have curvature singularities is that the matrix square root S = g −1 f is non-singular.
Thanks to an appropriate parametrisation of the metric functions, we have been able to detect pathological features of solutions. It also allowed us to find, for the first time, exact initial conditions for the metric fields at the Killing horizon, using Eddington-Finkelstein coordinates. The initial conditions lie on a 2-dimensional surface in the phase space of our system of ODE. Every possible bidiagonal static and spherically symmetric black hole solution must have the metric fields intersecting this surface at the Killing horizon and, if present, at the cosmological horizon.
We found that the sole BH solutions asymptotically converging to flat, de Sitter or Anti-de Sitter spacetimes have proportional metrics, equivalent to the corresponding GR solutions. Other theoretically consistent solutions exist, but they all diverge from flat, de Sitter or Anti-de Sitter space at large radii. We show that the former correspond to proper bidiagonal solutions, and the latter to improper bidiagonal solutions. The proper bidiagonal solutions are bidiagonal everywhere, whereas the improper bidiagonal solutions are bidiagonal inside and outside the Killing horizon, but not at the Killing horizon. Note, however, that all proper bidiagonal solutions do not necessarily correspond to GR solutions.
Our results show that black holes having very small differences in the metric fields at their Killing horizons, will have completely different asymptotic structures. In the light of the asymptotic exotic behaviour of the non-GR bidiagonal solutions, physical solutions would naïvely correspond to the GR solutions. However, assuming that S is non-singular, i.e., the solutions are theoretically consistent, we note that this behaviour start to become important only at radii similar to the Compton wavelength of the massive mode. This scale is many orders of magnitude outside the Killing horizon, typically even at cosmological length scales, and anyway far outside the region where the assumption of an isolated black hole breaks down. Therefore, in a real physical setting, the divergences may be either non-existing or non-observable.
Only studying static solutions, we can not answer the question of which solutions are realised in the process of gravitational collapse of matter. However, in [35][36][37][38] it was shown that the bidiagonal solutions are dynamically unstable, and as such cannot represent the end point of gravitational collapse. In addition, [37] showed that to have dynamically stable solutions, non-diagonal metric elements are needed, also at radii other then the horizon radius, r H . Therefore, in order to investigate what is the end point of gravitational collapse in HR bimetric theory, a full dynamical treatment of the process probably need to be performed. Also, the Lyapunov stability of non-bidiagonal BHs and electrically charged BHs should be explored.
where Q(r), N (r), a(r), Y (r), U (r) ∈ R. The determinants of these metrics are, The square root matrix S = g −1 f , after having chosen the principal branch, is, The absolute values can be removed assuming that the fields are always positive. The determinant and the trace of the square root are, The Ricci scalars for the two metrics are, where all fields are functions of the radial coordinate r. We are now going to discuss why this parametrisation is not optimal when studying BH solutions. First, this parametrisation is not defined inside the Killing horizon of the BH. Indeed, the norm of the translational Killing vector K µ = δ µ 0 is, In the coordinates adapted to the translational Killing vector, we cannot cross the Killing horizon since K becomes null at r = r H . On the other hand, using the Eddington-Finkelstein coordinates is problematic since K 2 = −Q 2 cannot be positive, so we cannot cross the Killing horizon. Yet, if we use the Eddington-Finkelstein coordinates, this parametrisation is defined in the interval r ∈ [r H , +∞), which, in principle, allows specification of the initial values on the Killing horizon. Nonetheless, the numerical solutions found in the literature (see, for example, [10,30,33,34]) are obtained by imposing initial conditions close to the Killing horizon, after having expanded the equations up to some (usually the first) order around r = r H . Namely, in order to define a Killing horizon at r H , one must impose Q(r H ) = a(r H ) = 0. In order for the determinants of the two metrics not to be zero at r H , we must also have N (r H ) = Y (r H ) = 0, with lim r→rH Q(r)/N (r) = const and lim r→rH a(r)/Y (r) = const. These conditions can be difficult to control when dealing with numerical solutions. Also, pathologies in the solutions are not straightforward to analyse. If only one of the metric fields becomes zero at finite r, the singularity introduced can be of different kinds, depending on the behaviour of more fundamental quantities.
For example, suppose that, at some finite valuer, Y (r) = 0, but the other metric fields and their derivatives are non-zero and finite. In this case det S, det f and f rr diverge, whereas the g-sector is regular. Therefore we have a determinant singularity at r =r, whose effects are studied in [64,65]. Suppose, instead, that for some finite r * , a(r * ) = 0, but the other metric functions and their derivatives are non-zero and finite. Then, det f = 0 and the Ricci scalar R f of f µν diverges, indicating a curvature singularity which is not hidden by a horizon, i.e. a naked singularity. In principle, one could think that pathologies are automatically excluded by the equations of motion, i.e. an actual solution would not have these problems. This is not the case, for example, in numerical solutions. We therefore set out to find a more convenient parametrisation for the metric functions in section III B.
We emphasise that detecting pathologies is equivalent to find the points where det(S) = 0, det(S −1 ) = 0, or where other elementary symmetric polynomials of S are not finite. With our new parametrisation of the metric fields, this corresponds to points with a singular behaviour of Σ(r), τ (r) and R(r), separately.
We now consider an explicit example of a pathological behaviour, already understood as such in [34]. The solution in question is plotted in Figure 11 and is a good example of the importance of Proposition 1. Note that this solution is defined in the interval (r H = 1, +∞) and has β 1 = 1, β 2 = β 3 = 0, β 0 and β 4 determined by the asymptotically flatness conditions in eq. (33) (so c = 1) and m g = 0.04 (m g is denoted with µ in [34]; the Compton wavelength of the massive graviton is then λ = 25 in geometrical units with the Planck constant h = 1). Hence m g m S = m g r H /2 = 0.02, with m S Schwarzschild mass of the BH. The second panel of Figure 1 in [34, p. 4] tells us what is the value of U (r H ) = u H to choose as initial condition, given the value of m g m S . Note that, since r H = 1 both in [34] and in our approach, u H corresponds to our R H . For m g m S = 0.02, u H ≈ 50.2.
In Figure 11, all relevant quantities for this solution are plotted: the metric functions and the first derivative U (r), which appears explicitly in f µν , det S and the Ricci and Kretschmann scalars for f µν . The curvature scalars for g µν diverge at r = 0 only, where the curvature singularity arises, and are not plotted. The third panel of Figure 11 shows det S; the grey solid lines highlight the points at which det S = 0. These lines are plotted also in the other panels, in order to understand what gives rise to the singular behaviour. In the first panel, we see all metric functions, except U (r) and U (r) which are plotted in the second panel. The solid grey lines correspond to the following 4 points: from left to right, (i) U (r) = 0, (ii) a(r) = 0, (iii) U (r) = 0, (iv) a(r) = 0. There are also large-dashed grey lines, highlighting the points at which U (r) = Y (r) = 0, keeping det S finite. Note that the second solid grey line and the first dashed grey line are almost overlapping, yet the points they represent are different. Now, we can look at the fourth panel; the curvature scalars for f are diverging when det S = 0 (whereas the g-sector is regular). This confirms the result of Proposition 1. Note that the change of sign of the fields is anyway not compatible with the selected principal branch of the square root matrix used in the equations of motion.
Having a singular square root S, this solution is not physically acceptable. In accordance with Proposition 1, it has four naked singularities at different finite r. In addition, it is not an asymptotically flat solution. Finally, U (r) is not a monotonic function; which is required to have a spherically symmetric spacetime [50]. We also studied this solution in our new parametrisation, finding the equivalent results with respect to the curvature scalars and det S.

Appendix B: Eddington-Finkelstein coordinates
Here we introduce Eddington-Finkelstein coordinates to avoid the coordinate singularity at the Killing horizon. In the coordinate chart x µ = (t, ξ, θ, φ), the metrics g, f and the square root matrix S are given by, The determinants of the metrics are, which are always non-zero except at r = 0, where the curvature singularity arises. This can be seen by looking at the Ricci and Kretschmann scalars for the g metric, The Ricci and Kretschmann scalars for the f metric are more complicated and there is no need to write them, since Proposition 1 implies that we can have induced curvature singularities in f sectors whenever det(S) = 0 or det(S −1 ) = 0. To avoid the coordinate singularity at the Killing horizon, we introduce the ingoing Eddington-Finkelstein coordinatesx β = (v, ξ, θ, φ) for g defined by the following Jacobian, where x α = (t, ξ, θ, φ) are the old Schwarzschild coordinates. This defines the ingoing null coordinate v for g (also known as 'advanced time'), dv = dt + dξ e q(ξ)/2 F (ξ) .
Applying the coordinate transformation (B5) to the metrics (B1) yields eq. (13). We could equally consider the outgoing null coordinate u; this would not affect any conclusions, however, since the equations are invariant under the transformation u = −v. Also, in principle, we could introduce the null coordinate adapted for f rather than for g.
Note that, contrary to the GR case, we cannot explicitly integrate the differential dv because we do not know analytically the functions q(ξ) and F (ξ). Nonetheless, we can still define the tortoise coordinate in an implicit form, where dξ = dr. Note that ξ * is a monotonic function of ξ separately inside and outside the event horizon (or, for a SdS solution, between the horizons, inside the event horizon and outside the cosmological horizon), Appendix C: On bidiagonality condition in [44] In the following we make Proposition 1 in [44] more precise. Consider a coordinate patch where two static and spherically symmetric metrics take the form, Suppose further that the Killing vector K = ∂ t is null with respect to g at r = r H . Proposition 1 in [44] states that, if both metrics describe smooth geometries and are diagonal at r = r H , the Killing vector K must also be null with respect to f at r = r H . The proof presented in [44] is elegant in its simplicity. It starts from the presuppositions A(r H ) = 0, B(r) = 0, and that the trace Tr(g −1 f ) = J/A+K/C +2r 2 /D is finite by assumption. Then, the individual terms in the trace cannot cancel since they necessarily have the same sign as both diagonal metrics have Lorentzian signature when the translation Killing vector K is timelike for r > r H . Therefore J(r H ) = 0, otherwise J/A would diverge at r = r H . Now, we further evolve the proof. Let us only consider the (t, r)-block of the coordinates, in which case det g = −AC and det f = −JK, which are non-zero by the premise of having smooth geometries. From the assumptions of Proposition 1, besides the trace, the determinant det(g −1 f ) = (J/A)(K/C) is also regular. Hence, because J/A is finite, K/C must also be finite, and they must be both non-zero. To see what happens at r = r H , we introduce the ingoing Eddington-Finkelstein coordinates (v, r) adapted for g, by dt 2 = dv 2 − 2 C/A dvdr + (C/A)dr 2 at r > r H . Assuming Lorentzian signature of the metrics, we have (up to square root signs, not shown for readability), Clearly, the metrics are regular at r = r H in this coordinate system. Since both J/A and K/C are finite and non-zero, it is also c := (J/A)(C/K) and we can write, At r = r H we have A(r H ) = J(r H ) = 0; thus, We can apply the theorem on canonical pair forms [55] on (C5). Thus, the metrics which describe smooth geometries cannot be both diagonal at r = r H unless K − J A C = 0 or equivalently J/A = K/C, that is, unless the (t, r)-blocks of g and f are conformal at r H (this corresponds to the proper bidiagonality condition in Proposition 3). In other words, strictly, Proposition 1 in [44] is stated for the proper bidiagonal metrics, but it is clear that the statement is also valid for the improper bidiagonal case for which 0 = K − J A C < ∞ (crossing condition in Proposition 2), i.e., when the metrics cannot be simultaneously diagonalised at r H ∈ {r | A(r) = 0} but they are bidiagonal elsewhere by assumption. In this appendix, we describe in more detail the numerical procedure used to obtain the solutions in section IV B 4. Moreover, we motivate the accuracy of the numerical solutions by studying their residuals, of which analysis is shown in Figures 12-15.
For solving ODEs numerically, we used Wolfram Mathematica 11 [66]. In particular, NDSolve`LSODA method was used, which is based on LSODA [67,68]. The numerical integration in Mathematica proceeds in two main steps. The first is the actual integration, which returns the values of the independent variable after each step and the values of the solution at each step, defining a grid of numerical points. The second step is to interpolate these points with a continuous function which is the final output of the numerical integration. This interpolation function should approximate the exact analytical solution.
In checking the validity of our numerical solutions, we first checked that our results are independent from the value of the WorkingPrecision parameter, which sets how many digits of precision are maintained in internal calculations. Ranging precision from 20 to 50, approximation and truncation errors did not affect the results. All the results presented in this paper were obtained with WorkingPrecision 50.
In addition, to verify that our numerical solutions are approximating the analytical solutions, we performed a backward error analysis, in agreement with the procedure described in [69,70]. Hence, our analysis was focused on the final interpolated function returned by Mathematica. We used the definition of residuals for a numerical solution as introduced in [69]. Suppose we have a first-order system of ODEs, y (x) = f (x, y(x)), with y, f ∈ R n and n ∈ N number of ODEs. We solve for y(x) numerically, obtaining its approximation u(x). Then, the residuals of the numerical solution u(x) are defined as the defect, r(x) := u (x) − f (x, u(x)). (D2) Clearly, if the numerical solution was equal to the exact solution, u(x) = y(x), this quantity would be identically zero. We then ask whether the equation, is close to the original one or not. This is equivalent to asking if residuals are small compared to f (x, u(x)). In [69,70], it is pointed out that, if u(x) is the approximation of y(x) in a given interval of the independent variable x, then, r(x) = u (x) − y (x) + higher order terms.
We will use eqs. (D2)-(D4) to show that the obtained solutions are accurate. In the following, we present the analysis only for the F metric function, but the same conclusions hold for all the metric fields. First, Figures 12 and 13 show that residuals are extremely small compared to the vector field component F F in eq. (26) for all solutions presented in this paper. This means, according to eq. (D3), that our numerical solutions satisfy equations which are very close to the original equations.
In order to study in more detail the properties of residuals, we consider GR solutions first. Knowing their analytic form, we can directly compute the residuals by using not only eq. (D2), but also eq. (D4). They are plotted in  We can now apply the same analysis to perturbations around GR solutions. The purpose is to show that the results presented in Figures 3-5 are not affected by numerical errors. Suppose that perturbations around GR solutions tend to them for large radii. Then, according to eq. (D4), after some finite r the difference |F num pert − F num GR | must be due only to numerical residuals. Conversely, if this difference is not close to the residuals calculated through eq. (D2), the numerical solution is not approximating an exact GR solution. We can then trust Figures 3-5, since the difference between the perturbations around GR solutions and GR solutions is not due to numerical errors.
In Figure 15, we see that the differences |F num pert − F num GR | and |F num pert − F F [V num pert ]| for perturbations of GR solutions differ by many order of magnitudes over almost the whole range of integration. For r > 10 14 , the two quantities start to be comparable. This strongly indicates that perturbations around GR solutions are always different from the GR solutions, independent of numerical errors. In addition, if we calculate the residuals by using the numerical points calculated by Mathematica before doing the interpolation, they are always many order of magnitude smaller than |F num pert − F num GR | (calculated numerically on the grid), even for r > 10 14 . This shows that the interpolation done by Mathematica is another source of error. Therefore, the fact that the red and blue curves in Figure 15 are approaching each other for r > 10 14 , is mostly due to an interpolation error and without physical significance. This is also suggested by the fact that residuals in Figure 15 do not exhibit the clean oscillatory behaviour of |F num pert − F num GR |.