Dynamically and thermodynamically stable black holes in Einstein-Maxwell-dilaton gravity

We consider Einstein-Maxwell-dilaton gravity with the non-minimal exponential coupling between the dilaton and the Maxwell field emerging from low energy heterotic string theory. The dilaton is endowed with a potential that originates from an electromagnetic Fayet-Iliopoulos (FI) term in $\mathcal{N}=2$ extended supergravity in four spacetime dimensions. For the case we are interested in, this potential introduces a single parameter $\alpha$. When $\alpha\rightarrow 0$, the static black holes (BHs) of the model are the Gibbons-Maeda-Garfinkle-Horowitz-Strominger (GMGHS) solutions. When $\alpha \rightarrow \infty$, the BHs become the standard Reissner-Nordstr\"om (RN) solutions of electrovacuum General Relativity. The BH solutions for finite non-zero $\alpha$ interpolate between these two families. In this case, the dilaton potential regularizes the extremal limit of the GMGHS solution yielding a set of zero temperature BHs with a near horizon $AdS_2\times S^2$ geometry. We show that, in the neighborhood of these extremal solutions, there is a subset of BHs that are dynamically and thermodynamically stable, all of which have charge to mass ratio larger than unity. By dynamical stability we mean that no growing quasi-normal modes are found; thus they are stable against linear perturbations (spherical and non-spherical). Moreover, non-linear numerical evolutions lend support to their non-linear stability. By thermodynamical stability we mean the BHs are stable both in the canonical and grand-canonical ensemble. In particular, both the specific heat at constant charge and the isothermal permittivity are positive. This is not possible for RN and GMGHS BHs. We discuss the different thermodynamical phases for the BHs in this model and comment on what may allow the existence of both dynamically and thermodynamically stable BHs.


