Charged black holes with axionic-type couplings: Classes of solutions and dynamical scalarization

Pedro G. S. Fernandes , Carlos A. R. Herdeiro, Alexandre M. Pombo , Eugen Radu, and Nicolas Sanchis-Gual 1,∥ Centro de Astrofísica e Gravitação CENTRA, Departamento de Física, Instituto Superior Técnico IST, Universidade de Lisboa UL, Avenida Rovisco Pais 1, 1049-001 Lisbon, Portugal Center for Research and Development in Mathematics and Applications (CIDMA), Campus de Santiago, 3810-183 Aveiro, Portugal Departamento de Matemática da Universidade de Aveiro, Campus de Santiago, 3810-183 Aveiro, Portugal


I. INTRODUCTION
Three types of bosonic fields, each consistent on its own as a classical relativistic field theory, are used in the physical description of nature: spin-0, -1, and -2 fields. According to the accepted models, the last two fields are realized, in their massless version, as the electromagnetic and gravitational field, whereas the scalar field is realized, in a massive, self-interacting version, as the Higgs boson. Beyond these concrete realizations, it has been a recurrent speculation in theoretical physics, ranging from quantum gravity to cosmology, that scalar fields play other roles. One particular example is a pseudoscalar field known as the axion. This (yet unobserved) particle was suggested in order to solve the strong CP problem, by Peccei and Quinn [1] (see also Refs. [2][3][4]). But it was later understood that, besides fulfilling its original purpose, the axion could have far-reaching implications in cosmology, namely, as a candidate for dark matter [5][6][7][8]. More recently, ultralight axionlike fields have been suggested to arise naturally from string theory compactifications [9,10], providing yet another possible origin in the context of a fundamental theory. The theoretical soundness of axions motivated experiments, both proposed and conducted, to detect axionic imprints; see, e.g., Refs. [11][12][13].
In the context of strong gravity, namely, black hole (BH) spacetimes, the interplay between the gravitational, electromagnetic, and nonstandard model scalar fields has a long history. If the scalar field is minimally coupled to the electromagnetic and gravitational field, remarkably, it cannot exist as an equilibrium configuration around spherical BHs [14]. This is a realization of the celebrated "no-hair" property of BHs; see, e.g., Refs. [15,16]. If nonminimal couplings are allowed, however, hairy BH solutions exist. In this paper, we shall be interested in such BHs when an axionic-type nonminimal coupling between the scalar and the electromagnetic field are considered, leading to charged BHs with axionic-type hair.
Nonminimal couplings between the electromagnetic and a scalar field, and the corresponding BH solutions, have long been considered in the context of, e.g., Kaluza-Klein theory and supergravity [17,18]. More recently, new families of such couplings have been under scrutiny in the context of BH spontaneous scalarization. Spontaneous scalarization is a strong gravity phase transition. It occurs when two phases (classes of solutions) coexist, one of them becoming dynamically preferred. First proposed by Damour and Esposito-Farèse in the context of neutron stars [19], it was later shown that BHs can undergo spontaneous scalarization in some models. The most recent interest, for asymptotically flat spacetimes, has been triggered by Refs. [20][21][22], centered in a class of models dubbed extended scalar-tensor Gauss-Bonnet. As pointed out in Ref. [23], however, in what concerns BH spontaneous scalarization, the latter class of models belongs to a wider universality class that also contains Einstein-Maxwellscalar (EMS) models. EMS models use a nonminimal coupling between the scalar and Maxwell invariant to induce scalarization and so require the presence of electric or magnetic charge. Using EMS models, it has been possible to establish that spontaneous scalarization of charged BHs occurs dynamically, leading to scalarized, perturbatively stable BHs [23][24][25][26][27].
In Ref. [28], a classification of EMS models was suggested in two main classes, according to the criterion that they admit or do not admit the electrovacuum Reissner-Nordström (RN) BH as a solution. The same can be done for the augmented EMS models, containing also an axionic-type nonminimal coupling, which will be the focus here. In the first class, the RN BH is not a solution, and the charged BHs all have nontrivial scalar profiles. For the augmented EMS models, this is the case of the canonical axion coupling. BHs in this model have been studied before [29][30][31][32][33][34], but we shall provide a fresh look at these solutions, exhibiting their full domain of existence. In the second class of models, the RN BH is a solution. We shall illustrate them in the augmented EMS models with a particular choice of nonminimal coupling that allows for the spontaneous scalarization of the RN BH. We shall construct the scalarized solutions and show they indeed arise dynamically from the evolution of the scalar-free RN BH, when a scalar perturbation is seeded. This paper is organized as follows. In Sec. II, we present the model, obtain the field equations for the ansatz that describes the class of solutions of interest, and provide some details on the construction of the solutions. We shall also review the criteria for the class of models wherein spontaneous scalarization can occur and describe how the emergence of the scalarized solutions is computed, in linear theory. Section III describes two physical relations used in assessing the accuracy of the obtained solutions. Section IV describes the effective potential for spherical perturbations, a simple tool that can establish perturbative stability and diagnose possible instabilities. Section V includes the bulk of our results, including a toy model, and the analysis of the two main examples considered. We close with a discussion of the obtained results.

