Kerr black holes with synchronised axionic hair

We construct and analyse Kerr black holes (BHs) with synchronised axionic hair. These are the BH generalisations of the recently constructed rotating axion boson stars [1]. Such BHs are stationary, axially symmetric, asymptotically flat solutions of the complex Einstein-Klein-Gordon theory with a QCD axion-like potential. They are regular everywhere on and outside the event horizon. The potential is characterised by two parameters: the mass of the axion-like particle, ma and the decay constant fa. The limit fa → ∞ recovers the original example of Kerr BHs with synchronised scalar hair [2]. The effects of the non-linearities in the potential become important for fa . 1. We present an overview of the parameter space of the solutions together with a study of their basic geometric and phenomenological properties, for an illustrative value of the coupling that yields a non-negligible impact of the self-interactions.


Introduction
Testing the Kerr hypothesis is a central goal of current strong gravity research [3,4]. This is the hypothesis that astrophysical black holes (BHs), when near equilibrium, are well described by the Kerr metric [5]. This working assumption is intimately connected with the no-hair conjecture [6], which states that the dynamically formed equilibrium BHs have no other macroscopic degrees of freedom beyond those associated with Gauss laws (and hence gauge symmetries) -see e.g. [7][8][9][10] for reviews.
Over the last few years it has been realized that the Kerr hypothesis can be challenged even within General Relativity (GR) and even with physically simple and reasonable energy-matter contents, due to the discovery that free, massive complex scalar or vector fields can endow the Kerr BH with synchronised bosonic "hair" [2,[11][12][13], see, e.g. [14][15][16][17][18][19][20][21][22] for generalizations. Moreover, if these fields are sufficiently ultra-light, this "hair" could occur in the mass range of astrophysical BH candidates, spanning the interval from a few solar masses, ∼ M (stellar mass BHs), to ∼ 10 10 M (supermassive BHs).
The existence of such hairy BHs circumvents various no-scalar (and no-Proca) hair theorems -see e.g. [7,23,24]. Notwithstanding, the key test to the no-hair conjecture is a dynamical one. Are these hairy BHs dynamically robust? That is, can they form dynamically and be sufficiently stable? The current understanding is that: 1) there are dynamical formation channels, namely: a) the superrading instability [25] of the Kerr solution [26,27] and b) mergers of bosonic stars [28]; 2) these hairy BHs may, themselves, still be afflicted by superradiant instabilities [29,30] but these can be very long-lived, with a time scale that can exceed a Hubble time [31]. The evidence so far, therefore, indicates BHs with synchronised ultralight bosonic hair are an interesting challenge to the Kerr hypothesis even in GR. This has led to different phenomenological studies, including the study of their shadows [32][33][34] and their X-ray phenomenology, namely the iron Kα line [35,36] and QPOs [37].
An important outstanding issue is the fundamental physics nature of the putative ultralight bosonic field that could endow BHs with "hair". Possible embeddings in high energy physics have been suggested, such as the string axiverse [38], which proposes a landscape of ultralight axion-like particles emerging from string compactifications. If these particles have sufficiently small couplings with standard model particles, they are dark matter, only observable via their gravitational effects, one of which would be endowing spinning BHs

