End point of nonaxisymmetric black hole instabilities in higher dimensions

TL;DR

This study investigates the nonlinear stability and end states of higher-dimensional rotating black holes, revealing stability in five dimensions and complex instability behaviors in six and seven dimensions, including potential violations of cosmic censorship.

Contribution

It provides the first nonlinear analysis of nonaxisymmetric instabilities in higher-dimensional black holes, showing how these black holes evolve and settle or develop singularities depending on their spin.

Findings

All five-dimensional singly spinning Myers-Perry black holes are stable.

In six and seven dimensions, black holes are stable only at low spins.

High-spin black holes develop instabilities leading to potential naked singularities.

Abstract

We report on the end state of nonaxisymmetric instabilities of singly spinning asymptotically flat Myers-Perry black holes. Starting from a singly spinning black hole in D=5,6,7 dimensions, we introduce perturbations with angular dependence described by m=2, m=3, or m=4 azimuthal mode numbers about the axis of rotation. In D=5, we find that all singly spinning Myers-Perry black holes are stable, in agreement with the results from perturbation theory. In D=6 and 7, we find that these black holes are nonlinearly stable only for sufficiently low spins. For intermediate spins, although the m=2 bar mode becomes unstable and leads to large deformations, the black hole settles back down to another member of the Myers-Perry family via gravitational wave emission; surprisingly, we find that all such unstable black holes settle to the same member of the Myers-Perry family. The amount of energy radiated into gravitational waves can be very large, in some cases more than 30% of the initial total mass of the system. For high enough spins, the m=4 mode becomes the dominant unstable mode, leading to deformed black holes that develop local Gregory-Laflamme instabilities, thus forming a naked singularity in finite time, which is further evidence for the violation of the weak cosmic censorship conjecture in asymptotically flat higher-dimensional spacetimes.

Tables

Table 1. Table 1: Real and Imaginary parts of the leading quasinormal modes in the ( ℓ 4 , m ) = ( 2 , 2 ) subscript ℓ 4 𝑚 2 2 (\ell_{4},m)=(2,2) sector of gravitational perturbations. Left: from perturbing the initial black hole. For a / μ 1 / 3 = 0.7 𝑎 superscript 𝜇 1 3 0.7 a/\mu^{1/3}=0.7 , black holes are stable, while for a / μ 1 / 3 ≥ 0.8 𝑎 superscript 𝜇 1 3 0.8 a/\mu^{1/3}\geq 0.8 black holes are unstable. Right : leading QNM in the ring-down phase. For all initial values of a / μ 1 / 3 𝑎 superscript 𝜇 1 3 a/\mu^{1/3} that end up settling back to an MP BH, we find that the end state appears to be the same black hole. For a / μ 1 / 3 = 1.5 𝑎 superscript 𝜇 1 3 1.5 a/\mu^{1/3}=1.5 , the end point of the instability is a naked singularity. We estimate an error of ± 0.02 plus-or-minus 0.02 \pm 0.02 in the real and imaginary parts of the frequencies (see footnote 5 ). For the a h / μ h 1 / 3 = 0.8 subscript 𝑎 h superscript subscript 𝜇 h 1 3 0.8 a_{\rm h}/\mu_{\rm h}^{1/3}=0.8 case, which is the lowest-spin case we consider in which the spinning black hole is unstable, the timescale of the growing mode is the longest. In practice, we find that the decay to the final state does not become evident until around t / μ 1 / 3 ∼ 200 similar-to 𝑡 superscript 𝜇 1 3 200 t/\mu^{1/3}\sim 200 , by which time the signal is already spoiled by accumulated numerical error at the typical resolutions with which we run. We found extracting the leading QNM for this case to be prohibitively expensive and have left its values as “xxx”.

Initial state		Final state		spin loss
$a_h/{\mu }_h^{1/3}$	${\omega }_{22}{\mu }_h^{1/3}$	$a_h/{\mu }_h^{1/3}$	${\omega }_{22}{\mu }_h^{1/3}$	$\Delta \left(a_h/{\mu }_h^{1/3}\right)$
$0.7$	$1.07-0.03i$	$0.7$	$1.07-0.03i$	—
$0.8$	$1.03+0.03i$	xxx	xxx	xxx
$0.9$	$1.03+0.07i$	$0.66$	$1.04-0.03i$	26.6%
$1.0$	$1.00+0.12i$	$0.64$	$1.08-0.05i$	36.0%
$1.1$	$0.95+0.17i$	$0.64$	$1.06-0.05i$	41.8%
$1.2$	$0.91+0.21i$	$0.63$	$1.05-0.05i$	47.5%
$1.3$	$0.88+0.24i$	$0.63$	$1.05-0.05i$	51.5%
$1.4$	$0.82+0.26i$	$0.63$	$1.05-0.05i$	55.0%
$1.5$	$0.74+0.28i$	$0.62$	$1.03-0.05i$	58.6%

Table 2. Table 2: Real and imaginary parts of the subleading ( ℓ 4 , m ) = ( 4 , 2 ) subscript ℓ 4 𝑚 4 2 (\ell_{4},m)=(4,2) and ( ℓ 4 , m ) = ( 4 , 4 ) subscript ℓ 4 𝑚 4 4 (\ell_{4},m)=(4,4) quasinormal mode for m = 2 𝑚 2 m=2 perturbation of a MP BH with the initial configuration a / μ 1 / 3 = 1.3 𝑎 superscript 𝜇 1 3 1.3 a/\mu^{1/3}=1.3 and a / μ 1 / 3 = 1.5 𝑎 superscript 𝜇 1 3 1.5 a/\mu^{1/3}=1.5 . We estimate an error of ± 0.04 plus-or-minus 0.04 \pm 0.04 in the real and imaginary parts of the frequencies (see footnote 5 ).

Initial state			Final state
$a_h/{\mu }_h^{1/3}$	${\omega }_{42}{\mu }_h^{1/3}$	${\omega }_{44}{\mu }_h^{1/3}$	$a_h/{\mu }_h^{1/3}$	${\omega }_{42}{\mu }_h^{1/3}$	${\omega }_{44}{\mu }_h^{1/3}$
$1.3$	$1.70+0.36i$	$1.61+0.40i$	$0.63$	$1.68-0.12i$	$1.83-0.13i$
$1.5$	$1.45+0.40i$	$1.57+0.45i$	$0.62$	$1.66-0.13i$	$1.87-0.17i$

Equations

$$d{s}^2=$$

$$+\Sigma d{\theta }^2+(r^2+a^2){\sin }^2\theta d{\varphi }^2+r^2{\cos }^2\theta d{\Omega }_{(D-4)}^2,$$

$$r=\rho {\left[1+\frac{1}{4}{\left(\frac{r_h}{\rho }\right)}^{D-3}\right]}^{\frac{2}{D-3}}$$

$$x=\rho \sin \theta \cos \varphi ,$$

$$y=\rho \sin \theta \sin \varphi ,$$

$$z=\rho \cos \theta \cos {\phi }_1,$$

$$w_1=\rho \cos \theta \sin {\phi }_1\cos {\phi }_2,$$

$$\dots$$

$$w_{D-4}=\rho \cos \theta \sin {\phi }_1\dots \sin {\phi }_{D-3}.$$

$$\chi =$$

$$\displaystyle\times\exp\left[-\left(\frac{\textrm{atanh}(\chi_{0})}{\textrm{atanh}(\chi_{H})}-\frac{\textrm{atanh}(\chi_{H})}{\textrm{atanh}(\chi_{0})}\right)^{2}\right]\Bigg{\}}\,,$$

$$m=2:f_2(x,y)=x^2-y^2,$$

$$m=3:f_2(x,y)=x^3-3x{y}^2,$$

$$m=4:f_4(x,y)=x^4-6x^2{y}^2+y^4.$$

$$({\partial }_t-{\beta }^i{\partial }_i)\alpha =c_{\alpha }\alpha \left(K-2\Theta \right),$$

$$({\partial }_t-{\beta }^j{\partial }_j){\beta }^i=c_{\beta }\left({\hat{\Gamma }}^i-f(t,\rho ){\bar{\Gamma }}^i\right)-\eta {\beta }^i,$$

$$f(t,\rho )=\exp \left[-{\left(\frac{t}{{\mu }^{1/(D-3)}}\right)}^2\left({\delta }_1{\left(\frac{{\rho }_h}{\rho }\right)}^2+{\delta }_2\right)\right],$$

$$h_{+}=\frac{{\tilde{\gamma }}_{xx}-{\tilde{\gamma }}_{yy}}{2}{\left(\frac{z}{{\mu }^{1/(D-3)}}\right)}^{\frac{D-2}{2}},$$

$$h_{\times }={\tilde{\gamma }}_{xy}{\left(\frac{z}{{\mu }^{1/(D-3)}}\right)}^{\frac{D-2}{2}}.$$

$$\displaystyle\Omega^{\prime}_{(\hat{\imath})(\hat{\jmath})}=~{}\frac{1}{2}\Big{[}R_{0k0l}m^{k}_{(\hat{\imath})}m^{l}_{(\hat{\jmath})}-R_{mk0l}m^{m}_{(\hat{1})}m^{k}_{(\hat{\imath})}m^{l}_{(\hat{\jmath})}$$

$$\displaystyle~{}-R_{0kml}m^{k}_{(\hat{\imath})}m^{m}_{(\hat{1})}m^{l}_{(\hat{\jmath})}+R_{mknl}m^{m}_{(\hat{1})}m^{k}_{(\hat{\imath})}m^{n}_{(\hat{1})}m^{l}_{(\hat{\jmath})}\Big{]}\,,$$

$$\displaystyle\Omega^{\prime}_{(\hat{a})(\hat{b})}=~{}\delta_{(\hat{a})(\hat{b})}\frac{1}{2\gamma_{ww}}\Big{[}R_{w0w0}-2\,R_{w0wk}m^{k}_{(\hat{1})}$$

$$\displaystyle-R_{wkwl}m^{k}_{(\hat{1})}m^{l}_{(\hat{1})}\Big{]}\,,$$

$$m_{(1)}$$

$$m_{(2)}$$

$$m_{(3)}$$

$$m_{(\hat{1})}$$

$$m_{(\hat{2})}$$

$$m_{(\hat{3})}$$

$$m_{(\hat{a})}=\frac{1}{\sqrt{{\gamma }_{ww}}}{\partial }_{w_a}.$$

