Stability of generalized Einstein-Maxwell-scalar black holes

We study the stability of static black holes in generalized Einstein-Maxwell-scalar theories. We derive the master equations for the odd and even parity perturbations. The sufficient and necessary conditions for the stability of black holes under odd-parity perturbations are derived. We show that these conditions are usually not similar to energy conditions even in the simplest case of a minimally coupled scalar field. We obtain the necessary conditions for the stability of even-parity perturbations. We also derived the speed of propagation of the five degrees of freedom and obtained the class of theories for which all degrees of freedom propagate at the speed of light. Finally, we have applied our results to various black holes in nonlinear electrodynamics, scalar-tensor theories and Einstein-Maxwell-dilaton theory. For the latter, we have also calculated the quasinormal modes.


Introduction
Stability is one of the most important measures in physics, telling us how likely a solution could exist for a long period. It discriminates cosmological solutions from transient situations. Electrovacuum black holes (BHs) in general relativity (GR) are stable. Indeed after the seminal work by Regge and Wheeler [1], where they proved partially the stability of the Schwarzschild black hole, Zerilli [2] provided the stability equations for the even-parity perturbations. Later, following the same approach, Moncrief [3,4] derived the stability of the Reissner-Nordstrom black hole. But because of the impossibility to separate variables in the rotating case, Teukolsky [5,6] developed a new approach which helped to study the linear stability of the Kerr black hole. Even if the Newman-Penrose formalism used by Teukolsky turned out to be extremely difficult to implement for other rotating black holes, such as the Kerr-Newman, the metric approach to non-rotating solution has been widely used in various beyond GR models, see e.g. [7][8][9]. We see that stability is an essential analysis of any solution. On the other hand, thermodynamic stability is also widely studied. It is important to mention that (mechanical) stability should be a condition before any thermodynamic analysis. Indeed, if the black hole is mechanically unstable, any Hawking radiation would destabilize it and therefore render any thermodynamic analysis impossible. Finally we should also mention the interesting aspects of instabilities in nature which provide for example phase transitions. In this direction, any unstable black hole in an asymptotically AdS spacetime could be associated via the AdS/CFT to a phase transition. From a more astrophysical aspect, BHs turn out to be extremely useful to study models beyond GR, in a strong gravitational regime. As we know, BHs are extremely simple objects. They are described in the vacuum by mainly two parameters, the mass and the angular momentum (the electric charge being usually negligible). These parameters, which describe entirely the BH, can be measured, in particular via the frequency emitted during the last moments of BH mergers, the so-called quasinormal modes (QNMs). Any hairy BH could be discriminated against through these modes. QNMs are easily computed if their equation is known. The analysis performed in this paper will provide these equations.
As we mentioned, in this paper, we address the problem of black hole stability for a generic class of theories beyond GR. Modifying Einstein's theory of gravity is an easy game but building a new physical theory turns out to be complicated. Many modified gravity theories turn out to be highly constrained by the data, from the Brans-Dicke model to massive gravity à la dRGT. We should follow a simpler road. Could we define theories with sufficient generality? For that, some of the simple elements that we will consider are 1. Ostrogradsky instability free theory, and therefore we will assume second order differential equations for the fields 1 2. Well-posed Cauchy problem. Interestingly, this condition reduces drastically the form of the models in presence of a scalar field. Indeed, only K-essence non-minimally coupled survives to these restrictions [13][14][15].
Of course, beyond GR theories can't be reduced to only the presence of a scalar field. For example, the dimensional reduction or Kaluza Klein compactification decomposes the higher 1 Even if higher order differential equation for a scalar field turn out to be possible, see beyond Horndeski [10,11], these theories haven't been very successful in describing the Universe. More generically, degenerate theories [12] can be constructed which do not propagate more degrees of freedom even if higher derivatives are present, by considering an appropriate coupling to the metric which satisfies the degenerate condition.
dimensional metric into at least a vector field along with the scalar field. It is therefore natural to consider a generic model with these two additional degrees of freedom, keeping in mind that K-essence types of theories are strongly hyperbolic. The models studied in this paper can be regarded as a non-linear extension of the hairy black hole studied by Gibbons and Maeda [16] in the context of Kaluza-Klein theories, but also rediscovered by Garfinkle, Horowitz and Strominger [17] in the context of string theory.
In this work, we therefore investigate the linear stability of generalized Einstein-Maxwellscalar black holes assuming that they are static and without magnetic charge. We will derive some generic results and apply them to specific models. This paper is organized as follows.
In section 2, we describe the model and the background equations. In section 3, we summarize the black hole perturbation formalism that we will apply in sections 4 for the odd-parity perturbations and section 5 for the even-parity perturbations. Finally, we will apply our results to various known solutions in the literature, before final comments and conclusions.
To facilitate the use of our results, a Mathematica ® notebook is available online [18].

Background equations of motion
We are interested in generalized Einstein-Maxwell-dilaton theories described by the following action where R is the Ricci scalar, f 1 (φ) is a function of the scalar field φ, representing the nonminimal coupling to gravity, and f 2 (φ, X, F ) is a function of the scalar field, its kinetic energy X = − 1 2 ∇ µ φ∇ µ φ and F = − 1 4 F µν F µν where the field strength is defined as We will study the linear stability of static spherically symmetric spacetimes. For that, we will consider a background spacetime described by the following line element 2 Notice that we kept a general function C(r) for the area of the sphere of constant r and t, because various solutions in the literature are written in this form.
To satisfy the background symmetries of our spacetime, we will consider that the scalar field only depends on the radial coordinate,φ = φ(r) and the vector field takes the nongeneric formĀ This function could also be decomposed into transverse and longitudinal modes using Helmholtz's theorem [19]. But because we do not consider any mass term à la Proca A µ A µ , the componentĀ 1 will play no game in our equations. This background is not the most generic spherically symmetric spacetime 3 . In fact, we should only impose that the energy-momentum tensor T µν is invariant under the lie derivative along the generators of the group SO(3), or we could consider a less restrictive situation where the electromagnetic part of the energy-momentum tensor is spherically symmetric. Therefore we could consider A =Ā 0 (r)dt +P cos θdϕ whereP is the magnetic charge. We will consider in this paper only electric charge. It is a genuine restriction. Using our background metric, the kinetic terms for the scalar and vector fields take the simple form,X = −Bφ 2 /2 andF = BĀ 2 0 /(2A), where a prime represents a derivative with respect to r.
Replacing the metric (2.2), the scalar and vector fields in the action (2.1), we obtain the equations of motion (EOM) by varying the action with respect to the functions X = {A, B, C, φ,Ā 0 }. We obtain the equations of motion E X = 0 where indicates a derivative wrt r. From the last equation follows that Q(r) is constant; this expresses the existence of a conserved charge arising from the U (1) symmetry. It is useful to introduce the following variables to simplify the notations hence, the equation of motion for the scalar field is written as E φ = J + S = 0.

