Plausible scenario for a generic violation of the weak cosmic censorship conjecture in asymptotically flat four dimensions
Felicity C. Eperon Bogdan Ganchev Jorge E. Santos

	



fce21@cam.ac.uk
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK 

bvg25@cam.ac.uk
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK 

jss55@cam.ac.uk
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK 


	We present a plausible counterexample to the weak cosmic censorship conjecture in the four-dimensional Einstein-Scalar theory with asymptotically flat boundary conditions. Our setup stems from the analysis of the massive Klein-Gordon equation on a fixed Kerr black hole background. In particular, we construct the quasinormal spectrum numerically, and analytically in the WKB approximation, then go on to compute its backreaction on the Kerr geometry. In the regime of parameters where the analytic and numerical techniques overlap we find perfect agreement. We give strong evidence for the existence of a nonlinear instability at late times.


Plausible scenario for a generic violation of the weak cosmic censorship conjecture in asymptotically flat four dimensions
    Jorge E. Santos
    December 30, 2023
==========================================================================================================================



 Introduction – Quantum gravity remains terra incognita, largely because it is hard to access experimentally. One might wonder why, since singularities are known to form under a variety of circumstances <cit.>. However, general relativity stubbornly conceals these regions of high curvature, where quantum gravity is likely to play a leading role, by hiding them behind an absolute event horizon. This phenomenon is observed for such large classes of initial data that it was promoted to a conjecture in <cit.>. This is the weak cosmic censorship conjecture (WCCC), which forbids the formation of naked singularities, i.e. singularities in causal contact with future null infinity, starting from generic initial data.

Here we propose a four-dimensional asymptotically flat counterexample to the WCCC motivated by the superradiant instability <cit.>, afflicting massive scalar perturbations around Kerr black holes (BHs) <cit.> - rotating, with spherical topology and considered the most general BH solutions of the vacuum Einstein equation <cit.>.

The angular dependence of such perturbations around Kerr is parametrized by two integers {m,ℓ}: m counts the number of nodes in the azimuthal direction and ℓ-|m| the number of zeros in the polar direction. We show that for any scalar field of mass μ and any nonzero value of the BH spin, for sufficiently large values of ℓ=m,  these perturbations herald instabilities around the BH, extracting energy and angular momentum. Furthermore, the time scales associated to these instabilities grow parametrically as e^4 ℓ logℓ, indicating that each of the ℓ modes decouples from the rest, evolving independently. 

As time progresses, modes with smaller values of ℓ stabilize one by one, forming scalar clouds around the BH, similar to those in <cit.>. However, these BHs were shown to be unstable to higher m-modes <cit.>, giving rise to the expectation of a cascade towards larger values of ℓ. This corresponds to a transfer of energy from lower ℓ-modes to higher ones, indicating an evolution towards smaller scales - a phenomenon akin to turbulence in nonrelativistic 3+1 fluids.

A possible stabilizing mechanism is the emission of gravitational waves (GWs) by the scalar clouds <cit.>. Were they to dissipate energy faster than superradiance creates them, the above scenario would not be possible. We numerically compute the GW emission for fixed gravitational coupling, M μ, and spin parameter, J/M^2, as a function of ℓ=m and find that it leads to energy and angular momentum dispersion that would not be able to counter the efficiency of superradiance.

Our paper is organized as follows: first we present our setup and provide analytic and numerical data for the instability time scales at large ℓ. We then compute, numerically, the energy radiated towards future null infinity in this process as well as the backreaction of the scalar on one of the components of the Weyl tensor. We see that modes with higher ℓ radiate less, implying that energy is accumulated at small scales more efficiently for larger values of ℓ. Finally we end with discussion of the results.



 Setup – We work with the Einstein-Hilbert action minimally coupled to a real massive scalar field ψ

    S=∫_ℳd^4x√(-g)(R/16 π G-∇_aψ∇^aψ-μ^2ψ^2),

where μ is the scalar field mass, g_ab the spacetime metric and R its Ricci scalar. The equations of motion are

	
    R_ab-R/2g_ab= 8π G T_ab

	where R_ab is the Ricci tensor of g_ab,
	
    ψ =μ^2 ψ ,

	and
	
    T_ab=2∇_a ψ∇_b ψ-g_ab∇^c ψ∇_c ψ-μ^2ψ^2 g_ab .

	

An important solution to these equations is the Kerr BH <cit.>, with ψ=0 and

    ds^2 = -Δ/Σ^2(dt-a sin^2θdϕ)^2+Σ^2(dθ^2+dr^2/Δ)
    
    
    +sin^2θ/Σ^2[a dt-(r^2+a^2)dϕ]^2 ,

where Δ = r^2+a^2-2 M r, Σ^2=r^2+a^2cos^2θ, ϕ∈(0,2π) is a periodic coordinate and θ∈(0,π) is a polar coordinate.  The BH event horizon is a null hypersurface with r=r_+≡ M+√(M^2-a^2), angular velocity Ω_K = a/(a^2+r_+^2) and surface gravity κ_K=(r_+^2-a^2)/[2 r_+(r_+^2+a^2)]. The constant M is the BH mass and a parametrizes its angular momentum via J = M a. The absence of naked singularities demands |a|≤ M with the inequality being saturated at extremality, when the Kerr BH event horizon becomes degenerate with κ_K=0.

We study Eq. (<ref>) on a fixed Kerr background (<ref>), which is stationary and axisymmetric with respect to the Killing fields ∂ /∂ t and ∂/∂ϕ, respectively. We consider perturbations of the form

    ψ(t,r,θ,ϕ) = e^-i ω t+i m ϕψ_ω m(r,θ)

and assume that ψ_ω m(r,θ) is separable, ψ_ω m(r,θ)=R_ωℓ m(r)S_ωℓ m(θ), with the label ℓ anticipating that the separation constant will be parametrized by an integer ℓ. Not all solutions to Eq. (<ref>) are separable, but we are interested in those composed of the sum (possibly infinite) of such separable solutions. ω∈ℂ is a complex frequency, determined by imposing appropriate boundary conditions. We are interested in finding unstable mode solutions for which Im(ω)>0.

Inserting the ansatz (<ref>) into Eq. (<ref>) yields a system of two second order ordinary differential equations, coupled via the separation constant Λ,

	
    Δ [Δ R_ωℓ m_,r]_,r+V(r)R_ωℓ m=0,

	
    1/sinθ[sinθ S_ωℓ m_,θ]_,θ
    
    	-[a^2 k^2cos^2θ+m^2/sin^2θ-Λ]S_ωℓ m=0,

	where
	
    V(r)   =-k^2r^4+2M μ ^2 r^3-(Λ+a^2 k^2)r^2+
       (2MΛ -4a m M ω+2Ma^2 ω^2)r-a^2(Λ-m^2) ,


with k≡√(μ^2-ω^2). Finding unstable modes amounts to finding the values of ω for which ψ has ingoing boundary conditions at the event horizon (consistent with the equivalence principle) and finite energy on a partial Cauchy surface t= const.