$${\Omega }_{\ell \dots }^{\prime }=\underset{r\to \infty }{\lim }{r}^{\frac{D-2}{2}}\int d{\Omega }_{(n)}{Y}_{\ell \dots }^{(A)(B)\ast }{\Omega }_{(A)B)}^{\prime }$$

$$\dot{M}(u)=-\underset{r\to \infty }{\lim }\frac{r^{D-2}}{8\pi }\int d{\Omega }_{(n)}{I}_{(A)(B)}{I}^{(A)(B)}$$

$$I_{(\alpha )(\beta )}={\int }_{-\infty }^u{\Omega }_{(A)(B)}^{\prime }d\tilde{u}.$$

$$\dot{M}(u)=-\frac{1}{8\pi }\sum _{\ell \dots }{\left({\int }_{-\infty }^u{\Omega }_{\ell \dots }^{\prime }(\tilde{u},r)d\tilde{u}\right)}^2.$$

$$r_{\mathrm{h}}={\left(\frac{A}{C}\frac{2\pi }{{\Omega }_{D-2}}\right)}^{\frac{1}{D-3}},{\mu }_{\mathrm{h}}=\frac{C}{2\pi }{r}_{\mathrm{h}}^{D-4},a_{\mathrm{h}}=r_{\mathrm{h}}\sqrt{\frac{{\mu }_{\mathrm{h}}}{r_{\mathrm{h}}^{D-3}}-1}$$

$$E_{\mathrm{avail}}=1-\frac{M_{\mathrm{irr}}}{M},$$

$$\frac{M_{\rm irr}}{M}=\frac{r_{+}^{d-3}}{\mu}\Bigg{(}1+\left(\dfrac{a}{r_{+}}\right)^{2}\Bigg{)}^{\frac{d-3}{d-2}}\,,$$

$$d{s}^2=-\mathcal{A}{e}^{\mathcal{B}}d{u}^2-2e^{\mathcal{B}}dudr+r^2{h}_{\alpha \beta }(d{\varphi }^{\alpha }+{\mathcal{C}}^{\alpha })(d{\varphi }^{\beta }+{\mathcal{C}}^{\beta }),$$

$$h_{\alpha \beta }(u,r,{\varphi }^{\gamma })={\omega }_{\alpha \beta }({\varphi }^{\gamma })+\sum _{s\ge 1}\frac{h_{\alpha \beta }^{(s)}(u,{\varphi }^{\gamma })}{r^{D/2+s-2}}.$$

Full text

I Introduction

Black holes play a central role in general relativity (GR), the currently accepted classical theory of gravity. The recent direct detections of gravitational waves from black hole binary mergers Abbott et al. (2016) together with the image of the shadow of the black hole at the center of the M87 galaxy by the Event Horizon Telescope Akiyama et al. (2019), have changed the perception of these objects from the purely mathematical to the “tangible”. All these observations are compatible with an equilibrium (or quasiequilibrium) state that is described by the Kerr black hole Kerr (1963). However, the astrophysical relevance of the Kerr black hole hinges on whether it is nonlinearly stable or not. All the evidence indicates that it is indeed stable, but a rigorous proof of the nonlinear stability of the Kerr black hole is not yet available.

Following the discovery of the Gregory-Laflamme (GL) instability Gregory and Laflamme (1993) of black strings and black branes in dimensions 5 and higher, black holes have become useful as laboratories to study dynamical instabilities, particularly due to their simplicity and their fundamental role in GR. The study of the evolution of black hole instabilities in the fully nonlinear regime has been a fruitful area of research; indeed, the pioneering work of Lehner and Pretorius Lehner and Pretorius (2010), with numerical simulations of the evolution of the GL instability of black strings in five dimensions, provided evidence that the weak cosmic censorship conjecture is false in asymptotically Kaluza-Klein spaces.

In higher-dimensional asymptotically flat spaces, the situation is qualitatively similar. Black rings Emparan and Reall (2002); Pomeransky and Sen’kov (2006) are asymptotically flat rotating black holes with horizon topology $S^{1}\times S^{n}$. In the limit of large angular momentum along the $S^{1}$, black rings resemble thin (boosted) black strings, and hence they ought to be unstable under a GL type of instability. This was confirmed in Ref. Santos and Way (2015). Reference Figueras et al. (2016) used numerical general relativity to simulate the nonlinear evolution of black ring instabilities and showed that thin enough black rings evolve into naked singularities in finite time, thus potentially violating the weak cosmic censorship in higher-dimensional asymptotically flat spaces.111There are several issues that still need to be addressed, in particular, the structure of null infinity. In these dynamical spacetimes that develop naked singularities, it is not clear whether null infinity is incomplete. This situation is not unique to black rings. In fact, a general picture of the stability/instability properties of higher-dimensional black holes has emerged. Emparan and Myers Emparan and Myers (2003) noticed that rapidly rotating black holes in higher dimensions, i.e., Myers-Perry (MP) black holes (BHs), resemble black membranes and hence should also be GL unstable. Higher-dimensional black holes generically admit a regime of large angular momentum with highly deformed horizons described by largely separated length scales, captured by the so-called black folds Emparan et al. (2009), which are expected to be dynamically unstable. To complete this general picture, it is thus important to determine the end point of these black hole instabilities in the rapidly spinning regime.

The instability of MP BHs was confirmed in Ref. Dias et al. (2009) for linearized axisymmetric perturbations on top of a MP BH background.222MP BHs with equal spins on all rotation planes are also known to be unstable Dias et al. (2010a). For these black holes, unlike the singly spinning ones, the total angular momentum is bounded and hence the endpoints of their instabilities can potentially be quite different. We will not consider MP BHs with equal spins in this article. A fully nonlinear evolution of this axisymmetric instability of MP BHs revealed a sequence of concentric rings connected by black brane segments that became progressively thinner, eventually leading to a naked singularity in finite asymptotic time Figueras et al. (2017). In the MP case, the dynamics of the horizon is not self-similar, unlike what was seen in the black string. The study of nonaxisymmetric instabilities was initiated by Shibata and Yoshino in Refs. Shibata and Yoshino (2010a, b). For a MP BH in $D$ dimensions with mass parameter $\mu$ and spin parameter $a$, these nonaxisymmetric instabilities set in at smaller dimensionless spin $a/\mu^{1/(D-3)}$ than the axisymmetric instabilities.

Rotating black holes become dynamically unstable for values of $a/\mu^{1/(D-3)}\sim O(1)$. Indeed, from thermodynamic considerations of the case in which a MP BH fragments and expels two nonspinning BHs, Emparan and Myers Emparan and Myers (2003) estimated that nonaxisymmetric instabilities should set in at around $a/\mu^{1/(D-3)}\approx 1$. References Shibata and Yoshino (2010a, b) found that MP BHs with $a/\mu^{1/(D-3)}\lesssim 0.7$ are stable in $D=5,6,7$, but for higher spins, the MP BHs are linearly unstable to a deformation of which the azimuthal angle dependence is $e^{im\phi}$ for $m=2$, i.e., a bar-mode deformation. In six and seven dimensions Shibata and Yoshino (2010b), they found that, due to gravitational wave emission, this bar-mode instability saturates and eventually damps, and the black hole settles down to a MP BH with a lower spin. A study of quasinormal modes of MP BHs Dias et al. (2014) corroborates these results, except in the five-dimensional (5D) case; while Ref. Shibata and Yoshino (2010a) found an instability, Ref. Dias et al. (2014) did not find an exponentially growing mode in the linear regime. These two results could be compatible with each other if in five dimensions the instability were nonlinear. In this article, we resolve the apparent tension between these linear and the nonlinear results. Analogous to the GL instability of black strings, for which higher harmonics become unstable for sufficiently thin strings, one would expect that for sufficiently large $a/\mu^{1/(D-3)}$ modes with $m>2$ become the dominant unstable modes. This is indeed the case for black rings Figueras et al. (2016), and we show here that the same happens for MP BHs.

In this paper, we investigate the nonlinear evolution of MP BHs in $D=5,6,7$ dimensions with dimensionless spins of $0.7\leq a/\mu^{1/(D-3)}\leq 1.5$, and perturbed by nonaxisymmetric deformations described by the $m=2$, $m=3$, and $m=4$ azimuthal mode numbers. We use the cartoon method Pretorius (2005); Shibata and Yoshino (2010b); Cook et al. (2016) to impose an SO($D-3$) symmetry that still captures the nonaxisymmetric instability while allowing us to restrict to $3+1$ dynamics in $D$ dimensions. We believe that this symmetry assumption should still be enough to capture the essential physics of the nonaxisymmetric instabilities and their end points.

The rest of this article is organized as follows: in Sec. II, we provide an overview of the numerical methods that we used in our simulations. Section III constitutes the bulk of the article, and there we present our results for the nonlinear evolution of unstable MP BHs in various spacetime dimensions and for different values of the dimensionless spin parameter. We summarize our results and conclude in Sec IV. Technical results are collected in the Appendixes. In Appendix A we review the notion of asymptotic flatness in higher dimensions introduced in Ref. Godazgar and Reall (2012) and how the matrix of Weyl scalars captures gravitational radiation. In Appendix B we collect several properties of tensor spherical harmonics. Appendices C and D contain the tensor harmonics on the $S^{3}$ and the $S^{4}$, respectively, that we have used in our extractions. In Appendix E we review the transformation properties of the multipoles of the Weyl tensor under changes of basis. In Appendix F we compare the contours of $\chi$ with the apparent horizon, in Appendix G we display contours of $\chi$ from the evolution of six dimensions with an $m=3$ perturbation, and Appendix H contains some convergence tests.

II Numerical Methods

In this section, we summarize the numerical methods that we have employed in our simulations. In Sec. II.1 we describe our evolution scheme together with the our choice of initial data and gauge evolution equations. In Sec. II.2 we provide some details about our approach to extracting gravitational waves in higher dimensions.

II.1 Evolution scheme

The results presented here are obtained by solving the Einstein field equations in the conformal and covariant Z4 (CCZ4) formulation Alic et al. (2012); Weyhausen et al. (2012). We use Cartesian coordinates to solve for the evolution of asymptotically flat black hole solutions in $D=5,6,7$ dimensions while imposing an SO($D-3$) symmetry using the modified cartoon method Pretorius (2005); Shibata and Yoshino (2010b); Cook et al. (2016). We redefine the constraint damping parameter $\kappa_{1}\rightarrow\kappa_{1}/\alpha$ as in Ref. Alic et al. (2013), in which $\alpha$ is the lapse function, and we typically use constraint damping values $\kappa_{1}=0.4,\kappa_{2}=0$. As in Ref. Figueras et al. (2017), we use initial data for a singly spinning MP BH,