Perturbation formalism
In this section, we will briefly summarize the formalism for perturbations around a static spherically symmetric spacetime. We consider a background metricḡ µν and a perturbed metric g µν =ḡ µν + h µν . The background metric will be described by the line element (2.2). Under a reparametrization of the angles (θ, ϕ) into R 2 (θ, ϕ), R 3 (θ, ϕ), the 10 components of the perturbed metric g µν transform as scalar, vector or tensor. For example, the component g 00 transforms as a scalar while g 0i , where i = (2, 3), transforms as a vector. In fact, g 0i dtdx i → g 0i ∂ j R i dtdx j . In summary, we have three scalars (g 00 , g 01 , g 11 ), two vectors of dimension two (g 0i , g 1i ) and a rank two tensor (g ij ), where (i, j) take the values (2, 3). Each vector field has an orthogonal decomposition and the rank 2 tensor can be decomposed as from which, we can conclude easily that in four dimensions τ ij = 0 because it has 3 different components (τ 22 , τ 23 , τ 33 ) and 3 restrictions (∇ i τ i j = 0, τ i i = 0). In D = 4, perturbations around spherically symmetric spacetime can be described by scalar and vector pertubations, which are also known as odd or axial for vector perturbations [1] and even or polar for scalar perturbations [2]. As we have seen, in higher dimensions, we have an additional tensor mode [20].
Considering the symmetries of our background spacetime, all scalars can be decomposed in the basis of (scalar) spherical harmonics as and all vectors can be decomposed into a basis of vector spherical harmonics 4 where E ij = √ det γ ij , γ ij being the two-dimensional metric on the sphere and ij being the Levi-Civita symbol with 23 = 1. The functions (∇ i Y m , E j i ∇ j Y m ) define the basis of our vector spherical harmonics.
In summary, we have 3 vector perturbations and 7 scalar perturbations These perturbations are not independent. We can eliminate four of them using the coordinate transformation x µ → x µ + ξ µ where ξ µ are infinitesimal. This transformation can also be decomposed into scalar and vector perturbations, we have for the scalar part 14) and for the vector part Under this transformation, the perturbed metric transforms according to h µν → h µν − 2∇ (µ ξ ν) and therefore, we have where a dot indicates a derivative wrt t. We see that for vector perturbations and for ≥ 2, the only full gauge fixing condition is the Regge-Wheeler gauge defined by h From which we see that 3 gauges (in the gravity sector) are possible, defined by . We will use in our calculations, the first gauge.
These perturbations should be supplemented by the perturbation of the scalar field and the vector field. The scalar field φ has only a scalar perturbation 5 while the vector field perturbation (A µ =Ā µ + δA µ ) can be decomposed into a scalar and a vector part. The vector perturbation is while the scalar perturbation can be written as Because of the gauge freedom associated to the vector field 6 , we can set A (2) m to zero. We will see later that for = (0, 1) some of the perturbations are identically zero which allows us to use other gauges and simplify the problem. The formalism summarized in this section applies to scalar and vector perturbations when ≥ 2. 5 which transforms under a gauge transformation as δφ m → δφ m − Bφ R m . We will later use this condition. 6 We define the transformation m , that we will redefine as A In summary, we have two type of modes which can be studied separately at the linear order unless we have parity-violating terms [21] such as Pontryagin density * RR. Also for any spherically symmetric background, the final equations of perturbations are independent of the order m [22]. Therefore, without any loss of generality, we will assume m = 0.

Odd-parity perturbations
Over the static and spherically symmetric background (2.2), we will consider small perturbations of the fields defined by eqs. (3.5, 3.6, 3.7, 3.29) Since modes with different evolve independently, we focus on a specific mode and omit the index.

Second-order action for higher multipoles, ≥ 2
In this section we will focus in higher multipoles ( ≥ 2), since dipole mode = 1 require a special treatment. Expanding the action (2.1) to second order in perturbations 7 and performing integration over the sphere (θ, ϕ), we find where the coefficients are given by and λ = ( +1). On-shell, E A = E B = EĀ 0 = 0 which allows us to simplify these coefficients. Notice that for = 1, we have a 2 = 0 on-shell, this is why this mode has to be studied separately. The term inside square brackets in (4.2) can be rewritten in a more convenient way such that the Lagrangian is simplified to The first order perturbations vanish on-shell. 8 We have renamed for simplicity of notations h where we have defined, b 1 = a 1 − (C a 3 ) /C and b 2 = a 5 + 2a 3 a 4 C /C ≡ λEĀ 0 /C. On-shell, the last coefficient vanishes, b 2 = 0.
Since no time derivative of h 0 appears in the Lagrangian, variation with respect to it, yields a constraint equation. However, because of the presence of h 0 in the action, this constraint results in a second-order ordinary differential equation which cannot be immediately solved for h 0 . To overcome this obstacle, we follow the procedure described in [23,24], we introduce an auxiliary field q(t, r) and define the following Lagrangian We can easily check that by substituting the EOM for q(t, r) into Eq. (4.5), we recover the original Lagrangian. Now, varying (4.5) with respect to h 0 and h 1 leads to these expressions relate the metric elements h 0 and h 1 to the auxiliary field q. Once q(t, r) is known, h 0 and h 1 are easily obtained, and thus all vector perturbations in the Regge-Wheeler gauge are obtained. Substituting these expressions into the Lagrangian and performing integration by parts one finds which explicitly shows that taking into account the vector sector, we only have two degrees of freedom, one associated with gravitational perturbation and another related to vector perturbation of the electromagnetic field. The coefficients of the above Lagrangian, in terms of the variables (4.3), are We see that in the absence of electromagnetic perturbations, the Lagrangian (4.8) becomes L (2) odd = α 1q 2 + β 1 q 2 + γ 1 q 2 , from which we should at least impose the no-ghost condition α 1 > 0 translating into a 2 < 0 and therefore f 1 > 0. Similarly, in the absence of the gravitational perturbations, we have L from which we impose the condition α 2 > 0 which means Q/Ā 0 > 0 and using eq.(2.8) gives f 2,F > 0.

