End point of nonaxisymmetric black hole instabilities in higher dimensions

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.


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 [1] together with the image of the shadow of the black hole at the center of the M87 galaxy by the Event Horizon Telescope [2], 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 [3]. 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 [4] 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 [5], 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 [6,7] are asymptotically flat rotating black holes with horizon topology S 1 × 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. [8]. Reference [9] used numerical general relativity to simulate the nonlinear evolution of black ring instabilities and showed that thin enough black rings evolve into * Until August 2018 naked singularities in finite time, thus potentially violating the weak cosmic censorship in higher-dimensional asymptotically flat spaces. 1 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 [10] 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 socalled black folds [11], 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. [12] for linearized axisymmetric perturbations on top of a MP BH background. 2 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 [14]. 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. [15,16]. For a MP BH in D dimensions with mass parameter µ and spin parameter a, these nonaxisymmetric instabilities set in at smaller dimensionless spin a/µ 1/(D−3) than the axisymmetric instabilities. 1 There 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. 2 MP BHs with equal spins on all rotation planes are also known to be unstable [13]. 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.
Rotating black holes become dynamically unstable for values of a/µ 1/(D−3) ∼ O(1). Indeed, from thermodynamic considerations of the case in which a MP BH fragments and expels two nonspinning BHs, Emparan and Myers [10] estimated that nonaxisymmetric instabilities should set in at around a/µ 1/(D−3) ≈ 1. References [15,16] found that MP BHs with a/µ 1/(D− 3) 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φ for m = 2, i.e., a bar-mode deformation. In six and seven dimensions [16], 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 [17] corroborates these results, except in the five-dimensional (5D) case; while Ref. [15] found an instability, Ref. [17] 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/µ 1/(D−3) modes with m > 2 become the dominant unstable modes. This is indeed the case for black rings [9], 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 ≤ a/µ 1/(D−3) ≤ 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 [16,18,19] 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. [20] 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 χ with the apparent horizon, in Appendix G we display contours of χ 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 A we describe our evolution scheme together with the our choice of initial data and gauge evolution equations. In Sec. II B we provide some details about our approach to extracting gravitational waves in higher dimensions.

A. Evolution scheme
The results presented here are obtained by solving the Einstein field equations in the conformal and covariant Z4 (CCZ4) formulation [21,22]. 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 [16,18,19]. We redefine the constraint damping parameter κ 1 → κ 1 /α as in Ref. [23], in which α is the lapse function, and we typically use constraint damping values κ 1 = 0.4, κ 2 = 0. As in Ref. [14], we use initial data for a singly spinning MP BH, where Σ = r 2 + a 2 cos 2 θ, ∆ = r 2 + a 2 − µ r 5−D , and µ 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 ρ, where r h is the largest real root of ∆(r h ) = 0. Then, the Cartesian coordinates are given by x = ρ sin θ cos φ , y = ρ sin θ sin φ , z = ρ cos θ cos ϕ 1 , Imposing the SO(D − 3) symmetry corresponds to working on the slice w 1 = . . . = w D−4 = 0; see Refs. [16,18,19].
For values of a/µ 1/3 ∼ 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 χ: where χ 0 is the conformal factor computed from the analytic initial data (1), χ 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: Introducing an m = 2 perturbation immediately at t/µ 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/µ 1 D−3 = 0. For some simulations, we also introduce this deformation in χ after the gauge adjustment period at around t/µ 1 D−3 ∼ 10. For higher spins, this amounts to choosing a different branch of the dynamical evolution. Perturbing after the gauge adjustment period at around t/µ 1 D−3 ∼ 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 B 3.
This method of perturbing the black hole violates the constraints, so even though the CCZ4 formulation quickly 3 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 [21], where c α is a freely adjustable coefficient. In our runs we found that c α = 1.5 worked well for D = 5, 6, 7. 4 To evolve the shift, we used the modified Gamma-driver condition introduced in Ref. [9], whereΓ i is the usual CCZ4 evolution variable,Γ i is the contracted Christoffel symbol computed from the (conformally 3 In practice this timescale is exponential and much faster than the physical timescale of the instabilities. 4 The preferred choice in the typical four-dimensional astrophysical setting is cα = 2. rescaled) initial data (1), and and where δ 1,2 are two adjustable parameters that control the decay of the initial gamma function near the horizon. Since the initial data haveΓ i = 0, the role of the extra term in (7) is to drive the gauge towardsΓ 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 δ 1 = 0.2 and δ 2 = 0.075. We can freely adjust the coefficient c β in (8); c β = 0.6 works well for us (note that the typical value in four-dimensional simulations is c β = 0.75). Following Ref. [9], we introduce diffusion terms (well inside the apparent horizon) on the righthand 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. [9] for more details. We use Kreiss-Oliger dissipation [24] to damp unphysical high-frequency modes that can arise at grid boundaries or during regridding, with typical dissipation values of σ = 0.4. We numerically solve the CCZ4 equations using the GRCHOMBO libraries [25,26], 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 ≤ z ≤ 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 fourthorder finite differences and integrate in time using 4th-order Runge-Kutta method (RK4). We obtain approximately thirdorder convergence, see Appendix H for convergence tests.
We also track the apparent horizon (AH) during several stages of the evolution, (see Ref. [9]); 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. [14] or more recently Ref. [27]. We leave this for future work. However, as in Ref. [9], we verify that in our working gauge certain contours of χ track the AH to within a few percent. Thus, we can use χ 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 χ.