II. MODEL
The model herein is a generalization of the EMS model discussed in Refs. [23,24]. It describes a real scalar field ϕ minimally coupled to Einstein's gravity and nonminimally coupled to the two relativistic and gauge invariants that can be built from the Maxwell tensor, F μν , in four spacetime dimensions: the Maxwell invariant, F μν F μν , and the parity violating invariant, F μνF μν . Here, where ϵ μνρσ denotes the Levi-Civita tensor density and g is the determinant of the spacetime metric. The generic model is described via the action (hereafter, we shall take units with 4πG The functions fðϕÞ and hðϕÞ determine the nonminimal couplings between the scalar field and the electromagnetic field. For hðϕÞ ¼ 0, the BH solutions have been constructed and discussed extensively in the literature, e.g., Refs. [17,18,23,24,28]. Here, we shall focus on constructing BH solutions and discussing some of their physical properties, in the case with hðϕÞ ≠ 0. Some previous work in this direction was considered, e.g., in Refs. [29,30].
Depending on the choices of the nonminimal coupling functions fðϕÞ, hðϕÞ, two types of models will be considered: models wherein all charged BHs have scalar hair and models wherein charged BHs can be either scalarized or scalar free. A generic, spherically symmetric line element to describe both the scalar-free and scalarized solutions is where NðrÞ ≡ 1-2mðrÞ=r and mðrÞ is the Misner-Sharp mass function [35]. Spherical symmetry requires the scalar field ϕðrÞ to have a radial dependence only, and an electromagnetic 4-potential ansatz of the following type, which allows a possible magnetic charge P, Integrating the trivial angular dependence, one obtains the effective Lagrangian, from which the equations of motion may be derived, with the prime denoting a radial derivative. Functions N, δ, V, ϕ have radial dependence only; for ease of notation, this dependence will be omitted henceforth.
To write the equations of motion, it is convenient to observe the existence of a first integral, where Q is an integration constant interpreted as the electric charge measured at infinity. Then, the equations of motion read ð2:7Þ where the dot denotes differentiation with respect to the scalar field, e.g., _ fðϕÞ ≡ df=dϕ. To solve this set of coupled, nonlinear ordinary differential equations, we have to implement suitable boundary conditions for the desired functions and corresponding derivatives. We assume the existence of an event horizon at r ¼ r H > 0 and that the solution possesses a power series expansion in ðr − r H Þ, ð2:9Þ Plugging these expansions in the field equations, the lowerorder coefficients are determined to be ð2:10Þ One observes that only two of the six parameters introduced in the expansions (2.9) are independent, which we choose to be ϕ 0 and δ 0 , the remaining being derived from these ones. The solutions in the vicinity of the horizon are determined by these two parameters, together with ðr H ; Q; PÞ. Some physical horizon quantities, such as the Hawking temperature T H , the horizon area A H , the energy density ρðr H Þ, and the Kretschmann scalar Kðr H Þ, are then determined by these five parameters as follows: To obtain the boundary conditions at spatial infinity, one performs an asymptotic approximation of the solution in the far field. Then, the equations of motion yield which introduce three new parameters: the scalar charge Q s , the electrostatic potential difference between the horizon and infinity Φ e , and the Arnowitt-Deser-Misner (ADM) mass M. From these asymptotic expansions, one collects a set of eight independent parameters, therefore, ðr H ; Q; P; ϕ 0 ; δ 0 ; Q s ; Φ e ; MÞ. As we shall see below, the full integration of the field equations relates these parameters, and for each choice of the coupling functions, the solutions of interest actually form a family of solutions with only three (continuous) parameters, typically taken to be the global charges ðM; P; QÞ, but possible with nonuniqueness. For later use, we collect the results and definitions where ξ represents any of the global charges of the model, while fðϕÞ and hðϕÞ remain unchanged.