Master equation
To arrive at our final result, it is convenient to rescale the variables as where (V g , V e ) represent respectively the vector perturbations associated to the gravitational and electromagnetic sector. Notice that we have assumed the previous conditions f 1 > 0 and f 2,F > 0. Finally, introducing the tortoise coordinate, dr = √ AB dr * we find where Ψ = (V g , V e ) t , and V ij are the components of a 2×2 symmetric matrix and are given, after substituting relations (4.3), by 14) The introduction of the functions (S 1 , S 2 ) will be more transparent in the stability analysis section. Variation of the action (4.11) with respect to Ψ gives us a wave-like equation Both modes propagate at the speed of light. Since matrix V is not diagonal, the previous equation is a set of two coupled differential equations. But the system can be diagonalized because the eigenvalues are independent of the radial coordinate 9 .

Stability analysis
Stability means that no perturbation grows unbounded in time. For that, we will Fourier transform our variables ( Ψ → e −iωt Ψ) such as eq.(4.15) is recast as where H = −∂ 2 r * + V. The frequency ω 2 appears as the eigenvalues of the operator H. Unstable modes are equivalent to purely imaginary modes ω 2 < 0 (see e.g. [24]), therefore the stability of the spacetime is related to the positivity of the operator H, namely that H has no negative spectra. To prove the stability, let us define the inner product where ψ = (ψ 1 , ψ 2 ) T and ξ = (ξ 1 , ξ 2 ) T . Stability means that the operator H is a positive self-adjoint operator in L 2 (r * ), the Hilbert space of square integrable functions of r * . Therefore, we need to prove the positivity defined as This condition will imply that given a well-behaved initial data, of compact support, χ remains bounded for all time. This is a sufficient condition. The rigorous and complete proof of the stability related to equations of the form (4.15) can be found in [25,26] using spectral theory. We have whereχ i (i = 1, 2) represent the complex conjugate of χ i and should not be confused with background quantities. In the second line, we have neglected the boundary term χ 1 ∂ r * χ 1 +χ 2 ∂ r * χ 2 coming from the integration by parts, because we assumed χ 1 and χ 2 to be smooth functions of compact support, while in the third line we have neglected the boundary term S 1 |χ 1 | 2 + S 2 |χ 2 | 2 . Therefore, we conclude that the black hole is stable under vector perturbations for ≥ 2 if the no-ghost condition is satisfied, namely 10 (4.20) Notice that these conditions are not equivalent to energy conditions (see Appendix B).

Stability of dipole perturbation, = 1
We should first mention that for = 0, spherical harmonics are constant and therefore the three vector perturbations of the metric and the vector perturbation of the electromagnetic field are identically zero as we can see from their definition (3.5, 3.6, 3.7, 3.29).
For the dipole perturbation, = 1 and assuming m = 0 because the final result is independent of the azimuthal angle, we have Y 0 We don't need to use the Regge-Wheeler gauge to eliminate this perturbation. We have seen that 10 It was also noticed that in Horndeski theory, the black hole is stable if the no-ghost and hyperbolicity conditions are satisfied [24]. Notice that for Lovelock black holes, the no-ghost and hyperbolicity conditions could be satisfied while the black hole is unstable [27].
Therefore, we can use this gauge freedom to fix h 1 = 0 but we have a residual gauge freedom defined as Λ → Λ + f (t)C(r). Let us consider the Lagrangian (4.4) which is valid also for = 1. After imposing the background equations, we get a 2 = b 1 = 0; and therefore 11 Variation of this action with respect to h 1 and h 0 gives uṡ where we have defined which solution is given by E = J /C(r), where J is an integration constant 12 . Considering our gauge freedom, we fix h 1 = 0, and integrating the previous equation, we obtain where F (t) is a constant of integration. This last term can be eliminated with the help of the residual gauge freedom, giving finally The variation of (4.23) wrt A v gives Defining a new variable a v = f 2,F A v and using the tortoise coordinate, we find where V 22 is the coefficient defined in eq.(4.13) for = 1, i.e. λ = 2.
The solution of eq.(4.29) is the sum of a particular solution that we can consider as a function of r only and a homogeneous solution. In order to understand this solution, let us consider the simple case where we have perturbed Reissner-Nordström in general relativity. For that we take f 1 = 2, f 2 = 4F and Q = 4q, from which we get A = B = 1−2M/r+q 2 /r 2 and C = r 2 . For this spacetime, a particular solution of eq.(4.29) is 30) 11 We will fix h1 = 0 after using the constraint derived from its variation because the gauge is not totally fixed by this condition. 12 We hope that this will not confuse the reader because of the function defined in eq.(2.9).
in other words, the Kerr-Newman solution at the first order in J if we define J = −12J. Therefore, the particular solution of eq.(4.29) should describe a slow rotating black hole while the homogeneous solution describes the propagation of the electromagnetic field in our spacetime [28]. This propagation is stable because of the potential V 22 which can be easily deformed to a positive potential using the function S 2 (r). In fact, as we have shown in Sec. (4.3), the operator associated is definite positive.
In conclusion, we have shown that for any theory described by the action (2.1) and considering the conditions f 1 > 0 and f 2,F > 0, the static spherically symmetric spacetime described by a generic metric (2.2) is stable under vector perturbations.