This problem can be tackled numerically (for any values of the parameters) and analytically (in certain regions of moduli space). We will first compute the modes using a WKB expansion in m, which we detail next.



 WKB Expansion and Numerical Validation– Our WKB expansion is valid for any spin parameter |a|<M, and only assumes m to be large. For small values of a and μ, it reproduces the results in <cit.>[The large m limit does not commute with the extremal limit. It can be shown that at any finite value of m, Im(ω)=0 at extremality. However, the region of the black hole moduli space where this behaviour occurs scales inversely with m, and is thus absent in the strict m→+∞ limit. The behaviour in this region could be investigated by switching the order of limits.] (up to a factor of 2, see <cit.>). For large m, Eq. (<ref>) can be approximated by the usual equation for spherical harmonics on a 2-sphere so that Λ=ℓ(ℓ+1)+O(ℓ^-1) and ℓ, m are integers with ℓ≥0, |m|≤ℓ. In the following we set m=ℓ and take the limit ℓ≫ 1.

To determine the large ℓ limit at fixed μ and a, we used the same method as in <cit.>, which combines a matched asymptotic and WKB type approach. We quote the final result and leave the details to the Supplemental Material,

	
    Re(ω M)=μ̂(1-μ̂^2/2ℓ^2)+𝒪(ℓ^-3) ,

	
    Im(ω M)=ℓ^-4 ℓ-9/2+p /2^2 ℓ+1-p √(π) p ! μ̂^4 ℓ+5sinh[π(ℓΩ _K-μ)/κ _K]×
    exp[-2/κ _K(ℓΩ _K-μ̂/r_+) arctan(Ω_K/κ_K)-2 (1-ℓ+p )]×
    [1+𝒪(ℓ^-1)] ,


where μ̂≡μ M and p∈ℕ_0 is a radial overtone. First note that if we set a=0, the argument of the sinh becomes negative, and the instability disappears. Furthermore, its onset sits precisely at the onset of superradiance, namely ℓΩ_K=Re(ω). More importantly for our purposes, in the limit ℓ→+∞, the growth rate of the instability scales as e^-4 ℓlogℓ, and no matter what the value of μ or a, one can always find a value of ℓ=ℓ_⋆≡⌈μ/Ω_K⌉ above which the instability sets in. This shows that all Kerr black holes are unstable to massive scalar field perturbations, irrespective of their initial spin |J|<M^2 and of the mass μ of the scalar perturbation.

One can test the regime of validity of our approximation, by comparing with data obtained by solving numerically, without approximations, the full equations. To this end, we perform a change of variables so that we only solve for smooth functions in their integration domain. This necessarily involves a choice of normalization, which we describe next. For numerical aid we define

	
    R_ωℓ m(r) = (1-r_+/r)^-i ω-m Ω_K/2 κ_Ke^-k r(r_+/r)^γ q_r(r) ,

	
    S_ωℓ m(θ) = sin^mθ q_θ(θ) ,

	

with q_r(r_+)=q_θ(0)=1 and γ≡ 1+2 M k-M μ ^2/k. We now use the methods of <cit.> to solve for the eigenpair (ω,Λ) using a Newton-Raphson routine [The only additional complication is that we are searching for extremely small growth rates when m increases, so using extended precision is mandatory.].

As seen in Fig. <ref>, our numerical data agrees excellently with (<ref>). Furthermore, one can measure deviations of our WKB expression to the exact numerical result, and it agrees with the error given in Eq. (<ref>).




 Backreaction – We want the GWs emitted by a scalar cloud around a Kerr BH and its leading order backreaction on the geometry. In the vector field case <cit.> it has been shown that the system evolves adiabatically; the emergence of the cloud due to superradiance, and the consecutive saturation of the vector mode responsible, due to the spinning down of the BH, proceed on a much faster time scale than the dispersion of energy and angular momentum due to GW emission from the cloud.

We proceed using nonlinear perturbation theory and declare

    ψ=∑_i=0^+∞ψ^(2i+1)ε^2i+1 ,  and   g=g_K+∑_i=1^+∞ g^(2i)ε^2i ,

where g_K is given by the Kerr metric (<ref>). We expand the equations of motion (<ref>) in a power series in ε. To first order in ε we solve eq. (<ref>) subject to a choice of initial data. For the case at hand, we choose ψ to be given by the real part of one of the unstable modes we have determined above. These are labelled by a given value of m. Furthermore, since  Im(ω M)≪Re(ω M), we take ω to be purely real. We then proceed to second order and attempt to compute the leading order backreaction on the metric, g^(2), and its associated curvature. The standard approach to the linearized Einstein equation presents us with a daunting task; however, Kerr BHs are algebraically special, allowing us to bypass computing g^(2), and directly calculate certain gauge invariant scalars built out of the Weyl tensor. These do not couple amongst themselves and we focus on the Newman-Penrose scalar ψ_4 since it also allows us to efficiently compute the GWs emitted by the scalar cloud. ψ_4 obeys the Teukolsky equation <cit.>

    [(Δ+3γ-γ̅+4μ+μ̅)(D+   4ϵ-ρ)-3ψ_2-
    
    -(δ̅+3α+β̅+4π-τ̅)(δ+   4β-τ)]ψ_4=4π T_4 ,

whereby all quantities appearing in Eq. (<ref>) are given in the Supplemental Material.  The source term T_4 is reconstructed directly from Eq. (<ref>) and thus depends on ψ^(1) and its gradient only.

The lhs of (<ref>) can be separated into angular and radial parts as in the vacuum case <cit.>, allowing us to solve for ψ_4 as an infinite sum of separable solutions using Green's method.

From ψ_4, Teukolsky <cit.> showed us how to compute the rate of gravitational radiation at future null infinity

    d^2E_s/dtdΩ=lim_r→∞r^2/4 π ω̂^2|ψ_4|^2,

where ω̂=2 ω and dΩ is the induced volume on a unit 2-sphere.

We work with the scaled expression

    P_E=dE_s/dt(M/ℳ_s)^2 ,

where

    ℳ_s=∫_r_+^+∞∫_0^π∫_0^2π√(-g) T^t_t dϕ dθ dr

is the total scalar field energy, - the energy of the perturbed initial data. P_E is independent of the scalar field amplitude and measures the energy radiated per m mode in the initial data. See the Supplemental Material on how to compute P_E.

As a measure of curvature, we monitored the following time independent quantity as a function of m:

    χ≡max_r,θ (|ψ_4|^2/ℳ_s^2).

The radial and azimuthal location of the maximum of (<ref>), (r,θ)=(r_⋆,θ_⋆), track the maximum of [ψ^(1)]^2.

Our results for the GW emission are shown in Fig. <ref>. The radiated angular momenta in this process is P_J=m/Re(ω M)P_E, in accordance with <cit.>. The fact that both P_E and P_J appear to decrease rapidly with increasing m shows that the evolution occurs, to very good approximation, at fixed energy and angular momentum. This is akin to the time evolution of the superradiance instability with anti-de Sitter asymptotics <cit.>, simulated recently in <cit.>, and shows hints of turbulent behaviour.


