# Hamiltonian unboundedness vs stability with an application to Horndeski theory

Eugeny Babichev, Christos Charmousis, Gilles Esposito-Farèse, Antoine Lehébel

► To cite this version:

Eugeny Babichev, Christos Charmousis, Gilles Esposito-Farèse, Antoine Lehébel. Hamiltonian unboundedness vs stability with an application to Horndeski theory. Physical Review D, 2018, 98 (10), pp.104050. 10.1103/PhysRevD.98.104050. hal-01768013

HAL Id: hal-01768013

https://hal.science/hal-01768013v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hamiltonian unboundedness vs stability with an application to Horndeski theory

E. Babichev,<fnref n="1" />,²,* C. Charmousis,<fnref n="1" />,† G. Esposito-Farèse,²,‡ and A. Lehébel<fnref n="1" />,§

¹Laboratoire de Physique Théorique, CNRS, Université Paris-Sud,
Université Paris-Saclay, F-91405 Orsay, France

²Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris,
GRÉCO, 98bis boulevard Arago, F-75014 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 10 September 2018; published 28 November 2018)

A Hamiltonian density bounded from below implies that the lowest-energy state is stable. We point out, contrary to common lore, that an unbounded Hamiltonian density does not necessarily imply an instability: Stability is indeed a coordinate-independent property, whereas the Hamiltonian density does depend on the choice of coordinates. We discuss in detail the relation between the two, starting from k-essence and extending our discussion to general field theories. We give the correct stability criterion, using the relative orientation of the causal cones for all propagating degrees of freedom. We then apply this criterion to an exact Schwarzschild-de Sitter solution of a beyond-Horndeski theory, while taking into account the recent experimental constraint regarding the speed of gravitational waves. We extract the spin-2 and spin-0 causal cones by analyzing respectively all the odd-parity and the ℓ = 0 even-parity modes. Contrary to a claim in the literature, we prove that this solution does not exhibit any kinetic instability for a given range of parameters defining the theory.

DOI: 10.1103/PhysRevD.98.104050

## I. INTRODUCTION

General relativity (GR) is an effective classical theory of gravity which is experimentally verified for a wide range of physical distance and gravitational strength scales. The former range up to 30 astronomical units or so, while the latter are performed from weak gravity tabletop experiments to very strong gravity environments of coalescing black hole binaries, involving recent gravitational wave detections (see [1] for a review, and for the latter [2,3]). GR is furthermore a theoretically robust theory, as it is unique when one imposes standard mathematical and physical assumptions: essentially the presence of a Levi-Civita connection for a sufficiently regular four-dimensional spacetime manifold equipped with a metric mediating gravitational interactions. These dictate that GR with a cosmological constant is the unique covariant theory with second-order field equations [4]. Since no mass term is associated with the rank-2 tensor mediating gravity, GR has two massless spin-2 degrees of freedom (d.o.f.). A second and independent assumption of GR is that all matter fields universally couple to this metric, in order to satisfy the weak equivalence principle. These postulates imply that

photons and gravitons propagate on the same causal cones, i.e., with the same speed.

The presence of dark energy at cosmological distance scales (but also the yet elusive dark matter) has opened up in recent years the possibility that GR may be an effective theory not only at UV scales but also at the deep IR: cosmological—but also lack of astrophysical—observations have raised questions concerning the viability of GR at large distances (see e.g., [5]). Given the aforementioned uniqueness, if we want to go beyond GR, we have to introduce novel d.o.f. This is true if we remain within the realm of Riemannian geometry, keeping the postulated geometrical structure of spacetime. The simplest of these d.o.f. is a scalar field, giving scalar-tensor theories of gravity (see [6–15]), but one can also consider vector(s) (see e.g., [16–20]), an additional metric field (see [21] for a review), etc. One should point out that even more general considerations, such as the presence of nontrivial torsion (e.g., [22]) or spacetime with extra dimensions (see [23] and references therein), also lead effectively to the addition of extra d.o.f.. The prototype example is that of the Dvali-Gabadadze-Porrati (DGP) five-dimensional<fnref n="1" /> braneworld model [25], which upon going to the decoupling limit [26] reduces to a particular Horndeski scalar-tensor theory (of the type studied in [27]).

<footnote>*eugeny.babichev@th.u-psud.fr</footnote>

<footnote>†christos.charmousis@th.u-psud.fr</footnote>

<footnote>‡gef@iap.fr</footnote>

<footnote>§antoine.lehebel@th.u-psud.fr</footnote>

<footnote n="1">¹See also [24] for relations between higher dimensional and 4-dimensional scalar-tensor gravity theories.</footnote>
Modified gravity d.o.f.—including gravitons—propagate in an effective metric which can be different from that of GR. The notion of causal cones and effective metric will be essential to our stability arguments, so let us dwell on this point before entering details in the body of the paper. Consider some background solution of a modified gravity theory. Linear perturbations of the background can be found by expanding the action up to the second order. According to their effective action, modes obey some second-order differential equation.<fnref n="2" /> We will concentrate on the kinetic part of the equations of motion (i.e., the principal part of the differential equation), since this part defines whether the most dangerous pathologies of the theory are absent. The kinetic operator is encoded in an effective metric, and in general it depends on the background solution. Provided that the equation of motion is hyperbolic, the kinetic operator (i.e., the effective metric) defines then a causal cone of propagation associated to the perturbative d.o.f. This causal cone is inherently different for different helicities—scalar, vector, or tensor—and its structure determines whether the modes are healthy or not.

Furthermore, matter is assumed to couple universally to a single metric in order to pass stringent fifth-force experiments. This introduces the matter causal cone, in addition to gravity cones, and the physical metric, associated to geodesic free-fall, which matter couples to. The prototype example is Brans-Dicke (BD) scalar-tensor theory (see e.g., [44]). There, the Jordan frame is the physical frame, whose associated metric gives geodesic free fall. The more calculation-friendly Einstein frame is related to the physical frame *via* a specific conformal transformation involving the metric and the scalar d.o.f.. Note that for Horndeski theory and beyond, modes of different species (in our example scalar and tensor) mix together and it is in general only for the most symmetric backgrounds that one manages to demix them.

The effect of multiple causal cones and mixings is that species—such as gravitons—can now have subluminal or superluminal propagation with respect to photons (matter light cone). However, these multiple possibilities have been recently constrained by a single observation. The simultaneous detection of gravity wave event GW170817 [2] and light emanating from the same source [45] at 40 Mpc distance strongly restricts the graviton (spin 2) causal cone and that of light to be essentially identical, just like in GR. This restricts the variety of modified gravity theories to a subclass of theories with gravitons propagating at the speed of light, i.e., such that $c_g = 1$ in the physical frame [46–52]. On the other hand, in some cases, it is technically easier to

work in a nonphysical frame, where the metric is disformally related to the physical metric. For Horndeski theory, e.g., one can find a disformal transformation of the metric that brings it to a unit propagation speed [46] ($c_g = 1$) theory of EST/DHOST type [38,39]. Note that this is not possible for any Horndeski theory; e.g., $G_5$ interactions involving the Gauss-Bonnet term are excluded [46,47].<fnref n="3" /> In a recent paper [53], considering a vacuum black hole background, we showed how the right choice of disformal transformation can ensure that the graviton speed is the same as that of light. Although we started from a shift-symmetric $G_4$ Horndeski theory, which is excluded in the physical frame, there exists a specific disformal transformation [46] upon which the graviton causal cone is identical to that of light as demanded by observation. The physical metric is a disformed metric of the Horndeski action, which mixes scalar and metric perturbations and brings us to a strictly $c_g = 1$ theory. In this case the initial Horndeski theory plays the role of the Einstein frame, whereas the target EST/DHOST theory [38,39] plays the role of the Jordan frame, in analogy to the familiar BD theory cited above.

In this paper we will see how, starting from the causal cone structure of propagating d.o.f., one can infer if the modes in question are healthy modes—in other words that they are not modes generating ghost or gradient instabilities. In particular, the sign of the determinant of the effective metric defines the hyperbolicity condition, which if satisfied, means that a particular solution is safe from imaginary speeds of propagation and therefore gradient instabilities. On the other hand, the local orientation of the cone tells us about absence/presence of ghost modes. Both requirements allow the local definition of a causal cone of propagation which then guarantees a healthy associated mode.

A complementary way to find the good or sick nature of propagating modes is often described *via* the associated Hamiltonian density of the modes in question. Once the effective action for the mode is known, one defines the conjugate momentum and writes down the Hamiltonian density of the associated field. It is known that if the Hamiltonian density is bounded from below, then the ground state is of finite energy and necessarily stable. The contrary is often assumed to be true: If a Hamiltonian is unbounded from below, then the system is unstable and admits ghost or gradient instabilities.

One of the main aims of our paper is to explicitly show that the above inverse statement is not always true. In other words, if a Hamiltonian density is unbounded from below,

<footnote n="2">²We do not consider here Lorentz-breaking theories [28,29] where equations of motion can be of higher order. Also, in extensions of Horndeski theory [30–43], one may get Euler-Lagrange equations of third order, but by manipulating these equations their order is reduced.</footnote>

<footnote n="3">³By excluded, throughout this paper, we refer to the cases where the extra mode is a dark energy field, giving an effective acceleration to the universe at late times. If the extra mode, say a scalar, is not varying at cosmological scales, but only locally, it may not influence gravitational waves in their 40 Mpc journey to Earth detectors.</footnote>
this does not necessarily signify that the mode in question generates a ghost or gradient instability. The reasoning is simple although it goes against standard lore originating from particle physics or highly symmetric backgrounds associated to Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology. The Hamiltonian is not a scalar quantity and therefore depends on the coordinate system it is associated with. As such we will explicitly see that Hamiltonian densities can be unbounded by below but under a coordinate transformation can be transformed to a bounded density. The key point will be the coordinate system on which the Hamiltonian is to be defined in relation to the effective causal cones.

In fact, we will see that the coordinate system will have to be of a certain “good” type in order for the Hamiltonian density to be conclusive. For our purposes, we will restrict ourselves to configurations where essentially the problem is mathematically 2-dimensional and involves the definition of a good timelike and spacelike direction. This includes the case for planar, cylindrical or spherical symmetry, for example. A “good” coordinate system will involve the existence of a common timelike direction for all causal cones. Second, it will involve the existence of a common spacelike direction exterior to all causal cones.<fnref n="4" /> If such a coordinate system exists, then we will show that the Hamiltonian density is bounded from below and the system is stable. If such a coordinate system does not exist, on the contrary, then the Hamiltonian density is always unbounded from below. The relevant criteria emerging from the causal cones will inevitably lead to the knowledge of ghost or gradient instabilities present in the system.

We will explicitly show all this for a general situation with two propagating d.o.f. with different causal structures in Sec. II, and apply it to the known stability criteria of k-essence. We will see explicitly how stability criteria are satisfied for well-defined causal cones and how on the contrary, the unboundedness of Hamiltonian densities can lead to wrong conclusions if not associated to a “good” coordinate system. We will then move on, in Sec. III, to apply our causal cone criteria to investigate the stability of a specific Horndeski theory admitting a nontrivial background black-hole solution [54]. The family of strongly gravitating solutions admit a time-dependent scalar field which is asymptotically a dark energy field allowing de Sitter acceleration inherently different from the vacuum cosmological constant. The mixed combination of space and time dependence for the scalar, as well as the higher order nature of the theory, leads to causal scalar and tensor cones which are quite complex. This is the reason why the

Hamiltonian analysis of [55] gave the wrong conclusion,<fnref n="5" /> stating the generic instability of such black holes for any coupling constants of the theory. Although the Hamiltonian associated with the graviton is unbounded by below in Schwarzschild coordinates, we will see that it is bounded by below in an appropriate coordinate system. The graviton and matter causal cones indeed keep compatible orientations, and can actually be chosen to exactly coincide at any spacetime point in order to satisfy the gravity speed constraint imposed by the GW170817 event. We will additionally complete this analysis by deriving the scalar causal cone, and showing that it also has a compatible orientation with the two previous cones for a certain range of parameters of the model. We compute this last cone by studying the $\ell = 0$ perturbations. We will conclude in Sec. IV.

## II. CAUSAL CONES AND HAMILTONIAN IN A GENERAL COORDINATE SYSTEM

Standard theories minimally coupling all fields to one metric tensor $g_{\mu\nu}$ possess a single causal cone, defined by $ds^2 = g_{\mu\nu}dx^\mu dx^\nu = 0$, or equivalently by $g^{\mu\nu}k_\mu k_\nu = 0$ for a wave vector $k_\mu$, $g^{\mu\nu}$ denoting as usual the inverse of $g_{\mu\nu}$. This is no longer the case when at least two fields are coupled to different metric tensors which are not proportional—even though they are generally related to each others. For instance, all Galileon models [aside from the simplest $(\partial_\mu\varphi)^2$ Lagrangian for a scalar field $\varphi$] and their generalizations called beyond-Horndeski theories [8–15,30–41,60], predict that the spin-2 and spin-0 d.o.f. propagate in different effective metrics, which depend on the background solution. And these two effective metrics actually also generically differ from $g_{\mu\nu}$, to which one may (or may not) choose that matter fields are universally coupled.