A. Models with spontaneous scalarization
Depending on the choice of the functions fðϕÞ and hðϕÞ, the model may accommodate either only BH solutions with a nontrivial scalar field profile or both scalar-free and scalarized charged BHs. In the latter case, the scalar-free BH solutions may become unstable, within some region of the parameter space, dynamically evolving to the scalarized solutions. This is the phenomenon of spontaneous scalarization of charged BHs [23,24]. In this subsection, we identify the conditions for scalarization to occur and study its onset, at linear level.

Spontaneous scalarization conditions
For spontaneous scalarization of charged BHs to occur, the coupling functions (and their derivatives) must satisfy two key conditions [20][21][22][23][24]. These read as follows: (1) The system must accommodate a scalar-free solution. Since the Klein-Gordon equation of motion is the existence of a scalar-free solution requires (for both Q, P nonzero) (2) Spontaneous scalarization occurs if the scalar-free solution is unstable against scalar perturbations δϕ. These obey (neglecting second-order terms) The instability occurs if the mass is tachyonic, i.e., μ 2 eff < 0. This constrains the second derivatives of the coupling functions. Besides these conditions, one may require, for convenience, that Maxwell's theory should be recovered near spatial infinity. Thus, No condition is imposed on the coupling function hðϕÞ, since the term F μνF μν is topological and hence nondynamical when the corresponding scalar coupling becomes constant.

Bifurcation of solutions: The existence line
Let us now consider the onset of spontaneous scalarization. We assume that the model under consideration admits the dyonic RN BH of Einstein-Maxwell theory as the scalar-free solution, that is, The scalarization phenomenon is assessed by considering scalar perturbations of the RN solution within the considered model. Following Refs. [23,24], we take a spherical harmonics decomposition of the scalar field perturbation: With this ansatz, the scalar field equation of motion (2.16) simplifies to Once the coupling functions are fully fixed, solving (2.23) is an eigenvalue problem: for a given l, requiring an asymptotically vanishing, regular at the horizon, smooth scalar field, a discrete set of BHs solutions is selected, i.e., a discrete set of RN solutions, each with a certain reduced charge q. These are the bifurcation points from the scalar-free solution. They are labeled by an integer n ∈ N 0 ; n ¼ 0 is the fundamental mode, whereas n ≥ 1 are excited states (overtones). The RN solutions with a smaller (larger) q than that of the bifurcation point are stable (unstable) against the corresponding scalar perturbation. In particular, the first bifurcation point, i.e, the one with the smallest q, which corresponds to the mode l ¼ 0 and n ¼ 0, marks the onset of the scalarization instability. Only RN BHs with q smaller than the first bifurcation point are stable against any sort of scalar perturbation. At each bifurcation point, a new family of (fully nonlinear) scalarized BH solutions emerges from the RN family, as static solutions of the equations of motion of the full model. In this paper, we shall consider only the first bifurcation point and the corresponding new family of spherically symmetric scalarized BHs that bifurcate from the RN family. To be concrete, let us take fðϕÞ ¼ 1 and hðϕÞ ¼ −αϕ 2 (with α > 0), a case we shall consider in detail later. Then, μ 2 eff ¼ −2αQP=r 4 , which confirms the tachyonic instability. 2 For l ¼ 0, one finds the exact solution where LP u is a Legendre function. The function UðrÞ approaches a constant nonzero value as r → ∞, ð2:25Þ thus, finding the l ¼ 0 unstable mode of the RN BH amounts to a study of the zeros of the Legendre functioncf. Fig. 1 (left panel). Bifurcation requires α above a minimum value given by α min ¼ ð1 þ β 2 Þ=ð8βÞ. Some illustrative values can be found in Table I. For β ≤ 1, a smaller β requires a larger α for bifurcation, and bifurcation occurs for a larger q, in contrast to the dyonic case [28]. Observe that for P ¼ 0 the Legendre function trivializes, since FF vanishes. Varying α, the set of bifurcation points constitutes the existence line, i.e., the set of RN configurations for which the zero mode exists. This line is shown in Fig. 1 (right panel) for some illustrative values of β in a q vs α plot.

III. PHYSICAL RELATIONS
Let us now briefly consider two physical relations that, besides their physical content, are used to test the numerical accuracy of the solutions found numerically. These are a Smarr-type law and a virial-type relation.