[TABLE]

where $\Sigma=r^{2}+a^{2}\cos^{2}\theta$, $\Delta=r^{2}+a^{2}-\mu\,r^{5-D}$, and $\mu$ and $a$ are the mass and spin parameters of the MP BH. In analogy with the transformation from Schwarzschild coordinates to isotropic coordinates in $D$ dimensions, we define a new radial coordinate $\rho$,

[TABLE]

where $r_{h}$ is the largest real root of $\Delta(r_{h})=0$. Then, the Cartesian coordinates are given by

[TABLE]

Imposing the SO$(D-3)$ symmetry corresponds to working on the slice $w_{1}=\ldots=w_{D-4}=0$; see Refs. Pretorius (2005); Shibata and Yoshino (2010b); Cook et al. (2016).

For values of $a/\mu^{1/3}\sim 1.5$, numerical noise is enough to trigger the instability in the $m=4$ sector. Note that although the initial (rapid) gauge adjustment induces a small burst of radiation that in practice contains modes with different $m$, this initial burst is induced by truncation error and hence is under control and small. For smaller values of the dimensionless spin, we trigger the instability by hand through an $m=2$ or $m=4$ deformation of the conformal factor $\chi$:

[TABLE]

where $\chi_{0}$ is the conformal factor computed from the analytic initial data (1), $\chi_{H}$ is some value of the conformal factor close to the apparent horizon, $A$ is the amplitude of the perturbation, and $f_{m}(x,y)$ is the function that induces the desired deformation:

[TABLE]

Introducing an $m=2$ perturbation immediately at $t/\mu^{\frac{1}{D-3}}=0$ results in a dumbbell configuration that eventually settles down to a MP BH. Unless specifically stated, simulations described here feature perturbations introduced at $t/\mu^{\frac{1}{D-3}}=0$. For some simulations, we also introduce this deformation in $\chi$ after the gauge adjustment period at around $t/\mu^{\frac{1}{D-3}}\sim 10$. For higher spins, this amounts to choosing a different branch of the dynamical evolution. Perturbing after the gauge adjustment period at around $t/\mu^{\frac{1}{D-3}}\sim 10$ results in an elongated configuration that eventually develops sharp edges at the ends. This last case is discussed in detail near the end of Sec. III.2.3.

This method of perturbing the black hole violates the constraints, so even though the CCZ4 formulation quickly333In practice this timescale is exponential and much faster than the physical timescale of the instabilities. takes us back to the constraint surface, we do not have control over where on the constraint surface we land. Therefore, to make sure that the mass and angular momentum of the perturbed data are not too far from those of the original MP BH we keep $A$ small (a typical value that we use is $A=0.02$).

We evolve the lapse with the standard $1+\log$ slicing condition Alic et al. (2012),

[TABLE]

where $c_{\alpha}$ is a freely adjustable coefficient. In our runs we found that $c_{\alpha}=1.5$ worked well for $D=5,6,7$.444The preferred choice in the typical four-dimensional astrophysical setting is $c_{\alpha}=2$. To evolve the shift, we used the modified Gamma-driver condition introduced in Ref. Figueras et al. (2016),

[TABLE]

where $\hat{\Gamma}^{i}$ is the usual CCZ4 evolution variable, $\bar{\Gamma}^{i}$ is the contracted Christoffel symbol computed from the (conformally rescaled) initial data (1), and

[TABLE]

with $\rho_{h}=r_{h}/(4^{\frac{1}{D-3}})$ and where $\delta_{1,2}$ are two adjustable parameters that control the decay of the initial gamma function near the horizon. Since the initial data have $\tilde{\Gamma}^{i}\neq 0$, the role of the extra term in (7) is to drive the gauge towards $\tilde{\Gamma}^{i}=0$ while making sure that the right-hand side of the equation remains relatively small throughout the evolution. Typical values of these parameters in our runs are $\delta_{1}=0.2$ and $\delta_{2}=0.075$. We can freely adjust the coefficient $c_{\beta}$ in (8); $c_{\beta}=0.6$ works well for us (note that the typical value in four-dimensional simulations is $c_{\beta}=0.75$). Following Ref. Figueras et al. (2016), we introduce diffusion terms (well inside the apparent horizon) on the right-hand side of the equations of motion for those variables that appear with second spatial derivatives. This improves the stability of our code, especially in the rapidly spinning regime. See Ref. Figueras et al. (2016) for more details. We use Kreiss-Oliger dissipation Kreiss et al. (1973) to damp unphysical high-frequency modes that can arise at grid boundaries or during regridding, with typical dissipation values of $\sigma=0.4$. We numerically solve the CCZ4 equations using the GRChombo libraries Clough et al. (2015); Adams et al. (2015), using up to nine levels of refinement, with a 2:1 refinement ratio, and typically $140^{3}$ points on the coarsest grid. We consider a cubic computational domain, with side length $L=200$. The $x$ and $y$ are standard Cartesian coordinates, with infinite range, which we cut off at some finite values; we typically impose periodic boundary conditions (BCs) for simplicity, but we have also experimented with Sommerfeld BCs, and the results are essentially the same. Needless to say, our choice of BCs limits the duration of our simulations. The $z$ coordinate has range $0\leq z\leq L$; at $z=0$, we impose regularity as derived from the cartoon method, and at $z=L$, we impose either periodic or Sommerfeld BCs. We discretize the equations using fourth-order finite differences and integrate in time using 4th-order Runge-Kutta method (RK4). We obtain approximately third-order convergence, see Appendix H for convergence tests.

We also track the apparent horizon (AH) during several stages of the evolution, (see Ref. Figueras et al. (2016)); our apparent horizon finder assumes that the AH is a star-shaped surface. While this is sufficient to describe the AH for the low-spin runs, this assumption breaks down when the deformations are too large. In this case, one would need to use a more general parametrization of the AH along the lines of Ref. Figueras et al. (2017) or more recently Ref. Pook-Kolb et al. (2019). We leave this for future work. However, as in Ref. Figueras et al. (2016), we verify that in our working gauge certain contours of $\chi$ track the AH to within a few percent. Thus, we can use $\chi$ to get a good approximation of the location and shape of the AH. In Appendix F we verify this claim with cases in which the AH satisfies our star-shaped assumption and can thus be calculated and compared to contours of $\chi$.

II.2 Gravitational wave extraction

One of the aims of this article is to provide a general picture of higher-dimensional black hole instabilities and their end points. To this end, the gravitational waves that are emitted during the evolution of these instabilities are a valuable source of insight. We are able to do this in two ways. The first is to extract the $h_{+}$ and $h_{\times}$ components of the metric perturbations at a certain radius along the $z$ axis Shibata and Yoshino (2010a, b),

[TABLE]

As noted in these references, $h_{+}$ and $h_{\times}$ basically contain the same information. This method has the advantage that it is very easy to implement and it accurately captures the $m=2$ modes, allowing us to compare our results with the existing results in the literature Shibata and Yoshino (2010b); Dias et al. (2014). However, this approach misses the higher $m$ modes and in particular the $m=4$ mode that becomes the dominant one at larger spins. Note that the perturbative calculations of Ref. Dias et al. (2014) only consider modes with $m=2$. Therefore, for this work, we have implemented a completely general approach to gravitational wave extraction in higher dimensions based on the higher-dimensional analogues of the Newman-Penrose scalars introduced in Ref. Godazgar and Reall (2012). We follow the implementation of Ref. Cook and Sperhake (2017) to calculate the nonvanishing components of the Weyl tensor along the outgoing null rays, $\Omega^{\prime}_{(A)(B)}$. For the class of spacetimes that we consider, the nonvanishing components of $\Omega^{\prime}_{(A)(B)}$ are Cook and Sperhake (2017)

[TABLE]

where the indices $\imath,\jmath=2,3$ run along the spatial angular dimensions of the computational domain and the indices $\hat{a},\hat{b},\ldots$ run along the transverse sphere. Here $R_{0klm}$, etc.; denote the components of the spacetime Riemann tensor, $m^{i}_{(\hat{1})}$ is the unit radial vector and $m^{k}_{(\hat{\imath})}$ are the angular unit vectors. To construct an orthonormal basis of vectors, we start from

[TABLE]

and use the standard Gram-Schmidt orthornormalization procedure. Changing to spherical coordinates as in (29), one can show that as $r\to\infty$, the basis orthonormalized vectors approach

[TABLE]

The unit vectors along the transverse sphere directions are simply given by,

[TABLE]

We extract the individual modes in the gravitational wave signal by projecting $\Omega^{\prime}_{(A)(B)}$ onto a basis of tensor harmonics on the sphere at infinity:

[TABLE]

where $\ell\ldots$ denotes the set of quantum numbers that identify each of the harmonics, and $d\Omega_{(n)}$ is the volume element on the round unit sphere $S^{n}$. In this article we focus on the $n=3,4$ cases but the procedure is completely general. In practice, we carry out the extraction at several radii to verify that we are in the wave zone. While a basis for tensor harmonics on $S^{3}$ is known, see, e.g., Ref. Lindblom et al. (2017), to the best of our knowledge, our results for tensor harmonics on $S^{4}$ are new. See Appendix D. Needless to say, for the $m=2$ modes, both methods give the same frequencies and growth/decay rates within our numerical errors.555To calculate the QNMs, we first extract the exponential growth/decay rates by fitting a straight line over the envelope of local extrema in the log-linear plot of the given data. The errors associated with this step are mainly due to the particular time interval chosen. The early and late data are contaminated by junk radiation and accumulated numerical error, respectively, whereas the signal shows a transition from growing to decaying regimes for intermediary times. Next, we use the information of the growth/decay rates to obtain a purely sinusoidal signal, from which we read the oscillatory frequency via a Fourier transformation. Here, the accuracy is restricted to the wavelengths available in the signal.

II.2.1 Energy flux

Following Ref. Cook and Sperhake (2017) (and references therein), the energy emitted in form of gravitational waves is given by

[TABLE]

with

[TABLE]

In the expressions above, $u$ is the retarded time coordinate associated to the line element expressed in the Bondi-Sachs form (see Appendix A).

Taking into account the decomposition in the tensor spherical harmonics (16), Eq. (17) is expressed as

[TABLE]

In practice, we neglect the gravitational waves before the start of the simulation and we perform the time integral in terms of the “computational” time $t$, starting at $t=0$. Total energy radiated, $E_{\text{rad}}$, can be computed by performing a final time integration.

III Results