The simplest example is k-essence [61–64], i.e., a Lagrangian for the scalar field given by a nonlinear function $f$ of the standard kinetic term, $\mathcal{L} = -\frac{1}{4}f(X)$, where $X \equiv g^{\mu\nu}\partial_\mu\varphi\partial_\nu\varphi$. If one writes the scalar field as $\varphi = \bar{\varphi} + \chi$, where $\bar{\varphi}$ denotes the background solution and $\chi$ a small perturbation, one finds that the second-order expansion of this Lagrangian reads $\mathcal{L}_2 = -\frac{1}{2}S^{\mu\nu}\partial_\mu\chi\partial_\nu\chi$, where

$$S^{\mu\nu} = f'(\bar{X})g^{\mu\nu} + 2f''(\bar{X})\nabla^{\mu}\bar{\varphi}\nabla^{\nu}\bar{\varphi}, \quad (1)$$

$f'$ and $f''$ being the first and second derivatives of function $f$ with respect to its argument $\bar{X}$ [62,65–68]. This means

<footnote n="4">⁴A causal cone represents an open set whose interior is bounded by the characteristics of the cone. The complementary of this set with boundary is an open set which is the exterior of the cone.</footnote>

<footnote n="5">⁵References [56–58] used similar arguments, and accordingly obtain too restrictive conditions for stability. Reference [59] also uses these arguments, but it proves the stability of the odd-parity modes outside neutron stars, and this is correct.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Possible relative orientations of two causal cones, in a coordinate system such that the grey cone with solid lines appears at ±45°. We do not plot the equivalent configurations exchanging left and right, and do not consider the limiting cases where some characteristics coincide. The first row (a)–(d) are safe cases in which the two metrics can be diagonalized simultaneously by an appropriate choice of coordinates—corresponding then to panels (b) or (c). Although the kinetic contribution to their Hamiltonian density is unbounded by below in cases (a) and (d), it is positive in (b) and (c). The second row (e)–(h) are again safe cases, for which the kinetic contribution to the Hamiltonian density can be proven to be positive in an appropriate coordinate system, actually corresponding to case (e), but the two metrics cannot be simultaneously diagonalized. [Let us recall that two quadratic forms can always be simultaneously diagonalized if at least one of them is positive (or negative) definite. Here both of our metrics have hyperbolic signature, and this is the reason why the nonsimultaneously diagonalizable cases (e)–(h) are possible.] The third row (i)–(l) are unstable cases, for which the two metrics can be simultaneously diagonalized as in (j) and (k), but they have then opposite signatures in this (t, x) subspace. Their total Hamiltonian density remains unbounded by below in all coordinate systems.</figcaption></figure>

that the spin-0 d.o.f. $\chi$ propagates in an effective metric$^6$ $S^{\mu\nu}$, which is not proportional to $g^{\mu\nu}$ as soon as $f''(\bar{X}) \neq 0$ and the background solution has a nonvanishing gradient $\partial_\mu \bar{\phi}$, and they define thus different causal cones. In this simple case, one can show that the spin-2 d.o.f. (the gravitons, perturbations of the metric tensor) do propagate in the initial metric $g^{\mu\nu}$. To simplify this example even further, one may actually consider it in flat spacetime, i.e., without any graviton, while universally coupling matter to $g_{\mu\nu}$. Then we still have at least two fields (matter and the

spin-0 d.o.f. $\chi$) which propagate in different metrics, defining two different causal cones.

The conditions for such a k-essence theory to be stable and have a well-posed Cauchy problem have been written several times in the literature [62,65–69], and we shall rederive them at the end of the present section from our general analysis. They read$^7$ $f'(\bar{X}) > 0$ and $2\bar{X}f''(\bar{X}) + f'(\bar{X}) > 0$. When the background scalar gradient $\partial_\mu \bar{\phi}$ is timelike with respect to $g^{\mu\nu}$, the causal cones can be represented as panels (a), (b), (c), or (d) of Fig. 1, where the grey cone (with solid lines) is defined by $g^{\mu\nu}$ and

<footnote>$^6$We use the notation $S^{\mu\nu}$ for the effective metric in which the Scalar d.o.f. propagates, while $G^{\mu\nu}$ will be used in Sec. III A to denote the effective metric in which spin-2 d.o.f. (Gravitons) propagate.</footnote>

<footnote>$^7$We choose the mostly-plus signature convention for the metric $g_{\mu\nu}$.</footnote>
the blue one (with dashed lines) by $S^{\mu\nu}$. Panel (a) is actually transformed into (b), and (d) into (c), if one chooses a coordinate system such that the spatial gradients $\partial_i \bar{\phi}$ vanish, the vector $\partial_\mu \bar{\phi}$ pointing then exactly in the time direction.

Panels (a) and (b) correspond to $f''(\bar{X}) < 0$, and mean that the spin-0 d.o.f. $\chi$ propagates slower than light (which is the fastest matter field). Panels (c) and (d) correspond to $f''(\bar{X}) > 0$, and describe a superluminal scalar field, but this does not lead to any causality problem as soon as this dashed cone remains always a cone, with a nonempty exterior where one may define Cauchy surfaces to specify initial data. This has already been discussed in detail in the literature [65–69]. Paradoxes only occur when one wants to specify initial data on the $t=0$ hypersurface in the situation of panel (d): This is forbidden because this hypersurface is not spacelike with respect to the dashed cone. Note that panels (b) and (c) are actually equivalent if one exchanges the meaning of the colors. If one chooses a coordinate system such that the scalar causal cone is at $\pm 45^\circ$ (grey cone with solid lines), then matter propagates within the dashed blue cone, and the case of a superluminal scalar perturbation $\chi$ now corresponds to panels (a)–(b). Causality becomes then more obvious than in panel (d).

Independently of the specific form (1) taken by the effective metric $S^{\mu\nu}$ in the case of k-essence, let us now consider any possible $S^{\mu\nu}$ in which a field $\chi$ propagates, to discuss all the cases of Fig. 1. To simplify, we shall assume that the standard metric $g_{\mu\nu}$ (to which matter is assumed to be universally coupled) is flat. If the Lagrangian defining the dynamics of $\chi$ reads as before $\mathcal{L}_2 = -\frac{1}{2}S^{\mu\nu}\partial_\mu\chi\partial_\nu\chi$, where we focus only on the kinetic term, then the conjugate momentum is defined as

$$p = \frac{\partial \mathcal{L}_2}{\partial \dot{\chi}} = -S^{00} \dot{\chi} - S^{0i} \partial_i \chi, \quad (2)$$

and the contribution of this field $\chi$ to the Hamiltonian density reads

$$\begin{aligned} \mathcal{H}_2 &= p \dot{\chi} - \mathcal{L}_2 \\ &= -\frac{1}{2S^{00}}(p + S^{0i} \partial_i \chi)^2 + \frac{1}{2} S^{ij} \partial_i \chi \partial_j \chi. \end{aligned} \quad (3)$$

Note that its positiveness depends only on $S^{00}$ and $S^{ij}$, but not on the mixed components $S^{0i}$, although we shall see that they are actually crucial for the stability analysis. Stability is indeed a physical (observable) statement, which should be coordinate independent, whereas the Hamiltonian density is not a scalar and depends thus on the coordinate system.

To simplify even further the discussion, let us assume that $S^{\mu\nu}$ is of the form

$$\begin{pmatrix} S^{00} & S^{01} & 0 & 0 \\ S^{01} & S^{11} & 0 & 0 \\ 0 & 0 & S^{22} & 0 \\ 0 & 0 & 0 & S^{33} \end{pmatrix}, \quad (4)$$

with $S^{22} \ge 0$ and $S^{33} \ge 0$, and let us focus on the $(t, x)$ subspace as in Fig. 1. [In the neighborhood of a spherical body, for instance, it is natural to choose spherical coordinates where $S^{0\theta} = 1/r^2$ and $S^{\phi\phi} = 1/(r^2 \sin^2 \theta)$ or similar in generalized Galilean and beyond-Horndeski theories, the difficulties being restricted to the $(t, r)$ subspace.] In order for this metric to define a cone, with nonempty interior and exterior, it is necessary that its determinant be negative:

$$D = S^{00}S^{11} - (S^{01})^2 < 0. \quad (5)$$

Note that this hyperbolicity condition does depend on the off-diagonal component $S^{01}$, contrary to the sign of Hamiltonian (3) above. The inverse $S_{\mu\nu}^{-1}$ of matrix (4) [its exponent $-1$ being explicitly written in order not to confuse it with $g_{\mu\lambda}g_{\nu\rho}S^{\lambda\rho}$] reads in the $(t, x)$ subspace

$$\begin{pmatrix} S^{11} & -S^{01} \\ -S^{01} & S^{00} \end{pmatrix} / D. \quad (6)$$

We can thus conclude that when $S^{\mu\nu}$ indeed defines a cone, then $S_{00}^{-1}$ has the opposite sign of $S^{11}$, and $S_{11}^{-1}$ the opposite sign of $S^{00}$.

Let us now consider the various cone orientations of Fig. 1. In the situation of panel (a), the time axis is outside the dashed (blue) cone defined by $S^{\mu\nu}$. This means that $S_{00}^{-1} dtdt > 0$, and therefore $S^{11} < 0$. This implies that the Hamiltonian density (3) is unbounded by below because of the contribution of $\frac{1}{2}S^{ij}\partial_i\chi\partial_j\chi$, when $\partial_i\chi$ is large enough (and $p$ is chosen to compensate $S^{0i}\partial_i\chi$). This conclusion remains the same for all panels of this Figure in which the time axis is outside the dashed cone, namely (f), (h), (i), (j), and (k). On the contrary, when the time axis is within the dashed cone (in all other panels of Fig. 1), this corresponds to $S^{11} > 0$, and the second term of Hamiltonian (3) is thus positive.

Similarly, in the situation of panel (d), the $x$ axis is within the dashed cone, therefore $S_{11}^{-1}dxdx < 0$, which implies $S^{00} > 0$. In this case, the Hamiltonian density (3) is unbounded by below because of the contribution of its first term $-(p + S^{0i}\partial_i\chi)^2/(2S^{00})$. This conclusion remains the same for all panels in which the $x$ axis is inside the dashed cone, namely (g), (h), (j), (k), and (l). In all other panels, the $x$ axis is outside the dashed cone, therefore $S^{00} < 0$ and the first term of Hamiltonian (3) is thus positive.
Note that panels (h), (j) and (k) have both their time axis outside the dashed cone and their x axis within it. This means that the Hamiltonian density (3) is always negative, while that corresponding to matter (coupled to $g_{\mu\nu}$ and propagating thereby in the solid grey cone) is always positive. It thus naively seems that any coupling between matter and $\chi$, or any indirect coupling via another field (for instance gravity), will lead to deadly instabilities. This is indeed the case for panels (j) and (k), but not for panel (h). Indeed, if one chooses another coordinate system such that the new time $t'$ lies within the intersection of both cones (superposition of the grey and blue regions), and the new spatial direction $x'$ is outside both cones (white region), then one gets simultaneously the four conditions $g'^{00} < 0$, $g'^{11} > 0$, $S'^{00} < 0$, and $S'^{11} > 0$. Therefore, both the Hamiltonian density (3) for the spin-0 d.o.f. $\chi$ and its analogue for matter are positive in this coordinate system. This suffices to prove that no instability can be caused by the kinetic terms in the situation of panel (h).

Let us recall that when a total Hamiltonian density (including all interacting fields) is bounded by below, then the lowest-energy state is necessarily stable. It is indeed impossible to reach a higher energy state (for any field) without violating energy conservation. But note that the converse theorem does not exist, as underlined by the reasoning above: A Hamiltonian density which is unbounded by below does not always imply an instability. In panel (h) of Fig. 1, this Hamiltonian was the sum of the positive contribution of matter and of the (always) negative contribution of the spin-0 field $\chi$, but we saw that there exist other coordinate systems in which both contributions are simultaneously positive.

To understand this better, let us just consider the boosts of special relativity in flat spacetime, instead of the arbitrary coordinate transformations allowed in GR. Then the metric $g^{\mu\nu}$ in which matter propagates always reads diag(-1, 1, 1, 1), and it defines the solid (grey) cones of Fig. 1. In the simple cases of panels (b) and (c), the components of $S^{\mu\nu}$ in the $(t, x)$ subspace read $k^2\text{diag}(-1/c_s^2, 1)$, where $k$ is a nonvanishing constant and $c_s$ is the velocity corresponding to the characteristics of the dashed (blue) cone.<fnref n="8" /> Indeed, the wave equation for the spin-0 field $\chi$ reads $S^{\mu\nu}\partial_\mu\partial_\nu\chi = 0$, and it admits as solutions arbitrary functions of $(x \pm c_s t)$. Panel (b) corresponds to $c_s^2 < 1$ while panel (c) to $c_s^2 > 1$. If we now perform a boost of velocity $-v$, we find that the components of $S'^{\mu\nu}$ in the new coordinate system read

$$ \frac{k^2}{c_s^2(1-v^2)} \begin{pmatrix} -1 + v^2 c_s^2 & v(1-c_s^2) \\ v(1-c_s^2) & c_s^2 - v^2 \end{pmatrix}. \quad (7) $$