A. Smarr-type law
The Smarr law [36,37] provides a relation between the total mass of the spacetime and other measurable quantities, like the horizon temperature and area. Its information complements that of the equations of motion, making it an interesting test to assess the accuracy of BH solutions obtained numerically.
The Smarr law can be obtained via the integral mass formula, which for our model reads where k a is the Killing vector field associated to staticity and T is the trace of the energy-momentum tensor. The full energy-momentum tensor is which is insensitive to the F μνF μν term in the action, as the latter is a topological invariant. One can thus arrive at the Smarr-type law where we refer to as the electrostatic and magnetostatic potential differences, respectively, and is interpreted as the axion-related electromagnetic energy. A nonlinear Smarr-type relation can also be established for this family of models, following Ref. [23], which reads Scaling arguments, initiated by the work of Derrick [38], are a powerful tool to establish no-go theorems for solitonic solutions (see, e.g., Ref. [39]) and no-hair theorems for BH solutions [15] as well as to provide a physical relation that must be obeyed by solutions of a given model. These relations are generalizations of the canonical virial theorem, which states an energy balance, and are often described as virial relations. They are typically independent from the equations of motion; thus, again, they are useful in assessing the accuracy of numerically generated solutions.
Consider the effective action and assume that a charged BH solution with scalar hair exists, described by the functions ϕðrÞ, δðrÞ, VðrÞ, NðrÞ, with suitable boundary conditions at the event horizon and at infinity. Next, consider the one-parameter family of configurations described by the scaled functions with F ∈ fϕ; δ; V; Ng. If the initial configuration had indeed been a solution, then the effective action for the scaled configurations must possess a critical point at λ ¼ 1: One can show that the left-hand side integrand is strictly positive. Thus, the virial identity shows that a nontrivial scalar field requires a nonzero electric/magnetic charge so that the right-hand side is nonzero. As an immediate corollary, neutral BHs cannot be hairy in this model.

IV. EFFECTIVE POTENTIAL FOR SPHERICAL PERTURBATIONS
Let us also introduce a diagnosis analysis of perturbative stability, against spherical perturbations, that shall be applied to the solutions derived and discussed in the next section.
Following a standard technique, see e.g., Ref. [24], we consider spherically symmetric, linear perturbations of an equilibrium solution, keeping the metric ansatz but allowing the functions N, δ, ϕ, V to depend on t as well as on r: The time dependence enters as a Fourier mode with frequency Ω, for each of these functions: Nðr; tÞ ¼ NðrÞ þ ϵN 1 ðrÞe −iΩt ; δðr; tÞ ¼ δðrÞ þ ϵδ 1 ðrÞe −iΩt ; ϕðr; tÞ ¼ ϕðrÞ þ ϵϕ 1 ðrÞe −iΩt ; Vðr; tÞ ¼ VðrÞ þ ϵV 1 ðrÞe −iΩt : ð4:2Þ From the linearized field equations around the background solution, the metric perturbations and V 1 ðrÞ can be expressed in terms of the scalar field perturbation, thus yielding a single perturbation equation for ϕ 1 . Introducing a new variable ΨðrÞ ¼ rϕ 1 , the scalar-field equation of motion may be written as which, by introducing the "tortoise" coordinate x as dx=dr ¼ e δ =N [40], can be written in the standard onedimensional Schrödinger-like form: The effective potential that describes spherical perturbations U Ω is defined as ð4:7Þ An unstable mode would have Ω 2 < 0, which for the asympotic boundary conditions of our model is a bound state. It follows from a standard result in quantum mechanics (see, e.g., Ref. [41]), however, that Eq. (4.5) has no bound states if U Ω is everywhere larger than the lowest of its two asymptotic values, i.e., if it is positive in our case. 3 Thus, an everywhere positive effective potential proofs mode stability against spherical perturbations. We remark that the existence of a region of negative potential is a necessary but not sufficient condition for instabilities to be present. In fact, for the fundamental, spherically symmetric scalarized solutions in Ref. [23], this region occurs for some solutions, which are, nonetheless, stable [25].

V. RESULTS FOR SPECIFIC EXAMPLES
WITH f ðϕÞ = 1 Let us now apply the foregoing formalism to particular cases. We shall focus on cases for which scalarization is sourced by the nonminimal coupling to the F μνF μν term, keeping fðϕÞ ¼ 1. Cases for which fðϕÞ ≠ 1 and hðϕÞ ¼ 0 have been considered more extensively in the literature; see, e.g., Refs. [17,18,23,24,28].

