Spinning black holes with axion hair

In this work we construct and analyze non-perturbative stationary and axially-symmetric black hole solutions in general relativity coupled to an electromagnetic and an axion field. The axion field is coupled to the electromagnetic field, which leads to hairy solutions in the presence of an electric charge and rotation. We investigate the existence and characteristics of these solutions for different values of the spin, charge and coupling constant. Our analysis shows the presence of violations of the Kerr–Newman bound, solutions with large positive and negative values of the gyromagnetic ratio, and the existence of multiple branches of solutions with distinct properties, demonstrating that black hole uniqueness does not hold in this scenario. The code used in this study is publicly available, providing a valuable tool for further research on this model.


Introduction
The existence of dark matter, which is thought to make up around 85% of the universe's mass, is one of the greatest mysteries in modern astrophysics.A promising dark-matter candidate are axions, first proposed by Peccei and Quinn as a potential solution to the strong CP problem [1].They can be extremely light and weakly interacting [2].Additionally, axion-like fields have been observed to arise naturally in string theory constructions [3,4], further adding to the motivation for studying their potential consequences.
In the strong gravity regime, massive scalar fields, such as axions, are expected to trigger superradiant instabilities around spinning black holes [5][6][7][8][9][10][11].This causes the black hole to lose some of its rotational energy, which is then transferred to a dense cloud of axions.Over time, the cloud will emit gravitational waves and slowly lose energy, causing it to shrink and eventually disappear.These gravitational waves, as well as a plethora of other unique signatures, can be detected by current and future instruments [12][13][14][15][16][17][18][19][20][21], making these systems an interesting subject for study.Moreover, it has been proposed that when photons interact with superradiant axion clouds, they can produce sudden, intense bursts of light that are similar to the emissions from a laser at a quantum level [22][23][24][25].
No-hair theorems state that electrovacuum, stationary black holes in Einstein-Maxwell theory are described by the Kerr-Newman solution [26] (see also Refs.[27,28] for reviews).However, if couplings are allowed between the scalar and electromagnetic fields, black holes can have additional properties, known as "hair".These hairy black hole solutions have long been studied in contexts such as Kaluza-Klein theory and supergravity [29,30].More recently, couplings of this type and respective black hole solutions have been explored in the context of black hole spontaneous scalarization [25,[31][32][33][34][35] (see Ref. [36] for a recent review).
Previous studies have investigated the effects of axions on magnetically charged static and spherically symmetric black holes 5 [34,37,38] and perturbative solutions away from the slowly-rotating Kerr-Newman black hole [25,39].In this work, we go beyond the perturbative/magnetically charged framework and construct fully non-linear stationary and axiallysymmetric solutions of electrically charged and rotating black hole spacetimes with axionic hair.Unlike uncharged and/or spherically symmetric cases without magnetic charge, the axion is sourced by a non-zero electromagnetic field in this scenario, resulting in equilibrium hairy configurations.
The structure of this paper is as follows.In Section 2, we present the model, along with the equations of motion and the ansatz used for solving them.We also explain the boundary conditions applied and our numerical methodology.Section 3 contains the main outcomes of our study.Subsection 3.1 deals with the solutions obtained for small spin values, while subsection 3.2 presents the results for higher spin values.We summarize our findings in Section 4. The source code utilized to construct the solutions examined in this work is accessible to the public [40], enabling further investigations into this model.

The model and setup
We will explore Einstein's gravity coupled to one of the best-motivated extensions of the Standard Model, which is described by the following action where totally anti-symmetric Levi-Civita symbol, and the axion a is a real pseudoscalar field with mass m a that couples to electromagnetism.The axionic coupling constant g aγγ , has dimensions of inverse mass.The currently allowed parameter space for the axion model described by action (1) is shown in Fig. 1, in the (m a , g aγγ ) plane.In the absence of a mass term, the axion enjoys the shift-symmetry a → a + c, for constant c.
Varying the action (1) with respect to the metric, we obtain the field equations where the total stress-energy tensor is given by  The figure shows the constraints on the axion coupling constant g aγγ and mass m a from various experiments.The QCD axion is represented by the yellow band.The magenta band represents constraints from collider physics experiments, the blue from dark matter experiments, the green from stellar experiments, and the gray from astrophysics.The dashed lines indicate projections from future observations.Adapted from Ref. [41].
which is unaffected by the axionic coupling.The axion equation of motion is and the Maxwell-axion equations are From now on, for convenience, we work in units of √ 2M Pl by performing the rescalings which amount to replacing 1/M 2 Pl → 2 in Eq. ( 2), while Eqs.( 4) and (5) remain unchanged.Furthermore, in this work we consider only a massless axion for simplicity.However, we have performed initial simulations for a massive axion and found that the results are qualitatively similar.
We aim to study black hole solutions to the equations of motion (2), ( 4), (5), that are regular on and outside of the event horizon, stationary and axially-symmetric.These solutions possess two commuting Killing vector fields, ξ = ∂ t and η = ∂ ϕ , in an adapted coordinate system.As such, we consider the following metric ansatz in quasi-isotropic coordinates: where and where f , g, h, W are dimensionless functions of r and θ, and r H is the coordinate location of the event horizon.This ansatz for the metric is motivated by discussions such as those in Ref. [42], and the Kerr-Newman solution can be written in this form (see Appendix A).For the electromagnetic four-potential, we consider the following ansatz: where A t and A ϕ depend on both r and θ, as does the axion a.To solve the system, we follow Ref. [43] and use the following combination of field equations (2) (written here in the form E = 0): together with the axion equation ( 4) and the two non-trivial components of the Maxwell equations (5).It is worth noting that the third equation in Eq. ( 9) is equivalent to an equation that, for a massless axion (m a = 0), is exactly the same as in vacuum General Relativity and depends only on the metric function g.Using an ansatz of the form g = g r (r)g θ (θ), the equation is separable and, by imposing regularity, we find that the angular part is trivial, leaving , as is the case for a Kerr-Newman black hole.Although we know a closed-form solution for the metric function g, we will not use it in the numerical method and will instead use it as a test for the code.
To solve the problem numerically we compactify our radial coordinate as . We impose the following boundary conditions.Axial symmetry and regularity on the axis of symmetry θ = 0, π, impose The involved functions possess definite parity with respect to θ = π/2 (all are even parity, except a which has odd parity), and therefore we need only consider the range θ ∈ [0, π/2].Parity considerations then imply Asymptotic flatness is imposed by the boundary conditions for x = 1, where are the dimensionless spin and charge, respectively, with M , J, and Q being the ADM mass, angular momentum, and electric charge respectively.These can be extracted from the asymptotic behaviour discussed below in Eq. ( 16).At the horizon we impose for x = −1, where we have used the gauge freedom of the electromagnetic field to impose the condition on A t .To avoid conical singularities we also impose h = 1 on the symmetry axis, as explained in Ref. [43].
Asymptotically we have the fall-offs where M , J and Q were defined above, D is the dipole moment of the axion, Φ is the electrostatic potential difference between infinity and the horizon, and µ M is the magnetic dipole moment.
The gyromagnetic ratio g (also known as g-factor ) measures how the magnetic dipole moment is induced by the total angular momentum and charge for a given mass For a Kerr-Newman black hole, g is always equal to 2. The Hawking temperature T H [44] and event horizon area A H are given by with the entropy of the black hole given by S = A H /4. The angular velocity of the horizon is To check the accuracy of our solutions we use a Smarr-type relation [45,46] that black hole solutions should obey which relates global charges at infinity and horizon quantities.
To solve the partial differential equations resulting from the field equations, we will use the code described in Ref. [43], which employs a pseudospectral method together with a Newton-Raphson root-finding algorithm to solve the non-linear system.We expand each of the functions in a spectral series with resolution N x and N θ in the radial and angular coordinates x and θ, respectively.The spectral series we use for each of the functions f, g, h, W, A t , A ϕ (collectively denoted by F) is given by  24) and our numerical results for the axion for a solution with χ = 0.001, g aγγ = 0.5, and q = 0.2.The initial guess provided to the code to obtain this solution was a Kerr-Newman black hole with the same χ and q, and a vanishing axion profile.
where T i (x) denotes the i th Chebyshev polynomial.For the axion, due to its odd parity properties, we use With these spectral expansions, all the angular boundary conditions are automatically satisfied and need not be imposed explicitly in the numerical method.
In our setup, we have five input parameters: (r H , χ, q, g aγγ , m a ) and typically use a resolution of N x × N θ = 40 × 8.For the vast majority of solutions we fix r H = 1 in the code.The error on our solutions is typically less than O (10 −12 ), as estimated by the Smarr-type relation (20), although errors increase for large charges and/or spins.As an initial guess for the solver, we usually use a comparable Kerr-Newman black hole or a previously obtained and similar black hole with axion hair.We focus on studying fundamental solutions (solutions without nodes in the scalar field radial and/or angular profiles), as these are expected to be the (most) stable solutions.