The model
The theory for which we shall obtain hairy BH solutions is the Einstein-Klein-Gordon model, that describes a massive complex scalar field, Ψ, minimally coupled to Einstein's gravity. The action of the theory is, using units such that G = c = = 1, where R is the Ricci scalar, Ψ is the complex scalar field ('*' denotes complex conjugation) and V is the scalar self-interaction potential. The equations of motion resulting from the variation of the action with respect to the metric, g µν , and scalar field, are, where is the energy-momentum tensor associated with the scalar field. A global U (1) transformation Ψ → e iχ Ψ, where χ is a constant, leaves the above action invariant; thus it is possible to write a scalar 4-current [2], j µ = −i (Ψ * ∂ µ Ψ − Ψ∂ µ Ψ * ) which is conserved: D µ j µ = 0. The existence of this symmetry and conserved current implies the existence of a conserved quantity -the Noether charge -that can be computed by integrating the timelike component of the 4-current, This quantity is interpreted as the number of scalar particles in a given solution, albeit this relation only becomes rigorous after field quantisation. For solitonic solutions (without an event horizon), moreover, Q is related with the total angular momentum as [44,45] This is a generic relation for rotating boson stars, already observed in other models with a self-interactions potential -see, e.g. [50]. As BH generalisations of the work done on rotating axion boson stars [1], we are interested in stationary, regular on and outside the event horizon, axi-symmetric and asymptotically flat solutions of the above equations of motion. The spacetime generated by these solutions possesses two commuting Killing vector fields, ξ and η, which, in a suitable coordinate system, can be written as ξ = ∂ t and η = ∂ ϕ corresponding to stationarity and axi-symmetry, respectively. With these Killing vector fields, we can define the following ansatz for the metric and the scalar field, written in terms of the spheroidal coordinates, {t, r, θ, ϕ}, where {F 1 , F 2 , F 0 , W ; φ} are ansatz functions that depend exclusively on the (r, θ) coordinates; r H is the radial coordinate of the event horizon; ω is the angular frequency of the scalar field; and m = ±1, ±2, . . . is the azimuthal harmonic index of the scalar field. In this work, we aim to study BHs with a surrounding axion-like scalar field; thus we specify the scalar selfinteraction potential as describing axion interactions. Following [1,47], we will use the QCD axion potential [51] to describe the axion interactions, added to a constant term in order to obtain asymptotically flat solutions. The potential in question has the form, 48 is the ratio between the up and down masses, m u and m d , respectively; µ a and f a define the axion-like particle (ALP) mass and quartic self-interaction coupling, respectively, and are two free parameters. This can be understood by performing an expansion of the potential around φ = 0, In this way, we confirm µ a is the ALP mass and f a is the quartic self-interaction coupling, In this work, we shall refer to µ a and f a as the ALP mass and decay constant, respectively. We note that the above expansion is only valid on the regime where φ f a . In fact, to leading order, only the mass term remains. Therefore, as already noted in [1] for the case of boson stars, for sufficiently large decay constant f a , we expect the BHs solutions to become very similar to the original Kerr BHs with synchronised scalar hair, obtained in [2].
In a similar fashion to the family of Kerr BHs with synchronised scalar hair, the solutions obtained in the work are possible due to the so-called synchronisation condition. This condition can be interpreted as a synchronisation between the angular velocity of the event horizon, Ω H , corresponding to the value of the metric function W at the horizon, Ω H = W (r H ), and the phase angular velocity of the scalar field, ω/m, cf. (8): Finally, let us mention that, as remarked in [1], the potential (9) allows for the existence of solutions even in the absence of the gravity term in the action (1). The simplest case corresponds to (non-gravitating) Q-ball-like solitons in a flat spacetime background. As expected, these solutions possess generalizations on a Kerr BH background. These bound states are in synchronous rotation with the BH horizon, i.e. they still obey the condition (12) and share most of the properties of the non-linear Q-clouds in [52]. A discussion of these aspects will be reported elsewhere.

Boundary conditions
To obtain numerical BH solutions appropriate boundary conditions must be imposed that enforce the sough physical behaviours. Such boundary conditions will now be summarised.
• Asymptotically boundary conditions: Asymptotic flatness implies that all ansatz functions must go to zero asymptotically, lim • Axial boundary conditions: Axial symmetry, together with regularity on the symmetry axis, implies that, Furthermore, we require the absence of conical singularities, thus, F 1 = F 2 on the axis. Since we focus on even parity solutions, which typically correspond to the fundamental solutions, they are symmetric w.r.t. the equatorial plane. Hence, one only needs to solve the equations of motion in the range 0 θ π/2 and impose the following boundary conditions, • Event horizon boundary conditions: To simplify the study of these boundary conditions, let us introduce a radial coordinate transformation, x = r 2 − r 2 H . With this coordinate transformation, we can perform a series expansion of the ansatz functions at the horizon, x = 0, and find that, where i , φ (0) are functions of θ. With these series expansions, we can naturally impose the following boundary conditions,