A. Flat spacetime toy model
Following Ref. [23], we start by presenting a flat spacetime toy model that illustrates the spontaneous scalarization phenomenon with an axionic-type coupling and, moreover, shows that gravity is not fundamental for the phenomenon to occur. 3 A simple proof is as follows. Write Eq. (4.5) in the equivalent After integrating from the horizon to infinity, it follows that Consider the particular case in which the background is the flat Minkowski spacetime, i.e., Eq. (2.2) with NðrÞ ¼ 1 and δ ¼ 0. The corresponding Maxwell-scalar model allows a scalar-free configuration, which is simply the dyonic Coulomb configuration, given by The scalar field equation of motion (2.7) reduces to Introducing γðϕÞ ≡ ðQ þ PhðϕÞÞ 2 and a new radial coordinate x ¼ 1=r, one arrives at which takes the form of a one-dimensional mechanical problem of a particle in a potential. The corresponding motion is described by the Lagrangian L ¼ T − U, with T ¼ ðdϕ=dxÞ 2 and U ¼ −γðϕÞ. We notice the existence of a first integral E 0 ¼ T þ U ¼ ðdϕ=dxÞ 2 − γðϕÞ, from which analytical solutions for the radial profile of the scalar field may be obtained by solving the integral and then inverting the solution to obtain ϕðxÞ and ϕðrÞ.
For an explicit example, consider the coupling function which respects all the previously discussed conditions for the occurrence of spontaneous scalarization. For this particular choice of coupling function, which does not have any particular motivation rather than illustrating the phenomenon with a simple solution, it is possible to obtain a closed-form spherically symmetric scalarized solution of the Maxwell-scalar system, where Ell denotes the elliptic integral of the second kind. This solution coexists with the scalar-free dyonic Coulomb solution, and one may inquire which one is preferred.
To address this point, and since the Coulomb solution has infinite energy due to the singularity at the location of a point charge, we consider a cutoff in the form of a conducting sphere located at r ¼ r 0 . The energy of any of the configurations can be computed as Denoting by E ϕ≠0 and by E ϕ¼0 the energy of the scalarized and scalar-free configurations, respectively, one obtains hence, the scalarized solution is energetically favored in a set of bands given by where n ∈ N 0 , besides labeling the bands, counts the number of nodes exterior to r 0 of the scalar field profile. As advertised, this toy model shows spontaneous scalarization through an axionic-type coupling; it is not exclusive of gravitational theories, and it may lead to favored configurations.

B. BHs with axionic hair (hðϕÞ = − αϕ)
Depending on the choice of the coupling function hðϕÞ, two sorts of models are possible: models in which the dyonic RN BH of electrovacuum is a solution and models in which it is not. In Ref. [28], a classification of EMS models-without an axionic coupling-was performed based on this criterion. Here, we shall illustrate these two sorts of models with two examples of hðϕÞ. In this subsection, we choose hðϕÞ ¼ −αϕ, which is the standard axionic coupling, for which the dyonic RN BH is not a solution and there are new BHs, all of which have axionic hair. In this case, the scalar field ϕ is dubbed "axion." The BHs of the model with hðϕÞ ¼ −αϕ were first constructed in Ref. [29]; our analysis will complement the study therein. Since all BH solutions are scalarized, no bifurcation points from the scalar-free RN BH exist. We will study the domain of existence of such BHs with axionic hair, which has not been presented before. For particular solutions, the radial profiles will be compared with the results in Ref. [29], obtaining agreement. We shall also obtain, for these BHs with axionic hair, the effective potential for spherical perturbations.