Small values of spin
Ref. [25] studied the stability of Reissner-Nordström black holes for the theory (1) and concluded that for small electric charge (q 1) and axion mass (M m a 1), the system is unstable against axial perturbations with multipole l if q q crit , where We therefore expect at least two branches of solutions to exist: one that exists for any value of q (labeled as branch 1) and another branch that exists for values q q crit (labeled as branch 2)6 .
χ=0.001, q =0.1, θ=π/4 A φ KN ×100 The first branch's solutions for a massless axion were also constructed perturbatively in Ref. [25] to first order in spin χ and axion coupling g aγγ around a Kerr-Newman black hole.To this order, non-trivial corrections appear only at the level of the axion and they read where r BL is the Boyer-Lindquist radial coordinate, related to the one in our coordinate system via Eq.(A.7), and r ± BL = M 1 ± 1 − q 2 − χ 2 .In this section, our focus is on constructing non-perturbative solutions and examining the existence of two solution branches.Specifically, we will concentrate on the small spin limit, where we take χ = 0.001 to ensure that Eq. ( 23) holds.As another test to the code we have compared our numerical results with the perturbative solution in Eq. ( 24), observing great agreement as shown in Fig. 2.
By taking q = 0.1, we can observe the presence of the two solution branches.The first branch exists for all values of g aγγ , while the second branch is only attainable for values of g aγγ 49, in agreement with Eq. ( 23).In all of the numerical solutions we obtained, we found that the metric function g(r, θ) matched the analytical results from the previous section with high accuracy.However, the function A ϕ , which is related to the magnetic field of the system, is significantly amplified compared to a comparable Kerr-Newman black hole.The two branches of solutions exhibit drastic differences in the profiles of A ϕ and a, having different signs.This has implications for the gyromagnetic ratio of the solutions defined in Eq. ( 17).An example of this can be seen in Fig. 3, where the branch 1 solution (represented by the red line) has g ≈ 575.02 and the branch 2 solution (represented by the blue line) has a negative gyromagnetic ratio of g ≈ −44.93 (recall that for a Kerr-Newman black hole, g = 2 always).
To gain a deeper understanding of the solution branches and their properties, we traced each branch through its domain of existence while monitoring various properties, including the gyromagnetic ratio, entropy, Hawking temperature, and the Gauss-Bonnet curvature scalar on the horizon.Specifically, we looked at the g aγγ = 50 case with χ = 0.001, and varied the parameter q.Based on Eq. ( 23), we anticipated that the second branch of solutions would emerge for values q 0.098.Fig. 4 displays the results for the gyromagnetic ratio, where the emergence of the second branch of solutions is manifest.We observe that the first branch starts at values close to g ≈ 2 for small values of q, similar to a Kerr-Newman black hole.However, once q increases, we see a significant departure from this behavior, with a sharp The figure shows the behavior of the gyromagnetic ratio g for the two branches of solutions.
The left panel displays a sharp increase in the gyromagnetic ratio for first branch solutions around q ≈ 0.098, coinciding with the appearance of a second branch of solutions.A more detailed view of this feature is provided in the right panel, which is a zoomed-in version of the plot on the left.These solutions were obtained with χ = 0.001 and g aγγ = 50. ( (2) increase occurring at around q ≈ 0.098 when the second branch of solutions emerges.This second branch of solutions has negative values of g, with very large magnitudes observed for q ≈ 0.098, and drives the gyromagnetic ratio to small values as q becomes larger, until the branch of solutions ends.However, we have noticed some features in the gyromagnetic ratio of the second branch for certain discrete values of the charge q.Although these features are too small to be visible in a plot for the small spin case, we will discuss them in more detail in the next section for higher spins.
In Figure 5, we can see the plots of entropy and Hawking temperature for solutions with small spin χ = 0.001.It is observed that the second branch of solutions behaves similarly to a Kerr-Newman black hole, whereas the first branch differs significantly.In particular, the first branch does not approach an extremal black hole as we increase q, because the Hawking temperature does not approach zero.Instead, the numerical monitoring of the Gauss-Bonnet curvature scalar on the horizon suggests the presence of a curvature singularity on the horizon as we approach the maximum value of q.Additionally, in the domain where both branches of solutions co-exist, the first branch is always entropically preferred over a comparable second branch solution.Moreover, it is possible for first branch solutions to violate the Kerr-Newman bound, which states that q should be less than or equal to 1 − χ 2 , which, for small spins, is approximately q 1. Fig. 6 presents the radial profiles of the involved functions and fields for two solutions belonging to the first branch that violate the Kerr-Newman bound, one of which is very close to the singular solution.It should be noted that while the analytical examples presented in Ref. [47] and Ref. [48] for related theories feature endpoint solutions that are singular, those solutions have a vanishing horizon area, which is not observed in our solutions.
We have observed another intriguing characteristic of the first branch solutions, namely the existence of counterrotating black holes where the horizon angular velocity Ω H and the total angular momentum J have opposite signs.This can be seen in Fig. 7, where the black hole solutions initially have a positive Ω H for small values of q, but undergo a transition from positive to negative values of Ω H at around q = q crit .This feature is also observed in the solutions of Fig. 6.Counterrotating black holes have been previously noted in works such as Refs.[49][50][51].
Fig. 8 shows that the maximum charge q max allowed by the first branch solutions strongly depends on the coupling g aγγ .Specifically, we observe that the violations of the Kerr-Newman The plot shows the maximum charge q max allowed by the first branch solutions as a function of the coupling constant g aγγ , with a fixed spin value of χ = 0.001.The maximum value of q max is approximately 1.525, which is obtained for a value of g aγγ around 25.It can be observed that q max experiences a significant growth around g aγγ ≈ 22.
bound are only possible for values of g aγγ less than or around 80, for a fixed spin value of χ = 0.001.For higher values of the coupling, charge to mass ratios are always less than unity.

Higher spins
In this section, we examine axionic solutions with higher spin values χ.The full domain of existence for first branch solutions with g aγγ = 50 (g aγγ = 100) is presented in Fig. 9 on the left, shaded in red (blue), respectively.These solutions can violate the Kerr-Newman bound (dashed line) throughout their domain, with values as high as q ≈ 1.31 observed for spins χ ≈ 0.2 in the g aγγ = 50 case.Approaching the red and blue lines, we observed a divergent behavior of the Gauss-Bonnet scalar on the horizon.As in the case of small spins, the maximum allowed value of charge q depends heavily on the coupling g aγγ , as demonstrated by comparing the g aγγ = 50 and g aγγ = 100 cases.In a manner similar to the case of small spins, we have also observed the presence of solutions that do not belong to the first branch.As we discussed earlier, for a fixed value of coupling, these solutions exist only above a critical value of the charge q.While in the case of small spins, we can predict the critical charge value from Eq. ( 23), numerical solutions are needed to determine the bifurcation points for higher spins.The existence line for these solutions for g aγγ = 50 was computed up to spins χ ≈ 0.96 and is presented in the right panel of Fig. 9.It is noteworthy that for higher spins, the critical value of the charge increases.In the case of χ = 0.35, at least four distinct branches of solutions were observed, including the first branch, whose domain of existence was presented in Fig. 9, as well as three others with markedly different properties, as shown in Fig. 10.The second branch corresponds to the second branch χ=0.35, g aγγ =50

Branch
(1)   seen in the small spin limit, but with the added observation of a resonance phenomenon on the gyromagnetic ratio of the solutions at certain discrete values of the charge q.In the specific case shown in Fig. 10, resonances were observed at q ≈ (0.141, 0.230, 0.396, 0.568, 0.710, 0.804).
A third new branch was found to bifurcate from the second branch and did not violate the Kerr-Newman bound, while a fourth new branch began at the existence line and did violate it.In the domain of co-existence, for the same charge and spin, the entropy of the solutions in the branches (normalized by 4πM 2 ) was observed to be ordered as S 1 > S 4 > S 3 > S 2 .All solutions we studied exhibit ergoregions, light rings, and innermost stable circular orbits that have the same topology as those of a comparable Kerr-Newman black hole.

Conclusions
In this study, we investigated hairy black hole solutions in General Relativity minimally coupled to electromagnetic and axion fields.When electric charge and rotation is considered, the coupling of the axion to the electromagnetic field results in hairy solutions due to F µν * F µν = 0, which sources the scalar field equation ( 4).We solved the field equations of the theory (1) numerically, using a pseudospectral method [43].
For small spin values, we verified the existence of two branches of solutions predicted in Ref. [25].The critical charge value at which a second branch of solutions appears agreed remarkably well with their prediction.The two branches of solutions have opposite signs for the magnetic function A ϕ and the axion, leading to gyromagnetic ratios with opposite signs.The first branch of solutions was always found to be preferred based on entropic considerations, and violations of the Kerr-Newman bound were observed for these solutions.Additionally, counterrotating black holes, where the angular velocity at the horizon and total angular momentum have opposite signs, were observed.
For higher spin values, we found at least four distinct branches of solutions, all with unique properties.Again, the first branch was preferred based on entropic considerations and was present for all charge/coupling values.We computed the domain of existence for first branch solutions for two distinct coupling values and the existence line at which solutions for other branches start to emerge.
Avenues for future research include a deeper study of the massive axion case, and of phenomenological properties of these solutions, such as ergoregions, light rings, innermost stable circular orbits, and shadows.It would also be interesting to investigate the polarization effects of light rays passing through our axion profiles [52][53][54].
Our study provides the first non-perturbative numerical exploration of spinning hairy black hole solutions in the axion model presented in Eq. ( 1).These solutions exhibit intriguing properties that are different from those of typical Kerr-Newman black holes.Our publicly available code [40] can be used to further investigate these solutions in future studies.

Figure 1 :
Figure 1:The figure shows the constraints on the axion coupling constant g aγγ and mass m a from various experiments.The QCD axion is represented by the yellow band.The magenta band represents constraints from collider physics experiments, the blue from dark matter experiments, the green from stellar experiments, and the gray from astrophysics.The dashed lines indicate projections from future observations.Adapted from Ref.[41].

Figure 2 :
Figure2: Comparison between the perturbative solution in Eq. (24) and our numerical results for the axion for a solution with χ = 0.001, g aγγ = 0.5, and q = 0.2.The initial guess provided to the code to obtain this solution was a Kerr-Newman black hole with the same χ and q, and a vanishing axion profile.

Figure 3 :
Figure3: The radial profiles of the function A ϕ and field a are shown for two solutions with g aγγ = 50, χ = 0.001, and q = 0.1, both plotted for an illustrative value of θ = π/4.The solution belonging to branch 1 is labeled (1), while the solution belonging to branch 2 is labeled(2).The profiles of A ϕ and a for the two solutions are markedly different.

Figure 4 :
Figure 4: The figure shows the behavior of the gyromagnetic ratio g for the two branches of solutions.The left panel displays a sharp increase in the gyromagnetic ratio for first branch solutions around q ≈ 0.098, coinciding with the appearance of a second branch of solutions.A more detailed view of this feature is provided in the right panel, which is a zoomed-in version of the plot on the left.These solutions were obtained with χ = 0.001 and g aγγ = 50.

Figure 5 :
Figure 5: The left plot shows the behavior of the entropy and the right plot shows the behavior of the Hawking temperature for both branches of solutions, with a dashed gray line indicating the profile for a Kerr-Newman black hole.The values of χ = 0.001 and g aγγ = 50 were used to obtain these solutions.

Figure 6 :
Figure 6: Radial profiles of two solutions belonging to the first branch that violate the Kerr-Newman bound, one of which is in close proximity to the singular solution.The value of θ is fixed at π/8.

Figure 7 :Figure 8 :
Figure7: Angular velocity of the horizon as a function of the charge to mass ratio q for solutions from the first branch.The vertical dashed lines indicate the critical values q crit obtained using Eq.(23) for each coupling constant.It is evident that the horizon angular velocity undergoes a transition from positive to negative values around q = q crit .

1-χ 2 -q 2 =0Figure 9 :
Figure 9: (Left) The region in parameter space where solutions of the first branch exist is shown for g aγγ = 50 and g aγγ = 100.(Right) Solutions outside the first branch exist for values of the charge to mass ratio q greater than those shown in the pink line for a coupling g aγγ = 50.The dot-dashed gray line indicates the value predicted by Eq. (23) for small spin, showing excellent agreement.

Figure 10 :
Figure 10: The behavior of the entropy (top left), Hawking temperature (top right), and gyromagnetic ratio (bottom) are shown for the four solution branches.The plot in the bottom right corner displays a zoomed-in region of the plot in the bottom left corner.The dashed gray line shows the profiles for a Kerr-Newman black hole for reference.The solutions were obtained with χ = 0.35 and g aγγ = 50.