Quasinormal modes of hot, cold and bald Einstein-Maxwell-scalar black holes

Einstein-Maxwell-scalar models allow for different classes of black hole solutions, depending on the non-minimal coupling function $f(\phi)$ 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(\phi)=1+\alpha\phi^4$ [1]. Besides the bald Reissner-Nordstr\"om 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\"om 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\"om 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 [2]. 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 (φ) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].
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 [10].
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 [1]. 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 [1], 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 nontrivial 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 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, 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 whereḟ (φ) = df (φ)/dφ, and T φ µν and T EM µν are the scalar stress-energy tensor and the electromagnetic stress-energy tensor, respectively.
To construct static spherically symmetric EMs black holes, we employ the line element 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 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 (3)-(5), we obtain the following set of ODEs for the functions g, m, a 0 and φ 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 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 (2). First of all, the RN black hole, which is also a solution of the EMs equations, is given by The scalarized EMs solutions are obtained numerically [1]. Their asymptotic behavior yields their global charges 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) [34]. Close to the horizon r = r H , the global charges have the expansion 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 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 : 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. 1, where we exhibit their reduced area a H (left) and reduced temperature t H (right) versus their charge to mass ratio q [1]. 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 [1]. 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 a0 , and employ the Ansatz For the metric, we introduce the perturbations functions F t and F r and employ the Ansatz 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 (17)- (19) into the general set of field equations (3)-(5) and utilizing the set of equations for the background functions (10) leads to the Master equation for the spherical perturbations. This is a single Schrödinger-like ODE for the function Z = rφ 1 : which involves the tortoise coordinate R * , where and the potential U 0 , reading 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., for r → ∞. As the perturbation propagates toward the horizon, we have to impose an ingoing wave behavior, i.e., 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 : where Y lm are the standard spherical harmonics. The corresponding Ansatz for the metric introduces the perturbation functions h 0 and h 1 . Inserting this Ansatz into the field equations (3)-(5) 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 where Ψ A denotes the perturbation functions in the form 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 For the electromagnetic field, we employ the perturbation functions a 1 , W 1 and V 1 and the Ansatz Finally, for the metric, we introduce the perturbation functions N , H 1 , L and T , and employ the Ansatz Inserting this Ansatz into the field equations (3)-(5) 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, The vector now contains 6 components: 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., [35][36][37][38]), 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 (10), subject to the boundary conditions following from the expansions at infinity (13) and at the horizon (14). For this purpose, we employ the solver COLSYS [39], 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 [32,33,[40][41][42]. 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. 2, showing the scaled real part of the frequency, ω R /M , in Fig. 2(a) and the scaled imaginary part −ω I /M in Fig. 2(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. 1, 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. 2, 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 scalarled 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 [35,36].

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. 3. 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. 4 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. 5, 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 6. 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. 7. 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 [35,36]. As seen above and once more highlighted in Fig. 8, 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. 1, 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. 8.

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. 1. 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.

A.1 Axial perturbations
By substituting the Ansatz (26) and (25) into the field equations (3)-(5), we obtain a set of differential equations for the axial perturbations: 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 In addition, close to the horizon, a perturbation with an ingoing wave behavior has to satisfy 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 .
On the other hand, polar perturbations with an ingoing wave behavior at the horizon satisfy