The data in Fig. (<ref>) is for a fixed value of the dimensionless spin parameter a/M. However, during the aforementioned cascade the BH will be gradually spinning down, hence the gravitational radiation for each value of m should ideally be computed by accounting for BH's loss of energy due to the superradiant modes active prior to the one under consideration. Nevertheless, using the superradiant condition Re(ω)>m Ω_K, one sees that Δ(a/M) for successive superradiant modes ∼ℓ^-1 as ℓ→∞, implying that in the regime of interest the dimensionless spin will be approximately constant and Fig. (<ref>) represents accurately the qualitative behaviour.

In Fig. <ref> we show the dependence of χ on the initial data, here labelled by m. A WKB-type analysis reveals that the leading behaviour for χ is power law in 1/m

    χ_WKB = μ̂^10/m^16π^2 2^6(μ/r_+)^6 .

One can consistently correct this approximation to next-to-leading order (see the Supplemental Material).

Note that P_E is harder to compute numerically than χ, hence why we have extended results for χ up to m=5.




 Gedanken experiment – We are now ready to present our possible counterexample to WCCC with asymptotically flat boundary conditions. Consider generic initial data for the Einstein-Scalar system. These data are controlled by a large functional freedom coming from the fact that we can choose the initial metric on a constant time slice, as well as the extrinsic curvature (so long as the Hamiltonian and momentum constraints are satisfied). In addition, we can also control the initial profile for the scalar field and its first time derivative on a constant time slice. We are going to choose our initial data to be close to that of the Kerr BH, so that deviations from the Kerr metric only occur at order 𝒪(ψ^2). This condition can be relaxed by considering initial data for the purely gravitational sector that is small in some norm. This essentially means that all the dynamics are being generated by the scalar field.

For generic scalar field initial data, we expect the scalar field profile to have some support on the unstable modes of the preceding sections, i.e. to excite unstable modes. Since all other modes decay with time [This decay can be very complicated to determine and is not exponential with time.], we expect the late time evolution to be dominated by the leading unstable modes and their backreaction. For each value of m there is an infinite number of such modes labelled by ℓ≥|m|. However, all of these modes stop being unstable as soon as the condition Re(ω)>m Ω_K is no longer satisfied. The dynamics of this change in angular momentum and energy is entirely controlled by the ℓ=m modes, implying that, after some suitably long time, the dynamics of the Einstein-Scalar system can be well approximated by restricting our attention to scalar profiles of the form

    ψ(t,r,θ,ϕ)= Re[∑_ℓ=0^+∞ a_ℓe^-iω_ℓ+i ℓϕ R_ω_ℓℓℓ(r) S_ω_ℓℓℓ(θ)] ,

and determining its leading order backreaction on the spacetime curvature. The coefficients a_ℓ are determined by our choice of initial data: for finite Sobolev norm initial data we expect a_ℓ to exhibit polynomial behaviour in 1/ℓ, whereas for C^∞ initial data we expect the coefficients a_ℓ to decay faster than any polynomial in 1/ℓ. Note that for real analytic initial data one can show that a_ℓ≈ e^-αℓ, for α>0.

Given that each ℓ mode evolves on an exponentially different timescale, as shown by Eq. (<ref>), they effectively decouple from each other, allowing us to study each term in Eq. (<ref>) and its backreaction on the metric separately. Eventually, a given ℓ=ℓ_⋆ mode becomes stable, but the system remains unstable to higher values of ℓ>ℓ_⋆. This cascading happens slowly, since the time scales for this effect are exponentially large. One might worry that the energy contained in these high ℓ modes is radiated away as time passes by, but we have seen in Fig. (<ref>) that this is not the case. In fact, the larger the value of ℓ, the smaller its radiative power is.

Finally, we have seen in Fig. (<ref>) that χ decays as ℓ^-16. This in turn implies that a mode with weigh a_ℓ will depend on ℓ as  a_ℓ^4 ℓ^-16. The reason for this is simple: the Teukolsky scalar ψ_4 is sourced by [ψ^(1)]^2, and χ is related to |ψ_4|^2, which translates into the overall scaling mentioned above. The curvatures are thus suppressed with time, and the evolution continues until all of the angular momentum is deposited into the scalar cloud, and the central black hole becomes Schwarzschild. However, one can show <cit.> that massive scalar field perturbations decay extremely slow around Schwarschild black holes, with a dependence as weak as 1/log(log t), for large t. This suggests that the hypothetical end point is itself nonlinearly unstable through a mechanism similar to the one reported in <cit.>. 
On the other hand, this slow decay warrants the question whether the gravitational radiation might not just disperse the clouds, before any nonlinearities become problematic. We cannot answer this in our analysis, as our approximation for radiation emission breaks down at late times. Therefore the outcome of our thought experiment depends on which of the competing processes - the GW emission or the nonlinear issues resulting from slower than logarithmic decay - wins over. If it is the latter, then the lack of a possible stationary end point leads us to conjecture that the spirit, if not the letter, of weak cosmic censorship is violated. Whether the curvature will be infinite in finite time is a question that we cannot settle with our current methods.

Note also that r_⋆(ℓ) increases with ℓ, posing a problem from a numerical perspective: 1) the time scales involved are enormous; 2) the cascading towards high ℓ values makes this problem dependent on high frequency modes (as the simulation of turbulence in 3+1 nonrelativitic fluids); and 3) the integration domain must extend to spatial infinity to observe this effect.



 Conclusions – We have seen (Fig. <ref>) that the efficiency of superradiance cannot be counteracted by GW emission, implying that the system will continue advancing to higher values of m with curvatures decreasing appropriately (as shown by Fig. <ref>) until a configuration with a central Schwarzschild black hole is reached. However, Schwarzschild is likely to be nonlinearly unstable due to the very slow decay of perturbations induced by massive scalar fields. Reaching this troublesome regime, given that effects we cannot account for in our analysis do not prevent this, will involve time scales much longer than the age of our Universe, of course, as one will have to go to large values of ℓ=m . Nevertheless our scenario provides the first plausible example of a system with asymptotically flat boundary conditions, where WCCC is violated.




	 Acknowledgments – We would like to thank O. J. C. Dias and H. S. Reall for reading an earlier version of this manuscript. We also want to thank W. E. East for many useful discussions that led to this paper. F.C.E. and B.G. are supported by STFC. J.E.S is supported in part by STFC Grants PHY-1504541 and ST/P000681/1. This work used the DIRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.ukwww.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure Capital Grant ST/J005673/1, STFC Capital Grant ST/H008586/1, and STFC DiRAC Operations Grant ST/K00333X/1. DiRAC is part of the National e-Infrastructure.