Extracting physical quantities
The main physical quantities of interest for our analysis are encoded in the metric functions evaluated either at the horizon or at spacial infinity. For the former, we obtain the horizon angular velocity by using the horizon boundary condition mentioned in the previous sections, Ω H = W | r H , together with Hawking temperature, T H , as well as the horizon area, A H , through the following expressions, The entropy of the computed BH follows from the Bekenstein-Hawking formula, S = A H /4. For the latter, we can compute the ADM mass, M , and angular momentum, J, through the asymptotic behaviour of g tt and g tϕ , The above quantities are mutually related through a Smarr-type formula [53], where we have introduced two new quantities, which can not be computed directly from the horizon or spacial infinity data, corresponding to the scalar field mass and angular momentum, respectively 2 . In their definition, Σ is a spacelike surface, bounded by the 2-sphere at infinity, S 2 ∞ , and the spatial section of the horizon, H. Moreover, the angular momentum of the scalar field is related to the Noether charge which arises from the global U (1) symmetry of the scalar field, Q, as J Ψ = mQ. For a solution composed entirely of axionic scalar hair (no horizon), the ADM angular momentum must be equal to the angular momentum of the scalar field; likewise, for a solution without axionic scalar hair, the ADM angular momentum of the solution equals the horizon angular momentum. Thus, in order to evaluate how hairy a given BH solution is, we define the following dimensionless parameter: One easily sees that, when q = 0, we have a bald BH, corresponding to a Kerr BH. On the other end, when q = 1, we have a solution entirely compose of axionic hair, which corresponds to a rotation axion boson star [1]. For any other value 0 < q < 1, we have a rotating BH surrounded by a non-trivial, backreacting, massive rotating axion scalar field.

Numerical approach
To perform the numerical integration of the equations of motion resulting from Eqs. (2) and (3) with the ansatz Eqs. (7) and (8), it is useful to rescale key quantities by µ a , together with f a → f a √ 4π. This leads to the disappearance of the ALP mass constant from the equations of motion numerically solved, but all global quantities will be express in terms of µ a .
In our approach, by expanding the equations of motion, we get a set of five coupled, non-linear, elliptic partial differential equations for the ansatz functions, F a = {F 0 , F 1 , F 2 , W ; φ}. They are compose of the Klein-Gordon equation, Eq. (3), together with the following combination of the Einstein equations, Eq. (2), The remaining equations, E r θ = 0 and E r r − E θ θ = 0 are not solved directly, instead, they are used as constraint equations to evaluate the accuracy of the numerical solution. Typically they are satisfied at the level of the overall numerical accuracy.
Our numerical treatment can be summarised as follows. We restrict the domain of integration to the region outside the horizon. Using the aforementioned radial coordinate transformation x = r 2 − r 2 H , we transform the radial region of integration from [r H , ∞) to [0, ∞). Then, we introduce a new radial coordinate that maps the semi-infinite region [0, ∞) to the finite region [0, 1]. Such map can be defined in several ways, but in this work we choose the new coordinate,x = x/(x + 1). After this, the equations F a are discretised on a grid inx and θ. Most of the results presented here were obtained on an equidistant grid with 251 × 30 points. The grid covers the integration region 0 x ≤ 1 and 0 θ ≤ π/2.
The equations of motion have been solved subject to the boundary conditions introduced above by using a professional package, entitled FIDISOL/CADSOL [54], which employs a Newton-Raphson method with an arbitrary grid and consistency order. This code uses the finite difference method, providing also an error estimate for each unknown function. For the solutions in this work, the maximal numerical error for the functions is estimated to be on the order of 10 −3 . The Smarr relation, Eq. (22), provides a further test of the numerical accuracy, leading to error estimates of the same order.
In our scheme, there are four input parameters: i) the decay constant f a in the potential, Eq. (9); ii) the angular frequency of the scalar field ω; iii) the azimuthal harmonic index m; and iv) the radial coordinate of the event horizon r H . The number of nodes n of the scalar field, as well as all other quantities of interest mentioned before, are computed from the numerical solution. For simplicity, we have restricted our study to the fundamental configurations, i.e. with a nodeless scalar field, n = 0 and with m = 1. Also, from the results presented in [1], we shall illustrate the effect 3 of the axion potential on the hairy BHs by performing a thorough study of the solutions with the specific decay constant f a = 0.05.

