Classification and asymptotic structure of black holes in bimetric theory

TL;DR

This paper classifies static, spherically symmetric black holes in Hassan-Rosen bimetric theory and finds that only solutions matching GR black holes converge to familiar spacetimes at large distances.

Contribution

It provides a new classification scheme for bidiagonal black holes and demonstrates the asymptotic behavior of solutions in bimetric theory.

Findings

Classified black holes based on horizon behavior.

Only GR-like solutions match standard spacetimes asymptotically.

Identified unique solutions converging to Minkowski, AdS, and dS.

Abstract

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.

Tables

Table 1. Table 1: The default global parameter values used in this paper. The asymptotic flatness condition and normalisation condition (see text below) are used to constrain three β 𝛽 \beta parameters, namely β 0 subscript 𝛽 0 \beta_{0} , β 1 subscript 𝛽 1 \beta_{1} and β 4 subscript 𝛽 4 \beta_{4} .

$m_g$	$\kappa ={\left(\frac{M_f}{M_g}\right)}^2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_4$
${10}^{-1}$	$\tan {\left(1\right)}^{-2}$	$-3+3{\beta }_2+2{\beta }_3$	$1-2{\beta }_2-{\beta }_3$	$-\frac{1}{2}$	$-\frac{2}{5}$	$-1-{\beta }_2-2{\beta }_3$

Equations

$$\mathcal{S}=\int {\mathrm{d}}^{4}x\left[\frac{1}{2}{M}_g^2\sqrt{-\det g}{R}^g+\frac{1}{2}{M}_f^2\sqrt{-\det f}{R}^f-m^2{M}_g^2\sqrt{-\det g}\sum _{n=0}^4{\beta }_n{e}_n\left(S\right)\right],$$

$$\displaystyle e_{0}(S)\coloneqq 1,\qquad e_{n}(S)\coloneqq?{\mbox{$S$}}^{[\mu_{1}}{}_{\mu_{1}}?\cdots?{\mbox{$S$}}^{\mu_{n}]}{}_{\mu_{n}}?,\ n\geq 1,$$

$$e_1(S)=\mathrm{Tr}(S),e_2(S)=\frac{1}{2}\left((\mathrm{Tr}S{)}^2-\mathrm{Tr}(S^2)\right),$$

$$e_3(S)=\frac{1}{6}\left((\mathrm{Tr}S{)}^3-3\mathrm{Tr}(S^2)\mathrm{Tr}(S)+2\mathrm{Tr}(S^3)\right),e_4(S)=\det (S).$$

$$R_{\mu \nu }^g-\frac{1}{2}{g}_{\mu \nu }{R}^g+m^2{g}_{\mu \lambda }\sum _{n=0}^3\left[{\beta }_n(-1{)}^n{{Y_{(n)}}^{\lambda }}_{\nu }\left(S\right)\right]$$

$$R_{\mu \nu }^f-\frac{1}{2}{f}_{\mu \nu }{R}^f+\frac{m^2}{\kappa }{f}_{\mu \lambda }\sum _{n=0}^3\left[{\beta }_{4-n}(-1{)}^n{{Y_{(n)}}^{\lambda }}_{\nu }\left(S^{-1}\right)\right]$$

$$?{\mbox{$Y_{(n)}$}}^{\mu}{}_{\nu}?\left(S\right)\coloneqq\sum_{k=0}^{n}(-1)^{k}\;?{\mbox{$(S^{n-k})$}}^{\mu}{}_{\nu}?\;e_{k}(S),\qquad\kappa\coloneqq\left(\dfrac{M_{f}}{M_{g}}\right)^{2}.$$

$${\nabla }_{\mu }^{(g)}\left\{\sum _{n=0}^3\left[{\beta }_n(-1{)}^n{{Y_{(n)}}^{\mu }}_{\nu }\left(S\right)\right]\right\}$$

$${\nabla }_{\mu }^{(f)}\left\{\sum _{n=0}^3\left[{\beta }_{4-n}(-1{)}^n{{Y_{(n)}}^{\mu }}_{\nu }\left(S^{-1}\right)\right]\right\}$$

$$\mathrm{d}{s}_g^2:=-g_{tt}(\xi )\mathrm{d}{t}^2+g_{\xi \xi }(\xi )\mathrm{d}{\xi }^2+r(\xi {)}^2\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2\right),$$

$$S:=\mathrm{diag}\left({\pm }_1|\tau (\xi )|,{\pm }_2|\Sigma (\xi )|,{\pm }_3|R(\xi )|,{\pm }_4|R(\xi )|\right),$$

$$\mathrm{d}{s}_f^2=-\tau (\xi {)}^2{g}_{tt}(\xi )\mathrm{d}{t}^2+\Sigma (\xi {)}^2{g}_{\xi \xi }(\xi )\mathrm{d}{\xi }^2+R(\xi {)}^{2}r(\xi {)}^2\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2\right).$$

$$\mathrm{d}{s}_g^2:=-{\mathrm{e}}^{q(\xi )}F(\xi )\mathrm{d}{t}^2+F(\xi {)}^{-1}\mathrm{d}{\xi }^2+r(\xi {)}^2\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2\right),$$

$$S:=\mathrm{diag}\left({\pm }_1|\tau (\xi )|,{\pm }_2|\Sigma (\xi )|,{\pm }_3|R(\xi )|,{\pm }_4|R(\xi )|\right),$$

$$\mathrm{d}{s}_f^2=-\tau (\xi {)}^2{\mathrm{e}}^{q(\xi )}F(\xi )\mathrm{d}{t}^2+\Sigma (\xi {)}^{2}F(\xi {)}^{-1}\mathrm{d}{\xi }^2+R(\xi {)}^{2}r(\xi {)}^2\left(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2\right).$$

$$g_{\mu \nu }=\left(\begin{array}{cccc}-{\mathrm{e}}^{q}F& {\mathrm{e}}^{q/2}& 0& 0\\ {\mathrm{e}}^{q/2}& 0& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2{\sin }^2\theta \end{array}\right),f_{\mu \nu }=\left(\begin{array}{cccc}-{\mathrm{e}}^{q}F{\tau }^2& {\mathrm{e}}^{q/2}{\tau }^2& 0& 0\\ {\mathrm{e}}^{q/2}{\tau }^2& \frac{{\Sigma }^2-{\tau }^2}{F}& 0& 0\\ 0& 0& R^2{r}^2& 0\\ 0& 0& 0& R^2{r}^2{\sin }^2\theta \end{array}\right),\text{vspace}-0.3em$$

$$?{\mbox{$S$}}^{\mu}{}_{\nu}?=\begin{pmatrix}\pm_{1}\,|\tau|&~{}~{}\mathrm{e}^{-q/2}\,\dfrac{\pm_{2}\,|\Sigma|-(\pm_{1})\,|\tau|}{F}&0&0\\ 0&\pm_{2}\,|\Sigma|&0&0\\ 0&0&\pm_{3}\,|R|&0\\ 0&0&0&\pm_{4}\,|R|\end{pmatrix}\!.$$

$$f_{\xi \xi }({\xi }_{\mathrm{H}}):=\underset{\xi \to {\xi }_{\text{mbox}\text{H}}}{\lim }\frac{\Sigma (\xi {)}^2-\tau (\xi {)}^2}{F(\xi )}=\mathrm{const}.$$

$$g=\left(\begin{array}{cc}-F& 1\\ 1& 0\end{array}\right),S=\left(\begin{array}{cc}\tau & a\\ b& \Sigma \end{array}\right),$$

$$f=\left(\begin{array}{cc}-F{\tau }^2+2b\tau & {\tau }^2+b\frac{\Sigma -\tau }{F}\\ {\tau }^2+b\frac{\Sigma -\tau }{F}& \frac{{\Sigma }^2-{\tau }^2}{F}\end{array}\right).$$

$$g_{\mu \nu }=\left(\begin{array}{cccc}0& {\mathrm{e}}^{q_{\mathrm{H}}/2}& 0& 0\\ {\mathrm{e}}^{q_{\mathrm{H}}/2}& 0& 0& 0\\ 0& 0& r_{\mathrm{H}}^2& 0\\ 0& 0& 0& r_{\mathrm{H}}^2{\sin }^2\theta \end{array}\right),f_{\mu \nu }=\left(\begin{array}{cccc}0& {\mathrm{e}}^{q_{\mathrm{H}}/2}{\tau }_{\mathrm{H}}^2& 0& 0\\ {\mathrm{e}}^{q_{\mathrm{H}}/2}{\tau }_{\mathrm{H}}^2& \boxed{f_{\xi \xi }({\xi }_{\mathrm{H}})}& 0& 0\\ 0& 0& R_{\mathrm{H}}^2{r}_{\mathrm{H}}^2& 0\\ 0& 0& 0& R_{\mathrm{H}}^2{r}_{\mathrm{H}}^2{\sin }^2\theta \end{array}\right).$$

$$A=\left(\begin{array}{cccc}0& \bullet & 0& 0\\ \bullet & 0& 0& 0\\ 0& 0& \bullet & 0\\ 0& 0& 0& \bullet \end{array}\right),B=\left(\begin{array}{cccc}0& \bullet & 0& 0\\ \bullet & \bullet & 0& 0\\ 0& 0& \bullet & 0\\ 0& 0& 0& \bullet \end{array}\right).$$

$$f_{\xi \xi }({\xi }_{\mathrm{H}})=0.$$

$$g_{\mu \nu }=\left(\begin{array}{cc}{g}_{\mathrm{TR}}& 0\\ 0& g_{\mathrm{ANG}}\end{array}\right),f_{\mu \nu }=\left(\begin{array}{cc}{\tau }_{\mathrm{H}}^2{g}_{\mathrm{TR}}& 0\\ 0& R_{\mathrm{H}}^2{g}_{\mathrm{ANG}}\end{array}\right),$$