Even-parity perturbations
In this section, we will study the scalar (even, polar) perturbations of the action (2.1). We follow the formalism that we have summarized in Sec. (3) 5.1 Second-order action for higher multipoles, ≥ 2 Following exactly the same procedure, we arrive at where the second-order Lagrangian is In what follows, we will reduce this Lagrangian and rewrite it in terms of three variables, representing the remaining three DOF of the theory. The variation of the action with respect to A 0 and A 1 gives Introducing the new variable 13 defined as

5)
13 which will be associated to the scalar part perturbation of the electromagnetic field.
we obtain Let us notice from Lagrangian (5.2) that H 1 has no derivatives, so we can integrate out this non-propagating field by using its own equation of motion Notice that the coefficient b 1 is zero on-shell, so we have dropped this term in the previous equation. At this point, the Lagrangian (5.2) only depends on the variables H 0 , H 2 , α, S e and δφ. Variation of the action with respect to H 0 gives us Taking the combination E 1 − a 11 /(2d 1 )S e (t, r), and using the identity 2a 1 − a 2 11 /2d 1 = 0, we obtain a relation that sets a constraint for the other fields, but this is not an algebraic constraint; however, in order to resolve this issue, we perform a field redefinition and use a new variable 14 S g defined by Using this relation, the constraint becomes an algebraic equation for α, which can be solved where T and P are given by Replacing this definition into E 1 , we can rewrite H 0 in terms of the other variables and their derivatives (this is a large expression, and at this point, it is unnecessary to write it). Substituting this relations into (5.2) we obtain a Lagrangian that only depends on variables S g (t, r), S e (t, r) and δφ(t, r) given by 14 which will be associated to the scalar part perturbation of the gravitational field.
The above Lagrangian can be rewritten in matrix form where χ = (S g , S e , δφ) T and (K ij , L ij , M ij ) are symmetric matrices, while D ij is antisymmetric 15 . The no-ghost condition requires the matrix K to be positive definite; for this purpose, we employ the Sylvester's criterion, which gives We obtain that the first two conditions reduce to where we defined Imposing the conditions that we obtained in the odd-parity sector, namely f 1 > 0 and f 2,F > 0, we get that the last term in (5.15) is always positive, so the previous two relations are reduced to a single constraint: λP 1 − f 1 > 0. Finally, the third condition is given by We recover the no-ghost condition computed in [29] for a general scalartensor theory.
Using the background equations, we find and assuming stability of odd-parity perturbations, viz. f 1 > 0, we obtain

Speed of propagation of scalar perturbations
We know that our model consists of five degrees of freedom, where two are related to gravity, two others to the electromagnetic field and finally one degree of freedom is associated to the scalar field. They are decomposed around a spherically symmetric spacetime as two vector perturbations and three scalar perturbations. As we have seen, vector perturbations propagate at the speed of light. To complete the analysis, we need to obtain the speed propagation of the even perturbations. In that direction, let us consider that the solution is of the form 16 χ ∝ e i(ωt−kr) . Considering the small scale limit, the dispersion relation obtained from (5.13) can be written as det ω 2 K ij − k 2 L ij = 0. The propagation speed c r in proper time can be derived by substituting ω = √ ABc r k into the dispersion relation. Solving for c 2 r , we obtain the propagation of the three degrees of freedom: The propagating mode different than the speed of light is related to the scalar field sector. Indeed, in the case where the vector field is absent, we have that P 3 = 2Σ, and c 2 r3 reduces to the propagation speed of the scalar field given in [29].
Using previous relations, we get 17 .

(5.26)
It is important to emphasize that these velocities are very weakly constrained by gravitational waves. Certainly, between the merger of two compact objects and us, the wave propagates mostly over a FLRW spacetime than a spherically symmetric background. But, we should also consider that a model described by our action is just an effective low energy description of some more fundamental theory. For that, we should impose standard conditions such as Lorentz invariance, unitarity, analyticity. Even if not proved generically, it was shown in various situations and around different backgrounds that these conditions imply nonsuperluminal propagation (see e.g. [31][32][33]). Hereafter, we will not consider this 16 We are calculating only the radial speed of propagation. 17 Notice that if we consider a theory for which f2,F + 2F f2,F F = 0, viz.
one of the degrees of freedom propagates at infinite speed, which could be similar to the Cuscuton [30], for which this perturbation do not carry any microscopic information.
condition but the much more restrictive, possibly more interesting, condition c 2 r 3 = 1, which translates from eq.(5.26) into We conclude that any action given by the form below will propagate at the speed of light five degrees of freedom where (f 1 , f 2 , f 3 ) are generic functions.

Master equation
From Lagrangian (5.13), we obtain the equations of motion where we have used χ(t, r) → e −iωt χ(r). Making a change of variable χ i (r) → S ij (r)Φ j (r), and changing to tortoise coordinate dr = √ ABdr * we get where the matrix S is solution of 18 and We see from eq.(5.30) that not all degrees of freedom propagate at the speed of light. But if we consider the condition (5.27), we have L ij = ABK ij , which implies where V is the matrix potential, which expression can be read off directly from (5.30) Given the matrix S, the potential V can be easily derived and the stability studied. Unfortunately, in the generic case, we were not able to obtain the matrix S but we will see in the following sections, various examples where these calculations can be performed easily.
Following the same procedure as in the vector perturbation sector, we define the operator H ij = −δ ij d 2 dr * + V ij , and we obtain where we have as usual performed an integration by parts and eliminated the boundary term. Even if a generic stability condition can't be obtained, we can discuss sufficient stability conditions. In fact, a sufficient but not necessary condition of the stability is obtained if the potential is definite positive. On the contrary, if there is a trial function 19 Φ 0 , such that (Φ 0 , HΦ 0 ) < 0, the spacetime is unstable [34].
In the case 20 where f 2,XF = 0, we found for large λ Therefore we conclude that if the integral is negative for sufficiently large λ which implies the instability.

