Axionic instabilities and new black hole solutions

The coupling between scalar and vector fields has a long and interesting history. Axions are one key possibility to solve the strong CP problem and axion-like particles could be one solution to the dark matter puzzle. Given the nature of the coupling, and the universality of free fall, nontrivial important effects are expected in regions where gravity is strong. Here, we show that i. A background EM field induces an axionic instability in flat space, for large enough electric fields. Conversely, a homogeneous harmonic axion field induces an instability in the Maxwell sector. When carried over to curved spacetime, this phenomena translates into generic instabilities of charged black holes (BHs). ii. In the presence of charge, BH uniqueness results are lost. We find solutions which are small deformations of the Kerr-Newman geometry and hairy stationary solutions without angular momentum, which are `dragged' by the axion. Axion fields must exist around spinning BHs if these are immersed in external magnetic fields. The axion profile can be obtained perturbatively from the electro-vacuum solution derived by Wald. iii. Ultralight axions trigger superradiant instabilities of spinning BHs and form an axionic cloud in the exterior geometry. The superradiant growth can be interrupted or suppressed through axionic or scalar couplings to EM. These couplings lead to periodic bursts of light, which occur throughout the history of energy extraction from the BH. We provide numerical and simple analytical estimates for the rates of these processes. iv. Finally, we discuss how plasma effects can affect the evolution of superradiant instabilities.

The coupling between scalar and vector fields has a long and interesting history. Axions are one key possibility to solve the strong CP problem and axion-like particles could be one solution to the dark matter puzzle. Extensive experimental and observational efforts are actively looking for "axionic" imprints. Given the nature of the coupling, and the universality of free fall, nontrivial important effects are expected in regions where gravity is strong. Rotating black holes (immersed, or not in magnetic fields) are a prime example of such regions. Here, we show that (i) A background electromagnetic field induces an axionic instability in flat space, for electric fields above a certain threshold value. Conversely, a homogeneous harmonic axion field induces an instability in the Maxwell sector. When carried over to curved spacetime, this phenomena translates into generic instabilities of charged black holes. We describe the instability and its likely final state, new black hole solutions.
(ii) In the presence of charge, black hole uniqueness results are lost. We find solutions which are small deformations of the Kerr-Newman geometry and hairy stationary solutions without angular momentum but which are "dragged" by the axion. Axion fields must exist around spinning black holes if these are immersed in external magnetic fields. The axion profile can be obtained perturbatively from the electro-vacuum solution derived by Wald.
(iii) Ultralight axions trigger superradiant instabilities of spinning black holes and form an axionic cloud in the exterior geometry. The superradiant growth can be interrupted or suppressed through couplings such as E·B (typical axionic coupling) but also more generic terms such as direct couplings to the invariant E 2 − B 2 . These couplings lead to periodic bursts of light, which occur throughout I. INTRODUCTION Peccei and Quinn first introduced the axion, a pseudo-Goldstone scalar field, in order to solve the strong CP problem [1]. More recently it was shown that a plenitude of ultralight axion-like bosons might arise from moduli compactification in string theory. In this "axiverse" scenario, a landscape of light axion-like fields can populate a mass range down to the Hubble scale m H ∼ 10 −33 eV, with possible implications at astrophysical and cosmological scales [2]. In particular, axions and axion-like particles are also strong candidates for both cold and non-cold dark matter [3][4][5][6].
The photon-axion mixing in the presence of an external magnetic or electric field, can be used to impose strong constraints on axion-like particles due to intergalactic magnetic fields (see e.g. Ref. [7] for a review) or can lead to a detectable signature in the spectra of high-energy gamma ray sources [8]. In addition, even in the absence of electromagnetic (EM) fields, superradiant instabilities can be triggered, through which the axion grows and "condensates" around massive, spinning black holes (BHs) [9][10][11][12][13][14]. The instability extracts rotational energy away from the spinning BH and deposits it into an axion cloud with high occupation number [14]. Eventually, gravitational-wave (GW) emission dominates over the superradiant growth, leading to a secular spin-down and decay of the cloud. Such systems are a promising source of GWs that can be detected with current and future detectors [15][16][17][18][19][20][21][22].
In an astrophysically relevant situation, BHs are often surrounded by a plasma in an accretion disk, which generates its own EM field. In addition, galactic magnetic fields and background EM radiation is present. The presence of magnetic fields in regions where gravity is strong may give rise to new phenomena, for example the triggering of instabilities or the induction of nontrivial BH hair. In the context of axions or axion-like particles, such scenarios have hardly been considered. One of the purposes of this work is to understand possible new BH configurations in the context of axionic physics.
It has also been argued recently that the coupling of superradiant axion clouds with photons can lead to bursts of radiation which in the quantum version resemble laserlike emission [23,24]. Thus, the evolution of superradiant instabilities would produce periodic emission of light. These arguments are order-of-magnitude, highly approximate and partially inconsistent, but have very recently been put on a firmer ground through the full numerical solution of the relevant equations [25]. More generally, the study of axion electrodynamics in curved spacetime has been the topic of some but few studies, with some results in the Schwarzschild background in the context of Pulsar magnetospheres [26] and polarization of EM waves passing through the scalar clouds around BHs [27]. The other purpose of this work is to explore thoroughly the possible instabilities associated to axionic or other scalar-type couplings to the Maxwell sector occuring in some dark matter models.
In summary, we find that the presence of electric charge or rotation leads to the appearance of new BH solutions with nontrivial axion hair. Axionic of scalartype couplings to the Maxwell sector are also found to trigger strong instabilities that affect specially superradiant clouds around BHs.
We adopt geometric units (G = c = 1) throughout, and a "mostly plus" signature.

II. SETUP
There are many possible and viable DM candidates [28]. We will focus our attention on what are perhaps the best motivated extensions of the Standard Model and of General Relativity, which include a massive (and real) scalar Ψ with possible axionic couplings to a vector (through the coupling constant k a ) and scalar couplings to the Maxwell invariant through a coupling constant k s , with p = 1, 2 being popular choices [29,30]. Thus, depending on the parity transformation of the (pseudo)scalar, coupling to the Maxwell sector is realized through E · B (pseudo-scalar) or E 2 − B 2 (scalar) invariant. We will study these couplings separately. Other couplings, such as ∇ µ Ψ∇ ν ΨF µν and ∇ µ Ψ∇ ν ΨA µ A ν are possible, but we will not consider these here. The mass of the scalar Ψ is given by m S = µ S , F µν ≡ ∇ µ A ν − ∇ ν A µ is the Maxwell tensor and * F µν ≡ 1 2 µνρσ F ρσ is its dual. We use the definition µνρσ ≡ 1 √ −g E µνρσ where E µνρσ is the totally anti-symmetric Levi-Civita symbol with E 0123 = 1. The quantities k a , k s are constants. We get the following equations of motion for the theory above: If Ψ is the QCD axion, √ k a ∼ 10 15 10 −5 eV µ S GeV .
The inverse energy parameter k s is tightly constrained for p = 1 but much less so for p ≥ 2 [29]. Astrophysical BHs can probe scalar fields with masses in the range ∼ [10 −20 , 10 −11 ] eV [21], but larger scalar field masses can in principle also be probed by hypothetical subsolar mass primordial BHs. Depending on the formation scenario, the QCD axion with masses in the range ∼ [10 −12 , 10 −2 ] eV is a strong dark matter candidate (see e.g. [31]). Current experiments are especially sensitive to this mass range but most of the relevant parameter space is not yet ruled out (see e.g. Fig. 1 in [32]). For masses below 10 −10 eV and arbitrary coupling constants the parameter space remains largely unconstrained. This range is especially interesting for stringy ultralight axions which could also explain part or all of the dark matter for axions with masses down to ∼ 10 −23 eV [5]. Here, we consider arbitrary coupling constants to keep the discussion as general as possible.
For all practical purposes, the right-hand-side of Eq. (2c) can be set to zero when the axionic coupling constant and the BH charge are is small: In natural units, the strength of a magnetic field around a source of mass M can be measured defining the characteristic magnetic field B M = 1/M associated to a spacetime curvature of the same order of the horizon curvature. In physical units this is given by B M ∼ 2.4 × 10 19 (M /M ) Gauss. We will use this approximation when looking for new BH solutions and also when performing dynamical simulations of superradiant clouds. Although our results are only valid when B/B M 1 (when backreaction on the metric is small), this limit includes the most interesting region of the parameter space. Indeed, for astrophysical BHs, a reference value for the largest magnetic field that can be supported in an accretion disk is given by B ∼ 4 × 10 8 (M/M ) −1/2 Gauss [33] so that the approximation B B M is well justified. We also neglect the backreaction of the scalar field onto the geometry, an approximation which is justified both perturbatively and numerically [17,21].

III. FLAT-SPACE INSTABILITIES
The theory above shows interesting aspects even in flat space. In fact, most of the strong-field effects that we will deal with could have been guessed from a flat-space analysis. We find that non-vanishing background EM or axion fields both may trigger instabilities, but the nature and details of such instability depends on which background field one discusses: a static background electric field gives rise to instabilities of flat-space. However, to have a similar effect for background axions, one needs a time-varying axion (or scalar).

A. Homogeneous background EM field
We start by exploring EM and scalar or axion fluctuations, determined by Eqs. (2a) and (2b), in the background of a homogeneous EM field.

Axionic couplings
Let us turn off the scalar coupling k s for now. A vanishing scalar field and constant background electric (in standard cartesian coordinates) E = (0, 0, E z ) and magnetic field B = (0, B y , 0) solves the equations of motion. Let's perturb this background solution by allowing a small but nonvanishing field Ψ ∼ η e −i(ωt−pix i ) and fluctuations in the vector field, where X µ , µ = 0, . . . , 3 are constants and is a small book-keeping parameter. The Klein-Gordon equation yields where p 2 ≡ p i p i . Maxwell's equations can be used to obtain X 0 , X 1 , X 2 , Finally, substituting back in (5) one finds the dispersion relation This equation can be solved for ω, with the result that an instability appears at a threshold value of the electric field. For example, for B y = 0, p z = 0 one gets the dispersion relation Thus, for E > E crit an instability (i.e., negative ω 2 , such that ω has a positive imaginary component and the fluctuations diverge exponentially in time) appears, with The instability requires a non-homogeneous fluctuation in the axion, something that will prove important in the discussion of BH instabilities. These results are in agreement with recent studies on axion-like phenomena in magnetic materials [34], and predict acritical electric field scaling like the inverse of the coupling constant k a .

Scalar couplings
Let us now focus on the scalar coupling only, i.e., we set k a = 0, k s = 0. We take the same ansatz for the scalar field and vector potential given by Eq. (4).
For the scalar coupling with p = 1, a background EM field also requires a background scalar. If B y = E z , a background scalar can be avoided and we find (turning off k a ) that the Klein-Gordon and Maxwell equations yield the dispersion relation Thus, instabilities can indeed be triggered by background EM fields. For p = 2 we find The dispersion relation for the scalar indicates that instabilities exist generically. Note, that only k 2 s enters the dispersion relation, and its sign is not fixed a priori. Depending on the sign of k 2 s there are indeed instabilities triggered by either the electric or the magnetic sector.

B. Homogeneous background axion field
The same analysis as above immediately shows that a constant background axion or scalar Ψ (and vanishing background EM field), triggers no instability. We thus turn to time-dependent background scalars. From the structure of the Klein-Gordon equation, the configurations of interest (dark halos, boson stars, superradiant clouds, etc) are indeed time dependent. The time dependence is harmonic and, in the non-relativistic approximation, set by the boson mass, Ψ ∼ e ±iµSt [14,35,36].

The instability: numerical results
FIG. 1. Time evolution of |α1(p, t)| rescaled by its initial value, for kapΨ0 = 0.5µS (upper panel) and kapΨ0 = 0.01µS (lower panel), and for selected frequencies. Notice that the p = µS/2 is generically excited and dominates the evolution at low couplings. Similar results hold for α2.
Consider therefore a uniform, coherent oscillating background axion state described by (this analysis closely parallels the discussion in Ref. [24]) where Ψ 0 (which we consider to be purely real, Im{Ψ 0 } = 0) determines the amplitude of the axion's oscillations, and complex conjugation is denoted with a * . Maxwell equations (we set k s = 0 here) can be analyzed using the following ansatz, where ω = |p| ≡ p and α µ , α * µ = α µ (p, t), α * µ (p, t) . When the coupling k a = 0, the solutions of Maxwell equations are a sum of plane waves, and α µ is timeindependent. There is no instability in such cases. It is the coupling to the axion or scalar that causes a possible time dependence for α µ . We note that there is U (1) symmetry gauge redundancy. To fix such gauge degree of freedom, we impose the Lorenz gauge condition Under this gauge condition, Maxwell equations are (latin letters stand for spatial indices) , where dots stand for derivatives with respect to time, and ijk is the totally anti-symmetric Levi-Civita symbol with xyz = 1. Let us denote the basis orthogonal to p by e (1) and e (2) . Then, the transverse components α (1) (t, p) and α (2) (t, p) obey the following system, Fig. 1 shows the time evolution of the strength |α (1) | of the vector potential for the following initial data, We consider both small and moderate coupling k a Ψ 0 . For small couplings, our numerical results show that there is an instability, α (I) ∼ e λt whose rate peaks at ω ∼ µ S /2, and for which (see Fig. 2 These results are compatible with previous predictions [24], and are not qualitatively affected by the initial conditions. On the other hand, for strong couplings, other (unstable) modes may also be excited.

FIG. 2.
Relation between the growth rate λ/µS and coupling kaΨ0 for the ω = 0.5µS mode. Red dots are numerical values, the dashed green line is the prediction from a small-coupling expansion, described in Appendix B 2 a, and which extends the lower order prediction in the main text.

Connection to the Mathieu equation
Let us go back to Eq. (19) and introduce the ansatz Now, we project the equations to a circular e ± , instead of a linear polarization basis, with ip × e ± = pe ± and y = y ω e ± . With a trivial time translation and re-scaling to a dimensionless time coordinate T = µ S t, the wavefunction y ω obeys the equation 1 In other words, we find that our problem is completely reduced to the well-known Mathieu equation! Thus, the complexity of having to deal with coupled (albeit ordinary) equations in the Lorenz gauge is no longer present. Mathieu equation predicts instabilities whenever ω 2 /µ 2 S = n 2 /4, n ∈ N and Ψ 0 k a p = 0 [37] or in terms of the relation between wavenumber and axion frequency, in agreement with the previous numerical results (see plot) for small values of the coupling. The Fouriertransformed vector potentials with wave numbers in the sequence (24) will dominate the other ones. One expects (based on the properties of Mathieu functions) that the first few of these modes will dominate and that -for small couplings -such instabilities will be significantly pronounced [38]. Perturbative investigation of unstable solutions demonstrate that the dominant rate of the instability is given by λ * /µ S = |Ψ 0 k a (ω * /µ S )| (see Appendix B for derivation and extension to higher orders), which for our problem reduces to Eq. (21), in accordance with previous work [24,39,40] and our numerical results. We thus recover in a considerably different way the same scaling for the instability rates of the coupled axion-EM system. In summary, a homogenous background of axions is an unstable configuration with a growth rate that scales linearly with the axion strength, for small couplings. We will find below in Section VI that this instability has a curved spacetime analog, even when the axion field is strongly inhomogeneous, and has the form of a condensate around spinning BHs.

C. Homogeneous background scalars
Let us apply the above analysis to scalar coupling case. Maxwell equations with scalar coupling are

Numerical results for the instability
We take again the ansatz (15) and (17) in the Lorenz gauge, Eq.18), to render the final equations sufficiently simple.
For p = 1, we find In order to extract the transverse mode, we introduce the projection operator P ij = δ ij − pipj p 2 . Then, the evolution equation for transverse mode is where α i = P ij α j , and θ is the phase of Ψ 0 . The time evolution of the vector potential α (I) is shown in Fig. 3 for the initial data (20). We find results similar to those of the axionic case: a homogeneous background scalar field and a vanishing EM field is an unstable configuration.
An instability is triggered whereby the vector grows exponentially, and ω = 0.5µ S seems to be -visibly at small couplings -the mode driving the state when p = 1.
The growth rates for the dominant ω = µ S /2 mode are shown in Fig. 4 for different coupling strength. Our results are consistent with for sufficiently small couplings. We can apply the same analysis to the p = 2 case. We find, The solutions of this equation are depicted in Fig. 3, again for the initial data (20). The dominant mode is now ω = µ S and for this mode the growth rates are

Analytical results
For p = 1 and small couplings k s , the substitution α = e iωt √ 1−ksΨ0 cos µ S t y ω allows us to re-write Eq. (26) as where again T = µ S t and we kept only the leading-order term in the coupling. Thus, we recover again Mathieu's equation. In this regime, we then easily get that the dominant mode is 2 and a growth rate λ = 0.25k s µ S Ψ 0 , in very good agreement with the numerics. y ω allows us to rewrite the relevant equation as where now T = 2µ S t. The same analysis therefore gives us the small-coupling dominant mode and the rate λ = Ψ 2 0 k 2 s µ S /4, in good agreement with the numerical results.

A. Perturbative framework
It is easy to see that Reissner-Nordström (RN) BHs are a solution to the field equations (2) when the scalar field vanishes and p ≥ 2 3 . We now wish to show that, as expected from the flat-space analysis of the previous section, this equilibrium solution of the field equations is unstable at large enough electric fields. We provide details in the axionic case. The scalar coupling follows through in the same way and we discuss the results below.
We consider the vector field A µ and massive scalar field Ψ propagating in a static, charged BH background and 2 These results are also derived in the Coulomb gauge in Appendix A 2, and we compute perturbatively the rate estimates in Appendix B. 3 For p = 1, RN is not a solution of the field equations, since the term Fµν F µν = 0 sources the scalar field equation (2a). This case is analogous to the more familiar Einstein-Maxwell-dilaton theories. We provide a perturbative solution in the next section.
coupled through the axionic coupling defined in the Lagrangian (1). The BH background is described by the metric where f (r) = 1 − 2M/r + Q 2 /r 2 , with M and Q the BH mass and charge, respectively, and the vector field satisfies A µ dx µ = Q/rdt. This spacetime has an event horizon and a Cauchy horizon located at r ± = M ± M 2 + Q 2 , respectively. We now linearize the field equations (2) around this background, and expand the perturbation functions in a complete basis of spherical harmonics. In particular we linearize the metric as and decompose the metric perturbations h µν in the Regge-Wheeler gauge: are the axial vector harmonics, and H l 0,1,2 , h l 0,1 , K l are functions of (t, r). Here a sum over the harmonic indices l and m is implicit.
The EM potential can be linearized in a similar way as: are again the axial vector harmonics, and u l (1,2,3,4) are functions of (t, r) and where again a sum over the harmonic indices l and m is implicit.
Finally, the scalar field is expanded as as In the following we omit the indices l and m, because in a spherical background different l and m modes completely decouple. This system can be studied in Fourier space by performing the decomposition: where δX(t, r) denotes schematically any perturbation function. The perturbation variables δX can be classified as "polar" or "axial" depending on their behavior under parity transformations. In a spherically symmetric background polar and axial perturbations always decouple. The background electric charge induces a coupling between gravitational and EM (grav-EM) perturbations, while due the pseudo-scalar nature of the axionic coupling, the scalar field only couples to the axial sector of the EM perturbations. Thus, only the axial sector is affected by the presence of the scalar field, while the polar sector is fully described by the usual grav-EM perturbations of the RN metric, which have been widely studied in the literature (see e.g. Ref. [41]).
Let us focus then on the axial sector. We define the Regge-Wheeler function through h 1 = rψ RW f (r) −1 and h 0 = i f (r) (rψ RW ) /ω and perform the linear transformations (from ψ RW , u 4 to Z ± ): with λ = 9M 2 + 4Q 2 (l − 1)(l + 2). After some algebra we find that axial perturbations can be described by a system of tree coupled ordinary differential equations, given by: where the "tortoise" coordinate r * is defined through the relation dr/dr * = f (r) and For spherically symmetric perturbations, l = 0, the scalar field decouples from the EM and gravitational perturbations, and one can conclude that the RN geometry is stable against radial perturbations. The same conclusion can be drawn when k a = 0 or Q = 0.
On the other hand, for l ≥ 1, in analogy to the flat space analysis of the previous section 4 , one expects the 4 As mentioned in the context of Eq. (9), the flat-space ax-existence of unstable modes above some critical value for the axion coupling constant k a . To find these modes we look for purely decaying modes at infinity 5 : while at the horizon, r = r + , we impose regular boundary conditions, i.e., purely ingoing waves described by where ψ 0 and Z ± 0 are constants. Since the system (40) -(42) is linear, one can fix the value of, e.g., ψ at the horizon and the problem becomes a 3-dimensional eigenvalue problem for ω and Z ± (r = r + ). We solved this system by using a 3-parameter shooting method: we shoot on the parameters ω and Z ± (r = r + ) and integrate the field equations starting from r = r + until the boundary conditions at infinity (44) are satisfied. Using the method outlined above, we find that the systems admits purely imaginary modes, ω = iω I , for which ionic instability around background electric field requires nonhomogeneous axion fluctuations. Thus, the instability only sets in for non-symmetric modes should not come as a surprise. 5 This system also allows for more generic boundary conditions describing quasibound states. For completeness we also discuss these modes in Appendix C.

Axionic couplings
ω I becomes positive above a critical value for k a . Given the ansatz (37) this therefore means that these modes grow exponentially in time and the system is unstable. An analogous instability was in fact shown to occur for asymptotically Anti-de Sitter RN BHs in Ref. [42]. For a given l, Q and µ S there exist an infinite family of solutions characterized by the number of nodes and k crit a . This tower of solutions signals the existence of new families of BH solutions, which we will discuss shortly. For the fundamental mode, i.e., the mode with no nodes, the results are summarized in Fig. (5) where we plot ω I as a function of k a for different values of Q, µ S and l.
We note that this instability is present even for massless scalars (M µ S = 0) as shown in Fig. 5. We also remark that it sets in for smaller couplings k a as we increase the BH charge and, hence, electric field. Fig. 5 furthermore indicates that the instability timescale does not depend monotonically on the multipole number l for sufficiently large coupling constants k a . For example, if the coupling is k a ≥ 50, Q = 0.1M and M µ S = 0 the l = 2 mode is the dominant one, i.e. it is the mode that grows faster.
For small BH charges Q M and small mass couplings M µ S 0.1, we find that the critical value above which the system becomes unstable for the nodeless modes is approximately given by The dependence of the threshold value of k a on the BH charge Q (or in other words, on the electric field close to the horizon) is in agreement with the results from the flat-space analysis (11). Near this critical value we find

Scalar couplings
The results for scalar couplings arise from a similar analysis. Focusing on the p = 2 case and on the spherically symmetric modes, we can decompose the scalar field as in (36) to reduce the analysis to a single ODE, For concreteness we show only the case with µ S = 0, although similar conclusions can be reached for generic values of µ S . We find that there is a critical value of k s beyond which an instability arises in a RN geometry. At this precise value of the coupling constant, k crit s , a new solution of the field equations is possible. In fact, we find that, for fixed M, Q there is an infinite set of constants k s for which a zero mode appears. This tower of solutions signals the existence of new families of BH solutions, analogous to the ones discussed in Refs. [43,44], which we will discuss at full nonlinear level shortly. We focus our attention on the smallest value of the coupling constant for which a zero mode exists, and for which the scalar has no nodes. For p = 2, our results are summarized in Fig. 6.
The threshold value is well described by k crit s = 2.4/Q, for small charge. This dependence on the charge (or, conversely, on the electric field) is the same as the one predicted from a flat-space analysis, (13). We should also stress that the scalar Ψ decays as 1/r at large distance, thus these BH solutions have a primary scalar hair and scalar charge, as we confirm in the next section when constructing these solutions at the full nonlinear level.

V. NEW BLACK HOLE SOLUTIONS
A. New BH solutions for axionic couplings

End-state of the RN instability
We showed that RN BHs are quite generically unstable when the EM field is coupled to a scalar field. The existence of a zero mode with frequency ω = 0 at the threshold of these instabilities suggests the existence of new BH solutions, branching-off from the RN solution for values of the coupling parameter at this threshold. As we discuss now, solutions with non-trivial scalar field indeed exist for the couplings considered here.
As shown in the previous section, for the axionic coupling the instability only exists for non-spherically symmetric axial perturbations, therefore due to the backreaction of these perturbations on the metric the end-state of this instability will most likely be an axisymmetric and stationary BH spacetime.
These solutions can then be obtained by solving a system of elliptic PDEs, which can only be solved numerically. We leave the construction of the full solutions for future work but we note that in asymptically AdS spacetimes an analogous instability and possible end-state solutions were constructed in Refs. [42,[45][46][47].
To gain insight into these solutions, we instead construct them perturbatively by making use of the fact that at k a = k crit a there should exist another stationary solution besides RN. Linearizing the field equations (2) and focusing on time-independent perturbations we find that, for l ≥ 2, we must have h 1 = 0, while h 0 , ψ and u 4 are given by the system of equations: where a prime denotes differentiation with respect to r. For l = 1 this can be simplified to: We solved these equations by imposing regular boundary conditions at the horizon and at infinity and using the shooting method outlined above. As expected, regular solutions only exist for non-zero BH electric charge and for k a = k crit a . Specific solutions are shown in Fig. 7. At spatial infinity the scalar field decays as Ψ ∝ 1/r l+1 , the vector field component decays as u 4 ∝ 1/r l while the metric component decays as h 0 ∝ 1/r l for l ≥ 2 and h 0 ∝ 1/r 2 for l = 1. Therefore these solutions have the peculiar property that frame-dragging effects are present even though their ADM angular momentum vanishes. The non-zero frame-dragging is instead due to the backreaction of the magnetic field induced by the axion on the metric. As a final remark, we should note that when obtaining equations (51) for l = 1 we have integrated out the solution that decays as h 0 ∝ 1/r and u 4 ∝ 1/r at spatial infinity. This additional solution, which exist for any k a , can be instead obtained from the system (50) with l = 1. For k a = 0, this solution is none other than the Kerr-Newman solution expanded up to first-order in the BH spin, while for k a = 0 this corresponds to a modified Kerr-Newman geometry as we discuss below.

Axions around rotating charged BHs
When the BH is spinning, one expects that the scalar acquires a non-trivial profile for any value of k a , since * F µν F µν = 0 for non-spherically symmetric spacetimes. In fact, it is easy to check that * F µν F µν = 0 when F µν = 0 or when the spacetime is spherically symmetric. Thus, the Kerr geometry together with F µν = Ψ = 0 is a full nonlinear solution of the theory above. In addition, for k a = 0, Kerr-Newman is also a solution. It is thus natural to look for rotating, charged BHs as perturbations of the Kerr-Newman geometry. At first order in k a 6 , we find that the only nontrivial correction appears at the level of the scalar, and to first order in rotation is described by (for µ S = 0 and k s = 0) 7 : Regular solutions exist for any nonzero value of the charge and can be found analytically after requiring regularity at the horizon and spatial infinity: These solutions decay as 1/r 2 at infinity and up to second-order in the BH charge are given by To summarize, when rotation is turned on, one expects two distinct families of solutions with axion hair. One family of solutions exists for any value of k a and thatto first order in the BH spin and in the axion coupling -is given by Eq. (52) with radial profile (54); a second family which is the rotating generalization of the solution discussed above, which we expect to branch-off the solution (52) at k a ∼ k crit a for small rotation.

An extended Wald solution
Aside from charged solutions, one should expect nontrivial axion configurations around non-spherically symmetric BHs with non-zero magnetic fields, since for this case one also has * F µν F µν = 0. In fact, it has been shown by Wald [48] that, neglecting backreaction, Kerr BHs immersed in a homogeneous magnetic field B aligned with the BH axis of symmetry allows for an exact analytical solution of Maxwell's equations: where k µ = (1, 0, 0, 0) and m µ = (0, 0, 0, 1) are the two Killing vectors that the Kerr metric admits. This field would lead to the BH accreting surrounding charge in the accretion disk and the interstellar medium. Therefore BHs would acquire a charge in those enviroments and be described by a Kerr-Newman spacetime, with a total vector potential given by [48] A with q = Q/M and Q is the accumulated BH charge. At equilibrium, the BH charge-to-mass-ratio is given by q = 2Ba. We can therefore analyse two different cases: (i) the BH is uncharged, and there is a net flow of charge from the surrounding medium, (ii) the BH is charged, but there is no net flow of charge from the surrounding medium. Recent estimates for supermassive BH in the Galactic center suggest that rotationally-induced charge is stable with respect to the discharging processes from the surroundings of an astrophysical plasma [49]. Let us then focus on the second (equilibrium) case and estimate the importance of the induced charge on the background geometry. Using the limit for a maximal astrophysically realistic magnetic field from Section II, we find q ≤ 10 −11 a/M , i.e. the geometry is still well described by the Kerr metric. Hence, we here consider a Kerr spacetime with the vector potential of the form 8 where F is a metric function given in Eq. (D2) and we refer the reader to Appendix E for details of the Wald solution.
Let us now consider, instead of Maxwell's equations, the generalized axionic equations (2). For k a = 0, Wald's solution is a solution to the problem, together with a vanishing scalar field. Thus, we are interested in a first-order (in k a B 2 M 2 ) production of axions, as a consequence of the EM background. The dominant term describing the axionic field is the equation where F µν denotes the Maxwell tensor corresponding to Wald's solution. Using Eq. (57) we find, to fifth order in the spinã = a/M , similar approach was recently considered in the context of pulsar magnetospheres [26]. 7 This equation can be obtained from the system (50) with l = 1 by considering the expansion parameter to be proportional to the dimensionless BH spin a/M . 8 In this subsection we work in Boyer-Lindquist coordinates (Appendix D) and suppress BL subscript.
One can now expand the left-hand side of Eq. (58) order by order in the spin, with Ψ = Φ 1ã + Φ 2ã 2 + . . .. To first order in rotation (and for µ S = 0) one gets and similar equations for higher order terms, each of which can be solved with an expansion in spherical harmonics. Finally, to first order in k a M 2 B 2 and fifth order in the spin for massless "axions" we find This field will in turn contribute to the background EM field, via Eq. (2b), but as a second order (in k a ) effect.
Let us now briefly comment on the scenario in which conditions for the superradiant instability are met (see Appendix F) and the conditions for the EM field instability (discussed in Section VI) are not satisfied. Then, Ψ (0) can be found by the expansion in α 2 g = (M µ S ) 2 [50]. The contribution of the magnetic field will appear, to first order in k a , at the fine-structure level (in α g ), since the dominant contribution of the driving term in Eq. (58) is aM/r 2 ∼ãα 4 g . Even though it is a subdominant effect, it could potentially be important for the consistent calculation of the level mixing in binary systems [50], for large magnetic fields and/or axion-photon coupling constant.
B. New BH solutions for scalar couplings

End-state of the RN instability
As already shown for the axionic couplings, the instability of the RN geometry discussed in Section IV suggests the existence of new BH solutions, branching-off from the RN solution 9 . We now construct those solutions for the scalar coupling described by the action (1) with k a = 0 but k s = 0.
Similar solutions were recently constructed for theories with a coupling of the form e −αΨ 2 F µν F µν [43], with α a coupling constant. In particular in Ref. [43] it was shown that the scalarized BH solutions are stable against 9 We note that the fact that the theory we consider has hairy BH solutions was already pointed out in Ref. [51], although they were not explicitly constructed. In Ref. [51] only BHs with magnetic charge were constructed.
spherically symmetric perturbations and found strong evidence that the solutions are the end-state of the instability by performing fully non-linear evolutions of the instability. In fact, when α 1 and k s 1 the coupling of Ref. [43] is equivalent to the quadratic coupling that we here study, hence all the conclusions should carry through. Here we only consider spherically symmetric spacetimes, and solutions for which the scalar field is nodeless. Solutions with nodes, and non-spherically symmetric (but static) solutions -which correspond to unstable polar modes with l > 0 -have also been constructed in Ref. [43].
The scalar field is given by Ψ ≡ ψ(r) while the most general spherically symmetric metric can be written as where N (r) = 1 − 2m(r)/r. The vector field is given by A µ dx µ = V (r)dt. After substitution in the field equations (2), we get (for generic values of p and µ S ): By imposing the existence and regularity of the solution across an event horizon at r = r H , in addition to regularity at infinity these equations can be solved using a standard shooting method. We refer the reader to Ref. [43] for more details. Our results for p = 2 and µ S = 0 are summarized in Figs. 8 and 9. Similar solutions can be constructed for p > 2 and µ S = 0. The left panel of Fig. 8 shows the dimensionless horizon area a H = A H /(16πM 2 ), with A H = 4πr 2 H , of the scalarized solutions as a function of the BH's charge-to-mass ratio Q/M . The right panel shows part of the domain of existence of the scalarized solutions. As can be seen, for a given k s , scalarized solutions only exist above a critical Q/M . The value at which these solutions branch-off the RN geometry agrees with the onset of the RN instability shown in Fig. 6. In agreement with Ref. [43] we find that solutions can be overcharged Q > M and for a given k s exist up to a maximum value of Q/M at which the numerics indicate that the horizon area vanishes.
In addition to the BH electric charge Q and the mass M these solutions have a scalar charge Q s . This can be seen in Fig. 9 where we show a specific scalarized solution. The scalar field can only be supported if Q = 0, i.e. Q s → 0 when Q → 0, in agreement with the expectation that a Schwarzschild BH is a stable solution of the field equations. Importantly, these results show that, for part of the parameter space, RN and scalarised solutions coexist with the same global charges.

Charged BH solutions for p = 1
For p = 1 RN is not a solution of the field equations (2). Charged BHs in this theory necessarily carry scalar hair for any value of k s since a non-vanishing F µν F µν will source the scalar field equation (2a). This case is analogous to the more familiar Einstein-Maxwell-dilaton theories. In fact, Einstein-Maxwell-dilaton theories have a coupling of the form e −2αΨ F µν F µν , with α as a coupling constant, and so are equivalent to the Lagrangian (1) with p = 1, k a = 0 and µ S = 0 when α 1 and k s 1. For Einstein-Maxwell-dilaton theory, closed exact analyt- ical BH solutions have been found [52]. For our specific coupling we were unable to find exact analytical solutions. However, since RN is a solution for k s = 0, perturbative solutions around RN can be found in an expansion in k s . For spherically symmetric solutions, to first order in k s and generic BH electric charge Q corrections to the RN solution appear only at the level of the scalar field and are given by solving the ODE (for µ S = 0 and k a = 0): Imposing regularity at the horizon and vanishing scalar field at infinity the solution is given by where r − = M − M 2 − Q 2 . One may easily find the first corrections to the metric and vector potential by expanding the solution around Q = 0. In particular up to second order in k s and fourth order Q the vector potential is given by: while the metric components are given by We note that the scalar field decays as 1/r at large distances and therefore these solutions carry scalar charge Q s given by Q s = k s Q 2 /(4M ).

VI. BURSTS OF LIGHT FROM SCALAR CLOUDS
Axion and axion-like light particles -even with negligible initial abundance -trigger superradiant instabilities around massive, spinning BHs [9][10][11][12][13][14]. The instability extracts rotational energy away from the spinning BH and deposits it into a cloud of scalars, with a spatial extent ∼ 1/(M µ 2 S ) [14]. Over long timescales, when the mass of the cloud is sufficiently large, GW-emission becomes important, and leads to a secular spin-down of the cloud (and BH), and a consequent cloud decay. Such systems are a promising source of GWs, both as resolvable and as a stochastic background, that can be detected with current and future detectors [15][16][17][18][19][20][21][22].
Our current understanding of the evolution of superradiant instabilities and accompanying GW emission neglects the coupling to matter, expected to be very weak. Arguments based on flat-space calculations similar to those worked out in Section III, suggest that when the axion strength exceeds a critical value (or in other words when the cloud extracts too much energy) an instability is triggered that might give rise to large amounts of EM radiation being emitted from BH systems [23]. Recently, these conjectures were shown to be true by some of us, through the evolution of Maxwell's field equations coupled to an axion field in a Kerr background. In particular it was shown that for critical values of the coupling k a Ψ 0 , EM fields are spontaneously excited in such environments, even at a classical level [25]. These instabilities can be indeed completely understood in the context of classical field theory, owing to the bosonic nature of axions and photons, that allows buildup of macroscopic numbers of particles.
Here, we will extend the results of Ref. [25] to scalar couplings and provide some analytical understanding of the mechanism. We refer the reader to Appendix G for the numerical formulation and initial data construction. We will always focus on a background axion or scalar field which is the product of the evolution of superradiant instabilities around spinning BHs. The growth of the scalar or axion due to superradiance is extremely slow to perform in a full nonlinear evolution; but see Refs. [13,53,54]. Thus, our setup is that of an axion fully grown by superradiance to some value, at which point we start monitoring the coupled system of Maxwell-Klein Gordon equations, in a fixed background geometry 10 following the approach of Ref. [57]. The initial data for the vector consists of a small azimuthal electric field of the form (see Appendix G 3 for further details) and The Θ profile is shown in Appendix G 3, and is used for completeness, since all our results are initial-data independent at the qualitative level.

A. Flat space analysis
Our main purpose is to show that EM instabilities may arise when axions or scalars exist, and they couple strongly to the Maxwell field. We will do this by order of complexity. Here, we artificially use a Minkowski background and we fix the scalar field to have the profile appropriate for clouds around BHs [17], Here, A 0 is amplitude of the field. For further details we refer the reader to Appendix F, and Refs. [17] and [25]. The purpose is to show that an instability exists even in this setup, but now with a critical threshold. In other words, the results worked out in Section III in a Our results are consistent with the existence of a critical coupling below which no instability is triggered, and well described by our analytical estimates in the small coupling regime.
Minkowski background generalize, except that the nonhomogeneous nature of the scalar or axion results in a critical coupling below which no instability is triggered. These results were reported for axions in our recent Letter [25].
We solved the (Maxwell) evolution equation of the EM field with a fixed scalar field for p = 1, 2 and mass couplings µ S M = 0.1, 0.2, 0.3, where M is the BH mass that supports the solution (75). Our results are summarized in Figs. 10 -11, and are consistent with the results we obtained for axion couplings recently [25].
The novel feature with respect to the homogeneous background scalar of Section III B is the existence of a critical coupling k s A 0 below which no instability occurs. This is apparent in Fig. 10 for both p = 1 and p = 2 (we stress that the results for axionic couplings can be found in Ref. [25]). The critical value is estimated below with simple analytical arguments. If the coupling k s A 0 is smaller than the threshold, the EM field dissipates away. This feature is induced by finite-size effects of the scalar cloud, as we argue below.
At large enough couplings, all initial conditions lead eventually to an instability (and exponential growth of the EM field), examples are shown in the bottom panels of Fig. 10. The growth rate depends very weakly on the initial data and on the coordinate at which the EM field is extracted. A closer inspection of instability rates allows us to estimate the critical coupling value. The rates are shown in Fig. 11 for different couplings, which also strongly indicates the existence of a critical threshold.

B. Kerr BHs
We have just discussed the exponential growth of an EM field around a "frozen" axionic or scalar cloud in a flat space background. These results suggest that when the effective coupling is larger than a threshold value, the EM field may grow exponentially -fed by the axionic cloud, which itself grew through super-radiance and extracted its energy from the spinning BHs. Here, we confirm this scenario with a fully numerical simulation around Kerr BHs (the geometry is kept fixed, but the coupled Maxwell-scalar system is evolved; following the approach of [57]). We refer the reader to Appendix D for notation on the Kerr metric representations, and summarize the formulation and initial data in Appendix G.
We solved the evolution equations for M µ S = 0.2, a = 0.5M (we also studied higher spins, the results are qualitatively the same), the results are summarized in Fig. 12 for scalar couplings with p = 1 and p = 2. As expected from the previous flat-space analysis, for small enough couplings any small EM disturbance dissipates away, and the profile of the axionic or scalar cloud is basically undisturbed. On the other hand, when the coupling is larger than a threshold, the EM field grows exponentially. As shown in Fig. 12, for large couplings an instability is indeed triggered. Because the instability acts to produce k ∼ µ S /2 vector fluctuations (for p = 1), at the nonlinear level these backreact on the scalar field, producing transient clumps of scalar field on these scales. This translates into an increase of the scalar, when observed sufficiently close to the BH, as seen in the upper panels of Fig. 12. On long timescales, the instability extracts energy from the scalar cloud and eventually lowers the effective coupling to sub-threshold values, leading to a now stable cloud. On even longer timescales, superradiance will grow the scalar to super-threshold values and the cycle begins again, as demonstrated in the axion scenario in Ref. [25] (see footnote 10). Similar effects have been found for scalar condensates with a self-interacting potential, but in the absence of couplings to the EM sector [58,59].
We would like to highlight a potential issue with the scalar couplings in general, and that clearly shows up when p = 1. When the effective coupling k s Ψ is of order unity, the kinetic term (left hand side of Eq. (2b)) for the vector field can vanish and the system becomes strongly coupled. The evolution in such case is ill-defined. In particular, we find for example that we cannot evolve E2 (see definition in caption of Fig. 12) in Fig. 12 past t = 520M , for this reason. It is possible that the dynamics of the gravity sector (neglected in this work) cure such anomalies, for example by producing BHs close to the threshold. Another possibility is that coupling to fermions will ensure that Schwinger-type creation works to prevent the EM field to ever approach such large values. The calculation of the time evolution near the strong coupling is beyond the purpose of our paper.

C. A simple analytic description of the results
Compared to the flat space analysis from Section III B, the (localized) axion "cloud" configuration introduces one more timescale in the problem, that of the time d needed for photons to leave the axion configuration, where d is a measure of the configuration size. Thus, there is another rate in the problem, If λ γ > λ * , with λ * being the estimate of the EM field instability rate for the homogeneous condensate, photons leave the configuration before the instability ensues and the effective rate of the instability is zero. In the other extreme, we can approximate the rate of the dominant FIG. 13. Instability rates for axionic couplings, in the presence of a background axion described by the cloud (75), for a Minkowski background. The analytical estimates for the instability rate λ for axionic couplings, as given to first order by Eq. (80) (dashed lines, full expansion is described in Appendix B 2 a), are compared with the numerical results of Ref. [25] (crosses). We find good agreement between our analytic estimates and numerical data for small mass (M µS ∼ 0.1) or small axion couplings.
instability by [40,60] λ where Ψ is some estimate of the average value of the axion field, to be implemented in the expression obtained for the homogeneous case (see Eq. (21)): Although the rate estimates were derived in flat space, the system under consideration is mostly in a weak field regime and we expect that they will provide a good description even in the context of instabilities around Kerr BHs. Let us therefore consider the dominant mode in the gravitational atom (75), frozen and embedded in Minkowski spacetime, and estimate the instability rate. For the measure of d we use the full-width-at-halfmaximum (FWHM) of the function (75) For the estimate of the field value we take the radial mean of the field on the FWHM and maximal contribution from the harmonic part of the function These estimates are in very good agreement with the results from the Letter [25] as can be seen in Fig. 13. In addition, a cutoff coupling below which no instability arises shows up naturally, explaining all the previous numerical results. In summary, a very simple and elegant analytic formula explains most of the results that we observe numerically.
The results above were worked out for the axionic coupling, but the underlying physics and mechanism remains the same for scalar-type couplings. Accordingly, we expect that the rate of the dominant instability is given by Eq. (77). Using Eq. (27), we find, for p = 1. Similarly, using Eq. (28) we find for p = 2, These estimates are shown together with numerical data in Figs. 11. Notice how such a simple estimate agrees very well with the full numerical evolution in the small coupling regime where the perturbative approximation is valid.

D. A note on plasma effects
FIG. 14. Time evolution of the massive scalar -massive vector field system around a Kerr BH with p = 2 and coupling k 2 s A 2 0 = 1.0, in which the initial data is an extended profile with (r0, w, E0) = (40M, 5M, 0.001). Here the scalar mass is µSM = 0.2 and the scalar cloud is evolving around a Kerr BH with a = 0.5M . The notation for the y−axis is explained in Appendix G 4.
Thus far, the system was assumed to evolve in a vacuum environment, when in reality the universe is filled with matter. We will approximate all this matter by plasma. The influence of plasma on axion-photon conversion has been discussed for superradiant axions [23,24], but also in other contexts [40,61]. EM wave propagation through plasma is described by the modified dispersion relation [62] where ν is the collision frequency and is the plasma frequency; m e , e and n e are the mass, charge and the concentration of the free electrons, respectively. Conceptually, it is helpful to consider two limiting cases -collisional (ω ν; appropriate in the context of plasma in the accretion disks) and collisionless (ω ν; in the context of the interstellar matter or a thin accretion disk).

Collisionless plasma
EM waves in this limit have a modified dispersion relation that is equivalent to providing a photon with a mass µ V = ω plasma . For high µ V ≥ (1/2)µ S , decay processes become kinematically forbidden 11 . For interstellar matter [61] ω plasma = n e 0.03cm −3 (6.4 × 10 −12 eV), the plasma frequency is below the range of the QCD axion mass and some of the ultra-light axions [see comments around Eq. (F7)]. Hence, one can expect EM instability not to be quenched on the primordial and lower range of the stellar BHs mass spectrum. We have numerically modeled this scenario as an scalar-Proca system. The time domain study is summarized in Fig. 14, and confirms this picture: when the (effective) mass of the vector field is larger than the axion mass, the burst is suppressed. We have checked that this suppression effect occurs in all our initial data and also for axionic couplings.
We will now provide some analytical control over these results by concisely reproducing some of the results of Ref. [24] and also expanding them to the scalar coupling case. Consider the dimensionless Mathieu equation, that describes parametric resonances of the EM field in the background of the homogeneous axion/scalar condensate (Section III B and III C) in the form where y is defined in Eq. (22) for axions, and in Section III C 2 for scalar couplings. Here, Υ = (p 2 + µ 2 V )/µ 2 S for axions and p = 1 scalars, and Υ = (p 2 + µ 2 V )/(2µ S ) 2 for p = 2 scalars. Also, = −Ψ 0 k a (p/µ S ), (1/2) (k s Ψ 0 /2) p for axions and scalars, respectively. As we showed already, the dominant instabilities occur for Υ = 1/4. It is 11 The decay process is a → γ + γ (for a Lagrangian with a ΨF 2 term), so if the photon has an (effective) mass µ V , in order for the decay to be energetically favourable, we should have µ S ≥ 2µ V .
well known that the critical stability curves on the Υ − (Ince-Strutt) diagram are given by [37] Inserting the appropriate Υ and , we can find the values of the parameters for which (87), a quadratic equation in p, has real solutions. First, consider the axion case where = (p). We find the critical plasma frequency in agreement with Ref. [24]. For scalar couplings we find Our results are in qualitative agreement with this prediction. One can also straightforwardly find corrections to the instability rate, induced by the effective mass. For example, in the axion case 12 The preceding analysis neglects the time-dependence of the plasma distribution, in particular it neglects also the backreaction by the cloud on the plasma. Although the full problem is outside the scope of this work, we note that the arguments of Ref. [40] suggest that nonharmonic time dependence would not jeopardize parametric resonances as long as the scalar mass is much smaller than the plasma frequency. However, the timeperiodic background of real scalars can drive matter resonantly in peculiar configurations oscillating with the (multiple of) scalar mass [36,63]. This behaviour would manifest itself as an additional harmonic term in the Mathieu equation, similar to the ones in Appendix B, and is expected to modify, but not eliminate the instability. The driving of plasma by the cloud could also deplete the plasma from the central regions [63], leading to a smaller effective photon mass and therefore to a more efficient instability.

Collisional plasma
Using estimates from Ref. [24] we find for the collision frequency and for the plasma frequency in the inner rims of the accretion disk around BHs. For BHs larger than M ∼ 10 −3 M , ω plasma > 10 −5 eV and the axion decay is forbidden in all of the parameter range interesting in a BH-superradiance context. However, one should also consider the geometry of the problem. Accretion disks are planar structures (when thin), immersed in a "spheroidal" scalar cloud. The EM field enhancement can originate in the space external to the accretion disk (there is a limitation from interstellar matter there, discussed in the previous subsection). Such waves can lead to Ohmic heating of the disk or disperse it through the radiation pressure. The quantitative analysis of this would probably depend on the geometry of the initial fluctuations. We should also note that the estimates of the peak luminosity from Ref. [23] (which are even lower than the ones estimated numerically in Ref. [25]; see below) indicate that the radiation pressure (if EM instability ensues) would blow away the surrounding plasma.
3. e + e − plasma Besides astrophysical plasma large electric fields can lead to Schwinger e + e − pair production. It was argued, at the order-of-magnitude level, that such plasma can indeed be created and reach large enough densities (and consequently critical ω plasma ) to block EM bursts [23]. Subsequently, e + e − annihilations would drive the plasma density down and restart the process again.
An adequate treatment of this phenomenon is beyond the scope of this work, but here we show that our numerical estimates are consistent with such high EM fields. We work with estimates from Ref. [25] for the axion coupling, µ S M = 0.2 and k a A 0 ≈ 0.3 − 0.4, where the peak luminosity was given by with M S representing the cloud mass (for the relation between A 0 , M and M S see Ref. [17]). The Maxwell invariant is where V ∼ r 3 ∼ 125 r 3 c (r c is the first Bohr radius, see Appendix F) is the cloud volume and t p ∼ 500M [25] is the time when the luminosity plateau develops. Taking |B| ∼ |E|, one finds where E c = 1.3 × 10 18 V/m is the critical Schwinger field.

VII. DISCUSSION
Extensions of the standard model of particle physics where a (ultralight) scalar or pseudoscalar Ψ couples to photons, are very "natural" to consider. These include terms of the form Ψ * F µν F µν or Ψ p F µν F µν , with p an integer. Such terms have been considered widely in the literature in the context of dark matter physics and cosmology, and in fact axionic-type couplings are elegant resolutions of the strong CP problem. Our purpose here was to discuss what looks like an important gap in our understanding: effects of axionic or scalar couplings in strong-gravity situations.
We have shown that such couplings lead to instabilities of homogeneous configurations and in fact seem to lead to violations of BH uniqueness results in General Relativity (see Ref. [64] for a discussion on these results and observational tests). In particular, for large enough couplings, charged RN BHs are unstable and spontaneously acquire a nontrivial scalar profile. In fact, unlike in Yang-Mills theory (where coloured BHs are unstable), it is plausible that the hairy solutions that we discussed at a perturbative level are stable (given that GR solutions are unstable, and a time evolution will in principle drive them to a hairy one). Thus, there are two possible solutions for a fixed set of physical parameters, despite the instability of one such solution.
Ultralight scalars or axions generically lead to superradiant instabilities of spinning BHs: on relatively short timescales the scalar field extracts a sizeable amount of the BH rotational energy and deposits it in a non axisymmetric "condensate" [14]. Since the scalar field is time-dependent, this massive condensate emits gravitational waves, a mechanism which can be used to place impressive constraints on the scalar mass [21,65,66]. In the presence of axionic or scalar-type couplings, a new channel onto which the system can radiate exists. Thus, one might rightfully worry that the superradiant instability responsible for the growth of the condensate might be quenched. What we showed, here and in a previous Letter [25] is that indeed such couplings will -if large enough and in the absence of plasma -give rise to periodic bursts of low-frequency EM radiation (but see footnote 10). These bursts deplete the condensate of a fraction of the scalar until superradiance becomes effective again. Thus, they effectively limit the maximum amount of energy that superradiance can extract before the EM burst. Our numerical results are supported by simple analytical estimates which cover in a concise and simple way all the phenomena we study, from flat to curved spacetime.
In practice, the interstellar plasma can influence the realization of the bursts. For collisionless plasma (interstellar medium and thin accretion disks), we find that the plasma would not affect significantly the existence of bursts for axions or scalars more massive than µ S ∼ 10 −12 eV.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for useful comments. The authors acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant "Matter and strong-field gravity: New frontiers in Einstein's theory" grant agreement no. where we defined We see that W p (t), i.e. the harmonic term that drives the instability, scales to leading order as k p s . To the lowest order in k s this equation reduces to the Mathieu equation for both p = 1: and the p = 2 case: Appendix B: Timescales for the instabilities of the flat-space homogeneous configurations

First order results
Here we derive the instability timescales, using a perturbative expansion in a small coupling parameter (k a , k s , Ψ 0 ), to be defined latter. We expect that for specific values of ω, given by Eqs. (24), (30) and (31), regular perturbation theory breaks down. However, one can start from the regular perturbation theory, see how instabilities buildup and use multi-scale [37] or dynamical renormalization group (DRG) [68] methods to regularize the problem. We will here use the later approach. As instabilities from a homogeneous axion field in flat space are exactly described by the Mathieu equation and for scalar couplings approximately, to leading order in k s , we will here derive the first-order result for the leading instability rate of the Mathieu equation.
We will consider Mathieu equation in the form of Eq. (86), with µ V = 0. Instabilities arise when Υ = n 2 /4, n ∈ N. We here focus on the dominant instability √ Υ = 1/2 and denoth the subscript of the function that governs dominant time-dependence of the vector potential with √ Υ. At zeroth order in , the solution is given by y 1/2 = Ae i(1/2)T + c.c. The differential equation for the first order correction is (B1) and the full solution is given by where T 0 is some arbitrary time where we imposed initial conditions. Higher-order terms will build up a secular terms of the form (T − T 0 ) m , where m is the order of the expansion, in the limit m → ∞ giving exponential growth. However, this behaviour invalidates our perturbative expansion. The DRG approach is based on the insight that the invalidation of the regular perturbation theory is a consequence of the big interval between T 0 and T [68,69]. In order to remedy this problem, we declare the parameters of the solution in Eq. (B2) as "bare" and rewrite them as the renormalized ones: Next, we expand T − T 0 = T − τ + τ − T 0 and choose a 1 ("counter-term") in such way to cancel secular terms ∝ (τ − T 0 ). The renormalized solution has the form . Arbitrariness of τ leads to the RG equation Working consistently at the 1 order and decomposing A(τ ) = X(τ ) + iY (τ ), we find i.e.
with ∂ τ X = Y . Finally, we choose τ = T as the "observational" time and conclude that the instability rate, to first order in is λ = , for (dimensionless) Mathieu equation. Substituting appropriate and rescaling back to the physical time [see definitions of √ Υ and below Eq. (86)], we obtain the estimates consistent with the numerical results in Eqs. (21), (27) and (28).

a. Axion coupling
Instabilities for the homogeneous configurations with the axion couplings are exactly described by the Mathieu equation. Hence, we will here obtain solution of the Mathieu equation to the second order, using DRG. Differential equation for the second order contribution is of the form 1/2 = −(2 cos T )y DRG procedure is analogous to the first order case. We will consider only leading order harmonics with T /2 as they will give dominant contribution to the instability. Note that the term of the form (T − τ )(τ − T 0 ) will selfconsistently cancel, sign of the renormalizability of the differential equation [69]. Second-order coefficient in Eq. (B3) is From the RG equation (B5) we obtain λ (R,2) a = √ 1 − 2 . As we should trust this solution to the order of O( 3 ), we perform Padé resummation of the results [37]. As the perturbative result λ (R,2) a is even function, we used the first non-trivial approximant (2, 1) and the final rate estimate is with appropriate defined below below Eq. (86). This results gives very good description of the numerical data in both Sections III B and VI, as shown on Fig. 2 and Fig. 13.

b. Scalar coupling
As will become clear later, we will first consider p = 2 case. Equation for the second order correction 15 is For p = 1 case, second order correction is governed by the equation For the rate estimate we obtain the same results as for the axion case. This result is clearly not a good description as the numerical results indicate (Fig. 4) that the function 15 We use Y 1/2 label for the higher order p = 2 corrections with 1/2 and Y 1/2 ≡ y 1/2 . For p = 1 we use J 1/2 mutatis mutandis. λ p=1 ( ) is divergent. Therefore, we go to the third-order contribution After Padé resummation at the order (1, 2) we obtain Besides the unstable modes discussed in the main text, massive fields around BHs can also form long-lived quasibound states (see e.g. Ref. [14] and references therein). Here, we discuss the effect that the axionic coupling has on these modes. In particular, due to this coupling, we expect that the long-lived scalar cloud will in turn trigger the excitation of a long-lived signal of EM waves. We consider the system of eqs. (40)- (42) and focus on the l = 1 case for simplicity. For this case, the field equation for Z − (42) completely decouples and to compute the modes of the system one can set Z − = 0. At the horizon we impose the usual regular boundary conditions (45). On the other hand one can check that asymptotically the most general solution behaves as: where k S = µ 2 S − ω 2 . If (ω) < µ S , one can then find regular modes by imposing B = D = 0. These modes are spatially localized states that slowly leak EM radiation to infinity due to the axionic coupling. On the other hand, if (ω) > µ S , the condition A = C = 0 yields purely outgoing waves at infinity, and allow us to compute the quasinormal modes of the system [70].
To find the quasi-bound state modes we employ two complementary methods. We used a direct extension of the two-parameter shooting method explained in the main text: fixing the scalar field at the horizon, we shoot for ω and the value of the vector at the horizon (45), such that the quasibound-state boundary conditions are met. We also used a direct integration method that allow us to reduce the problem to a one-parameter shooting method [66,[71][72][73], that we here describe. This method allows to compute modes of systems with an arbitrary number of coupled equations but for concreteness let us consider the system of coupled equations (40) and (41) for l = 1 with ingoing wave boundary condition at the horizon (45).
By imposing this boundary condition we obtain a family of solutions at infinity characterized by 2 parameters, corresponding to the 2-dimensional vector of the coefficients {ψ 0 , Z + 0 }. The quasibound-state spectrum can be computed by choosing a suitable orthogonal basis for the space of initial coefficients {ψ 0 , Z + 0 }. To do so we perform two integrations from the horizon to infinity and construct the 2 × 2 matrix where B and D are obtained from the boundary conditions at infinity (C1) and the superscripts denote a particular vector of the chosen basis, for example, B (1) corresponds to {ψ 0 , Z + 0 } = {1, 0} and B (2) corresponds to {ψ 0 , Z + 0 } = {0, 1}. The eigenfrequency ω = ω R + iω I will then correspond to the solutions of which in practice corresponds to minimizing det S m in the complex plane. The eigenfrequencies for M µ S = 0.4 as a function of k a are shown in Fig. 15. A generic conclusion of our analysis is that the coupling does not significantly affect the quasibound states, however the time-scale, 1/ω I , over which these states decay slightly increases the larger the coupling k a .
Appendix D: Kerr spacetime -coordinates and notation

Boyer-Lindquist coordinates
Kerr spacetime in Boyer-Lindquist coordinates (t BL , r, θ, ϕ BL ) is given by [74] where Σ =r 2 + a 2 cos 2 θ , and a = J/M , while J and M are the BH's angular momentum and mass, respectively. It is also useful to define the dimensionless rotational parameterã = a/M and note that 0 <ã < 1 (consequence of the Cosmic Censorship Conjecture). The Kerr spacetime admits two horizons and the angular velocity of the outer one is

Kerr-Schild coordinates
For numerical purposes it is often more useful to use Kerr-Schild coordinates. The two coordinates systems are related via [75] In these "spherical"-type Kerr-Schild coordinates the line element is Here we summarize the solution to Einstein-Maxwell theory derived by Wald [48] that represents a rotating BH encompassed by an originally uniform magnetic field. The EM field F in the background of a Kerr BH is written as [48] F =F 10 ω 1 ∧ ω 0 + F 13 ω 1 ∧ ω 3 where the orthonormal tetrad is and In BL coordinates the Maxwell tensor is given by 16 with components The corresponding 4-vector potential is given by 2aM r 1 + cos 2 θ , 0, 0, sin 2 θ F − 4a 2 M r .
(E6) 16 Note that we suppress subscripts BL to deflate the notation, and simply remind the reader that we refer to (t BL , r, θ, ϕ BL ).
The explicit expressions for the electric and magnetic fields are rather involved, but they asymptote to lim r→∞ E i =(0, 0, 0) , lim r→∞ B i = B 0 (cos θ, 0, 0) , (E7) at spatial infinity. In the limit of a non-rotating BH background the fields reduce to where α is the lapse function. The invariants of the EM field are The latter vanishes in limit of zero spin.

Appendix F: Gravitational atom
Superradiant instability around Kerr BHs happens as long as the superradiant condition is satisfied [14] ω m < Ω + , where ω is a bosonic wave frequency and m is a spherical harmonic "quantum" number. Since the field can "leak" through the horizon, solutions of the Klein-Gordon equation are quasi-bound and the field frequency is a complex number.
In the weak-field regime (scenario dubbed as a gravitational atom) velocities and densities are small, as the cloud spreads over a large volume around the BH [17]. Particle dispersion relation is (for now we omit the imaginary part) ω 2 = µ 2 + p 2 and in the weak-field regime we expect ω ∼ µ, so that, up to the second order in wave number p, Dimensionless expansion parameter here is the group velocity of the field v = p/µ. In order to understand long-range behaviour, note that the field is trapped by BH gravity (p 2 < 0), so that we expect exponential tail Ψ ∼ e ipr ∼ e −|p|r . Length scale associated with the particles momentum is De Broglie wavelength λ D = 2π/|k| and in the nearhorizon limit important scale is the gravitational radius r g = 2M . We will use virial theorem 17 (2T ∼ V ) to understand dependence of the typical size of the cloud r c on µ and M : De Broglie wavelength of the wave on radius r c depends on the number of modes excited as nλ D = 2r c π. We find where α g is the fine structure constant and λ c = 1/µ is the (reduced) Compton wavelength. Finally, we see that the behaviour of the real part of the spectrum is the same as for the Hydrogen atom, mutatis mutandis: This equation describes leading-order behaviour of the real part of the frequency. For higher-order corrections see Ref. [50]. Imaginary part (decay width) was analytically calculated in the weak-field limit in Ref. [9] and it implies that the dominant growth mode is |2 1 1 , where the field is described by Eq. (75). In this state the density peak is located at r = 5r c .

Appendix G: Formulation as Cauchy problem
In this Appendix, we summarize the Cauchy problem for our system, which we use in Section VI.

3+1 decomposition
The equations of motion of this system are given by Eq. 2. In this work, we ignore the dynamics of gravity, and solve the Klein-Goldon equation Eq.2a, and Maxwell's equations coupled to a scalar field according to Eq.2b. In order to calculate the time evolution of this system, we apply the 3+1 decomposition to these equations with Lorenz gauge. In this decomposition, the metric function is written as where α is a lapse function, β i is a shift vector, and γ ij is a spatial metric.
The vector potential A µ is written as where A φ = −n µ A µ , and A i = γ µ i A µ . By substituting these expressions in the field equations, we get the following evolution equations for this system: where Π is the momentum conjugate of the scalar field, E i and B i are the electric field and magnetic field, defined as E i = γ i µ F µν n ν , B i = γ i µ * F µν n ν = − ikl D k A l ; (G8) A i is a spatial component of vector field; D i is a covariant derivative with respect to γ ij ; K is the trace of extrinsic curvature; Z is a auxiliary field, which is introduced to stabilize the constraint damping mode; ijk = − 1 √ γ E ijk and E ijk is the totally anti-symmetric Levi-Civita symbol with E xyz = 1. The above equations of motion include both the axion and the scalar coupling (note that they are implemented mutually exclusively) as well as a massive photon (set at µ V = 0, except in Section VI D). We also get the following constraint equations, Our numerical time-evolution code is based on this formalism, expressed in Kerr-Schild type coordinates discussed below,

Background spacetime
For the background spacetime, we use Kerr-Schild coordinates (see Appendix D 2). But, in order to avoid the coordinate singularity, we use Cartesian type coordinates, which are defined by the following coordinate transformations: x = r cos ϕ KS sin θ − a sin ϕ KS sin θ , (G10) y = r sin ϕ KS sin θ + a cos ϕ KS sin θ , (G11) z = r cos θ . (G12) In these coordinates, the metric can be written as where H and l µ are defined as l µ = 1, rx + ay r 2 + a 2 , −ax + ry r 2 + a 2 , z r .

Initial Data
In order to construct the initial data, one must solve the constraint equations (G9). For the scalar field, we use following a simple axion cloud profile as initial data: Ψ(t, r, θ, φ) = A 0 rM µ 2 e −rM µ 2 /2 cos(φ − ω R t) sin θ , (G20) where A 0 is the amplitude of the cloud, and ω R µ is a bounded-state frequency.
For axionic couplings, the constraint equation is given by The initial data that we use is given by where E 0 (r, θ) is an arbitrary function of r and θ. One can show that this profile satisfies the constraint equations. For E 0 (r, θ) we use a Gaussian profile 18 .
where E 0 , r 0 , and w are the amplitude, the peak radius and the width of the initial electric field. For scalar couplings, the constraint equation is the following, As the initial profile, we use the following solution of the constraint equations, where F (r, θ) is an arbitrary function of r and θ. In this work, we use where E 0 , r 0 and w are constants, which characterize the strength, the radius, and width of the Gaussian profile of the electric field. Θ characterizes the θ-dependence of the profile. In our study, we use two different profiles for F (r, θ). The first is a simple constant value, which we call "extended profile" since it is directionindependent. The second profile we use is Θ(θ) = sin 4 4θ for 0 < θ < π 4 0 for π 4 < θ < π. (G30) We term this a "localized profile" since it is sharply peaked along some directions only.

Numerical code
FIG. 16. Time evolution of the norm of constraint. This shows second-order convergence, because of the integration scheme used to compute the norm.
We developed a numerical code which solves the time evolution problem for the scalar and EM fields under the above formulation. Our numerical code is written in C++; we adopt a 4th order time-integration Runge-Kutta method, and 4th order discretization for the spatial direction. To calculate long term simulation, we use fixed mesh refinement. The grid structure is layered around BHs. Furthermore, to avoid the physical singularity inside the horizon, the metric in the horizon is replaced with a smooth regular function. The numerical domain of our simulation is −600M ≤ x ≤ 600M , −600M ≤ y ≤ 600M , and −600M ≤ z ≤ 600M . Refinement level is typically 4, and the ratio of resolution between adjacent refinement level is 2. To avoid high frequency modes which comes from boundary between adjacent refinement region, a Kreiss-Oliger artificial dissipation is added.
As a test simulation, we calculate the time evolution of extended initial data whose (r 0 , w, E 0 ) = (40M, 5M, 0.001), with p = 1 scalar coupling, and a = 0.5M . The time evolution of the norm of the constraint is depicted in Fig.16. This evolution shows 2nd-order convergence, simply because we used a second-orderaccurate integration scheme to compute the norm.