We thus immediately see that $S'^{11} < 0$ (with $S'^{00}$ still negative) when we choose $|c_s| < |v| < 1$ in the case of panel (b), i.e., that we obtain the situation of panel (a), as described below Eq. (6). Although we started from the stable situation of panel (b), in which the total Hamiltonian density is positive, we thus find that the contribution of the spin-0 d.o.f. is no longer bounded by below in this boosted frame corresponding to panel (a). This is the main lesson: The unboundedness by below of the Hamiltonian density is a mere coordinate effect in the present situation, and it has no physical meaning. The model is stable, but one is not computing the “right” quantity in the boosted frame of panel (a). [We shall come back to this “right” quantity below.]

Note that a negative value of $S'^{11}$ in the boosted frame of panel (a) always comes together with a significant nonzero value of $|S'^{01}| = |S'^{10}| > \sqrt{-D}$, where $D$ is the determinant (5). The reason is that this determinant must remain negative in all coordinate systems—and actually remains strictly equal to $D$ when one considers only special-relativistic boosts as here. These nonzero off-diagonal components of $S'^{\mu\nu}$ are crucial for the existence of an inverse boost taking us back to the situation of panel (b), where the total Hamiltonian density is positive. If they were absent, then the metric diag($S'^{00}, S'^{11}$) would be negative definite, it would not define any causal cone, and the Cauchy problem would be ill-posed. Note also that the magnitude of these off-diagonal components of $S'^{\mu\nu}$ is also crucial. For instance, panel (i) of Fig. 1 corresponds to $S'^{00} < 0$ and $S'^{11} < 0$ like panel (a), and it does satisfy $|S'^{01}| > \sqrt{-D}$, but also the inequality $|S'^{01}| < |S'^{00} + S'^{11}|/2$ which leads to the situation of panel (j) when diagonalizing $S^{\mu\nu}$ by an appropriate boost. In this case (j), the two metrics $g^{\mu\nu}$ and $S^{\mu\nu}$ have opposite signatures in the $(t, x)$ subspace, so that the spin-0 d.o.f. $\chi$ behaves as a ghost in this subspace, and the model is unstable as soon as $\chi$ is somehow coupled to matter (including indirectly, e.g., via gravity).

Let us now apply a boost to the case of panel (c) of Fig. 1. If we choose $|c_s|^{-1} < |v| < 1$, then we find from Eq. (7) that $S'^{00} > 0$ (with $S'^{11}$ still positive), i.e., we obtain the situation of panel (d). Here again, as described above, we thus find that the contribution of the spin-0 d.o.f. to the Hamiltonian density is no longer bounded by below in this boosted frame, whereas is was positive in the initial frame corresponding to panel (c). [The fact that the first term of (3), proportional to $\dot{\chi}^2$, becomes negative is related to the wrong time-orientation<fnref n="9" /> of the null vector $N^\mu$ with respect

<footnote n="8">⁸We set $c = 1$ for the velocity of light, corresponding to the solid (grey) cone.</footnote>

<footnote n="9">⁹On the other hand, the possible negative values of Hamiltonian (3) in the previous case of panel (a) is less obvious, since the null vectors $N^\mu$ (with respect to $S'^{-1}_{\mu\nu}$) always remain future-oriented. In that case, negative values are caused by the second term of (3) involving the spatial derivative $\partial_1\chi$, and they are thus caused by a specific spatial dependence of the initial data.</footnote>
to $S'^{-1}_{\mu\nu}$ in the boosted frame of panel (d): When it points toward positive values of $x'$, it seems to go backwards with respect to time $t'$. As underlined above, the hypersurface $t' = 0$ cannot be consistently used to specify initial data in this case, since it is not spacelike with respect to $S'^{-1}_{\mu\nu}$, therefore the sign of Hamiltonian (3) at $t' = 0$ does not have much meaning anyway.] The conclusion is the same as before: The unbounded by below Hamiltonian in the boosted frame of panel (d) is a mere coordinate effect, without any physical meaning, and the model is actually stable, as proven by the positive total Hamiltonian density in the frame of panel (c).

It is also instructive to compute the energy of a system in a boosted frame (still in flat spacetime, to simplify the discussion). Although it differs from $g^{\mu\nu}$, the effective metric $S^{\mu\nu}$ is a tensor; see for instance Eq. (1) for the particular case of k-essence. Therefore, the Lagrangian $\mathcal{L}_2 = -\frac{1}{2}S^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi$ is diffeomorphism invariant, and this implies that four Noether currents are conserved. They read

$$-T_{\mu}^{\nu} = \frac{\delta \mathcal{L}_2}{\delta(\partial_{\nu}\chi)} \partial_{\mu}\chi - \delta_{\mu}^{\nu} \mathcal{L}_2, \quad (8)$$

where the index $\mu$ specifies which of the four currents is considered, $\nu$ denotes its components, and $\delta_{\mu}^{\nu}$ is the Kronecker symbol. [A global minus sign is introduced in definition (8) so that the mixed component $T_0^0$ denotes the opposite of the energy density, like in general relativity.] The current conservation reads as usual $\partial_{\nu}T_{\mu}^{\nu} = 0 \Leftrightarrow \partial_{0}T_{\mu}^{0} + \partial_{i}T_{\mu}^{i} = 0$. When integrating this identity over a large spatial volume $V$ containing the whole physical system under consideration, the spatial derivatives become vanishing boundary terms, and one gets the standard conservation laws for total energy and momentum, $\partial_{t}P_{\mu} = 0$, with $P_{\mu} = -\int_{V}T_{\mu}^{0}d^{3}x$. For $\mu = 0$, the energy density $-T_{0}^{0}$ coincides with the on-shell value of the Hamiltonian density (3). As recalled above, if it is bounded by below, then the lowest-energy state must be stable. But it should be underlined that the three components of the total momentum $P_{i}$ are also conserved, and that the components $-T_{i}^{0} = p\partial_{i}\chi$ [with $p$ still given by Eq. (2)] have no preferred sign, since there is no privileged spatial direction. When changing coordinates, the total 4-momentum of the system becomes $P'_{\lambda} = (\partial x^{\mu}/\partial x'^{\lambda})P_{\mu}$, and in particular, the energy gets mixed with the initial 3-momentum, $P'_{0} = (\partial x^{\mu}/\partial x'^{0})P_{\mu}$, or simply $P'_{0} = (P_{0} + vP_{1})/\sqrt{1-v^2}$ for a mere boost of velocity $v$ in the $x$ direction. Of course, $g^{\mu\nu}P_{\mu}P_{\nu}$ (as well as $S^{\mu\nu}P_{\mu}P_{\nu}$) is a scalar quantity, and it remains thus invariant under coordinate transformations. However, it is not always negative, contrary to the standard “minus rest mass squared” in special relativity, therefore the magnitude of the spatial components $P_{i}$ is not always bounded by $P_{0}$. For instance, in panels (c) or (d) of Fig. 1, a scalar field