$$\mathrm{d}{s}_g^2$$

$$\mathrm{d}{s}_f^2$$

$$\Phi (\xi ):=\frac{\Sigma (\xi )-\tau (\xi )}{F(\xi )}.$$

$$\sigma (\xi ):=\frac{\Sigma (\xi )}{{\left(R(\xi )r(\xi )\right)}^{\prime }}=\frac{\Sigma (\xi )}{R^{\prime }(\xi )r(\xi )+R(\xi )r^{\prime }(\xi )}.$$

$$q^{\prime }r$$

$$F^{\prime }r$$

$${\sigma }^{\prime }r$$

$${\tau }^{\prime }r$$

$$\frac{\Sigma }{2R\sigma F}\left[\frac{m^2{r}^2}{\kappa }{\sigma }^2{\alpha }_1+\tau \left(F+{\sigma }^2\left(\frac{m^2{r}^2}{\kappa }{\alpha }_2-1\right)\right)\right],$$

$${\tau }^{\prime }r$$

$$r^{\prime }$$

Full text

Classical Theories of Gravity, Modified Theories of Gravity, Bimetric Theory, Black Holes

1. I Introduction

2. II The Hassan-Rosen bimetric theory

3. III Bidiagonal ansatz and field equations

4. III.1 Static and spherically symmetric bidiagonal BHs

5. III.2 Choice of the metric fields and behaviour at the Killing horizon

6. III.3 Field equations

7. III.4 Exact initial conditions at the Killing horizon

8. IV Black hole solutions

9. IV.1 Global parameters

10. IV.2 Numerical solutions and phase space analysis

11. IV.2.1 Choice of initial conditions at the Killing horizon

12. IV.2.2 Lyapunov stability and the conjugate system of ODEs

13. IV.2.3 Proper bidiagonality and GR solutions

14. IV.2.4 Lyapunov instability of GR solutions

15. IV.3 The causal structure of the improper bidiagonal solutions

16. V Conclusions

17. A Old parametrisation for the metric fields

18. B Eddington-Finkelstein coordinates

19. C On bidiagonality condition in [44]

20. D Numerical details

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 lrr-2014-4 ; PhysRevLett.116.061102 . Yet, GR cosmology cannot give a satisfactory explanation for the physical origin of the accelerated expansion of the Universe Ade:2015xua . This issue contributes to the so-called “cosmological constant problem” RevModPhys.61.1 , 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 Bergstrom:2000pn ; Bertone:2004pz ; Feng:2010gw . 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 Aoki:2014cla ; Blanchet:2015bia ; Enander:2015kda ; Babichev:2016bxi .

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 Lovelock:1971yv ; Schwartz:2013pla . 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 Fierz:1939 ; Fierz211 . 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 Boulware1972227 ; Boulware:1973my . However, more recently a ghost-free nonlinear theory of a massive spin-2 field — de Rham-Gabadadze-Tolley (dRGT) massive gravity — was proposed deRham:2010ik ; deRham:2010kj . 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 generalisation of the dRGT model to an arbitrary reference metric was given in Hassan:2011tf . That the dRGT massive gravity indeed does not suffer from the BD ghost was proven in Hassan:2011hr .

A further generalisation of the dRGT massive gravity was done in Hassan:2011zd by giving dynamics to the reference metric. For this model, the absence of the BD ghost was proven in Hassan:2011ea . The resulting ghost-free, nonlinear theory of two dynamical interacting spin-2 tensor fields — Hassan-Rosen (HR) bimetric theory — describes the interaction between a massless and a massive spin-2 field, together with their self-interactions and interactions with non-gravitational fields. This theory contains seven propagating modes, two of them associated with the massless spin-2 field and five of them with the massive spin-2 field.

Having constructed these consistent theories of gravity, their viability as alternative theories to GR should be tested by studying their phenomenology. In this respect, the effort has been dedicated to finding black hole (BH) solutions both in dRGT massive gravity and in HR bimetric theory.

The simplest solutions one can consider are static and spherically symmetric. In GR, the no-hair theorems Israel:1967wq ; Israel:1967za ; Carter:1971zc guarantee that a static and spherically symmetric BH without electrical charge is described by only one constant parameter, its mass. The validity of these theorems is not proved in dRGT massive gravity and HR bimetric theory. Then, it is natural to ask whether there are asymptotically flat, static and spherically symmetric solutions, not completely determined by their mass.

Unlike in GR, in HR bimetric theory, on which we will be focusing in the present paper, the static and spherically symmetric field equations cannot, in general, be integrated analytically. Nevertheless, in some regions suitable approximations can be made, allowing for an analytical study of the behaviour of solutions.

In Comelli:2011wq , the authors performed an analytical study of spherically symmetric BH solutions, finding that there are two principal branches of solutions, later denoted bidiagonal and non-bidiagonal.111 Note that the two branches (non-bidiagonal and bidiagonal) were introduced by Salam:1976as in the context of strong-gravity. For the non-bidiagonal solution, there is no coordinate system in which the two metrics can be simultaneously diagonalised, whereas, for the bidiagonal solutions, such coordinate frame exists. The authors show that, in the former case, the equations can be integrated analytically and the solutions are equivalent to those in GR: Schwarzschild, Schwarzschild-de Sitter (SdS) and Schwarzschild-Anti-de Sitter (SAdS). In the latter case, the equations of motion cannot, in general, be solved analytically. However, the authors derived linearised solutions, valid outside the Vainshtein radius Vainshtein:1972sx of the system (inside the Vainshtein radius, by definition, non-linearities must be taken into account). Linearising the equations decouples the massless mode from the massive one, and the metric functions are sums of a Newtonian potential (massless mode) and a Yukawa potential (massive mode).

Bidiagonal, static and spherically symmetric solutions were also studied in Babichev:2013pfa . The authors argue that approximate analytical solutions can be found by solving an algebraic equation, valid inside a region between a point reasonably close to the event horizon (the point reasonably close to the horizon being defined as the point when the gravitational fields become large) and the Compton wavelength of the massive graviton. In this way, the analytical behaviour of the solutions can be also understood inside the Vainshtein radius.

Restricting ourselves to exact, bidiagonal solutions, in Babichev:2014tfa it was shown for the first time that, under suitable conditions, a Kerr $g$-metric and a flat reference metric, $f$, is an exact solution in dRGT massive gravity. Also, the bi-Kerr case with both $g$ and $f$ being Kerr metrics is an exact solution in HR bimetric theory. Since Schwarzschild geometry can be obtained from Kerr, setting the angular momentum to zero, this implies that the Schwarzschild solution with a flat reference metric $f$ is an exact solution in dRGT massive gravity and that the bi-Schwarzschild solution is an exact solution in HR bimetric theory. Charged BH solutions on the Reissner-Nordström-de Sitter form were also found for dRGT massive gravity and HR bimetric theory in Babichev:2014fka . To summarise, all known exact BH solutions of HR bimetric theory correspond to the GR solutions.

In addition to exact and approximate analytical solutions, numerical solutions have been found in the interval $r\in(r_{\mathrm{\scriptscriptstyle H}},+\infty)$, with the Killing horizon radius, $r_{\mathrm{\scriptscriptstyle H}}$, not being included in the domain. A comprehensive analysis of static and spherically symmetric numerical solutions was performed in Volkov:2012wp , in which bidiagonal solutions different from the GR solutions were obtained numerically. Small perturbations around GR solutions were added to the initial conditions at the Killing horizon after which numerical integration was performed, using the new perturbed initial conditions. The results indicated that only SAdS solutions admit perturbations which have AdS-type asymptotic in the leading order. On the other hand, asymptotically flat perturbations of Schwarzschild solutions were not found and SdS solutions were found not to admit non-compact perturbations in the sense that perturbations diverge at finite radii. Numerical studies were also performed in Brito:2013xaa , in which hairy, static and spherically symmetric solutions, being asymptotically flat were explicitly found. The authors proposed that they may represent the final stable state of perturbed Schwarzschild BHs, which are in general unstable Babichev:2013una ; Brito:2013wya ; Babichev:2014oua ; Babichev:2015zub (see below for more details about stability of BH solutions).

Astrophysical and cosmological properties of bidiagonal static and spherically symmetric BHs were studied in Enander:2013kza ; Enander:2015kda . In Enander:2013kza , the bidiagonal solutions found in Comelli:2011wq were rederived using isotropic coordinates to study their lensing properties. The relationship between bidiagonal static and spherically symmetric BHs and stars in HR bimetric theory was discussed in Enander:2015kda . The authors noted that BH solutions and star solutions differ in their asymptotic properties at large radii, indicating that stars cannot gravitationally collapse into bidiagonal BHs.

Stability is of large importance when discussing the physical relevance of BH solutions. Two types of stability are of importance for BH solutions: The first is the stability of a solution against a small change in its initial conditions, the second is the stability of a solution on external perturbations. In the former case, perturbations are made in the phase-space of the system of differential equations, whereas in the latter, perturbations are dynamical (e.g., a radial oscillation for a Schwarzschild BH). In the following, we will refer to the former case as Lyapunov stability and to the latter as dynamical stability. Note that the study of Lyapunov stability is standard practice when considering differential equations in GR, see for example (rendall2008partial, , sections 10.2–10.5).