Second order Lagrangian for = 0
Similarly to the vector perturbations, the generic analysis we have performed in the previous section is valid only for higher modes, ≥ 2. For scalar perturbations, we will see that for = 0, we have only one perturbation related to the scalar field, the breathing mode, which is the spherically symmetric perturbation and for = 1 we will have 2 perturbations related to the scalar field and the electromagnetic field. Using the relations derived in Sec.3, we have, for = 0, without defining any gauge, h 0i = h 1i = 0 and only K 00 survives in the expression of h ij . We can choose the freedom on ξ 1 to fix K 00 = 0. We still have the 19 Because the lowest eigenvalue ω0 of the spectrum of H gives ∀ Φ, ω 2 0 ≤ (Φ, HΦ)/(Φ, Φ). Therefore if a trial function Φ0 is such that (Φ0, HΦ0) < 0, we have ω 2 0 < 0, for a normalizable trial function. 20 Notice that this condition includes most, if not all models studied in the literature for generic Einstein-Maxwell-dilaton theories. freedom on choosing ξ 0 . Using these relations, we obtain from eq.(5.2) using = 0 (λ = 0) L (2) even, =0 = a 1 H 2 0 + H 0 a 2 H 2 + a 4 H 2 + a 7 δφ + a 9 δφ + a 10 δφ + a 11 (A 0 −Ȧ 1 ) Variation of the action with respect to A 0 and A 1 gives us the relations ∂ t S e (t, r) = ∂ r S e (t, r) = 0, where S e (t, r) is defined in (5.5); thus, S e (t, r) ≡ C 1 is constant.
Variation wrt H 0 , and using the previous relation gives Therefore, we conclude that where C 2 is a new integration constant. This equation will be used to eliminate the variable H 2 .
We can use the freedom on ξ 0 , to partially fix the gauge by considering H 1 = 0 which would also gives us an additional freedom viz. H 0 → H 0 − f (t). Considering now the variation wrt to H 2 and fixing this gauge, we obtain an expression for H 0 which can be integrated and the integration constant g(t) eliminated by the remaining gauge freedom. These expressions are not specially illuminating and therefore will not be written but we will demonstrate the effect of this mode on a particular solution. Finally, the equation wrt to δφ gives the equation where we have defined From this equation, we can see that the no-ghost condition, K 0 > 0, is the same as the one for higher multipoles, i.e., equation (5.19). On the other hand, in the limit of large wavenumber k, we get that the velocity of the propagation is which also coincides with c 2 r 3 given in Eq. (5.22). Since only the scalar wave is excited in the monopole perturbation, this result allows us to interpret c 2 r 1 and c 2 r 2 as the propagation speed of gravitational and vector perturbations, and c 2 r 3 as the propagation of the scalar waves.
The eq.(5.44) admits a particular and homogeneous solution. The homogeneous solution will describe the propagation of the spherical scalar wave in the fixed background metric while the particular solution should modify the constants of the metric. For that, let us consider the electric GM-GHS black hole (see section 6.4) with Using the previous equations, we found that a particular solution of eq.(5.44) is a constant which gives us (H 0 , H 2 ) and therefore . That gives us A = B = 1 − 2M/r − α/r. Therefore, only the mass has been modified. Also we have which shifts the electric charge F 01 → F 01 − (M C 1 + 4αq)/8M r 2 . In conclusion, we expect that for any spacetime, the particular solution shifts the constants of the spacetime while the homogeneous solution describes a scalar wave propagating in a fixed background.

Second order Lagrangian for = 1
Considering now the dipole perturbation, we see that for = 1 and therefore Y 0 1 ∝ cos θ, we have h 22 = (K(t, r) − G(t, r)) cos θ , (5.52) h 33 = (K(t, r) − G(t, r)) cos θ sin 2 θ , (5.53) Therefore, the metric perturbation h ij depends on K and G only through the combination K − G. We can use one function of the coordinate transformation, ξ 2 , to set K − G = 0 (ξ 3 = 0). We also use the transformation ξ 0 to define β(t, r) = 0. Thus, we have a remaining degree of freedom that we can use to set δφ(t, r) = 0. As we have seen, δφ → δφ − Bφ R.
Using the transformation of coordinate ξ 1 sets the scalar field to zero. In this case, the gauge is totally fixed and therefore we can proceed as usual. The Lagrangian for the dipole mode is obtained taking λ = 2 ( = 1) in (5.2) Performing similar calculations to those we did in the section of higher multipoles, after taking δφ = 0, we get a second-order Lagrangian similar to (5.12). Again, after redefining the variables, we obtain that the Lagrangian takes the canonical form In this case, χ is a two-component vector representing the two propagating degrees of freedom, K 1 , L 1 and M 1 are 2 × 2 symmetric matrices and D 1 is antisymmetric. The no-ghost conditions, positivity of the matrix K 1 ij , are the same as in the higher multipole case, i.e., equations (5.15) and (5.16), which are reduced to the stability condition (5.19) after taking λ = 2. Also, from this Lagrangian, we obtain that the propagation speeds along the radial direction are given by , (5.57) which again coincide with the propagation speed in the case of higher multipoles Eq. (5.22). We conclude that the degree of freedom associated to the scalar field travels at the same speed, regardless of the multipole order.
6 Application to specific models Because f 2,F = 4 > 0, we conclude that the vector perturbations are stable. For completeness, let us write the perturbation potential for this sector which corresponds exactly to the potential derived 21 originally in [3]. Considering the scalar perturbation sector, we need first to eliminate one row and one column of our matrices because of the absence of the scalar field. Then, we need to solve the eq. (5.31) which gives r + c 4 (6.5) 21 Our perturbations (Vg, Ve) correspond respectively to (πg,π f ) of [3].
where (c 1 , c 2 , c 3 , c 4 ) are constants of integration, which should be chosen freely as soon as c 1 c 4 −c 2 c 3 = 0 in order for S to be invertible. Finally, the calculation of the scalar potential is trivial and we found that our potential V is similar to the original potential 22 derived in [4] with a change of a constant basis matrix. Both matrices describe the same problem. In fact, if we have we can always define a transformation Φ = P Ψ and the similar potential W = P VP −1 such that if P is a constant matrix. From the potential, the stability can be easily obtained because of the positivity of the matrix.

