title: Quasinormal modes of hot, cold and bald Einstein-Maxwell-scalar black holes
authors: Jose Luis Blázquez-Salcedo Carlos A. R. Herdeiro Sarah Kahlen Jutta Kunz Alexandre M. Pombo Eugen Radu
text: 
  
  
  
  
  Quasinormal modes 
  of hot, cold and bald
  
   
   Einstein-Maxwell-scalar black holes
  
  
  
  
  Jose Luis Blázquez-Salcedo^,
  Carlos A. R. Herdeiro^,  
  
  Sarah Kahlen^, Jutta Kunz^,
  Alexandre M. Pombo^ and 
  Eugen Radu^ 
  
  
  
  
  
  ^Institut für  Physik, Universität Oldenburg, Postfach 2503,
  D-26111 Oldenburg, Germany
  
  
  
  
  ^Departamento de Matemática da Universidade de Aveiro and  
   Centre for Research and Development  in Mathematics and Applications (CIDMA), 
   Campus de Santiago, 3810-183 Aveiro, Portugal
  
  
  
  
  
  
   
  
  
  
  Einstein-Maxwell-scalar models allow for different classes of black hole solutions, depending on the non-minimal coupling function f(ϕ) employed, between the scalar field  and the  Maxwell invariant.
  
  Here, we address the linear mode stability of the black hole solutions obtained
  recently for a quartic coupling function, f(ϕ)=1+αϕ^4 <cit.>.
  Besides the bald Reissner-Nordström solutions, this coupling allows for two branches of
  scalarized black holes, termed cold and hot, respectively.
  For these three branches of black holes we calculate the spectrum of quasinormal modes.
  It consists of polar scalar-led modes, polar and axial electromagnetic-led modes,
  and polar and axial gravitational-led modes. We demonstrate that the only unstable mode
  present is the radial scalar-led mode of the cold branch. 
  Consequently, the bald Reissner-Nordström branch and the hot scalarized branch are both mode-stable.
  The non-trivial scalar field in the scalarized background solutions leads to
  the breaking of the degeneracy between axial and polar modes present for Reissner-Nordström solutions.
  This isospectrality is only slightly broken on the cold branch, but it is strongly broken on the hot branch.
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  The phenomenon of spontaneous scalarization of black holes has received much interest in recent years.
  In certain scalar-tensor theories, 
  
  the well-known black holes of General Relativity (GR) remain solutions of the
  field equations, while in certain regions of the parameter space, additional branches of black hole solutions arise,
  that are endowed with scalar hair.
  Spontaneous scalarization can, for instance, be charge-induced,
  when a scalar field is suitably coupled to the Maxwell invariant F^2
  <cit.>.
  Whereas the onset of the instability is universal, the
  
   properties of the resulting scalarized black holes and the branch structure of the solutions then depend 
  significantly on the coupling function f(ϕ)
  
  <cit.>. 
  
  
  
  Einstein-Maxwell-scalar (EMs) models include two distinctive classes, depending on the choice of f(ϕ): models that have black holes with scalar hair only (and do not allow the GR solutions) or models that allow both  the GR solutions and new hairy black holes. The latter splits into  two further sub-classes: models wherein the GR black  holes are unstable, in  some region of parameter space, against the tachyonic instability  that promotes scalarization and models wherein the GR black holes are never afflicted by this instability. The last two sub-classes  have been labelled as classes IIA and IIB, respectively,  in <cit.>.  
  
  
  
  
  Many previous studies on EMs BHs in the last few years have been motivated by charge-induced spontaneous scalarization, and thus have focused on black holes of class IIA 
  
  <cit.>. 
  In the static case, the branch of Reissner-Nordström (RN) black holes then develops 
  a zero mode at a critical charge to mass ratio q_cr(α),
  that depends on the coupling strength α. 
  At this bifurcation point, a new branch of scalarized black holes emerges. 
  The zero mode of the RN black hole branch 
  then turns into an unstable mode, while the new scalarized branch may or may not be stable,
  depending on the coupling function
  <cit.>.
  
  Of course, unstable excited scalarized branches exist as well.
  
  Let us note at this point that this observed pattern of branches and instabilities
  is very similar to the one seen in curvature-induced spontaneously scalarized black holes
  
  <cit.>. 
  
  Recently, we have considered an  EMs model 
  
  employing the quartic 
  coupling function f(ϕ)=1+αϕ^4, which qualifies as class IIB,
  since the RN branch does not become unstable against scalar perturbations
  <cit.>.
  Instead, we have observed the following interesting pattern:
  Close to the extremal RN solution, corresponding to the mass to charge ratio q=1,
  a first branch of scalarized black holes emerges.
  This first branch then exists in the range q_min(α) ≤ q ≤ 1.
  At q_min(α), it bifurcates with a second branch of scalarized black holes,
  which extends throughout the interval q_min(α) ≤ q ≤ q_max(α),
  and ends in an extremal singular solution at
  q_max(α) > 1.
  Considering the properties of these three branches, we have termed the RN branch
  as the bald branch, since it does not carry scalar hair, while we have termed the first
  and the second branches as the cold and the hot branches, respectively, 
  according to the black hole horizon temperatures.
  
  In this first study <cit.>, 
  we have also addressed the radial modes of the black holes for these three branches,
  since the radial modes signal instabilities with respect to scalar perturbations
  and thus the onset of scalar hair.
  While our analysis has shown that there are no unstable radial modes on the RN branch 
  and on the hot scalarized branch, we have found that
  the cold scalarized branch develops an unstable mode close to the point q=1. This instability is present throughout the interval q_min(α) < q < 1 and ends with a zero mode
  at the bifurcation point q_min(α) with the hot scalarized branch.
  Thus, the cold scalarized branch is clearly unstable. However, it remained an open issue
  
  whether there are really two stable coexisting branches, the bald RN branch and the hot scalarized branch.
  This question has motivated our present investigation.
  
  Here, we study linear mode stability of the black holes on all three branches.
  To that end, we calculate the lowest quasinormal modes for each type of mode.
  In particular, since the theory has scalar and vector fields coupled to gravity,
  we have to consider the perturbations of all these fields.
  In the presence of non-trivial background fields, 
  i.e., when the black hole solutions carry scalar hair and electromagnetic charge,
  the different types of perturbations generically couple to each other,
  leading to scalar-led, electromagnetic-led and gravitational-led modes
  instead of pure scalar, electromagnetic or gravitational modes,
  which one would find for a Schwarzschild black hole, for instance.
  Since parity even (polar) and parity odd (axial) perturbations decouple,
  we arrive at the following set of modes when expanding in spherical harmonics:
  polar scalar-led l ≥ 0 modes, axial and polar electromagnetic-led l ≥ 1 modes,
  and axial and polar gravitational-led l ≥ 2 modes.
  
  In Section II, we present the EMs theory studied
  and the general set of equations. We specify the Ansatz for the spherically symmetric
  background solutions and give the resulting set of ordinary differential equations (ODEs).
  We discuss the asymptotic expansions for asymptotically flat black hole solutions 
  with a regular horizon in Section  III. Here, we also recall basic properties of the three branches of
  black hole solutions. In Section IV, we formulate the sets of perturbation equations
  for spherical, axial and polar perturbations, deferring more details to the Appendix.
  We present our numerical results for the quasinormal modes in Section V.
  Here, we show that no further unstable modes arise on any of the three branches.
  Moreover, we discuss the breaking of isospectrality,
  i.e., the splitting of the axial and polar electromagnetic-led and gravitational-led modes,
  caused by the presence of the scalar hair. We conclude in Section V.
  
  
  
  
  § EMS THEORY
  
  
  We consider EMs theory described by the action
  
      𝒮=∫ d^4 x √(-g)[R-2∂ _μϕ∂ ^μϕ -f (ϕ) F_μν F^μν] ,
  
  where R is the Ricci scalar, ϕ is a real scalar field, and F_μν is the Maxwell field strength tensor.
  The coupling between the scalar and Maxwell fields is determined by the coupling function f(ϕ),
  for which we assume a quartic dependence,
  
      f(ϕ) = 1 +αϕ^4  .
  
  For a positive coupling constant α,
  the global minimum of the coupling function is at ϕ=0.
  
  The set of coupled field equations is obtained via the variational principle and reads
  
      R_μν-1/2g_μνR=
      T_μν^ϕ+T_μν^EM , 
      ∇_μ(√(-g) f(ϕ) F^μν )
      =0  , 
      1/√(-g)∂_μ
       (√(-g)g^μν∂_νϕ )
       = ḟ(ϕ) F_μνF^μν ,
  
  where ḟ(ϕ) = d  f(ϕ) /d ϕ, and
  T_μν^ϕ and T_μν^EM
  are the scalar stress-energy tensor
  and the electromagnetic stress-energy tensor, respectively.
  
      T_μν^ϕ≡1/2∂_μϕ∂_νϕ
      -1/2g_μν( 1/2(∂_αϕ)^2+V(ϕ)  )  , 
      
      T_μν^EM≡2 f(ϕ)
      (F_μαF_ν^  α
      -1/4g_μνF^2) .
  
  
  To construct static spherically symmetric EMs black holes, 
  we employ the line element 
  
      ds^2=-g(r) dt^2 + dr^2/1-2m(r)/r +r^2(dθ^2+sin^2 θ dφ^2)  ,
  
  where the metric functions g and m depend on the radial coordinate r.
  The black holes are supposed to carry electric charge and, in the scalarized case, also scalar charge.
  We therefore parametrize the gauge potential and the scalar field by
  
      A   =     a_0(r) dt  , 
      ϕ   =   ϕ_0(r)  ,
   
  where a_0 and ϕ_0 are the electric and the scalar field function,
  respectively, which depend on the radial coordinate r.
  
  Inserting this Ansatz into the general set of EMs equations (<ref>)-(<ref>),
  we obtain the following set of ODEs
  for the functions g, m, a_0 and ϕ_0: 
  
      ∂_rδ   =   -r/4 ( ∂_rϕ_0)^2 , 
      ∂_r m    =   1/2  f(ϕ_0) e^2 δr^2( ∂_ra_0  )^2+r/8(r-2m)  (∂_rϕ_0)^2 , 
      ∂^2_rϕ_0   =   e^2δ r/r-2m ( r(∂_rϕ_0)  f(ϕ_0) -2ḟ(ϕ_0) )
      ( ∂_r a_0)^2
             +2 m-r/(r-2m)r (∂_rϕ_0)  , 
      ∂^2_r a_0    =   (1/4 (∂_rϕ_0) ^2r
      -ḟ(ϕ_0)/f(ϕ_0) (∂_rϕ_0)-2/r) (∂_r a_0)  ,
  
  where we have introduced the function δ(r) via g=(1-2m/r)e^-2δ.
  Inserting the first equation into the last yields a
  first integral for the electromagnetic field 
  
      ∂_r a_0 = Q/f(ϕ_0) e^δr^2 ,
  
  where Q is the electric charge of the black holes.
  With this first integral, the above set of equations can be simplified.
  
  
  
  § EMS BLACK HOLES
  
  
  We here briefly recall the properties of static spherically symmetric electrically charged
  EMs black hole solutions with the quartic coupling function (<ref>).
  First of all, the RN black hole, which is also a solution of the EMs equations, is given by
  
      g=1-2M/r+Q^2/r^2 , 
      m=M-Q^2/2r , 
      a_0 = -Q/r , 
      ϕ_0 = 0  .
   
  
  The scalarized EMs solutions are obtained numerically <cit.>.
  Their asymptotic behavior yields their global charges
  
      g    =   1-2M/r+ O(r^-2)  , 
      
      m    =   M-(Q^2+Q_s^2/4)1/2r+O(r^-2)  , 
      ϕ_0   =   Q_s/r + (M Q_s - Q^2ḟ(0)/f(0))1/r^2 + O(r^-3)  , 
      
      a_0    =   -Q/r +  Q Q_s/2r^2ḟ(0)/f(0) + O(r^-3)  ,
  
  with black hole mass M, electric charge Q, and scalar charge Q_s.
  Note that there is no conservation law for the scalar field, and that
  the existence of a horizon imposes a non-trivial relation Q_s = Q_s(M,Q)
  <cit.>.
  
  Close to the horizon r=r_H, the global charges have the expansion
  
      g    =   g_1 (r-r_H) + O((r-r_H)^2)  , 
      
      m    =   r_H/2+Q^2/2r_H^2 f(ϕ_H)(r-r_H) + O((r-r_H)^2)  , 
      ϕ_0    =   ϕ_H-2 Q^2 ḟ(ϕ_H)/r_H f(ϕ_H) (r_H^2 f(ϕ_H) -Q^2)(r-r_H) +O((r-r_H)^2)  , 
      
      a_0    =   -Ψ_H+Q√(g_1/r_H f(ϕ_H) (r_H^2 f(ϕ_H) -Q^2))(r-r_H)+O((r-r_H)^2)  ,
  
  where ϕ_0(r_H)=ϕ_H and a_0(∞)-a_0(r_H)=Ψ_H
  are the value of the scalar field and the electrostatic potential at the horizon.
  Further physically relevant horizon properties
  are, for instance, the temperature T_H and the horizon area A_H, respectively given by
  
      T_H   =   1/4π√(g_1/r_H^3  f(ϕ_H) (r_H^2   -Q^2 f(ϕ_H))) , 
      
      A_H   =    4π r_H^2  .
  
  
  
  
  
  The EMs black hole solutions can be characterized by three dimensionless quantities:
  by their charge to mass ratio q, their reduced horizon area a_H and their
  reduced horizon temperature t_H:
  
      q≡Q/M ,       a_H≡A_H/16π M ^2 ,     t_H≡ 8π M T_H   .
  
  
  In the following we will focus on the case with α=200, since for other values of the coupling constant the general properties of the black holes are qualitatively similar.
  
  We illustrate the EMs black hole branches in Fig. <ref>,
  where we exhibit their reduced area a_H (left) and reduced temperature t_H (right) versus their
  charge to mass ratio q <cit.>.
  The RN black holes are shown in black. The first scalarized branch emerges close to the
  extremal RN black hole solution. Along the branch, the mass to charge ratio q decreases,
  while the reduced area a_H and temperature t_H increase.
  At a minimal value of q, the first branch bifurcates with the second branch.
  Along the second branch, a_H decreases while t_H increases monotonically
  with increasing q. Thus, the first branch represents the cold (blue) branch,
  and the second branch the hot (red) branch.
  
  
  
  
  
  § LINEAR PERTURBATION THEORY
  
  
  Previously, we have addressed the radial stability of these EMs black holes,
  by looking for radially unstable modes <cit.>.
  Our analysis has revealed radial stability for the RN branch and for the hot scalarized branch.
  In contrast, for the cold scalarized branch, we have found an unstable radial mode.
  Here, our goal is more ambitious, since we want to fully clarify the
  
   linear mode stability of the
  RN and the hot scalarized branch. As we shall see these are both stable.
  
  To show linear mode stability, we will evaluate the quasinormal mode spectrum
  of the black holes on these branches. For completeness, we will consider the spectrum on all 
  three branches, including the unstable cold scalarized branch.
  The symmetry of the background solutions suggests to consider
  perturbations for the following three cases: 
  purely spherical perturbations (l=0), axial or odd-parity ((-1)^l+1) perturbations,
  and polar or even-parity ((-1)^l) perturbations.
  
  
  
   §.§ Spherical perturbations
  
  
  To study the spherical perturbations, we have to perturb all the fields: the scalar field,
  the electromagnetic field and the metric. 
  For the scalar and the electromagnetic field,
  we introduce the perturbation functions ϕ_1 and F_a_0,
  and employ the Ansatz
  
      ϕ=ϕ_0(r) + ϵ e^-i ω tϕ_1(r)  , 
      
      A=   a_0(r) (1+ϵ e^-i ω t F_a_0(r)) dt  .
   
  For the metric, we introduce the perturbations functions F_t and F_r
  and employ the Ansatz
  
      ds^2=-g(r) (1+ϵ e^-i ω t F_t(r)) dt^2 + 1+ϵ e^-i ω t F_r(r)/1-2m(r)/rdr^2 +r^2(dθ^2+sin^2 θ dφ^2)  .
  
  We use ϵ as the control parameter in the linear expansion. 
  The complex quantity ω=ω_R+iω_I   
  corresponds to the sought-after eigenvalue of the respective quasinormal mode.
  Its real part is the oscillation frequency; its imaginary part is the damping rate.
  
  
  Inserting the Ansatz (<ref>)-(<ref>)
  into the general set of field equations 
  (<ref>)-(<ref>) 
  and utilizing the set of equations for the background functions (<ref>) 
  leads to the Master equation for the spherical perturbations.
  This is a single Schrödinger-like ODE
  for the function Z=rϕ_1:
  
      d^2Z/d(R^*)^2   =    (U_0(r) - ω^2)Z  ,
  
  which involves the tortoise coordinate R^*, where
  
      ∂_r R^*    =   1/√(g(1-2m/r)) ,
  
  and the potential U_0, reading
  
      U_0(r)   =   r-2m/2r^5 f(ϕ_0) e^2δ[  (Q^2-r^2 f(ϕ_0) )  (r∂_rϕ_0)^2
      -4ḟ(ϕ_0)/f(ϕ_0) Q^2 r ∂_r ϕ_0 
      . 
         +   . 
      2Q^2 (2(ḟ(ϕ_0)/f(ϕ_0))^2-1  )
      +4rm f(ϕ_0) ]  .
  
  
  To determine the quasinormal modes, we need to impose
  an adequate set of boundary conditions.
  As the perturbation propages toward infinity, 
  we have to impose an outgoing wave behavior,
  i.e., 
  
      Z =A_ϕ^+e^iω R^*(1+iM/2ω1/r^2 + O(r^-3))    ,
  
  for r→∞.
  As the perturbation propagates toward the horizon,
  we have to impose an ingoing wave behavior,
  i.e., 
  
      Z    =   A_ϕ^-e^-iω R^*
            ( 1 + 2r_H f(ϕ_H)
      (i (ḟ(ϕ_H)/f(ϕ_H))^2
                      +e^δ_Hω r_H  ) Q^2  
                      -f(ϕ_H) e^δ_Hω r_H^3/( f(ϕ_H)r_H^2-Q^2 ) 
          ( -iQ^2 + f(ϕ_H) r_H^2
           (i+2r_He^δ_Hω) ) (r-r_H) 
      + O((r-r_H)^2) )   ,
  
  for r → r_H.
  In these expansions, the quantities A_ϕ^± denote arbitrary amplitudes
  for the perturbation, while all other terms are 
  fixed by the background solution.
  
  
  
  
   §.§ Axial perturbations
  
  
  When considering axial perturbations, the scalar field does not get affected
  due to symmetry. Here, only the electromagnetic field and the metric enter.
  An appropriate Ansatz for the electromagnetic field involves
  the perturbation function W_2:
  
      A= a_0(r) dt -ϵ W_2(r)e^-iω t∂_φY_lm(θ,φ)/sinθdθ
      + ϵ W_2(r)e^-iω tsinθ∂_θY_lm(θ,φ) dφ ,
  
  where Y_lm are the standard spherical harmonics. 
  The corresponding Ansatz for the metric
  
      ds^2   =   -g(r) dt^2 
      + 1/1-2m(r)/rdr^2 +r^2(dθ^2+sin^2 θ dφ^2)
      
             + 2 ϵ h_0(r)e^-iω t∂_φY_lm(θ,φ)/sinθ dtdθ +
      2 ϵ h_0(r)e^-iω tsinθ∂_θY_lm(θ,φ) dtdφ
             + 2 ϵ h_1(r)e^-iω t∂_φY_lm(θ,φ)/sinθ drdθ 
      +
      2 ϵ h_1(r)e^-iω tsinθ∂_θY_lm(θ,φ) drdφ  ,
  
  introduces the perturbation functions h_0 and h_1.
  
  
  Inserting this Ansatz into the field equations 
  (<ref>)-(<ref>)
  leads to a set of coupled differential equations,
  which is shown in the Appendix.
  It involves first order equations for 
  the metric perturbation functions h_0 and h_1,
  and a second order equation for 
  the electromagnetic perturbation function W_2.
  This system can be put into the form
  
      ∂_r Ψ_A = M_AΨ_A ,
  
  where Ψ_A denotes the perturbation functions in the form
  
      Ψ_A = [ [     h_0;     h_1;     W_2; ∂_r W_2 ]]  ,
  
  and M_A denotes a 4× 4 matrix
  which contains the background functions,
  the angular number l, 
  and the eigenvalue ω of the quasinormal mode.
  To determine the axial quasinormal modes,
  as for the radial ones, we have to impose the proper set
  of boundary conditions at infinity (outgoing) and at
  the horizon (ingoing).
  These can be found in the Appendix.
  
  
  
  
   §.§ Polar perturbations
  
  
  Since the radial perturbations are a set of polar perturbations (l=0),
  it is clear that the general polar perturbations involve again all the fields.
  For the scalar field, we introduce the perturbation function ϕ_1
  and the Ansatz
  
      ϕ   =   ϕ_0(r) + ϵ e^-i ω tϕ_1(r) Y_lm(θ,φ)  .
   
  For the electromagnetic field, we employ the perturbation functions a_1, W_1 and V_1
  and the Ansatz
  
      A   =    ( a_0(r)+ϵ e^-i ω t a_1(r) Y_lm(θ,φ)) dt 
      + ϵ W_1(r) e^-i ω t Y_lm(θ,φ) dr
       
            + ϵ V_1(r) e^-i ω t∂_θY_lm(θ,φ) dθ + ϵ V_1(r) e^-i ω t∂_φY_lm(θ,φ) dφ .
   
  Finally, for the metric, we introduce the perturbation functions N, H_1, L and T,
  and employ the Ansatz 
  
      ds^2   =   -g(r) (1+ϵ e^-i ω t N(r) Y_lm(θ,φ)) dt^2 
      - 2ϵ e^-iω t H_1(r) Y_lm(θ,φ) dtdr
      + 1-ϵ e^-i ω t L(r)Y_lm(θ,φ)/1-2m(r)/rdr^2
      
            +(r^2-2ϵ e^-i ω t T(r)Y_lm(θ,φ))(dθ^2+sin^2 θ dφ^2)  .
  
  
  Inserting this Ansatz into the field equations 
  (<ref>)-(<ref>) 
  then again results in a set of coupled differential equations,
  which are shown in the Appendix.
  In fact, this system of equations can be simplified
  when one introduces the new functions F_0, F_1 and F_2,
  that are defined in terms of the perturbation functions W_1, V_1 and a_1
  (as discussed in the Appendix).
  The resulting Master equations
  
  can again be put in vectorial form,
  
      ∂_r Ψ_P = M_P Ψ_P  .
   
  The vector now contains 6 components:
  
      Ψ_P = [ [     H_1;       T;     F_0;     F_1;     ϕ_1; ∂_r ϕ_1 ]]  .
  
  The 6 × 6  matrix M_P
  depends again on the background functions,
  the angular number l and the eigenvalue ω of the quasinormal mode.
  
  We remark that
  for a RN black hole, the equations for the scalar perturbations are decoupled.
  However, the equations for the electromagnetic perturbations couple with the equations
  for the metric perturbations.
  When the charged black hole also carries scalar hair,
  all equations are coupled to each other.
  
  Again, we must impose the proper boundary conditions
  at infinity (outgoing) and at the horizon (ingoing)
  to determine the quasinormal modes.
  As in the axial case, these are shown in the Appendix.
  
  
  
  
  
  § QUASINORMAL MODE SPECTRUM
  
  
  In the following, we briefly discuss the method used to extract the quasinormal mode
  spectrum of the EMs black holes and recall the nomenclature for the modes.
  Then, we present our numerical results for the l=0, l=1 and l=2 cases
  and discuss isospectrality breaking.
  
  
  
  
  
   §.§ Numerics and nomenclature
  
  
  We calculate the quasinormal modes of the black holes on all three branches,
  bald, cold and hot,
  making sure that
  the respective spectra match,
  when the background solutions get close.
  Since the RN quasinormal modes are well-known
  
  (see e.g., <cit.>), 
  this provides an independent check of the numerics used, as well as a first guess for the spectrum of the cold branch, which is typically close to the RN branch for large values of the electric charge.
  All calculations are performed for a coupling constant α=200.
  
  In a first step, we obtain the background solutions with high precision,
  solving numerically the set of  ODEs (<ref>),
  subject to the boundary conditions following from the expansions
  at infinity (<ref>) and at the horizon (<ref>).
  For this purpose, we employ the solver COLSYS <cit.>,
  which uses a collocation method for boundary-value ODEs together with 
  a damped  Newton method of quasi-linearization. 
  The problem is linearized and solved at each iteration step, 
  employing a spline collocation at Gaussian points. 
  The solver features an adaptive mesh selection procedure, 
  refining the grid until the required accuracy is reached.
  
  Once the background solutions are known, we follow the
  procedure that is analogous to the one we have used before
  to calculate the quasinormal modes of hairy black holes
  <cit.>. 
  We split the space-time into two regions,
  the inner region r_H+ϵ_H ≤ r ≤ r_J, 
  and the outer region r_J ≤ r ≤ r_∞. 
  In the inner region, we impose the respective ingoing wave behavior; 
  
  in the outer region, we impose the respective outgoing wave behavior.
  
  Then, we calculate sets of linearly independent solutions numerically
  and match them at the common border r_J of the two regions.
  The eigenvalue ω of the quasinormal modes is found
  when the functions and their derivatives are continuous at the matching point r_J.
  
  
  We follow a common nomenclature for quasinormal modes, that is used
  when scalar and electromagnetic fields are coupled in the background solutions.
  Without such a coupling in the background solutions,
  one would simply obtain scalar or electromagnetic or gravitational modes
  by solving the respective scalar, electromagnetic or gravitational perturbation equations, 
  since the different types of perturbations would not be coupled to each other.
  However, when all fields are already present in the background solutions,
  the different types of perturbations couple.
  By taking the respective charges, Q_s and Q, to zero, the perturbations
  decouple again.
  Therefore, we employ a nomenclature that reveals this decoupling limit.
  
  
    i.
  Modes that are connected with purely scalar perturbations 
  are called scalar-led modes. 
  Typically they have dominant amplitude A_ϕ^±.
  
    ii.
  Modes that are connected with purely electromagnetic perturbations
  are called EM-led modes. 
  Typically they have dominant amplitude A_F^±.
  
    iii.
  Modes that are connected with purely gravitational perturbations 
  are called grav-led modes. 
  Typically they have dominant amplitude A_g^±.
  
  Scalar-led modes exist for l ≥ 0, EM-led modes for l ≥ 1
  and grav-led for l ≥ 2.
  
  In the following, we present our results for the quasinormal modes 
  for the cases l=0, 1 and 2, successively.
  
  
  
   §.§ Spectrum of l=0 quasinormal modes
  
  
  
  
  The l=0 perturbations are most interesting, since they 
  are the ones that discriminate between stable and unstable EMs black hole solutions.
  In particular, they feature an unstable mode on the cold black hole branch,
  whereas the bald and the hot black hole branches possess only stable modes,
  as our study shows.
  
  
  We exhibit the lowest scalar-led l=0 modes vs. the charge to mass ratio q
  in Fig. <ref>,
  showing the scaled real part of the frequency, ω_R/M, in Fig. <ref>(a)
  and the scaled imaginary part -ω_I/M in Fig. <ref>(b).
  A positive imaginary part signals instability, whereas a negative imaginary part
  yields the damping time of the mode.
  
  Following the color coding of Fig. <ref>,
  the bald RN branch is shown in black, the cold EMs branch in blue, and the hot EMs branch in red. 
  
  As seen in Fig. <ref>,
  the fundamental RN branch has only little dependence on q,
  deviating only slightly from the Schwarzschild mode all the way up to the 
  extremal black hole. The full RN branch is also free of unstable modes.
  
  
  The cold EMs branch features an unstable scalar-led mode
  throughout its domain of existence.
  The mode tend to vanish for the largest charge to mass ratio of the cold branch, q=1.
  
  The purely imaginary frequency grows as q decreases. At a certain value of q it reaches a maximum, from where it decreases again to the end point of the branch,
  where it reaches a zero mode.
  Here the bifurcation with the hot branch is encountered, and thus the minimal
  charge to mass ratio q_min of both EMs branches is reached.
  Continuity then requires that at q_min, also the hot EMs branch
  has a zero mode.
  
  Besides the unstable mode, the cold branch features a stable
  scalar-led mode, which is close to the scalar-led mode of the RN black hole.
  
  Along this branch, from q=1 to q_min the frequency ω_R/M
  first decreases slowly and then increases, becoming larger than the frequency
  of the RN branch, while the damping rate |ω_I/M| increases monotonically.
  
  The corresponding stable scalar-led modes of the hot EMs branch extend from q_min
  to q_max, where an extremal singular EMs solution is reached.
  The fundamental branch starts at the zero mode at q_min,
  from where both the frequency ω_R/M and the damping rate |ω_I/M|
  rise at first almost vertically. They continue to rise monotonically,
  reaching final values above those of the extremal RN branch.
  The first overtone branch starts at the stable scalar-led mode at the bifurcation point.
  Its frequency ω_R/M rises also monotonically, but its damping rate |ω_I/M|
  exhibits an overall but not monotonic decrease.
  
  The reason, we exhibit not only the fundamental stable mode for the hot EMs branch 
  but also the first overtone is to 
  demonstrate continuity of the modes at the bifurcation with the cold EMs branch.
  The zero mode of the cold EMs branch turns into the fundamental mode of the
  hot EMs branch, whereas the fundamental mode of the RN branch can be followed
  via the first stable mode of the cold EMs branch to the first stable overtone
  of the hot EMs branch. 
  
  
  Of course, all three classical branches feature sequences of (further) overtones,
  
  not studied here in detail. In the RN case, they are well known. 
  There, the rapidly damped modes possess several peculiar features. For instance,
  the higher modes of the non-extremal black holes have been observed to spiral towards the modes
  of the extremal black hole with increasing q
  <cit.>.
  
  
  
  
   §.§ Spectrum of l=1 quasinormal modes
  
  
  
  
  
  The l=1 modes consist of polar scalar-led modes and both axial and polar EM-led modes.
  We start the discussion with the scalar-led modes, shown in Fig. <ref>.
  Again, the lowest RN mode changes smoothly with increasing charge to mass ratio q.
  The frequency ω_R/M increases monotonically, while the damping rate |ω_I/M|
  remains almost constant, showing a slight decrease towards the extremal endpoint q=1.
  
  The lowest scalar-led l=1 mode of the cold EMs branch changes smoothly 
  from the q=1 point to the bifurcation point
  with the hot EMs branch at q_min, exhibiting a monotonic 
  decrease of the frequency ω_R/M and a monotonic increase of 
  damping rate |ω_I/M|.
  Along the hot EMs branch, the change of the eigenvalue with increasing q
  is no longer monotonic. As compared to the mode of the extremal RN solution,
  the frequency ω_R/M of the extremal EMs solution is lower
  and the damping rate |ω_I/M| is higher.
  
  In Fig. <ref> we exhibit the axial and polar EM-led l=1 modes.
  The RN black holes are known to exhibit isospectrality, i.e., the axial and polar 
  EM-led modes are degenerate. The EM-led RN modes have an analogous
  q-dependence to the scalar-led RN modes, but their absolute values differ
  somewhat. 
  
  Starting from the q=1 bifurcation point, the axial and polar modes
  of the cold EMs branch follow closely the RN modes at first.
  The frequencies ω_R/M of both axial and polar modes are slightly higher
  than the RN frequencies, but do not deviate strongly even at 
  the bifurcation point q_min with the hot EMs branch.
  The damping rates |ω_I/M| start to deviate from the RN
  damping rate earlier, showing an opposite behavior for the axial and polar modes.
  For the axial modes, the damping rates increase monotonically,
  whereas for the polar modes, a more sinusoidal pattern is seen.
  
  Along the hot EMs branch from q_min≤ q ≤ q_max,
  the frequencies  ω_R/M of both axial and polar EM-led modes
  rise monotonically, with the polar frequencies rising almost twice
  as much compared to the axial ones.
  The damping rates |ω_I/M| change in  a non-monotonic manner
  again, first rising from the bifurcation point, and then exhibiting some
  oscillation.
  At the extremal EMs solution, the damping rates  reach similar values,
  corresponding to about twice the extremal RN value.
  
  
  
  
  
   §.§ Spectrum of l=2 quasinormal modes
  
  
  
  
  
  
  
  We now turn to the l=2 modes, which consist of 
  polar scalar-led modes, axial and polar EM-led modes, and additionally
  axial and polar grav-led modes. 
  Thus, we now have five types of modes for generic EMs black holes.
  For the RN case, however, there is again isospectrality of the axial and polar EM-led modes,
  and there is also isospectrality of the axial and polar grav-led modes.
  Thus, basically only three types of modes are left,
  which all show a very similar pattern.
  It also resembles the pattern seen for l=0 and l=1.
  We will therefore now focus the discussion on the more interesting 
  and varied behavior of the modes of the EMs black holes. 
  
  As seen in Fig. <ref>,
  the scalar-led l=2 modes change monotonically on the cold EMs branch.
  The frequency ω_R/M decreases, and the damping rate |ω_I/M|
  increases toward the bifurcation point q_min.
  Along the hot EMs branch, the frequency ω_R/M rapidly reaches a minimum
  and then increases, while the damping rate |ω_I/M| exhibits a sinusoidal
  behavior. At the extremal EMs solution, both the frequency and the damping rate
  are higher than that at the extremal RN solution.
  
  The EM-led l=2 modes are exhibited in Fig <ref>.
  They show a pattern that is very similar to the pattern of the EM-led l=1 modes,
  although the numerical values of the frequencies and damping rates differ, of course.
  The grav-led l=2 modes are exhibited in Fig. <ref>.
  The frequencies ω_R/M on the axial and polar cold EMs branches
  follow very closely the RN branch almost up to 
  the bifurcation point q_min.
  This also holds for the damping rate |ω_I/M| on the axial cold EMs branch.
  Only the damping rate along the polar EMs branch starts to deviate from the
  RN one somewhat earlier.
  Along the hot EMs branch, the frequencies ω_R/M show opposite
  behavior with the frequency increasing monotonically on the axial branch,
  while decreasing monotonically along the polar branch.
  The damping rate |ω_I/M| also increases monotonically along
  the axial branch, whereas it exhibits a strong sinusoidal behavior
  along the polar branch.
  
  
  
  
   §.§ Isospectrality
  
  
  
  
  As addressed in the above discussion, for RN black holes, 
  isospectrality of axial and polar quasinormal modes holds. Here, axial and polar modes
  coincide for EM-led modes as well as for grav-led modes, for any allowed
  angular number l.
  Moreover, in the extremal case, the EM-led modes with angular number l
  agree closely with the grav-led modes with angular number l+1
  <cit.>.
  As seen above and once more highlighted in Fig. <ref>,
  this isospectrality survives to a large extent along the cold EMs branch,
  since the EMs modes follow the RN modes closely over a large range of their existence.
  In view of Fig. <ref>, this may, however, not be too surprising,
  since the cold EMs branch itself follows closely the RN branch
  almost up to the bifurcation point q_min.
  
  In contrast, the hot EMs branch reaches far beyond the RN branch
  and thus also far beyond the cold EMs branch.
  Its modes are therefore expected to show a generic behavior on their own.
  In particular, this branch features a significant background scalar field.
  The presence of a background scalar field, however, leads to 
  a coupling of all the different perturbations, scalar, vector and gravitational
  in the polar modes, whereas the axial modes remain free of
  scalar perturbations.
  Therefore it should not be surprising that isospectrality gets strongly
  broken along the hot EMs branch, as illustrated in Fig. <ref>.
  
  The breaking of isospectrality in the presence of a scalar field
  has been observed before, of course, notably in Einstein-Maxwell-dilaton
  black holes <cit.>, 
  and Einstein-scalar-Gauss-Bonnet black holes
  <cit.>.
  
  
  
  
  
  § CONCLUSION
  
  
  EMs theories represent an interesting simple setting to study 
  generic properties of spontaneous scalarization of black holes  and, more generically, the interplay between bald and hairy black holes.
  Although the presence of scalar hair
  
   is charge-induced,
  many similarities with the more realistic curvature-induced
  spontaneous scalarization and scalar hairy black holes have been observed,
  lending weight to EMs studies beyond an intrinsic theoretical interest.
  
  Here, we have investigated the linear mode stability of EMs
  black hole solutions with a quadratic coupling function,
  representing type IIB characteristics:
  the presence of scalarized black black holes,
  while the RN branch remains stable throughout.
  In particular, the system features besides the bald RN branch
  a cold EMs branch and a hot EMs branch.
  The cold EMs branch exists 
  in the interval q_min(α) ≤ q ≤1,
  and the hot EMs branch 
  in q_min(α) ≤ q ≤ q_max(α).
  
  The EMs black holes on the cold branch
  has an unstable 
  scalar-led mode, present for all q_min(α) ≤ q ≤1.
  At q_min(α), this instability becomes a zero mode that is shared
  with the hot branch. By continuity, the hot branch then
  exhibits a stable scalar-led mode, starting at q_min(α).
  
  In our mode analysis, we have calculated the lowest quasinormal modes
  of all types of perturbations, scalar-led modes, 
  axial and polar EM-led modes, and axial and polar grav-led modes
  for all three (bald, cold, hot) branches.
  
  For the stable modes on the cold EMs branch,
  we have found close similarity with the respective RN modes,
  almost up to the bifurcation point with the hot EMs branch q_min(α).
  Since the cold EMs branch itself follows closely the RN branch,
  this behavior could have been anticipated.
  In contrast, the modes of the hot EMs branch show 
  a wide variation and large deviations from the modes of RN branch.
  But again, the hot EMs branch reaches far beyond the RN branch
  (for sufficiently large α).
  
  RN black holes possess degenerate axial and polar modes
  in the EM-led and grav-led sectors. As expected, this isospectrality 
  gets broken in the presence of a non-trivial background scalar field,
  since scalar perturbations contribute in the polar case but not in the axial case.
  Not surprisingly, the breaking of isospectrality is very limited on the cold EMs branch,
  but becomes very strong on the hot EMs branch, away from the bifurcation point.
  
  The analysis here has focused on an illustrative value of the coupling α=200 
  and for ℓ=0,1,2, but the results concerning instability are generic for all values of α; 
  moreover, it is not to be expected that higher l modes introduce more instabilities.
  Thus our analysis has shown that there are no further unstable modes
  apart from the unstable scalar-led l=0 mode of the cold EMs branch.
  This implies that there are two linearly mode-stable branches in the system,
  the bald RN branch and the hot EMs branch.
  While both branches have large regions in parameter space,
  where their black holes are the only existing black holes,
  there is also an overlap region q_min(α) ≤ q ≤1.
  Here, both branches coexist and both are mode-stable.
  
  However, when these branches are considered from a
  thermodynamical point of view, their reduced horizon areas differ
  except at one critical point q_th(α), where the two curves cross,
  as seen in Fig. <ref>.
  This might suggest that the RN black holes represent the physically preferred state
  for 0 < q < q_th(α),
  whereas the hot EMs black holes represent the physically preferred state
  for q_th(α) < q ≤ q_max(α).
  Here, dynamical calculations of the EMs evolution equations might
  give further insight into the interesting question of which 
  type of black hole will represent the end state of collapse in the overlap region.
  
  
  
  
  
  
  
  
    
  Acknowledgments
  JLBS, SK and JK gratefully acknowledge support by the
  DFG Research Training Group 1620  Models of Gravity
  and the COST Action CA16104 GWverse.
  JLBS would like to acknowledge support from
  the DFG Project No. BL 1553. 
  AMP is supported by the FCT grant PD/BD/142842/2018. This work was supported by  the Center for Research and Development in Mathematics and Applications (CIDMA) through the Portuguese Foundation for Science and Technology (FCT - Fundação para a Ciência e a Tecnologia), references UIDB/04106/2020 and UIDP/04106/2020, by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6
  of the article 23, of the Decree-Law 57/2016, of August 29,
  changed by Law 57/2017, of July 19 and by the projects PTDC/FIS-OUT/28407/2017 and CERN/FIS-PAR/0027/2019. This work has further been supported by  the  European  Union's  Horizon  2020  research  and  innovation  (RISE) programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. We would like to acknowledge networking support by the COST Action GWverse CA16104.
  
  
  
  
  
  § PERTURBATION EQUATIONS
  
  
  In this appendix we present the explicit form of the perturbation equations
  for the axial and polar case, and also the first terms in the asymptotic expansion of the perturbation functions.
  
  
  
   §.§ Axial perturbations
  
  
  By substituting the Ansatz (<ref>)
  and (<ref>) into the field equations 
  (<ref>)-(<ref>), 
  we obtain a set of differential equations for the axial perturbations:
  
      ∂_r h_0 = - iω h_1
       -4f(ϕ_0)(∂_ra_0)W_2  +2/r h_0-i/2ω r^3g   h_1[ -( ∂_rϕ_0) ^2 g ^2r^2(r-2m) 
      -2rg^2l(l+1)
       +2rg( ∂_rg   )( r∂_rm  + m - r ) 
       -2r^2g( ∂_r^2 g   )(r-2m)   +4r^2f(ϕ_0)g     ( ∂_ra_0) ^2(r-2m) 
      
        +( ∂_rg   ) ^2r^2(r-2m)+4 ( ∂_rm   ) r g^2  +4g^2 (r-m) ]  ,
  
  
      ∂_r h_1 =
      -iω r h_0 /g  ( r-2m )  +h_12g   ( ∂_r m   ) r-r( ∂_rg   ) (r-2m) -2g  m/2rg  ( r-2m  )  ,
  
  
      ∂_r^2W_2  =
      (  4/gf(ϕ_0)(∂_r a_0)^2-ω^2r/g(r-2m)+l(l+1)/r(r-2m)) W_2  -[ ∂_r g/2g - ∂_rm/r-2m + ḟ(ϕ_0)/f(ϕ_0)∂_rϕ_0 + m/r(r-2m)] ∂_rW_2
        +i(∂_r a_0) h_1/2ω r^3 g^2[ 4g^2(r-m) + 4rg^2∂_rm - r^2g^2(r-2m)(∂_rϕ_0)^2 - 2rg^2l(l+1)
      
       
      +4gr^2f(ϕ_0)(r-2m)(∂_ra_0)^2 -2gr^2(r-2m)(∂_r^2g)+r^2(r-2m)(∂_rg)^2 + (2r^2g∂_rm+2rg(m-r))∂_rg]
      
       
      +h_0[(∂_rg/2g^2+∂_rm/m(r-2m)-ḟ(ϕ_0)/f(ϕ_0)∂_rϕ_0/g+3m-2r/gr(r-2m))∂_ra_0-∂_r^2a_0/g]  .
  
  
  This system of coupled differential equations consists of
  two first order differential equations for h_0 and h_1, plus a second order differential equation for W_2.
  
  A perturbation with an outgoing wave behavior satisfying this system of equations has to behave like
  
      h_0 = r e^iω R[A_g^+(-ω+O(r^-1))+A_F^+(Q/r^2+O(r^-3))]  ,  
      
      h_1 = rω e^iω R[A_g^+(1+O(r^-1))+A_F^+(-Q/ω r^2+O(r^-3))]  ,  
      
      W_2 = e^iω R[A_g^+(Q(l-1)(l+2)/4ω r^2+O(r^-3))+A_F^+(1+O(r^-1))]  .
  
  
  In addition, close to the horizon, a perturbation with an ingoing wave behavior has to satisfy
  
      h_0    =    e^-iω R[A_g^-(1+O(r-r_H)) + A_F^-(O(r-r_H))]  ,  
      
      h_1    =   ω/r-r_He^-iω R[A_g^-(r_H^3/2√(f(ϕ_H))/ω√(g_1(r_H^2f(ϕ_H)-Q^2))+O(r-r_H)) + A_F^-(O(r-r_H))]  , 
      
      W_2    =    e^-iω R[A_g^-(O(r-r_H))+A_F^-(1+O(r-r_H))]  .
  
  
  The asymptotic expansion is determined by two independent amplitudes, one related to
  the space-time perturbation A_g^± 
  and another related to the electromagnetic perturbation A_F^±.
  
  
  
   §.§ Polar case
  
  
  
  By substituting the Ansatz for the polar perturbations
  (<ref>)-(<ref>)
  into the field equations (<ref>)-(<ref>), we obtain the following set of equations
  
  
      L+N=0  ,
  
  
  
      ∂_rH_1  =4f(ϕ_0)( ∂_ra_0) (W_1-∂_rV_1)  +(∂_rδ)r(r-2m)+2(∂_rm   )r-2 m
      /r ( r-2 m  ) H_1 
      -2 iω r^2/r ( r-2 m  )(T+L)  ,
  
  
  
      ∂_r(N + T)  = r^3( ∂_r ϕ_0   ) ^2(r-2m)f(ϕ_0)-4Q^2/8	( r-2m   ) r^2f(ϕ_0) (L-N) 
      +r-3m/r(r-2m)N
      
      
      +r-m/r(r-2m)L
      - 2e^2δ/r(r-2m)(Q(iω V_1+a_1)-i/4ω r^2H_1)+1/2ϕ_1  (∂_r ϕ_0)  ,
  
  
  
      ∂_rN =-r/4( ∂_r ϕ_0   ) (∂_rϕ_1) 
      - 1/r-2m[ ( -( ∂_rδ)r(r-2m)
      - ( ∂_rm ) r-m  +r ) (∂_rT)   -e^2δf(ϕ_0)r^2(∂_ra_0   )(∂_ra_1)
      
      -i e^2δω rH_1-L-1/2(2e^2δf(ϕ_0)r^2(∂_ra_0)^2+l(l+1) ) N
      -1/2e^2δḟ(ϕ_0)(∂_ra_0)^2r^2ϕ_1
      
      +2r^3e^2δω^2-l(l+1)(r-2m)+2(r-2m)/2(r-2m)T  -i ( ∂_ra_0   ) e^2δf(ϕ_0)		ω r^2W_1]  ,
  
  
  
      ∂_rT  =i/8ω r^2( 4e^2δf(ϕ_0)r^2( ∂_r a_0 )^2+r( ∂_r ϕ_0 ) ^2(r-2m)+2l(l+1)-8(∂_rm) ) H_1 
      
      +L/r+1/4ϕ_1  (∂_r ϕ_0)-1/r( r-2m)[r(∂_rδ) (r-2m)+r((∂_rm) +1) -3m] T ,
  
  
  
      ∂_r^2T  =1/r(∂_rL)+ 1/4 ( ∂_r ϕ_0   ) (∂_rϕ_1)+1/4r(r-2m)[ 4( 		r( ∂_rm   ) +5 m		 -3 r ) (∂_rT) + 4 ( ∂_ra_0   )   e^			2δf(ϕ_0)  r^2(∂_ra_1)  
      
      +2r^2e^2δḟ(ϕ_0)(∂_ra_0   )^2ϕ_1 +( 4 			e^2δf(ϕ_0)r^2( ∂_ra_0   ) ^2+r ( ∂_r ϕ_0   ) ^2(r-2m)+2l(l+1)-4 (2(∂_rm)  		-1) ) L  
      
      -(r ( ∂_r		 ϕ_0   ) ^2(r-2m)  -8		 ∂_rm   ) N +2 ( l+2 )  ( l-1 ) T  +4i ( ∂_ra_0   )   e^2δf(ϕ_0)				 ω r^2W_1]  ,
  
  
  
      ∂_r^2(T+N)=
      iω r e^2δ/r-2m(∂_rH_1)
      -2-3r(∂_rδ)-4(∂_rm)+6m(∂_rδ)/r-2m(∂_rN) 
      + 1/2(∂_rϕ_0)(∂_rϕ_1)
      
      
      + r(r-2m)(∂_rδ)+2m-2r+2r(∂_rm)/r-2m(∂_rT)
      - 2re^2δf(ϕ_0)/r-2m(∂_ra_0)(∂_ra_1) 
      - iω e^2δ(m-r+r∂_rm)/(r-2m)^2H_1
      
      
      -2riω e^2δf(ϕ_0)(∂_ra_0)/r-2mW_1
      -r e^2δḟ(ϕ_0)(∂_ra_0)^2/r-2mϕ_1
      +e^2δω^2 r^2/(r-2m)^2(N-T)+2re^2δf(ϕ_0)(∂_ra_0)^2/r-2mT 
      
      
      +[
      2∂_r^2δ+2/r-2m∂_r^2m-1/2(∂_rϕ_0)^2-2(∂_rδ)^2+2(∂_rδ)(m+r-3r(∂_rm))/r(r-2m)](N+T)  ,
  
  
  
      ∂_rV_1  =ie^2δωr^3/l(l+1)(r-2m)[(∂_ra_1)  +( ∂_ra_0   )(L+N-2T+ḟ(ϕ_0)/f(ϕ_0)ϕ_1)] + [1-ω^2r^3e^2δ/l(l+1)(r-2m)]W_1  ,
  
  
  
      ∂_r^2V_1  = 
      -r^2 e^2δω/(r-2m)^2[ω V_1-ia_1] 
      +∂_rW_1
      
       
      -1/r(r-2m)[2(m-r∂_r m)+r(r-2m)(ḟ(ϕ_0)/f(ϕ_0)(∂_r ϕ_0)-(∂_r δ))  ] [∂_rV_1-W_1]  ,
  
  
  
      ∂_r^2a_1  =  -iω(∂_rW_1)  - (∂_ra_0) [∂_r(L+N-2T+ḟ(ϕ_0)/f(ϕ_0)ϕ_1)]  
       
      
       +l(l+1)/r ( r-2m ) (iω V_1+a_1)-(ḟ(ϕ_0)/f(ϕ_0)( ∂_r ϕ_0) +(∂_rδ)+2/r)[(∂_r a_1)+iω W_1]
      
      
      -[( ∂_r a_0   )  ( ∂_r ϕ_0   ) ḟ(ϕ_0)/f(ϕ_0)+ ( ∂_r δ) (∂_r a_0)  +(∂_r^2a_0)  +2/r(∂_ra_0)] [ḟ(ϕ_0)/f(ϕ_0) ϕ_1+2(L+N)]  ,
  
  
  
      ∂_r^2ϕ_1  =( ∂_r ϕ_0   ) [∂_r(N+2T-L)]
      -irω e^2δ∂_rϕ_0/r-2mH_1-4r e^2δḟ(ϕ_0)/r-2mN-4riω  e^2δḟ(ϕ_0)/r-2mW_1 
      
      -4r e^2δḟ(ϕ_0) (∂_ra_0)/r-2m (∂_ra_1)+2r(∂_rm)+2m-2r+r(r-2m)(∂_rδ)/r(r-2m)(∂_rϕ_1) 
      
      +[2∂_rϕ_0∂_rδ-2∂_r^2ϕ_0-4r e^2δḟ(ϕ_0)(∂_ra_0)^2/r-2m+4(∂_rϕ_0)(m-r+r∂_rm)/r(r-2m)]L 
      
      +[l(l+1)/r(r-2m)-r^2ω^2e^2δ/(r-2m)^2-(ḟ(ϕ_0)/f(ϕ_0))^22re^2δf(ϕ_0)(∂_ra_0)^2/r-2m]ϕ_1  .
  
  
  An useful redefinition of the electromagnetic perturbations is 
  
      F_0(r)    =    -iω W_1(r) - dW_1(r)/dr , 
      
      F_1(r)    =    -iω V_1(r) - a_1(r)  , 
      
      F_2(r)    =    -W_1(r) + dV_1(r)/dr ,
  
  which allows to simplify the system of equations in a way that is more convenient for numerical calculations. After some algebra, the minimal set of differential equations, or
  Master equations, can be written in vectorial form by defining a 6× 6 matrix (see equation (<ref>)). 
  
  Polar perturbations with an outgoing wave behavior at infinity satisfying the previous equations behave like
  
      H_1 = r e^iω R[A_g^+(-2+O(r^-1))+A_F^+(2iQ/ω r^2+O(r^-3))+A_ϕ^+(-i/4ω^2 r^4(ω.ḟ/f|_∞ Q^2-iQ_S)+O(r^-5))]  ,  
            T = e^iω R[A_g^+(1+O(r^-1))+A_F^+(-iQ/ω r^2+O(r^-3))+A_ϕ^+(iQ_S/4ω r^3+O(r^-4))]  ,  
            F_0 = 1/r^2 e^iω R[A_g^+(-2Q+O(r^-1))+A_F^+(-il(l+1)/ω+O(r^-1))+A_ϕ^+(.ḟ/f|_∞Q/r+O(r^-2))]  ,  
            F_1 = e^iω R[A_g^+(Q/r+O(r^-2))+A_F^+(1+O(r^-1))+A_ϕ^+(-.ḟ/f|_∞Q/4r^2+O(r^-3))]  ,  
            ϕ_1 = 1/re^iω R[A_g^+(2Q_SM/r-.ḟ/f|_∞2 Q^2/r+O(r^-1))+A_F^+( Q l(l+1)/ω^2r^2-i QQ_S/ω r^2+O(r^-3))+A_ϕ^+(1+O(r^-1))]  .
  
  Note that, since the background solutions we have considered satisfy ϕ(∞)=0, and the theory coupling is such that .ḟ/f|_∞=0, some of the terms in this expansion vanish.
  
  On the other hand, polar perturbations with an ingoing wave behavior at the horizon satisfy
  
      H_1 = 1/r-r_H e^-iω R[A_g^-(
      2ir_H^2(2ie^δ_Hf(ϕ_H)ω r_H^3+f(ϕ_H) r_H^2-Q^2)/(2r_H e^δ_Hω+il(l+1))(f(ϕ_H)r_H^2-Q^2) + O(r-r_H))
      . 
      .
      +A_F^-(	4e^δ_Hf(ϕ_H)Q r_H^3/(r_H^2f(ϕ_H)-Q^2)l(l+1)(r-r_H)+O((r-r_H)^2))
      . 
      .
      +A_ϕ^-(ḟ(ϕ_H)/f(ϕ_H)-4 r_H Q^2 r_H^3/(r_H^2f(ϕ_H)-Q^2)l(l+1)(r-r_H)+O((r-r_H)^2))]  ,  
      
      T = e^-iω R[A_g^-(
      1 + O(r-r_H))
      .                                                           
      .
      +A_F^-(-2e^δ_Hf(ϕ_H)Q r_H(2r_H e^δ_Hω+il(l+1))/(2r_H^3e^δ_Hf(ϕ_H)ω+if(ϕ_H)r_H^2-iQ^2)l(l+1)(r-r_H)+O((r-r_H)^2))
      . 
      .
      +A_ϕ^-(ḟ(ϕ_H)/f(ϕ_H)- Q^2(2r_H e^δ_Hω+il(l+1))/(2r_H^3e^δ_Hf(ϕ_H)ω+if(ϕ_H)r_H^2-iQ^2)l(l+1)(r-r_H)+O((r-r_H)^2))]  ,  
      
      F_0 = e^-iω R[A_g^-(O(r-rH))+A_F^-(1 + O(r-r_H))+A_ϕ^-(O(r-rH))]  ,  
      
      F_1 = e^-iω R[A_g^-(-2iω/f(ϕ_H)l(l+1)+O(r-r_H)))+A_F^-(-ie^δ_Hω r_H^2/l(l+1)+O(r-rH))
      . 
      .
      +A_ϕ^-(ḟ(ϕ_H)/f(ϕ_H)^2iω Q/l(l+1)+O(r-rH))]  ,  
      ϕ_1 = e^-iω R[A_g^-(O(r-rH))+A_F^-(O(r-rH))+A_ϕ^-(1 + O(r-r_H))]  .
  
  
  The asymptotic expansion of the polar perturbations is characterized by three undetermined amplitudes: A_g^± (space-time perturbation amplitude), A_F^± (electromagnetic perturbation amplitude)
  and A_ϕ^± (scalar perturbation amplitude).
  
  
  
  
  ieeetr
  
  
  