Introduction
Stable configurations are the preferred configurations in (thermo)dynamics. For BHs, the classical concept of stability is that of dynamical (or mechanical) stability. Is the BH robust against small, mechanical perturbations? Perturbative stability of BHs is typically assessed by computing the spectrum of quasinormal modes (QNMs) of the isolated BH. The absence of growing modes establishes mode stability. For the only static vacuum BH of General Relativity, the Schwarzschild solution [1], its mode stability was established by the seminal works of Regge and Wheeler [2] and Zerilli [3]. In this sense, the Schwarzschild BH is a preferred configuration.
The advent of BH thermodynamics in the 1970s [4,5], as a consequence of a semi-classical treatment of gravity, yields a different angle on BH stability. Is a BH in thermodynamical equilibrium with its environment robust against small fluctuations of, say, energy, or another of its defining parameters? This question is, a priori, different, because the BH is interacting with an environment, or reservoir, rather than being isolated. Different reservoirs, or statistical ensembles, can be considered. For the case of a Schwarzschild BH, whose only defining parameter is its mass, the question simplifies. Fixing the reservoir's temperature, one simply asks how does the BH temperature responds to a small fluctuation of the BH energy. It turns out that the BH heats up/cools down when it loses/absorbs energy. That is, it has a negative specific heat. Thus, under a small energy exchange, Schwarzschild BHs run away from thermal equilibrium when placed in a reservoir at fixed temperature. They are (locally) thermodynamically unstable. Consequently, in this sense, these BHs are not preferred configurations.
Consider now the addition of electric charge. In electrovacuum General Relativity, the only static electrically charged BH solution with a connected horizon is the Reissner-Nordström (RN) BH [6]. Dynamically, it is still perturbatively stable, since its QNMs decay in time [7,8]. Thermodynamically, however, one can now consider different interactions. Firstly, consider the BH can exchange energy, but not electric charge, with the reservoir. Thus, the BH charge is fixed. Is the heat capacity (at constant charge) still negative? It must be for small charge, as the solution reduces to the Schwarzschild BH. For sufficiently large charge to mass ratio, q ≡ Q e /M > √ 3/2, however, the specific heat (at constant charge) becomes positive [9]. Thus, preventing any charge exchanges, RN BHs with sufficiently large q are (locally) thermodynamically stable, oscillating around the reservoir's temperature when small exchanges of energy occur. This is the canonical ensemble. Under these conditions, RN BHs with q > √ 3/2 are preferred. Considering no electric charge exchanges is, however, non-generic. In a generic situations such exchanges do occur. The reservoir is now not only a reservoir of energy, kept at constant temperature, but also of (charged) particles, kept at constant "chemical" potential, which in this case corresponds to the electrostatic potential. This is the grand-canonical ensemble. Now, the assessment of stability must also consider a possible run away mode triggered by the wrong evolution of the BH's chemical (electrostatic) potential. As the BH absorbs (releases) positive charge, if its chemical potential decreases (increases), this promotes more charge absorption (release), and hence a run-away mode. The response function monitoring this effect is the isothermal permittivity (at constant temperature). It so happens that for RN BHs it becomes negative precisely for q > √ 3/2, exactly when the heat capacity (at constant charge) becomes positive. This can be understood as follows. Fixing the temperature, an increase of charge implies an (even larger) increase of the BH mass, for q > √ 3/2. This implies the BH size increases sufficiently so that, despite the charge increase, the electrostatic potential decreases. The bottom line is that RN BHs are always locally unstable in the grand canonical ensemble. For small charges, exchange of energy at constant charge (or at constant electrostatic potential) promotes a run-away mode. For large charges, exchange of charge at constant temperature promotes the instability. Thus, under these conditions, RN BHs are not preferred configurations. Intriguingly, however, the isothermal permittivity diverges as extremality (q → 1) is approached, suggesting the RN family is on the verge of another thermodynamical phase.
The advent of supergravity and string theory naturally led to considering Einstein-Maxwell models with an extra scalar field, a dilaton, non-minimally coupled to the Maxwell field with a particular exponential coupling. In fact, such models naturally occur also in the context of Kaluza-Klein theories. Gibbons [10], subsequently also with Maeda [11], considered the charged BH solutions in these models. They are charged BHs that possess scalar (dilaton) "hair". The BH solutions of these Einstein-Maxwell-dilaton models were later reobtained by Garfinkle, Horowitz and Strominger in the context of string theory [12]. We shall refer to them as GMGHS solutions.
The GMGHS BHs, with the particular coupling emerging from string theory, are sometimes regarded as a sort of generalisation of the RN BH of electrovacuum, albeit they do not reduce to the latter, except in the uncharged limit. Moreover, GMGHS BHs introduce three qualitative new physical aspects. Firstly, the excited dilaton around the BH creates an effective medium with electric properties, which can be faced as, say, inducing a varying magnetic permeability. Secondly, in the standard electrovacuum-like description (the Einstein frame), the GMGHS BHs do not have a smooth extremal limit; they are singular in that limit. Thirdly, these BHs allow a charge to mass ratio greater than unity.
The GMGHS BHs can be made more RN-like, in particular concerning the second property of the previous paragraph, by augmenting the model with a particular dilaton potential that has been shown to emerge in N = 2 supergravity in four spacetime dimensions, extended with vector multiplets and deformed by a Fayet-Iliopoulos term [13]. In [14], Anabalon, Astefanesei and Mann obtained exact charged BH solutions in this model, that we shall refer to as AAM BHs. The potential introduces one single extra parameter, α. Then, as α → 0, the stringy GMGHS BHs are recovered. On the other hand, as α → ∞, the potential confines the scalar field to vanish but the electric charge may remain non-zero, yielding the RN electrovacuum family. In this sense, the AAM family of BHs interpolates between the RN and the GMGHS families. Now, it turns out that the dilatonic potential regularises the extremal limit, yielding a family of extremal BHs with a near horizon geometry of Robinson-Bertotti type, AdS 2 × S 2 , analogue to the extremal RN solution. Nonetheless, the AAM solution retains the GMGHS property that it still allows overcharged BHs with q = 1. The AAM family of solutions, therefore, presents itself as an arena to test the hypothesis that the q → 1 RN BHs are on the verge of another thermodynamical phase.
In this work we therefore investigate the dynamical and thermodynamical stability of the AAM BHs. It was recently pointed out [15] that some of these BHs are thermodynamically stable in both the canonical and grand-canonical ensemble, unlike RN BHs. In this paper we perform a thorough scanning of the domain of existence, precisely identifying the subset of thermodynamically stable AAM BHs. Our analysis makes clear that this only occurs in the overcharged regime, q > 1. Indeed, the AAM solutions have three thermodynamical phases: i) a Schwarzschild-like phase, which is unstable since the specific heat (at constant charge or at constant electrostatic potential) is negative; ii) a near-extremal RN-like phase, which is unstable since the isothermal permittivity is negative; and iii) a new stable phase, for which a necessary, but not sufficient, condition is that the BH is overcharged. In this sense the AAM BHs are an extension of RN BHs into the overcharged regime, which is made possible by the conjugation of the dilaton non-minimal coupling and potential.
Our work also clarifies that, like RN BHs, GMGHS BHs are never thermodynamically stable in both the canonical and grand-canonical ensemble. Moreover, we show that, like RN BHs, the AAM family is dynamically perturbatively stable, in the sense that no growing quasi-normal modes are found, building upon the recent work [16], see also [17]. The dynamical robustness of AAM solutions is furthermore confirmed by fully non-linear numerical simulations with the Einstein-Maxwell-dilaton system. Even starting with initial data far off from the equilibrium solutions, the evolutions relax to the latter. These simulations were performed using a similar code and setup to the ones reported in [18][19][20]. All this put together means that there is a subset of AAM solutions that are preferred, both dynamically and thermodynamically. To the best of our knowledge, this is the first such example for asymptotically flat BHs, without using artefacts such as box boundary conditions. This paper is organised as follows. In Section 2 we describe the model and the AAM family of solutions. In particular we consider the extremal limit showing the existence of AdS 2 × S 2 geometries, unlike in the GMGHS sub-family. In Section 3 we present a linear analysis of dynamical stability. We consider separately spherical, axial and polar perturbations. Then we discuss the spectrum of QNMs showing that, within our scanning, they always decay in time, providing strong evidence of mode stability. In Section 4 we consider a non-linear analysis of dynamical stability, showing that fully non-linear dynamical evolutions within the Einstein-Maxwell-dilaton model, starting from initial data which can be considered a highly perturbed AAM BH, the evolution converges to the latter. In Section 5 we discuss thermodynamical stability in both the canonical and grand-canonical ensemble, spelling out the conditions and identifying the region of the domain of existence of AAM BHs where thermodynamical stability holds. Finally, in Section 6 we summarise our results and provide a discussion of their significance.
2 Einstein-Maxwell-dilaton BHs with a supergravity potential