The domain of existence
At the end of the previous section, we fixed two of the four input parameters of the problem (f a , m). Thus, the full domain of existence is obtained by varying the remaining two input parameters: the angular frequency of the scalar field, ω; and the radial coordinate of the event horizon, r H . Since it is impossible to obtain all possible BH solutions, we obtained a very large number of them (∼ 30000) and we have extrapolated this large discrete set of solutions into the continuum, which defines the region where one can find the BH solutions with axionic hair.
Such a region can be expressed and plotted in various ways. In the left panel of Fig. 1 (main panel), we show it in an ADM mass M µ a vs. angular frequency ω/µ a plane. We can observe that most (but not all) of the numerical solutions region is bounded by two specific lines 4 , • The axion boson star line -corresponding to the solitonic limit, in which both the event horizon radius and area vanish, r H = 0 and A H = 0, and the solutions have no BH horizon; therefore q = 1. Such line is represented in both panels of Fig. 1 as a red solid line.
• The existence line -corresponding to specific subset of vacuum Kerr BHs which can support stationary scalar clouds (with an infinitesimally small φ), first discussed by Hod [55,56], thus having q = 0. These solutions are obtained by linearising the theory, and since, on that regime, the axion self-interacting potential reduces to the mass potential -cf. Eq. (10) -the existence line will be the same as the one obtained for the family of Kerr BHs with scalar hair. Such line is represented in both panels of Fig. 1 as a blue dotted line.  . The consequence of this observation is that for certain frequencies, the ADM mass is not maximised by a boson star, but rather by a hairy BH.
By changing the representation, however, and plotting the domain of existence in the ADM angular momentum Jµ 2 a vs. angular frequency ω/µ a plane -right panel of Fig. 1 -, we see that, for the angular momentum, the boson star line bounds all axionic BHs with decay constant f a = 0.05. Thus for all frequencies, the angular momentum is maximised by a boson star.
Another distinctive feature of the domain of solutions of the axionic BHs is the existence of a local maximum for the mass and angular momentum at ω ∼ 0.9, which is not the global maximum. For solutions with smaller angular frequency, it is possible to have BHs with more mass and angular momentum than the ones near the local maximum; in fact, these quantities are maximised for the solutions with the smallest possible value of angular frequency. Such is not the case in the absence of the axionic potential (inset of both panels in Fig. 1).
In Fig. 2 both the ADM angular momentum, J and the dimensionless spin, j, defined as j ≡ J/M 2 , are exhibited vs. the ADM mass, M , on the left and right panels, respectively. In the former, we can see that the axionic BHs can have a higher mass and angular momentum than their f a → ∞ counterparts (inset). We can also see a zigzag of the boson star line. This behaviour is explained by the sudden drop on the ADM mass and angular momentum around ω/µ a ≈ 0.84 -cf. Fig. 1. In the right panel, we see a considerable violation of the Kerr bound, j 1 in part of solution space. This already occurred for the f a → ∞ limit (inset). Again, it is possible to visualise the zigzag behaviour of the boson star line.
Let us now study the horizon geometry of the axionic BHs. Their event horizon has a spherical topology but a spheroidal geometry, similarly to Kerr BHs. This can be seen by studying the spatial cross-section of the horizon, through the induced metric, Due to the rotation of the solutions, the horizon is squashed at the poles. To show this, we compute the horizon  Fig. 1 the insets show the free scalar field limit f a → ∞, where the extremal line is also shown.
circumference along the equator, L e , and along the poles, L p , We define the sphericity as the ratio of both circumferences above [57], For values in which the sphericity is greater (lower) than 1, the horizon will be squashed (elongated) at the poles, leading to an oblate (prolate) spheroid. From the left panel in Fig. 3 we see that all solutions have a sphericity larger than the unity; thus all solutions have an oblate horizon, as expected.  Another physical quantity of interest associated with the horizon is its linear velocity v H [21,57,58]. Such quantity measures how fast the null geodesics generators of the horizon rotate relatively to a static observer at spatial infinity. Its definition is quite simple, only taking into account the perimetral radius of the circumference located at the equator, R e ≡ L e /2π, and the horizon angular velocity, The horizon linear velocity is presented in the right panel of Fig. 3. The central feature in this plot is the fact that all solutions have horizon linear velocity smaller than the unity, which, in the units we are using, corresponds to the speed of light. Therefore, null geodesics generators of the horizon never rotate relatively to the asymptotic observer at superluminal speeds, even though some solutions strongly violate the Kerr bound j 1.
A final insight about the horizon geometry of the axionic BHs is obtained from investigating whether an isometric embedding of the spatial sections of the horizon is possible in Euclidean 3-space E 3 . For a Kerr BH, such embedding is possible iff its dimensionless spin obeys j j (S) [59], where j (S) ≡ √ 3/2 was dubbed the Smarr point [57]. For j > j (S) the Gaussian curvature of the horizon becomes negative in the vicinity of the poles [59], which, prevents the embedding (due to occurring at a fixed point of the axi-symmetry). As expected, this feature also occurs for the axionic BHs. Due to the existence of scalar hair around the BH, we have an extra degree of freedom, which converts the Smarr point into a Smarr line. Such line is represented in both panels of Fig. 3 as a solid black line. One observes that, for both the sphericity and the horizon linear velocity, the Smarr line is not constant. This contrasts with the behaviour for f a → ∞; in that case, the sphericity of the Smarr line was constant and equal to the value of the Smarr point in Kerr [57]. Thus, the axion potential destroys the constancy of the sphericity along the Smarr line.