B. 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 × components of the metric perturbations at a certain radius along the z axis [15,16], As noted in these references, h + and h × 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 [16,17]. 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. [17] 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 higherdimensional analogues of the Newman-Penrose scalars introduced in Ref. [20]. We follow the implementation of Ref. [28] to calculate the nonvanishing components of the Weyl tensor along the outgoing null rays, Ω (A)(B) . For the class of spacetimes that we consider, the nonvanishing components of Ω (A)(B) are [28] where the indices ı,  = 2, 3 run along the spatial angular dimensions of the computational domain and the indicesâ,b, . . . run along the transverse sphere. Here R 0klm , etc.; denote the components of the spacetime Riemann tensor, m i (1) is the unit radial vector and m k (î) are the angular unit vectors. To construct an orthonormal basis of vectors, we start from and use the standard Gram-Schmidt orthornormalization procedure. Changing to spherical coordinates as in (B2), one can show that as r → ∞, the basis orthonormalized vectors approach The unit vectors along the transverse sphere directions are simply given by, We extract the individual modes in the gravitational wave signal by projecting Ω (A)(B) onto a basis of tensor harmonics on the sphere at infinity: where . . . denotes the set of quantum numbers that identify each of the harmonics, and dΩ (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. [29], 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. 5

Energy flux
Following Ref. [28] (and references therein), the energy emitted in form of gravitational waves is given bẏ with 5 To 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.
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 aṡ 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 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 A, 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 B to study the nonlinear stability properties of 6D MP BHs of various spins: in Sec. III B 1, we study the gravitational waveforms and extract the quasinormal modes (QNMs), in Sec. III B 2 we study the physics of the AH, and in Sec. III B 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 (20) with A the horizon's area and C its equatorial circumference. In the dynamical regime, Eq. (20) defines r h , µ h , and a h . Note that our results scale with the chosen length scale µ, related to the black hole mass of the initial configuration. As the evolution proceeds, we dynamically rescale the physical quantities according to the length µ h . For instance, nonlinear dynamics that settles to a final MP BH in six dimensions results in a final, dimensionless spin of (a/µ 1/3 ) final = a h /µ 1/3 h . Before we proceed, we remind the reader of some wellknown properties of singly spinning MP BHs in D ≥ 5 spacetime dimensions. See Ref. [30] 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 ≥ 6, the angular momentum is unbounded from above, and for any finite angular momentum, the black hole horizon is nondegenerate and has finite area.
In Ref. [15], 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/µ 1/2 ∼ 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/µ 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. [17] studied perturbations of singly spinning MP BHs in D ≥ 5 with m ≥ 2. In the D = 5 case, the results of Ref. [17] 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. [15] and [17] would be compatible if the instability were nonlinear.
Here, we revisit the results from Ref. [15] 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/µ 1/2 for which Ref. [15] observed an instability. When comparing our results to those in Ref. [15], 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. [15], 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/µ 1/2 = 0.89, for instance, Fig. 1 shows the gravitational waveform extracted on the z axis [15] at z 0 /µ 1/2 = 21.5. The plot below shows the ( 3 , m) = (2, 2) scalar multipole of the gravitational waveform extracted from Ω (A)(B) . From the damped oscillatory decay, we obtain a dominant QNM with frequency ωµ 1/2 = 1.283 − 0.032 i, which agrees with the results from perturbation theory of [17], i.e, ω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 Ω (A)(B) , and they agree within the truncation error.
In five dimensions, we find that gravitational waveform is essentially in the ( 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 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. [17] 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/µ 1/2 that are arbitrarily close to extremality. As pointed out in Ref. [31], 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/µ 1/2 → 1. In this subsection we consider the nonlinear evolution of perturbed 6D MP BHs of different values of the spin parameter a/µ 1/3 . The results in seven dimensions are qualitatively similar, and we will not show them here.
Recall that in D ≥ 6 the dimensionless spin parameter a/µ 1/3 of singly spinning MP BHs can be arbitrarily large. Reference [10] noted that the thermodynamics of spinning black holes exhibits a qualitative change in behavior as the angular momentum increases beyond a certain (dimensiondependent) critical value of a crit /µ 1/(D−3) ; for values of a/µ 1/(D−3) smaller than this critical value, MP BHs behave similarly to the Kerr black hole, and for values of a/µ 1/(D−3) greater than this critical value, they behave like black membranes. In six dimensions, one finds a crit /µ 1/3 1.29, and the ultraspinning instability kicks in at a/µ 1/3 1.57 [32]. The end points of certain ultraspinning instability were worked out in Ref. [14]. On the other hand, the nonaxisymmetric instabilities of 6D MP BHs kick in at a/µ 1/3 0.74 [16,17]. In this article we study the endpoint of nonaxisymmetric instabilities for MP BHs with 0.7 ≤ a/µ 1/3 ≤ 1.5. The results presented here, taken in conjunction with those in Ref. [14], constitute a complete picture of the dynamics of singly spinning MP BHs.

Gravitational waves
In Table I, 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 [17]. In the left column of Table I 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/µ 1/3 = 0.7 are stable, while for a/µ 1/3 ≥ 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/µ 1/3 that we have considered.
In Fig. 2 we display the gravitational waveforms corresponding to black holes with initial spin a/µ 1/3 = 1.3 (top) and a/µ 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 4 = 2 (top) and 4 = 4 (bottom) modes in the scalar-derived sector of tensor harmonics for a/µ 1/3 = 1.5. Most of the signal is in the ( 4 , m) = (2, 2) mode, but the 4 = 4 also displays a sizeable signal. How-

Initial state
Final state spin loss = 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, 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". ever, for MP BHs with a/µ 1/3 = 1.3 the ( 4 , m) = (2, 2) mode governs the nonlinear evolution, with qualitatively similar results for a/µ 1/3 = 1.5. Table II displays the quasinormal frequencies in the 4 = 4 sector for a/µ 1/3 = 1.3 and a/µ 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/µ 1/3 for the initial black hole. For small enough perturbations and a small enough initial value of a/µ 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/µ 1/3 ) 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.

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/µ 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/µ 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/µ 1/3 , so it is conceivable that one can achieve even higher emission percentages for larger initial values a/µ 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/µ 1/3 or a sufficiently large initial deformation is no longer another MP BH but is a naked singularity.

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 [9] or in the ultraspinning instability of MP black holes  , 2) and ( 4, m) = (4, 4) quasinormal mode for m = 2 perturbation of a MP BH with the initial configuration a/µ 1/3 = 1.3 and a/µ 1/3 = 1.5. We estimate an error of ±0.04 in the real and imaginary parts of the frequencies (see footnote5).
Gravitational wave (m = 2 mode) at z0/µ 1/3 = 21.5, together with snapshots of the apparent horizon for initial spin parameter a/µ 1/3 = 1.3 (top) and a/µ 1/3 = 1.5 (bottom). The dashed and dotted lines provide the fit for the QNMs in the growing and decaying phases. Table I displays the frequencies normalized with respect to the initial/final horizon scales. For larger spins, the system stays longer in the transient regime between the initially growing phase and the final ring down. [14]. The latter reference managed to find non-star shaped AH using an intrinsic parametrization of the surface, while Ref. [27] 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. [9,14], we have verified that certain contours of the conformal factor χ provide a very accurate description of Gravitational waveforms for initial spin parameter a/µ 1/3 = 1.3 obtained from the higher-dimensional Weyl scalars projected on a sphere at finite radius: ( 4, m) = (2, 0) and the real part of the ( 4, m) = (2, 2) modes (top) and ( 4, m) = (4, 0) and real parts of the ( 4, m) = (4, 2), (4, 4) modes (bottom). Most of the signal lies on the ( 4, m) = (2, 2) mode, but higher harmonics are also excited due to the nonlinearities. the AH, within less than a few percent. See Appendix F for quantitative comparison between the contours of χ and the actual AH for various runs in which the latter can be found. Therefore, we will use contours of χ 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 a/µ 1/3 = 1.5 a/µ 1/3 = 1.3 a/µ 1/3 = 1.1 a/µ 1/3 = 0.9 Horizon spin parameter a h /µ end points. 6 6 Note that since contours of χ provide accurate approximations of the AH, in principle one could use them as suitable reference surfaces to find the AH using the method of [27].  Fig. 2. This suggests that at higher spins gravitational wave emission becomes less efficient and instabilities have time to develop and eventually form a naked singularity. This run crashes close to the singularity but we artificially froze the evolution around the black hole to extract the gravitational waves. Therefore, the part of the waveform after (t − r0)/µ 1/2 ∼ 23 is unphysical.
For MP BH's with an initial a/µ 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. 7 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/µ 1/3 = 1.5 in Fig. 6 with their counterpart waveforms for the m = 2 instability for a/µ 1/3 = 1.3 in Fig. 2, notice that the amplitude for the m = 2 case is much larger, even though the AH defor- FIG. 7. Snapshots of the χ = 0.5 contour during the evolution of an unstable MP BH with initial a/µ 1/3 = 1.5 and no initial perturbation. This contour tracks the AH very closely; see Appendix F. The evolution is dominated by an m = 4 present in the numerical noise. During the evolution, the black hole develops a square shape, and the tips of the square eventually grow into long arms. The latter eventually become GL unstable. The last snapshot shows the second generation bulge propagating along the arm, a characteristic feature of GL dynamics. mation 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/µ 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/µ 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 χ = 0.5 contours from the evolution of the MP BH with a/µ 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. 8 However, the local evolution of the instability along each of the arms appears to be analogous to the ultraspinning case of Ref. [14]. 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 8 Recall that the partial differential equations (PDE) that we are solving are effectively (3 + 1) dimensional, while in Ref. [5] or [14], the system is (2+1) dimensional, thus allowing one to get much closer to the singularity. In analogy with the m = 4 case, the appearance of pointy tips precedes the formation of long and thin arms, which eventually become GL unstable. We suspect that this is the most likely fate of this evolution.
instability of black strings [5], 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/µ 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/µ 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/µ 1/3 = 0 only results in a dumbbell 9 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/µ 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.

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 4 = 2 and 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/µ 1/3 = 1.3 and a/µ 1/3 = 1.5. The peak of emission occurs around (t − r 0 )/µ 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 B 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/µ 1/3 = 1. 3 (continuous red) and a/µ 1/3 = 1.5 (dotted blue). The emission reaches its peak during the initial growing phase, followed by a rather constant emission during the transient period and finally dropping to zero during the ring down -cf. Fig.2. Bottom: total radiated energy, corresponding to approximately 23.54% (a/µ 1/3 = 1.3) and 31.15% (a/µ 1/3 = 1.5) of the BH's mass. a/µ 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 [34][35][36] (in D = 6, the fraction of mass radiated is (8.19 ± 0.05) × 10 −4 in the head-on collisions and (1.99 ± 0.05) × 10 −1 in the binary merger), or in the first black hole binary merger recorded by LIGO [1], 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 with M irr the irreducible mass 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 rad /E avail ∼ 0.60 and E rad /E avail ∼ 0.64 for a/µ 1/3 = 1.3 and a/µ 1/3 = 1.5, respectively.
Finally, we compare the energy emission for the MP BH with spin a/µ 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. [15] and the perturbation theory results of [17]. In the 6D case we found that slowly spinning black holes are stable. As the initial spin increases, at around a/µ 1/3 ≈ 0.8 MP BHs become unstable under perturbations with an m = 2 azimuthal mode number (i.e., bar mode), as Ref. [16] had previously observed and Ref. [17] 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/µ 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/µ 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/µ 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/µ 1/3 ≈ 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 [9]. 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. [37]. 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 [9], 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/µ 1/3 = 1.572 in D = 6 [32]; the evolution and end points of the latter were studied in Ref. [14]. 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. [14], provides a complete picture of the dynamics of some of the most notable instabilities that afflict MP black holes and their respective end points. 10 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 higherdimensional 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 10 The superradiant instability is another important instability that affects these black holes, but we have not considered it here. 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/µ 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. [37] 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. [37] is realized in finite dimensions. D = 6 should suffice.

ACKNOWLEDGMENTS
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 Mid-lands+ 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 We assume the spacetime to be asymptotically flat in the sense of Ref. [20]. Thus, in the wave zone, the line element is expressed in the Bondi-Sachs form [38,39] ds 2 = −Ae B du 2 −2e B du dr +r 2 h αβ (dφ α +C α )(dφ β +C β ), (A1) with u a retarded time and, A, B, C α and h αβ functions on the Bondi coordinate (u, r, φ α ). Moreover, det h αβ = det ω αβ , where ω αβ is the unit metric on the unit n-sphere S n . Here, we are interested in the asymptotic expansion [40] h αβ (u, r, αβ the Bondi news function. We consider the higher-dimensional generalization of the Geroch-Held-Penrose formalism [41] and introduce a tetrad basis ( , k, m A ) satisfying µ k µ = 1, m µ (A) m (B)µ = δ (A)(B) , and with vanishing further contractions. Asymptotically, one particular choice for the tetrad is Gravitational waves are extracted from the projection of the Weyl tensor C µνρσ onto the tetrad basis (A3). Specifically, the news function h (1) αβ is given by the leading contribution of Ref. [20], From the symmetries of C µνρσ , one verifies that Ω 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!C and!D we specialize to the cases of interest for us, namely, S 3 and S 4 .

(B1)
Here, θ 1 ∈ [0, 2π] and θ a ∈ [0, π] ∀a = 2 . . . , n. In the three-dimensional computational domain, we parametrize the numerical Cartesian coordinates (x, y, z) by new spherical coordinates (r, θ, ϕ) via x = r cos(θ n ) , y = r sin(θ n ) cos(θ n−1 ) , z = r sin(θ n ) sin(θ n−1 ) , 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. [29] 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 with 0 ≤ θ 1 ≤ 2π and 0 ≤ θ i ≤ π, ∀i = 2, . . . , 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 , 12 ∇ a ∇ a S n ,..., 1 = − n ( n + n − 1) S n,..., 1 , where the integers i satisfy n ≥ n−1 ≥ . . . ≥ 2 ≥ | 1 |. The scalar harmonics are normalized so that δ n n ..
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 and they satisfy where g is the metric on the round unit S n . On the other hand, divergence-free vector harmonics V n ,..., 1 a on the S n satisfy (B12) 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 (B6), form a complete basis of vectors on the S n , see Refs. [42][43][44]. From the scalar harmonics (B4), one can obtain two types of symmetric tensor harmonics: S n,..., 1 These satisfy: On the other hand, we have ∇ c ∇ c S n,..., 1 (2)ab = [− n ( n + n − 1) + 2 n] S n ,..., 1 (2)ab (B19) ∇ a S n ,..., 1 g ab S n,..., 1 δ n n ... 1 1 = S n g ac g bd S n ,..., 1 (2)ab S n ,..., 1 (2)cd . (B22) Recall that the matrix of Weyl scalars, Ω AB , encoding the gravitational radiation is symmetric and traceless. Therefore, the family 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, 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: 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 n−2 = . . . = 1 = 0 for the quantum numbers labeling the harmonics. We shall restrict to this class of harmonics from now on. (C3) • ( 3 , 2 ) = (2, 1): (C4) 13 Since tensor harmonics are symmetric Y ... 21 = Y ... 12 so we do not need to list all components.
b. 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 c. 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 6 sin 2 (θ 2 ) cos(2θ 3 ) + 3 cos(2θ 2 ) + 1 . (C8)

Adapted basis
Given the symmetries of the spacetimes that we consider, it seems more natural to write the metric on the S 3 as ds 2 = dχ 2 + sin 2 χ dφ 2 + cos 2 χ dψ 2 , with χ ∈ [0, π/2] and φ, ψ ∈ [0, 2π]. Indeed, in this form, the rotation axis in the spacetime coincides with the φ direction in (C9). Written in these coordinates, scalar harmonics on the S 3 are given by Given the symmetries of spacetimes that we are considering, we are interested in harmonics withm = 0. We consider the obvious basis of vectors on the S 3 given the form of the metric in (C9), In this basis, the scalar-derived tensor harmonics on the S 3 with 3 = 2 are given by the following (C13) The ( 3 , m) = (2, −2) harmonics are given by the complex conjugate of the ( 3 , m) = (2, 2) harmonics.
As discussed, for example, in Ref. [45], for a fixed n one can relate the tensor harmonics on given basis to those on a different basis by a simple linear transformation. 14 For the case at hand, we can relate the scalar derived tensor harmonics in the basis (C12) to the harmonics in the other basis (C2) as follows where in this expression (χ, φ) should be understood as functions of (θ 3 , θ 2 ) and R b a is a rotation matrix given by (C16) For the 3 = 2 scalar-derived tensor harmonics, we find that the transformation matrix A m 2 is given by: where the rows correspond to the values of m = 0, 2, −2, respectively, and similarly the columns correspond to the values of 2 = 0, 1, 2. Furthermore, as shown in Ref. [45] reviewed in Appendix E, the inverse of this matrix (C17), determines the transformation of the coefficients in the multipolar expansion of the Weyl scalars in different bases.
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. [44]. 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 , (D1). These (normalized) scalar harmonics are given by: where P m (x) are the associated Legendre polynomials, is the normalization constant, and 4 ≥ 3 ≥ 0. In this basis, the projected components of the tensor harmonics with 2 = 1 = 0 are as follows 15 a. Scalar-derived tensor harmonics  (D7) 15 Note that for this particular class of harmonics, Y  The main advantage of using the standard form of the metric on the S 4 (D1) 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 ( 4 , 3 ) = (2, 2) multipole of the Weyl scalars in the scalar derived, vector derived and transverse traceless tensor harmonics sectors for the same a/µ 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 (D1) 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 4 but different 3 mix and one cannot accurately extract the frequencies and growth/decay rates of the various modes. 16  (D24) 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 (D23) can written as [46] Y 4,m, 2,m = N (sin χ) |m| (cos χ) 2 e imφ Ym 2 (θ, ψ) 2 F 1 2 + |m| − k, k + 3 2 , 2 + 3 2 ; cos 2 χ (D25) where N is a normalization constant, Ym 2 (θ, ψ) 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 4 by 4 = 2k − ( 2 + |m|) . (D26) Since we are considering spacetimes with a manifest SO(3) symmetry, we only need to consider harmonics with 2 = 0 (and consequentlym = 0) and hence we will drop the corresponding labels from now on. Writing the S 4 as in (D23) and in the basis (D24), the projected scalar derived tensor harmonics obtained from (D25) are: •  (D33) For the 4 = 2 scalar-derived tensor harmonics, the transformation matrix is given by where the rows correspond to m = 0, 2, −2 and the columns correspond to 3 = 0, 1, 2. Similarly, for the 4 = 4 harmonics, we find 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 χ ∼ 0.4 contour, and the gauge appears to evolve in such a way that the AH at late times is described well by the χ ∼ 0.5 contour. Thus, we continue to use the χ ∼ 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). as the m = 4 perturbation discussed sec. III B 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 appear-ance of these sharp features precedes the formation of long and thin arms, which eventually become GL unstable.
FIG. 16. Snapshots of the χ = 0.5 contour during the evolution of an unstable MP BH with initial a/µ 1/3 = 1.3 and a m = 3 perturbation. During the evolution, the black hole develops a triangle shape, and its tips eventually grow into long arms. The latter eventually become GL unstable, as discussed in Sec. III B 3

Appendix H: Convergence Test
Here we present a numerical test taken from a simulation of a 6D MP BH perturbed at t/µ 1/3 = 10 by an m = 2 deformation in χ 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/µ 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.