Nonlinear electrodynamics
In this section, we generalize the previous section by considering nonlinear electrodynamics (NED) black holes, given by the action where L is an arbitrary function of F = −F µν F µν /4. From the background equations, it is easy to find that which can be easily integrated to where α is a constant of integration. In the case where we would be interested only to the background solution, we could always reparametrize our time coordinate such that α = 1 but because the perturbations are time dependent we will keep that constant. The vector perturbations are stable if the condition (4.20) is satisfied (f 1 (φ) = 1/2 > 0), viz.
(6.11) 22 Where our variables (Sg, Se) correspond respectively to (Q, H) of [4] For completeness of the vector perturbations, the potential (4.12,4.13,4.14) is given by In [28], the perturbations for this model were studied in a gauge invariant formalism. We find exactly the same potential when C(r) = r 2 .
Focusing on the even-parity sector, since we have only two degrees of freedom, the stability conditions reduce to (for simplicity of the expression, we have eliminated the obvious positive factors) giving the same condition as the odd-parity sector, namely L ,F > 0. Following with the even-parity sector, we compute the matrix S (5.31) where as in the previous case, the constants (c 1 , c 2 , c 3 , c 4 ) can be freely defined as soon as c 1 c 4 −c 2 c 3 = 0. The expression of the potential in the generic case is long and not necessarily instructive. It can be easily obtained from the eq.(5.35). It is anyway interesting to derive the most used expression in the literature, namely A(r) = B(r) and C(r) = r 2 which implies from eq.(6.10), α = 1/4. We will take without any restriction c 2 = c 3 = 0 and c 4 = 2Q(λ − 2)c 1 . We get , where we have defined a(r) = rA − 2A, B = λ − 2 + 4Q 2 /r 2 L ,F and κ = −rĀ 0 /(2Ā 0 ) = L ,F /(L ,F + 2F L ,F F ). This potential is the same as in [28]. From which we can easily derive the stability condition for a specific model or calculate the QNMs. It is important to notice that the stability will depend only on the metric potential A and its derivatives. In fact, using the equations of motion, we can eliminate the electric potentialĀ 0 and the Lagrangian L, reducing the potential matrix to , (6.21) From which we conclude that different theories (Lagrangians) having the same black hole solution will have the same stability. In that direction, we found that considering A = B = 1 − 2M/r + q 2 /r 2 , the stability potential reduces to the Moncrief potential 23 which implies the stability of the Reissner-Nordström metric independently of the nonlinear electrodynamics theory considered.

Bardeen black hole
As an application, let us consider one of the first singularity free black hole proposed by Bardeen [35]. This metric was shown to be an exact solution of Einstein equations with a nonlinear magnetic monopole [36]. It is described by the metric 2M r 2 (r 2 + q 2 ) 3/2 (6.23) and C = r 2 . Because we have focused in this paper on electric source, we know from the FP duality 24 , that a solution can be generated by two different Lagrangians, one theory with a magnetic field and the other with an electric field 25 . It is easy to find that a Lagrangian 23 Our potential V can be written as V = P W P −1 with W the Moncrief potential and P = 0 1 −1 0 , therefore they are similar. 24 The Lagrangian L(F ) can also be written in terms of the field P = Pµν P µν /4 where Pµν = Fµν L,F [37]. This is the so-called P framework. The equations in the P framework are equivalent to the equations in the F framework by performing a transformation while the metric remains unchanged [38]. Therefore, the FP duality connects theories with different Lagrangians but similar metric. A purely magnetic solution in the P framework will correspond to a purely electric solution in the F framework or vice versa. This transformation exists [28] if L,F (L,F + 2F L,F F ) = 0. 25 It would be interesting to see if a magnetically charged BH follows the same stability condition as the electrically charged BH. Indeed, we found previously that in the electrically charged case, the stability condition is independent of the Lagrangian.
with electric field generating the Bardeen spacetime is A 0 (r) = α r 5 (r 2 + q 2 ) 5/2 , (6.24) where α is a constant. From which, we could obtain L(F ) but because the expression is numerical and not analytical, we will just mention that such Lagrangian exists 26 . We can see from the electric potential, that we will not recover the Coulomb potential at large distances, which could also be seen from the metric that reduces to 1 − 2M/r + 3M q 2 /r 3 for large r. Therefore, even at large distances, the black hole is not similar to Reissner-Nordström spacetime. But interestingly, long distance quantum corrections to the Schwarzschild black hole goes like 1/r 3 [39,40] which is similar to the Bardeen solution if we identify q 2 to 62 /45π Looking now to the stability of this black hole, we found that odd-parity perturbations are stable if f 2,F = 3M/αr 2 is positive which implies the trivially satisfied condition M α > 0.
Considering finally the even parity sector, we have analysed the stability numerically from the potential (6.20,6.21,6.22). For that, we redefined the variables r →rq, M →M q such thatM remains the only free parameter. Considering the normalized mass in the range 0 <M < 10 and for = 2, ..., 10, we found that the potential is definite positive and therefore the black hole is stable in this range.