We now present results for the nonlinear evolution of nonaxisymmetric perturbations of singly spinning black holes in five and six dimensions. The dynamics of MP black holes in seven dimensions is qualitatively similar to the six dimensional (6D) case. We start by discussing the 5D case in Sec. III.1, in which we resolve the existing tension in the literature regarding the stability/instability of MP BHs with a sufficiently large spin. We move on in Sec. III.2 to study the nonlinear stability properties of 6D MP BHs of various spins: in Sec. III.2.1, we study the gravitational waveforms and extract the quasinormal modes (QNMs), in Sec. III.2.2 we study the physics of the AH, and in Sec. III.2.3 we discuss the end point of the instabilities for large angular momentum.

In the following, we will use several geometric measures of the AH to estimate the mass and spin of the black hole throughout the evolution, even though the relations between the geometry of the horizon and the physical quantities of the black hole are only valid for equilibrium configurations. The following relations between the MP parameters and the horizon are valid in the stationary configuration

[TABLE]

with $A$ the horizon’s area and $C$ its equatorial circumference. In the dynamical regime, Eq. (20) defines $r_{\rm h}$, $\mu_{\rm h}$, and $a_{\rm h}$. Note that our results scale with the chosen length scale $\mu$, related to the black hole mass of the initial configuration. As the evolution proceeds, we dynamically rescale the physical quantities according to the length $\mu_{\rm h}$. For instance, nonlinear dynamics that settles to a final MP BH in six dimensions results in a final, dimensionless spin of $(a/\mu^{1/3})_{\rm final}=a_{\rm h}/\mu^{1/3}_{\rm h}$.

Before we proceed, we remind the reader of some well-known properties of singly spinning MP BHs in $D\geq 5$ spacetime dimensions. See Ref. Emparan and Reall (2008) and references therein for more details. For fixed total mass, in five dimensions, the angular momentum of a MP BH is bounded from above by an extremal bound; configurations saturating this bound have a naked ring singularity. This situation should be contrasted with the Kerr family of solutions, for which the solution saturating the Kerr bound corresponds to an extremal black hole with a degenerate horizon. On the other hand, for $D\geq 6$, the angular momentum is unbounded from above, and for any finite angular momentum, the black hole horizon is nondegenerate and has finite area.

III.1 Case $D=5$

In Ref. Shibata and Yoshino (2010a), the nonlinear dynamics of rapidly spinning MP BHs in five spacetime dimensions were studied for the first time, under perturbations that break the axial symmetry along the rotation plane while preserving the transverse U$(1)$. They found that for values of $a/\mu^{1/2}\sim 0.87$, there exists an exponentially growing mode that induces a bar-shaped deformation on the AH, i.e., an unstable $m=2$ mode. The authors found that the growth of the instability appeared to be stronger for larger values of $a/\mu^{1/2}$. However, the nonlinear evolution of the instability was not followed, so determining the end point of the instability remained an open problem. On the other hand, Ref. Dias et al. (2014) studied perturbations of singly spinning MP BHs in $D\geq 5$ with $m\geq 2$. In the $D=5$ case, the results of Ref. Dias et al. (2014) indicate that all linear modes decay exponentially, which would lead one to conclude that MP BHs are linearly stable in five dimensions, at least for the values of the spin parameters that were studied numerically. The results of Refs. Shibata and Yoshino (2010a) and Dias et al. (2014) would be compatible if the instability were nonlinear.

Here, we revisit the results from Ref. Shibata and Yoshino (2010a) by carrying out analogous fully nonlinear simulations of perturbed MP BHs with various values of the spin parameter, and in particular for values of $a/\mu^{1/2}$ for which Ref. Shibata and Yoshino (2010a) observed an instability. When comparing our results to those in Ref. Shibata and Yoshino (2010a), one should note that, although our gauge choice is different, our method of perturbing the black hole is essentially the same as the one used in this reference.

Contrary to the results in Ref. Shibata and Yoshino (2010a), we find that all rapidly spinning black holes in five dimensions that we are able to simulate are nonlinearly stable. For a perturbed black hole with $a/\mu^{1/2}=0.89$, for instance, Fig. 1 shows the gravitational waveform extracted on the $z$ axis Shibata and Yoshino (2010a) at $z_{0}/\mu^{1/2}=21.5$. The plot below shows the $(\ell_{3},m)=(2,2)$ scalar multipole of the gravitational waveform extracted from $\Omega^{\prime}_{(A)(B)}$. From the damped oscillatory decay, we obtain a dominant QNM with frequency $\omega\mu^{1/2}=1.283-0.032\,i$, which agrees with the results from perturbation theory of  Dias et al. (2014), i.e, $\omega r_{+}=0.585-0.015\,i$. For a further check, we computed the frequency and decay rate of the leading QNM from both $h_{+}$ and $\Omega^{\prime}_{(A)(B)}$, and they agree within the truncation error.

In five dimensions, we find that gravitational waveform is essentially in the $(\ell_{3},m)=(2,2)$ scalar-derived sector: see Fig. 1. The projection onto the other tensor harmonics, i.e., vector-derived and transverse traceless tensors, gives waveforms with amplitudes that are smaller than those in the $\ell_{4}=2$ scalar sector by 1 order of magnitude or more. Higher multipoles again give much smaller waveforms.

Our results together with the results of Ref. Dias et al. (2014) suggest that 5D MP BHs are both perturbatively and nonperturbatively stable under perturbations that break the axial symmetry on the rotation plane. In this article, we do not study the stability properties of MP BHs with values of $a/\mu^{1/2}$ that are arbitrarily close to extremality. As pointed out in Ref. Yang et al. (2015), the dynamics of the Kerr black holes near extremality can exhibit turbulence, so it is conceivable that similar dynamics will arise for MP BHs in the limit $a/\mu^{1/2}\to 1$.

III.2 Case $D=6$

In this subsection we consider the nonlinear evolution of perturbed 6D MP BHs of different values of the spin parameter $a/\mu^{1/3}$. The results in seven dimensions are qualitatively similar, and we will not show them here.

Recall that in $D\geq 6$ the dimensionless spin parameter $a/\mu^{1/3}$ of singly spinning MP BHs can be arbitrarily large. Reference Emparan and Myers (2003) noted that the thermodynamics of spinning black holes exhibits a qualitative change in behavior as the angular momentum increases beyond a certain (dimension-dependent) critical value of $a_{\textrm{crit}}/\mu^{1/(D-3)}$; for values of $a/\mu^{1/(D-3)}$ smaller than this critical value, MP BHs behave similarly to the Kerr black hole, and for values of $a/\mu^{1/(D-3)}$ greater than this critical value, they behave like black membranes. In six dimensions, one finds $a_{\textrm{crit}}/\mu^{1/3}\simeq 1.29$, and the ultraspinning instability kicks in at $a/\mu^{1/3}\simeq 1.57$ Dias et al. (2010b). The end points of certain ultraspinning instability were worked out in Ref. Figueras et al. (2017). On the other hand, the nonaxisymmetric instabilities of 6D MP BHs kick in at $a/\mu^{1/3}\simeq 0.74$ Shibata and Yoshino (2010b); Dias et al. (2014). In this article we study the endpoint of nonaxisymmetric instabilities for MP BHs with $0.7\leq a/\mu^{1/3}\leq 1.5$. The results presented here, taken in conjunction with those in Ref. Figueras et al. (2017), constitute a complete picture of the dynamics of singly spinning MP BHs.

III.2.1 Gravitational waves

In Table 1, we summarize the QNM frequencies in the $m=2$ sector that we have extracted from our fully nonlinear simulations. To obtain the results displayed on this table, we considered sufficiently small perturbations, with amplitudes $A=0.02$. Although our simulations are fully nonlinear, thus making it difficult to obtain accurate results for linear waves, our results agree within the truncation error with the perturbative calculations of Dias et al. (2014). In the left column of Table 1 we display the normalized frequencies of the leading mode obtained from perturbing the initial black hole with an $m=2$ mode. This shows that black holes with $a/\mu^{1/3}=0.7$ are stable, while for $a/\mu^{1/3}\geq 0.8$, the sign of the imaginary part of the frequency indicates a linear instability. In the right column, we display the normalized frequencies of the leading QNM in the ring-down phase for those runs that settle back onto another MP BH. We measured the AH area and the length of the circumference of the AH on the equatorial plane of the final state to estimate the angular momentum parameter of the final MP black hole using (20). These results indicate that all runs settle down to what appear to be the same black hole (within the numerical error), independent of the initial value of $a/\mu^{1/3}$ that we have considered.

In Fig. 2 we display the gravitational waveforms corresponding to black holes with initial spin $a/\mu^{1/3}=1.3$ (top) and $a/\mu^{1/3}=1.5$ (bottom), extracted using the metric perturbations. Together with the waveforms, we also display snapshots of the AH at different stages of the evolution as an inset above the figure. From these AH snapshots, it is clear that the dynamics is governed by an $m=2$ mode which deforms the AH into a bar shape. The black hole radiates mass and angular momentum until it spins down and settles down to an equilibrium MP BH, with lower mass and angular momentum. As we increase the spin parameter, the systems stays for a longer period in the transient regime between the initial dynamics described by perturbation theory and the final ring-down phase.

Figure 3 shows the gravitational waveforms extracted from the Weyl scalars for both $\ell_{4}=2$ (top) and $\ell_{4}=4$ (bottom) modes in the scalar-derived sector of tensor harmonics for $a/\mu^{1/3}=1.5$. Most of the signal is in the $(\ell_{4},m)=(2,2)$ mode, but the $\ell_{4}=4$ also displays a sizeable signal. However, for MP BHs with $a/\mu^{1/3}=1.3$ the $(\ell_{4},m)=(2,2)$ mode governs the nonlinear evolution, with qualitatively similar results for $a/\mu^{1/3}=1.5$. Table 2 displays the quasinormal frequencies in the $\ell_{4}=4$ sector for $a/\mu^{1/3}=1.3$ and $a/\mu^{1/3}=1.5$.

The end point of the evolution of these instabilities depends, in a nontrivial way, on the amplitude of the initial perturbation and the value of $a/\mu^{1/3}$ for the initial black hole. For small enough perturbations and a small enough initial value of $a/\mu^{1/3}$, the end state of the nonaxisymmetric instabilities is another MP BH. Surprisingly, the perturbed black holes of this type that we have been able to simulate all appear to settle down to the same MP BH, the one with $(a/\mu^{1/3})_{\textrm{final}}=0.63$. We are not aware of the physical reason that singles out this particular member of the MP family, but our simulations suggest that it behaves as an attractor. It would be worth investigating this point further.