Radial profiles
Some typical solutions of the various functions that define scalarized BHs obtained from numerical integration can be found in Fig. 2 for two illustrative values of the coupling constant α ¼ 15, 35, while keeping Q, P, and r H constant. Some characteristic quantities can be found in Table II.
All BHs with axionic hair we have constructed have an axion profile that is nodeless. The axion is a monotonically decreasing function of the radius; hence, ϕ 0 is always the maximum value of the axion field. Data reveal that (as expected) for larger couplings α there is a higher degree of scalarization, as seen from the value of ϕ 0 in Fig. 2 and in the values of the scalar charge (Table II).
In Ref. [29], two approximations for the axion radial profile were derived from Eq. (2.7) as ϕðrÞ ≈ α QP r þ r − ln r r − r − þ Oðα 2 Þ; for small α; ð5:10Þ p denote the radial coordinate of the horizons of the scalar-free RN BH. A comparison between these analytical approximations and the numerical results can be seen in Fig. 3.
We have introduced two improved analytical approximations for large and small couplings (in the same spirit as in Ref. [29]) that besides taking into account the value ϕ 0 for the axion at the event horizon in the small coupling limit does not neglect α 2 terms. Both approximations can be obtained starting from (2.7), which for this particular choice of coupling functions takes the form In the small coupling limit, the axion field is considered as perturbations around the dyonic Reissner-Nordström metric (2.21) (N ¼ 1-2Mx þ ðQ 2 þ P 2 Þx 2 and δ ¼ 0). By solving the differential equation (5.13) in this limit, one obtains for small α; ð5:14Þ where LP and LQ are Legendre functions of first and second kinds, respectively, and v ¼ 1 2 ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi and B are integration constants chosen to guarantee an asymptotically vanishing solution; the horizon value of the axion is ϕ 0 . The opposite approximation is also obtained from (5.13), by considering the analogous flat space limit (N ¼ 1 and δ ¼ 0) [29]. The solution of the differential equation (5.13) is then for large α: ð5:15Þ These new approximations are compared with the previous ones and the numerics in Fig. 4. Data reveal that, indeed, the new approximations hold better than the previous ones, with the differences being more relevant in the small coupling case. Although less obvious from Fig. 4, numerics indicate that the new large coupling approximation also gives an overall better approximation. For code tests on the numerical solutions, we report differences of 10 −9 for the Virial relation and 10 −8 for the Smarr law.

Domain of existence
The domain of existence, in the ðα; qÞ plane, for BHs with axionic hair is presented in Fig. 5 (left panel). It is delimited by a set of critical solutions that we call the critical line. At the critical line, numerics suggest a divergence of the Kretschmann scalar and of the horizon temperature, together with a vanishing of the horizon area. Simultaneously, M and Q s remain finite and nonzero. However, for sufficiently small values of the coupling, for which the axionic hair is almost negligible, some physical quantities behave similarly to those for an extremal RN BH, e.g., the event horizon temperature (cf. Fig. 5, right panel). It is interesting to notice that, along α ¼ constant branches, q increases beyond unity; therefore, BHs with axionic hair can be overcharged.
For a given α, as q increases, so does the axion's initial amplitude ϕ 0 . Thus, the global maximum of the axion, which always occurs at the horizon, increases, for fixed α, with q. This means one may take ϕ 0 as a measure of q and vice versa. As another feature, fixing Q, there is a wider domain of existence for larger P. This behavior contrasts with the one observed for the dyonic scalarized BHs in Ref. [28] [which has hðϕÞ ¼ 0 and P ≠ 0]. This wider domain of existence, for larger β, may be associated to the fact that the coupling hðϕÞF μνF μν is directly proportional to P; thus, for larger values of β, the axion couples more strongly to the source term.

Effective potential for spherical perturbations
The effective potential for spherical perturbations U Ω for a sample of BHs with axionic hair is plotted in Fig. 6. It is not positive definite in all cases, but it is regular in the entire range −∞ < x < ∞. For other values of the coupling α, the potential always behaves in a similar way: for sufficiently small values of β, the potential is always positive (and hence, free of instabilities), until β reaches a value at which the potential starts to have a negative region. The potential vanishes at the horizon and at infinity.

C. Spontaneous scalarization with an axionic-type coupling [hðϕÞ = − αϕ 2 ]
We now turn to the second sort of models, for which the dyonic RN BH of electrovacuum is a solution but new, dubbed scalarized BHs, are also possible. In these models, the phenomenon of spontaneous scalarization may occur. To illustrate it, we consider the coupling function hðϕÞ ¼ −αϕ 2 [and fðϕÞ ¼ 1 as before], hereby dubbed "power-law coupling." The coupling constant α is a dimensionless positive constant. This coupling function obeys all the necessary conditions for the occurrence of spontaneous scalarization, and it is the natural generalization of the axionic case. Again, we shall present radial profiles and the domain of existence of the scalarized BHs. Moreover, since the RN BH is a solution of the model (the scalar-free solution), an entropic analysis that will show the scalarized solutions are entropically favored will be performed. Then, a fully nonlinear dynamical analysis will show that the scalar-free solutions indeed dynamically evolve to form scalarized BHs.