Ergo-regions
An ergo-region is a part of a spacetime, outside the event horizon, wherein the norm of the asymptotically timelike Killing vector ξ = ∂ t becomes positive and thus the vector becomes spacelike. Ergoregions are associated with the possibility of energy extraction from a spinning BH, via the Penrose process [60,61], or superradiant scattering [25]. In the context of BHs with synchronised hair, the superradiant instability of vacuum Kerr BHs in the presence of ultralight scalar fields is one of the possible channels of formation of these hairy BHs [26,27]. Thus, it is of relevance to analyse ergo-regions for the axionic BHs.
Kerr BHs possess an ergo-region whose boundary has a spherical topology and touches the BH horizon at the poles -such surface is called an ergo-sphere. For the axionic BHs in the f a → ∞ limit, the ergo-regions can be more complicated [29]. Some solutions have a Kerr-like ergo-region; but others have a more elaborate ergo-region topology, with two disjoint parts, one Kerr-like and another of toroidal topology. The latter were dubbed ergo-Saturns. The toroidal ergo-region is inherited from the mini boson star environment around the horizon, since these stars develop such ergo-regions, when sufficiently compact. Such ergo-torii also occur for spinning axion boson stars [1]. Thus, we expect some axionic BHs to develop an ergo-Saturn. This is confirmed in Fig. 4. We have found that the ergo-region of axionic BHs follows a qualitatively similar distribution to that of their f a → ∞ limit: there are solutions which possess an ergo-sphere and others that develop an ergo-Saturn. The latter occur on the far left of the domain of existence in Fig. 4, where the most massive BHs exist, having the lowest possible values of the angular frequency of the scalar field.