III.2.2 Apparent horizon

In Fig. 4, we display the evolution of the estimated black hole spin parameter from the geometry of the AH (20). This estimate only corresponds to the spin parameter of a black hole in equilibrium, and therefore it is exact only at the initial and final stages of the evolution. However, it indicates that highly spinning unstable black holes can radiate an enormous amount of angular momentum due to the large deformation of the horizon during the highly nonlinear stages of the evolution. For instance, black holes with an initial $a/\mu^{1/3}=1.5$ can radiate up to $58\%$ of the initial spin. In Fig. 5 we display four snapshots of the AH during the evolution of a perturbed MP with initial $a/\mu^{1/3}=1.3$. As can be seen from these snapshots, the AH develops a strong bar-shaped type of deformation, which leads to a very strong emission of gravitational radiation. As the plots in Fig. 4 show, the amount of angular momentum radiated increases with the initial value of $a/\mu^{1/3}$, so it is conceivable that one can achieve even higher emission percentages for larger initial values $a/\mu^{1/3}$. However, as we will discuss next, there may be a (dimension-dependent) dynamical upper bound due to the fact that the generic end point of the evolution of the instability of MP BH with a sufficiently large initial value of $a/\mu^{1/3}$ or a sufficiently large initial deformation is no longer another MP BH but is a naked singularity.

III.2.3 Large initial spin

Given the large deformations observed at moderate spins and small perturbations, e.g., see Fig. 5, one may wonder what happens when the spin of the unperturbed black hole or the size of the deformation is increased. In this article, we are interested in the end points of perturbative instabilities of black holes, so we will always keep the initial deformation small enough to remain in the perturbative regime initially and vary the spin of the MP BH. As we now show, the inevitable end point of the instabilities for sufficiently rapidly spinning black holes is a naked singularity. Before we proceed, note that as the deformations of the AH beyond spherical symmetry become increasingly large during the highly nonlinear stages of the evolution the AH ceases to be a star-shaped surface, which is an assumed property underlying the AH finder we use. The same situation was shown to occur in unstable black rings Figueras et al. (2016) or in the ultraspinning instability of MP black holes Figueras et al. (2017). The latter reference managed to find non-star shaped AH using an intrinsic parametrization of the surface, while Ref. Pook-Kolb et al. (2019) proposed the use of a reference surface. We will leave the interesting problem of implementing one of these methods to the present situation for future work. However, as in Ref. Figueras et al. (2016, 2017), we have verified that certain contours of the conformal factor $\chi$ provide a very accurate description of the AH, within less than a few percent. See Appendix F for quantitative comparison between the contours of $\chi$ and the actual AH for various runs in which the latter can be found. Therefore, we will use contours of $\chi$ as approximations of AH for situations where the latter is not a star-shaped surface, and hence get intuition of the physics of the instabilities and their end points.666Note that since contours of $\chi$ provide accurate approximations of the AH, in principle one could use them as suitable reference surfaces to find the AH using the method of Pook-Kolb et al. (2019).

For MP BH’s with an initial $a/\mu^{1/3}=1.5$, we find that numerical noise is enough to excite the $m=4$ mode, which is the one that dominates the subsequent nonlinear evolution.777The same happens for black rings: moderately thin rings are unstable under a GL mode with $m=2$, while for sufficiently thin rings the fastest growing instability is in $m=4$ sector. Presumably, for even thinner rings, higher $m$ modes are more unstable, just as in the GL instability of black strings.

Modes with $m=2$ that are present in the noise also get excited but have much smaller amplitude. See Fig. 6. Comparing the waveforms of the $m=4$ instability for $a/\mu^{1/3}=1.5$ in Fig. 6 with their counterpart waveforms for the $m=2$ instability for $a/\mu^{1/3}=1.3$ in Fig. 2, notice that the amplitude for the $m=2$ case is much larger, even though the AH deformation is quite extreme in both cases. This suggests that the efficiency of the gravitational wave emission during the evolution of the instabilities decreases as $a/\mu^{1/3}$ increases, hence giving more time for these instabilities to grow and eventually form a naked singularity. The simple physical reason for this is that the structures that form at large values of $a/\mu^{1/3}$, as deformed as they may be, contain only a small amount of the total mass.

In Fig. 7, we show representative snapshots of $\chi=0.5$ contours from the evolution of the MP BH with $a/\mu^{1/3}=1.5$. The initial deformation of the black hole horizon into a square shape is characteristic of the $m=4$ mode. Subsequently, as the instability develops, the corners of the square grow into four arms. These arms become longer and thinner as time goes by, until at some point a local GL instability kicks in along the arms, signaled by the appearance of a bulge in the central part of each of the arms. These bulges travel outward towards the tip of the arms due to centrifugal force, which leads to a further thinning of arms, and mass accumulates at the tips, forming nearly spherical bulges. It becomes very difficult to continue the simulation beyond this point since resolving the new generations of bulges that form due to the local GL instability becomes computationally very expensive.888Recall that the partial differential equations (PDE) that we are solving are effectively $(3+1)$ dimensional, while in Ref. Lehner and Pretorius (2010) or Figueras et al. (2017), the system is $(2+1)$ dimensional, thus allowing one to get much closer to the singularity. However, the local evolution of the instability along each of the arms appears to be analogous to the ultraspinning case of Ref. Figueras et al. (2017). In this case, it was observed that the interaction between the centrifugal force and mass accretion by the big bulges accelerates the formation of a naked singularity compared to the GL instability of black strings Lehner and Pretorius (2010), and the process is no longer self-similar. Note that the central bulge becomes more spherical as the evolution proceeds and in fact contains most of the mass of the system. This explains why the gravitational wave emission is less efficient at this stage of the evolution, at sufficiently large values of $a/\mu^{1/3}$. We believe that this is a generic aspect of the dynamics of unstable black holes at large angular momentum. The most likely end point of this instability is that quantum gravity effects will govern the breakup of the arms, and the resulting four black holes will either fly off to infinity or be recaptured by the central black hole. We also performed simulation with $m=3$ perturbations, in which we observe essentially the same qualitative dynamics, but with the black hole developing three arms instead of four — see Appendix G.

If we start with a black hole with initial $a/\mu^{1/3}=1.5$ and perturb it with an $m=2$ mode, we can in principle explore the nonlinear evolution of this sector of perturbations at large angular momentum. Doing so by perturbing at $t/\mu^{1/3}=0$ only results in a dumbbell 999We thank Roberto Emparan for telling us about the existence of two sectors of instabilities for the $m=2$ perturbations Andrade et al. .configuration that eventually settles back down to a MP BH, as shown in the top panel of Fig.8 for an initial $m=2$ perturbation with amplitude of $A=0.02$. However, by perturbing after the gauge adjustment period at around $t/\mu^{1/3}=10$, the perturbation turns out to be in a sector where unstable modes dominate the nonlinear regime. In this case, we find that during the highly nonlinear regime of the subsequent evolution, the AH becomes more elongated and develops very sharp features at the edges; see the bottom panel in Fig. 8. The appearance of these sharp features is not a numerical artifact and it makes the numerical simulation of the system very difficult. While we were unable to continue this particular simulation beyond this point due to the cost of the computation, one is tempted to conjecture that the end point of the evolution may well be a naked singularity. Indeed, the spacetime curvature at these sharp tips grows large, as spatial gradients diverge and the timescales become very fast. In fact, in the $m=4$ case, the appearance of the sharp tips precedes the formation of long and thin arms, which eventually become GL unstable. We suspect that the same will happen in the present $m=2$ case: it is quite possible that the evolution will continue by forming two long arms, joined at a central region that is nearly spherical and contains most of the mass. These long arms will become thinner over time due to the centrifugal force and eventually will develop local GL instabilities. Therefore, the end point of the evolution is likely the formation of a naked singularity in finite asymptotic time.

III.2.4 Energy flux

We end this section by calculating the energy emitted in the form of gravitational waves according to Eq. (19). As discussed in the previous sections (see also Appendix §D), the main contribution to wave signal comes from the $\ell_{4}=2$ and $\ell_{4}=4$ projections of the Weyl tensor into the scalar-derived tensor harmonics. The top panel of Fig. 9 compares the normalized energy flux for the $m=2$ perturbation of MP BH’s with initial spin $a/\mu^{1/3}=1.3$ and $a/\mu^{1/3}=1.5$. The peak of emission occurs around $(t-r_{0})/\mu^{1/3}=15-20$, which corresponds exactly to the period dominated by the growing QNM phase. Then, the emission rate stays relatively constant during the transient period, before dropping to zero in the final ring down phase. These three phases of the emission are evident in the gravitational wave signal depicted in Fig.2.

The orders of magnitude in the energy flux are essentially the same for both systems. However, as noticed in Sec. III.2.1, the larger the initial spin, the longer the transient period. Hence, the system emits energy for a longer time, before settling down to the final BH. Indeed, the bottom panel of Fig. 9 shows the normalized radiated energy, and we observe and emission of approximately $23.54\%$ and $31.15\%$ of the BH’s mass with initial spin $a/\mu^{1/3}=1.3$ and $a/\mu^{1/3}=1.5$, respectively. This should be contrasted with the head-on collisions or the merger of binary black holes in higher dimensions, for which less than $1\%$ of the total mass is radiated Witek et al. (2010); Cook et al. (2017, 2018) (in $D=6$, the fraction of mass radiated is $(8.19\pm 0.05)\times 10^{-4}$ in the head-on collisions and $(1.99\pm 0.05)\times 10^{-1}$ in the binary merger), or in the first black hole binary merger recorded by LIGO Abbott et al. (2016), which radiated roughly $4.6\%$ of the initial mass into gravitational waves.

It is instructive to compare the energy emitted by the gravitational waves against the total available energy given by

[TABLE]

with $M_{\rm irr}$ the irreducible mass

[TABLE]

where $r_{+}$ is the horizon location of the MP BH in the Boyer-Lindquist type of coordinates (1). The energy carried out by gravitational waves corresponds to $E_{\rm rad}/E_{\rm avail}\sim 0.60$ and $E_{\rm rad}/E_{\rm avail}\sim 0.64$ for $a/\mu^{1/3}=1.3$ and $a/\mu^{1/3}=1.5$, respectively.