Radial profiles
Some typical profiles for the various functions that define scalarized BHs, in the power-law coupling model, obtained from numerical integration, can be found in Fig. 7, for two illustrative values of the coupling constant, α ¼ 20, 35, while keeping Q, P, and r H constant (some characteristic quantities can be found in Table III). Two excited solutions (with n ¼ 1 and n ¼ 2) are also presented.
Some qualitative features observed for the BHs with axionic hair remain here. For the nodeless solutions, ϕ 0 is the maximum of the scalar field radial profile. Larger couplings imply stronger scalarization and smaller electrostatic potentials. The excited states of Fig. 7 (bottom panel) have larger δ 0 and electrostatic potentials when compared to the fundamental solutions.
From Eq. (2.7), it is also possible to derive two approximations, for small and large couplings, in a similar fashion as it was done for the case of the axionic coupling, where the integration constants were chosen such that the solution is regular and asymptotically vanishing. A comparison between these analytical approximations and the numerical results can be found in Fig. 8. The approximations are qualitatively worse than in the axionic case. Likely, this is due to the nonlinearities induced by the quadratic scalar field terms. The small coupling approximation neglects α 2 terms, and the large coupling one does not possess the benefit of imposing the scalar field value at the event horizon, since one of the integration constants was used to impose regularity of the approximation. We remark that, once again, tests to the code reveal relative differences 10 −9 for the Virial relation and 10 −7 for the Smarr law.

Domain of existence
In the case studied in the previous section, dyonic BHs with axionic hair, the domain of existence was, generically, bounded by a critical line: a set of singular solutions in which the horizon shrinks to a point. This is also the case for the purely electrically charged scalarized BHs studied in Refs. [23,24]. A different behavior was observed for the dyonic BHs in Ref. [28]. The domain of existence is then bounded by an extremal line: a set of solutions for which the horizon temperature tends to zero, while the Kretschmann scalar and the horizon area remain finite and nonzero. We will now see that for the dyonic scalarized BHs with the power-law coupling herein the latter behavior also occurs.
To obtain the extremal BH solutions, a different nearhorizon expansion that accounts for a degenerate horizon must be used (see, e.g., Refs. [28,42,43]). A double zero for NðrÞ at the horizon is accommodated with the following expansion: Again, the field equations related these parameters. We take ϕ 0 and r H as the independent parameters. Then, In this way, we obtained directly extremal solutions, which make up a boundary of the domain of existence. This domain is presented in Fig. 9, in the ðα; qÞ plane. The boundary of the domain of existence is composed of several extremal BH solutions, each characterized by a single value of p that is unique and exclusive of that solution. The domain of existence reveals a region of nonuniqueness where, for the same charge-to-mass ratio q < 1, RN BHs and scalarized BHs coexist. Unlike in the domain of existence of the axionic case, for larger β, there is a narrower domain of existence. This trend, however, is not universal. For small enough α, the domain of existence can be wider for larger β- Fig. 10 (left panel). Figure 10 (right panel) shows, in fact, a more subtle behavior. The solid blue lines are extremal solution for α ¼ constant. One can observe that there is an optimal value of β for overcharging. As α increases, however, this optimal value tends to the minimum value of β.

Effective potential for spherical perturbations
The effective potential U Ω for the scalarized BHs in the power-law coupling model is plotted in Fig. 11 for some illustrative examples. Again, it vanishes at the BH event horizon and at infinity, but it is not positive definite. For other values of the coupling and β ratios, the potential always has a similar form. The existence of a region of negative potential does not imply instability. Thus, conclusions about linear stability require a study of  quasinormal modes, similarly to what was done in Ref. [25] for the purely electrical case. Below, however, we shall present evidence, using fully nonlinear numerical simulations that the scalarized solutions are not only stable but, in fact, form dynamically.