The model
The model under consideration is an Einstein-Maxwell-dilaton model endowed with a particular dilaton potential. It is described by the action: Here, R is the Ricci scalar, F µν ≡ ∂ µ A ν − ∂ ν A µ is the Maxwell field and φ is the dilaton field, which couples non-minimally to the Maxwell field with a coupling constant γ. In the following, we use units where 8πG = 1 = c. The dilaton is endowed with the potential discussed in [15], which has the form: where α is a dimensionful constant, with dimensions length −2 . This potential is plotted in Fig. 1. For small φ, this potential behaves V αφ 5 /30; thus it contains no mass term. We notice that V is invariant under the discrete symmetry α → −α, φ → −φ. Thus, taking α > 0 without any loss of generality, the requirement V > 0 imposes that the physical solutions necessarily have a positive φ.
The equations of motion that follow from (2.1) are where are the dilaton and electromagnetic energy-momentum tensors. The potential (2.2) was first considered in [14]. It was originally engineered to obtain exact solutions in Einstein-Maxwell-dilaton gravity. Subsequently, it was shown that the Einstein-dilaton sector of (2.1) is a consistent truncation of N = 2 supergravity in four spacetime dimensions, coupled to a vector multiplet and deformed by a Fayet-Iliopoulos term [13]. Augmenting the latter model by introducing the standard Maxwell term with the dilatonic coupling that emerges in low energy heterotic string theory [12] leads to the action (2.1). It is possible, but it has not been explicitly shown, that the full model (2.1) emerges from supergravity. That is an interesting open question. But both its F = 0 and α = 0 truncations emerge from supergravity models.

The solutions in Schwarzschild-like coordinates
The charged, spherically symmetric BH solutions of the model (2.1) are the standard RN BH for γ = 0 = α, the GMGHS solution [12] for α = 0 and arbitrary γ, and the AAM BHs [14] for γ = 1 and arbitrary α. The last case will be the focus of our work. Obviously, the AAM solutions reduce to the subset of the GMGHS solutions with γ = 1 when α = 0. But, as shown below, they turn out also to reduce to RN BHs when α → ∞.
As discussed in [15], these solutions form two branches distinguished by the sign of the parameter α in the dilaton potential (2.2). In what follows we shall take α > 0. Then, following [15], there is a family of BH solutions with strictly positive φ, thus only probing the strictly positive region of the dilaton potential. Employing Schwarzschild-like coordinates 1 , the geometry of these BH solutions takes the form where the Misner-Sharp mass function m(r) [21] and the redshift function σ(r) read 1 In [15] the solution was considered in a different coordinate system, which makes its properties less transparent. The relation between the radial x−coordinate in [15] and the r-coordinate herein is x = (1 + 2η 2 r 2 + 1 + 4η 2 r 2 )/(2η 2 r 2 ), with η = 1/Qs. Also, the parameter q in Ref. [15] should not be confused with q = Qe/M cf. (2.20) in this work. where while the matter fields take the form: A 0 is an arbitrary constant. Besides A 0 , this solution contains two arbitrary parameters Q e and Q s , which correspond to the electric charge and scalar charge respectively, as read e.g. from the 1/r terms in the asymptotic expansion of the scalar and gauge field. Observe, however, that Q e is a gauge charge associated to a gauge symmetry, whereas Q s is not. 2 The ADM mass is read off from the asymptotic value of m(r); it can be expressed as a function of Q e , Thus, the scalar charge is not independent from the mass and gauge charge. We conclude that the scalar hair of the solution is of secondary type -see [23] for a discussion of scalar hair for asymptotically flat BHs. For a certain parameter range the above solution describes BHs. Then, there is a BH horizon at r = r H > 0. The relation between the electric charge, 3 scalar charge and r H reads: For these BHs, the Hawking temperature T H and the event horizon area A H read Also, of relevance for the analysis below, working in a gauge with A t (r H ) = 0 we find the electrostatic potential at infinity, which corresponds to the chemical potential in the thermodynamical analysis, is (2.14) For any choice of the dilaton potential (and a vanishing dilaton at infinity), one can verify that the solutions satisfy the first law of BH thermodynamics in the form Also, one can prove that solutions satisfy the following Smarr-law 16) 2 We also emphasise that unlike [22], due to the existence of the potential, to obtain asymptotically flat solutions the asymptotic value of the scalar field should be kept fixed. 3 The electric charge can be computed by the covariant form of Gauss' law Since there are no sources to Maxwell's equations (2.3) outside the horizon, the value of the charge is independent of the choice of the surface S, with At ∼ A 0 − Qe/r asymptotically. The scalar charge Qs is computed from the asymptotics of the scalar field, φ ∼ Qs/r. where the space integral is taken over a domain bounded by the event horizon and the sphere at infinity 4 . Moreover, as usual in Einstein gravity, the entropy is given by the Bekenstein-Hawking formula The model possesses the scaling symmetry where λ > 0 is a constant. Under this scaling symmetry, all other quantities scale appropriately, e.g. M → λM , Q e → λQ e and Q s → λQ s . We frame the physical discussion using quantities which are invariant under this transformation. These are the following reduced quantities The domain of existence of these solutions in the (q, a H ) and (t H , a H )-planes is shown in Fig. 2, where the parameter α spans the whole positive real line. Since α is dimensionful, a dimensionless parameter is obtained as αM 2 . One can see that the solution (2.6)-(2.9) interpolate between the dilatonic GMGHS BHs, occurring for α = 0 and RN BHs, which are approached as αM 2 → ∞ in which case φ → 0. The latter limit can be shown as follows. First define a new constant Q as introducing in the solution (2.6)-(2.9), the limit α → ∞ yields The expression (2.16) holds for a generic model (2.1). For the AAM solution, the explicit form of M φ reads This corresponds to RN BH with charge Q and horizon radius r H . The metric and gauge functions can be checked to have the correct limit. In Fig. 3 we plot the profile functions defining the solution for a RN (α = ∞), a AAM (α = 100) and a GMGHS (α = 0) BH, all with q = 0.9. One observes, in particular, that the behaviour of the AAM solutions is "in between" the RN and the GMGHS solutions. Setting the electric charge to zero, the AAM solutions (2.6)-(2.9) trivialise. This does not, by itself, exclude the existence of Einstein-scalar field BH solutions of the model (2.1), without electric charge. However, one can easily prove that no such solutions exist. The proof is based on a Bekenstein-type argument [23] and uses the equation for the scalar field, which, for (2.6) takes the form Multiplying this equation by φ, a simple manipulation yields When integrating the above relation between r H and infinity, the left hand side vanishes for a regular configuration, since N (r H ) = 0 and φ(∞) = 0; the integrand of the right hand side, on the other hand, is strictly positive (since dV /dφ > 0). Thus, only φ ≡ 0 is possible. The presence of an electromagnetic field results in a strictly negative extra-term in the right hand side of (2.24), which allows circumventing this argument.