Finally, we compare the energy emission for the MP BH with spin $a/\mu^{1/3}=1.5$ perturbed with $m=2$ and $m=4$ deformations. Figure 10 confirms that the efficiency of the emission is lower by 2 orders of magnitudes in the $m=4$ dynamics, which gives more time for the instabilities to grow, as discussed in the previous section. In particular, the $m=4$ perturbation depicted in Fig. 7 evolves towards a configuration with a large central bulge that contains most of the black hole’s mass and becomes progressively more spherical as the evolution proceeds. This mostly spherical configuration is thus expected to be less efficient at emitting gravitational waves than the configuration that results from the $m=2$ perturbation depicted in the top panel of Fig. 8.

IV Conclusion

In this article we describe the evolution and end points of nonaxisymmetric perturbations of singly spinning MP BHs in five and six dimensions. In the 5D case, we find that all black holes that we have considered are nonlinearly stable, thus resolving the tension between the numerical simulations of Ref. Shibata and Yoshino (2010a) and the perturbation theory results of Dias et al. (2014). In the 6D case we found that slowly spinning black holes are stable. As the initial spin increases, at around $a/\mu^{1/3}\approx 0.8$ MP BHs become unstable under perturbations with an $m=2$ azimuthal mode number (i.e., bar mode), as Ref. Shibata and Yoshino (2010b) had previously observed and Ref. Dias et al. (2014) confirmed. We have been able to evolve MP BHs with significantly larger initial spins than previous works, and the picture that has emerged is the following: in the nonlinear regime, unstable black holes develop a bar shape which leads to a strong emission of gravitational waves, through which the black hole eventually settles back down to another member of the MP family with a lower angular momentum. See Fig. 5. As the initial spin of the black hole grows, the amount of mass and angular momentum emitted by the gravitational waves during the evolution also grows, reaching a staggering $31.15\%$ of radiated mass for a black hole with initial spin $a/\mu^{1/3}=1.5$. This is much larger than previously observed in head-on collisions and merger of black holes in higher dimensions. It would be interesting to explore whether these high emission rates have implications for scenarios with large extra dimensions. Another surprising fact uncovered in this work is that all unstable MP BHs that settle back down to another member of the MP family in fact seem to settle onto the same black hole, the MP BH with $a/\mu^{1/3}=0.63$. Therefore, this particular member of the MP family seems to be an attractor, at least for the evolution of unstable BHs under a bar-mode type of deformation. We do not have a physical explanation for this behavior since the $a/\mu^{1/3}=0.63$ black hole does not to have any particular physical properties that single it out among nearby solutions. It would be very interesting to understand why this particular solution seems to be special.

For initial spins of $a/\mu^{1/3}\approx 1.5$, it is the $m=4$ mode that has the fastest growth rate and will generically dominate the subsequent evolution, even though lower $m$ modes are still present. The same situation was observed in the case of the black ring in 5D Figueras et al. (2016). The evolution of these instabilities is different than those at lower spins; the $m=4$ mode results in a square-type deformation, which then evolves into four arms (see Fig. 7). These arms grow longer and thinner over time, until they eventually develop a local GL instability. In this article we have been able to follow the evolution of this local instability beyond the formation of the first generation of bulges. After this point, the simulation becomes too computationally expensive. However, at this point the fate of these black holes is sealed: the development of the GL instability is at least exponential in time, while the emission of gravitational waves only follows a power law. Therefore, once the local GL instability kicks in, a naked singularity will inevitably form in finite asymptotic time before the unlucky black hole has had time to radiate away enough mass and angular momentum. Since this mechanism is generic, one may interpret this as further evidence that the weak cosmic censorship conjecture is false in higher-dimensional asymptotically flat spacetimes. Note that it is not clear whether the type of singularity that forms is visible to external faraway observers and hence whether it constitutes an honest violation of the conjecture. To settle this issue, one would need to analyze the (in)completeness of null infinity, which remains an open question for this class of spacetimes. However, arbitrarily large curvatures near the horizon do become visible. Since the mass contained in these regions of increasingly large curvature becomes negligible over time, the potential violation of the weak cosmic censorship conjecture would be of the mildest possible type in the language of Ref. Andrade et al. (2019). For these rapidly spinning black holes, one can also trigger an instability in the $m=2$ sector of perturbations. Even though we were not able to follow the evolution to the point of seeing the appearance of local GL instabilities and hence naked singularities, it seems that this is the most likely end point.

The instabilities discussed in this article are of the “elastic” type Figueras et al. (2016), as opposed to the GL type of instabilities. Intuitively, elastic instabilities deform the black hole horizon without introducing substantial thickness variations. GL instabilities in singly spinning MP BHs kick in at slightly larger spins, $a/\mu^{1/3}=1.572$ in $D=6$ Dias et al. (2010b); the evolution and end points of the latter were studied in Ref. Figueras et al. (2017). Whether the elastic type of instabilities or the GL ones dominate the evolution of a given black hole is a question of suitably preparing the initial conditions. We do not see any obstruction in being able to construct open sets of initial conditions for which either the elastic or GL instability dominates the subsequent evolution by exciting the desired mode with a sufficiently large initial amplitude. Therefore, this article, together with Ref. Figueras et al. (2017), provides a complete picture of the dynamics of some of the most notable instabilities that afflict MP black holes and their respective end points.101010The superradiant instability is another important instability that affects these black holes, but we have not considered it here.

One important result in higher-dimensional black hole physics that emerges from our work is that sufficiently rapidly rotating black holes in higher-dimensions are unstable and evolve into a naked singularity in finite time as a result of local GL dynamics, regardless of the nature of the original instability. Both the elastic and GL instabilities are quite generic in the sense that they should afflict other types of higher-dimensional black holes apart from MP BHs and black rings. This leads us to the following conjecture. The GL instability is the only mechanism that higher-dimensional vacuum GR has to change the topology of a black hole horizon in dynamical spacetimes. If true, this conjecture would make the detailed understanding of the GL instability an even more pressing issue. The fact that the GL instability in black strings appears to be self-similar and perhaps controlled by an attractor offers hope that some analytic progress is possible in this case.

The dynamics of singly spinning MP BHs in seven dimensions is qualitatively and quantitatively very similar to the 6D ones. Namely, BHs are nonlinear stable for sufficiently small angular momentum. At larger spins, they first become unstable under a bar-mode type of deformation. The evolution of this instability again ends on another member of the MP family, with lower mass and spin. For larger angular momentum, the $m=4$ mode becomes the fastest-growing unstable mode, and the end point is a naked singularity, which forms in finite asymptotic time. The only difference between six and seven dimensions is that naked singularities appear to form at lower values of the angular momentum. For instance, in seven dimensions, the $a/\mu^{1/4}=1.3$ MP BH is unstable under an $m=4$ mode, and the end point of the instability is a naked singularity. As a rule of thumb, it seems that forming naked singularities is “easier” as $D$ increases.

One may object that unstable BHs should not be regarded as “generic” and consequently the potential violations of the weak cosmic censorship conjecture in this and previous works should not be considered generic. First, the question of whether it is possible or not to construct open sets of initial conditions sufficiently close to the unstable black holes considered in these works remains open. Second, in this article we have shown that black holes with sufficiently large angular momentum develop local GL instabilities which inevitably end up forming naked singularities. Last but not least, using the large $D$ limit of GR, Ref. Andrade et al. (2019) provided compelling evidence that it should be possible to generically form single black holes with large angular momentum by colliding black holes with a nonzero impact parameter. Moreover, this reference also suggested that these black holes should evolve into naked singularities via local GL instabilities. In this article, we have shown that these local GL instabilities do indeed kick in when one expects them to do so. It now remains to show that indeed the scenario put forward in Ref. Andrade et al. (2019) is realized in finite dimensions. $D=6$ should suffice.

Acknowledgements.

We would like to thank Tomas Andrade, Jay Armas, Roberto Emparan, Mahdi Godazgar, Harvey Reall, and Ulrich Sperhake for discussions and the referee for a careful reading of the manuscript. We have benefited from being part of the GRChombo Collaboration (www.grchombo.org), our special thanks go to the whole collaboration. HB, PF, MK and RPM were supported by the European Research Council Grant No. ERC-2014- StG 639022-NewNGR “New frontiers in numerical general relativity”. PF is also supported by a Royal Society University Research Fellowship (Grant No. UF140319). We acknowledge that the results of this research have been achieved using the DECI resource Cartesius based in the Netherlands with support from the PRACE DECI-14-Tier-1 call. Some of the computations reported in this article were also performed at the MareNostrum4 supercomputer at the Barcelona Supercomputing Centre, Grants No. FI-2019-1-0010, FI-2018-3-0004, FI-2017-3-0012, FI-2017-2-0018, FI-2017-1-0042. We also acknowledge the use of Athena at HPC Midlands+, which was funded by the EPSRC on Grant No. EP/P020232/1, in this research, as part of the HPC Midlands+ consortium. Part of this work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital Grants No. ST/P002307/1 and ST/R002452/1 and STFC operations Grant No. ST/R00689X/1. DiRAC is part of the National e-Infrastructure.

Appendix A Asymptotic flatness and gravitational waves

We assume the spacetime to be asymptotically flat in the sense of Ref. Godazgar and Reall (2012). Thus, in the wave zone, the line element is expressed in the Bondi-Sachs form Bondi et al. (1962); Sachs (1962)

[TABLE]

with $u$ a retarded time and, $\mathcal{A}$, $\mathcal{B}$, ${\mathcal{C}}^{\alpha}$ and $h_{\alpha\beta}$ functions on the Bondi coordinate $(u,r,\phi^{\alpha})$. Moreover, $\det h_{\alpha\beta}=\det\omega_{\alpha\beta}$, where $\omega_{\alpha\beta}$ is the unit metric on the unit $n$-sphere $S^{n}$. Here, we are interested in the asymptotic expansion Tanabe et al. (2011)

[TABLE]

with $h_{\alpha\beta}^{(1)}$ the Bondi news function.

We consider the higher-dimensional generalization of the Geroch-Held-Penrose formalism Durkee et al. (2010) and introduce a tetrad basis $(\ell,k,m_{A})$ satisfying $\ell^{\mu}k_{\mu}=1$, $m_{(A)}^{\mu}m_{(B)}{}_{\mu}=\delta_{(A)(B)}$, and with vanishing further contractions. Asymptotically, one particular choice for the tetrad is

[TABLE]

Gravitational waves are extracted from the projection of the Weyl tensor $C_{\mu\nu\rho\sigma}$ onto the tetrad basis (25). Specifically, the news function $h^{(1)}_{\alpha\beta}$ is given by the leading contribution of Ref. Godazgar and Reall (2012),

[TABLE]