Hayward black hole
An other interesting singular free black hole is defined as [41] This black hole is known as the Hayward black hole. It reproduces the Schwarzschild spacetime at large distances with corrections O(1/r 4 ) which are not similar to loop quantum contribution as noticed in the Bardeen solution. But similarly, the black hole is regular at the origin and has a de Sitter core, A 1 − r 2 /l 2 .
As in the previous case, the solution can be derived from various theories. We will consider nonlinear electrodynamics. It is easy to find that a Lagrangian exists, numerically, 26 It was shown in [38] that a Lagrangian having a Maxwell asymptotic "does not admit a static, spherically symmetric solution with a regular center and a nonzero electric charge", which seems to contradict our result. In fact, the Lagrangian we found is a multivalued function, therefore we would need different Lagrangians for different ranges of the radial coordinate r. But it is single-valued if we restrict our analysis to the exterior region, which is the interesting area for our analysis. Of course, these models can't be continued until the singularity and therefore loose their interest as singular free models. which gives this metricĀ where α is a constant. Within this theory, we would like to know if the black hole is stable. The odd-parity sector is trivially stable because f 2,F = (r 3 + 2M l 2 ) 3 /18α 2 M 2 l 2 r 7 > 0 As for the Bardeen black hole, the even-parity perturbations will be studied numerically. Interestingly this spacetime is invariant under a scaling transformation [42]. Therefore, we can renormalize our variables, r → l √ M r, M → M 3/2 l, which reduces the solution to only one parameter. We have checked the stability for 0 < (M/l) 2/3 ≤ 10 and 2 ≤ l ≤ 10. We found that the potential is definite positive and therefore the black hole is stable in this range.

Born-Infeld gravity
Considering open superstring theory, loop contributions lead to a Lagrangian of the Born-Infeld type where b is related to the tension of the D-brane [43]. In four dimensions, the determinant can be expanded into [44] where F is the dual of the electromagnetic tensor, ( F ) µν = 1 2 µνρσ F ρσ . Focusing from now on an electric source and adding gravity, the action becomes The static spherically symmetric electrically charged black hole was obtained in [45] where g(r) = − r 4 + b 2 Q 2 dr which can be expressed in terms of hypergeometric functions.
This solution has interesting features such as recovering the Reissner-Nordström solution at large distances, with highly suppressed corrections πb . Also notice that around the singularity, at r = 0, the metric takes the form A −2M/r but because M is not the ADM mass, it can take any sign causing a time-like or space-like singularity according to the sign of m − Γ( 1 4 ) 2 Q 3/2 /6 √ πb. But contrary to [46], we conclude that the singularity is unavoidable, even if M = 0, because the curvature scalar is R = 4Q/br 2 around the singularity.
It is trivial to see that odd-parity perturbations are stable because f 1 = 2 > 0 and f 2,F = 4/ 1 + b 2 F/2 > 0. Turning now, to the even-parity perturbations, we have checked for a large range of parameters and found that the eigenvalues of the matrix V (6.20,6.21,6.22) are positive outside the horizon, from which we can conclude the stability of this black hole.

Scalar-tensor theory
After studying spacetimes with the presence of an electric charge, we will consider another particular case of our analysis for which only the scalar field is present, viz.
where we have redefined the functions of the Lagrangian in a more standard notations of a K-field non-minimally coupled to gravity also known as K-inflation [47] or K-essence [48]. Since Nordström and his attempt to write a scalar theory of gravity [49], scalar fields remained a very dynamical field of research from phenomenological models to more advanced and motivated theories. It has been motivated as early as Kaluza and Klein where a scalar field appears in the compactification of a fifth dimension. Similarly, string theory behaves in the low-energy limit as general relativity coupled to a scalar field, the dilaton, and a totally antisymmetric field strength H µνρ [50]. Also, we should mention that the dynamics of the tachyon field near the minimum of its potential is given by [51,52] K(φ, X) = −F (φ) 1 + (∂ µ φ) 2 , or the attempts to describe dark energy by the ghost condensation scenario [53] or inflation [54]. It would be impossible to list all the different models with a scalar field, but we should mention that these models since Brans and Dicke have been largely used as a model testing of general relativity [55].
As in the previous cases, the odd-parity sector is trivial. Indeed, we only need to satisfy the condition F (φ) > 0 outside the black hole. For example, any theory expressed in the string frame [56] where F (φ) = e −φ would be trivially stable under odd-parity perturbations.
A generic sufficient condition for the stability of even-parity perturbations has not been found in closed form but we will work out particular examples. This analysis can be easily generalized to any BH within this class of theories with the Mathematica ® notebook file provided in [18].
Among all these theories, it was shown that Horndeski theories and Lovelock are weakly hyperbolic but might fail to be strongly hyperbolic [13][14][15]. The only sub-class which is strongly hyperbolic are the K-field models defined in eq. (6.37). Unfortunately for most of these models, a non trivial black hole seems to be impossible. For example, it was shown that if the model is shift symmetric, viz. f 1 = 1 and f 2 = f 2 (X), the only solution is Schwarzschild [57]. For that, it was assumed that the area of constant r-spheres should neither be infinite nor zero at the horizon. Violating this condition, the so-called cold black holes (with zero surface gravity) can be constructed [58]. For that, we need to violate the null energy condition by considering f 1 = 1 and f 2 = −X. We see that the scalar field has the "wrong" sign for the kinetic term. The solution is where m is the ADM mass and k is related to the scalar charge. Indeed, we have at infinity φ √ 2m 2 −2k 2 r . Even if very interesting, this black hole violates the no-ghost condition (5.14). In fact, it is easy to show that det(K) < 0. Maybe a stable black hole could be constructed in that direction by considering a non-linear kinetic term. We should here mention that various papers on wormholes consider a ghost-like scalar field. While performing a stability analysis at the equation level and not the action level, we could reach the erroneous conclusion of a normal behaviour of that solution. But because the Hamiltonian is unbounded from below, any coupling to matter such as a detector or any nonlinear regime of the theory would couple to the negative energy modes and render the theory unstable. These theories are pathological.
Another way to escape the no-go theorem previously mentioned is to brake the shift invariance and assume a potential. We know that important no-hair theorems restrict the existence of non-trivial black holes, see e.g. [59][60][61][62]. One particular solution is given by a scalar field non-minimally coupled to the scalar curvature, the so-called BBMB black hole [63][64][65]. The action is given by f 1 (φ) = 1 2 − φ 2 6 and f 2 (φ, X) = X while the line element is given by and φ(r) = √ 3M/(r − M ), where M is the ADM mass. The horizon is located at r = M where the scalar field diverges but with a regular geometry of the horizon. It is also a cold black hole which connects smoothly to the Minkowsky space M = 0 but never to the Schwarzschild spacetime. In fact, the black hole has the same causal structure than the extremal Reissner-Nordstrom black hole. We conclude easily that the black hole is unstable because it violated the condition f 1 (φ) > 0 in the region M < r < 2M , which is consistent with the radial perturbations performed in [66].