The extremal limit and near horizon geometry
Unlike the α = 0 case (the GMGHS solutions), the AAM BHs with α = 0 possess a regular extremal limit, with a nonzero event horizon area, mass, and electric charge. For GMGHS solutions, the effective potential [24] does not have a critical point at the horizon and so the extremal solution is a naked singularity. This can be also understood from the fact that the inner horizon is singular, the Kretschmann scalar is divergent there, and so in the extremal limit the singularity is pushed to the outer horizon. In the case of AAM solution, the effective potential has a new contribution coming from the dilaton potential, yielding a regular near horizon geometry of the extremal solution.
The extremal limit is found by imposing the constraint T H = 0, which results in an involved relation between r H and Q s (or, equivalently, r H and Q e , cf. (2.13)). Therefore, similarly to the RN case, the extremal AAM solutions are characterized by a single parameter, which is more conveniently taken as the event horizon radius, r H . A different route to study this limit is to impose that the general model (2.1) possesses an AdS 2 × S 2 geometry as a solution and obtain the horizon data by using the entropy function formalism [25][26][27]. 5 This configuration describes the near horizon geometry of an extremal BH, and has a Robinson-Bertotti-type line element [29,30], with The matter fields ansatz is Here, v 0 , v 1 ,ẽ, φ H are constants that determine the near geometry and the values of the electric field and dilaton at horizon. Under this ansatz, the equations of the model reduce to the following three algebraic The constants v 0 and v 1 determine the AdS 2 'radius' and S 2 radius, respectively. The radius of the horizon can be computed, as usual, from the relation v 1 = r 2 H and the constantẽ determines the electric charge Q e via the conservation of the flux, Thus, one finds that the parameters v 0 and v 1 , which enter the near horizon geometry, as well as the electric charge Q e , are fixed by the value of the scalar field at the horizon φ H , yielding a continuum of solutions with Equivalently, one can obtain all horizon data as a function of the physical electric charge, but since we can not solve analytically the last equation to obtain the value of the dilaton at the horizon, we prefer to work with φ H rather than with Q e . We can explicitly check that all the expressions (2.29) are finite at the horizon for the potential (2.2). Since the entropy, S, of the BHs with this near horizon geometry is given by the Bekenstein-Hawking formula, it can be computed as S = πv 1 . Consequently, (2.29) together with (2.2) imply a relatively simple expression for the entropy of extremal AAM solutions as a function of the scalar field at the horizon: 6 If α = 0 (the GMGHS limit) there is no regular AdS 2 × S 2 geometry. Following the corresponding analysis in the Section 2.2, to recover the RN limit we first define with Q a new constant. Then, the limit α → ∞ yields v 0 = v 1 = Q 2 and S = πQ 2 , which are the results for an extremal RN solution with Q e = Q.