Light rings and timelike innermost stable circular orbits
A phenomenological aspect of importance is the structure of circular orbits (COs) of both massless and massive particles around a BH. In particular, the (timelike) innermost stable circular orbit (ISCO) and the (null) light rings (LRs) are of special relevance. The former is associated with a cut-off frequency of the emitted synchrotron radiation generated from accelerated charges in accretion disks; the latter is related to the real part of the frequency of BH quasi-normal modes [62], as well as to the BH shadow [63]. The structure of COs can be obtained as follows. Given the geometry in Eq. (7), we can compute the effective Lagrangian for equatorial, θ = π/2, geodesic motion as, where all ansatz functions depend only on the radial coordinate r, the dot,˙, denotes the derivative w.r.t the proper time, and = {−1, 0} for massive (timelike) particles and massless (lightlike) particles, respectively. Due to the existence of two Killing vector fields, we can writeṫ andφ in terms of the energy E and angular momentum L of the particle, Inverting the above system of equations and replacing the result into the effective Lagrangian, Eq. (34), we can obtain an equation forṙ, which defines a potential V (r), In order to obtain COs, both the potential and its derivative must be zero, i.e. V (r) = V (r) = 0. Depending on whether we are considering massless or massive particles, these two equations will yield different results.
For massless particles, the first equation, V (r) = 0, will give two algebraic equation for the impact parameter of the particle, b + ≡ L + /E + and b − ≡ L − /E − , corresponding to co-and counter-rotating orbits, respectively. The second equation, V (r) = 0, together with the impact parameters, will give the radial coordinate of the coand counter-rotating LRs. Whenever it is possible to obtain a real solution for the radial coordinate, the BH possesses LRs.
In Fig. 5, we show the distribution of hairy BHs with different number of LRs in the angular frequency, ω/µ a vs. angular momentum, Jµ 2 plane. The left panel shows the counter-rotating case, whereas the right panel shows the co-rotating case. In both cases, it is always possible to have at least one LR, as for the Kerr BH. In the counter-rotating case, however, if the surrounding scalar field is compact enough, an extra pair of LRs emerge, leading to a hairy BH with 3 LRs. This is the case of a large set of hairy BHs with low values of ω/µ a . In fact, we see the same behaviour for the free scalar field case (inset plot); even the solid black line separating the two regions has a qualitatively similar shape. The radial coordinate of the several LRs is shown in Figs. 6 and 7 as a blue dashed line. For massive particles, V (r) = V (r) = 0 will yield two algebraic equations for the energy and angular momentum of the particle, {E + , L + } and {E − , L − } corresponding to co-and counter-rotating orbits, respectively. Their stability can be verified by analysing the sign of the second derivative of the potential V (r). Given a BH solution it will have three distinct regions concerning (timelike) COs: • No circular orbits (No COs) -Whenever we obtain an imaginary solution for the energy and angular momentum of the massive particle. This is a region in which no COs exist. Such region is shaded in light red, in Figs. 6 and 7.
• Unstable circular orbits (UCOs) -In this region, it is possible to obtain real solutions for the energy and angular momentum, but the sign of the second derivative of the potential V (r) is positive, implying that COs are unstable. Such region is shaded in light yellow, in Figs. 6 and 7.
• Stable circular orbits (SCOs) -In this region, it is also possible to obtain real solutions for the energy and angular momentum, but now the second derivative of the potential V (r) has a negative sign, thus implying that COs are stable. Such region is shaded in light green, in Figs. 6 and 7.
The ISCO, as the name entails, is the innermost stable circular orbit, i.e., the stable CO with the smallest radial coordinate. At this orbit, the second derivative of the potential vanishes, as it corresponds to the transition between the SCOs region and the UCOs region. The ISCO is represented as a solid red line in Figs. 6 and 7.
In Fig. 6 we present the structure of counter-rotating COs for four sets of hairy BHs with constant q = {0.5, 0.8, 0.9, 0.999}. For each set of BHs, we have constructed a plot showing the maximal value of the scalar field, φ max , as a function of the radial coordinate x = r 2 − r 2 H . This ensures that every horizontal line (lines with constant φ max ) corresponds to a unique BH solution making it easier to analyse the structure of COs.
The first observation from Fig. 6 is that, by increasing q, it is possible to obtain solutions with larger φ max . This is reasonable since by increasing q the BHs become hairier. Let us now analyse each individual plot. For the case where q = 0.5 (right bottom panel), the structure of the counter-rotating COs for all BHs is identical to the one for Kerr. There is only one LR and one ISCO, which separates the No COs region from the UCOs region, and separates the UCOs region from the SCOs region, respectively. At large enough x we only have SCOs, but as we approach the horizon, we reach the ISCO and consequently enter the UCOs region. If we continue towards the horizon, we eventually cross the LR and enter the No COs region. The radial coordinate of the LR and ISCO increases monotonically with the increase of the maximal value of the scalar field.
For the case where q = 0.8 (left bottom panel) most of the structure is similar, but a new feature emerges. Starting at the regime where there are BHs with a dilute scalar field (φ max < 0.11), their structure is the same as for Kerr, but for BHs with φ max > 0.11 and φ max < 0.15 an extra region of UCOs appears besides the one already present between the LR and the ISCO, and disconnected from the latter. The new region appears as a single point around xµ a ≈ 4.8 and increases in size as φ max increases until its inner boundary merges with the ISCO. For φ max > 0.15, the structure is again the same as that of Kerr, but now the radial coordinate of the ISCO is significantly larger than for BHs with a more dilute scalar field. In fact, there is a discontinuity when we study the evolution of the radial coordinate of the ISCO as it changes with the maximum value of the scalar field, as one can see in the left bottom panel of Fig. 6.
In the q = 0.9 case (right top panel), we have a structure similar as for the q = 0.8 when we consider BHs with φ max < 0.21, but, for the remaining BHs, the scalar field environment is compact enough to develop a pair of extra LRs. Such pairs give rise to a new region of No COs, disconnected from the already existing one between the event horizon and innermost LR. It starts as a single point for a BH with φ max ≈ 0.21; then it grows in size until its inner boundary (one of the LRs) merges with the innermost LR, connecting both regions of No COs. At the same time, the ISCO also merges with the LRs, meaning that such BH has a degenerate point in which two LRs and the ISCO converge. For BHs with a larger value of φ max , a Kerr-like structure again emerges, but now both the LR and the ISCO appear at a larger radial coordinate than for small φ max .
Lastly, in the q = 0.999 case (left top panel), the most complex structure is observed. This case inherits the several features discussed in the previous cases, such as the existence of two different and disconnected UCOs and No COs regions, but also presents a new one. For BHs with φ max > 0.18, a third region of UCOs between the two already existing ones. This third region starts as a single point and increases slowly in size as φ max increases until its inner boundary connects with the ISCO, and the innermost region of UCOs merges with this new region. Although it is not visible in the left top panel of Fig. 6, one can draw the conclusion that there is a BH with larger φ max for which the innermost two LRs together with the ISCO converge to the same single point, akin to what happened in the q = 0.9 case. Now we turn to the co-rotating COs. In Fig. 7 we followed the same idea that we used to study the counterrotating COs and plotted four sets of BHs solutions with constant q = {0.5, 0.8, 0.9, 0.999} in the φ max vs. x plane. Fig. 7 manifests that the structure of co-rotating COs is much simpler than the counter-rotating one; only the q = 0.999 case presents significant differences from the other ones; the structure of the latter is the same as for a Kerr BH. Moreover, the radial coordinate of the co-rotating LR and ISCO is always smaller than the one of the counter-rotating LR and ISCO, as it is for the Kerr BH [64]. Let us comment on the qualitatively different q = 0.999 case (first panel of Fig. 7). For BHs with a dilute scalar field, a Kerr-like structure is observed; but above φ max ≈ 0.14 a second region of UCOs emerges. This region exists until we reach a BH with φ max ≈ 0.30, where both regions of UCOs merge together. For BHs with slighter larger φ max a second region of UCOs again emerges, but now this region occurs at large radii, quickly decreasing to smaller radii as φ max increases. Although not shown in this panel, we can predict that this new region of UCOs will converge, as the previous one, with the already existing and closer to the event horizon UCOs region.