Linear, dynamical stability analyses of BH solutions were performed in Babichev:2013una ; Brito:2013wya ; Babichev:2014oua ; Babichev:2015zub . In Babichev:2013una it was shown that linear perturbations around bidiagonal bi-Schwarzschild solutions, both in dRGT massive gravity and HR bimetric theory, are unstable, consistent with the results in Brito:2013wya for linear, spherically symmetric perturbations of a Schwarzschild BH. In Brito:2013wya , it was also showed that Kerr BHs are unstable against radial perturbations and superradiant perturbations (a recent review about superradiance can be found in Brito:2015oca ). The dynamical instability of the bidiagonal Schwarzschild solution was confirmed in Babichev:2014oua , where it was also shown that radial perturbations around static and spherically symmetric non-bidiagonal solutions do not exhibit unstable modes. These results were again obtained in Babichev:2015zub , where the authors perform an extensive study of the quasi-normal modes of both bidiagonal and non-bidiagonal solutions. Note that for parameter values corresponding to the partially massless case (for partially masslessness, see Hassan:2012gz and reference therein), the bidiagonal GR solutions (Schwarzschild, SdS and SAdS) are dynamically stable Brito:2013yxa .

General properties of the horizon structure of static and spherically symmetric BHs and their thermodynamics in modified theories of gravity can be found in, e.g., Deffayet:2011rh ; Banados:2011hk . Recent reviews about BH solutions in dRGT massive gravity and HR bimetric theory can be found in Volkov:2013roa ; Volkov:2014ooa ; Babichev:2015xha .

Structure of the paper.

Figure 2 at page 23, apart from describing the new proposed classification for BH solutions, shows the structure of the paper until that point. The caption of the figure contains a detailed description of what we do in each section. In any case, a more general description is reported in the following.

In section II we briefly review the HR bimetric theory formulation. In section III, we describe the theoretical background and clarify the meaning of the bidiagonal ansatz. We also present an original proposition stating that, if one metric is regular and without curvature singularities, to have a non-singular square root matrix $S=\sqrt{g^{-1}f}$ is a necessary condition for not inducing curvature singularities into the other metric. We further introduce a new parametrisation for the two metrics, allowing for a description of the solutions both outside, inside and at the Killing horizon, as well as detecting a possible pathological behaviour of the solutions. This new parametrisation allows us, for the first time, to study the behaviour of the solution a the Killing horizon, and leads us to present two propositions regarding the behaviour of the metrics at their common Killing horizon. Taking into account these results, we propose a new classification for the solutions within the bidiagonal ansatz and discuss the fields equations and initial conditions at the Killing horizon. In section IV, we present our BH solutions and study their convergence properties when $r\rightarrow\infty$, both numerically and analytically through a Lyapunov stability analysis.

Executive summary.

Using a combination of numerical and analytical techniques, we conclude that the only solutions converging to Minkowski, de Sitter (dS) and Anti-de Sitter (AdS) spacetimes for large radii in HR bimetric theory, are the GR solutions: Schwarzschild, Schwarzschild-de Sitter (SdS) and Schwarzschild-Anti-de Sitter (SAdS). For these solutions, the $g$ and $f$ metrics are conformal, i.e., proportional to each other.

Non-GR solutions exist, but they diverge from Minkowski, dS and AdS, even if they differ only infinitesimally from the GR solutions at the Killing horizon. Also, all the found perturbations around Schwarzschild and SdS solutions display singularities in the interaction potential between the two metrics at finite radii. For this reason, we do not consider them as acceptable solutions.

These results, concerning the asymptotic behaviour of the solutions, are compatible with those in Volkov:2012wp , where the leading order asymptotic properties of the solutions are studied for $r\rightarrow\infty$. Therein, the relative differences between non-GR solutions and GR ones are taken into account, whereas, in this paper, we focus on the absolute differences between them (which is the relevant property for Lyapunov stability). We note that the behaviour of the relative differences does not tell anything about the behaviour of the absolute difference (i.e., about the Lyapunov stability); therefore, the results of Volkov:2012wp do not imply ours. Therefore, the main motivation of this work is to extend the results of Volkov:2012wp by studying the Lyapunov stability of the solutions.

Moreover, we study for the first time the causal structure of BH spacetimes in HR bimetric theory at the common Killing horizon by using Eddington-Finkelstein coordinates. In the light of this analysis, we propose a new classification for the BH solutions and point out that choosing a bidiagonal ansatz outside or inside the Killing horizon, does not guarantee bidiagonality exactly at the Killing horizon.

As already pointed out in Volkov:2012wp , 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 $\Lambda=0$, then we will have two Schwarzschild solutions. Similarly for SAdS when we have different $\Lambda<0$ and for SdS when we have different $\Lambda>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:

1. (i)

the only bidiagonal BH solutions in HR bimetric theory asymptotically converging to Minkowski, dS and AdS spacetimes are GR solutions with conformal metrics, 2. (ii)

bidiagonal BH solutions display dynamical linear instabilities, 3. (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 Hassan:2011zd ,

[TABLE]

where $M_{g},M_{f}$ are the reduced Planck masses and $R^{g},R^{f}$ are the Ricci scalars respectively for $g_{\mu\nu}$ and $f_{\mu\nu}$. The kinetic terms for $g_{\mu\nu}$ and $f_{\mu\nu}$ have the standard Einstein-Hilbert form. The interaction potential is given in terms of the square root matrix $S\coloneqq\sqrt{g^{-1}f}$ (so that ?{\mbox{S}}^{\mu}{}_{\rho}??{\mbox{S}}^{\rho}{}_{\nu}?=g^{\mu\rho}f_{\rho\nu}) through the elementary symmetric polynomials $e_{n}(S)$. These polynomials are the scalar invariants of $S$ and can be expressed in the following way,

[TABLE]

which can be expanded (noting that $e_{n}(S)\equiv 0$ for $n>4$),

[TABLE]

The $\beta_{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_{\mu\nu}$ and $f_{\mu\nu}$ separately, yielding,

[TABLE]

where $R^{g}_{\mu\nu}$ and $R^{f}_{\mu\nu}$ are the Ricci tensors of $g_{\mu\nu}$ and $f_{\mu\nu}$, respectively, and,

[TABLE]

Finally, the Bianchi constraints read,

[TABLE]

III Bidiagonal ansatz and field equations

This section is organised as follows. In subsection III.1, 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 subsection III.2, 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 subsection III.3, we expand the field equations as nonlinear ODEs and discuss their properties. Finally, in subsection III.4, we consider the exact initial conditions at the Killing horizon (not determined in previous works).

@fb@secFB

III.1 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 Comelli:2011wq . Whether this is the case also for bidiagonal solutions is the main question we want to address in this paper.

As proved in Deffayet:2011rh ; Banados:2011hk , 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 $\mathcal{K}$. Since the geometry of $g$ is static and spherically symmetric, assuming that the vector field $\mathcal{K}$ is unique, we can adopt coordinates as prescribed in wald1984 : We begin by selecting a spacelike surface $\mathcal{S}$ orthogonal to the orbits of the one-parameter isometry generated by $\mathcal{K}$. We use $t$ to denote this parameter and refer to it as the ‘time coordinate’. The image of this isometry are the surfaces $\mathcal{S}_{t}$ also orthogonal to $\mathcal{K}$, which we label by $t$ so that all the points at $\mathcal{S}_{t}$ have the same $t$. Now, we select a two-sphere on the spacelike surface $\mathcal{S}_{t}$ and choose angular coordinates $(\theta,\phi)$ where the metric on each 2-sphere takes the form $r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2}\right)$. 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 $\mathcal{S}_{t}$ along the radial geodesics orthogonal to the selected sphere. Assuming a non-vanishing $\nabla_{\mu}r\neq 0$, we select $(r,\theta,\phi)$ as the coordinates on $\mathcal{S}_{t}$, resulting in the spacetime coordinates $(t,r,\theta,\phi)$. Allowing further for the reparametrisation of $r$ in terms of a new coordinate $\xi$ (with the unitary gauge $r=\xi$), we arrive at the coordinate chart $x^{\mu}=(t,\xi,\theta,\phi)$. The introduction of $\xi$ 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 $\xi$ as the independent variable. In these coordinates, the metric $g$ takes the form,

[TABLE]