Mode stability -linear analysis
Let us now consider the dynamical stability of the AAM BHs described in section 2. Following the analysis in [16], we first consider mode stability against linear perturbations of the metric and fields. We perform the analysis first for spherically symmetric perturbations and, subsequently for generic non-spherical perturbations.

Spherical perturbations
Following the standard method, we consider linear spherically symmetric perturbations of the AAM BHs (see for example [31]). This can be implemented using the following ansatz: Here N (r), f (r) = N (r)σ 2 (r), φ 0 (r) and a 0 (r) correspond to the unperturbed solutions (2.6)-(2.9); P 1 (r), S 1 (r), φ 1 (r), V 1 (r) are the perturbation functions all associated to a Fourier mode with frequency ω, and is an infinitesimal parameter. The frequency is in general a complex number, ω = ω R + iω I . The real part ω R is related with the oscillation frequency of the perturbation. If the imaginary part ω I is negative, then the perturbation is exponentially damped with damping time 1/ω I . Mode instabilities would occur if QNMs with ω I > 0 exist. A straightforward computation shows that both the metric and the matter fields perturbations are determined by φ 1 . As such, the study of the system reduces to a single equation for the scalar field perturbation. This equation can be written in the standard 1D Schrödinger form: where we have defined the 'tortoise' coordinate x, and the new function Ψ by The perturbation potential U ω has the unenlightening form Expressed in this 1D Schrödinger form, the diagnosis of an unstable mode solution of (3.34) would be ω 2 < 0 (or more explicitly, ω R = 0 and ω I > 0 with Ψ| r H = Ψ| ∞ = 0). Since the potential is regular in the entire x−range and it vanishes at the BH event horizon and at infinity, this would be a bound state. A standard result in quantum mechanics is that (3.34) will have no bound states if the potential U ω is everywhere greater than the lower of its two asymptotic values. Although the potential is not manifestly positive definite, scanning the space of solutions, this positivity is indeed satisfied. This means that ω 2 is positive and the AAM BHs are stable against spherical perturbations.

Non-spherical perturbations
After a decomposition using tensor spherical harmonics, the generic non-spherical perturbations can be differentiated in two decoupled channels, axial perturbations and polar perturbations, depending on how they transform under reflection of the angular coordinates [32][33][34][35].
The axial channel only perturbs the metric and the gauge field; not the scalar field [16]. The ansatz for this sort of perturbations introduces three perturbation functions, h 0 , h 1 and W 2 . It reads: where Y lm (θ, φ) are the usual spherical harmonics. The system of equations obtained from the linearised Einstein-Maxwell-dilaton equations consists of two first order differential equations for h 0 and h 1 (the metric perturbations) coupled to a second order differential equation for W 2 (the electro-magnetic perturbation). The polar channel, on the other hand, perturbs all the fields: the metric, the gauge field and the scalar field. The ansatz in this case introduces eight perturbation functions, H 0 (r), H 1 (r), L(r), T (r), a 1 (r), W 1 (r), V 1 (r) and φ 1 (r). Now it reads: It is convenient to define Introducing this ansatz on the field equations and performing a number of algebraic manipulations, it is possible to see that the polar perturbations are described by: two first order differential equations for H 1 and T (metric perturbations), another two first order differential equations for F 0 and F 1 (electro-magnetic perturbations), and a second order differential equation for φ 1 (scalar perturbation). These six functions determine the rest of functions (H 0 , L and F 2 ) via some algebraic equations.
As a summary, the minimal system of equations for the non-spherical perturbations can be written like [16] ∂ r Ψ j = M j Ψ j , (3.42) where j = {Axial , Polar } and we have the perturbation functions for each of these channels: We have separated with semicolons the space-time perturbations, electromagnetic perturbations and scalar perturbations (only present in the polar channel). The coefficients of the 4 × 4 matrix M Axial and the 6 × 6 matrix M Polar depend on the unperturbed metric and field functions, the l harmonic index and the mode complex frequency ω.

