Spontaneous scalarisation of charged black holes

Extended scalar-tensor-Gauss-Bonnet (eSTGB) gravity has been recently argued to exhibit spontaneous scalarisation of vacuum black holes (BHs). A similar phenomenon can be expected in a larger class of models, which includes e.g. Einstein-Maxwell-scalar (EMS) models, where spontaneous scalarisation of electrovacuum BHs should occur. EMS models have no higher curvature corrections, a technical simplification over eSTGB models that allows us to investigate, fully non-linearly, BH scalarisation in two novel directions. Firstly, numerical simulations in spherical symmetry show, dynamically, that Reissner-Nordstr\"om (RN) BHs evolve into a perturbatively stable scalarised BH. Secondly, we compute the non-spherical sector of static scalarised BH solutions bifurcating from the RN BH trunk. Scalarised BHs form an infinite (countable) number of branches, and possess a large freedom in their multipole structure. Unlike the case of electrovacuum, the EMS model admits static, asymptotically flat, regular on and outside the horizon BHs without spherical symmetry and even without any spatial isometries, which are thermodynamically preferred over the electrovacuum state. We speculate on a possible dynamical role of these non-spherical scalarised BHs.

Introduction. The newborn field of gravitational wave astronomy [1][2][3][4][5][6] will test the nature of astrophysical black holes (BHs) in an unprecedented way. It is therefore of the utmost importance to theoretically scrutinise all physically reasonable alternatives to the Kerr BH [7] hypothesis, a paradigm supported by the uniqueness theorems for vacuum General Relativity (GR) [8][9][10].
A reasonable alternative BH model requires, conservatively, being a solution of a consistent physical theory and having robust dynamical properties. The latter include a formation mechanism and sufficient stability. A novel class of such models has recently emerged in extended scalar-tensor-Gauss-Bonnet (eSTGB) gravity [11][12][13]. In eSTGB the Schwarzschild BH may become spontaneously scalarised, since linear theory reveals it becomes unstable against scalar perturbations for sufficiently small BHs. Moreover, BHs with scalar hair exist in the model which are thermodynamically preferred over the vacuum solution [11] and stable against spherical perturbations [14]. This phenomenon is akin to the spontaneous scalarisation found long ago for neutron stars in the context of scalar-tensor theories [15], but with the scalarisation induced by strong spacetime curvature, rather than matter (see [16] for a discussion of matter induced BH scalarisation).
The studies of BH scalarisation in [11][12][13][14] were both non-dynamical and, for the non-linear scalarised solutions, restricted to the spherically symmetric sector. However, the instability of the vacuum solutions observed in linear theory contains non-spherical modes. This raises the issues of what other scalarised, static, BH solutions are admitted in eSTGB and what is the dynamical evolution and endpoint of the instability of vacuum BHs.
A complete understanding of these evolutions requires non-linear numerical simulations; these are challenging, particularly if no symmetries are imposed, in view of the higher curvature terms in eSTGB gravity. In this letter we tackle these questions by observing that the eSTGB model belongs to a larger universality class, that in particular includes the technically simpler Einstein-Maxwellscalar (EMS) models [44]. In an illustrative EMS model, we perform, in spherical symmetry, non-linear numerical simulations exhibiting, dynamically, the scalarisation of the Reissner-Nordström (RN) BH. The process indeed forms a perturbatively stable scalarised BH. In the same EMS model, we study the non-spherical, static, scalarised solutions, and show they are thermodynamically preferred over the (electrovacuum) spherical BH. Since in cousin models to EMS, non-linear numerical evolutions of binary BHs have been performed [17], even non-spherical, non-linear evolutions of the scalarisation instability in EMS models are within reach.
Besides exhibiting, dynamically, the scalarisation instabilty in the spherical sector, our investigation of the static non-spherical sector constructs the first examples of static, asymptotically flat, regular on and outside the event horizon BHs without spatial isometries. This is a maximal violation of Israel's uniqueness theorem [18], stating that static BHs in electrovacuum must be spherical. The non-minimal couplings in EMS models circumvent in a radical way the uniqueness and simplicity of (electro)vacuum BHs [10], which holds even if a minimal coupling to a real scalar field is included [19,20]. It endows BHs with multipole freedom. Universal conditions and a toy model. Consider the following scalar field action on a non-trivial background (e.g. a curved spacetime and/or electromagnetic field): (1) f (φ) is the coupling function and I is a source term which generically depends on matter field(s), ψ, and the metric tensor, g µν . The scalar field equation of motion is φ = f I/4. f (φ) must allow the existence of a scalarfree solution, with φ = 0. This requires f (0) = 0. Spontaneous scalarisation occurs if the scalar-free solution is unstable against scalar perturbations δφ. These obey ( − µ 2 eff )δφ = 0, where µ 2 eff ≡ f (0)I(background)/4. If µ 2 eff < 0, a tachyonic instability sets in, driving the system away from the scalar-free solution.
Consider, as a field theory illustration of the instability, where F = dA is the Maxwell tensor, in Minkowski spacetime. The scalar-free configuration is the Coulomb solution (in standard spherical coordinates): Consider also (in appropriate units) These choices obey f (0) = 0 and µ 2 eff < 0. Thus (3) is unstable against scalarisation. This model admits a simple closed form spherically symmetric scalarised solution: where the integration constant obeys |ζ| < 1, to avoid divergences in the coupling. The total energy of the configurations is E = 4π ∞ r0 ρdr, where ρ = −T t t is the energy density. We consider a conducting sphere at r = r 0 and solutions (3)-(4) only for r r 0 , to regularise the total energy. The energies of these exterior configurations are, respectively, E (φ=0) and E (φ =0) . Then, r0 . The scalarised solution is energetically favoured in a set of bands, defined as Q/r 0 ∈ π]n+1/2, n+1[. The integer n ∈ N 0 , labelling the bands also counts, via (4), the number of nodes exterior to r 0 of the scalar field profile. Thus, within these bands, the instability (likely) evolves dynamically (3) into (4). Spontaneous scalarisation of BHs. The toy model shows that spontaneous scalarisation (i) is not exclusive of gravitational models and (ii) can be supported by an electromagnetic non-minimal coupling. We shall now focus on the case of BHs and consider the gravitational model (G = 1 = c): where R is the Ricci scalar. On a spherical, scalar-free BH solution, with a generic line element performing a (real) spherical harmonics decomposition of the scalar field φ(r, θ, ϕ) = m Y m (θ, ϕ)U (r), the scalar field equation becomes e δ r 2 d dr Eq. (7) is an eigenvalue problem: for a given , requiring an asymptotically vanishing, smooth scalar field selects a discrete set of BHs. These are the bifurcation points of the scalar-free solution. The (test) scalar field profiles they support -hereafter scalar clouds -are distinguished by the node number of U (r), n ∈ N 0 , besides ( , m).
In [11][12][13] the model (5) was studied with I = L GB , where L GB is the Gauss-Bonnet invariant, and various different coupling functions, satisfying f (0) = 0 and µ 2 eff < 0 (see also [16,21]). In this case the spontaneous scalarisation of the Schwarzschild solution is induced by this higher curvature correction, which is non-trivial for those BHs. Here, we shall take the source (2). Then, scalarisation of the Reissner-Nordström (RN) BH is induced without the need of higher curvature corrections. Scalarisation in EMS models. We are interested in models (5) with (1) and (2) which: (i) admit the scalarfree RN solution. This rules out the usual Einstein-Maxwell-dilaton model [22], where f (φ) = e −αφ ; (ii) approach the standard Einstein-Maxwell system in the far field, i.e. φ → 0 and f (φ) → 1 as r → ∞.
From the scalar equation of motion, it is possible to derive two Bekenstein type identities [23] that set the following constraints on f (for a purely electric field F 2 < 0): f ,φφ > 0 and φf ,φ > 0 for some range of the radial coordinate r. A simple potential compatible with the above requirements, that we shall use hereafter, is The coupling α is a dimensionless constant and the conditions on f imply α < 0 for F 2 < 0. The RN (scalar-free) solution is given by (3) and (6), with δ = 0, N (r) = 1 − 2M/r + Q 2 /r 2 . The scalar perturbations on this background are given by (7) with µ 2 eff = αQ 2 /r 4 < 0, hence exhibiting the instability. For = 0, one finds an exact test field solution where u ≡ ( √ 4α + 1 − 1)/2, r H ≡ M + M 2 − Q 2 and P u is a Legendre function. For generic parameters (α, Q, r H ), finding the = 0 bifurcation points from RN reduces to studying the zeros of this function as r → ∞. Examples are shown in Fig. 4 below. Bifurcation requires α below a maximal value, α max . For fundamental modes (n = 0) we obtain α max = −1/4 ( = 0), α max −2.784 ( = 1) and α max −7.087 ( = 2). Spherical sector (no nodes): domain of existence. The spherical scalarised BHs bifurcate from the RN solution for any α < −1/4. They are the non-linear realisations of the (n, , m) = (0, 0, 0) clouds. These non-linear solutions were obtained by using the ansatz (6), φ = φ(r), A = V (r)dt, and standard numerical techniques to solve coupled non-linear ODEs -see Appendix A for details. Varying α we have obtained the domain of existence of spherical scalarised BHs (with n = 0) shown in Fig. 1.
For each α, a branch of scalarised solutions bifurcates from a RN BH with a particular charge to mass ratio q ≡ Q/M . As α varies, these RN BHs define an existence line. Each constant α-branch ends at a critical, (likely) singular, configuration: the numerics indicate the Kretschmann scalar and the horizon temperature diverge, the horizon area vanishes, whereas the mass and the scalar "charge" (defined as Q s ≡ − lim r→∞ r 2 dφ/dr) remain finite. Along the α=constant branch, q increases beyond unity. Thus, scalarised BHs can be overcharged.
This potential vanishes both at the BH event horizon and at infinity, being regular everywhere in between. It follows that the Schrödinger eq. will have no bound states if the potential is everywhere positive. For all n = 0 solutions analysed, this positivity is indeed satisfied. The absence of such bound states guarantees the stability of the scalarised BHs against this class of perturbations. A final piece of dynamical evidence that the scalarisation of the RN BH leads to the solutions in Fig. 1 comes from performing fully non-linear numerical evolutionssee [24] and Appendix B for details. We start with a RN BH and a small Gaussian perturbation for φ and monitor the scalar field energy, E φ = ∞ r H n α n β T SF αβ dV , where T SF is the φ stress-energy tensor and n the 4-velocity of the Eulerian observer in the 3+1 spacetime decomposition [25]. Fig. 3 shows the time evolution of E φ for an initial RN BH with q = 0.2 and different couplings. One observes an initial growth of E φ followed by a saturation and equilibrium. For the largest |α| saturation is faster and it is preceded by an overshooting of the equilibrium value. The inset shows the value of the scalar field at the horizon (areal) radius, both after equilibrium is reached in each simulation (crosses) and from the scalarised BHs with the same q = 0.2 computed as static solutions (line). The values agree within 1%. Thus, the end point of the evolutions are the perturbatively stable scalarised BHs with the same q as the initial unstable RN BH. A larger set of numerical evolutions, exploring the (q, α) space, confirms scalarisation is generic and leads to the BHs in Fig. 1, but for larger q values of the initial RN BH, the evolution does not preserve (even approximately) q. A discussion of these evolutions will be presented elsewhere. Non-spherical sector: multipole freedom. Generic (n, , m)-scalar clouds also have non-linear realisations, yielding branches of non-spherical scalarised BHs bifurcating from the RN trunk. We have constructed such solutions using the Einstein-De Turck approach [26,27] -see [28] for a review and Appendix C for details.
All configurations constructed are regular on and outside a topologically (but not geometrically) spherical horizon and asymptotically flat. Solutions bifurcating from a scalar cloud with m = 0 (m = 0) are (are not) axially symmetric. The latter, in fact, have no spatial isometries. They provide the first explicit example of static, asymptotically flat BHs without any continuous (spatial) symmetries (see [29][30][31] for related work).
In Fig. 4 we exhibit results concerning solutions with (n, , m) = (0, 1, 0); (0, 1, 1); (0, 2, 0); (0, 2, 1); (0, 2, 2) and (n, 0, 0), for an illustrative value of the coupling α. The latter are for n = 1, 2, 3, corresponding to excited spherical scalarised BHs. The top panel shows that as either or n increases, the bifurcation point moves towards extremality of the scalar-free RN solution. This bifurcation point does not depend on m (Fig. 4, bottom panel) as anticipated from the linear theory analysis. Also, the relative location of the bifurcation point for the first few values of , n can be seen in the inset of the bottom panel. We remark that when continuing scalar clouds with , m = 0 into the non-linear regime, the corresponding scalar field does not remain a pure , m mode; nonlinearities excite all modes with the same m and parity.
All scalarised solutions are entropically preferred over a comparable RN solution (with the same total charge and mass). Within the scalarised solutions, the preferred one is the fundamental state with = 0 = m, i.e. spherically symmetric. Also, for all studied cases, for constant m, the entropy is maximised by the branch of solutions emerging from the = m zero mode. These entropic arguments are, however, inconclusive for dynamical considerations. General dynamics: the charged drop analogy. Consider a generic (i.e. non-spherical) perturbation of the RN BH with sufficiently large q. The existence of the non-spherical branches of scalarised BHs allows nonspherical perturbations to grow rather than damp (as they would in electrovacuum). A suggestive parallelism with charged drops in fluid dynamics can be drawn.
A self-gravitating, uncharged, isolated liquid mass is spherical and stable against small perturbations [32]. But if the liquid is electrically charged, conducting and surrounded by an insulator (e.g. a gas), the competition between the cohesive tension and the electric repulsion makes the spherical drop (which remains a solution) unstable, for charge beyond the Rayleigh limit [33]. New branches of non-spherical solutions of the fluidelectrostatic equations, associated to each spherical harmonic, emerge, bifurcating from the trunk family of spherical fluid balls. Such structure is analogous to that shown in Fig. 4 -see, e.g. Fig. 2 in [34]. Both experiments [35] and simulations [36] show that the unstable charged drop evolves towards a non-spherical shape.
BHs in the EMS model (unlike charged drops) have a = 0 mode instability which may always dominate and lead to a spherical scalarised BHs. If, however, there are initial conditions such that = 0 modes dominate, a dynamical symmetry breaking may occur. Fully nonlinear numerical evolutions in the EMS model with current technology can probe this possibility.
Consider the metric ansatz (6), with N (r) = 1 − 2m(r)/r, where m(r) is the Misner-Sharp mass function [37]. The gauge connection ansatz is purely electric, A = V (r)dt and the scalar field is a function of r only. This ansatz yields the following effective Lagrangian: the corresponding equations of motion being The electric potential can be eliminated from the above equations noticing the existence of a first integral, where Q is the electric charge. We assume the existence of a horizon located at r = r H > 0. In its vicinity one finds the following approximate solution to the field equations where This approximate solution contains two essential parameters, φ 0 and δ 0 , which are found matching the above expansion with the following asymptotics of the solutions in the far field: Apart from M and Q, the essential parameters here are Φ (the electrostatic potential) and Q s (the scalar "charge"). The solutions satisfy the virial identity: Since the term in square brackets, in the first line is always non-negative, this shows that a nontrivial scalar field can only be supported if Q = 0. The data at infinity is specified by the ADM mass M , electric charge Q, electrostatic potential Φ and Q s . The horizon data is the Hawking temperature and horizon We have noticed that the solutions satisfy the following interesting relation Remarkably, one can shown that (A8) holds for any choice of the coupling function f (φ) (with φ → Q s /r as r → ∞).
One observes that the solutions do not change when considering the scaling r → λr, Q → λQ, while α is fixed. q, a H are invariant under these transformations.

Appendix B: Numerical evolutions
To address numerical evolutions in the EMS system, we have modified the code described in [24,38] and adapted the evolution equations in [17] to the coupling described above. The 3+1 metric split reads ds 2 = −(α 2 0 + β r β r )dt 2 + 2β r dtdr + e 4χ a dr 2 + b r 2 dΩ 2 , where the lapse α 0 , shift component β r , and the (spatial) metric functions, χ, a, b depend on t, r. The electric field E µ = F µν n ν has only a radial component and the magnetic field B µ = F µν n ν vanishes, where n µ is the 4velocity of the Eulerian observer [39].
As in [24], we take as initial state a RN BH with ADM mass M , and charge Q, together with a vanishing scalar field. We set, as the scalar field initial data, a Gaussian distribution of the form with A 0 = 3 × 10 −4 , r 0 = 10M and λ = √ 8. As the geometrical initial data, we choose a conformally flat metric with a = b = 1 together with a time symmetry condition, i.e. vanishing extrinsic curvature, K ij = 0. The conformal factor is given by At t = 0, we choose a "pre-collapsed" lapse α 0 = ψ −2 and a vanishing shift β r = 0. Initially, the electric field is given by E r = Q r 2 ψ 6 . The logarithmic numerical grid extends from the origin to r = 10 3 M and uses a maximum resolution of ∆r = 0.0125M . A broader survey of the parameter space will be presented elsewhere.
In spherical symmetry, the evolution equations for the electric field and an extra variable, Ψ, to damp dynamically the constrains, take the form where K is the trace of K ij , we have taken κ 1 = 1 and Π ≡ −n a ∇ a φ.
To solve the Klein-Gordon equation we evolve the following system of first-order equations: 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.
νρ is the Levi-Civita connection associated to the spacetime metric g that one wants to determine, and a reference metricḡ is introduced, (Γ µ νρ being the corresponding Levi-Civita connection). Solutions to (C1) solve the Einstein equations iff ξ µ ≡ 0 everywhere. To achieve this, we impose boundary conditions which are compatible with ξ µ = 0 on the boundary of the domain of integration. Then, this should imply ξ µ ≡ 0 everywhere, a condition which is verified from the numerical output.
We use a metric ansatz with seven functions, F 1 , F 2 , F 3 , F 0 , S 1 , S 2 , S 3 which depend on (r, θ, ϕ): The reference metric is that of the RN BH, with D 1 = F 2 = F 3 = F 0 = 1, S 1 = S 2 = S 3 = 0 and N (r) where r H 0 is the event horizon radius andq < r H another input constant.
The ansatz for the Maxwell field is A = V (r, θ, ϕ)dt; for the scalar field it has the most general expression compatible with a static spacetime, φ = φ(r, θ, ϕ). The solutions have a nodeless scalar field n = 0, being indexed by two integers ( , m), which classify the corresponding scalar zero mode. The particular case of the axially symmetric (or even spherically symmetric) BHs can also be studied within this framework for m = 0. For example, the axially symmetric metric has S 2 = S 3 = 0, and all remaining functions depend on (r, θ) only.
In this approach, the problem reduces to solving a set of nine PDEs with suitable boundary conditions (BCs). The BCs are found by constructing an approximate form of the solutions on the boundary of the domain of integration compatible with the requirement ξ µ = 0, regularity of the solutions and asymptotic flatness.
In numerics, we have imposes the following BCs. At infinity: where Q is an input parameter fixing the electric charge; at θ = 0: together with ∂ θ φ = 0 for axially symmmetric BHs (m = 0) and φ = 0 in rest; at θ = π/2: together with φ = 0, except if + m is an even number, in which case we impose ∂ ϕ φ = 0; at ϕ = 0: at ϕ = π/2: together with φ = 0 for odd m and ∂ ϕ φ = 0 for even m. Finally, the numerical treatment of the problem is simplified by introducing a new (compact) radial coordinate x, such that r = r H 1−x 2 and 0 x 1, such that the horizon is located at x = 0. This results in the following boundary conditions at the horizon: A general large-r expression of the solutions can be constructed as a series in 1/r (e.g. F 0 = 1 + c t /r + . . . , where c t is a constant). The mass and the electric charge of the solutions are read off from the far field asymptotics where Φ is the electrostatic potential. Also, the solutions have an horizon area We have studied in a systematic way the solutions with boundary data corresponding to a = 1, m = 0, 1 and = 2, m = 0, 1, 2, for several values of the coupling constant α. The field equations are first discretized on a (r, θ, ϕ) grid with N r × N θ × N ϕ points. The grid spacing in the r direction is non-uniform, while the values of the grid points in the angular directions are uniform. Typical grids have sizes around 150 × 24 × 24. The resulting system is solved iteratively until convergence is achieved. All numerical calculations for non-spherical configurations are performed with a professional software [40], based on the iterative Newton-Raphson method. The code automatically provides an error estimate for each unknown function, which is the maximum of the discretization error divided by the maximum of the function. For most of the non-spherical solutions reported in this work, the typical numerical error is estimated to be around 10 −3 .