with $g_{tt}(\xi)=-\mathcal{K}^{2}$. Note that $r(\xi)$ is necessarily monotonic choquet2008general . Observe also that this chart breaks down at points where $\mathcal{K}$ or $\nabla_{\mu}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\left(U,S(S(V)\right)$, we have $g\left(S(U),V\right)=g\left(U,S(V)\right)$ and $f(U,V)=g\left(S(U),S(V)\right)$, 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 higham2008 .

The equivalent expressions of the above statements in matrix notation are: $f=gS^{2}$ ($S$ is the square root), $\left(gS\right)^{\mathsf{{\scriptscriptstyle T}}}=S^{\mathsf{{\scriptscriptstyle T}}}g=gS$ ($S$ is self-adjoint), and $f=S^{\mathsf{{\scriptscriptstyle 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,

[TABLE]

where $\tau(\xi)$, $\Sigma(\xi)$ and $R(\xi)$ are real fields. The set of independent $\pm_{n}$ in (8) indicates the square root branch (giving 16 branches in total for $n=1,...,4$). The branch will be selected later in subsection III.3.

Consequently, we have $f_{\xi\xi}(\xi)=\Sigma(\xi)^{2}\,g_{\xi\xi}$ and $f_{tt}(\xi)=\tau(\xi)^{2}\,g_{tt}(\xi)=-\mathcal{K}^{\prime 2}$, where $\mathcal{K}^{\prime 2}$ is the norm of the translational Killing vector $\mathcal{K}$ taken with respect to $f$; hence,

[TABLE]

The radial coordinate for $f$ is effectively $R(\xi)r(\xi)$, which must be monotonic. The reparametrisation $U(\xi)=R(\xi)r(\xi)$ will cast (9) into a standard form, as in (7). For a non-vanishing $\tau(\xi)$, $\mathcal{K}^{\prime 2}=0$ if and only if $\mathcal{K}^{2}=0$. Hence, the ansatz (7)–(8) ensures a common Killing horizon.

In the above, we required a non-vanishing $\tau(\xi)$ 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(\cdot,\cdot)=g(S(\cdot),S(\cdot))$. 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 kupeli1987 ; kupeli1996 . 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 Volkov:2012wp , 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 straumann2013general . 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, $\tau(\xi)\Sigma(\xi)$ could be finite with $\tau(\xi)\to\infty$ and $\Sigma(\xi)\to 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, $\tau(\xi)\neq 0$, $\Sigma(\xi)\neq 0$ and $R(\xi)\neq 0$ for all $\xi$.

@fb@secFB

III.2 Choice of the metric fields and behaviour at the Killing horizon

The fields $g_{tt}(\xi)$ and $g_{\xi\xi}(\xi)$ in (7) can be denoted as $q(\xi)$ and $F(\xi)$, so that $F(\xi)=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,

[TABLE]

[TABLE]

with the induced,

[TABLE]

From now on, we employ the following naming conventions:

• •

The fields at an arbitrary point $\xi$ are denoted just by their name.

Example: $r\coloneqq r(\xi)$, $q\coloneqq q(\xi)$, $R\coloneqq R(\xi)$ etc.

• •

The derivatives of the fields (denoted with primes) are always with respect to $\xi$.

• •

The fields at some designated point, $\xi=\xi_{\mathrm{label}}$, inherit the same ‘label’.

Example: At $\xi=\xi_{\mathrm{\scriptscriptstyle H}}$, we have $r_{\mathrm{\scriptscriptstyle H}}\coloneqq r(\xi_{\mathrm{\scriptscriptstyle H}})$, $q_{\mathrm{\scriptscriptstyle H}}\coloneqq q(\xi_{\mathrm{\scriptscriptstyle H}})$, $R_{\mathrm{\scriptscriptstyle H}}\coloneqq R(\xi_{\mathrm{\scriptscriptstyle H}})$ etc.

As earlier pointed out, the coordinate system $(t,\xi,\theta,\phi)$ is not suitable to go through the Killing horizon where $\mathcal{K}^{2}$ vanishes. For this purpose we introduce the ingoing Eddington-Finkelstein coordinates $x^{\mu}=(v,\xi,\theta,\phi)$ adapted for $g_{\mu\nu}$. In these coordinates, the metrics $g$ and $f$ are (see Appendix B for more details),

[TABLE]

with the square root,

[TABLE]

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_{\xi\xi}$ 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(\xi_{\mathrm{\scriptscriptstyle H}})=0$ is that the following limit is finite,

[TABLE]

Proof.

The proof simply follows from (13).222The form of $S$ does not influence the result. Consider the $2\times 2$ blocks of $g$ and $S$ in the form,

$g=\begin{pmatrix}-F&1\\ 1&0\end{pmatrix}\!,\qquad S=\begin{pmatrix}\tau&a\\ b&\Sigma\end{pmatrix}\!,$

so that $f=gS^{2}$. The symmetrisation $f=f^{\mathsf{{\scriptscriptstyle T}}}$ implies $a=(\Sigma-\tau)/F$ without any restriction on $b$. Then (with $b=0$ in our case),

$f=\begin{pmatrix}-F\tau^{2}+2b\tau&~{}\tau^{2}+b\frac{\Sigma-\tau}{F}\\ \tau^{2}+b\frac{\Sigma-\tau}{F}&~{}\frac{\Sigma^{2}-\tau^{2}}{F}\end{pmatrix}\!.$

∎

We refer to the condition (15) as the crossing condition.333This crossing condition is similar to the GR one for $g=F/(UV)\,\mathrm{d}U\mathrm{d}V$ in Kruskal-Szekeres coordinates, where $F/(UV)$ should remain finite at $U=0$ or $V=0$, which correspond to the Killing horizons. 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 subsection IV.3.

From now on, we assume that the crossing condition is satisfied. In the limit $\xi\rightarrow\xi_{\mathrm{\scriptscriptstyle H}}$ when $F(\xi)\rightarrow 0$, we have,

[TABLE]

Now, denoting non-zero components by $\bullet$, two real symmetric matrices $A$ and $B$ cannot be simultaneously diagonalised if they have the following structure,

[TABLE]

This follows from the theorem on canonical pair forms (see Uhlig:1973 ; Hassan:2017ugh 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(\xi_{\mathrm{\scriptscriptstyle H}})=0$, if and only if,

[TABLE]

We refer to the condition (18) as the proper bidiagonality condition. A similar statement can be made by extending Proposition 1 in Deffayet:2011rh (see Appendix 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_{\mu\nu}$ are the Eddington-Finkelstein coordinates for $f_{\mu\nu}$ too, and they are a null frame for both metrics. This is obvious from (16) after setting $f_{\xi\xi}=0$. In such case, having also $g_{vv}=f_{vv}=0$, the $(v,\xi)$ coordinates become null for both metrics.

Proposition 3 contains two pieces of information: The first is that, if (18) and consequently $\Sigma_{\mathrm{\scriptscriptstyle H}}^{2}=\tau_{\mathrm{\scriptscriptstyle H}}^{2}$ 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 $\tau_{\mathrm{\scriptscriptstyle H}}^{2}$ and for the angular part, $R_{\mathrm{\scriptscriptstyle H}}^{2}$. We refer to the time-radial 2-block of the metric with $g_{\mathrm{\scriptscriptstyle TR}}$ and to the angular 2-block with $g_{\mathrm{\scriptscriptstyle ANG}}$. At the Killing horizon,

[TABLE]

hence, the metrics are block-conformal with conformal factors $\tau_{\mathrm{\scriptscriptstyle H}}^{2}$ and $R_{\mathrm{\scriptscriptstyle H}}^{2}$. 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.444The Jacobian (72) is ill-defined at the Killing horizon independently of the proper bidiagonality condition.

An example of coordinates that diagonalises the metrics, if the proper bidiagonality condition is true, is, $T\coloneqq(v+r)/\sqrt{2}$, $X\coloneqq(-v+r)/\sqrt{2}$ (which is the $\pi/4$ clockwise rotation of the two basis vectors $\partial_{v}$ and $\partial_{\xi}$). In these coordinates, we obtain the following two block-conformal metrics at horizon,

[TABLE]

If we use the same coordinates without imposing the proper bidiagonality condition, $f$ will be manifestly non-diagonal with a non-zero $\mathrm{d}T\mathrm{d}X$ component.

This different behaviour at the horizon suggests a new classification for the solutions:

1. (i)

proper bidiagonal solutions, for which the proper bidiagonality condition holds, and 2. (ii)

improper bidiagonal solutions, for which only the crossing condition holds. The improper bidiagonal solutions are bidiagonal everywhere except at the Killing horizon.

@fb@secFB

III.3 Field equations

To expand the field equations (4), we need to specify the branch of the square root. We select the principal branch where $\pm_{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,

[TABLE]

According to eq. (15), $\Phi(\xi)$ must be finite at the horizon. For increased readability, we also introduce the auxiliary field $\sigma(\xi)$ by,

[TABLE]

Note that $\sigma$ is not to be considered an independent field since its dynamics is governed by $\Sigma$. Also, given that $rR$ must be monotonic choquet2008general and $\Sigma$ must be non-zero and non-diverging due to the regularity of the elementary symmetric polynomials of $S$, $\sigma$ must also be non-zero and finite.

With the fields given in (10)–(12), gauge-fixing $r(\xi)=\xi$, and restricting ourselves to the principal square root branch, the field equations are,

[TABLE]

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 $\alpha_{n}[R]$ defined as,

[TABLE]

which come out from the interaction between the two metrics; they contain the $\beta$-parameters and are always multiplied by the energy scale $m$ of the interaction potential.555Note that our $\alpha$-utility functions are the analogous of the same-named functions defined in Volkov:2012wp , but shifted for index 1 and written with our choice of the metric fields. Since all the fields in (23) are only functions of $\xi$ and there is no explicit dependence on $\xi$ in the equations, the system is an autonomous system of nonlinear ODEs with $\xi$ 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.666This scaling property can be used to make the radial coordinate dimensionless by setting the Compton wavelength to 1. By simply replacing $r\to 1$ and $m\to m\,\mathrm{e}^{\tilde{r}}$, the system will be given in terms of $\tilde{r}=\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^{\prime}=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 $\tau$ and $F$ have a linear dependence. From eqs. (23d) and (23e), we get,

[TABLE]

Taking the derivative of either $\tau$ or $F$ from (25), and plugging it back into (23), enables us to solve for $R^{\prime}$ (or $\Sigma$).777Following a procedure similar to that described in Volkov:2012wp . Although $\tau$ and $F$ could be used on equal footing, the remnant equations are simpler for $\tau$. Thus, we select $F$ to be dynamical and $\tau$ 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^{\prime}$, we also know $\Sigma$ and $\Phi$ algebraically (again noting that $\sigma$ took over dynamics from $\Sigma$). 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, $\mathrm{e}^{q}\mathrm{d}t\to\mathrm{d}t$. As a consequence, we can integrate eq. (23a) and determine $q$ after the other fields are known (up to an integration constant).

To summarise: $\tau$, $\Sigma$ and $q$ can be determined once the fields $r$, $F$, $\sigma$ and $R$ are known, and for the latter we have an autonomous system of four nonlinear ODEs, where three of them are coupled,

[TABLE]

Here, we introduced two vector fields, $\mathscr{V}$ and $\mathscr{F}$, to simplify expressions (these vector fields will be extensively used later on in the paper). The fields $\mathscr{F}_{F}$ and $\mathscr{F}_{\sigma}$ are given by (23b) and (23c), respectively. The expression for $R^{\prime}=\mathscr{F}_{R}$ will not be quoted here as its total number of indivisible sub-expressions is around 2600.

@fb@secFB

III.4 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).888Since we are dealing with rational functions, the existence and the uniqueness of the solutions to (26) are guaranteed by the continuous differentiability of $\mathscr{F}$ in all its variables and parameters, unless at poles. 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\in(0,+\infty)$.

Denoting one of the Killing horizons $\xi_{\mathrm{\scriptscriptstyle H}}$, we can start integrating from some $r_{\mathrm{\scriptscriptstyle H}}=r(\xi_{\mathrm{\scriptscriptstyle H}})$.999Here we assume the existence of at least one Killing horizon. Observe that a solution to (23) can have none, one or many Killing horizons satisfying $F(\xi)$=0. In particular, some solutions can have no Killing horizons, if we start integration from $F(\xi)\neq 0$. This is our first of four initial conditions. The second exact initial condition at the horizon is, of course, $F(\xi_{\mathrm{\scriptscriptstyle 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,

[TABLE]

at the horizon. The algebraic expression for $\tau$ is known and, as pointed out earlier, it is possible to derive an equation for the combination $R+rR^{\prime}$, which is inside $\Sigma$. It turns out that both $\tau$ and $\Sigma$ depend on $F$ as rational functions,

[TABLE]

where $P$ indicates the polynomial in $F$ at the numerator and $Q$ the polynomial in $F$ at the denominator. Therefore, we can write,

[TABLE]

The coefficients $p_{i}$ and $q_{j}$ are functions of $R$, $\sigma$ and $r$. In order to have a finite limit when $F\rightarrow 0$, the term $p_{0}$ must vanish at the horizon,

[TABLE]

After fixing $r_{\mathrm{\scriptscriptstyle H}}$, eq. (30) is fourth-order in $R_{\mathrm{\scriptscriptstyle H}}$ and third-order in $\sigma_{\mathrm{\scriptscriptstyle H}}$ and we can solve it for either. We consider $R_{\mathrm{\scriptscriptstyle H}}$ as a free parameter, and solve eq. (30) for $\sigma_{\mathrm{\scriptscriptstyle H}}$ given a value of $R_{\mathrm{\scriptscriptstyle H}}$. It turns out that there are always only two independent non-zero solutions for $\sigma_{\mathrm{\scriptscriptstyle 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,

[TABLE]

Therefore, after having solved for $\sigma_{\mathrm{\scriptscriptstyle 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 $\Phi_{\mathrm{\scriptscriptstyle H}}=\left.(p_{1}/q_{0})\right|_{\xi=\xi_{\mathrm{\scriptscriptstyle H}}}$, then $p_{1}(r_{\mathrm{\scriptscriptstyle H}},R_{\mathrm{\scriptscriptstyle H}},\sigma_{\mathrm{\scriptscriptstyle H}})=0$ if and only if we have a proper bidiagonal solution. Solving the last equation is equivalent to finding proper bidiagonal solutions.

To summarise, besides $r_{\mathrm{\scriptscriptstyle H}}$, we only have one free parameter, $R_{\mathrm{\scriptscriptstyle H}}$. Once specified, and after choosing the branch for $\sigma_{\mathrm{\scriptscriptstyle 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.1 is about the choice of global parameters, i.e., parameters appearing in the action (1). In subsection IV.2, we describe the detailed analysis of found BH solutions. In particular, in subsection IV.2.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 subsection IV.3, we study the causal structure of improper bidiagonal solutions both inside and outside the Killing horizon.

@fb@secFB

IV.1 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 Fierz:1939 ; Fierz211 is $m_{\mathrm{g}}=10^{-1}$. To facilitate comparison with the results of Volkov:2012wp , we use the default set of global parameters given in Table 1. 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^{\mu}=(t,r,\theta,\phi)$, we demand that the bi-Minkowski solution,

[TABLE]

is a solution of eq. (26), and obtain the so-called asymptotic flatness conditions,101010These conditions are actually the ‘flatness’ conditions, because the bi-Minkowski solution is valid everywhere and not only asymptotically.

[TABLE]

These conditions can also be deduced from cosmological solutions vonStrauss:2011mq . In addition, we use the normalisation condition,

[TABLE]

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}\rightarrow m^{2}\delta$. By using the normalisation condition, we effectively trade one of the free parameters $\beta_{1},\beta_{2},\beta_{3}$ in favour of $m_{\mathrm{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 $\beta$ parameters by rescaling $\beta_{n}\to\beta_{n}c^{n}$.

@fb@secFB

IV.2 Numerical solutions and phase space analysis

IV.2.1 Choice of initial conditions at the Killing horizon

By definition, in the 4-dimensional phase space $\left(r,F,\sigma,R\right)$ 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 $\mathcal{A}$ embedded in the hyperplane $F=0$. To be more precise, eq. (30) defines an algebraic variety hartshorne1977alg , i.e., the set of solutions of a system of polynomial equations, in our case eq. (30). The surface $\mathcal{A}$ is shown in Figure1a. The red and blue regions represent the two independent branches for $\sigma_{\mathrm{\scriptscriptstyle H}}$ (which are the solutions of eq. (30)), as functions of $R_{\mathrm{\scriptscriptstyle H}}$ and $r_{\mathrm{\scriptscriptstyle H}}$. Setting the Killing horizon to some $r_{\mathrm{\scriptscriptstyle H}}=r_{*}$, defines a 1-dimensional curve $\mathcal{C}\subset\mathcal{A}$, as the intersection of the plane $r_{\mathrm{\scriptscriptstyle H}}=r_{*}$ and $\mathcal{A}$ (see Figures1a and1b, where $r_{*}=1$).111111 Note that the distinction between the branches is artificial. We can always introduce a change of coordinates $(r,\sigma,R)\to(\tilde{r},\tilde{\sigma},\tilde{R})$ so that any point becomes centred in the new coordinates with the tangent at the algebraic variety being horizontal.

In Figure1b, the vertical solid lines indicate the values of $R_{\mathrm{\scriptscriptstyle H}}$ for which $\Phi_{\mathrm{\scriptscriptstyle H}}=0$, i.e., the corresponding solutions are proper bidiagonal. The vertical dashed lines represent improper bidiagonal solutions close to the proper bidiagonal ones. Once we have $R_{\mathrm{\scriptscriptstyle H}}$, we must choose one branch for $\sigma_{\mathrm{\scriptscriptstyle H}}$ (red or blue in the top panel). In this example, we choose the blue branch for all the five solutions considered. $\circ$ indicates the Schwarzschild solution (being a proper bidiagonal solution), $\square$ indicates an improper bidiagonal solution close to the Schwarzschild one, $\bullet$ indicates the SAdS solution (also a proper bidiagonal solution) and $\blacksquare$ an improper bidiagonal solution close to SAdS. The SdS solution is not considered in Figure 1 since it has $R_{\mathrm{\scriptscriptstyle H}}<0$ for this parameter set.121212The meaning of $R_{\mathrm{\scriptscriptstyle H}}<0$ is clarified in subsection IV.2.4. In subsubsection IV.2.2, we will explain in more detail why GR solutions correspond to these initial conditions. $\diamond$ indicates a proper bidiagonal solution which is non-GR (we come back to this solution in subsubsection IV.2.3).

Summarising, for a BH to have a Killing horizon at some $r_{*}$, the metric fields must assume values lying on $\mathcal{C}$ at $r=r_{*}$. Starting from initial conditions belonging to $\mathcal{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_{\mathrm{\scriptscriptstyle H}}$ belonging to the domain of $\mathcal{C}$ are allowed (that is, for which eq. (30) has real solutions).

We choose the Killing horizon to be at $r_{\mathrm{\scriptscriptstyle H}}=1$. Given that $m_{\mathrm{g}}=10^{-1}$, the Compton wavelength of the massive mode is $\lambda=10$, and $r_{\mathrm{\scriptscriptstyle H}}/\lambda=10^{-1}$. The Fierz-Pauli mass is observationally constrained to be very small, $m_{\mathrm{g}}\lessapprox 10^{-30}-10^{-33}$eV deRham:2014zqa . 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_{\mathrm{\scriptscriptstyle H}}$ becomes smaller (see Figure1a), 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.

IV.2.2 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 Lyapunov1992General :

Definition** (Lyapunov stability).**

Consider a generic autonomous system of ODEs with independent variable $t$ and let $\mathscr{S}$ denote the set of all its solutions. A solution $x(t)\in\mathscr{S}$ is called stable if, for every $\epsilon>0$, there exist a $\delta_{\epsilon}>0$ and a $t_{0}$ such that the following holds,

[TABLE]

If, in addition, there exist a $\delta>0$ and a $\bar{t}$ such that the following holds,

[TABLE]

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 (rendall2008partial, , 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 arrowsmith1990 . 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 $\mathcal{V}^{\prime}(\xi)\equiv 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 $C^{1}$-diffeomorphism, $x(t)\rightarrow x(t)-x_{0}(t),$ where $x(t)$ is another generic solution of the system of ODEs Lyapunov1992General . 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^{\infty}$-diffeomorphism given by,

[TABLE]

defined in $r\in(0,+\infty)$ and where $\Lambda\in\mathbb{R}$ is constant. Here, $\{\rho,\bar{F},\bar{\sigma},\bar{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 $\rho=0$.131313Note that $\rho=0$ is outside of the original domain, but it can be included by compactification.

Plugging the shifted functions into eqs. (23b), (23c), (23f) and into the equation for $R^{\prime}$ in eq. (26), we obtain the conjugate ODEs system,

[TABLE]

where all the fields are functions of $\xi$. Here, $\Sigma$ is a function of $\bar{\mathscr{F}}_{\bar{R}}(\rho,\bar{F},\bar{\sigma},\bar{R})$ rather than $\bar{R}^{\prime}$. Also, $\alpha_{n}$ and $b_{n}$ now stand for $\alpha_{n}[\bar{R}+R_{0}]$ and $b_{n}[\bar{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 $\left(\bar{F},\bar{\sigma},\bar{R},\rho\right)\equiv\left(\bar{F}_{0},\bar{\sigma}_{0},\bar{R}_{0},\rho_{0}\right)$, with the subscript 0 indicating a constant value. If $\bar{F}_{0}\neq 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 arrowsmith1990 , we conclude that we can have a fixed point of the conjugate system only if $\bar{F}_{0}=0$. Besides, $\rho^{\prime}\equiv 0$ if and only if $\rho\equiv 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 $\left(\bar{F},\bar{\sigma},\bar{R},\rho\right)\equiv\left(0,\bar{\sigma}_{0},\bar{R}_{0},0\right)$ 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 $\rho$ in the denominators),

[TABLE]

where $R_{0}+\bar{R}_{0}=R_{\mathrm{\scriptscriptstyle H}}$ from (37). In particular, condition (a) makes $\bar{F}^{\prime}$ vanish identically, whereas conditions (b) and (c) combined make $\bar{\sigma}^{\prime}$ vanish identically, and condition (d) is equivalent to $\bar{R}^{\prime}=0$ identically. By imposing these necessary conditions, we are considering all the fixed points of the conjugate system of the form,

[TABLE]

Keeping $\rho$ arbitrary and substituting $\left(\rho,\bar{F},\bar{\sigma},\bar{R}\right)=\left(\rho,0,0,\bar{R}_{0}\right)$ in eq. (37), together with the conditions in eq. (39), we obtain (coming back to the original system),

[TABLE]

Finally, using eqs. (25a) and (41), we obtain the following quartic equation in $R$ which yields all possible $R_{\mathrm{\scriptscriptstyle H}}$,

[TABLE]

with the resulting (see eq. (39a)),

[TABLE]

If we require that $R_{\mathrm{\scriptscriptstyle H}}=1$ is a flat solution, i.e., that $0=\Lambda(1)=\alpha_{0}[1]+\alpha_{1}[1]=\alpha_{1}[1]+\alpha_{2}[1]$, we get the asymptotic flatness conditions (33). Conversely, imposing (33) ensures that $R_{\mathrm{\scriptscriptstyle 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\rightarrow\infty$. Second, this shows that GR solutions are exact solutions of eq. (26). Note that the cosmological constant of the GR solutions depend on $R_{\mathrm{\scriptscriptstyle H}}$, the free parameter in our initial conditions. Therefore, depending on $R_{\mathrm{\scriptscriptstyle H}}$, we will have Schwarzschild, SdS or SAdS solutions (see Figure1b). The asymptotic structure is thus determined by the initial conditions at the Killing horizon. For example, for $R_{\mathrm{\scriptscriptstyle H}}=1$, $\Phi\equiv 0$ and $\Lambda(R_{\mathrm{\scriptscriptstyle H}})=0$, i.e., the Schwarzschild solution. Note again that the asymptotic flatness conditions (33) guarantee that $R_{\mathrm{\scriptscriptstyle H}}=1$ is always a possible initial value for $R(\xi)$. Substituting the fields in eq. (41) into the metrics, the metrics become conformal, with conformal factor $R_{\mathrm{\scriptscriptstyle H}}^{2}$. These results are all compatible with previous results (see, for example, Volkov:2012wp ). After the considerations made in this subsection, we can better understand Figure 1, in which GR solutions have $\sigma_{\mathrm{\scriptscriptstyle H}}=1$ and $\Phi_{\mathrm{\scriptscriptstyle H}}=0$.

We note that, $\Lambda$ being determined by $R_{\mathrm{\scriptscriptstyle H}}$ through eq. (43), all values of $R_{\mathrm{\scriptscriptstyle H}}$ for which $\Lambda(R_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle H}}$ has real solutions for $\sigma_{\mathrm{\scriptscriptstyle H}}$ and $\Lambda(R_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle H}}$. In Volkov:2012wp 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.

IV.2.3 Proper bidiagonality and GR solutions

Equation (41) shows that GR solutions have $\Phi=0$ identically.

[TABLE]

The constraint $\Phi_{\mathrm{\scriptscriptstyle H}}=0$ defines a finite number of points on $\mathcal{C}$, hence a finite number of proper bidiagonal solutions (we saw three of them in Figure1b). Proper bidiagonality does not imply GR solutions, as we can see in Figure1b. The solution indicated by $\diamond$ has $\Phi_{\mathrm{\scriptscriptstyle H}}=0$, $\sigma_{\mathrm{\scriptscriptstyle H}}=1,R_{\mathrm{\scriptscriptstyle H}}=3/2$ but $R^{\prime}_{\mathrm{\scriptscriptstyle H}}\neq 0$ at the Killing horizon; therefore, it is not a fixed point (i.e., a GR solution).

However, we can prove the following,

[TABLE]

Assuming $\Phi\equiv 0$ implies $\Sigma\equiv\tau$ and

[TABLE]

Substituting $\sigma^{\prime}$ and $\tau^{\prime}$ into $\Phi^{\prime}$ by using the equations of motion, and substituting $\tau\equiv\Sigma=\sigma\left(rR^{\prime}+R\right)$, we obtain,

[TABLE]

where all the metric fields are functions of $\xi$. By hypothesis, $\Phi^{\prime}$ is identically zero. The second term depends on $\beta$ parameters which are inside the utility functions $\alpha_{1}$, $b_{1}$ and $b_{2}$, whereas the first term does not depend on them. Therefore, in order for $\Phi^{\prime}$ to be identically zero, independently of the $\beta$ parameters, $R$ must be constant (i.e., $R\equiv R_{\mathrm{\scriptscriptstyle H}}$) and $\sigma\equiv 1$. This implies $\Sigma\equiv\tau\equiv R\equiv R_{\mathrm{\scriptscriptstyle H}}$, i.e., $g_{\mu\nu}$ and $f_{\mu\nu}$ 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.

IV.2.4 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 Appendix D for more details about the accuracy of these solutions). Using the parameter values given in Table 1, the set of the fixed points (42) gets the values,

[TABLE]

This is consistent with Volkov:2012wp . When used as the initial values, these correspond to GR solutions. First, note that the value $R_{\mathrm{\scriptscriptstyle H}}=1$ is exact, guaranteed by the asymptotic flatness conditions in eq. (33). Second, there are two negative values for $R_{\mathrm{\scriptscriptstyle H}}$. For GR solutions, $R\equiv R_{\mathrm{\scriptscriptstyle H}}$; For non-GR solutions, instead, given that $R$ cannot change sign (since this implies $\det S=0$), if $R_{\mathrm{\scriptscriptstyle H}}$ is negative, then $R$ will be negative for all $\xi$ (i.e., all $r$). The same holds for $\sigma$ which is positive (because $\sigma_{\mathrm{\scriptscriptstyle H}}=1$). Then, from eq. (14), we can deduce that different signs of $R_{\mathrm{\scriptscriptstyle 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 $\sigma_{\mathrm{\scriptscriptstyle H}}$. A negative $R$ is equivalent to retaining the absolute value of $R$ in eq. (14) together with $\pm_{3}=\pm_{4}=-1$. We emphasize that the equations of motions (23) are valid only for $R>0$. However, the transformation,

[TABLE]

is equivalent to changing the square root branch, but for different $\beta$ parameters and different asymptotic flatness conditions (i.e., in a different phase space).

Choosing $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{Schw}}\equiv 1$, we obtain the Schwarzschild solution. Perturbations around it can be obtained by considering $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{Schw}}(1+\epsilon)$ and selecting the same branch for $\sigma_{\mathrm{\scriptscriptstyle H}}$ as the Schwarzschild solution (Figure1b). Hence, these perturbations are improper bidiagonal (like the perturbations around SAdS and SdS, treated later). All found perturbations around Schwarzschild solutions have a diverging $\sigma(r)$ at some finite radial coordinate $r=\hat{r}$; the bigger $\epsilon$, the smaller $\hat{r}$. A diverging $\sigma(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 $\sigma(r)\rightarrow s(r)\coloneqq\sigma(r)^{-1}$ and solve for $s(r)$, which is 0 when $\sigma(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,

[TABLE]

with $R_{\mathrm{\scriptscriptstyle H}}=1$ and the numerical perturbation around it, $F_{\mathrm{pert}}$, with $R_{\mathrm{\scriptscriptstyle 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 Brito:2013xaa where asymptotically flat solutions other than the Schwarzschild were found, our result shows that Schwarzschild solution is the unique asymptotically flat solution (see also Appendix 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\rightarrow\infty$, independent of the global parameter values.

The solutions with $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{SAdS}}\approx 2.633$ and $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H},2}^{\mathrm{SAdS}}\approx-0.6459$ correspond to SAdS solutions. Perturbations around these two solutions behave very similarly, so we describe only $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{SAdS}}$, belonging to the same phase space as the Schwarzschild solution. Plotting the metric functions, the SAdS solution naïvely appears asymptotically stable. However, the study of the difference $F_{\mathrm{pert}}-F_{\mathrm{SAdS}}$ at large radii, which matters for Lyapunov stability, shows that it is unstable. Here, $F_{\mathrm{SAdS}}$ indicates the analytic SAdS solution,

[TABLE]

which is 0 when $r=r_{\mathrm{\scriptscriptstyle H}}$. The difference between the improper bidiagonal perturbation around SAdS with $R_{\mathrm{\scriptscriptstyle H}}=3$ and the analytic SAdS solution with $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{SAdS}}$ is plotted in Figure 4. $F_{\mathrm{pert}}-F_{\mathrm{SAdS}}$ oscillates with increasing amplitude at radii larger than the Compton wavelength of the massive mode ($\lambda=10$ in our case). A similar asymptotic behaviour around SAdS was found in Sushkov:2015fma .

Such a result shows that SAdS is an unstable solution. We emphasise that all the found perturbations around SAdS have a non-diverging $\sigma(r)$ and a monotonically increasing $R(r)r$ (which is a necessary condition to have a spherically symmetric spacetime choquet2008general ). 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_{\mathrm{\scriptscriptstyle H}}=R_{0}^{\mathrm{SdS}}\approx-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 $\sigma(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 Volkov:2012wp , 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\rightarrow\infty$. On the other hand, integrating inwards from the cosmological horizon, $\sigma(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 around SdS solution, we integrate numerically from the cosmological horizon, imposing initial conditions lying on the algebraic variety. Note that $R_{\mathrm{\scriptscriptstyle C}}=R_{\mathrm{\scriptscriptstyle H}}$ for the SdS solution (because $R$ is constant), of course corresponding to different $r_{\mathrm{\scriptscriptstyle H}}$ and $r_{\mathrm{\scriptscriptstyle C}}$.

The SdS solution with Killing horizon $r_{\mathrm{\scriptscriptstyle H}}=1$ has the cosmological horizon at $r=r_{\mathrm{\scriptscriptstyle C}}\approx 3.128$. In Figure 5, we show the difference between the analytical SdS solution in eq. (51) with $R_{\mathrm{\scriptscriptstyle H}}=R_{0}^{\mathrm{SdS}}\approx-8.557$ and $\Lambda(R^{\mathrm{SdS}}_{\mathrm{\scriptscriptstyle H}})>0$, and a perturbation around it, with $R_{\mathrm{\scriptscriptstyle H}}=R_{\mathrm{\scriptscriptstyle H}}^{\mathrm{SdS}}(1+\epsilon)$, with $\epsilon=10^{-3}$. Having the same qualitative features as 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_{\mathrm{\scriptscriptstyle 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 $\rho=0$, corresponding to infinite radii. When $\rho=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 $\rho$, 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\rightarrow\infty$, 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,

[TABLE]

Einstein equations reduce to,

[TABLE]

which can be formally rewritten as an autonomous system of two ODEs,

[TABLE]

Following the procedure described for example in jordan2007nonlinear , we introduce variables which compactify the phase space,

[TABLE]

Since $\varrho\geq 0$ and $\varrho^{2}+\varphi^{2}\leq 1$, with the last equality valid only in each of the limits $r\rightarrow\infty$ and $F\rightarrow\infty$, the whole phase space is compactified onto the unitary right half disk by stereographic projection.

Depending on the value of $\Lambda$, 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_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle H}}$ keeping $R_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle 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_{\mu\nu}$ and $g_{\mu\nu}$, i.e., $R_{\mathrm{\scriptscriptstyle 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 $\left(\rho,\bar{F},\bar{\sigma},\bar{R}\right)$. Using a color code for one dimension, we can make 3-dimensional plots. Here, the $\rho$ 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 $\left(0,0,0\right)$.

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 $\left(\bar{R},\bar{\sigma}\right)$ projections in which we are looking at the phase space “from the top”, the second column contains the $\left(\bar{R},\bar{F}\right)$ projections in which we are looking the phase space “from the front”.

The value of $\rho=1/r$ is indicated by the colour of the trajectories: the bluer the trajectory, the larger the $\rho$ (smaller $r$), the redder the trajectory, the smaller the $\rho$ (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 Figures7a and7b, we see that the Schwarzschild solution is Lyapunov unstable. In Figure7a 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 Figure7b, we see that the phase space flow is departing from the Schwarzschild solution also along the $\bar{F}$ axis when $r\rightarrow\infty$. 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_{\mathrm{\scriptscriptstyle 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 Figures7c and7e, 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 $\left(0,0,0\right)$ when $r\rightarrow\infty$ (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 Figures7d and7f, it is evident that SAdS and SdS solutions are unstable. Although trajectories get closer to $\left(0,0,0\right)$ in the $\left(\bar{R},\bar{\sigma}\right)$ plane, they are oscillating around the SAdS and SdS solutions with a growing amplitude along the $\bar{F}$ 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, $\kappa\to\infty$, 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 Volkov:2012wp . This means that perturbations oscillate around it with constant amplitude, not being asymptotically flat.

@fb@secFB

IV.3 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^{\mu}=(v,r,\theta,\phi)$ adapted for $g$. We can construct a complex null frame field (vielbein) of the metric $g$ which suits the ingoing null radial geodesics as,

[TABLE]

It is easy to verify that $g^{\mu\nu}=-2\ell_{g}^{(\mu}n_{g}^{\nu)}+2m_{g}^{(\mu}\bar{m}_{g}^{\nu)}$ 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,

[TABLE]

where the re-parametrisation invariance is ensured by the einbein $e(\lambda)$. Moreover, $e(\lambda)$ acts as a Lagrangian multiplier imposing the mass-shell constraint so that $\mu^{2}=1,0,-1$ corresponds to time-, null- or space-like geodesics, respectively. For the null radial geodesics, besides $\mu^{2}=0$, we have $\mathrm{d}\theta/\mathrm{d}\lambda=\mathrm{d}\phi/\mathrm{d}\lambda=0$. Note also that the vielbein (56) solely depends on the radial coordinate $x^{1}\equiv r$. Under these conditions, the variation of (57) with respect to $x^{\mu}(\lambda)$ and $e(\lambda)$ gives the following two geodetic equations, up to a suitably gauged $e(\lambda)$,

[TABLE]

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 $\ell_{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,

[TABLE]

for which $f^{\mu\nu}=-2\ell_{f}^{(\mu}n_{f}^{\nu)}+2m_{f}^{(\mu}\bar{m}_{f}^{\nu)}$ 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 $\ell_{f}$ and $n_{f}$, respectively, with the corresponding geodetic equations,

[TABLE]

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_{\star},r)=(v-r,r)$. In these coordinates, the null radial geodesics and the null cone fields for the improper bidiagonal solutions described in subsection IV.2, 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_{\star}$ 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 (Figures8b, 9b, 10b, 10c), since the translational Killing vector becomes null for both the metrics. Inside the horizon, both null cones shrink while $r\rightarrow 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_{\star}$, 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 Hassan:2017ugh , 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 (Figures8b, 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 (subsection III.2).

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=\sqrt{g^{-1}f}$ at finite $r$. These results are complementary to those in Volkov:2012wp 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=\sqrt{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 Brito:2013wya ; Babichev:2013una ; Babichev:2014oua ; Babichev:2015zub it was shown that the bidiagonal solutions are dynamically unstable, and as such cannot represent the end point of gravitational collapse. In addition, Babichev:2014oua showed that to have dynamically stable solutions, non-diagonal metric elements are needed, also at radii other then the horizon radius, $r_{\mathrm{\scriptscriptstyle 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.

Acknowledgements.

We want to thank Jonas Enander, Fawad Hassan, Bo Sundborg, Ingemar Bengtsson, Stefan Sjörs and Lars Bergström for helpful discussions. Our gratitude goes to Mikhail S. Volkov and Richard Brito for reading the paper and providing valuable remarks. We also thank Anders Lundkvist, Giacomo Monari, Julius Engelsöy, Sebastian Baum and Luca Visinelli for useful suggestions. Support for this study from the Swedish Research Council for EM is acknowledged.

Appendix A Old parametrisation for the metric fields

Given the coordinate system $x^{\mu}=(t,r,\theta,\phi)$, the most commonly used parametrisation of the metric fields for the bidiagonal ansatz is (see, for example, Volkov:2012wp ; Brito:2013xaa ; Babichev:2013pfa ; Enander:2015kda ),

[TABLE]

where $Q(r),N(r),a(r),Y(r),U(r)\in\mathbb{R}$. The determinants of these metrics are,

[TABLE]

The square root matrix $S=\sqrt{g^{-1}f}$, after having chosen the principal branch, is,

[TABLE]

The absolute values can be removed assuming that the fields are always positive. The determinant and the trace of the square root are,

[TABLE]

The Ricci scalars for the two metrics are,

[TABLE]

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 $\mathcal{K}^{\mu}=\delta^{\mu}_{0}$ is,

[TABLE]

In the coordinates adapted to the translational Killing vector, we cannot cross the Killing horizon since $\mathcal{K}$ becomes null at $r=r_{\mathrm{\scriptscriptstyle H}}$. On the other hand, using the Eddington-Finkelstein coordinates is problematic since $\mathcal{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\in[r_{\mathrm{\scriptscriptstyle H}},+\infty)$, which, in principle, allows specification of the initial values on the Killing horizon. Nonetheless, the numerical solutions found in the literature (see, for example, Volkov:2012wp ; Brito:2013xaa ; Babichev:2013pfa ; Enander:2015kda ) 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_{\mathrm{\scriptscriptstyle H}}$. Namely, in order to define a Killing horizon at $r_{\mathrm{\scriptscriptstyle H}}$, one must impose $Q(r_{\mathrm{\scriptscriptstyle H}})=a(r_{\mathrm{\scriptscriptstyle H}})=0$. In order for the determinants of the two metrics not to be zero at $r_{\mathrm{\scriptscriptstyle H}}$, we must also have $N(r_{\mathrm{\scriptscriptstyle H}})=Y(r_{\mathrm{\scriptscriptstyle H}})=0$, with $\lim\limits_{r\rightarrow r_{\mathrm{\scriptscriptstyle H}}}Q(r)/N(r)=\mathrm{const}$ and $\lim\limits_{r\rightarrow r_{\mathrm{\scriptscriptstyle H}}}a(r)/Y(r)=\mathrm{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 value $\bar{r}$, $Y(\bar{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=\bar{r}$, whose effects are studied in Gratia:2013gka ; Gratia:2013uza . 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_{\mu\nu}$ 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 subsection III.2.

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 $\Sigma(r)$, $\tau(r)$ and $R(r)$, separately.

We now consider an explicit example of a pathological behaviour, already understood as such in Brito:2013xaa . 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_{\mathrm{\scriptscriptstyle H}}=1,+\infty)$ and has $\beta_{1}=1,\beta_{2}=\beta_{3}=0$, $\beta_{0}$ and $\beta_{4}$ determined by the asymptotically flatness conditions in eq. (33) (so $c=1$) and $m_{\mathrm{g}}=0.04$ ($m_{\mathrm{g}}$ is denoted with $\mu$ in Brito:2013xaa ; the Compton wavelength of the massive graviton is then $\lambda=25$ in geometrical units with the Planck constant $h=1$). Hence $m_{\mathrm{g}}m_{\mathrm{S}}=m_{\mathrm{g}}r_{\mathrm{\scriptscriptstyle H}}/2=0.02$, with $m_{\mathrm{S}}$ Schwarzschild mass of the BH. The second panel of Figure 1 in (Brito:2013xaa, , p. 4) tells us what is the value of $U(r_{\mathrm{\scriptscriptstyle H}})=u_{\mathrm{\scriptscriptstyle H}}$ to choose as initial condition, given the value of $m_{\mathrm{g}}m_{\mathrm{S}}$. Note that, since $r_{\mathrm{\scriptscriptstyle H}}=1$ both in Brito:2013xaa and in our approach, $u_{\mathrm{\scriptscriptstyle H}}$ corresponds to our $R_{\mathrm{\scriptscriptstyle H}}$. For $m_{\mathrm{g}}m_{\mathrm{S}}=0.02$, $u_{\mathrm{\scriptscriptstyle H}}\approx 50.2$.

In Figure 11, all relevant quantities for this solution are plotted: the metric functions and the first derivative $U^{\prime}(r)$, which appears explicitly in $f_{\mu\nu}$, $\det S$ and the Ricci and Kretschmann scalars for $f_{\mu\nu}$. The curvature scalars for $g_{\mu\nu}$ 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^{\prime}(r)$ which are plotted in the second panel. The solid grey lines correspond to the following 4 points: from left to right, (i) $U^{\prime}(r)=0$, (ii) $a(r)=0$, (iii) $U^{\prime}(r)=0$, (iv) $a(r)=0$. There are also large-dashed grey lines, highlighting the points at which $U^{\prime}(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 choquet2008general . 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^{\mu}=(t,\xi,\theta,\phi)$, the metrics $g$, $f$ and the square root matrix $S$ are given by,

[TABLE]

The determinants of the metrics are,

[TABLE]

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,

[TABLE]

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 coordinates $\bar{x}^{\beta}=\left(v,\xi,\theta,\phi\right)$ for $g$ defined by the following Jacobian,

[TABLE]

where $x^{\alpha}=\left(t,\xi,\theta,\phi\right)$ are the old Schwarzschild coordinates. This defines the ingoing null coordinate $v$ for $g$ (also known as ‘advanced time’),

[TABLE]

Applying the coordinate transformation (72) to the metrics (68) 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 $\mathrm{d}v$ because we do not know analytically the functions $q(\xi)$ and $F(\xi)$. Nonetheless, we can still define the tortoise coordinate in an implicit form,

[TABLE]

where $\mathrm{d}\xi=\mathrm{d}r$. Note that $\xi^{*}$ is a monotonic function of $\xi$ separately inside and outside the event horizon (or, for a SdS solution, between the horizons, inside the event horizon and outside the cosmological horizon),

[TABLE]

Appendix C On bidiagonality condition in [44]

In the following we make Proposition 1 in Deffayet:2011rh more precise. Consider a coordinate patch where two static and spherically symmetric metrics take the form,

[TABLE]

Suppose further that the Killing vector $\mathcal{K}=\partial_{t}$ is null with respect to $g$ at $r=r_{\mathrm{\scriptscriptstyle H}}$. Proposition 1 in Deffayet:2011rh states that, if both metrics describe smooth geometries and are diagonal at $r=r_{\mathrm{\scriptscriptstyle H}}$, the Killing vector $\mathcal{K}$ must also be null with respect to $f$ at $r=r_{\mathrm{\scriptscriptstyle H}}$. The proof presented in Deffayet:2011rh is elegant in its simplicity. It starts from the presuppositions $A(r_{\mathrm{\scriptscriptstyle H}})=0$, $B(r)=0$, and that the trace $\operatorname{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 $\mathcal{K}$ is timelike for $r>r_{\mathrm{\scriptscriptstyle H}}$. Therefore $J(r_{\mathrm{\scriptscriptstyle H}})=0$, otherwise $J/A$ would diverge at $r=r_{\mathrm{\scriptscriptstyle 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_{\mathrm{\scriptscriptstyle H}}$, we introduce the ingoing Eddington-Finkelstein coordinates $(v,r)$ adapted for $g$, by $\mathrm{d}t^{2}=\mathrm{d}v^{2}-2\sqrt{C/A}\,\mathrm{d}v\mathrm{d}r+(C/A)\mathrm{d}r^{2}$ at $r>r_{\mathrm{\scriptscriptstyle H}}$. Assuming Lorentzian signature of the metrics, we have (up to square root signs, not shown for readability),

[TABLE]

Clearly, the metrics are regular at $r=r_{\mathrm{\scriptscriptstyle H}}$ in this coordinate system. Since both $J/A$ and $K/C$ are finite and non-zero, it is also $c\coloneqq(J/A)(C/K)$ and we can write,

[TABLE]

At $r=r_{\mathrm{\scriptscriptstyle H}}$ we have $A(r_{\mathrm{\scriptscriptstyle H}})=J(r_{\mathrm{\scriptscriptstyle H}})=0$; thus,

[TABLE]

We can apply the theorem on canonical pair forms Uhlig:1973 on (80). Thus, the metrics which describe smooth geometries cannot be both diagonal at $r=r_{\mathrm{\scriptscriptstyle H}}$ unless $K-\frac{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_{\mathrm{\scriptscriptstyle H}}$ (this corresponds to the proper bidiagonality condition in Proposition 3). In other words, strictly, Proposition 1 in Deffayet:2011rh 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\neq K-\frac{J}{A}C<\infty$ (crossing condition in Proposition 2), i.e., when the metrics cannot be simultaneously diagonalised at $r_{\mathrm{\scriptscriptstyle H}}\in\{r\,|\,A(r)=0\}$ but they are bidiagonal elsewhere by assumption.

Appendix D Numerical details

In this appendix, we describe in more detail the numerical procedure used to obtain the solutions in subsubsection IV.2.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 mathematica . In particular, NDSolve`LSODA method was used, which is based on LSODA ref60 ; Brown:1989:RSM:2613951.2613987 . 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 Enright:1989:NEI:2613951.2614000 ; Shampine:2005:SOD:1052070.1052077 . 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 Enright:1989:NEI:2613951.2614000 . Suppose we have a first-order system of ODEs,

[TABLE]

with $y,f\in\mathbb{R}^{n}$ and $n\in\mathbb{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,

[TABLE]

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,

[TABLE]

is close to the original one or not. This is equivalent to asking if residuals are small compared to $f(x,u(x))$. In Enright:1989:NEI:2613951.2614000 ; Shampine:2005:SOD:1052070.1052077 , it is pointed out that, if $u(x)$ is the approximation of $y(x)$ in a given interval of the independent variable $x$, then,

[TABLE]

We will use eqs. (82)–(84) 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 $\mathscr{F}_{F}$ in eq. (26) for all solutions presented in this paper. This means, according to eq. (83), 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. (82), but also eq. (84). They are plotted in Figure 14, where the red and blue curves are almost coinciding over the whole range of integration. These two curves are the residuals calculated using eqs. (82) (red) and (84) (blue). Being equal, our numerical GR solutions are approximating the analytical GR solutions.

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. (84), after some finite $r$ the difference $|{F^{\prime}}_{\mathrm{pert}}^{\mathrm{num}}-{F^{\prime}}_{\mathrm{GR}}^{\mathrm{num}}|$ must be due only to numerical residuals. Conversely, if this difference is not close to the residuals calculated through eq. (82), 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^{\prime}}_{\mathrm{pert}}^{\mathrm{num}}-{F^{\prime}}_{\mathrm{GR}}^{\mathrm{num}}|$ and $|{F^{\prime}}_{\mathrm{pert}}^{\mathrm{num}}-\mathscr{F}_{F}[\mathscr{V}_{\mathrm{pert}}^{\mathrm{num}}]|$ 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^{\prime}}_{\mathrm{pert}}^{\mathrm{num}}-{F^{\prime}}_{\mathrm{GR}}^{\mathrm{num}}|$ (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^{\prime}}_{\mathrm{pert}}^{\mathrm{num}}-{F^{\prime}}_{\mathrm{GR}}^{\mathrm{num}}|$.