From the symmetries of $C_{\mu\nu\rho\sigma}$, one verifies that $\Omega^{\prime}_{AB}$ is symmetric and traceless.

Appendix B Spherical harmonics on an $n$-sphere

In this Appendix we collect some well-known properties of spherical harmonics on an $S^{n}$. In Appendixes and we specialize to the cases of interest for us, namely, $S^{3}$ and $S^{4}$.

We start by setting up our notation and conventions.111111Uppercase latin indices $A,B,...$ run from [math] to $D-1$. Lowercase latin letters $i,j,...$ are spatial indices running from $1$ to $D-1$, and greek letters denote angular indices from $2$ to $D-1$. The sphere has dimension $K=D-2$, and the computation domain has $d=3$ spatial dimensions. Following Ref. Cook and Sperhake (2017), we consider spherical coordinates $Y^{A}=(t,r,\theta_{a})$ which relate to a Cartesian coordinate system $X^{a}=(t,w_{i})$ via

[TABLE]

Here, $\theta_{1}\in[0,2\pi]$ and $\theta_{a}\in[0,\pi]$ $\forall a=2\,\ldots,n$. In the three-dimensional computational domain, we parametrize the numerical Cartesian coordinates $(x,y,z)$ by new spherical coordinates $(r,\theta,\varphi)$ via

[TABLE]

with $\theta_{n},\theta_{n-1}\in[0,\pi]$.

Now we are in position to collect some properties of scalar harmonics and higher rank tensor harmonics on an arbitrary $n$-sphere $S^{n}$, with $n>2$. Ultimately we are interested in tensor harmonics on the $S^{3}$ and the $S^{4}$. While the former are well known, see e.g., Ref. Lindblom et al. (2017) and references therein, our results for the tensor harmonics on the $S^{4}$ are, to the best of our knowledge, new.

We write down the metric $ds^{2}_{n}$ on the unit round $S^{n}$ as

[TABLE]

with $0\leq\theta_{1}\leq 2\pi$ and $0\leq\theta_{i}\leq\pi$, $\forall i=2,\ldots,n$. We refer to this particular form of the metric on the $S^{n}$ as the “standard” form and the corresponding basis of unit vectors as the standard basis.

Recall that spherical harmonics on the unit $S^{n}$ are defined as eigenfunctions of the Laplace operator on the $S^{n}$, 121212Here, we work on the unit $S^{n}$. One can easily reinstate the radius in our formulas.

[TABLE]

where the integers $\ell_{i}$ satisfy $\ell_{n}\geq\ell_{n-1}\geq\ldots\geq\ell_{2}\geq|\ell_{1}|$. The scalar harmonics are normalized so that

[TABLE]

where the volume element of $S^{n}$ is omitted to avoid clutter. From the scalar harmonics, on can construct scalar-derived vector harmonics on the $S^{n}$. These are given by

[TABLE]

and they satisfy

[TABLE]

where $g$ is the metric on the round unit $S^{n}$. On the other hand, divergence-free vector harmonics $\mathbb{V}_{a}^{\ell_{n},\ldots,\ell_{1}}$ on the $S^{n}$ satisfy

[TABLE]

One can show that there are $n-1$ linearly independent, orthogonal, and divergence-free vector harmonics on the $S^{n}$ which, together with the scalar derived vector harmonics (33), form a complete basis of vectors on the $S^{n}$, see Refs. Chodos and Myers (1984); Rubin and Ordonez (1985); Higuchi (1987).

From the scalar harmonics (31), one can obtain two types of symmetric tensor harmonics:

[TABLE]

These satisfy:

[TABLE]

On the other hand, we have

[TABLE]

Recall that the matrix of Weyl scalars, $\Omega^{\prime}_{AB}$, encoding the gravitational radiation is symmetric and traceless. Therefore, the family $\mathbb{S}_{(1)ab}$ of scalar-derived tensor harmonics cannot contribute to the multipolar expansion of the Weyl scalars. Hence, from now on we will only consider the second family, $\mathbb{S}_{(2)ab}$, of scalar-derived tensor harmonics and we shall drop the subscript (2).

From each family of vector harmonics, one can similarly construct the corresponding family of vector-derived tensor harmonics:

[TABLE]

These satisfy,

[TABLE]

Finally, transverse traceless symmetric tensor harmonics are defined by

[TABLE]

Needless to say, the different families of tensor harmonics are orthogonal to each other and together form a complete basis of tensor harmonics on the $S^{n}$.

The class of dynamical spacetimes that we consider possesses a transverse round $S^{n-2}$. Therefore, to extract the gravitational waveforms form the components of the Weyl scalars we can restrict ourselves to tensor harmonics on the $S^{n}$ that preserve an internal SO$(n-1)$ symmetry. This amounts to assuming that $\ell_{n-2}=\ldots=\ell_{1}=0$ for the quantum numbers labeling the harmonics. We shall restrict to this class of harmonics from now on.

Appendix C Tensor harmonics on the $S^{3}$

In this Appendix we collect some results on tensor harmonics on the $S^{3}$.

C.1 Standard basis

We write the round metric on the unit $S^{3}$ as:

[TABLE]

with $\theta_{3},\theta_{2}\in[0,\pi]$ and $\theta_{1}\in[0,2\pi]$. We choose the following basis of unit vectors on the $S^{3}$:

[TABLE]

In this basis, the projected components of the harmonics with $\ell_{1}=0$ are:131313Since tensor harmonics are symmetric $\mathbb{Y}^{\ell\ldots}_{21}=\mathbb{Y}^{\ell\ldots}_{12}$ so we do not need to list all components.

Scalar-derived tensor harmonics

• •

$(\ell_{3},\ell_{2})=(2,0)$:

[TABLE]

• •

$(\ell_{3},\ell_{2})=(2,1)$:

[TABLE]

• •

$(\ell_{3},\ell_{2})=(2,2)$:

[TABLE]

Vector-derived tensor harmonics

There are two families of vector-derived tensor harmonics on the $S^{3}$, but given the class of spacetimes that we consider in this article, only one of them has a nonzero overlap with the Weyl scalars. The relevant vector harmonics for us are given by the following

• •

$(\ell_{3},\ell_{2})=(2,1)$

[TABLE]

• •

$(\ell_{3},\ell_{2})=(2,2)$

[TABLE]

Transverse traceless tensor harmonics

There are two families of transverse traceless tensor harmonics on the $S^{3}$, but, as with the vectors, only one of them has a nonzero overlap with the Weyl tensor. This is given by the following

• •

$(\ell_{3},\ell_{2})=(2,2)$

[TABLE]

C.2 Adapted basis

Given the symmetries of the spacetimes that we consider, it seems more natural to write the metric on the $S^{3}$ as

[TABLE]

with $\chi\in[0,\pi/2]$ and $\phi,\psi\in[0,2\pi]$. Indeed, in this form, the rotation axis in the spacetime coincides with the $\phi$ direction in (67). Written in these coordinates, scalar harmonics on the $S^{3}$ are given by

[TABLE]

where $m=m_{+}+m_{-}$, $\bar{m}=m_{+}-m_{-}$, $P_{n}^{(a,b)}(x)$ is the Jacobi polynomial and

[TABLE]

with $|m_{\pm}|\leq\frac{\ell_{3}}{2}$ and $\frac{\ell_{3}}{2}-m_{\pm}\in\mathbb{N}$.

Given the symmetries of spacetimes that we are considering, we are interested in harmonics with $\bar{m}=0$. We consider the obvious basis of vectors on the $S^{3}$ given the form of the metric in (67),

[TABLE]

In this basis, the scalar-derived tensor harmonics on the $S^{3}$ with $\ell_{3}=2$ are given by the following

• •

$(\ell_{3},m)=(2,0)$

[TABLE]

• •

$(\ell_{3},m)=(2,2)$

[TABLE]

The $(\ell_{3},m)=(2,-2)$ harmonics are given by the complex conjugate of the $(\ell_{3},m)=(2,2)$ harmonics.

As discussed, for example, in Ref. Gualtieri et al. (2008), for a fixed $\ell_{n}$ one can relate the tensor harmonics on given basis to those on a different basis by a simple linear transformation.141414From representation theory, it follows that harmonics with different $\ell_{n}$ do not mix. Likewise, scalar-derived tensor harmonics in one basis will only mix with the scalar-derived tensor harmonics in the other basis. For the case at hand, we can relate the scalar derived tensor harmonics in the basis (70) to the harmonics in the other basis (60) as follows

[TABLE]

where in this expression $(\chi,\phi)$ should be understood as functions of $(\theta_{3},\theta_{2})$ and $R_{a}^{\phantom{a}b}$ is a rotation matrix given by

[TABLE]

For the $\ell_{3}=2$ scalar-derived tensor harmonics, we find that the transformation matrix $A^{m}_{\phantom{m}\ell_{2}}$ is given by:

[TABLE]

where the rows correspond to the values of $m=0,2,-2$, respectively, and similarly the columns correspond to the values of $\ell_{2}=0,1,2$. Furthermore, as shown in Ref. Gualtieri et al. (2008) reviewed in Appendix E, the inverse of this matrix (75), determines the transformation of the coefficients in the multipolar expansion of the Weyl scalars in different bases.

Appendix D Tensor harmonics on the $S^{4}$

In this Appendix we collect some properties of tensor harmonics in $S^{4}$. Furthermore, we list the harmonics that we have used to produce the waveforms in the main text.

D.1 Standard basis

We first write the metric on the $S^{4}$ in the standard form:

[TABLE]

where $ds^{2}_{S^{3}}$ is the metric on the $S^{3}$ written as in (59) and $\theta_{4}\in[0,\pi]$. Furthermore, we choose the obvious basis of unit vectors on the $S^{4}$:

[TABLE]

The advantage of using this form of the metric on the $S^{4}$ is that we can easily construct tensor harmonics of any rank following the algorithm of Ref. Higuchi (1987). To do so, we start from scalar harmonics on the $S^{4}$. Since we are interested in spacetimes that possess an SO$(3)$ in six dimensions, we can restrict ourselves to scalar harmonics that are constant on the $S^{2}$ sitting inside the $S^{3}$, that in turn sits inside the $S^{4}$, (76). These (normalized) scalar harmonics are given by:

[TABLE]

where $P_{\ell}^{m}(x)$ are the associated Legendre polynomials,

[TABLE]

is the normalization constant, and $\ell_{4}\geq\ell_{3}\geq 0$. In this basis, the projected components of the tensor harmonics with $\ell_{2}=\ell_{1}=0$ are as follows151515Note that for this particular class of harmonics, $\mathbb{Y}^{(\ell_{4},\ell_{3})}_{44}=\mathbb{Y}^{(\ell_{4},\ell_{3})}_{33}$, where here $\mathbb{Y}$ denotes any tensor harmonic.

