Black holes with synchronised Proca hair: linear clouds and fundamental non-linear solutions

Recent studies have made key progress on the black hole/solitonic solutions of the Einstein-Proca system. Firstly, fully non-linear dynamical evolutions of the Kerr black hole superradiant instability, triggered by a Proca field, have shown the formation of a new equilibrium state, a spinning black hole with synchronised Proca hair. Secondly, non-linear evolutions of spinning Proca stars have established that they are dynamically stable, unlike their scalar cousins. Thirdy, separability of the Proca equation on the Kerr background has been achieved. Motivated by these results, in this paper we reconsider Kerr black holes with synchronised Proca hair. The separability of the Proca equation on the Kerr background allows us to examine the stationary Proca clouds in greater detail, in particular their dependence on the different quantum numbers. These stationary clouds occur at a set of existence lines in the Kerr parameter space, from which the black holes with synchronised Proca hair bifurcate. We construct the domain of existence of these black holes, comparing the fundamental states missed in the original study with the first excited states and with the cousin scalar model, giving illustrative examples of Kerr-like and non-Kerr-like BHs. In the vanishing event horizon limit, these hairy black holes connect to the fundamental states of spinning Proca stars, which include the dynamically stable solutions.


Introduction
Successfully tested up to the T eV scale, the Standard Model of particle physics turns out to describe but a tiny fraction of the matter-energy density of the Universe. It does not explain the phenomenological evidence for the existence of dark matter and dark energy, which are believed to make up about 95% of the Cosmos [1]. Many models have been put forward to explain the dark side of the Universe. Some, in particular, relate dark matter to hypothetical new, ultralight bosonic particles which are sufficiently weakly coupled to ordinary matter to have remained elusive to past and present experimental searches [2,3].
Bosonic particles have an interesting interaction with Kerr black holes (BHs). They can extract the BH's rotational energy through a radiation enhancement mechanism known as superradiance [4]. For Kerr BHs, superradiance occurs when the phase angular velocity of the boson, ω, fulfills the condition where m j is the boson's azimuthal total angular momentum and Ω H , r H = M + √ M 2 − a 2 are, respectively, the BH's horizon angular velocity and event horizon (Boyer-Lindquist) radial coordinate, in terms of the BHs's ADM mass M and total angular momentum J = M a. The enhancement is most efficient when the reduced Compton wavelength of the boson, λ C = /(µc), is comparable to the BH's gravitational radius, R G = GM/c 2 , i.e.
where µ is the boson's mass; α is the so-called gravitational fine-structure constant. For the known astrophysical BH masses, ranging between 1 − 10 10 M , this implies that the bosonic particles are ultralight, with a mass range of roughly 10 −10 − 10 −20 eV . The non-vanishing bosonic field mass plays the role of a mirror, trapping the bosons in the vicinity of the BH and creating a recurrent energy/angular momentum enhancement of the bosonic state. At the linear level, i.e. disregarding the bosons' backreaction on the background spacetime, the energy feeding of the particles fuels an exponential growth known as superradiant instability, or 'BH bomb' [5]. At the non-linear level, the exponential superradiant growth stalls when the inequality (1.1) saturates, i.e. One may say the (phase angular velocity of the) cloud and the (horizon angular velocity of the) BH synchronise. A simple entropic estimate shows that up to about 29% of the BH's energy could be mined, in an astrophysical timescale, by this process. The result is a classical condensate (often dubbed as cloud but also as BH 'hair ') which is stationary with respect to the slowed-down BH [6][7][8]: a Kerr BH with synchronised bosonic hair. These are stationary BH solutions of Einstein's gravity minimally coupled to complex bosons, first discussed in [9] for scalar and in [10] for vector bosons. They challenge the no-hair hypothesis [11] (see also [12,13]) even in General Relativity. According to this hypothesis, BHs that could form dynamically in the presence of astrophysically (potentially) relevant generic matter-energy are fully characterised by global charges associated with Gauss laws, such as M and J, and have no other degrees of freedom, broadly referred to as 'hair'. The domain of existence of BHs with synchronised bosonic hair has two important boundaries. Firstly, for vanishing horizon size it yields the set of spinning bosonic stars, which have long been known in the scalar case [14,15], but only recently constructed in the vector case, a.k.a. Proca stars [16]. Very recently, it has been shown that the spinning scalar stars suffer from a non-axisymmetric instability, whereas the spinning Proca stars are dynamically robust [17]. This suggests that the Proca case may be dynamically more interesting. Secondly, for vanishing bosonic field, the hairy BHs bifurcate from the Kerr family at the Kerr solutions that admit linear bound states of the corresponding massive bosonic field. These states exist at the threshold of superradience, i.e. when Eq. (1.3) holds, and are commonly known as stationary clouds.
Stationary clouds around Kerr BHs were first found in the scalar case and around an extremal (a = M ) BH [18]. Remarkably, in this particular case the radial function can be solved analytically in terms of confluent hypergeometric functions. This analysis has then be extended, typically using numerical methods, to other regimes and other BHs -see e.g. [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. Stationary clouds in the Kerr case are analogous to the atomic orbitals in the hydrogen atom [37] -see also [38]. They are finite on and outside the BH's event horizon, decay exponentially at spatial infinity and can be labeled by four quantum numbers: n, the number of nodes of the radial function; , the orbital angular momentum; j, the total angular momentum; and m j , the projection of the total angular momentum along the BH's axis of rotation. Similar configurations have also been obtained in analogue models of gravity such as the draining bathtub vortex [39,40].
Most studies of stationary clouds around rotating BHs have focused on the scalar case, whose equations of motion are separable on the Kerr spacetime. As for the massive vector bosons, it remained unclear for decades whether the Proca equation was separable or not on Kerr and the only study of clouds tackled the problem by solving the corresponding partial differential equations [10] -see also [41]. Recently, however, the separability of the Proca equation for a large family of spacetimes that includes the Kerr BH was established using a proper ansatz [hereafter the Frolov-Krtouš-Kubizňák-Santos (FKKS) ansatz] [42]. This development has allowed more detailed studies of the Proca superradiant instability -see e.g. [43,44]. The first goal of this paper is to make use of this development to determine and characterize the stationary vector clouds around Kerr BHs in terms of {n, , j, m j }.
The stationary (scalar or vector) clouds define an existence line on the Kerr parameter space from which the BHs with synchronised hair bifurcate [9,10]. There is a discrete set of families of BHs with synchronised hair, labelled by the parameters (n, m j ); the parameters ( , j) do not have significance when going from the linear to the non-linear theory. In the scalar case, the fundamental family of hairy BHs has nodeless scalar field profiles, corresponding to n = 0 and m j = 1 [45]; nodeful solutions, i.e. with n = 0, are excited states with higher energy [46]. The same holds for the solitonic limit. In the particular case of spherical, static scalar boson stars (m j = 0), it has been shown dynamically that the n > 0 excited states decay into the fundamental n = 0 ground state [47].
In the original study of Proca stars [16] it was proved that, for static spherical Proca stars, one of the profile functions of the Proca potential must have at least one node; there are no nodeless solutions. In consistency with this observation, the spinning Proca stars reported in the same paper had one node for the corresponding function. Subsequently, the original study of BHs with synchronised Proca hair constructed BH solutions that also have one node of the same function [10]. It was observed in [8] when interpreting the results in [7], however, that in the spinning case (but not in the static case) there are nodeless Proca stars, and also hairy BHs with Proca hair, and these are the true fundamental states. Nonetheless, the latter have not been studied in detail in the literature. The second goal of this paper is, therefore, to report a detailed study of the fundamental solutions of these hairy BHs. In particular, their solitonic limit corresponds precisely to the solutions that have been recently shown to be dynamically robust [17] -see also [48]. This paper is organised as follows. In Section 2 we consider the linear analysis of the stationary clouds on a fixed Kerr geometry. In Subsection 2.1 the relativistic quantummechanical description of vector bosons is briefly addressed. The notation introduced therein will be useful to label the stationary clouds. Subsection 2.2 reviews the Proca equation on a curved spacetime, introduces the FKKS ansatz for the Proca field and presents the radial and angular equations it yields for the Kerr case, in Boyer-Lindquist coordinates. Subsection 2.3 then sets the stage for the numerical integration of those equations and covers the results. Section 3 deals with the non-linear analysis. After briefly describing the setup in Subsection 3.1, the domain of existence of the fundamental BHs with Proca hair is discussed in Subsection 3.2 and compared with that of the first excited states of hairy BHs and the cousin scalar model. In Subsection 3.3 we analyse illustrative solutions of both hairy BHs and spinning Proca stars. Finally, a concise overview of the work is sketched in Section 4, together with some closing remarks on future prospects. Appendix A provides some illustrations of the vector spherical harmonics.
2 Linear analysis: stationary clouds on a fixed Kerr geometry

Vector bosons
In relativistic quantum mechanics, particles are described by the orbital angular momentum L and the intrinsic angular momentumŜ. The components of the individual operators satisfy the angular momentum commutation relations, i.e.
Vector bosons are characterized by s = 1, which means that the quantum number j can take a single value when = 0 (j = 1) and three different values when > 0 (j = − 1, , + 1). In this case, the eigenstates of the operatorĴ 2 in the spherical coordinate representation are the ('pure-orbital') vector spherical harmonics Y j,m j , which can be expressed in terms of the scalar spherical harmonics Y j,m j as The vector spherical harmonics have parityΠ = (−1) +1 . Thus, upon a parity transformation, Y j,m j acquires a factor of (−1) j , when j = ± 1, and of (−1) j+1 , when j = .
Y j±1 j,m j (Y j j,m j ) are said to have electric-type (magnetic-type) parity: they correspond to the magnetic (electric) field of electric multipole radiation and the electric (magnetic) field of magnetic multipole radiation [50,51]. The explicit form and a graphic representation of the first few 'pure-orbital' vector harmonics are provided in Appendix A.
In curved spacetimes, { , j, m j } are not in general legitimate quantum numbers, since curvature can break the conservation of angular momentum 1 . However, in Schwarzschild spacetime, the total angular momentumĴ is still conserved. This means that vector bosons only have definite total angular momentum. Such definiteness is broken in Kerr spacetime, in which the total angular momentum is no longer conserved. Nevertheless, choosingĴ z to be aligned with the symmetry axis of Kerr spacetime at spatial infinity, m j remains a conserved quantity.
It is convenient to use the quantum numbers { , j, m j } to identify vector bosons, always bearing in mind that they are only physically meaningful in Minkowski spacetime. In particular, in the following, vector states will be labelled with | , j, m j ≡ | , s = 1, j, m j .

Proca equation
The Lagrangian density of a massive complex vector boson A α reads where F αβ = 2A [β;α] . F αβ is the electromagnetic-field tensor, which is antisymmetric and gauge invariant, α, β = 0, 1, 2, 3 and µ is the boson's mass. The variation of the action integral S = V d 4 x √ −g L M with respect to the fieldĀ α leads to the Proca field equation [52]: Writing the equation in terms of the electromagnetic four-potential A α , its four-divergence reads Using the identities The dynamics of the divergenceless electromagnetic four-potential A α is thus encoded in a set of four Klein-Gordon equations, one per component. The non-trivial separability of the Klein-Gordon equation in Kerr spacetime was first unveiled via variable separation by Carter [53], shortly after noting the complete integrability of the Hamilton-Jacobi equation for Kerr geodesics [54]. Carter's seminal work broke down the original second-order partial differential equation (PDE) into two coupled second-order ordinary differential equations (ODEs), and paved the way for a thorough study of Kerr linear perturbations. Although the four equations of motion (2.11) are individually separable for a specific ansatz, the separability does not extend to the Lorenz condition (2.10). In fact, the separability of these five second-order PDEs is not trivial and was only achieved recently via the FKKS ansatz [42] (see also Ref. [55]), following [56]. This separability has been established for the Kerr-NUT-(A)dS family of spacetimes. The FKKS ansatz, which embodies the explicit and hidden symmetries of the metric, is where B αβ is the polarisation tensor 2 and Z is an auxiliary complex scalar function for which a multiplicative separation of variables will hold -cf. Eq.(2.17) below. The polarisation tensor B αβ is defined in terms of the principal tensor 3 h αβ as where λ plays the role of a separation constant.
To solve the Proca equation in the Kerr background with the ansatz (2.12) one proceeds as follows. In Boyer-Lindquist coordinates (t, r, θ, ϕ), the Kerr metric reads where Σ ≡ r 2 + a 2 cos 2 θ and ∆ ≡ r 2 − 2M r + a 2 . The Kerr spacetime is stationary and axisymmetric; it has an event horizon at r = r H , the largest root of ∆. In Boyer-Lindquist coordinates, the Killing vectors associated with these continuous symmetries are ξ t = ∂ t and ξ ϕ = ∂ ϕ , respectively. Its principal tensor reads Using the metric and principal tensor, one constructs the polarization tensor. Its symmetric and antisymmetric parts are Note that the r-dependent terms are decoupled from the θ-dependent terms, apart from the common factor 1/Σ. Additionally, we take the complex scalar function in Eq. (2.12) with the form where R and S are dubbed radial and angular functions, respectively, and ω and m j are the eigenvalues related to the aforementioned isometries. With this construction, the ansatz (2.12) reduces the Proca equation to the two separated equations As mentioned in Subsection 2.1, vector states can have electric-type or magnetic-type parity. The ansatz (2.12) encodes all the electric-type states and the magnetic-type states with j = |m j | [37,43]. It remains unclear, however, whether all the magnetic-type states are captured by the FKKS ansatz or not. If so, further study of the Proca equation in Kerr spacetime may shed some light on how to recover all of them. If not, a new question arises: whether it is possible to find a new ansatz which contains all the states.
In the following, all physical quantities will be expressed in terms of the boson's reduced Compton wavelength. It is therefore convenient to set λ C = µ −1 = 1.

Stationary vector clouds around Kerr black holes
In general, Eqs. (2.18) and (2.19) form a non-standard coupled eigenvalue problem with an eigenvalue pair {ω, λ}. To construct the stationary vector clouds around Kerr BHs, we need to find the bound states whose phase angular velocity fulfills the synchronization condition (1.3). Since ω is fixed a priori, the eigenvalue pair can be chosen to be either {M, λ} or {a, λ}. The existence of stationary clouds is only allowed for specific values of the background parameters M and a. Such quantization follows from the regularity of the bound states and results in an existence line in the two-dimensional Kerr parameter space defined by (M, a) or, alternatively, (M, Ω H ).
The next subsections summarise the algorithm to solve the radial equation (2.18) together with the synchronization condition and therefore determine the existence lines of stationary vector clouds around Kerr BHs. For convenience, Eqs. (2.18) and (2.19) will be considered as written in terms of r H and a instead of M and a. For this purpose, using the identity the function ∆ may be written as (2.22)

Angular equation
In the Minkowski limit, the phase angular velocity equals the inverse of the reduced Compton wavelength: ω = µ = 1. Hence, Eq. (2.19) can be written in the form 4 where Q 0 denotes the leading-order form of the function Q at spatial infinity and the superscript 'E' will become clear in the remainder of the present section.Ĵ 2 coincides with the square of the orbital angular momentum operator, (2.24) whose eigenfunctions are the well-known scalar spherical harmonics of degree j and order Since the Kerr spacetime is asymptotically flat, j may be defined as the total angular momentum at spatial infinity.
This means that, at leading order, the electromagnetic four-potential A α = (A t , A) takes the form is the leading-order form of the function R at spatial infinity. The spatial part of A α can be written as The angular dependence of A 0 is described by the 'pure-orbital' vector spherical harmonics in flat space. Equation (2.26) with eigenvalues λ E 0,∓ , corresponds to the j = ± 1 electrictype states of the vector field (cf. Subsection 2.1). This explains the superscript 'E'.
In the zero-angular-momentum (ZAMO) frame, characterized by the tetrad where Ξ ≡ (r 2 + a 2 ) 2 − a 2 ∆ sin 2 θ, the electric field E and the magnetic field B have the following components: is the four-dimensional Kronecker delta. The leading-order terms of E (a) and B (a) are given by where we used the fact thatê • (a) →ê (a) as r → +∞. The electric field E 0 depends on the difference j − , whereas the magnetic field B 0 does not.
When α 1, the angular eigenstates and eigenvalues for the electric-type states may be written as an expansion in α. The next-to-leading order corrections to the angular eigenfunctions induce couplings to vector spherical harmonics of the same parity and thus the angular functions Q have definite parity. For future reference, we present the expansion for the angular eigenvalues below [37]: where the first terms of the series are given by withn ≡ n + + 1.n ∈ N is the principal quantum number and n ∈ N 0 is the node number -see Subsubsection 2.3.3. Equation (2.23) allows us to recover the electric-type states of A a solely. The magnetictype states with j = = |m j |, the only ones which are known to be captured by the FKKS ansatz, can be recovered considering the limits where the superscript 'M' labels all quantities related to magnetic-type states with j = |m j |.
As first shown in [37], the leading-order form of A is proportional to the vector spherical harmonic Y j j,j . Unluckily, no expansion of λ M in powers of α is known. However, when considering marginally-bound states (ω 2 = µ 2 = 1), the angular eigenvalue yields

Radial equation
The integration of the radial equation is performed via the expansion for the radial function R in Eq. (2.18). N is the number of terms of the partial sum and the coefficients {c n } n=0,...,N are functions of 5 r H , a, µ, m j and λ, which, in turn, depends on and j. Plugging the expansion into Eq. (2.18) and equating coefficients order by order, it is possible to write {c n } n=1,...,N in terms of c 0 . The latter is usually set to 1. The choice of N should be a trade-off between computational time and accuracy. Once the coefficients {c n } n=0,...,N are defined, one fixes the numerical values of r H , µ, , j and m j , assigns a guess value to a and computes the corresponding guess value for λ.
The radial equation is then integrated from r = r H (1 + δ), with δ 1, to r = r ∞ , where r ∞ stands for the numerical value of infinity. The solution must satisfy the boundary conditions where the prime denotes differentiation with respect to r. The previous step is repeated for different guess values of a, until the solution satisfies the boundary conditions

Results
When scanning the parameter space in search of stationary vector clouds with fixed quantum numbers ( , j, m j ), solutions with different numbers of nodes n (n ∈ N 0 ) are found. Thus, each vector state may be labelled using the notation |n, , j, m j . Configurations with n = 0 (n ∈ N) are dubbed fundamental (excited) states. The greater the node number n, the more energetic the state. The frequency spectra of massive vector quasi-bound states, i.e. with a complex frequency, can be written in the form .
The frequencies and corresponding instability rates were computed analytically in [37] via matched asymptotic expansions, except for the magnetic-type (j = ) vector states. The expression in Eq. (2.35) for j = is a conjecture. Nonetheless, the authors of [37] confirmed that the conjectured frequencies do agree with those found numerically without relying on separability of the Proca equation. In fact, the analytic approximation is accurate when α 0.2, even for near-extremal Kerr BHs. Vector instability rates are proportional to the factor (ω − m j Ω H ) and thus vanish whenever the synchronization condition holds. In that case, the contour lines for which |n,j,m j denotes the frequency of a massive scalar quasibound state with quantum numbers {n, j, mj}. This suggests that the magnetic-type vector states are equivalent to the scalar states with the same total angular momentum. If so, it should be possible to show that Eq. (2.18) for magnetic-type states and its scalar counterpart are equivalent, at least in some limiting case.
constitute an analytical approximation to the existence lines of stationary vector clouds in the parameter space of Kerr BHs. For future reference, theses curves will be referred to as analytical existence lines (AEL), whereas those obtained via the numerical algorithm laid out in Subsubsection 2.3.2 will be named numerical existence lines (NEL). Additionally, all existence lines will be presented in a (M, Ω H )-plane normalized to the boson's mass µ, in which the domain of existence of Kerr BHs is shaded light green.
The mass spectrum of Kerr BHs which support stationary vector clouds may be derived by solving Eq. (2.36) for α. This yields The first two terms in Eq. (2.37) depend on n, and m j , but not on j. The next-toleading-order term, which depends on j through g V , must be taken into account to capture the leading-order behavior of stationary vector clouds. The existence lines for the vector states |0, , 1, 1 with 7 ∈ {0, 1, 2} are shown in Figure 1 (top panel). When m j = 1, the line corresponding to the lowest values of Ω H belongs to the electric-type state |0, 0, 1, 1 , which is therefore the fundamental mode with m j = 1. The analytical existence line for the electric-type state |0, 0, 1, 1 is in agreement with its numerical counterpart when α 1. As α increases to values near the extremal case a = M (black solid line), the two lines diverge from each other. This behavior appears to be a generic feature of existence lines corresponding to states for which j = m j and < j (see Figure 1 -bottom panel). The discrepancy, whose source remains unclear, suggests that higher-order corrections to the α-expansion in Eq. (2.35) are needed when describing clouds around rapidly-rotating Kerr BHs. On the other hand, the analytical and numerical existence lines for the electric-type state |0, 2, 1, 1 , for which j = m j but > j, appear to overlap over the full range of α.
When {n, j, m j } are fixed, the existence lines move towards greater values of Ω H as the orbital angular momentum increases. In fact, the larger the value of , i.e. the greater the energy of the state, the greater must the angular velocity Ω H be for stationary equilibrium. Moreover, the existence lines converge in the limit of vanishing mass, i.e. M → 0, which reflects the fact that the spacetime becomes insensible to the cloud's features. These trends were also found for stationary scalar clouds around Kerr BHs [21].
The variation of the node number n when { , j, m j } are fixed yields identical behavior. The existence lines for the vector states |n, 0, 1, 1 with n ∈ {0, 1, 2} are plotted in Figure 1 (bottom panel). The node number plays a similar role to that of the principal quantum number in the description of hydrogen's energy levels: the larger the node number n, the more energetic the state. Given two existence lines with the same { , j, m j }, the one with  the largest node number n lies to the right with respect to other in the (M, Ω H )-plane. Additionally, they converge in the limit of vanishing M .
The radial profile of the clouds I, II and III in Figure 1 (bottom panel) are displayed in Figure 2 (top panel). The function R is finite over the whole r domain outside the event horizon and vanishes (exponentially) as r → +∞, as required by asymptotic flatness. Besides, the local maximum closest to the event horizon decreases with increasing n. Finally, Figure 2 (bottom panel) shows the dependence of the radius of the cloud, hereafter denoted by r C , on the rotation parameter a for different vector stationary clouds with n = 0. r C is defined as the value of r closest to r H that locally maximizes the function 4π|R| 2 . Its value diverges in the Schwarzschild limit (a → 0), in accordance with the fact that Schwarzschild BHs cannot carry stationary vector clouds. Moreover, the minimum of r C , which occurs at a = M , is finite, which means that Kerr BHs do not support sufficiently tight clouds. Similar observations were already reported for stationary scalar clouds in [21].

Non-linear analysis: hairy black holes and Proca stars
We now address the fully non-linear solutions of the Einstein-complex-Proca model, described by the action where L M is the Proca Lagrangian density (2.4). Varying this action one obtains the Proca equations (2.5) and the Einstein equations where the Proca energy-momentum tensor is: We follow the conventions of [10]. More details on the formalism can be found therein. If one linearises the model (3.1) in the Proca field, one ends up with the vacuum Einstein equations and a test Proca field on a fixed curved background (that solves the vacuum Einstein equations). This corresponds precisely to the analysis of Section 2.

The ansatz
To find the hairy BHs that bifurcate from the linear clouds that were studied in Section 2 we use the metric ansatz 8 and F i , W are functions of the spheroidal coordinates (r, θ). The parameter r H is the radial coordinate of the event horizon, which is θ-independent. For the Proca potential, we use an ansatz that depends on four functions (V, H a ). All these functions depend on (r, θ). The ansatz has a harmonic time and azimuthal dependence, which introduces a (positive) frequency, ω > 0, and the azimuthal harmonic index, m ∈ Z: A = e i(mϕ−ωt) (iV dt + H 1 dr + H 2 dθ + iH 3 sin θdϕ) . Here, m should be identified with m j of Section 2 and we shall focus on m = 1. We follow closely [10], wherein all details can be found, namely: the explicit equations of motion for this ansatz (in Appendix B therein) and the boundary conditions at the horizon, spatial infinity and on the axis (in Section 4 therein). Details on the numerical method can be found in Section 3.3 of [45]. The key feature for the existence of these BHs is the synchronisation condition (1.3), where Ω H = W (r H ), the non-diagonal metric function in Eq. (3.4), which on the horizon is independent of θ, and m, ω are the parameters in the Proca ansatz (3.6).

Domain of existence
When finding solutions via a relaxation method, such as the Newton-Raphson method used for this work, the initial guess plays a key role to guarantee convergence to the desired solutions. In [45], the construction of hairy BHs started from the spinning Proca stars in [16], which have one node for the temporal component of the Proca potential, V . Consequently, the hairy BHs reported in [45] also have one node in V . At that point, it was found no evidence for nodeless solutions of either spinning Proca stars or BHs with Proca hair, even though it was stated in [45] that no proof for the inexistence of nodeless solutions could be established (except for spherical Proca stars). These results were reconsidered after the numerical evolutions of the Kerr superradiant instability have been reported [7]. The data describing the equilibrium points attained in these evolutions matched spinning BHs with Proca hair and a nodeless Proca potential temporal component V , first constructed in [8], wherein their domain of existence was exhibited. This domain of existence is shown in Figure 3, together with the domain of existence of the nodeful (n = 1) solutions reported in [45]. The hairy BHs exist in the blue shaded regions. In each case (n = 0 or n = 1) the domain of existence is bounded by the solitonic limit (red solid lines) wherein the hairy BHs become spinning Proca stars with the same n, and by the bald limit (blue dotted lines), wherein they meet the Kerr parameter space at the corresponding existence line, with m j = m, the same n and ( , j) = (0, m j ). Thus, the two blue dotted lines plotted in Figure 3 correspond to the blue (n = 0) and yellow (n = 1) numerical existence lines plotted in Figure 1 (bottom panel).  The existence line from which the fundamental non-linear solutions bifurcate follows a similar rationale to that observed for the scalar case [9]. For a given m = m j , the existence line with = 0 and j = m j is the leftmost one in the Kerr parameter space plotted in Figure 1. Thus it represents the threshold between the Kerr BHs that are stable against all modes with that m j and the ones that are unstable against at least one such mode. Since m j is the only of the three quantum numbers ( , j, m j ) that remains significant in the non-linear theory -it is associated to an isometry -, for each m j the existence line whence the hairy BHs bifurcate is the one with ( , j) = (0, m j ). BHs emerging from the other existence lines with m = 1 are likely to exist but are excited states, with either more radial or angular nodes. Inspection of Figure 3 reveals two main features. Firstly, as expected, the excited states (n = 1) can attain a larger ADM mass; secondly, the fundamental states of the BHs with Proca hair exist for a larger ω-range; this also seems intuitive: excited states require a larger minimum angular velocity. 9 The same trends are observed when comparing the fundamental states (n, m) = (0, 1) of the scalar and the Proca hairy BHs - Figure 4, with the scalar case playing the role of the excited Proca family, in this comparison. This had already been observed for the solitonic limit in [48]. In the bald limit, this means that the fundamental Proca existence line spans lower Ω H BHs. This is a manifestation of the well known fact that the superradiant instability is stronger for the vector case [5].  Unlike the scalar or the excited Proca case, in the case of n = 0 spinning Proca stars, that compose the (red solid line) boundary of the domain of existence, it was not possible to explore the domain of solutions after the backbending, i.e. when the minimum frequency is attained -see inset in Figure 4. The reason is that these solutions become rather compact and hence strong gravity configurations, making their computation numerically challenging. To assess this, we have used the same measure of compactness as, e.g. in [45,58], namely: where R 99 is the perimetral radius that contains 99% of the star's mass, M 99 . We recall that bosonic stars do not have a surface where a discontinuity of the energy density occurs; rather, they decay exponentially, vanishing only at infinity. The perimetral radius is a geometrically meaningful radial coordinate R: a circumference along the equatorial plane has perimeter 2πR. The inverse compactness of the n = 0, 1 Proca stars and n = 0 scalar boson stars, all with m = 1, is shown in Figure 5. One observes that the inverse compactness is always greater than unity, meaning that all these stars are less compact than a BH. Moreover, the fundamental Proca stars become the most compact ones, precisely at the backbending, where they attain an inverse compactness 1.1. Another token of strong field gravity is the formation of ergo-regions in spinning spacetimes. In the solitonic limit, both the scalar and vector spinning stars do not have ergoregions when ω → µ (see Figures 3 and 4), corresponding to the dilute regime where the stars are not compact and not strongly relativistic. Moving along the spiral and away from this dilute regime, in all cases the ergo-region appears in the first branch, i.e., before the first backbending and for quite compact stars. The first occurrence of an ergo-region along the sequence of bosonic stars is marked with a square in Figure 5. The comparison between the three different cases shows that compactness is not the only factor determining the existence of an ergo-region. For all three cases (Proca with n = 0, 1 and scalar with n = 0) the ergo-region of these stars is toroidal. In the family of the hairy BHs, this toroidal region adds up to the ergo-sphere around the spinning horizon. We have not scanned in detail the parameter space but one will get a rich ergo-region structure, including ergo-Saturns, analogous to those found for BHs with synchronised scalar hair (n = 0) [59], Proca hair (n = 1) [10] and other cousin models, e.g. [34,35,60,61].

Analysis of specific solutions
In order to get a better intuition on the impact of the node number on the solutions let us consider a comparative study between the profile functions of two spinning Proca stars, one with n = 0 and another with n = 1, and both with the same frequency ω/µ = 0.9. These two configurations are highlighted as two stars in Figure 3.
In Figure 6 we compare the metric functions of the two illustrative Proca stars in terms of a compactified radial coordinate, to have an overview of the whole radial domain, and for three different θ-values. Whereas in the fundamental state all metric functions are rather smooth and monotonic, in the first excited state there is some extra structure, mostly noticeable along the equatorial plane (θ = π/2). The metric function W , in particular, is no longer monotonic. The nodeless vs. nodeful structure of the n = 0 vs. n = 1 spinning Proca stars becomes evident in Figure 7. One observes, in particular, that all four Proca potential functions have the same number of nodes, n = 0 or n = 1, for each star. Moreover, the temporal and radial component of the potential have a trivial structure along the θ = 0 symmetry axis. Finally, the extra structure of the excited states becomes clear when analysing more invariant quantities, such as the Noether charge density, the Ricci curvature scalar and the Komar energy density, that are exhibited in Figure 8.
The Noether charge results from the global U (1) symmetry of Eq. (3.1), which is invariant under the global transformation A β → e iχ A β , where χ is a constant. Thus, a conserved 4-current exists The Noether charge, Q, which is interpreted as the particle number (indeed becomes the particle number upon quantisation), is obtained integrating the time component of this current on a spacelike hypersurface Σ: The Nother charge density is thus J t , which is plotted in the top panel of Figure 8. The Komar energy density results from the Komar mass computed at infinity. Using Gauss's law, one relates the latter with the horizon Komar mass M H (in the cases which have a horizon) and a volume integral on a spacelike hypersurface between the horizon and infinity. One obtains [45] (where k α is the asymptotic timelike Killing vector field):  Figure 7: Same as in Figure 6, but for the Proca potential functions. where M (P) is the energy contained in the Proca field (outside a horizon, in case there is one): The integrand is the Proca energy density, which is plotted in the second from bottom panels of Figure 8. A similar analysis can be done for the Komar angular momentum density, showing that the Komar angular momentum density is T t ϕ [45], plotted in the bottom panels of Figure 8.
All invariant quantities in Figure 8 demonstrate that whereas the n = 1 stars have a Saturn-like morphology -which was observed in [45] -, with the energy density or the particle number having a global maximum at the centre and a local maximum at some radial distance, the fundamental states are spheroidal. This contrasts with the toroidal shape of the fundamental spinning scalar boson stars [14]. This morphological difference was argued to be related to the different dynamical stability of the fundamental states of scalar/vector spinning bosonic stars [17].
We now turn to hairy BHs. The single most important observation concerning hairy BHs in this model is that they can be quite Kerr-like or strongly non-Kerr-like. This follows from the fact that the hairy BHs interpolate between the Kerr family and a solitonic limit (Proca stars) whose properties and phenomenology can be quite different from Kerr. So, here we shall focus on two illustrative solutions that exemplify this range of possibilities.
First, we consider an example of a fairly Kerr-like BH with Proca hair, with n = 0. It is chosen in the region where these BHs matched the endpoint of the dynamical of evolutions reported in [7] -see [8]. Moreover, within this region, it is chosen to be as hairy as those evolutions suggest a hairy BH can be, when forming dynamically from the superradiant instability of Kerr BHs. This hairy BH, labelled HBH 1 , has 10 [all quantities in Eqs. (3.12) and (3.13)  H are the dimensionless spin in terms of global and horizon quantities, respectively. Thus, HBH 1 has 9.5% of its energy and 39.3% of its spin outside the horizon. These were roughly the maximal values of extraction via superradiance observed in [7]. Also note that both j and j H are smaller than unity; thus the hairy BH obeys the Kerr bound, both in terms of horizon and asymptotic quantities. It is known that spinning BHs with synchronised hair can violate the Kerr bound -see e.g. [9,62].
Second, we consider an example of a fairly non-Kerr-like BH with Proca hair, with n = 0. This hairy BH, labelled HBH 2 , has (3.13) 10 Do not confuse the dimensionless spin j in this Section with the total angular momentum j of Section 2.
Thus, HBH 2 has 76.9% of its energy and 97.8% of its spin outside the horizon. Moreover, this BH violates the Kerr bound in terms of asymptotic quantities, since j = 1.56 > 1, but not in terms of horizon quantities. In this sense it behaves more like a star. Both these solutions are marked with triangles, and labelled with the corresponding numbers, in Figure 3.
In Figure 9 we exhibit the metric functions outside the horizon for the two hairy BHs. In the case of HBH 1 , the insets show a comparable Kerr BH, that is with the same total mass M and angular momentum J. In order to make the latter comparison, the Kerr metric is expressed in the gauge (3.4) -see Appendix A in [45]. No comparable Kerr BH exists for HBH 2 , as the latter violates the Kerr bound. Again we use a compactified radial coordinate. The metric functions of HBH 1 already show some qualitative differences relatively to those of the comparable Kerr BH. The latter has a considerably higher Ω H /µ 1.72; indeed HBH 1 and the comparable Kerr roughly correspond to the two points in Fig. 5 of [8], representing the longest migration. On the other hand, the metric functions of HBH 2 are a hybrid between the Kerr metric functions and those of a Proca star -see Figure 6 (left panels). Indeed, as can be seen from the physical parameters in Eq. (3.13), HBH 2 has over three quarters of the total mass and almost the totality of the angular momentum stored in the Proca field outside the horizon. Thus, it is more accurately described as a spinning Proca star with a BH horizon at its centre, than as a BH horizon surrounded by a Proca cloud. The latter is an appropriate description for HBH 1 .
In Figure 10 the Proca potentials are shown for the two hairy BHs. One can appreciate the difference in boundary conditions as compared to the Proca stars in Figure 7. V, H 2 , H 3 are non-zero on the BH horizon and zero at the origin, for stars; H 1 is the opposite. On the other hand, the most apparent differences between HBH 1 and HBH 2 are the larger magnitude of the Proca potential functions for the latter, together with a steeper behaviour. This is intuitive from the fact the second BH has a much larger fraction of its energy in the Proca field.
Finally, in Figure 11 we represent some physical quantities of the two hairy BH solutions. The Noether charge density is one order of magnitude larger for HBH 2 and with a steeper profile. This impacts on the Ricci scalar curvature, known to manifest the spacetime deformation due to matter, which has a clear lump outside the horizon for the hairiest solution. The Komar energy and angular momentum densities are also larger in magnitude and with sharper profiles, becoming asymptotically more similar to those of the Proca star exhibited in Figure 8 (left panels).

Conclusion
In this paper we have analysed linear vector clouds of a massive Proca field around a Kerr BH and the BHs with synchronised Proca hair that can be considered as the non-linear realisations of these clouds. Our analysis has been inspired by a series of fairly recent developments that motivates revisiting the Einstein-(complex)-Proca model and its BH and solitonic solutions. Notice, however, that the linear analysis in Section 2 does not depend on the fact that the Proca field is complex, unlike the analysis in Section 3, where  Figure 11: Some physical quantities for HBH 1 (left panel) and HBH 2 (right panel).
the existence of the stationary solutions describing hairy BHs and Proca stars relies on the field being complex.
Concerning the linear analysis of Section 2, the key physical property of these boundstate configurations is the synchronization of their phase angular velocity with the event horizon angular velocity. Furthermore, they resemble the hydrogen's atomic orbitals and can be described in terms of {n, , j, m j }. The quantum numbers label the existence lines of stationary vector clouds in the two-dimensional parameter space of Kerr BHs. These curves mark the bifurcation of the Kerr family towards the new family of BHs with Proca hair and constitute one of the boundaries of the domain of existence of the latter [10].
As for massive scalar bosons [21], the analysis of the vector clouds shows that, for a fixed value of the azimuthal total angular momentum m j , the cloud's energy, which is proportional to its phase angular velocity, is mainly determined by the node number n, and the orbital angular momentum . The bound-state |0, 0, 1, 1 has the lowest possible energy and higher values of n and/or correspond to higher-energy states. Thus, the existence line of this bound state is wherein the fundamental states of the hairy BHs bifurcate from. Moreover, despite not having a relevant impact on the cloud energy, the total angular momentum j allows for the existence of = 0 bound states, which is rooted in the nonvanishing intrinsic angular momentum of the bosons.
The existence lines obtained numerically were compared with analytical approximations recently reported in the literature [37]. In general, the agreement is excellent for all values of the Kerr BH's rotation parameter, except when j = m j and < j. In this case, a discrepancy arises for near-extremal Kerr BHs; the reason behind this observation remains to be clarified.
The analysis' starting point was the FKKS ansatz [42] for the separation of the Proca equation. This ansatz prevents the need for approximations or time-consuming numerical algorithms when studying massive vector bosons in Kerr-NUT-(A)dS spacetimes and has already been used to address quasi-bound states in the Kerr and Kerr-Newman backgrounds [37,[42][43][44]63]. This ansatz, however, does decouple and separate the torsionmodified Proca equation (known as Troca equation) in the Chong-Cvetič-Lü-Pope spacetime of D = 5 minimal gauged supergravity [44] and in the Kerr-Sen spacetime of lowenergy heterotic string theory [44]. Note that All these works focused on the dynamics of massive vector bosons in the frequency domain, in which the particles are described as monochromatic waves. Future research should then dive into a yet-to-be-explored timedomain analysis of these simplified equations of motion. Of particular interest would be to perform long-time evolutions of massive vector Gaussian wave packets in superradianceprone spacetimes.
Concerning the non-linear analysis in Section 3, here we have exhibited the domain of existence of the fundamental states of BHs with synchronised Proca hair and compared some of their properties with the first excited states, discussed in [10], and the cousin hairy BHs obtained in the scalar case [9,45]. Then, we have analysed some illustrative solutions of both the solitonic limit (Proca stars) and hairy BHs. We emphasise that all the solutions considered here have azimuthal harmonic index m (≡ m j ) = 1. Higher m solutions also exist, corresponding to another sort of excitation. 11 There are two main ideas to retain from our results.
Firstly, there are morphological differences between the cases compared herein, which may have various implications. This is summarised in Figure 12, where surfaces of constant scalar or Proca density for stars and hairy BHs are exhibited for the three families of solutions we have compared. Spinning scalar bosonic stars (n = 0) are toroidal; spinning vector stars with n = 0 are spheroidal and with n = 1 have a Saturn-like morphology. The corresponding hairy BHs are a non-linear bound state of such a star with a horizon, deforming the star's morphology. Figure 12: Morphology of the surfaces of constant energy density of scalar (n = 0) and vector/Proca (n = 0 and n = 1) spinning bosonic stars (top panels) and BHs with synchronised hair of the corresponding type (bottom panels). For the BHs, the black disk represents the horizon. All these solutions have m = 1. The hairy BHs bifurcate from clouds with = 0. Besides further radial excitations (with higher n), these families of solutions possess azimuthal excitations (with higher m) and also further polar-angular excitations (corresponding, in the linear limit, to higher ).
Secondly, in all these families of hairy BHs, in particular in the n = 0 BHs with Proca hair, there are Kerr-like solutions, but also rather non-Kerr-like examples. At the moment, at least one formation channel for the hairy Kerr-like solutions is known -superradiance; HBH 1 discussed in Subsection 3.3 belongs to this set. The solution shows already some interesting deviations from Kerr and it will be very interesting to analyse how these deviations impact on astrophysical observables. In this respect, one could reconsider some of the 11 Higher m increases the number of nodes in the azimuthal direction. An in-depth study of the higher m solutions in the scalar case is found in [34]. Some particular higher m solutions in the Proca case can be found in [17]. We remark that n = 0 spinning Proca stars with m > 1 possess surfaces of constant energy density with a toroidal morphology.
analysis done for the scalar case or for the excited BHs with Proca hair, namely of shadows [64,65], X-ray spectroscopy [66,67] or quasi-periodic oscillations [68]. Of course, one of the most interesting open questions concerns the dynamical properties of these BHs, including quasi-normal modes. The recently established dynamical robustness of spinning Proca stars [17] has paved the way to perform dynamical evolutions of these BHs, from which one could, in particular, extract waveforms for binary evolutions. Work in this direction is underway.

A Vector spherical harmonics
The components of the first few ('pure-orbital') vector spherical harmonics Y j,m j in terms of the spherical unit vectors {ê (r) ,ê (θ) ,ê (ϕ) } are Their real part are displayed in Figure 13.