Plausible scenario for a generic violation of the weak cosmic censorship conjecture in asymptotically flat four dimensions

TL;DR

This paper proposes a plausible counterexample to the weak cosmic censorship conjecture in four-dimensional Einstein-Scalar theory, using analysis of the Klein-Gordon equation on Kerr black holes and demonstrating potential nonlinear instability.

Contribution

It introduces a new scenario where the weak cosmic censorship conjecture may be violated, supported by numerical and analytical analysis of quasinormal modes and backreaction effects.

Findings

Numerical and analytical quasinormal mode spectra agree.

Evidence of nonlinear instability at late times.

Potential violation of weak cosmic censorship in this setting.

Abstract

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.

Equations

$$S={\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left(\frac{R}{16\pi G}-{\nabla }_a\psi {\nabla }^a\psi -{\mu }^2{\psi }^2\right),$$

$$R_{ab}-\frac{R}{2}{g}_{ab}=8\pi G{T}_{ab}$$

$$\square \psi ={\mu }^2\psi ,$$

$$T_{ab}=2{\nabla }_a\psi {\nabla }_b\psi -g_{ab}{\nabla }^c\psi {\nabla }_c\psi -{\mu }^2{\psi }^2{g}_{ab}.$$

$$\begin{gathered} \mathrm{d}{s}^2=-\frac{\Delta }{{\Sigma }^2}{\left(\mathrm{d}t-a{\sin }^2\theta \mathrm{d}\varphi \right)}^2+{\Sigma }^2\left(\mathrm{d}{\theta }^2+\frac{\mathrm{d}{r}^2}{\Delta }\right) \\ +\frac{{\sin }^2\theta }{{\Sigma }^2}[a\mathrm{d}t-(r^2+a^2)\mathrm{d}\varphi {]}^2, \end{gathered}$$

$$\psi (t,r,\theta ,\varphi )=e^{-i\omega t+im\varphi }{\hat{\psi }}_{\omega m}(r,\theta )$$

$$\Delta\,\big{[}\Delta\,{R_{\omega\ell m}}_{,r}\big{]}_{,r}+V(r)R_{\omega\ell m}=0,$$

$$\frac{1}{\sin\theta}\big{[}\sin\theta\,{S_{\omega\ell m}}_{,\theta}\big{]}_{,\theta}\\ -\left[a^{2}\,k^{2}\cos^{2}\theta+\frac{m^{2}}{\sin^{2}\theta}-\Lambda\right]S_{\omega\ell m}=0,$$

$$V(r)$$

$$(2M\Lambda -4amM\omega +2M{a}^2{\omega }^2)r-a^2(\Lambda -m^2),$$

$$\mathrm{Re}(\omega M)=\hat{\mu }\left(1-\frac{{\hat{\mu }}^2}{2{\ell }^2}\right)+\mathcal{O}({\ell }^{-3}),$$

$$\begin{gathered} \mathrm{Im}(\omega M)=\frac{{\ell }^{-4\ell -\frac{9}{2}+p}}{2^{2\ell +1-p}\sqrt{\pi }p!}{\hat{\mu }}^{4\ell +5}\sinh \left[\frac{\pi \left(\ell {\Omega }_K-\mu \right)}{{\kappa }_K}\right]\times \\ \exp \left[-\frac{2}{{\kappa }_K}\left(\ell {\Omega }_K-\frac{\hat{\mu }}{r_{+}}\right)\arctan \left(\frac{{\Omega }_K}{{\kappa }_K}\right)-2(1-\ell +p)\right]\times \\ \left[1+\mathcal{O}({\ell }^{-1})\right], \end{gathered}$$

$$R_{\omega \ell m}(r)={\left(1-\frac{r_{+}}{r}\right)}^{-i\frac{\omega -m{\Omega }_K}{2{\kappa }_K}}{e}^{-kr}{\left(\frac{r_{+}}{r}\right)}^{\gamma }{q}_r(r),$$

$$S_{\omega \ell m}(\theta )={\sin }^m\theta q_{\theta }(\theta ),$$

$$\psi =\sum _{i=0}^{+\infty }{\psi }^{(2i+1)}{\epsilon }^{2i+1},\text{and}g=g_K+\sum _{i=1}^{+\infty }{g}^{(2i)}{\epsilon }^{2i},$$

$$\displaystyle\Big{[}\big{(}\Delta+3\gamma-\bar{\gamma}+4\mu+\bar{\mu}\big{)}\big{(}D+$$

$$\displaystyle-\big{(}\bar{\delta}+3\alpha+\bar{\beta}+4\pi-\bar{\tau}\big{)}\big{(}\delta+$$

$$\frac{{\mathrm{d}}^2{E}_s}{\mathrm{d}t\mathrm{d}\Omega }=\underset{r\to \infty }{\lim }\frac{r^2}{4\pi {\hat{\omega }}^2}{\left|{\psi }_4\right|}^2,$$

$$P_E=\frac{\mathrm{d}{E}_s}{\mathrm{d}t}{\left(\frac{M}{{\mathcal{M}}_s}\right)}^2,$$

$${\mathcal{M}}_s={\int }_{r_{+}}^{+\infty }{\int }_0^{\pi }{\int }_0^{2\pi }\sqrt{-g}\text{tensor}T{}_t^t\mathrm{d}\varphi \mathrm{d}\theta \mathrm{d}r$$

$$\chi \equiv \underset{r,\theta }{\max }\left(|{\psi }_4{|}^2/{\mathcal{M}}_s^2\right).$$

$${\chi }_{\mathrm{WKB}}=\frac{{\hat{\mu }}^{10}}{m^{16}{\pi }^2{2}^6}{\left(\frac{\mu }{r_{+}}\right)}^6.$$

$$\psi (t,r,\theta ,\varphi )=\mathrm{Re}\left[\sum _{\ell =0}^{+\infty }{a}_{\ell }{e}^{-i{\omega }_{\ell }+i\ell \varphi }{R}_{{\omega }_{\ell }\ell \ell }(r)S_{{\omega }_{\ell }\ell \ell }(\theta )\right],$$

$$\displaystyle\frac{r-r_{+}}{r_{+}}\big{[}\Delta\,R_{\omega\ell m}^{\prime}$$

$$\displaystyle\frac{\big{[}\sin\theta\,S_{\omega\ell m}^{\prime}(\theta)\big{]}_{,\theta}}{\sin\theta}-$$

$$V(r)=a_1+$$

$$a_1=\Lambda +$$

$$a_2={\tilde{a}}^2{\tilde{\omega }}^2$$

$$a_4=$$

$$\tilde{a}=\frac{a}{r_{+}},\tilde{M}=\frac{M}{r_{+}},\tilde{\omega }=\omega r_{+},\tilde{\mu }=\mu r_{+},$$

$$\displaystyle(x+\tilde{a}^{2})\big{[}\tilde{\Delta}\,R_{\omega\ell m}^{\prime}(x)\big{]}_{,x}$$

$$V(x)=b_1+b_{2}x+b_3$$

$$b_1={\tilde{\omega }}^2{\left(1+{\tilde{a}}^2\right)}^2\left(1+4{\tilde{a}}^2\right)$$

$$\displaystyle-\tilde{a}^{2}\big{[}\ell($$

$$\displaystyle b_{2}=-\left(1+\tilde{a}^{2}\right)\big{[}\left(1+3\,\tilde{a}^{2}\right)$$

$$-\ell (\ell +$$

$$b_3=\left(3+4{\tilde{a}}^2\right){\tilde{\omega }}^2-$$

$$\tilde{k}{x}^3\ll x{m}^2,b_3{x}^2\ll x{m}^2,$$

$$(x+\tilde{a}^{2})\big{[}\tilde{\Delta}\,R_{H}^{\prime}(x)\big{]}_{,x}+\Big{[}b_{1}+b_{2}\,x+\frac{a_{4}}{1+x}\Big{]}R_{H}=0.$$

$$\displaystyle R_{H}(x)=A_{in}(-u)^{\delta}(1+u)^{\phi}\,\tensor{[}_{2}]{F}{{}_{1}}\Big{[}c_{-},c_{+};c;-u\Big{]}$$

Full text

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 Senovilla:2018aav . 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 Penrose:1969pc . 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 Bardeen:1972fi ; Starobinsky:1973aij ; Damour:1976kh ; Zouros:1979iw ; Detweiler:1980uk ; Dolan:2007mj ; Yoshino:2013ofa ; Brito:2014wla ; Brito:2015oca ; Shlapentokh-Rothman:2013ysa , afflicting massive scalar perturbations around Kerr black holes (BHs) Kerr:1963ud - rotating, with spherical topology and considered the most general BH solutions of the vacuum Einstein equation Robinson:2004zz .

The angular dependence of such perturbations around Kerr is parametrized by two integers $\{m,\ell\}$: $m$ counts the number of nodes in the azimuthal direction and $\ell-|m|$ the number of zeros in the polar direction. We show that for any scalar field of mass $\mu$ and any nonzero value of the BH spin, for sufficiently large values of $\ell=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\,\ell\,\log\ell}$, indicating that each of the $\ell$ modes decouples from the rest, evolving independently.

As time progresses, modes with smaller values of $\ell$ stabilize one by one, forming scalar clouds around the BH, similar to those in Chodosh:2015oma ; Herdeiro:2014goa . However, these BHs were shown to be unstable to higher $m$-modes Ganchev:2017uuo , giving rise to the expectation of a cascade towards larger values of $\ell$. This corresponds to a transfer of energy from lower $\ell$-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 Sasaki:2003xr ; Brito:2017zvb . 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\,\mu$, and spin parameter, $J/M^{2}$, as a function of $\ell=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 $\ell$. 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 $\ell$ radiate less, implying that energy is accumulated at small scales more efficiently for larger values of $\ell$. Finally we end with discussion of the results.

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

[TABLE]

where $\mu$ is the scalar field mass, $g_{ab}$ the spacetime metric and $R$ its Ricci scalar. The equations of motion are

[TABLE]

An important solution to these equations is the Kerr BH Kerr:1963ud , with $\psi=0$ and

[TABLE]

where $\Delta=r^{2}+a^{2}-2\,M\,r$, $\Sigma^{2}=r^{2}+a^{2}\cos^{2}\theta$, $\phi\in(0,2\pi)$ is a periodic coordinate and $\theta\in(0,\pi)$ is a polar coordinate. The BH event horizon is a null hypersurface with $r=r_{+}\equiv M+\sqrt{M^{2}-a^{2}}$, angular velocity $\Omega_{K}=a/(a^{2}+r_{+}^{2})$ and surface gravity $\kappa_{K}=(r_{+}^{2}-a^{2})/[2r_{+}(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|\leq M$ with the inequality being saturated at extremality, when the Kerr BH event horizon becomes degenerate with $\kappa_{K}=0$.

We study Eq. (2b) on a fixed Kerr background (3), which is stationary and axisymmetric with respect to the Killing fields $\partial/\partial t$ and $\partial/\partial\phi$, respectively. We consider perturbations of the form

[TABLE]

and assume that $\widehat{\psi}_{\omega m}(r,\theta)$ is separable, $\widehat{\psi}_{\omega m}(r,\theta)=R_{\omega\ell m}(r)S_{\omega\ell m}(\theta)$, with the label $\ell$ anticipating that the separation constant will be parametrized by an integer $\ell$. Not all solutions to Eq. (2b) are separable, but we are interested in those composed of the sum (possibly infinite) of such separable solutions. $\omega\in\mathbb{C}$ is a complex frequency, determined by imposing appropriate boundary conditions. We are interested in finding unstable mode solutions for which $\mathrm{Im}(\omega)>0$.

Inserting the ansatz (4) into Eq. (2b) yields a system of two second order ordinary differential equations, coupled via the separation constant $\Lambda$,

[TABLE]

with $k\equiv\sqrt{\mu^{2}-\omega^{2}}$. Finding unstable modes amounts to finding the values of $\omega$ for which $\psi$ has ingoing boundary conditions at the event horizon (consistent with the equivalence principle) and finite energy on a partial Cauchy surface $t={\rm 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 $\mu$, it reproduces the results in Detweiler:1980uk 111The large $m$ limit does not commute with the extremal limit. It can be shown that at any finite value of $m$, $\mathrm{Im}(\omega)=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\to+\infty$ limit. The behaviour in this region could be investigated by switching the order of limits. (up to a factor of $2$, see Furuhashi:2004jk ). For large $m$, Eq. (5b) can be approximated by the usual equation for spherical harmonics on a 2-sphere so that $\Lambda=\ell(\ell+1)+O(\ell^{-1})$ and $\ell,\,m$ are integers with $\ell\geq 0$, $|m|\leq\ell$. In the following we set $m=\ell$ and take the limit $\ell\gg 1$.

To determine the large $\ell$ limit at fixed $\mu$ and $a$, we used the same method as in Eperon:2016cdd , which combines a matched asymptotic and WKB type approach. We quote the final result and leave the details to the Supplemental Material,

[TABLE]

where $\hat{\mu}\equiv\mu M$ and $p\in\mathbb{N}_{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 $\ell\Omega_{K}=\mathrm{Re}(\omega)$. More importantly for our purposes, in the limit $\ell\to+\infty$, the growth rate of the instability scales as $e^{-4\ell\log\ell}$, and no matter what the value of $\mu$ or $a$, one can always find a value of $\ell=\ell_{\star}\equiv\lceil\mu/\Omega_{K}\rceil$ 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 $\mu$ 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

[TABLE]

with $q_{r}(r_{+})=q_{\theta}(0)=1$ and $\gamma\equiv 1+2Mk-M\mu^{2}/k$. We now use the methods of Dias:2015nua to solve for the eigenpair $(\omega,\Lambda)$ using a Newton-Raphson routine 222The 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. 1, our numerical data agrees excellently with (6b). Furthermore, one can measure deviations of our WKB expression to the exact numerical result, and it agrees with the error given in Eq. (6b).

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 East:2017ovw ; East:2017mrj 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

[TABLE]

where $g_{K}$ is given by the Kerr metric (3). We expand the equations of motion (2c) in a power series in $\varepsilon$. To first order in $\varepsilon$ we solve eq. (2b) subject to a choice of initial data. For the case at hand, we choose $\psi$ 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 $\mathrm{Im}(\omega M)\ll\mathrm{Re}(\omega M)$, we take $\omega$ 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 $\psi_{4}$ since it also allows us to efficiently compute the GWs emitted by the scalar cloud. $\psi_{4}$ obeys the Teukolsky equation Teukolsky:1973ha ; Press:1973zz ; Teukolsky:1974yv

[TABLE]

whereby all quantities appearing in Eq. (9) are given in the Supplemental Material. The source term $T_{4}$ is reconstructed directly from Eq. (2c) and thus depends on $\psi^{(1)}$ and its gradient only.

The lhs of (9) can be separated into angular and radial parts as in the vacuum case Teukolsky:1973ha , allowing us to solve for $\psi_{4}$ as an infinite sum of separable solutions using Green’s method.

From $\psi_{4}$, Teukolsky Teukolsky:1974yv showed us how to compute the rate of gravitational radiation at future null infinity

[TABLE]

where $\hat{\omega}=2\,\omega$ and $\mathrm{d}\Omega$ is the induced volume on a unit 2-sphere.

We work with the scaled expression

[TABLE]

where

[TABLE]

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$:

[TABLE]

The radial and azimuthal location of the maximum of (13), $(r,\theta)=(r_{\star},\theta_{\star})$, track the maximum of $[\psi^{(1)}]^{2}$.

Our results for the GW emission are shown in Fig. 2. The radiated angular momenta in this process is $P_{J}=\frac{m}{\mathrm{Re}(\omega M)}P_{E}$, in accordance with Teukolsky:1974yv . 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 Dias:2011at ; Niehoff:2015oga , simulated recently in Chesler:2018txn , and shows hints of turbulent behaviour.

The data in Fig. (2) 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 $\mathrm{Re}(\omega)>m\Omega_{K}$, one sees that $\Delta(a/M)$ for successive superradiant modes $\sim\ell^{-1}$ as $\ell\rightarrow\infty$, implying that in the regime of interest the dimensionless spin will be approximately constant and Fig. (2) represents accurately the qualitative behaviour.

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

[TABLE]

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 $\chi$, hence why we have extended results for $\chi$ 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 $\mathcal{O}(\psi^{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 333This 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 $\ell\geq|m|$. However, all of these modes stop being unstable as soon as the condition $\mathrm{Re}(\omega)>m\Omega_{K}$ is no longer satisfied. The dynamics of this change in angular momentum and energy is entirely controlled by the $\ell=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

[TABLE]

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

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

Finally, we have seen in Fig. (3) that $\chi$ decays as $\ell^{-16}$. This in turn implies that a mode with weigh $a_{\ell}$ will depend on $\ell$ as $a_{\ell}^{4}\,\ell^{-16}$. The reason for this is simple: the Teukolsky scalar $\psi_{4}$ is sourced by $[\psi^{(1)}]^{2}$, and $\chi$ is related to $|\psi_{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 futureWork 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 Eperon:2016cdd ; Keir:2014oka ; Keir:2016azt . 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_{\star}(\ell)$ increases with $\ell$, posing a problem from a numerical perspective: 1) the time scales involved are enormous; 2) the cascading towards high $\ell$ 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. 2) 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. 3) 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 $\ell=m$ . Nevertheless our scenario provides the first plausible example of a system with asymptotically flat boundary conditions, where WCCC is violated.

I Supplemental Material

WKB expansion – We present the details of the expansion in large $\ell$ that we use to determine the scalar instability growth rate. We work with $l-m=n$, where $n\sim\mathcal{O}(1)$, $n\in\mathbb{N}_{0}$. To this end we rewrite the Klein-Gordon equation from the main text in the following way

[TABLE]

where $\tilde{k}=\tilde{\mu}^{2}-\tilde{\omega}^{2}$, $\Lambda$ is the angular separation constant and we have introduced dimensionsless variables

[TABLE]

As explained in the main text, finding unstable modes amounts to determining the values of $\tilde{\omega}$ for which $\psi$ has ingoing boundary conditions at the event horizon (consistent with the equivalence principle), and finite energy on a partial Cauchy surface $t={\rm const}$ 444For unstable modes, this condition is equivalent to requiring the modes to be outgoing at spatial infinity..

In the eikonal limit $\ell\gg 1$ the spheroidal eigenvalue is $\Lambda=\ell(\ell+1)+O(\ell^{-1})$; this can easily be verified by a series expansion of the angular equation after prefactoring $\sin^{m}\theta$ from $S(\theta)$. In order to determine $\tilde{\omega}$ we divide the domain in two intersecting regions, solve (16a) in each of them and then match the solutions in the overlap. To this end we introduce the new variable $x=r/r_{+}-(1+\tilde{a}^{2})$

[TABLE]

with $\tilde{\Delta}=(1+x)(x+\tilde{a}^{2})$

Near-horizon region – Using quasimodes Holzegel:2013kna ; Keir:2014oka , one can show that $\tilde{k}=\mathcal{O}(\ell^{-2})$. We thus take the near-horizon region to be defined by $x\ll\ell^{2}$ and see that

[TABLE]

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

[TABLE]

The solution 555Multiply the equation by $(1+x)$, then change variables to $u=(x+\tilde{a}^{2})/(1-\tilde{a}^{2})$. is a linear combination of Gauss hypergeoemtric functions

[TABLE]

where $u=(x+\tilde{a}^{2})/(1-\tilde{a}^{2})$, $A_{in/out}$ are constants and

[TABLE]

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

[TABLE]

Looking at the near-horizon limit allows us to obtain

[TABLE]

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 (18), to ensure that we will get the correct asymptotic behaviour,

[TABLE]

We work with $x\gg 1$, which leaves us666First multiply (18) by $1+x$, then perform transformation (25), followed by dividing the resulting equation by $(1+x)(x+\tilde{a}^{2})(t_{0}+t_{1}\,x+t_{2}\,x^{2}+t_{3}\,x^{3}+2\,x^{4})/\left[2\,x\,(x+i_{0})^{2}\sqrt{x(x+i_{1})+i_{0}^{2}}\right]$, with $t_{0}=i_{0}^{3}(1+\tilde{a}^{2})-2\,\tilde{a}^{2}\,i_{0}^{2}+\tilde{a}^{2}\,i_{0}\,i_{1},\quad t_{1}=i_{0}(2\,i_{1}(1+\tilde{a}^{2})+2\,\tilde{a}^{2})-i_{0}^{2}(1+\tilde{a}^{2})+2\,i_{0}^{3}-\tilde{a}^{2}\,i_{1},\quad t_{2}=3\,i_{0}(1+\tilde{a}^{2}+i_{1}),\quad t_{3}=1+\tilde{a}^{2}+4\,i_{0}+i_{1},\quad i_{0}=d_{1}/(2\,\tilde{k}),\quad i_{1}=d_{2}/\tilde{k}$, and expand near spatial infinity ($x=1$) with

[TABLE]

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

[TABLE]

resulting in

[TABLE]

The general solution to (28) is then a linear combination of the Whittaker functions $W_{k,\nu}(z)$ and $W_{-k,\nu}(-z)$, thus

[TABLE]

with $B_{1/2}$ constants. Expanding for large $z$, and using

[TABLE]

we see that

[TABLE]

Imposing that the solution decays at infinity sets $B_{1}=0$.

Matching – We now perform the matching procedure in the overlapping region $1\ll x\ll\ell^{2}$. To this end, we take $x\sim\ell^{\frac{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

[TABLE]

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

[TABLE]

Afterwards we equate the coefficients in front of the equivalent terms in (32) and (33) in order to derive

[TABLE]

Furthermore, in the limit $l\rightarrow\infty$ the RHS above has very small real and imaginary parts, implying that $\Gamma[\frac{1}{2}+\nu-\kappa]$ in the denominator on the LHS must have a pole, thus

[TABLE]

with $N\in\mathbb{N}_{0}$, corresponding to the radial node of the scalar field. This allows us to deduce the real part of $\sqrt{\tilde{k}}$

[TABLE]

Next, with $\lvert\tilde{\omega}_{R}\rvert\gg\lvert\tilde{\omega}_{I}\rvert$ Yang:2012he ; Dolan:2010wr , implying $\lvert\tilde{k}_{R}\rvert\gg\lvert\tilde{k}_{I}\rvert$, we can expand (36) for large $\ell$ to obtain the real part of $\tilde{\omega}$

[TABLE]

confirming our expectations. Moreover, our calculation is accurate to $\mathcal{O}(1)$ only, hence we can replace $\tilde{\omega}$ with $\tilde{\mu}$ everywhere except inside $\tilde{k}$, where the leading order term cancels. This reduces (34) to

[TABLE]

For the imaginary part of $\tilde{\omega}$ we allow (35) to be complex

[TABLE]

with $\epsilon\ll 1$, $\epsilon\in\mathbb{C}$. We can look at the $\epsilon\rightarrow 0$ limit of (39), using (36), to derive

[TABLE]

which also enables us to find, in the limit $\ell\rightarrow\infty$,

[TABLE]

Next, we look at (38) for $\epsilon\rightarrow 0$ and obtain

[TABLE]

Taking $\ell$ large, we can rearrange for $\epsilon$, which allows us to derive an expression for $\tilde{\omega}_{I}$ via (41). It is valid for any spin parameter $\lvert a\rvert<M$ and scalar field mass $\mu$,

[TABLE]

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

[TABLE]

which is rewritten in the main text for compactness. Moreover, for $\tilde{a}=0$ the $\sinh$ changes sign and we reproduce the correct behaviour for Schwarzschild futureWork

[TABLE]

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

[TABLE]

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 (17), factor out the singular behaviour at the boundaries

[TABLE]

and finally apply Newton’s method. We then solve for $\mathscr{R}_{\omega\ell m}(x)$ and $\mathscr{S}_{\omega\ell m}(z)$.

No special treatment at the boundaries is needed as (47) force the system to pick the right solution. An initial guess is (43), $\Lambda=\ell(\ell+1)$, $\mathscr{R}_{\omega\ell m}(x)=1$, $\mathscr{S}_{\omega\ell m}(z)=1$.

Computation of $\psi_{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 - $\psi_{R}=\Re\,\psi$.

As explained in the main text, we treat the cloud as a perturbing source and use Teukolsky’s equation Teukolsky:1973ha ; Press:1973zz ; Teukolsky:1974yv

[TABLE]

to extract the NP scalar, $\psi_{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 Frolov:1998wf .

The tetrad projections we need in $T_{4}$ are $T_{nn}$, $T_{n\bar{m}}$ and $T_{\bar{m}\bar{m}}$. Moreover using a null tetrad implies we can ignore the $g_{\mu\nu}$ term. In a Kinnersley tetrad Kinnersley:1969zza we have

[TABLE]

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

[TABLE]

Next, we define the stress tensor of a scalar as

[TABLE]

and plug in $\psi_{R}=1/2(\psi+\bar{\psi})$ for $\Upsilon$ and $\Psi$, leading us to

[TABLE]

with bar indicating complex conjugation. The ansatz (50) implies that the terms above will contain exponential prefactors, which reveal that the first and second one source outgoing waves with $\hat{m}=2\,m$, $\hat{\omega}=2\,\omega$ and $\hat{m}=-2\,m$, $\hat{\omega}=-2\,\omega$, respectively, whereas the third one corresponds to the energy of the cloud, given by the gravitational mode with $\hat{m}=\hat{\omega}=0$.

The LHS of (48a) can be separated into angular and radial parts as in the vacuum case Teukolsky:1973ha . Thus take

[TABLE]

with spin coefficient $\rho=-1/(r-i\,a\,\cos\theta)$, and multiply both sides by $2\,\Sigma(r,\theta)/\rho^{4}$. $\mathcal{S}(\theta)$ can be identified with a spin-weighted spheroidal harmonic, satisfying

[TABLE]

with $c=a\,\hat{\omega}$, $\tensor{[}_{s}]{A}{{}_{\hat{\ell}\hat{m}\hat{\omega}}}$ the separation constant and $s=-2$. Multiplying by $\tensor{[}_{s}]{\bar{\mathcal{S}}}{{}_{\hat{\ell}\hat{m}}}$ and integrating over $\theta$ results in the following equation for the radial function

[TABLE]

where $K=\left(r^{2}+a^{2}\right)\hat{\omega}-a\,\hat{m}$ and the source term is the integrated over angles stress-energy tensor

[TABLE]

with the normalisation condition

[TABLE]

Equation (54) should be evaluated for $\hat{\omega}=\pm 2\,\omega$, $\hat{m}=\pm 2\,m$, $s=-2$ on the LHS with T_{\mu\nu}\big{(}\psi,\psi\big{)}, for the $+$, and T_{\mu\nu}\big{(}\bar{\psi},\bar{\psi}\big{)}, for the $-$, inside $T_{4}$ on the RHS, and for $m=\hat{\omega}=0$ on the LHS with T_{\mu\nu}\big{(}\psi,\bar{\psi}\big{)} in $T_{4}$ on the RHS.

However, the Teukolsky equation is invariant under complex conjugation followed by $m\rightarrow-m$ and $\omega\rightarrow-\omega$, hence for $\hat{m}\neq 0\neq\hat{\omega}$ we only need to determine the contribution of T_{\mu\nu}\big{(}\psi,\psi\big{)} and double the result.

Moreover, by looking at the asymptotic behaviour of the homogeneous radial equation ((54) with $\mathcal{T}_{\hat{\ell}\hat{m}\hat{\omega}}=0$) we can show that the $\hat{m}=\hat{\omega}=0$ mode is subleading to all the rest at spatial infinity. To this end, transform the homogeneous (54) with the help of

[TABLE]

and then take the limit of $r\rightarrow\infty$, leaving us with

[TABLE]

for finite777Due to the aforementioned symmetry of the Teukolsky equation, we can work, without loss of generality, with the case of $\hat{\omega}=2\,\omega$, $\hat{m}=2\,m$ only. frequency $\hat{\omega}$ and azimuthal number $\hat{m}$. The outgoing part of $Y_{\hat{\ell}\hat{m}}$ behaves as $e^{i\,r_{*}\,\hat{\omega}}r^{-s}$, leading to an $e^{i\,r_{*}\,\hat{\omega}}r^{-2s-1}$ asymptotic behaviour for $\mathcal{R}_{\hat{\ell}\hat{m}}$, which together with the definition (52) allows us to deduce that the outgoing contribution of $\psi_{4}$ near spatial infinity, for finite $\hat{\omega}$ and $\hat{m}$, behaves as

[TABLE]

However, were we to set $\hat{\omega}=\hat{m}=0$ before we take the limit of $r\rightarrow\infty$, instead of (58), we are left with

[TABLE]

where we have used $\tensor{[}_{s}]{A}{{}_{\hat{\ell}00}}=\ell(\ell+1)-s(s+1)$. The solution with finite energy at infinity corresponds to asymptotic behaviour of the form $r^{-\ell-3}$ for $\psi_{4}$, which is clearly subleading to non-zero $\hat{\omega}$ and $\hat{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 $\hat{\omega}=2\,\omega$, $\hat{m}=2\,m$ case for (54), with T_{\mu\nu}\big{(}\psi,\psi\big{)} inside $T_{4}$ on the RHS.

In this way, $\psi_{4}$ has been projected onto a basis of spin-weighted spheroidal harmonics

[TABLE]

where the definition of $\tensor{[}_{s}]{\mathcal{S}}{{}_{\hat{\ell}\hat{m}}}(\theta)$ requires $\hat{\ell}\geq\hat{m}$.

Thus, given a Kerr BH, for fixed $M\,\mu$ and $J/M^{2}$, and its $\ell=m$ scalar superradiant spectrum, we can use the latter as a source on the RHS of (54), which sets the values of the GW frequency $\hat{\omega}$ and mode-number $\hat{m}$, letting us solve (54) for any allowed value of $\hat{\ell}$.

Knowing $\psi_{4}$, we can compute the gravitational radiation at infinity and the maximum of the curvature component represented by the NP scalar $\psi_{4}$ with the formulae given in the main text. The argument for ignoring the $\hat{m}=\hat{\omega}=0$ mode in the computation of the gravitational radiation at infinity does not hold for the maximum of $\psi_{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 (53) and (54) as for (16a) and (16b).

The angular equation – (53) is an eigenvalue problem which we integrate using Newton’s method. We write

[TABLE]

where $\iota_{\pm}=\lvert\frac{m\pm s}{2}\rvert$ and solve for $\tensor{[}_{s}]{\mathbb{S}}{{}_{\ell m}}$. We drop the hats on $m$ and $\omega$, as we solve the equation generally. The eigenvalue $\tensor{[}_{s}]{A}{{}_{\ell m\omega}}$ is treated as an unknown and we normalise $\tensor{[}_{s}]{\mathbb{S}}{{}_{\ell m}}$ to $1$ using (56). Appropriate boundary conditions are imposed at the edges of the domain. These can be obtained by expanding the equation for $\tensor{[}_{s}]{\mathbb{S}}{{}_{\ell m}}$ in power series around the poles. At $z=0$ we require

[TABLE]

and similarly at $z=1$

[TABLE]

A starting point to our iteration is the series solution, provided in Berti:2005gp , up to twelfth order.

The radial equation – (54) 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 (54) for any $\hat{\ell}$. We take out the prefactors in (47) from the scalar field inside $T_{4}$, so that we can directly substitute our numerical solutions $\mathscr{R}_{\omega\ell m}(x)$ and $\mathscr{S}_{\omega\ell m}(z)$. We then factor out the singular behaviour of $\mathcal{R}_{\hat{\ell}\hat{m}}$ at infinity and the horizon, derived via Frobenius analysis, and transfer them to the RHS. Specifically

[TABLE]

with $\mathbb{R}_{\hat{\ell}\hat{m}}(x)$ dimensionless and

[TABLE]

Moreover, $x^{2\,\hat{\gamma}}$ gets cancelled by the $x^{\gamma}$ from the radial scalar functions (47) inside $T_{4}$. We then solve for $\mathbb{R}_{\hat{\ell}\hat{m}}(x)$. We need to impose boundary conditions at spatial infinity, since Teukolsky:1973ha reveals that the ingoing behaviour of the radial function (the one we do not want) is subleading

[TABLE]

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)}_{\hat{\ell}\hat{m}}=0$ is equivalent to demanding

[TABLE]

where $\hat{\varpi}=2\,\tilde{\omega}$.

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

[TABLE]

where we have used (56), the orthogonality of the spin-weighted spheroidal harmonics, and have multiplied by 2 to accound for the GWs with $\hat{m}=-2\,m$, $\hat{\omega}=-2\,\omega$.

Numerical convergence – As we cannot integrate (54) for infinitely many $\hat{\ell}$, we truncuate the sum until (69) converges. We will look at the $m=4$ case here ($\hat{m}=8$ for the GWs), as this was the hardest one to tackle numerically. We include 9 GW modes, $\hat{\ell}=8$ to $\hat{\ell}=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. (4) and (5)

Analytic approximation to $\chi$ – In this section we outline the steps for the WKB approximation of $\chi\equiv\max_{r,\theta}\;\left(|\psi_{4}|^{2}/M_{s}^{2}\right)$. We work with $\ell=m$ and $\ell\gg 1$.

First, we need a WKB expansion of the scalar field near its maximum - i.e. $R_{\omega\ell m}(r)$, $S_{\omega\ell m}(\theta)$, $\Lambda$ and $\tilde{\omega}$, as well as the location of the minimum of the potential of the radial KG equation expanded for large $\ell=m$.

We start with the angular equation (16b), which rewriten in our numerical coordinates (46), together with the redefinition (47), looks like

[TABLE]

Expand the variables in series in $m$

[TABLE]

plug them in and solve order by order, requiring regularity of the solution. The known behaviour of $\tilde{\omega}$ from before is utilised. In this way one can obtain $S_{i}(z)$ and $\lambda_{i}$ to any order as a function of the $\omega_{i}$.

Afterwards we transform the radial part (16a) using

[TABLE]

The equation will not be given here - the above transformation correponds to the first two terms in the expansion for large $\ell=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 $\Lambda=m(m+1)+\mathcal{O}(1)$, for which no knowledge of $\omega_{3}$ or higher is required. Afterwards, we take

[TABLE]

and solve order by order for $f(u)$, $R_{i}(u)$ and $\omega_{i}$, again demanding regularity of the solution. This can be done to any order, but we only give the first few here

[TABLE]

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 (54). The expressions are very lengthy and not illuminating at all, so we will give only the most important ones. Treating the angular equation (53) as the scalar one produces

[TABLE]

whereby we have accounted for $\hat{m}=2\,m$ and $d_{1}$ is a constant that can be chosen to make the normalisation of the function easier888Depending on how one does the numerics, $\tensor{[}_{s}]{\bar{\mathcal{S}}}{{}_{\hat{\ell}\hat{m}}}$ inside $T_{4}$ might also need to be normalised.. Afterwards we need to obtain a series expansion in $\ell=m$ for $\mathcal{T}_{\hat{\ell}\hat{m}\hat{\omega}}$ (55). The normalisation $\eta_{\hat{\ell}\hat{m}}$ has to be obtained by plugging (75) in (56) and integrating. The final result is

[TABLE]

With all ingredients present, we change variables in (54) with (72) and WKB expand $\mathcal{R}_{\hat{\ell}\hat{m}}$ for large $\ell=m$. Matching order by order resutls in

[TABLE]

allowing us to to rebuild $\psi_{4}$ for $\ell=m$ via (61).

We also need an expression for the scalar cloud energy

[TABLE]

Working with coordinates $z$ and $u$, we substitute our expansions for the scalar field in the definition of $\mathcal{M}_{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 $\ell=m$ and integrated over $u$ from $\frac{\tilde{\mu}^{2}(1+\tilde{a}^{2})}{2\,m\,(m+1)}$ to infinity, producing

[TABLE]

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

Having $\psi_{4}$ and $\mathcal{M}_{s}$ in the large $\ell=m$ limit allows one to derive the expression for $\chi_{WKB}$ in the main text by dividing them and expanding in series in $\ell=m$.