Scalar-derived tensor harmonics

• •

$(\ell_{4},\ell_{3})=(2,0)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(2,1)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(2,2)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,0)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,1)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,2)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,3)$:

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,4)$:

[TABLE]

Vector-derived tensor harmonics

There are three families of vector-derived tensor harmonics on the $S^{4}$, but with our symmetry assumptions, only one of them has a nonzero overlap with the Weyl scalars. The relevant vector harmonics for us are given by the following

• •

$(\ell_{4},\ell_{3})=(2,1)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(2,2)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,1)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,2)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,3)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,4)$

[TABLE]

Transverse traceless tensor harmonics

Given the symmetries of the spacetimes that we are considering, there is only one family of transverse traceless tensor harmonics on the $S^{4}$ that can contribute to the multipolar expansion of the Weyl scalars. The relevant tensor harmonics are given by:

• •

$(\ell_{4},\ell_{3})=(2,2)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,2)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,3)$

[TABLE]

• •

$(\ell_{4},\ell_{3})=(4,4)$

[TABLE]

D.2 Multipolar expansion of the Weyl scalars in the standard basis

The main advantage of using the standard form of the metric on the $S^{4}$ (76) and the associated basis of vectors to obtain the multipolar expansion of the Weyl tensor, is that one can systematically construct the required tensor harmonics of any rank. In this way, we can identify the sector of tensor harmonics that captures most of the signal.

In Fig. 11 we display the $(\ell_{4},\ell_{3})=(2,2)$ multipole of the Weyl scalars in the scalar derived, vector derived and transverse traceless tensor harmonics sectors for the same $a/\mu^{1/3}=1.3$ simulation reported in the main text, Fig. 2. As this figure illustrates, most the signal is in the scalar derived sector. We have checked that this is the case in all our simulations and this is why in the main text we only report on the multipoles from the scalar derived tensor harmonics. Perhaps it should not be surprising that, given the symmetries of the spacetimes that we are considering, most the waveforms are captured by the scalar derived tensor harmonics.

The drawback of using the basis (76) is that it is not aligned with the axis of rotation of the black hole spacetimes that we have considered. This has the consequence that the various harmonics with the same $\ell_{4}$ but different $\ell_{3}$ mix and one cannot accurately extract the frequencies and growth/decay rates of the various modes.161616This phenomenon was observed in Gualtieri et al. (2008) in 4$D$, for instance, when extracting gravitational waves from head on collisions of black holes and considering an extraction frame that is not suitably aligned with the collision axis. We would like to thank Ulrich Sperhake for pointing this out to us. This mixing of modes with the same $\ell_{3}$ and different $\ell_{2}$’s can be seen in Fig. 12. In this plot we display the $(\ell_{4},\ell_{3})=(2,0)$, $(2,1)$ and $(2,2)$ multipoles of the Weyl scalars for the same $a/\mu^{1/3}=1.3$ run as in the main text. As this plot suggests, the various modes appear to have similar frequencies and decay rates. To disentangle the various modes one would have to rotate the basis or consider an adapted basis. We chose the latter option.

D.3 Adapted basis

As in the 5$D$ case, it turns out to be convenient to write the metric on the $S^{4}$ in a form that makes the symmetries and the rotation plane of spacetimes that we consider manifest:

[TABLE]

with $\chi\in[0,\pi/2]$, $\theta\in[0,\pi]$ and $\phi,\psi\in[0,2\pi]$. In these coordinates, the angle $\phi$ on the extraction $S^{4}$ coincides with the angle on the rotation plane of the full spacetime. We choose the obvious basis of angular vectors on the $S^{4}$ written as in (98):

[TABLE]

In the previous subsection §D.2, we have shown that the leading waveforms are in the sector of scalar derived tensor harmonics. Therefore, to accurately extract the leading modes that govern the dynamics of the black holes that we are interested in, we only need to concentrate on this sector of tensor harmonics. Scalar harmonics on the $S^{4}$ in the coordinates (98) can written as Berti et al. (2006)

[TABLE]

where $N$ is a normalization constant, $\mathbb{Y}^{\bar{m}}_{\ell_{2}}(\theta,\psi)$ are the standard spherical harmonics on the $S^{2}$, ${}_{2}F_{1}$ is the ordinary hypergeometric function and $k$ is a positive integer related to the eigenvalue $\ell_{4}$ by

[TABLE]

Since we are considering spacetimes with a manifest SO$(3)$ symmetry, we only need to consider harmonics with $\ell_{2}=0$ (and consequently $\bar{m}=0)$ and hence we will drop the corresponding labels from now on.

Writing the $S^{4}$ as in (98) and in the basis (99), the projected scalar derived tensor harmonics obtained from (100) are:

• •

$(\ell_{4},m)=(2,0)$:

[TABLE]

• •

$(\ell_{4},m)=(2,2)$:

[TABLE]

• •

$(\ell_{4},m)=(4,0)$:

[TABLE]

• •

$(\ell_{4},m)=(4,2)$:

[TABLE]

• •

$(\ell_{4},m)=(4,4)$:

[TABLE]

The $(\ell_{4},m)=(2,-2)$, $(4,-2)$, $(4,-4)$ harmonics are given by the complex conjugate of the $(2,2)$, $(4,2)$, $(4,4)$ respectively.

As we have already seen in Appendix C.2, for a fixed $\ell_{4}$ we can relate the harmonics in this basis (99) to those in the standard basis (77) by a linear transformation,

[TABLE]

where in the lhs of this expression $(\chi,\phi)$ should be understood as functions of $(\theta_{4},\theta_{3})$, and $R_{a}^{\phantom{a}b}$ is a rotation matrix given by

[TABLE]

For the $\ell_{4}=2$ scalar-derived tensor harmonics, the transformation matrix is given by

[TABLE]

where the rows correspond to $m=0,2,-2$ and the columns correspond to $\ell_{3}=0,1,2$. Similarly, for the $\ell_{4}=4$ harmonics, we find

[TABLE]

where the rows correspond to $m=0,2,-2,4,-4$ and the columns correspond to $\ell_{3}=0,1,2,3,4$ respectively.

Appendix E Transformation of the Weyl multipoles under changes of basis

In this Appendix, we review the transformations of the Weyl multipoles under rotations of the basis vectors and relate them to the transformation properties of the tensor harmonics. This is a straightforward generalization of the results in Ref. Gualtieri et al. (2008) to the higher dimensions.

The expansion of the Weyl scalars in multipoles is given by

[TABLE]

where $\ell$ is a collective label that specifies each of the tensor harmonics on the $S^{n}$ at infinity, $\theta$ denotes the collection of angles on such a sphere and $\psi_{l}$ are the corresponding multipoles. Under a rotation, $\Omega_{AB}^{\prime}(\theta)$ transforms as

[TABLE]

where the $R_{A}^{\phantom{A}B}$ are rotation matrices. Given the transformation of the harmonics under rotations,

[TABLE]

it follows that

[TABLE]

This transformation rule allows us to compute the Weyl multipoles in the (99) from the multipoles computed in the standard basis (77).

Appendix F Apparent horizon and contours of $\chi$

In cases in which the AH is a star-shaped surface, we find that contours of the conformal factor $\chi$ closely track the AH obtained by our apparent horizon finder. Figure 13 demonstrates this for a 6D MP black hole of dimensionless spin $a/\mu^{1/3}=1.3$ an initial $m=2$ deformation in $\chi$ described by (4) and (5). For this simulation, the AH can be tracked by our apparent horizon finder for the entire evolution. The AH at early times is roughly followed by the $\chi\sim 0.4$ contour, and as the gauge evolution proceeds, the AH at late times is most closely followed by the $\chi\sim 0.6$ contour. In particular, at time $t/\mu^{1/3}\sim 30$ when the AH is most elongated, the AH remains described well by the $\chi\sim 0.6$ contour.

Following $\chi$ contours is particularly useful when the AH ceases to be star shaped, and can no longer be tracked by our apparent horizon finder (though, see Refs. Figueras et al. (2017); Pook-Kolb et al. (2019) for different apparent horizon finder implementations that get around the star-shaped requirement). Figures 14 and 15 show the $\chi$ contours for two such cases of 6D MP black holes of dimensionless spin $a/\mu^{1/(D-3)}=1.5$ with initial $m=2$ and $m=4$ deformations in $\chi$, respectively, along with their AH shapes for as long as they can be tracked.

Figure 14 corresponds to an initial $m=2$ deformation, and Fig. 15 corresponds to an initial $m=4$ deformation. In both cases, the AH at early times is roughly followed by the $\chi\sim 0.4$ contour, and the gauge appears to evolve in such a way that the AH at late times is described well by the $\chi\sim 0.5$ contour. Thus, we continue to use the $\chi\sim 0.5$ contour as a proxy for the AH shape even at late times when the AH can no longer be found. With this proxy, the MP BH with an initial $m=2$ deformation develops a thin bar shape (see the last panel of Fig. 14), and the MP BH with an initial $m=4$ deformation develops four thin arms (see the last panel of Fig. 15).

Appendix G $m=3$ perturbation

We display here for completeness the evolution of the $6D$ MP BH with initial spin $a/\mu^{1/3}=1.3$ perturbed with a $m=3$ deformation, see eq. (5). Fig. 16 shows the evolution of the $\chi=0.5$ contour, which has qualitatively the same features as the $m=4$ perturbation discussed sec. III.2.3. The initial deformation of the black hole horizon into a triangle shape is characteristic of the $m=3$ mode. The corners of the triangle grow into three arms, which become elongated and develop sharp features at the edges. As in the $m=4$ case, the appearance of these sharp features precedes the formation of long and thin arms, which eventually become GL unstable.

Appendix H Convergence Test

Here we present a numerical test taken from a simulation of a 6D MP BH perturbed at $t/\mu^{1/3}=10$ by an $m=2$ deformation in $\chi$ described by Eqs. (4) and (5). Figure 17 shows the resulting gravitational waves via the $h_{+}$ component of the metric perturbation given by (10), extracted on the $z$ axis at $z/\mu^{1/3}=29$. This is done at three different resolutions, for which each linear dimension is covered by $80$, $120$, and $160$ points, respectively. The pointwise difference in the waveform between subsequent resolutions yields a convergence factor of approximately $3$, which indicates a rate convergence that is between second and third order.