Conclusions and Remarks
In this work, we have constructed and analysed BHs with synchronised axionic hair, which are BH generalisation of the rotating axion boson stars recently found in [1] -see also [47]. These are stationary, axi-symmetric, asymptotically flat and regular on and outside the event horizon solutions of the Einstein-Klein-Gordon equations of motion with a QCD axion-like potential, Eq. (9). This family of axionic BHs is described by three parameters: the radial coordinate of the event horizon, r H , the angular frequency of the scalar field, ω and the decay constant of the QCD potential, f a . In this work, we have thoroughly scanned the space of solutions with decay constant f a = 0.05, since, from the results found in [1], this yields a case with considerable impact of the axion potential, and hence considerable differences from the free scalar field case, obtained as the f a → ∞ limit. The latter yields the original example of Kerr BHs with synchronised scalar hair [2]. Even larger deviations from this original example may occur for even smaller f a , but then the numerics to obtain such solutions becomes more challenging.
When comparing the axionic BHs with their free scalar field counterparts [2], there are both differences and similarities. Some key differences are: (i) axionic BHs can have more mass and angular momentum and lower values of the angular frequency of the scalar field; (ii) the existence of a local maximum for the mass and angular momentum, that does not exists for the f a → ∞ case; (iii) the presence of a small region of frequencies where the mass of axionic BHs is no longer bounded by the axion boson stars. In such region, we have a degeneracy of solutions with the same angular frequency and ADM mass, but such degeneracy is easily lifted by specifying any other physical quantity; (iv) the variation of the sphericity along the Smarr line, which, in the free scalar field case was constant and equal to the sphericity of the Smarr point for Kerr BHs, but varies for axionic BHs.
Concerning the similarities, we may emphasise: (i) the clear violation of the Kerr bound, j 1; (ii) the sphericity and horizon linear velocity are bounded by the same existence line. Thus, both families can only have the same values of sphericity and horizon linear velocity as the Kerr ones, which implies that all BHs have a horizon which is an oblate spheroid and the rotation of its null generators (relatively to a static observer at spatial infinity) never exceeds the speed of light; (iii) the topology of the ergo-regions is either an ergo-sphere or an ergo-Saturn. In both families, the ergo-Saturn only appears for the solutions with the lower values of angular frequency.
In this work, we also presented a study of the structure of COs. Counter-rotating COs show a more complex structure than the co-rotating COs counterpart. Concerning the former, we saw that solutions with low or moderate amounts of hair (e.g. q = 0.5) will have a structure similar to the Kerr one. For very hairy solutions (q 0.8), a new region of UCOs can emerge, as well as a new region of No COs when the axionic hair outside the event horizon is compact enough to yield an extra pair of LRs. In the extreme case, where most of the BH solution is composed of axionic hair (q = 0.999), it possible to have a third region of UCOs, leading to an intercalation of SCOs and UCOs regions.
Let us briefly comment on energy conditions. In [1], it was shown that the weak energy condition (WEC) and the dominant energy condition (DEC) always hold for the case of axion boson stars, whereas the strong energy condition (SEC) could be violated, and in fact, it is violated for a zero angular momentum observer (ZAMO). For the BH generalisation we can prove, following [1], that the WEC and DEC are never violated, but the SEC can be violated, for instance, for a ZAMO.
Due to the groundbreaking results presented by the LIGO-Virgo collaboration detecting several gravitational waves events, starting with [65], as well as the results by the EHT collaboration on the first image of a shadow of the M87 supermassive BH [66], two possible and interesting directions to follow up on this work are: (1) to study the possible gravitational waves generated by the collision of two axionic hairy BHs, both in the head-on scenario, as well as, in the more realistic scenario of an in-spiral binary system. Recently [67] such collisions were made for Proca stars [68] obtaining a suggestive agreement with the particular gravitational wave event GW190521 [69]; (2) and to study the shadow of axionic hairy BHs, to analyse the impact of the QCD axion-like potential in the BH shadows. This would be a generalisation of the analysis in [32].
Union's Horizon 2020 research and innovation (RISE) programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. The authors would like to acknowledge networking support by the COST Action CA16104.