BH in Einstein-Maxwell-dilaton theory
Finally, we consider the presence of all fields, namely the metric, the scalar and vector field. In this context, an interesting solution was derived in Einstein-Maxwell-dilaton (EMd) theory described by the Lagrangian [16,17] where a is the dilaton coupling constant. The model appears as a low energy limit of heterotic string theory for a = 1 and as the dimensionally reduced five dimensional vacuum Einstein action in the Einstein frame, for a = √ 3. The static and spherically symmetric black hole solution in the case where the Maxwell field is electric is given by

43)
Here r + is the position of the event horizon and r − corresponds to a curvature singularity R µνρσ R µνρσ ∝ (r − r − ) −2(1+3a 2 )(1+a 2 ) except for a = 0 where the solution reduces to Reissner-Nordström metric. The two parameters (r + , r − ) are related to the ADM mass M and the electric charge Q by [67] Given the background solution, we can perform the analysis of perturbation theory. As a first step, we focus on odd-parity perturbations. The no-ghost conditions (4.20) are trivially satisfied since f 1 (φ) = 2 > 0 and f 2,F (φ, X, F ) = 4e −2aφ > 0. We conclude that the oddparity perturbations are stable. From Eqs.(4.12) -(4.14), we get that the components of the potential associated with this type of perturbations are where we have defined . (6.50) It can be checked that the potential matrix has constant eigenvalues; therefore, we can diagonalize this matrix to decouple the system of wave-like equations. After this diagonalization, the potentials are As a specific case, setting a = 1, we recover the potential computed in [68]. Now, we turn to the problem of even-parity perturbations. The no-ghost condition (5.19) is trivially satisfied Therefore, even-parity perturbations are well defined. It can also be check that all perturbations propagate at the speed of light. In order to study completely the stability of these modes, we need to compute the matrix S (5.31) in order to obtain the potential associated to these perturbations. Unfortunately, in these type of models, the integration of eq. (5.31) is not easy 27 and therefore we will rely to a numerical analysis. For that, we integrate numerically eq.(5.31) assuming any initial conditions as soon as det(S) = 0. From which we can easily obtain the potential (5.35). We found that the eigenvectors of the matrix (5.35) −S −1 K −1 BS are constants and therefore the system has the same eigenvalues than the matrix −K −1 B. As a consequence, we do not need to calculate the matrix S, but only the eigenvalues of the matrix −K −1 B. Each eigenvalues will be the effective potential of one of the three type of perturbations. From these effective potentials, we have calculated the QNMs using the sixth order WKB method [69][70][71][72]. We see in the Fig.(1), the real and imaginary part of the QNMs associated to this black hole compared to the Reissner-Nordström solution. We see that contrary to the latter, the electric charge breaks the isospectrality 28 . From the numerical results, we found that Im(ω) < 0, for a large number of , and therefore we conclude that the black hole is linearly stable. As an example, we have represented, in Figs.(1,2), the real and imaginary parts of the QNM for a = 1 and a = 0.5 respectively. These results coincide with [68] for a = 1 and with [73,74] in the case a = 0.5 27 The matrices obtained are extremely large and therefore not very useful. 28 Same spectrum of quasinormal modes for even and odd perturbations.  We have studied the stability of static spherically symmetric black holes in generalized Einstein-Maxwell-scalar theories with an electrically charged Maxwell field. We have obtained the master equations in the odd and even parity sectors. We found that the ghost free conditions reduce to ,φ + 2f 1 f 2,X > 0 . (7.1) We found that four degrees of freedom propagate at the speed of light, while a degree of freedom associated to the scalar field could propagate faster or slower than the speed of light, its expression is given in eq.(5.26), from which we obtained the subclass of generalized Einstein-Maxwell-scalar theories where all degrees of freedom propagate at the speed of light eq. (5.28). Imposing the ghost free conditions, we found that odd-parity perturbations are always stable and we obtained explicitly the effective potential matrix. Considering the even-parity sector, we derived the master equation and the procedure to calculate the effective potential matrix for a given model. We also found that assuming ghost free conditions, the perturbations are unstable if f 2,F + 2F f 2,F F < 0 for models with f 2,XF = 0. Finally, we have applied this formalism for a large number of BH for which we could easily obtain the stability conditions and calculate the QNMs. We could recover all previously derived results in the literature which guarantees the correctness of our calculations. As an application, the paper can be used to study the QNMs of gravitational waves, electromagnetic radiation and scalar radiation. The presence of these three fields appears naturally in higher dimensional theories of gravity such as supergravity or string theory. Interestingly, we didn't find any stable hairy black hole in scalar tensor theories.
In order to facilitate future analysis of black hole stability in generalized Einstein-Maxwell-dilaton gravity, all the perturbed equations are listed in a Mathematica ® notebook available online [18].

A Even-parity coefficients
The coefficients which appear in the second-order Lagrangian (5.2) are

B Energy conditions
Considering the generalized Einstein-Maxwell-dilaton action (2.1), we obtain the GR-like equation G µν = T µν with Because the energy momentum tensor is diagonal we can interpret these components as the energy density and the three principal pressures in that frame in which the energymomentum tensor is diagonal. We have where ρ is the energy density, P r is the radial pressure and P t is the tangential pressure. The energy conditions can be geometric if related directly to the curvature tensor and physical if related to the fluid energy-momentum tensor. In our case, because we will assume that the additional fields can be interpreted as an additional fluid, the equation is reduced to a GR type equation and therefore the geometric and the physical energy conditions are equivalent and described by the following relations Null energy condition (NEC) ρ + P r ≥ 0 , ρ + P t ≥ 0 (B.3)