The quasinormal modes
The AAM BH are solutions of the model (2.1), with γ = 1 and arbitrary α, and as we have seen these solutions are known in closed form, cf. section 2. For the analysis of linear stability and QNMs spectrum, however, we have obtained numerically the BHs in the more general case of model (2.1) with both values of γ and α arbitrary. This allows us to better explore limiting cases. For example, varying the dilaton coupling γ from γ = 1 to γ = 0, provides a simple procedure to recover the electrovacuum model. In this way, we can continuously track the QNMs, connecting them to all known spectra: Schwarzschild, RN and GMGHS [36]. This is an excellent cross-check on the numerical calculation of QNMs. Following this methodology and the procedure described in [16] (see also [37]), we have computed the spectrum of QNMs of the AAM BHs. When the model reduces to electrovacuum γ = 0 = α, our results for the QNM of the RN BH reproduce the results in [38]. In the other well known limit when α = 0, our results reproduce the spectrum calculated in [36] for the GMGHS BHs. In these two limiting cases all the modes are stable.
After benchmarking the method with these two limiting cases, we have tackled the AAM BHs in [15], scanning for possible unstable modes, in particular for the thermodynamically stable solutions. In the illustrative case α = 1, r H = 2, we have scanned for unstable modes several thermodynamically stable solutions between Q e = 3.305929 (extremal), and Q e = 3.11784 (critical, separating thermodynamically stable from unstable solutions -cf. Section 5), as well as solutions close to the critical point (Q e = 3, 3.1). No unstable modes were found for the l = 1, 2 cases (note l = 0 corresponds to the previous spherically symmetric perturbations). It is unlikely that unstable modes may exist for larger values of l.
In addition, by slowly varying α and γ, we have tracked the modes that for γ = 0 = α correspond to the RN spectrum and for γ = 1, α = 0 correspond to the GMGHS spectrum. It is possible to see that all of these modes remain stable when getting to the AAM family with γ = 1 and (say) αr 2 H = 4. This is shown in Fig. 4, where the left (right) panels exhibits the scaled real (imaginary) part of ω versus q, for gravitational (top), electromagnetic (middle) and scalar (bottom) perturbations. All the curves in these plots are for l = 2 modes. With a solid thick grey line we show the RN modes. The GMGHS modes are shown with a dashed red curve (axial) and a dashed orange curve (polar). The AAM modes are shown with a dotted blue curve (axial) and a dotted cyan curve (polar).
In Fig. 4 we call grav-led modes to those that correspond to purely gravitational modes when the couplings are turn to zero (Schwarzschild). Typically, for small and intermediate values of q, the grav-led modes excite more strongly the gravitational perturbations; that is, the amplitude of the metric perturbation functions is larger than the amplitude of the electromagnetic and scalar perturbations. Similarly, EM-led modes and scalar-led modes correspond respectively to purely electromagnetic an scalar modes in the Schwarzschild limit [16].
In Fig. 4 we can observe how the GMGHS and AAM modes diverge from the RN modes as we increase q. In all cases, the AAM curves fit "in between" the RN and GMGHS curves. We can also appreciate how the spacetime and electromagnetic modes split into two (axial and polar) when the charge is increased, for  both the GMGHS and AAM cases. For the scalar modes, no such splitting occurs, since they only exist as polar modes. This is related with the breaking of isospectrality in these dilatonic BHs. It is well known that RN modes possess isospectrality, meaning axial and polar modes are equal. This isospectrality is typically broken by the dilaton, as already noted in [16,36] (see also [39,40]).
Moreover, the presence of the dilaton couples the scalar perturbations to the full polar equations. Thus, while the axial channel retains two families of modes -spacetime (Fig. 4 top) and electromagnetic modes (Fig. 4 middle), as for RN -, the polar channel acquires three families of modes -the previous two plus the scalar modes (Fig. 4 bottom). In the limit q = 0, all cases reduce to the Schwarzschild BH, and the modes converge to the Schwarzschild spectrum. Isospectrality is restored, and polar modes merge with axial modes of the same family. The dilatonic modes converge to the modes of a minimally coupled scalar field on the Schwarzschild background. Moreover, considering the two limits of the domain of existence shown in Fig. 2, in the limit α = 0 with γ = 1, the modes converge to the QNMs of the stringy GMGHS BHs, where isospectrality is also broken. In the limit α → ∞, the modes tend to converge to the RN modes, and again isospectrality is restored.
In this scanning, we have observed that all modes remain stable, although close to extremality there is a tendency to increase the damping time, corresponding to a smaller values of the imaginary part of ω. This trend, however, occurs already for RN BHs.
To conclude, the spectrum of QNMs of the AAM BHs is qualitatively very similar to the one studied recently for scalarised RN BHs in [16]: all QNMs are damped, although the spectrum is richer due to the broken isospectrality and the non-trivial scalar degree of freedom. All these features strongly indicate that the AAM BHs are mode stable for arbitrary values of α and q, in both the axial and polar channels, and for arbitrary l numbers.