Entropic and dynamical preference
Since the model under consideration is general relativity coupled to some matter, the Bekenstein-Hawking BH entropy formula holds. Thus, the entropy analysis reduces to the analysis of the horizon area. It is convenient to use the reduced event horizon area a H . Then, in the region where the RN BHs and scalarized BHs coexist-the nonuniqueness region-for the same q, the scalarized solutions are always entropically preferred, as seen in Fig. 12.
The entropic preference, together with the instability of the scalar-free RN BHs against scalarization, suggest that the latter evolve toward the former, when perturbed, at least if the evolutions are approximately conservative. We have performed fully nonlinear numerical evolutions to test this scenario, following our previous work [23,24]. The initial data are a dyonic RN BH, with ADM mass M, electric charge Q, and magnetic charge P, which was evolved using the numerical code described in Refs. [23,44,45], adapted to the power-law coupling hðϕÞ ¼ −αϕ 2 . The code uses spherical coordinates under the assumption of spherical symmetry employing the second-order partially implicit Runge-Kutta method developed by Refs. [46,47].
The 3 þ 1 metric is given by ds 2 ¼ −ðα 2 0 þ β r β r Þdt 2 þ 2β r dtdr þ e 4χ ½adr 2 þ br 2 dΩ 2 , where the lapse α 0 ; shift component β r ; and the (spatial) metric functions χ, a, and b depend on t, r. As in Ref. [23], the scalar (initial data) perturbation to trigger the instability is a Gaussian distribution of the form For the dyonic RN initial data, the conformal factor is given by : ð5:22Þ At t ¼ 0, we choose a "precollapsed" lapse α 0 ¼ ψ −2 and a vanishing shift β r ¼ 0. Initially, the electric field is given by E r ¼ Q r 2 ψ 6 and B r ¼ P r 2 ψ 6 . Since we are considering a RN BH with magnetic charge, we also have to take into account the evolution equation of the magnetic field, even in spherical symmetry. The evolution equations for the electric and magnetic fields and two extra variables, Ψ E and Φ B , to damp dynamically the constrains, take the form in spherical symmetry ∂ t E r ¼ β r ∂ r E r − E r ∂ r β r þ ðα 0 KE r − D r Ψ E Þ þ 2αα 0 ϕΠB r ; ð5:23Þ ∂ t B r ¼ β r ∂ r B r − B r ∂ r β r þ ðα 0 KB r þ D r Φ B Þ; ð5:24Þ where K is the trace of the extrinsic curvature K ij , Π ≡ −n a ∇ a ϕ, and we have chosen the damping terms κ 1 ¼ κ 2 ¼ 1.
The Klein-Gordon equation is solved by evolving the following system of first-order equations: The matter source terms in the Einstein equations are free from the axionic coupling.
The numerical simulations show that the scalar perturbation triggers the spontaneous scalarization of the RN BH. The horizon electric charge decreases as the energy of the field increases, while the horizon magnetic charge remains unchanged, until we reach equilibrium and a scalarized solution forms at the end point of the dynamical scalarization. The scalar cloud grows near the horizon and expands radially. The radial profile of the cloud decreases monotonically with increasing radii. In Fig. 13, we plot the scalar field value at the horizon ϕðr H Þ for both the dynamical evolutions and the static solutions with the same Q and P. We obtain quite good agreement with the static solutions described above. As in our previous works [23,24], the dynamical solutions better match the static ones for lower values of Q, showing that the evolutions are more approximately conservative in that case.

VI. CONCLUSIONS
This paper considers an augmented EMS sclar model and its BH solutions, in particular in the context of the spontaneous scalarization of charged BHs. The model's novelty consists of adding an axioniclike coupling hðϕÞF μνF μν to the EMS action, which has been previously studied in the context of BH spontaneous scalarization. Depending on the choice of hðϕÞ, the model can accommodate BHs with axioniclike and, possibly, also the standard RN BH of electrovacuum. In this case, the latter may become unstable against scalar perturbations and spontaneously scalarized, which we have shown to occur dynamically in one illustrative example.
Spontaneous scalarization provides a dynamical mechanism for new BHs to emerge under particular circumstances. Although it is commonly considered that electric charges are astrophysically irrelevant (but see e.g., Ref. [48]), the model herein can be considered as a toy model for spontaneous scalarization induced by higher curvature corrections, which may have astrophysical FIG. 13. Scalar field value at the horizon for different values of Q, while β ¼ 0.5, in the power-law coupling model. The solid lines are obtained from data of the scalarized BHs obtained as static solutions of the field equations. The individual points are obtained from the numerical evolutions, starting from a scalarfree RN BH with the same global charges M, Q, and P. The agreement is better for the lower charges, showing that scalarization only redistributes the electric charge and energy between the horizon and the scalar hair, with minor leaking toward infinity. This leaking appears to become more significant for larger charges. relevance. In this context, for instance, it was pointed out recently that the observed shadow of the BH in the center of the M87 galaxy [49] is compatible with spontaneously scalarized Kerr BHs within some range of the parameters of the model [50]. However, the spontaneous scalarization with higher curvature parity-violating terms, such as the Chern-Simons term, has not yet been considered, except in the academic model in Ref. [51]. 4

ACKNOWLEDGMENTS
This work has been supported by Fundação para a Ciência e a Tecnologia (FCT), within Project No.