§ SUPPLEMENTAL MATERIAL


 WKB expansion – We present the details of the expansion in large ℓ that we use to determine the scalar instability growth rate.  We work with l-m=n, where n∼𝒪(1), n∈ℕ_0. To this end we rewrite the Klein-Gordon equation from the main text in the following way

	
    r-r_+/r_+[Δ R_ωℓ m'   (r)]_,r+V(r) R_ωℓ m(r)=0,
    [sinθ S_ωℓ m'(θ)]_,θ/sinθ-   [ã^2 k̃ cos^2θ+m^2/sin^2θ-Λ]S_ωℓ m(θ)=0,
    
    	V(r)=a_1+   a_2 r/r_++a_3 r^2/r_+^2-k̃ r^3/r_+^3+a_4/r/r_+-ã^2,
    
    	a_1=Λ+   ã ω̃(1+ã^2)[ã ω̃(1+ã^2)-2 m],
    
    	a_2=ã^2 ω̃^2   (1+ã^2)-Λ,   a_3=μ̃^2+ã^2 ω̃^2,
    
    	a_4=   ã^2[m-ã ω̃(1+ã^2)]^2,


where k̃=μ̃^2-ω̃^2, Λ is the angular separation constant and we have introduced dimensionsless variables

    ã=a/r_+, M̃=M/r_+, ω̃=ω r_+, μ̃=μ r_+,

As explained in the main text, finding unstable modes amounts to determining the values of ω̃ for which ψ has ingoing boundary conditions at the event horizon (consistent with the equivalence principle), and finite energy on a partial Cauchy surface t= const [For unstable modes, this condition is equivalent to requiring the modes to be outgoing at spatial infinity.].

In the eikonal limit ℓ≫1 the spheroidal eigenvalue is Λ=ℓ(ℓ+1)+O(ℓ^-1); this can easily be verified by a series expansion of the angular equation after prefactoring sin^mθ from S(θ). In order to determine ω̃ we divide the domain in two intersecting regions, solve (<ref>) in each of them and then match the solutions in the overlap. To this end we introduce the new variable x=r/r_+-(1+ã^2)

    (x+ã^2)[Δ̃ R_ωℓ m'(x)]_,x   +Ṽ(x)R_ωℓ m(x)=0,
    
    V(x)=b_1+b_2 x+b_3    x^2-k̃ x^3+a_4/1+x,
    
    b_1=ω̃^2(1+ã^2)^2(1+4 ã^2)   -2 m ã ω̃(1+ã^2)
    
    -ã^2[ℓ(   ℓ+1)+(1+ã^2)^2μ̃^2],
    
    b_2=-(1+ã^2)[(1+3 ã^2)   μ̃^2-3 (1+2 ã^2)ω̃^2]
    
    -ℓ(ℓ+   1),
    
    b_3=(3+4 ã^2)ω̃^2-   (2+3 ã^2)μ̃^2,

with Δ̃=(1+x)(x+ã^2)

 Near-horizon region – Using quasimodes <cit.>, one can show that k̃=𝒪(ℓ^-2). We thus take the near-horizon region to be defined by x≪ℓ^2 and see that

    k̃ x^3≪ x m^2,   b_3 x^2≪ x m^2,

implying that we can drop the cubic and quadratic terms inside V(x), leaving us with

    (x+ã^2)[Δ̃ R_H'(x)]_,x+[b_1+b_2 x+a_4/1+x]R_H=0.

The solution [Multiply the equation by (1+x), then change variables to u=(x+ã^2)/(1-ã^2).] is a linear combination of Gauss hypergeoemtric functions

    R_H(x)=A_in(-u)^δ(1+u)^ϕ [_2]F_1[c_-,c_+;c;-u]
       +A_out(-1)^-2 δ(-u)^-δ(1+u)^ϕ [_2]F_1[λ_-,λ_+;2-c;-u],

where u=(x+ã^2)/(1-ã^2), A_in/out are constants and

    δ=   i Ω_0/1-ã^2[1+ã^4(1-ã^4)(2 μ̃^2-3 ω̃^2)/Ω_0^2]^1/2,
    ϕ=   -i m-ã ω̃ (1+ã^2)/1-ã^2,  Ω_0=m ã-ω̃ (1+ã^2)
    
    c_±=   1/2±√(1-4 b_2)/2+ϕ+δ,   c=1+2 δ,
    λ_±=   c_±-c+1.

Demanding ingoing waves only at the horizon (x=0) requires setting A_out=0. To see this, note that in the limit x→0 the hypergeometric functions take on a constant value to leading order, hence one just needs to know which of (-u)^±δ gives the correct behaviour there. This can most easily be deduced by performing a transformation from BL to Kerr coordinates

    dv=dt+r^2+a^2/Δdr,  dχ=dϕ+a/Δdr.

Looking at the near-horizon limit allows us to obtain

    e^-i ω te^i m ϕ=e^-i ω ve^i m χ(r-r_+)^i(ω̃1+ã^2/1-ã^2-m ã/1-ã^2),

which has opposite sign to what we expect at the horizon, as the transformation eliminates any singular behaviour.

 Far away region – Before zooming into spatial infinity we make the following transformation to (<ref>), to ensure that we will get the correct asymptotic behaviour,

    ℛ_ωℓ m(x)=   h(x) Q_∞(x), h(x)=√(4 k̃^2 x^2+4 k̃ x d_2+d_1^2)/2 k̃ x+d_1,
    
    d_1   =(1+2 ã^2)(2 μ̃^2-3 ω̃^2),
    
    d_2   =(1+3 ã^2)μ̃^2-(2+5 ã^2)ω̃^2.

We work with x≫1, which leaves us[First multiply (<ref>) by 1+x, then perform transformation (<ref>), followed by dividing the resulting equation by (1+x)(x+ã^2)(t_0+t_1 x+t_2 x^2+t_3 x^3+2 x^4)/[2 x (x+i_0)^2√(x(x+i_1)+i_0^2)], with t_0=i_0^3(1+ã^2)-2 ã^2 i_0^2+ã^2 i_0 i_1,   t_1=i_0(2 i_1(1+ã^2)+2 ã^2)-i_0^2(1+ã^2)+2 i_0^3-ã^2 i_1,   t_2=3 i_0(1+ã^2+i_1),   t_3=1+ã^2+4 i_0+i_1,   i_0=d_1/(2 k̃),   i_1=d_2/k̃, and expand near spatial infinity (x=1)] with

    x^2 Q”_∞(x)+2 x Q'_∞(x)   -(k̃ x^2+e_1 x-b_2)Q_∞(x)=0,
    
    e_1=(1+ã^2)   (μ̃^2-2 ω̃^2),

to solve. The above can be transformed into Whittaker's differential equation with the help of

    x=z/2 √(k̃), j(z)=Q_∞(z)/z,

resulting in

    j”(z)+(-1/4-e_1/2 √(k̃) z+b_2/z^2)j(z)=0.

The general solution to (<ref>) is then a linear combination of the Whittaker functions W_k,ν(z) and W_-k,ν(-z), thus

    R_∞=B_1 v(z) W_-κ,ν(-z)+B_2 v(z) W_κ,ν(z),
    κ=-e_1/2√(k̃), ν=1/2√(1-4 b_2), v(z)=h(z)/z,

with B_1/2 constants. Expanding for large z, and using

    z=2 √(k̃)[r/r_+-(1+ã^2)]2 √(k̃)r/r_+,

we see that 

    lim_r→∞[B_1 v(z) W_-κ,ν(-z)   +B_2 v(z) W_κ,ν(z)]=
    
    B_1 e^√(k̃) r/r_+r^-1-κ   +B_2 e^-√(k̃) r/r_+r^-1+κ.

Imposing that the solution decays at infinity sets B_1=0.

 Matching – We now perform the matching procedure in the overlapping region 1≪ x≪ℓ^2. To this end, we take x∼ℓ^3/2 and expand the near-horizon and far-away solutions for large and small variables respectively.

Looking at the near-horizon region first, we can straightforwardly use the standard asymptotic series of the Hypergeometric function, as all the parameters grow slower than the variable of the function, leaving us with

    lim_x→∞   R_H=A_in(-1)^δ Γ[c_+-c_-] Γ[c]/Γ[c_+] Γ[c-c_-][x/(1-ã^2)]^ν-1/2
       +A_in(-1)^δ Γ[c_–c_+] Γ[c]/Γ[c_-] Γ[c-c_+][x/(1-ã^2)]^-ν-1/2.


Next we consider the expansion of the far-away solution for small variable and again we can use the standard series in the literature to do so, as the parameters outgrow the argument of the Whittaker function, thus

    lim_x→ 0R_∞=-B_2   (2 √(k̃))^ν-1/2Γ[-2 ν]/Γ[1/2-ν-κ]x^ν-1/2
       -B_2(2 √(k̃))^-ν-1/2Γ[2 ν]/Γ[1/2+ν-κ]x^-ν-1/2.


Afterwards we equate the coefficients in front of the equivalent terms in (<ref>) and (<ref>) in order to derive

    Γ[1/2-ν-κ]/Γ[1/2+ν-κ]=(2 (1-ã^2)√(k̃))^2 νG(ℓ),
       G(ℓ)=Γ[c_–c_+] Γ[-2 ν] Γ[c_+] Γ[c-c_-]/Γ[2 ν] Γ[c_+-c_-] Γ[c_-] Γ[c-c_+].

Furthermore, in the limit l→∞ the RHS above has very small real and imaginary parts, implying that Γ[1/2+ν-κ] in the denominator on the LHS must have a pole, thus

    1/2+ν-κ=-N,

with N∈ℕ_0, corresponding to the radial node of the scalar field. This allows us to deduce the real part of √(k̃)

    Re(√(k̃))=-e_1/2N+2 ν+1.

Next, with |ω̃_R|≫|ω̃_I| <cit.>, implying |k̃_R|≫|k̃_I|, we can expand (<ref>) for large ℓ to obtain the real part of ω̃

    ω̃_R=μ̃-(1+ã^2)^2 μ̃^3/8 ℓ^2+𝒪(ℓ^-3),

confirming our expectations. Moreover, our calculation is accurate to 𝒪(1) only, hence we can replace ω̃ with μ̃ everywhere except inside k̃, where the leading order term cancels. This reduces (<ref>) to

    Γ[1/2-ν-κ]/Γ[1/2+ν-κ]=[μ̃^2(1-ã^4)/N+1+ν̂/2]^ν̂Ĝ(ℓ),
       Ĝ(ℓ)=Γ[-ν̂]^2 Γ[σ_++ϕ̂_+] Γ[σ_++ϕ̂_-]/Γ[ν̂]^2 Γ[σ_-+ϕ̂_+] Γ[σ_-+ϕ̂_-],
       ν̂=√((1+2 ℓ)^2-4(1+ã^2)(2+3 ã^2)μ̃^2),
       σ_±=1±ν̂/2+δ̂,  ϕ̂_±=±i ã μ̃ (m-ã(1+ã^2))/1-ã^2,
       δ̂=i ξ/1-ã^2[1-ã^4 μ̃^2 (1-ã^4)/ξ^2]^1/2,
       ξ=ã m-(1+ã^2)μ̃,

For the imaginary part of ω̃ we allow (<ref>) to be complex

    1/2+ν-κ=-N+ϵ,

with ϵ≪1, ϵ∈ℂ. We can look at the ϵ→0 limit of (<ref>), using (<ref>), to derive

    Im(√(k̃))=ϵ2 i e_1/(1+2 N+2 ν)^2,

which also enables us to find, in the limit ℓ→∞,

    ω̃_I=i Im(ϵ)(1+ã^2)^2 μ̃^3/4 ℓ^3.


Next, we look at (<ref>) for ϵ→0 and obtain

    (-1)^N N! ϵ=[μ̃^2(1-ã^4)/N+1+ν̂/2]^ν̂Ĝ(ℓ)/Γ[-N-ν̂].

Taking ℓ large, we can rearrange for ϵ, which allows us to derive an expression for ω̃_I via (<ref>). It is valid for any spin parameter | a|<M and scalar field mass μ,

	
    ω̃_I   =2^N-5-6 ℓ(1+ã^2)^3+2 ℓ((1-ã^2)^2ℓ^2+4 ã^2 m^2)^1/2+ℓμ̃^5+4 ℓ/√(π) N! ℓ^-N+11/2+6 ℓ sinh[2 π(ã m-(1+ã^2)μ̃)/1-ã^2]×
       exp[2 ℓ-2(N+1)+2 μ̃ (1+ã^2)^2arctan[2 ã m/ℓ(1-ã^2)]/1-ã^2-4 ã m arctan[2 ã m/ℓ(1-ã^2)]/1-ã^2].


Setting ℓ=m (n=0) we get the growth rate of the dominant modes in the spectrum

	
    ω̃_I,ℓ=m   =2^N-5-6 ℓ(1+ã^2)^4+4 ℓμ̃^5+4 ℓ/√(π) N! ℓ^-N+9/2+4 ℓsinh[2 π(ã ℓ-μ̃(1+ã^2))/1-ã^2]×
       exp[-2(N+1+μ̃ (1+ã^2)^2arctan[-2 ã/1-ã^2]/1-ã^2)+ℓ(2+4 ã arctan[-2 ã/1-ã^2]/1-ã^2)],


which is rewritten in the main text for compactness. Moreover, for ã=0 the sinh changes sign and we reproduce the correct behaviour for Schwarzschild <cit.>

    ω̃_I,ã=0=-2^N-5-6 ℓe^-2(N+1)+2 ℓμ̃^5+4 ℓ/√(π) N! ℓ^-N+9/2+4 ℓsinh[2 π μ̃].



 KG Numerically– We apply spectral collocation methods on a discretised Chebyshev grid with

    z=1+cosθ/2,   x=1-r_+/r,

where x=0 and x=1 correspond to the event horizon and spatial infinity, and z=0 and z=1 to the north and south poles, respectively. We use dimensionless variables (<ref>), factor out the singular behaviour at the boundaries

    R_ωℓ m(x)=(1-x)^β e^α/1-x x^γ ℛ_ωℓ m(x),
       S_ωℓ m(z)=z^m/2(1-z)^m/2 𝒮_ωℓ m(z),
       α=-√(μ̃^2-ω̃^2),  β=1+(1+ã^2)(μ̃^2-2ω̃^2)/2√(μ̃^2-ω̃^2),
       γ=-i(ω̃1+ã^2/1-ã^2-m ã/1-ã^2),

and finally apply Newton's method. We then solve for ℛ_ωℓ m(x) and 𝒮_ωℓ m(z).

No special treatment at the boundaries is needed as (<ref>) force the system to pick the right solution. An initial guess is (<ref>), Λ=ℓ(ℓ+1), ℛ_ωℓ m(x)=1, 𝒮_ωℓ m(z)=1.

 Computation of ψ_4 – We want to compute GWs emitted by a real scalar cloud around a Kerr BH. A complex scalar is not realistic and produces negligible gravitational radiation, since its stress-tensor is stationary - its use is a computational convenience for the determination of the field's QNMs, hence to find GWs we consider the real part of the superradiant modes - ψ_R= ψ.

As explained in the main text, we treat the cloud as a perturbing source and use Teukolsky's equation <cit.>

	
    [(Δ+3γ-γ̅+4μ+μ̅)(D+4ϵ-ρ)-3ψ_2-
       (δ̅+3α+β̅+4π-τ̅)(δ+4β-τ)]ψ_4=4π T_4,
       T_4=(Δ+3γ-γ̅+4μ+μ̅)×
       [(δ̅-2τ̅+2α)T_nm̅-(Δ+2γ-2γ̅+μ̅)T_m̅m̅]
       +(δ̅-τ̅+β̅+3α+4π)×
       [(Δ+2γ+2μ̅)T_nm̅-(δ̅-τ̅+2β̅+2α)T_nn].


to extract the NP scalar, ψ_4, encoding the information about the outgoing GWs. Our sign convention is (-,+,+,+), opposite to Teukolsky's. Nevertheless, the equation in terms of NP variables is unchanged. For their definitions, in our convention, we use <cit.>.

The tetrad projections we need in T_4 are T_nn, T_nm̅ and T_m̅m̅. Moreover using a null tetrad implies we can ignore the g_μν term. In a Kinnersley tetrad <cit.> we have

    T_nn=   e^-2 i ω te^2 i m ϕ/8 Σ^2[i K R_ωℓ m+Δ_k R_ωℓ m,r]^2S_ωℓ m^2,
    
    T_nm̅=   e^-2 i ω te^2 i m ϕ/4√(2) Σ(r-i a cosθ)[i K R_ωℓ m+Δ_k R_ωℓ m,r]×
       [S_ωℓ m,θ-(a ω sinθ-m/sinθ)S_ωℓ m]R_ωℓ m S_ωℓ m,
    
    T_m̅m̅=   e^-2 i ω te^2 i m ϕ/4(r-i a cosθ)^2 R_ωℓ m^2×
       [S_ωℓ m,θ-(a ω sinθ-m/sinθ)S_ωℓ m]^2.


Furthermore, in the main text it is stated that the GW frequency and mode number are related to the scalar field ones by ω̂=2 ω, m̂=2 m. To see this we first remind ourselves of our ansatz for the scalar field

    ψ(t,r,θ,ϕ)=e^-i ω t+i m ϕR_ωℓ m(r)S_ωℓ m(θ) .

Next, we define the stress tensor of a scalar as

    T_μν(Υ,Ψ)=2∇_μΥ∇_νΨ-g_μν(∇_σΥ∇^σΨ+μ^2ΥΨ),

and plug in ψ_R=1/2(ψ+ψ̅) for Υ and Ψ, leading us to

    T_μν(ψ_R,ψ_R)=1/4[T_μν(ψ,ψ)+T_μν(ψ̅,ψ̅)+2T_μν(ψ,ψ̅)],

with bar indicating complex conjugation. The ansatz (<ref>) implies that the terms above will contain exponential prefactors, which reveal that the first and second one source outgoing waves with m̂=2 m, ω̂=2 ω and m̂=-2 m, ω̂=-2 ω, respectively, whereas the third one corresponds to the energy of the cloud, given by the gravitational mode with m̂=ω̂=0.

The LHS of (<ref>) can be separated into angular and radial parts as in the vacuum case <cit.>. Thus take

    ψ_4=e^-i ω̂ te^i m̂ ϕ ρ^4 ℛ(r) 𝒮(θ),

with spin coefficient ρ=-1/(r-i a cosθ), and multiply both sides by 2 Σ(r,θ)/ρ^4. 𝒮(θ) can be identified with a spin-weighted spheroidal harmonic, satisfying

    1/sinθ[sinθ    [_s]𝒮_ℓ̂m̂,θ]_,θ+[(c cosθ)^2-2 c s cosθ+s
       +[_s]A_ℓ̂m̂ω̂-(m̂+s cosθ)^2/sinθ^2][_s]𝒮_ℓ̂m̂=0,

with c= a ω̂, [_s]A_ℓ̂m̂ω̂ the separation constant and s=-2. Multiplying by [_s]𝒮̅_ℓ̂m̂ and integrating over θ results in the following equation for the radial function

    Δ_k^2[Δ_k^-1ℛ_ℓ̂m̂,r]_,r+[K^2+4 i(r-M)K/Δ_k-8 i ω̂ r
    
    -a^2 ω̂^2+2 a m̂ ω̂-[_s]A_ℓ̂m̂ω̂]ℛ_ℓ̂m̂=𝒯_ℓ̂m̂ω̂,

where K=(r^2+a^2)ω̂-a m̂ and the source term is the integrated over angles stress-energy tensor

    𝒯_ℓ̂m̂ω̂=4π/η_ℓ̂m̂∫2 Σ(r,θ)/ρ^4[_s]𝒮̅_ℓ̂m̂ T_4 sinθ dθ,

with the normalisation condition

    ∫_0^π|[_s]𝒮_ℓ̂m̂(θ)|^2sinθ dθ=η_ℓ̂m̂.

Equation (<ref>) should be evaluated for ω̂=±2 ω, m̂=±2 m, s=-2 on the LHS with T_μν(ψ,ψ), for the +, and T_μν(ψ̅,ψ̅), for the -, inside T_4 on the RHS, and for m=ω̂=0 on the LHS with T_μν(ψ,ψ̅) in T_4 on the RHS.

However, the Teukolsky equation is invariant under complex conjugation followed by m→-m and ω→-ω, hence for m̂0ω̂ we only need to determine the contribution of T_μν(ψ,ψ) and double the result.

Moreover, by looking at the asymptotic behaviour of the homogeneous radial equation ((<ref>) with 𝒯_ℓ̂m̂ω̂=0) we can show that the m̂=ω̂=0 mode is subleading to all the rest at spatial infinity. To this end, transform the homogeneous (<ref>) with the help of

    ℛ_ℓ̂m̂(r)=Δ_K^-s/2/√(r^2+a^2)Y_ℓ̂m̂(r),  dr_*/dr=r^2+a^2/Δ_K,

and then take the limit of r→∞, leaving us with

    Y_ℓ̂m̂”(r_*)+[ω̂^2+2 i s ω̂/r]Y_ℓ̂m̂(r_*)=0,

for finite[Due to the aforementioned symmetry of the Teukolsky equation, we can work, without loss of generality, with the case of ω̂=2 ω, m̂=2 m only.] frequency ω̂ and azimuthal number m̂. The outgoing part of Y_ℓ̂m̂ behaves as e^i r_* ω̂r^-s, leading to an e^i r_* ω̂r^-2s-1 asymptotic behaviour for ℛ_ℓ̂m̂, which together with the definition (<ref>) allows us to deduce that the outgoing contribution of ψ_4 near spatial infinity, for finite ω̂ and m̂, behaves as

    lim_r→∞ψ_4∼ e^i ω̂ r_*/r  ( outgoing mode).

However, were we to set ω̂=m̂=0 before we take the limit of r→∞, instead of (<ref>), we are left with

    Y_00”(r_*)-ℓ(ℓ+1)/r^2Y_00(r_*)=0,

where we have used [_s]A_ℓ̂00=ℓ(ℓ+1)-s(s+1). The solution with finite energy at infinity corresponds to asymptotic behaviour of the form r^-ℓ-3 for ψ_4, which is clearly subleading to non-zero ω̂ and m̂ modes.

Therefore, for the computation of the gravitational radiation from the scalar clouds, which is performed at spatial infinity, we only need to look at the ω̂=2 ω, m̂=2 m case for (<ref>), with T_μν(ψ,ψ) inside T_4 on the RHS.

In this way, ψ_4 has been projected onto a basis of spin-weighted spheroidal harmonics

    ψ_4=e^-i ω̂ te^i m̂ ϕ ρ^4 ∑_ℓ̂ℛ_ℓ̂m̂(r) [_s]𝒮_ℓ̂m̂(θ),

where the definition of [_s]𝒮_ℓ̂m̂(θ) requires ℓ̂≥m̂.

Thus, given a Kerr BH, for fixed M μ and J/M^2, and its ℓ=m scalar superradiant spectrum, we can use the latter as a source on the RHS of (<ref>), which sets the values of the GW frequency ω̂ and mode-number m̂, letting us solve (<ref>) for any allowed value of ℓ̂.

Knowing ψ_4, we can compute the gravitational radiation at infinity and the maximum of the curvature component represented by the NP scalar ψ_4 with the formulae given in the main text. The argument for ignoring the m̂=ω̂=0 mode in the computation of the gravitational radiation at infinity does not hold for the maximum of ψ_4, as the latter is expected to be located at a finite distance from the black hole, near the peak of the scalar cloud, as we find in our numerics. Nevertheless, we have checked numerically that its contribution to the maximum is negligible compared to the rest.

 Numerical integration – We use the same numerical method for (<ref>) and (<ref>) as for (<ref>) and (<ref>).

 The angular equation – (<ref>) is an eigenvalue problem which we integrate using Newton's method. We write

    [_s]𝒮_ℓ m=z^ι_-(1-z)^ι_+[_s]𝕊_ℓ m,

where ι_±=|m± s/2| and solve for [_s]𝕊_ℓ m. We drop the hats on m and ω, as we solve the equation generally. The eigenvalue [_s]A_ℓ mω is treated as an unknown and we normalise [_s]𝕊_ℓ m to 1 using (<ref>). Appropriate boundary conditions are imposed at the edges of the domain. These can be obtained by expanding the equation for [_s]𝕊_ℓ m in power series around the poles. At z=0 we require

    [-4+(m-s)^2   -8 ι_–(2 ι_-)^2].∂_z[_s]𝕊_ℓ m(z)|_z=0=
    
    -[m^2-4 s+2    m s-4 [_s]A_ℓ m-8 ã s ω-4 ã^2ω^2+
    
    (2 ι_-)^2+4 ι_+   +4 ι_-(1+2 ι_+)].[_s]𝕊_ℓ m(z)|_z=1,

and similarly at z=1

    [-4+(m+s)^2   -8 ι_+-(2 ι_+)^2].∂_z[_s]𝕊_ℓ m(z)|_z=1=
    [m^2-4 s-2    m s-4 [_s]A_ℓ m+8 ã s ω-4 ã^2ω^2+
    
    (2 ι_+)^2+4 ι_+   +4 ι_-(1+2 ι_+)].[_s]𝕊_ℓ m(z)|_z=1.

A starting point to our iteration is the series solution, provided in <cit.>, up to twelfth order.

 The radial equation – (<ref>) is a sourced ODE, which we can invert once the RHS is known. As explained above, we need the superradiant modes of the scalar field on a fixed Kerr background, which we have numerically, so we can integrate (<ref>) for any ℓ̂. We take out the prefactors in (<ref>) from the scalar field inside T_4, so that we can directly substitute our numerical solutions ℛ_ωℓ m(x) and 𝒮_ωℓ m(z). We then factor out the singular behaviour of ℛ_ℓ̂m̂ at infinity and the horizon, derived via Frobenius analysis, and transfer them to the RHS. Specifically

    ℛ_ℓ̂m̂(x)=r_+^2(1-x)^β̂ e^α̂/1-x x^2+2 γ̂ ℝ_ℓ̂m̂(x),

with ℝ_ℓ̂m̂(x) dimensionless and

    α̂=2 i ω̃, β̂=-3-2 i(1+ã^2) ω̃,
    γ̂=-i(ω̃1+ã^2/1-ã^2-m̂ ã/1-ã^2).

Moreover, x^2 γ̂ gets cancelled by the x^γ from the radial scalar functions (<ref>) inside T_4. We then solve for ℝ_ℓ̂m̂(x).
We need to impose boundary conditions at spatial infinity, since <cit.> reveals that the ingoing behaviour of the radial function (the one we do not want) is subleading

    ℛ_ℓ̂m̂(x)∼ Z^(in)_ℓ̂m̂ e^-2 i ω̃/1-x(1-x)+Z^(out)_ℓ̂m̂e^2 i ω̃/1-x/(1-x)^3,

The boundary condition that selects the outgoing waves only can be deduced by expanding the homogeneous Teukolsky equation for large radial variable , which shows that setting Z^(in)_ℓ̂m̂=0 is equivalent to demanding

    ∂_xx.ℝ_ℓ̂m̂(x)|_x=1+1/4(1-ã^2)^2ϖ̂^2[(1-ã^2)^2[_s]A_ℓ̂m̂ϖ̂^2
       +2(1-ã^2)[_s]A_ℓ̂m̂ϖ̂[]2 ã m̂ ϖ̂+ã^4ϖ̂^2+2 ã^6ϖ̂^2
       -(i+ϖ̂)^2+ã^2(-1+4 i ϖ̂-3 ϖ̂^2))
       +2 ϖ̂(-3 i-12 ϖ̂+20 i ϖ̂^2+ã^8(8 i-11/2ϖ̂)ϖ̂^2
       +8 ϖ̂^3+2 ã^12 ϖ̂^3+2 ã^10ϖ̂^2(i+ϖ̂-2 ã^7 m̂ ϖ̂(i-2 ϖ̂)
       +4 ã^5 m̂(1+i/2ϖ̂+1/2ϖ̂^2)-8 ã^3 m̂(1-3 i/2ϖ̂+3/4ϖ̂^2)
       +4 ã m̂(1-3 i ϖ̂-2 ϖ̂^2)+ã^4(3 i-2 ϖ̂-20 i ϖ̂^2+1/2ϖ̂^3)
       +ã^2(3 i+19 ϖ̂+2 m̂^2 ϖ̂-4 i ϖ̂^2+12 ϖ̂^3)
       -ã^6(3 i+5 ϖ̂+6 i ϖ̂^2+11 ϖ̂^3)]]ℝ_ℓ̂m̂(x=1)=0,

where ϖ̂=2 ω̃.

In terms of the functions for which we solve numerically, the formula for the gravitational radiation given in the main text takes the following form

    dE_s/dt=1/(2 ω̃)^2∑_ℓ̂η_ℓ̂m̂|ℝ̂_ℓ̂m̂(x)|^2_x→1,

where we have used (<ref>), the orthogonality of the spin-weighted spheroidal harmonics, and have multiplied by 2 to accound for the GWs with m̂=-2 m, ω̂=-2 ω.

 Numerical convergence – As we cannot integrate (<ref>) for infinitely many ℓ̂, we truncuate the sum until (<ref>) converges. We will look at the m=4 case here (m̂=8 for the GWs), as this was the hardest one to tackle numerically. We include 9 GW modes, ℓ̂=8 to ℓ̂=16, to get good convergence for the radiated energy. Moreover, very high grid resolution was needed in the radial direction, in order to resolve the oscillating behaviour of the solution far away from the BH (exactly at spatial infinity the oscillating part is discarded by the boundary conditions). This is summarised in Figs. (<ref>) and (<ref>)



 Analytic approximation to χ – In this section we outline the steps for the WKB approximation of χ≡max_r,θ (|ψ_4|^2/M_s^2). We work with ℓ=m and ℓ≫1. 

First, we need a WKB expansion of the scalar field near its maximum - i.e. R_ωℓ m(r), S_ωℓ m(θ), Λ and ω̃, as well as the location of the minimum of the potential of the radial KG equation expanded for large ℓ=m.

We start with the angular equation (<ref>), which rewriten in our numerical coordinates (<ref>), together with the redefinition (<ref>), looks like

    (1-z) z 𝒮_ωℓ m”(z)+(m+1)(2z-1)𝒮_ωℓ m'(z)
       -[Λ-m(m+1)-ã^2 k̃ (2z-1)^2]𝒮_ωℓ m(z)=0.

Expand the variables in series in m

    𝒮_ωℓ m(z)   =1+S_1(z)/m+S_2(z)/m^2+𝒪(m^-3),
    Λ   =λ_1 m^2+λ_2 m+λ_3+𝒪(m^-1),
    ω̃   =μ̃-(1+ã^2)^2 μ̃^3/8 m^2+ω_3/m^3+𝒪(m^-4),

plug them in and solve order by order, requiring regularity of the solution. The known behaviour of ω̃ from before is utilised. In this way one can obtain S_i(z) and λ_i to any order as a function of the ω_i. 

Afterwards we transform the radial part (<ref>) using 

    r=2 r_+ m(m+1)/(1+ã^2)μ̃^2u.

The equation will not be given here - the above transformation correponds to the first two terms in the expansion for large ℓ=m of the position of the minimum of the radial potential. This is straightforward to get after deriving the leading order behaviour of the angular eigenvalue Λ=m(m+1)+𝒪(1), for which no knowledge of ω_3 or higher is required. Afterwards, we take

    R_ωℓ m(u)=e^m f(u)[1+R_1(u)/m+𝒪(m^-2)],

and solve order by order for f(u), R_i(u) and ω_i, again demanding regularity of the solution. This can be done to any order, but we only give the first few here

    S_1(z)   =c_2,   S_2(z)=c_3,
    
    f(u)   =c_0-u+log u,
    
    R_1(u)   =c_1-μ̃^2/4(1+ã^2)^2[u+1/u+3 log u],
    λ_1   =λ_2=1,  λ_3=0,  ω_3=(1+ã^2)^2 μ̃^3/4,

whereby c_i are constants that can for example be chosen so that the scalar field is unity at the maximum (u=1), which is what we do to compare with the numerics (where we can just divide by the numerical maximum).

Next, we move onto the Teukolsky equation (<ref>). The expressions are very lengthy and not illuminating at all, so we will give only the most important ones. Treating the angular equation (<ref>) as the scalar one produces

    [_s]𝒮_ℓ̂m̂   =z^|2m-s/2|(1-z)^|2m+s/2|[1-
       -2 ã μ̃ z(s+ã μ̃ (1-z))+d_1/m+𝒪(m^-2)],
    [_s]A_ℓ̂m̂ω̂   =2 m (2 m+1)-s(s+1)+𝒪(m^-1),

whereby we have accounted for m̂=2 m and d_1 is a constant that can be chosen to make the normalisation of the function easier[Depending on how one does the numerics, [_s]𝒮̅_ℓ̂m̂ inside T_4 might also need to be normalised.]. Afterwards we need to obtain a series expansion in ℓ=m for 𝒯_ℓ̂m̂ω̂ (<ref>). The normalisation η_ℓ̂m̂ has to be obtained by plugging (<ref>) in (<ref>) and integrating. The final result is

    𝒯_ℓ̂m̂ω̂=e^2 m(1-u+log u)16 √(2) π r_+^2 u^4/(1+ã^2)^4 μ̃^6m^9[1-4/m
       -μ̃^2 [3 (1+ã^4)(u-1)^2+ã^2(u (6 u-13)+6)]/6 u m
    
    -   2 i μ̃(1+ã^2)+3 μ̃^2(1+ã^2)^2log (u)/2 m+𝒪[1/m^2]].

With all ingredients present, we change variables in (<ref>) with (<ref>) and WKB expand ℛ_ℓ̂m̂ for large ℓ=m. Matching order by order resutls in

    ℛ_ℓ̂m̂=-e^2 m(1-u+log u)√(2) π r_+^2 u^2/(1+ã^2)^2 μ̃^4m^5[1-2/m
       -μ̃^2 [3 (1+ã^4)(u-1)^2+ã^2(u (6 u-13)+6)]/6 u m
    
    -   2 i μ̃(1+ã^2)+3 μ̃^2(1+ã^2)^2log (u)/2 m+𝒪[1/m^2]],

allowing us to to rebuild ψ_4 for ℓ=m via (<ref>).

We also need an expression for the scalar cloud energy

    ℳ_s=∫_r_+^+∞∫_0^π∫_0^2π√(-g) T^t_t dϕ dθ dr.

Working with coordinates z and u, we substitute our expansions for the scalar field in the definition of ℳ_s and after simplifications the angular integrals can be done analytically - with Mathematica or with a book on integrals. The resulting expression is expanded in ℓ=m and integrated over u from μ̃^2(1+ã^2)/2 m (m+1) to infinity, producing

    ℳ_s=e^2 mπ^3/2 Γ[3+2m,(1+ã^2)μ̃^2/m+1]/2^4 m-1m^2 m-5/2(1+ã^2)^3 μ̃^4[1+𝒪(m^-1)],

where we also have higher order terms, but the expressiosn are too cumbersome.

Having ψ_4 and ℳ_s in the large ℓ=m limit allows one to derive the expression for χ_WKB in the main text by dividing them and expanding in series in ℓ=m.

unsrt