perturbation propagating outside the solid (grey) cone obviously corresponds to a positive $g^{\mu\nu}P_{\mu}P_{\nu}$, i.e., a spacelike $P_{\mu}$ with respect to $g^{\mu\nu}$. It is thus clear that a negative value of $P'_{0} = (P_{0} + vP_{1})/\sqrt{1-v^2}$ is reachable for a large enough boost velocity $|v| < 1$. The fact that $P'_{0}$ can also become negative in the case of panel (a) is much less obvious, but it can be checked that it coincides with the spatial integral of the on-shell expression of Hamiltonian (3) with the boosted effective metric (7). In such a case, a large enough boost velocity $|c_s| < |v| < 1$ generates a negative $S'^{11}$, and thereby a possibly negative Hamiltonian (3), when initial data on the $t' = 0$ hypersurface are chosen with a large spatial gradient $\partial_{x'}\chi$ (but a small $\partial_{t'}\chi$). Up to now, we are merely rephrasing our previous conclusions with a slightly different viewpoint. But what is more interesting is to understand why situations like panels (a) or (d) of Fig. 1 are stable in spite of their Hamiltonian density (3) which is unbounded by below. The reason is simply that not only their total energy $P'_{0}$ is conserved, but also their 3-momentum $P'_{i}$. And it happens that the linear combination $(P'_{0} - vP'_{1})/\sqrt{1-v^2}$, which is thus also conserved, is bounded by below, since it obviously gives the positive expression of $P_{0}$ in the initial frame of panels (b) or (c). In other words, stability is not ensured by the boundedness by below of the Hamiltonian density, in the present case, but by that of the linear combination $-T'^{0}_{0} + vT'^{0}_{1}$. In more general situations involving arbitrary coordinate transformations, the initial energy $P_{0}$ which is bounded by below is again a linear combination of conserved quantities in the new frame, $P_{0} = (\partial x'^{\mu}/\partial x^{0})P'_{\mu}$.

In conclusion, although the Hamiltonian density is not bounded by below in the situations corresponding to panels (a), (d), (f), (g) and (h) of Fig. 1, there exists a choice of coordinates mapping them to panels (b), (c) or (e), where the new total Hamiltonian density is bounded by below. This suffices to guarantee the stability of the lowest-energy state, as computed in this new coordinate system. The only generically unstable cases correspond to the third row of Fig. 1, panels (i) to (l), because their total Hamiltonian density is never bounded by below in any coordinate system. They are such that the matrix $S^{\mu\lambda}g_{\lambda\nu}$ is diagonalizable and possesses two negative eigenvalues. Conversely, it is easy to write the inequalities needed on the components of the effective metric $S^{\mu\nu}$ to be in the eight safe cases corresponding to the first two rows, panels (a) to (h): In addition to the hyperbolicity condition (5), one just needs

$$S^{00} < S^{11} \quad \text{and/or} \quad |S^{00} + S^{11}| < 2|S^{01}|, \quad (9)$$

when focusing on the $(t, x)$ subspace in a coordinate system such that $g_{\mu\nu} = \text{diag}(-1, 1)$. But these inequalities are less enlightening than Fig. 1 itself, in which it is immediate to see whether the two causal cones have both a common exterior (when one should specify initial data) and a
common interior. When one chooses new coordinates such that time lies within the cone intersection, and space is outside both cones, then the total Hamiltonian density caused by kinetic terms becomes positive.

The above results can also be formulated in a covariant way. A given solution is stable if and only if all effective metrics (here $g^{\mu\nu}$ and $S^{\mu\nu}$, but also $G^{\mu\nu}$ introduced in Sec. III below) are of hyperbolic mostly plus signature, and it is possible to find a contravariant vector $U^\mu$ and a covariant vector $u_\mu$ such that

$$g_{\mu\nu} U^{\mu} U^{\nu} < 0, \qquad S_{\mu\nu}^{-1} U^{\mu} U^{\nu} < 0, \dots, \qquad (10)$$

$$g^{\mu\nu} u_{\mu} u_{\nu} < 0, \qquad S^{\mu\nu} u_{\mu} u_{\nu} < 0, \dots, \qquad (11)$$

and

$$T_{\mu}^{\nu} U^{\mu} u_{\nu} \geq 0, \qquad (12)$$

where $T_\mu^\nu$ denotes the Noether currents for all fields,<fnref n="10" /> including Eq. (8) for the scalar field. Equation (10) is the covariant way of formulating the existence of a common interior to all causal cones, where a “good” time direction may be chosen, namely $dx^0$ in the direction of $U^\mu$. Equation (11) expresses the existence of a spatial hyper-surface exterior to all causal cones, defined by $u_\mu dx^\mu = 0$, where “good” spatial coordinates may be chosen. Finally, Eq. (12) states that the Hamiltonian density is positive in such a “good” coordinate system. Let us underline that $U^\mu$ and $u_\mu$ are generically not related by lowering or raising the index with any of the effective metrics. This is the crucial difference with general relativity (with standard minimally coupled fields), in which a single metric $g_{\mu\nu}$ defines the causal cone of all d.o.f. In this simpler case of GR, the above conditions boil down to finding a single timelike vector $U^\mu$ for which the usual weak energy condition $T_{\mu\nu} U^\mu U^\nu \geq 0$ is satisfied. In particular, Eqs. (10) and (11) become then equivalent if one chooses $u_\mu = g_{\mu\nu} U^\nu$. One may also extend conditions (10)–(12) by imposing that all “future-oriented” contravariant and covariant vectors $U^\mu$ and $u_\mu$ satisfying (10) and (11) respect inequality (12). [By “future-oriented,” we mean here that these two vectors must have consistent orientations, i.e., that their scalar product $U^\mu u_\mu < 0$, otherwise one could change the sign of one of them without spoiling conditions (10) nor (11) but making (12) negative.] This would be the full generalization of the weak energy condition to our more subtle case involving several causal cones, and our previous discussion shows

that it would indeed be satisfied if the solution is stable. However, let us underline again that stability is actually ensured as soon as one pair of vectors $U^\mu$ and $u_\mu$ satisfies Eqs. (10)–(12).

As an application of the above results, let us rederive the stability conditions for the effective metric (1) corresponding to k-essence. Let us first choose a locally inertial frame such that $g_{\mu\nu} = \text{diag}(-1, 1, 1, 1)$. Then, if $\partial_\mu \bar{\phi}$ is timelike with respect to $g_{\mu\nu}$, it is always possible to boost this coordinate system such that $\partial_i \bar{\phi} = 0$. We thus get $S^{\mu\nu} = \text{diag}([-f' + 2\dot{\bar{\phi}}^2 f''], f', f', f')$. To be in the situation of panels (b) or (c) of Fig. 1, it is necessary to have $S^{00} < 0$ and $S^{xx} > 0$, therefore we need $-f' + 2\dot{\bar{\phi}}^2 f'' < 0$ and $f' > 0$. Since $\bar{X} = g^{\mu\nu} \partial_\mu \bar{\phi} \partial_\nu \bar{\phi} = -\dot{\bar{\phi}}^2$ in this specific coordinate system, the covariant expressions of these conditions are necessarily $f'(\bar{X}) > 0$ and $2\bar{X}f''(\bar{X}) + f'(\bar{X}) > 0$, as mentioned one paragraph below Eq. (1). Note that no condition is imposed on $f''(\bar{X})$ alone. The result remains the same when the background scalar gradient $\partial_\mu \bar{\phi}$ is spacelike (still with respect to $g_{\mu\nu}$). Then one may choose the x coordinate in its direction, so that its only non-vanishing component be $\bar{\phi}' \equiv \partial_1 \bar{\phi}$. In this coordinate system, the components of the effective metric read $S^{\mu\nu} = \text{diag}(-f', [f' + 2\dot{\bar{\phi}}^2 f''], f', f')$, while $\bar{X} = +\bar{\phi}'^2$, therefore we recover strictly the same covariant inequalities. Finally, when $\partial_\mu \bar{\phi}$ is a null vector (again with respect to $g_{\mu\nu}$, i.e., $\bar{X} = 0$), it is possible to choose a coordinate system in which $\partial_\mu \bar{\phi} = (\dot{\bar{\phi}}, \dot{\bar{\phi}}, 0, 0)$, and the nonvanishing components of the effective metric read $S^{00} = -f' + 2\dot{\bar{\phi}}^2 f''$, $S^{11} = f' + 2\dot{\bar{\phi}}^2 f''$, $S^{01} = S^{10} = -2\dot{\bar{\phi}}^2 f''$, and $S^{22} = S^{33} = f'$. We then find that one of the characteristics defined by $S^{\mu\nu}$ coincides with one of those defined by $g^{\mu\nu}$, corresponding to a velocity $-1$ for spin-0 perturbations. This is thus a limiting case of those plotted in Fig. 1. But when $f'(\bar{X}) > 0$, consistently with the same covariant inequalities as above, one finds that the causal cones defined by $g^{\mu\nu}$ and $S^{\mu\nu}$ have both a common interior and a common exterior, and the background solution is thus stable.

### III. STABLE BLACK HOLE SOLUTIONS IN A SUBCLASS OF (BEYOND) HORNDESKI THEORIES

Let us now illustrate our findings with a specific example, stemming from Horndeski theory. We will discuss certain solutions of the following action, which has been studied quite a lot due to its simple self-tuning properties:

$$S_J[g_{\mu\nu}, \varphi] = \int \sqrt{-g} d^4 x [\zeta(R - 2\Lambda_{\text{bare}}) + \beta G^{\mu\nu} \partial_\mu \varphi \partial_\nu \varphi \\ - \eta \varphi_\lambda^2], \qquad (13)$$

where we use the simplifying notation $\varphi_\lambda \equiv \partial_\lambda \varphi$, so that $\varphi_\lambda^2 = g^{\mu\nu} \partial_\mu \varphi \partial_\nu \varphi$ (which was also denoted as X in Sec. II).

<footnote n="10">¹⁰Of course, any other *conserved* tensor constructed from $T_\mu^\nu$ by adding the divergence of an antisymmetric Belinfante tensor is also allowed [70,71], and in particular the standard symmetric energy-momentum tensor $(2/\sqrt{-g})(\delta S_{\text{field}}/\delta g_{\mu\nu})$ defined as in general relativity (with one of its indices lowered with $g_{\mu\nu}$), where $S_{\text{field}}$ denotes the contribution of a given field to the action.</footnote>
$\zeta$ is the Planck mass squared divided by $16\pi$, and $\eta, \beta$ and $\Lambda_{\text{bare}}$ are some constants. In terms of standard Horndeski notation, this action corresponds to

$$G_4 = \zeta - \frac{\beta}{2} \varphi_{\lambda}^2, \quad G_2 = -2\zeta\Lambda_{\text{bare}} - \eta\varphi_{\lambda}^2. \quad (14)$$

Static and spherically symmetric black hole solutions of the above theory were first derived in [72] while they were extended in [54,73,74] to the case of nonvanishing $\Lambda_{\text{bare}}$. A new family of solutions, with a linearly time-dependent scalar field, was proposed in [54]. These solutions enjoy novel regularity properties thanks to the time dependence of the scalar. Some solutions have spacetime metrics that are identical to their GR counterparts (apart from the value of the cosmological constant). As a result, they are often referred to as stealth solutions. More importantly, time dependence of the scalar field qualifies the scalar to be a dark energy field responsible for late-time acceleration (as well as self-tuning properties). These solutions were claimed to be unstable under linear perturbations [55], and more recently the theory (13) was ruled out observationally. The aim of the forthcoming section is to show that the former result is in fact wrong, while the latter crucially depends on how the metric couples to matter. Put in other words, if the physical metric to which matter couples minimally is $g_{\mu\nu}$, then the above theory is ruled out (more precisely, the scalar field is ruled out as a dark energy candidate). Indeed, the speed of gravitons in this theory generically deviates from the speed of light [46,47] in inconsistency with the simultaneous observation of gravitational and electromagnetic waves from the same source, GW170817 [2,45]. However, it is easy to map the action (13) to a beyond Horndeski theory in which gravitational waves do travel at the speed of light in accordance with observations. This has been checked in weakly curved backgrounds [46,47] but also in strongly curved spherically symmetric backgrounds [53] (and Ref. [20] recently proved so for vector-tensor theories too). To make therefore the theory (13) viable, the matter action should be minimally coupled to $\tilde{g}_{\mu\nu}$, the physical metric:

$$\tilde{g}_{\mu\nu} = g_{\mu\nu} - \frac{\beta}{\zeta + \frac{\beta}{2}\varphi_{\lambda}^2} \partial_{\mu}\varphi\partial_{\nu}\varphi. \quad (15)$$

Of course, any metric proportional to this $\tilde{g}_{\mu\nu}$ would also be allowed, since it would not change the causal cone, even if the conformal factor depends on $\varphi_{\lambda}^2$. One should then work with the action

$$S_j[g_{\mu\nu}, \varphi] + S_m[\tilde{g}_{\mu\nu}, \Psi], \quad (16)$$

where $S_m$ is some given matter action with matter fields, collectively denoted as $\Psi$, universally coupled to the physical metric $\tilde{g}_{\mu\nu}$. In standard nomenclature for BD

gravity, the nonphysical $g_{\mu\nu}$ would be called the “Einstein frame” metric. However, its perturbations do not describe pure spin-2 d.o.f. in the present case, because of the kinetic mixing introduced by the $G^{\mu\nu}\varphi_{\mu}\varphi_{\nu}$ term of action (13). We will therefore call $g_{\mu\nu}$ the “Horndeski frame” metric rather than the “Einstein frame” one. On the other hand, we call $\tilde{g}_{\mu\nu}$ the “Jordan frame” physical metric. As in standard BD theory, it is easier to work in the nonphysical frame because the metric sector is simpler there. We should keep in mind that our analogy is to be taken with caution, because the frames of the higher order theories are related disformally (15), and not conformally as in BD theory. Indeed, the disformal factor (15) has been chosen in order to impose a unit speed for the gravitational waves in the physical (or Jordan) frame, at least in weakly curved backgrounds. We recently reported [53] that the black hole solutions found in [54] are again black holes with respect to the physical frame. This is not a trivial result, as a disformal transformation may change the nature of solutions, rendering them even singular upon going from one frame to the other. We will explicitly work out the physical disformed metric in the next section. Furthermore, we study the stability of some solutions of the theory (16). We do so in the Horndeski frame, as stability properties carry through upon field redefinitions (15) as long as these are not singular. A priori, three causal cones must be considered in our analysis, and must have compatible orientations for the solutions to be stable: the matter causal cone associated to $\tilde{g}_{\mu\nu}$, and the cones associated to scalar and gravitational perturbations (with their associated effective metrics). Quite remarkably, as we will see, the graviton perturbation cone will end up being identical to the matter light cone in the physical frame, demonstrating that gravity waves travel at same speed as light, even in a strongly curved region of spacetime (close to the event horizon). This will effectively reduce the number of causal cones under scrutiny from three to two. We will see in the next section how to construct these causal cones and effective metrics in a spherically symmetric background.

Regarding stability, Appleby and Linder examined action (13) with vanishing $\Lambda_{\text{bare}}$ in a cosmological framework [75]. From the study of scalar perturbations, they found that there always exists either a gradient instability or a ghost. This pathology can however be cured by the introduction of a bare cosmological constant $\Lambda_{\text{bare}}$, as we will see below. The stability of the black hole static solutions was discussed in [76–78], based on a more generic theory than (13). The authors employed the well-established Regge-Wheeler formalism: they decomposed the perturbations into odd and even modes, each mode being decoupled of all others at the linear level. Stable parameter regions were exhibited for the action (13). Then, Ogawa et al. tackled the case where the scalar field acquires time-dependence [55]. They claimed that the solutions were always unstable, whatever the coupling
parameters of the theory. However, their argument made use of the fact that the Hamiltonian is unbounded from below; as we argued in the former section, this cannot be a satisfactory criterion to decide on the stability of some solution. We show in the last paragraph of this section that there indeed exist stable black hole solutions for given parameters. We will first derive the effective metrics in which graviton and scalar perturbations respectively propagate.

### A. The effective metrics for graviton and scalar perturbations

We will focus our analysis on perturbation theory around spherically symmetric Schwarzschild-de Sitter solutions. It is indeed known that the action (13) allows for a “stealth” Schwarzschild black hole, as well as Schwarzschild-de Sitter metrics with a non-trivial scalar profile:

$$ds^2 = -A(r)dt^2 + \frac{dr^2}{B(r)} + r^2(d\theta^2 + \sin^2\theta d\phi^2), \quad (17)$$

$$A(r) = B(r) = 1 - \frac{2Gm}{r} - \frac{\Lambda_{\text{eff}}}{3} r^2, \quad (18)$$

$$\Lambda_{\text{eff}} = -\frac{\eta}{\beta}, \quad (19)$$

$$\varphi = q \left[ t \pm \int \frac{\sqrt{1 - A(r)}}{A(r)} dr \right], \quad (20)$$

$$q^2 = \frac{\eta + \beta \Lambda_{\text{bare}}}{\eta \beta} \zeta, \quad (21)$$

where $q$ parametrizes the linear time-dependence of the scalar field, and $m$ corresponds to the mass of the black hole. The constant $\Lambda_{\text{eff}}$ plays the role of an effective cosmological constant, and is a priori independent of $\Lambda_{\text{bare}}$ in Eq. (13), with the velocity integration constant $q$ playing the role of a tuning integration constant to $\Lambda_{\text{bare}}$ via relation (21). For consistency, the right-hand side of Eq. (21) should be positive; since $\zeta$ is always positive, we must therefore have

$$(\eta + \beta \Lambda_{\text{bare}}) \eta \beta > 0, \quad (22)$$

for this solution. A second background solution obtained in the case $\eta = 0$ and $\Lambda_{\text{bare}} = 0$ is the stealth Schwarzschild black hole solution, which reads

$$ds^2 = -A(r)dt^2 + \frac{dr^2}{B(r)} + r^2(d\theta^2 + \sin^2\theta d\phi^2), \quad (23)$$

$$A(r) = B(r) = 1 - \frac{2Gm}{r} \quad (24)$$

$$\varphi = q \left[ t \pm \int \frac{\sqrt{1 - A(r)}}{A(r)} dr \right], \quad (25)$$

and for which $q$ is a free parameter. This solution is characterized by an asymptotically flat metric, does not have self-tuning properties and is not a limit of the de Sitter black hole (17).

Before we proceed to the stability analysis let us apply the disformal transformation (15) to the above background solutions and examine the nature of the physical metric and scalar field background solution (17). The family of background solutions verifies the relation $\varphi_{\lambda}^2 = -q^2$ and this simplifies the disformal transformation and coordinate transformations thereof. To simplify notation, we can set

$$\mathcal{F} = - \frac{\beta q^2}{\zeta - \frac{\beta}{2} q^2}, \quad (26)$$

and note that the disformed metric acquires off diagonal terms in the original $(t, r)$ coordinates due to the $t$ and $r$ scalar field dependence. We then diagonalize the physical metric using

$$\tilde{t} = \sqrt{1 - \mathcal{F}} \left\{ t \mp \int \frac{\mathcal{F} \sqrt{1 - A(r)}}{A(r)[A(r) - \mathcal{F}]} dr \right\}, \quad (27)$$

where the minus sign corresponds to the plus one in Eqs. (20) and (25). Note that, for this coordinate transformation to be well defined, one needs

$$\mathcal{F} < 1. \quad (28)$$

For the solution (17)–(21), this bound reads

$$(3\eta + \beta\Lambda_{\text{bare}})(\eta - \beta\Lambda_{\text{bare}}) > 0. \quad (29)$$

We have to keep in mind this constraint for the upcoming stability analysis. The background solution (20) in the physical frame $\tilde{g}_{\mu\nu}$ then recovers the same form as the original background, namely:

$$d\tilde{s}^2 = -\tilde{A}(r)d\tilde{t}^2 + \frac{dr^2}{\tilde{B}(r)} + r^2(d\theta^2 + \sin^2\theta d\phi^2), \quad (30)$$

$$\tilde{A}(r) = \tilde{B}(r) = 1 - \frac{2G\tilde{m}}{r} - \frac{\tilde{\Lambda}_{\text{eff}}}{3} r^2, \quad (31)$$

$$\varphi = \tilde{q} \left[ \tilde{t} - \int \frac{\sqrt{1 - \tilde{A}(r)}}{\tilde{A}(r)} dr \right], \quad (32)$$

$$\tilde{q} = \frac{q}{\sqrt{1 - \mathcal{F}}}, \quad \tilde{m} = \frac{m}{1 - \mathcal{F}},$$

$$\tilde{\Lambda}_{\text{eff}} = \frac{\Lambda_{\text{eff}}}{1 - \mathcal{F}} = \left( \frac{\Lambda_{\text{eff}} + \Lambda_{\text{bare}}}{3\Lambda_{\text{eff}} - \Lambda_{\text{bare}}} \right) \Lambda_{\text{eff}}, \quad (33)$$
in the $(\tilde{t}, r)$ coordinate system with respect to the rescaled parameters of the solution<fnref n="11" /> $\tilde{\Lambda}_{\text{eff}}$, $\tilde{m}$, $\tilde{q}$. Before studying the stability, let us remark that the solution in the physical frame is asymptotically de Sitter only for positive $\tilde{\Lambda}_{\text{eff}}$. Since this solution is meant to describe our current Universe, we impose the positivity of $\tilde{\Lambda}_{\text{eff}}$. In terms of the Lagrangian parameters, this translates as:

$$ \eta\beta(\eta - \beta\Lambda_{\text{bare}})(3\eta + \beta\Lambda_{\text{bare}}) < 0, \quad (34) $$

to be combined with constraints (22) and (29). It is easy to check that the three conditions together imply that the solution was also asymptotically de Sitter in the original Horndeski frame, i.e., that $\Lambda_{\text{eff}} > 0$. The above transformation can be trivially extended to the stealth solution (23)–(25). We hence recover the announced result that the physical metrics are again black hole solutions.

$$ h_{\mu\nu} = \begin{pmatrix} A(r)H_0(t, r) & H_1(t, r) & 0 & 0 \\ H_1(t, r) & H_2(t, r)/B(r) & 0 & 0 \\ 0 & 0 & K(t, r)r^2 & 0 \\ 0 & 0 & 0 & K(t, r)r^2 \sin^2\theta \end{pmatrix}, \quad (37) $$

where the $H_i$ and $K$ are free functions. Inserting these perturbations into the action, we isolate the terms which are quadratic in $h_{\mu\nu}$ and $\chi$. This gives the second order perturbed action, that we can write as

$$ \delta_s^{(2)} S_J = \int dt dr 4\pi r^2 \mathcal{L}_s^{(2)}, \quad (38) $$

where the factor $4\pi r^2$ corresponds to the trivial angular integration, and $\mathcal{L}_s^{(2)}$ is the Lagrangian density from which we can extract the causal structure of the perturbations. The subscript “s” stands for scalar, since we choose to excite only a spherically symmetric mode. We can simplify the calculations using the diffeomorphism invariance generated by an infinitesimal vector $\xi^\mu$. In the new system of coordinates $\hat{x}^\mu = x^\mu + \xi^\mu$, the metric and scalar transform according to

$$ \hat{g}_{\mu\nu} = g_{\mu\nu} - 2\nabla_{(\mu}\xi_{\nu)}, \quad (39) $$

Let us now proceed with the perturbative analysis. Generically, the theory (13) has one scalar d.o.f., and two polarizations of a massless spin-2 d.o.f. We first want to obtain the effective metric in which the scalar mode propagates. To this end, we shall focus on a spherically symmetric perturbation. If such a dynamical breather mode exists, it necessarily corresponds to the scalar d.o.f. We perturb the metric and scalar field according to

$$ g_{\mu\nu} = \bar{g}_{\mu\nu} + h_{\mu\nu}, \quad (35) $$

$$ \varphi = \bar{\varphi} + \chi, \quad (36) $$

where a bar denotes the background solution, and $h_{\mu\nu}$ and $\chi$ depend only on $t$ and $r$ since we look for a spherically symmetric perturbation. Using the formalism developed by Regge and Wheeler [80], $h_{\mu\nu}$ can be written in spherical coordinates as:

$$ \hat{\varphi} = \varphi - \partial_\mu \varphi \xi^\mu. \quad (40) $$

With a well-chosen $\xi^\mu$, we can in fact set $K$ and $\chi$ to zero. This completely fixes the gauge. Explicitly,

$$ \xi^\mu = \left( \frac{1}{q} \left( \chi + \varphi' \frac{Kr}{2} \right), -\frac{Kr}{2}, 0, 0 \right), \quad (41) $$

a prime standing for a derivative with respect to $r$. In this gauge, $\mathcal{L}_s^{(2)}$ reads, after numerous integrations by parts and using the background field equations,

$$ \begin{aligned} \mathcal{L}_s^{(2)} = & c_1 H_0 \dot{H}_2 + c_2 H'_0 H_1 + c_3 H'_0 H_2 + c_4 H_1 \dot{H}_2 \\ & + c_5 H_0^2 + c_6 H_2^2 + c_7 H_0 H_2 + c_8 H_1 H_2. \end{aligned} \quad (42) $$

Here a dot represents a time derivative, and all $c_i$ are background coefficients with radial (but no time) dependence, the detailed expression of which can be found in Appendix. This three-field Lagrangian should boil down to a Lagrangian depending on a single dynamical variable. As a first step in this direction, it is easy to eliminate $H_2$ since the associated field equation is algebraic in $H_2$:

$$ H_2 = -\frac{1}{2c_6} (-c_1 \dot{H}_0 - c_4 \dot{H}_1 + c_3 H'_0 + c_7 H_0 + c_8 H_1). \quad (43) $$

Inserting back this expression in $\mathcal{L}_s^{(2)}$, we obtain,

<footnote n="11">¹¹Similarly to Eq. (20), there actually exist two branches for the scalar field, corresponding to a plus or minus sign in front of the $r$ integral. In the physical frame, we keep only the minus branch, so that this solution is mapped to a homogeneous and expanding one in FLRW coordinates, see [79]. This minus sign actually also corresponds to a minus sign in the Horndeski frame of Eq. (20).</footnote>
$$ \mathcal{L}_s^{(2)} = \tilde{c}_1 \dot{H}_0^2 + \tilde{c}_2 H_0'^2 + \tilde{c}_3 H_0' \dot{H}_0 + \tilde{c}_4 \dot{H}_1^2 + \tilde{c}_5 \dot{H}_0 H_1' \\ + \tilde{c}_6 \dot{H}_0 \dot{H}_1 + \tilde{c}_7 H_0' H_1 + \tilde{c}_8 \dot{H}_0 H_1 \\ + \tilde{c}_9 H_0^2 + \tilde{c}_{10} H_0 H_1 + \tilde{c}_{11} H_1^2, \qquad (44) $$

where the $\tilde{c}_i$ coefficients are again given in Appendix in terms of the $c_i$. A trickier step is to trade $H_1$ and $H_0$ for a single variable, since the associated field equations are differential equations, not algebraic ones. To this end, we introduce an auxiliary field $\pi_s$ as a linear combination of $H_0, H_1$ and their first derivatives:

$$ \pi_s = \dot{H}_0 + a_2 H_0' + a_3 \dot{H}_1 + a_4 H_1' + a_5 H_0 + a_6 H_1, \quad (45) $$

with some $a_i$ coefficients to be determined soon. The idea is to introduce $\pi_s$ at the level of the action, group all the derivatives inside $\pi_s$, and then to solve for the algebraic equations giving $H_0$ and $H_1$ in terms of $\pi_s$. Therefore, we rewrite the Lagrangian as

$$ \mathcal{L}_s^{(2)} = a_1[-\pi_s^2 + 2\pi_s(\dot{H}_0 + a_2H_0' + a_3\dot{H}_1 + a_4H_1' + a_5H_0 + a_6H_1)] \\ + a_7H_0^2 + a_8H_1^2 + a_9H_0H_1. \qquad (46) $$

Variation of (46) with respect to $\pi_s$ ensures Eq. (45). Now, a simple identification with Lagrangian (44) allows us to determine the $a_i$ in terms of the $\tilde{c}_i$. Again, these coefficients are given in Appendix. Variation of (46) with respect to $H_0$ and $H_1$ gives a system of two linear equations, which we can easily solve to write these two fields in terms of $\pi_s$ and its derivatives. We do not write down their expression here because of their consequent length, but the procedure is straightforward.<fnref n="12" /> At this point, we have obtained a Lagrangian density in terms of a single variable $\pi_s$. We will examine its kinetic part only, neglecting the potential associated to this d.o.f. and thereby focusing on the causal structure. This kinetic part reads

$$ \mathcal{L}_{s; \text{Kin}}^{(2)} = -\frac{1}{2}(S^{tt}\dot{\pi}_s^2 + 2S^{tr}\dot{\pi}_s\pi_s' + S^{rr}\pi_s'^2), \quad (47) $$

with

$$ S^{tt} = \frac{c_1^2 c_3^2 c_4^2}{4c_2 D} (-2c_4^2 c_5 + c_2 c_1 c_4' - c_2 c_4 c_1' - c_1 c_4 c_2'), \quad (48) $$

$$ S^{rr} = -\frac{c_1^2 c_3^2 c_4^2}{2D} (-c_3 c_8 + c_2 c_6), \quad (49) $$

$$ S^{rr} = -\frac{c_1^2 c_3^2 c_4^2}{4D} (-c_4 c_7 + c_1 c_8), \qquad (50) $$

$$ \mathcal{D} = c_6^2 \{2(-c_3 c_8 + c_2 c_6)(-c_2 c_4 c_1' + c_1 c_2 c_4' - c_1 c_4 c_2') \\ + [4c_3 c_8 c_5 + c_2(c_7^2 - 4c_6 c_5)]c_4^2 \\ - 2c_2 c_4 c_7 c_1 c_8 + c_2 c_1^2 c_8^2\}. \qquad (51) $$

Alternatively, we can remark that the scalar mode propagates to linear order in the given black hole background (17) with an effective two-dimensional metric $S_{\mu\nu}^{-1}$:

$$ \mathcal{L}_{s; \text{Kin}}^{(2)} = -\frac{1}{2} S^{\mu\nu} \partial_{\mu} \pi_s \partial_{\nu} \pi_s. \qquad (52) $$

We can read from Eq. (47) the inverse metric:

$$ S^{\mu\nu} = \begin{pmatrix} S^{tt} & S^{tr} \\ S^{tr} & S^{rr} \end{pmatrix}, \qquad (53) $$

and the metric itself:

$$ S_{\mu\nu}^{-1} = \frac{1}{S^{tt}S^{rr} - (S^{tr})^2} \begin{pmatrix} S^{rr} & -S^{tr} \\ -S^{tr} & S^{tt} \end{pmatrix}. \qquad (54) $$

From this last object, we can determine the hyperbolicity condition, the propagation speeds, and all the information we need for the causal structure of the scalar mode. The hyperbolicity condition for instance reads

$$ (S^{tr})^2 - S^{tt}S^{rr} > 0. \qquad (55) $$

The speed of a wave moving towards or away from the origin is then given by

$$ c_s^{\pm} = \frac{S^{tr} \pm \sqrt{(S^{tr})^2 - S^{tt}S^{rr}}}{S^{tt}}. \qquad (56) $$

The hyperbolicity condition ensures that these propagation speeds are well defined. At any given point, $c_s^+$ and $c_s^-$ generate the scalar causal cone. Finally, one needs to know where the interior of the cone is located. This can be easily determined by checking whether a given direction (for instance the one generated by the vector $\partial_t$) is time or spacelike with respect to the metric $S_{\mu\nu}^{-1}$.

A similar analysis must be carried out for the spin-2 mode. It was actually already realized by Ogawa et al. in [55]. They studied odd-parity perturbations, which cannot correspond to a scalar d.o.f.—the latter always has even parity. Hence odd-parity perturbations correspond to one of the two spin-2 polarizations. We checked the calculations of Ogawa et al., and we are in full agreement as to the quadratic Lagrangian derived in their paper. For brevity, we only reproduce the final result here, applied to the solution (17)–(21); the gravity perturbations propagate in a

<footnote n="12">¹²The case of the stealth Schwarzschild black hole [54] is more subtle. There, the two equations determining $H_0$ and $H_1$ become linearly dependent and the procedure cannot be applied. Section III D is devoted to this particular case.</footnote>
two-dimensional effective metric $G_{\mu\nu}^{-1}$, which essentially coincides with the physical metric $\tilde{g}_{\mu\nu}$, as given by Eq. (15):

$$G_{\mu\nu}^{-1} = \frac{\Lambda_{\text{eff}}}{\Lambda_{\text{bare}} + \Lambda_{\text{eff}}} \tilde{g}_{\mu\nu}. \qquad (57)$$

The two metrics are related by a constant conformal factor. Therefore, they have identical causal structure at any point of spacetime, provided that the conformal factor is positive. When this is the case, matter and gravitons propagate exactly the same way; it is enough to analyze the overlap conditions of two of the three causal cones, say the matter and scalar one. On the contrary, if the above conformal factor is negative, the cone of the graviton is exactly complementary to the matter one and they have no overlap nor common exterior. We are therefore led to impose that

$$\Lambda_{\text{eff}}(\Lambda_{\text{bare}} + \Lambda_{\text{eff}}) > 0, \qquad (58)$$

i.e., in terms of the Lagrangian parameters:

$$\eta(\eta - \beta\Lambda_{\text{bare}}) > 0, \qquad (59)$$

The hyperbolicity condition coming from $G_{\mu\nu}^{-1}$, or equivalently $\tilde{g}_{\mu\nu}$, reads:

$$(G^{ir})^2 - G^{tt}G^{rr} > 0, \qquad (60)$$

and the speeds of inwards/outwards moving gravitons are given by

$$c_{\text{g}}^{\pm} = \frac{G^{ir} \pm \sqrt{(G^{ir})^2 - G^{tt}G^{rr}}}{G^{tt}}. \qquad (61)$$

In a nutshell, we have found the effective metrics in which gravitons and massless scalar propagate and they are given by $G^{\mu\nu}$ and $S^{\mu\nu}$ respectively. The Schwarzschild-de Sitter solution is well defined in both Einstein and physical frames provided Eqs. (22), (29) and (34). Additionally, under condition (59), gravitons propagate in the exact same metric as matter. In particular, the speed of gravitational waves is identical to the speed of light in a strongly curved background.

## B. Homogeneous solutions: A stability window

We will first apply the above analysis to de Sitter solutions, that is solution (17)–(21) with $m = 0$. Of course, in this case, the analysis presented above is not strictly necessary, but it allows us to cross check our results with cosmological perturbation theory. In particular, we arrive at the same conclusion as [75] for the model (13) with vanishing $\Lambda_{\text{bare}}$: There is no stable homogeneous configuration. However, switching on a nontrivial $\Lambda_{\text{bare}}$, the hyperbolicity conditions (55) and (60) read respectively:

$$(3\beta\Lambda_{\text{bare}} + \eta)(\eta - \beta\Lambda_{\text{bare}}) < 0, \qquad (62)$$

$$(3\eta + \beta\Lambda_{\text{bare}})(\eta - \beta\Lambda_{\text{bare}}) > 0. \qquad (63)$$

These two conditions must be supplemented with the fact that the graviton and scalar cones have a nonempty intersection and a common exterior. Again, compatibility with the matter causal cone will follow automatically, since $\tilde{g}^{\mu\nu}$ and $G^{\mu\nu}$ are conformally related, with a positive factor provided Eq. (59). It is enough to check the orientation of the cones at $r = 0$, since the solution under analysis is homogeneous. If the cones have compatible orientations at $r = 0$, this will remain true everywhere else. The calculation is then particularly simple, since $S^{ir}$ and $G^{ir}$ vanish at $r = 0$, meaning that the cones are either aligned (and symmetric around the $t$ axis) or inclined at ninety degrees. As soon as the hyperbolicity condition (62) for the scalar and the constraint (59) are satisfied, the $t$ axis is contained in the cone associated to $S_{\mu\nu}^{-1}$. We therefore need the graviton cone to contain also the $t$ axis. This is the case if

$$\eta(3\eta + \beta\Lambda_{\text{bare}}) > 0. \qquad (64)$$

Thus, there are in total seven conditions to fulfill for stability and existence of the solution: Eqs. (22), (29), (34), (59), (62), (63) and (64). They actually define an nonempty subspace of the parameter space. The cosmological solution is stable if and only if

$$\text{either } \eta > 0, \beta < 0 \text{ and } \frac{\Lambda_{\text{bare}}}{3} < -\frac{\eta}{\beta} < \Lambda_{\text{bare}}, \qquad (65)$$

$$\text{or } \eta < 0, \beta > 0 \text{ and } \Lambda_{\text{bare}} < -\frac{\eta}{\beta} < 3\Lambda_{\text{bare}}. \qquad (66)$$

In the following section, we give an example of parameters that fulfill this criterion. Let us stress that the above restrictions prevent one from using the theory (13) as a self-tuning model. Indeed, the above equations tell us that the effective cosmological constant has to be of same magnitude as the bare one. Rewriting these conditions in terms of the observed $\tilde{\Lambda}_{\text{eff}}$, we obtain

$$\text{either } \eta > 0, \beta < 0 \text{ and } \Lambda_{\text{bare}} < \tilde{\Lambda}_{\text{eff}}, \qquad (67)$$

$$\text{or } \eta < 0, \beta > 0 \text{ and } \Lambda_{\text{bare}} < \tilde{\Lambda}_{\text{eff}} < \frac{3}{2}\Lambda_{\text{bare}}. \qquad (68)$$

Again, this means that self-tuning is impossible in this specific model, since the observed cosmological constant must always be larger than the bare one.
### C. Black holes in de Sitter: Example of a stable configuration

Our analysis is fully relevant when the solution no longer describes a homogeneous cosmology, but rather a black hole embedded in such a cosmology. The three conditions for the background solution to exist and the frame transformation to be well defined, Eqs. (22), (29), (34), do not depend on the presence of a mass $m \neq 0$. Therefore, they remain identical when a black hole is present. Additionally, the condition (59) for $\tilde{g}_{\mu\nu}$ and $G_{\mu\nu}^{-1}$ to have compatible orientations is unchanged. Provided Eq. (59), photons coupled to $\tilde{g}_{\mu\nu}$ and gravitons will travel with the exact same speed, even in a highly curved background. This is a very positive feature of the theory (13). Indeed, comparing the speed of gravitational and electromagnetic waves with arbitrary accuracy cannot rule out this model.

On the other hand, the expressions of $G_{\mu\nu}^{-1}$ and $S_{\mu\nu}^{-1}$ become very complicated with a nonvanishing black hole mass. It is still possible to prove that the hyperbolicity conditions for both $S_{\mu\nu}^{-1}$ and $G_{\mu\nu}^{-1}$ are not modified with respect to the de Sitter case. They are again given by Eqs. (62), (63). To ensure the compatibility of orientation between the scalar cone and the graviton one is however more tricky. We checked numerically that the condition (64) for these two cones to be compatible in the de Sitter case leads to compatible cones also when the mass parameter $m$ is switched on. That is, for parameters in the range (65) or (66), the scalar cone seems to have a compatible orientation with the graviton cone even close to the black hole horizon. This remains true for arbitrary mass of the black hole (as long as the black hole horizon remains smaller than the cosmological horizon). Figure 2 provides an illustrative example of this numerical check, for a given set of parameters that falls in the range allowed by the corresponding de Sitter solution. In this case, the cones have compatible orientations everywhere. Remarkably, the scalar causal cone entirely opens up when approaching the black hole horizon, without becoming pathological.

Let us stress here why Ref. [55] would have claimed that the situation exhibited in Fig. 2 is unstable, in the light of the discussion of Sec. II. In this paper, the graviton metric components $G^{tt}$ and $G^{rr}$ were required to be negative and positive respectively. It was proven, however, that the product $G^{tt}G^{rr}$ is always positive in the vicinity of a horizon. Figure 2 shows that, indeed, the $t$ axis “leaves” the causal cone of the graviton (red cone), close to the event and cosmological horizon, while the $r$ axis remains in the exterior of the cone. This makes the quantity $G^{tt}G^{rr}$ positive close to the horizon, and the associated Hamiltonian unbounded by below. However, our analysis so far clearly shows that it does not signal an instability in any way.

From Fig. 2, one can additionally draw conclusions on the case where matter couples to $g_{\mu\nu}$, rather than $\tilde{g}_{\mu\nu}$. Of course, matter was chosen to couple to $\tilde{g}_{\mu\nu}$ on physical

grounds. However, in a generic scalar-tensor theory where no relation is assumed between the propagation of gravity and light, one has to investigate these three different cones. The time coordinate $t'$ in Fig. 2 is rescaled with respect to $t$ in such a way that the causal cone associated with $g_{\mu\nu}$ is at $\pm 45^\circ$, with a timelike $t'$ direction (that is, the $t'$ axis lies inside the cone of $g_{\mu\nu}$). Even in this more restrictive situation, the plots of Fig. 2 show that the three cones would actually be compatible, and the solution would be stable.

### D. A special case: Stealth Schwarzschild black hole

As mentioned above, there exists an exact asymptotically flat Schwarzschild solution when $\eta$ and $\Lambda_{\text{bare}}$ vanish, with a nontrivial scalar profile (23)–(25). In this case, the parameter $q$ is no longer related to the coupling constants of the action and is in fact a free parameter. We should also emphasize that solution (23)–(25) is the unique static and spherically symmetric solution with a linearly time dependent scalar field and $\eta = \Lambda_{\text{bare}} = 0$. The procedure for determining the effective metric of scalar perturbations, described in Sec. III A, breaks down for this background. It is not possible to carry on after Eq. (44), and to express the fields $H_0$ and $H_1$ in terms of $\pi_s$. The reason is that for the stealth Schwarzschild solution $H_0$ and $H_1$ cannot be simultaneously expressed in terms of the master variable $\pi_s$ from the Lagrangian introduced in (46). Therefore, we need to find another way of extracting the scalar mode from the second-order Lagrangian (42) which now reads

$$ \mathcal{L}_s^{(2)} = c_1 H_0 \dot{H}_2 + c_2 H'_0 H_1 + c_3 H'_0 H_2 + c_4 H_1 \dot{H}_2 + c_6 H_2^2 + c_7 H_0 H_2 + c_8 H_1 H_2, \quad (69) $$

as $c_5 = 0$ for the relevant background. The equation of motion for $H_2$ following from (42) is algebraic in terms of $H_2$, so as before we can find $H_2$ in terms of $H_0$ and $H_1$:

$$ H_2 = -\frac{1}{2c_6} (-c_1 \dot{H}_0 - c_4 \dot{H}_1 + c_3 H'_0 + c_7 H_0 + c_8 H_1). \quad (70) $$

Substituting (70) in (69) and rearranging terms, we can write (69) as

$$ \mathcal{L}_s^{(2)} = a_1 (\dot{H}_0 + a_2 H'_0 + a_3 \dot{H}_1 + a_5 H_0 + a_6 H_1)^2 + a_7 H_0^2 + a_8 H_1^2 + a_9 H_0 H_1, \quad (71) $$

where the coefficients $a_i$ are given in the Appendix. We now introduce new variables $x(t, r)$ and $y(t, r)$ grouping together the time and space derivatives:
<figure><img src="image_4.png" /><figcaption>FIG. 2. The scalar and graviton/matter causal cones in Schwarzschild-de Sitter geometry, respectively in dashed blue and plain red. The parameters of the Lagrangian are chosen so that the associated cosmological solution is stable: $\eta = \frac{1}{2}$, $\beta = -1$, $\zeta = 1$, $\Lambda_{\text{bare}} = 1$ in Planck units. The radius $r$ varies between the black hole horizon located at $r \approx 9.4 \times 10^{-3}$ and the cosmological horizon at $r \approx 2.4$. For this set of parameters, the graviton cone always lies inside the scalar cone; as a consequence, they have compatible orientations. In this plot, the time coordinate $t'$ has been rescaled with respect to the original one, so that the causal cone associated to the unphysical metric $g_{\mu\nu}$ corresponds to lines at $\pm 45^\circ$.</figcaption></figure>

$$H_1 = \frac{c_1}{c_4} \left( x - \frac{c_1}{c_3} y \right), \qquad H_0 = \frac{c_1}{c_3} y. \qquad (72)$$

Indeed, the Lagrangian (71) then takes the form

$$\mathcal{L}_s^{(2)} = \mathcal{P}^2 + \mathcal{A}y^2 + Bxy + Cx^2, \qquad (73)$$

where

$$\mathcal{P} = \dot{x} - y' + \tilde{a}_1 x + \tilde{a}_2 y, \qquad (74)$$

$$\tilde{a}_1 = \frac{2c_2c_6 - c_3c_8}{c_3c_4},$$

$$\tilde{a}_2 = \frac{c_4c_1c'_3 - c_3c_4c'_1 + c_8c_1^2 - c_4c_7c_1}{c_1c_3c_4},$$

$$\mathcal{A} = \frac{c_1^2[c_4(c_2c'_1 + c_1c'_2 + 2c_4c_5) - c_1c_2c'_4]}{2c_3^2c_4^2},$$

$$\mathcal{B} = \frac{c_1^2c_2(c_1c_8 - c_4c_7)}{c_3^2c_4^2},$$

$$\mathcal{C} = \frac{c_1^2c_2(c_2c_6 - c_3c_8)}{c_3^2c_4^2}. \qquad (75)$$

and

Variation of (73) with respect to $y$ yields the constraint
$$2\mathcal{P}' + 2A\mathcal{y} + B\mathcal{x} = 0. \qquad (76)$$

The above constraint (76) contains $y''$, $y'$, $\dot{x}'$ and it may be seen as an equation which determines $\mathcal{y}$ in terms of $\mathcal{x}$ and its derivatives. To find $\mathcal{y}$ from (76), the use of nonlocal (in space) operators is in general required. For our purposes, however, we do not need to know the exact expression of $\mathcal{y}$ in terms of $\mathcal{x}$, since we are only interested in the absence of ghost and gradient instabilities. This means that we focus on higher derivative terms, i.e., we neglect $\sim \mathcal{x}$ with respect to $\sim \dot{\mathcal{x}}$ or $\sim \mathcal{x}'$ as well as $\sim \mathcal{y}$ with respect to $\mathcal{y}'$. With this approximation in mind, Eq. (76) becomes

$$\dot{x}' - \mathcal{y}'' = 0, \qquad (77)$$

which after integration over $r$ and setting to zero the integration constant yields

$$\dot{x} = \mathcal{y}'. \qquad (78)$$

By the same token, Eq. (76) shows that the term $\mathcal{P}^2$ in (73) is of lower order in derivatives in comparison with the last three terms, because from (76) one can see that $\mathcal{P}$ is of lower order compared to $\mathcal{x}$ and $\mathcal{y}$. As a consequence, to the leading order in derivatives, the Lagrangian (73) is

$$\mathcal{L}_{\text{s; Kin}}^{(2)} = A\mathcal{y}^2 + B\mathcal{xy} + C\mathcal{x}^2, \qquad (79)$$

One can also verify that

$$D = -\frac{c_1^4 c_2}{c_3^4 c_4^4 c_6^2} D, \qquad (86)$$

where $D$ is defined in (51). In terms of $D$, the hyperbolicity condition found in (55) reads,

$$\frac{c_1^8}{16c_6^4 D} < 0. \qquad (87)$$

As long as $D < 0$, i.e., the hyperbolicity condition (84) is satisfied for $\tilde{\pi}$, the hyperbolicity condition is also satisfied for $\pi$. Moreover, for $D < 0$ the variables $\tilde{\pi}$ and $\pi$ and the kinetic matrices for $\tilde{\pi}$ and $\pi$ are related as,

$$S^{ab} = -\frac{c_1^4}{4c_6^2 D} \tilde{S}^{ab}, \qquad \pi = \frac{2c_6}{c_1^2} \sqrt{-D} \tilde{\pi}. \qquad (88)$$

where indices $a$ and $b$ are either $t$ or $r$.

where the subscript “Kin” stresses that only higher order terms (kinetic part) are left in the Lagrangian. We then introduce $\tilde{\pi}$ as

$$x = \tilde{\pi}', \qquad (80)$$

and from (78), we easily obtain

$$y = \dot{\tilde{\pi}}, \qquad (81)$$

where we set to zero the integration constant. Finally, substituting (80) and (81) in (79), we find the kinetic part of the Lagrangian for the scalar perturbations:

$$\mathcal{L}_{\text{s; Kin}}^{(2)} = -\frac{1}{2}(\tilde{S}^{tt}\dot{\tilde{\pi}}^2 + 2\tilde{S}^{tr}\dot{\tilde{\pi}}\tilde{\pi}' + \tilde{S}^{rr}\tilde{\pi}'^2), \qquad (82)$$

where

$$\tilde{S}^{tt} = -2A, \qquad \tilde{S}^{tr} = -B, \qquad \tilde{S}^{rr} = -2C. \qquad (83)$$

One can obtain the same result for the de Sitter black hole by following the above method rather than (47). Indeed, first of all, the hyperbolicity condition for (82) reads

$$D \equiv \tilde{S}^{tt}\tilde{S}^{rr} - (\tilde{S}^{tr})^2 < 0. \qquad (84)$$

The explicit expression for $D$ in terms of $c_i$ is given by

$$D = -\frac{c_1^4 c_2 \{c_2(c_4 c_7 - c_1 c_8)^2 - 2(c_2 c_6 - c_3 c_8)[c_4(c_2 c_1' + c_1 c_2' + 2c_4 c_5) - c_1 c_2 c_4]\}}{c_3^4 c_4^4}. \qquad (85)$$

The advantage of the Lagrangian (82) obtained here is that it also allows us to treat the case of stealth Schwarzschild black hole, for which the method of Sec. III A fails. Indeed, for the stealth solution it turns out that $D = D = 0$ (in other words, $H_0$ and $H_1$ are linearly dependent). However, the kinetic matrix $\tilde{S}^{ab}$ remains finite, see (83), while the kinetic matrix $S^{ab}$ diverges, as it can be seen from (88).

For the Lagrangian (82), the vanishing determinant of the kinetic matrix means that the equation of motion is parabolic (for all $r$). Per se, this fact does not necessarily mean that the perturbations are pathological on the considered background. For instance, in the case of the k-essence Lagrangian $\mathcal{L} = G_2(X)$, for solutions where $dG_2/dX = 0$ with timelike $\nabla^\mu \varphi$, the perturbations behave as dust, i.e., they are governed by a wave equation with $c_s^2 = 0$. The determinant of the kinetic matrix in this case is also zero, since only the $tt$ component of the kinetic matrix is nonvanishing. For the stealth solution (23)–(25), the kinetic matrix reads
$$ \tilde{S}^{ab} \propto \begin{pmatrix} \frac{\mu r}{(\mu-r)^2} & \frac{\sqrt{\mu r}}{\mu-r} \\ \frac{\sqrt{\mu r}}{\mu-r} & 1 \end{pmatrix}, \qquad (89) $$

where we defined $\mu = 2Gm$. Notice that, for $r \gg \mu$, all the terms of (89) apart from $\tilde{S}^{rr}$ tend to zero. The global factor of Eq. (89) may have any sign, depending on the parameters of the model and the (arbitrary) value of $q$ in (25). When this global factor is negative, the dynamics of the perturbations indeed corresponds to dust (i.e., a vanishing velocity, similarly to the example of k-essence described above), but the infinitely thin cone of propagation tends towards the $r$ axis. This, together with the fact that the graviton cone has a “usual” behavior at $r \to \infty$, makes the stealth solution pathological. It corresponds to a limit of panels (i) and (j) of Fig. 1 when the dashed (blue) cone is infinitely thin. On the other hand, for a positive global factor in Eq. (89), the scalar dynamics corresponds to the limit of panels (c) and (d) of Fig. 1 when the dashed (blue) cone totally opens, i.e., its sound velocity is infinite. In that case, the scalar field is no longer a propagating d.o.f.

## IV. CONCLUSIONS

In this paper, we have studied stability criteria for solutions in modified gravity theories. We then applied these criteria to establish the stability of certain hairy black holes whose hair is supplied by a dark energy scalar field [54,79,81,82].

Throughout this study, we focused on scalar-tensor theories, but the tools we have developed, as well as the stability criteria concerning Hamiltonian densities, are generically applicable in modified gravity theories. The starting ingredient for the applicability of our tools are multiple gravitational modes, a clear characteristic of theories going beyond GR. In order to treat the problem consistently, we have formulated the notion of causal cones, each of which is associated to a healthy propagating d.o.f. Indeed, the local existence of well-defined causal cones permits us to determine the healthy propagation of modes about an effective background solution. We saw that, unlike standard lore, the Hamiltonian densities associated to each of the modes do not suffice to exhibit an instability. The failure of the Hamiltonian criterion, in more complex background metrics, is due to the fact that it is not a scalar quantity. Each Hamiltonian density, associated to a propagating mode, depends on the particular coordinate system one is using. So although a Hamiltonian density which is bounded from below signals that the mode is stable, the converse is not true. Namely, a Hamiltonian density found to be unbounded from below in some coordinate system is inconclusive on instability. One may find a coordinate transformation rendering the Hamiltonian bounded from below as we saw explicitly in Sec. II.

Standard lore is recovered only for a class of “good” coordinate systems that are defined with respect to all the causal cones present in the system: scalar, graviton, matter, etc. Namely, these “good” coordinate systems exhibit a timelike coordinate common to all causal cones and spacelike coordinates for all causal cone exteriors. If such a coordinate system exists, then the Hamiltonian is indeed bounded from below and the modes are well behaved, propagating in a timelike direction with a hyperbolic operator. If not, then indeed the Hamiltonian for at least one of the modes is always unbounded from below, and the said mode presents a gradient or ghost instability.

The subtlety arises due to the complexity of the background solution. Indeed, the key point for our examples here is that the background scalar is space and time dependent. Then the causal cones can tilt and open up as we approach the horizon (event or cosmological). As a result, the original time (space) coordinate of the background metric may “leave” the interior (exterior) of a causal cone associated to some mode. This can lead to a misinterpretation of the Hamiltonian density associated to the initial coordinates, which leave the causal cone of the mode in question. We should emphasize that the failing Hamiltonian stability criterion is not due to the mixing of modes, as illustrated by the simple k-essence example of Sec. II. We may also consider the case of a theory including a $G_3$ Horndeski term (the DGP term) where the mixing and demixing of modes has been completely resolved [27] for an arbitrary background.<fnref n="13" /> In such a model were found self-accelerating vacua for the so-called kinetic gravity braiding (KGB) model [83] (see also [84]). In a standard FLRW coordinate system, where the dark energy scalar depends purely on cosmological time, such KGB self-accelerating solutions generically give (depending on the coupling constants of the theory) a stable vacuum, with an associated Hamiltonian which is bounded from below. When one considers the precise same stable vacua in a spherical coordinate system, where the metric is static but the scalar field now depends both on space and time, the same Hamiltonian density can be found to be unbounded from below. This, as we emphasized, is an artifact of a bad use of a coordinate system (here spherical) whereas the FLRW coordinates are indeed “good” (satisfying the causal cone criteria, its cosmological time remaining notably within the causal cones). This example demonstrates that misinterpretations related to Hamiltonians are not due to mixing of modes but, crucially, to the background depending (or not) on multiple coordinates. It is for this reason that we do not encounter problems with the Hamiltonian in FLRW systems e.g., or with static black holes (with a static scalar field). One therefore expects our analysis to be relevant for backgrounds (in modified gravity) with lesser symmetry, e.g., stationary backgrounds involving rotating black holes.

<footnote n="13">¹³This result is not known for $G_4$ theories for example.</footnote>
For stationary backgrounds e.g., the $\theta$ and $r$ dependent background effective metrics may again tilt and open up as we approach horizons. Also clearly our analysis could be used in vector-tensor theories where similar black holes to the one studied here have been found [85–87]. In any case, let us emphasize that there always exist “bad” coordinates in which the Hamiltonian density of a stable solution appears unbounded by below. It suffices that its time direction be outside at least one of its causal cones, or that one of its spatial directions be inside one of them. For instance, a mere exchange of $t$ and $x$ creates such a spurious pathology, whereas the physics is obviously unchanged. As usual in GR, one should never trust coordinate-dependent quantities.

We also underlined that when there exists a good coordinate system in which the total Hamiltonian density is bounded by below, then it may also be computed in other coordinate systems, but it no longer corresponds to the mere Hamiltonian. It becomes a linear combination of the energy and momentum densities, whose spatial integrals over the whole system are all conserved. In other words, the stability of the solution is still guaranteed by the boundedness by below of a conserved quantity, but this is no longer the mere energy which plays this role. The conditions for the existence of such a good coordinate system, i.e., for stability, are written in a covariant form in Eqs. (10)–(12).

Using the above tools, we have corrected a misinterpretation [55–58] in the literature about the said instability of a class of hairy black holes. It is true, as stated in [55], that the Hamiltonian for the graviton is always unbounded by below in Schwarzschild’s coordinates when approaching a horizon. However, at the same time, the graviton causal cone remains compatible (and may even coincide) with the matter causal cone under some conditions on the parameters defining the model. In other words, there exists coordinate systems where the Hamiltonian for the graviton is bounded by below. We completed this stability analysis by computing the scalar causal cone, thanks to the study of $\ell = 0$ perturbations. Again, we found that there exists a domain of parameters where the three causal cones share a common time and a common spacelike hypersurface. Hence, the class of hairy black holes studied here, that is quite generically encountered in Horndeski and beyond Horndeski theories [79,82], is free of ghost and gradient instability pathologies for a given range of parameters of the model. This is an important result considering the absence of stable hairy black holes in gravitational physics (see e.g., [88–91] for two celebrated cases). We have demonstrated this result for a particular Horndeski theory but the result will be similar for other cases and we even expect some cases will allow for self-tuning properties [53]. These are among some of the subjects to be treated in future studies.

## ACKNOWLEDGMENTS

We wish to thank A. Fabbri, R. Parentani, S. Robertson, V. Rubakov and A. Vikman for interesting discussions. The

authors acknowledge support from the research program “Programme national de cosmologie et galaxies” of the CNRS/INSU, France, from the call Défi InFIniTi, and from PRC CNRS/RFBR (2018–2020) No. 1985 “Modified gravity and black holes: Consistent models and experimental signatures”.

## APPENDIX: MONOPOLE PERTURBATION

The aim of this Appendix is to display the explicit expressions of the various coefficients used in our analysis of Sec. III A. Those entering Eq. (42) and later read

$$c_1 = -\frac{\beta q}{r} \varphi' \sqrt{\frac{B}{A}}, \qquad (A1)$$

$$c_2 = 2Bc_1, \qquad (A2)$$

$$c_3 = -\frac{1}{2r} \sqrt{\frac{B}{A}} (-2\zeta A + \beta q^2 - 3\beta A B \varphi'^2), \qquad (A3)$$

$$c_4 = \frac{2}{A} c_3, \qquad (A4)$$

$$c_5 = \frac{q^2}{4r^2} \frac{1}{\sqrt{AB}} [2\beta(1 - B - rB') + \eta r^2], \qquad (A5)$$

$$c_6 = -\frac{1}{4Ar^2} \sqrt{\frac{B}{A}} \left\{ \frac{1}{2} A \varphi'^2 [(2\beta - 12\beta B + \eta r^2)A - 12\beta B r A'] + rA' (\beta q^2 - 2\zeta A) - 2\zeta A^2 - \beta q^2 A \right\}, \qquad (A6)$$

$$c_7 = \frac{q^2}{4Ar^2} \frac{1}{\sqrt{AB}} [2\beta BrA' + A(2\beta - 2\beta B + \eta r^2)], \qquad (A7)$$

$$c_8 = \frac{q}{2Ar^2} \varphi' \sqrt{\frac{B}{A}} [-6\beta BrA' + A(2\beta - 6\beta B + \eta r^2)]. \qquad (A8)$$

The coefficients entering Eq. (44) and later read

$$\tilde{c}_1 = -\frac{c_1^2}{4c_6}, \qquad \tilde{c}_2 = -\frac{c_3^2}{4c_6}, \qquad (A9)$$

$$\tilde{c}_3 = \frac{c_1 c_3}{2c_6}, \qquad \tilde{c}_4 = -\frac{c_4^2}{4c_6}, \qquad (A10)$$

$$\tilde{c}_5 = \frac{c_4 c_3}{2c_6}, \qquad \tilde{c}_6 = -\frac{c_1 c_4}{2c_6}, \qquad (A11)$$

$$\tilde{c}_7 = \frac{2c_2 c_6 - c_8 c_3}{2c_6}, \qquad \tilde{c}_8 = \frac{c_8 c_1 - c_4 c_7}{2c_6}, \qquad (A12)$$

$$\tilde{c}_9 = -\frac{c_6(c_7 c_3' + c_3 c_7' - c_7^2 + 4c_5 c_6) - c_7 c_3 c_6'}{4c_6^2}, \qquad (A13)$$
$$ \tilde{c}_{10} = -\frac{c_7 c_8}{2c_6}, \qquad \tilde{c}_{11} = -\frac{c_8^2}{4c_6}. \qquad (A14) $$

Finally, the coefficients entering Eq. (71) read

$$ a_1 = \tilde{c}_1, \qquad a_2 = \frac{\tilde{c}_3}{2\tilde{c}_1}, \qquad (A15) $$

$$ a_3 = \frac{\tilde{c}_6}{2\tilde{c}_1}, \qquad a_4 = 0, \qquad (A16) $$

$$ a_5 = \frac{2\tilde{c}_1\tilde{c}_7 - \tilde{c}_8\tilde{c}_3}{\tilde{c}_6\tilde{c}_3}, \qquad a_6 = \frac{\tilde{c}_7}{\tilde{c}_3}, \qquad (A17) $$

$$ a_7 = \tilde{c}_9 - a_1 a_5^2 + (a_1 a_2 a_5)', \qquad a_8 = \tilde{c}_{11} - a_1 a_6^2, \qquad (A18) $$

$$ a_9 = \tilde{c}_{10} - 2a_1 a_6 a_5. \qquad (A19) $$

[1] C. M. Will, Living Rev. Relativity **17**, 4 (2014).

[2] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Phys. Rev. Lett. **119**, 161101 (2017).

[3] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Phys. Rev. Lett. **116**, 061102 (2016).

[4] D. Lovelock, J. Math. Phys. (N.Y.) **12**, 498 (1971).

[5] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. **513**, 1 (2012).

[6] C. M. Will, Theory and Experiment in Gravitational Physics (Cambridge University Press, Cambridge, England, 1993).

[7] T. Damour and G. Esposito-Farèse, Classical Quantum Gravity **9**, 2093 (1992).

[8] G. W. Horndeski, Int. J. Theor. Phys. **10**, 363 (1974).

[9] D. B. Fairlie, J. Govaerts, and A. Morozov, Nucl. Phys. B**373**, 214 (1992); D. B. Fairlie and J. Govaerts, J. Math. Phys. (N.Y.) **33**, 3543 (1992); Phys. Lett. B **281**, 49 (1992).

[10] A. Nicolis, R. Rattazzi, and E. Trincherini, Phys. Rev. D **79**, 064036 (2009).

[11] C. Deffayet, G. Esposito-Farèse, and A. Vikman, Phys. Rev. D **79**, 084003 (2009).

[12] C. Deffayet, S. Deser, and G. Esposito-Farèse, Phys. Rev. D **80**, 064015 (2009).

[13] C. de Rham and A. J. Tolley, J. Cosmol. Astropart. Phys. **05** (2010) 015.

[14] C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, Phys. Rev. D **84**, 064039 (2011).

[15] C. Deffayet and D. A. Steer, Classical Quantum Gravity **30**, 214006 (2013).

[16] G. W. Horndeski, J. Math. Phys. (N.Y.) **17**, 1980 (1976).

[17] C. Deffayet, S. Deser, and G. Esposito-Farèse, Phys. Rev. D **82**, 061501 (2010).

[18] L. Heisenberg, J. Cosmol. Astropart. Phys. **05** (2014) 015.

[19] G. Tasinato, J. High Energy Phys. **04** (2014) 067.

[20] R. Kase, M. Minamitsuji, and S. Tsujikawa, Phys. Lett. B **782**, 541 (2018).

[21] C. de Rham, Living Rev. Relativity **17**, 7 (2014).

[22] J. Zanelli, arXiv:hep-th/0502193.

[23] C. Charmousis, E. Kiritsis, and F. Nitti, J. High Energy Phys. **09** (2017) 031.

[24] C. Charmousis, Lect. Notes Phys. **892**, 25 (2015).

[25] G. R. Dvali, G. Gabadadze, and M. Porrati, Phys. Lett. B **485**, 208 (2000).

[26] A. Nicolis and R. Rattazzi, J. High Energy Phys. **06** (2004) 059.

[27] E. Babichev and G. Esposito-Farèse, Phys. Rev. D **87**, 044032 (2013).

[28] P. Horava, Phys. Rev. D **79**, 084008 (2009).

[29] T. Jacobson, Proc. Sci. QG-Ph (2007) 020

[30] M. Zumalacárregui and J. García-Bellido, Phys. Rev. D **89**, 064046 (2014).

[31] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, Phys. Rev. Lett. **114**, 211101 (2015).

[32] C. Lin, S. Mukohyama, R. Namba, and R. Saitou, J. Cosmol. Astropart. Phys. **10** (2014) 071.

[33] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, J. Cosmol. Astropart. Phys. **02** (2015) 018.

[34] C. Deffayet, G. Esposito-Farèse, and D. A. Steer, Phys. Rev. D **92**, 084013 (2015).

[35] D. Langlois and K. Noui, J. Cosmol. Astropart. Phys. **02** (2016) 034.

[36] D. Langlois and K. Noui, J. Cosmol. Astropart. Phys. **07** (2016) 016.

[37] M. Crisostomi, M. Hull, K. Koyama, and G. Tasinato, J. Cosmol. Astropart. Phys. **03** (2016) 038.

[38] M. Crisostomi, K. Koyama, and G. Tasinato, J. Cosmol. Astropart. Phys. **04** (2016) 044.

[39] J. B. Achour, D. Langlois, and K. Noui, Phys. Rev. D **93**, 124005 (2016).

[40] C. de Rham and A. Matas, J. Cosmol. Astropart. Phys. **06** (2016) 041.

[41] J. B. Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui, and G. Tasinato, J. High Energy Phys. **12** (2016) 100.

[42] H. Motohashi and T. Suyama, Phys. Rev. D **91**, 085009 (2015).

[43] H. Motohashi, K. Noui, T. Suyama, M. Yamaguchi, and D. Langlois, J. Cosmol. Astropart. Phys. **07** (2016) 033.

[44] T. P. Sotiriou, Lect. Notes Phys. **892**, 3 (2015).

[45] B. P. Abbott et al. (LIGO Scientific and Virgo and Fermi-GBM and INTEGRAL Collaborations), Astrophys. J. **848**, L13 (2017).

[46] J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. Lett. **119**, 251304 (2017).

[47] P. Creminelli and F. Vernizzi, Phys. Rev. Lett. **119**, 251302 (2017).

[48] J. Sakstein and B. Jain, Phys. Rev. Lett. **119**, 251303 (2017).
[49] L. Lombriser and A. Taylor, J. Cosmol. Astropart. Phys. 03 (2016) 031.

[50] L. Lombriser and N. A. Lima, Phys. Lett. B **765**, 382 (2017).

[51] N. Bartolo, P. Karmakar, S. Matarrese, and M. Scomparin, J. Cosmol. Astropart. Phys. 05 (2018) 048.

[52] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Phys. Rev. Lett. **119**, 251301 (2017).

[53] E. Babichev, C. Charmousis, G. Esposito-Farèse, and A. Lehébel, Phys. Rev. Lett. **120**, 241101 (2018).

[54] E. Babichev and C. Charmousis, J. High Energy Phys. 08 (2014) 106.

[55] H. Ogawa, T. Kobayashi, and T. Suyama, Phys. Rev. D **93**, 064078 (2016).

[56] K. Takahashi, T. Suyama, and T. Kobayashi, Phys. Rev. D **93**, 064068 (2016).

[57] K. Takahashi and T. Suyama, Phys. Rev. D **95**, 024034 (2017).

[58] R. Kase, M. Minamitsuji, S. Tsujikawa, and Y. L. Zhang, J. Cosmol. Astropart. Phys. 02 (2018) 048.

[59] A. Maselli, H. O. Silva, M. Minamitsuji, and E. Berti, Phys. Rev. D **93**, 124056 (2016).

[60] D. Bettoni and S. Liberati, Phys. Rev. D **88**, 084020 (2013).

[61] J. Bekenstein and M. Milgrom, Astrophys. J. **286**, 7 (1984).

[62] C. Armendariz-Picon, T. Damour, and V. F. Mukhanov, Phys. Lett. B **458**, 209 (1999).

[63] T. Chiba, T. Okabe, and M. Yamaguchi, Phys. Rev. D **62**, 023511 (2000).

[64] C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, Phys. Rev. D **63**, 103510 (2001).

[65] E. Babichev, V. F. Mukhanov, and A. Vikman, J. High Energy Phys. 09 (2006) 061.

[66] J. P. Bruneton, Phys. Rev. D **75**, 085013 (2007).

[67] J. P. Bruneton and G. Esposito-Farèse, Phys. Rev. D **76**, 124012 (2007).

[68] E. Babichev, V. Mukhanov, and A. Vikman, J. High Energy Phys. 02 (2008) 101.

[69] Y. Aharonov, A. Komar, and L. Susskind, Phys. Rev. **182**, 1400 (1969).

[70] F. J. Belinfante, Physica **6**, 887 (1939).

[71] A. Bandyopadhyay, Int. J. Theor. Phys. **38**, 1531 (1999).

[72] M. Rinaldi, Phys. Rev. D **86**, 084048 (2012).

[73] A. Anabalón, A. Cisterna, and J. Oliva, Phys. Rev. D **89**, 084050 (2014).

[74] M. Minamitsuji, Phys. Rev. D **89**, 064017 (2014).

[75] S. Appleby and E. V. Linder, J. Cosmol. Astropart. Phys. 03 (2012) 043.

[76] T. Kobayashi, H. Motohashi, and T. Suyama, Phys. Rev. D **85**, 084025 (2012); **96**, 109903(E) (2017).

[77] T. Kobayashi, H. Motohashi, and T. Suyama, Phys. Rev. D **89**, 084042 (2014).

[78] A. Ganguly, R. Gannouji, M. Gonzalez-Espinoza, and C. Pizarro-Moya, Classical Quantum Gravity **35**, 145008 (2018).

[79] E. Babichev and G. Esposito-Farèse, Phys. Rev. D **95**, 024020 (2017).

[80] T. Regge and J. A. Wheeler, Phys. Rev. **108**, 1063 (1957).

[81] C. Charmousis and D. Iosifidis, J. Phys. Conf. Ser. **600**, 012003 (2015).

[82] E. Babichev, C. Charmousis, and A. Lehébel, Classical Quantum Gravity **33**, 154002 (2016).

[83] C. Deffayet, O. Pujolas, I. Sawicki, and A. Vikman, J. Cosmol. Astropart. Phys. 10 (2010) 026.

[84] E. Babichev, C. Charmousis, A. Lehébel, and T. Moskalets, J. Cosmol. Astropart. Phys. 09 (2016) 011.

[85] J. Chagoya, G. Niz, and G. Tasinato, Classical Quantum Gravity **33**, 175007 (2016).

[86] M. Minamitsuji, Phys. Rev. D **94**, 084039 (2016).

[87] E. Babichev, C. Charmousis, and M. Hassaine, J. High Energy Phys. 05 (2017) 114.

[88] J. D. Bekenstein, Ann. Phys. (N.Y.) **91**, 75 (1975).

[89] K. A. Bronnikov and Y. N. Kireev, Phys. Lett. **67A**, 95 (1978).

[90] T. J. T. Harper, P. A. Thomas, E. Winstanley, and P. M. Young, Phys. Rev. D **70**, 064023 (2004).

[91] M. S. Volkov and D. V. Gal'tsov, Phys. Rep. **319**, 1 (1999).