Dynamical stability -non-linear analysis
The linear stability analysis, as discussed in Section 3, does not rule out that large perturbations can cause instabilities. One way to assess this possibility is by performing fully non linear numerical simulations, within the framework of numerical relativity, starting with a highly perturbed configuration. With this goal in mind, we have performed non-linear evolutions of the model (2.1) with γ = 1.
To test the stability of an AAM BH against large perturbations we consider two scenarios. Firstly, we have started with a RN BH with some small scalar field profile around it. We face this as a highly perturbed AAM BH. This initial data has the advantage of being readily accessible within the framework of numerical relativity. It is, however, constraint violating initial data, in the sense that it does not solve the constraint equations obtained from (2.1). Nonetheless, as it often happens with constraint violating data, the evolution converges to a true, stable, solution of the model, which in this case is an AAM BH. The latter, however, has q < 1, similarly to the initial RN configuration. Thus, although such evolutions confirm the stability of the AAM solutions, one cannot, in this way, probe the overcharged regime, which is the most interesting one if one is interested in the thermodynamically stable BHs. This issue can be solved in our second scenario, where we start with a GMGHS BH configuration as initial data, but in a model with α = 0. Again, this is constraint violating initial data. Unlike the RN case, however, there are GMGHS BHs with q > 1; thus, we can start with such overcharged configurations and the evolutions confirm an overcharged AAM configurations forms. In this way we show that even overcharged AAM solutions are stable in a non-linear sense. We remark, however, that the simulations we performed did not form AAM solutions precisely in the thermodynamically stable region. Reaching this region is numerically challenging within our approach, as it is very close to extremality. Nonetheless, the results herein establish that, in the sense we discussed, both undercharged and overcharge AAM solutions are non-linearly preferred. We expect this extends to the thermodynamically stable region.
Accordingly, in the first scenario we take as initial data a RN BH configuration with charge to mass ratio q = 0.2, with the following dilaton field initial Gaussian distribution as an illustrative example, we have taken A 0 = 3 × 10 −4 , r 0 = 10M and λ = √ 8. We have evolved this system for the coupling γ = 1 and taken α = 0.01.
The framework for this numerical evolutions is the 3+1 spacetime decomposition. For the metric, this split is given by ds 2 = −(α 2 0 + β r β r )dt 2 + 2β r dtdr + e 4χ a dr 2 + b r 2 dΩ 2 , (4.44) where the lapse α 0 , shift component β r , and the (spatial) metric functions, χ, a, b depend on t, r. A conformally flat metric with a = b = 1 is chosen together with a time symmetry condition, i.e. vanishing extrinsic curvature, K ij = 0. A description of the code to perform the evolutions and previous numerical studies of dynamical scalarisation of RN BHs can be found in [18,19,[41][42][43]. The evolution are performed in spherical coordinates under the assumption of spherical symmetry. The time integration uses a second-order Partially Implicit Runge-Kutta (PIRK) method developed by [44,45]. Let us briefly describe the evolution equations for the models (2.1) that are used in our numerical evolutions, under the 3+1 split and assuming spherical symmetry. For the electric field E r and an extra variable, Ψ, introduced to damp dynamically the constrains, they take the form where K is the trace of K ij , we have taken κ 1 = 1 and Π ≡ −n a ∇ a φ. The Klein-Gordon equation is given by: (4.46) Figure 6: Snapshots of the time evolution of the scalar field profile on the equatorial plane in the evolutions that do not impose spherical symmetry. One observes the initial RN configuration with a small dilaton perturbation evolving towards a dilatonic AAM BH.
The matter source terms for the scalar field, to be used in the Einstein equations, read and for the electric field The momentum density j em r vanishes because there is no magnetic field in spherical symmetry. Numerical evolutions are made under a spacetime discretisation. The logarithmic numerical grid extends from the origin to r = 1500M and uses a maximum resolution of ∆r = 0.0125M .
Let us now describe the results obtained within this setup starting with the first sort of initial data described above. In Fig. 5 (left panel) we exhibit the radial profile of the scalar field at late times (t = 600), and we compare it with the corresponding analytic profile of the AAM solution described in Section 2 with the corresponding value of q and r H . In the right panels we show the time evolution of the BH charge Q e and the amplitude of the scalar field extracted at radius r 0 = 11.09. These quantities clearly stabilise reaching an equilibrium configuration, which matches the AAM solution. At the end of the evolution, the final charge is Q e = 0.196. To assess the dynamical stability of the AAM without imposing spherical symmetry, we carried out numerical simulations in a 3D cartesian grid using the Einstein Toolkit [46,47]. To perform the evolutions we have used a numerical grid with 11 refinement levels with where the first set of numbers indicates the spatial domain of each level and the second set indicates the resolution. The evolution is identical to the spherically symmetric case. In Fig. 6 we plot snapshots of the evolution of the scalar field. Similar diagnostics as for the 1+1 evolutions attest the convergence of the initial data towards an AAM BH. The second sort of initial data is dealt with in a very similar way. We now start with a GMGHS solution with q = 1.1. We consider three models with α = 0.1, 1, 10. The evolution changes the value of the scalar field and increases slightly the value of the charge to mass ratio of the BH. We obtain, respectively, q = 1.1001, 1.1001, 1.101. The radial profile of the scalar field can be matched with a AAM BH with, respectively, Q s = 1.10, 1.03, 0.8. The profiles of the scalar field at t = 1000 and the evolution (left panel) of the evolution of the scalar field at some extraction radius (right panel) can be seen in Fig. 7.
These evolutions agree with our expectations. Starting from RN, but forcing the system to evolve to an AAM solutions, due to the dilaton coupling, the BH grows a scalar charge. Starting from GMGHS, but forcing the system to evolve to an AAM BH due to the non-trivial potential, the BH loses some scalar charge. The first/second case establishes the dynamical formation of an undercharged/overcharged AAM BH.

Thermodynamical stability
As pointed out in [15], a unique property (within asymptotically flat spacetime BHs) of AAM BHs is that they possess a subset which is locally thermodymically stable. We shall now examine precisely when this occurs in the domain of existence.
Thermodynamical stability can be local or global. Moreover, different thermodynamic ensembles represent physically different situations and may not lead, in general, to the same conclusions regarding the thermodynamical stability. Mathematically, (local) thermodynamical stability is equated with the sub-additivity of the entropy function. In the canonical ensemble, this is equivalent to the positivity of the specific heat at constant electric charge In the grand canonical ensemble, one requires instead the positivity of the specific heat at constant electric potential, and the positivity of the isothermal permittivity In fact, if C Q and T are positive, the identity implies C Φ > 0. Thus, local thermodynamical stability follows from C Q 0 and T 0. For the Schwarzschild BH there are no electric variables and the specific heat is negative whereas the electric permittivity is (5.55) Thus, RN BHs exhibit two phases, depending on the sign of r 2 H − Q 2 e , which vanishes for q = √ 3/2. For q < √ 3/2, C Q < 0 and T > 0. This is the Schwarzschild-like phase. For q > √ 3/2, C Q > 0 and T < 0. This is the near extremal RN-like phase. RN BHs are stable in the canonical ensemble in the near extremal RN-like phase. But they are always unstable in the grand canonical ensemble. Notice, however, that C Q is vanishing at extremality, wherein T is diverging. Thus, the RN family seems to be approaching a new phase at extremality, wherein BHs cease to be possible.
For the GMGHS family with γ = 1, parameterised by M, r H , the response functions take the simple form The negativity of the specific heats implies local thermodynamical instability.
In the AAM case, the expression of C Q , C Φ , and T are long and not enlightening [15], thus we do not include them here. All the response functions can, nonetheless, be presented as functions of r H and Q e . A study of these quantities show the existence of a region in the domain of existence where the thermodynamic stability is satisfied in both canonical and grand-canonical ensembles. This region is bounded by the set of extremal solutions T H = 0 and a set of critical configurations where C Φ and T diverge (and change sign afterwards) while C Q remains positive and finite. The critical configurations are found numerically, the relations providing a good fit for the corresponding curves in Fig. 8. Let us also remark that all thermodynamically stable AAM BHs have 1 < q < √ 2, while √ 2/2 < Φ < 1. Moreover, the set of locally thermodynamically stable solutions are also globally stable. That is, in a grand canonical ensemble, (i.e. for the same T H , Φ) they minimise the Gibbs free energy (5.57) The generic picture is summarised in Fig. 9 (left panel). For any value of 1/ √ 2 < Φ < 1, the G(T H ) curve consists in two parts. The branch minimising the free energy G starts at T H = 0 and ends at some maximal T H (Φ), consisting in configurations which are locally stable. The situation changes for Φ < 1/ √ 2, in which case, similar to the RN or GMGHS cases, one single branch of locally unstable solutions is found 7 , with G > 0.
When considering instead a canonical ensemble, one finds the existence of two branches, for any (fixed) value of the electric charge Q e . The solutions minimising the Helmholtz free energy F = M − T H S , (5.59) are located on the lower branch, which starts with the T H = 0 extremal BHs. These configurations have also a positive specific heat, C Q > 0, while a part of them are also stable in a grand canonical ensemble, C Φ > 0.
To summarise, we conclude that a set of AAM solutions which are overcharged, q > 1, and with a large enough chemical potential, Φ > 1/ √ 2, are thermodynamically stable, both locally and globally.

Conclusions
In this paper we have discussed the set of BH solutions found in [14] within Einstein-Maxwell-dilaton theory with a certain dilaton potential. Such AAM BHs can be thought of as a family of solutions that interpolates between the standard RN electrovacuum BHs and the GMGHS solutions of low energy heterotic string theory in four dimensions, retaining some of the features of both these limits. In particular, these BHs have a regular extremal limit and no electric charge outside the horizon, analogously to the RN BH; on the other hand, 7 One finds  they can be overcharged, i.e. to have a charge to mass ratio exceeding unity, as GMGHS. The combination of these properties allows in particular the exceptional feature in BH physics of exhibiting thermodynamical stability in both the canonical and grand canonical ensemble. In this sense, the overcharged AAM BHs can be faced as an extension of the RN family beyond extremality. Although there is a subset of AAM BHs that have both dynamical and thermodynamical stability, they are still afflicted by the decay induced by quantum effects, that is, Hawking radiation, except for the extremal solutions. The extremal AAM BHs are then stable also against Hawking evaporation. One cannot exclude, however, if these solutions are not supersymmetric, that non-perturbative effects may destabilise them.
We would like also to point that similar results in the thermodynamical analysis have been found for a second model discussed in [15], still described by (2.1), but with γ = √ 3 and a different V (φ). Clearly, some of the analysis herein could be repeated for that model.
Finally, despite the ingredients we have identified, we cannot pinpoint exactly the mechanism behind the existence of these dynamically and thermodynamically stable BHs. In particular the properties of the potential that permit them to exist. It is well known that AdS BHs can become thermodynamically stable. It is then tempting to think the dilaton potential is inducing AdS-like features. There is, however, an important difference between these two cases. In the AdS case, large BHs are thermodynamically stable; in the case analysed herein, the stable BHs are the smallest ones. In this respect it is worth remarking that the potential (2.2) is decaying towards the spatial infinity in the AAM BH solutions. Thus, even if induces a box-like effect, such effect may be more effective for small BHs.