Collapse and Nonlinear Instability of AdS with Angular Momenta

TL;DR

This study numerically explores the nonlinear stability and collapse dynamics of rotating scalar fields in AdS$_5$, revealing how angular momentum influences collapse times and stability of boson stars, and identifying endpoints as black holes.

Contribution

It provides new insights into how angular momentum affects collapse times and stability of boson stars in AdS, including the delay of collapse and the nature of instability outcomes.

Findings

Angular momentum delays gravitational collapse in AdS.

Boson stars near AdS are stable, far from AdS are unstable.

Collapse endpoints include Myers-Perry and hairy black holes.

Abstract

We present a numerical study of rotational dynamics in AdS$_5$ with equal angular momenta in the presence of a complex doublet scalar field. We determine that the endpoint of gravitational collapse is a Myers-Perry black hole for high energies and a hairy black hole for low energies. We investigate the timescale for collapse at low energies $E$, keeping the angular momenta $J\propto E$ in AdS length units. We find that the inclusion of angular momenta delays the collapse time, but retains a $t\sim1/E$ scaling. We perturb and evolve rotating boson stars, and find that boson stars near AdS appear stable, but those sufficiently far from AdS are unstable. We find that the dynamics of the boson star instability depend on the perturbation, resulting either in collapse to a Myers-Perry black hole, or development towards a stable oscillating solution.

Tables

Table 1. Table 1: Comparison of quasinormal modes (QNM) extracted using Prony analysis from late-time evolution of Gaussian initial data with λ = 1 𝜆 1 \lambda=1 , ϵ = 0.5 italic-ϵ 0.5 \epsilon=0.5 ( J ≈ 0.0136 𝐽 0.0136 J\approx 0.0136 , E ≈ 0.0873 𝐸 0.0873 E\approx 0.0873 ) and from linear perturbation theory of the Myers-Perry black hole with the same E 𝐸 E and J 𝐽 J .

QNM	Prony Analysis	Perturbation Theory
$⟨B⟩$	$5.829-5.263\times {10}^{-3}i$	$5.828-5.303\times {10}^{-3}i$
$⟨{\Pi }_{\Re }⟩$	$4.7528-0.03130i$	$4.7523-0.03127i$
$⟨{\Pi }_{\Im }⟩$	$4.7522-0.03126i$	$4.7523-0.03127i$

Equations

$$S=\frac{1}{16\pi G_5}\int {\mathrm{d}}^{5}x\sqrt{-g}\left(R+12-2|\nabla \Pi {|}^2\right).$$

$$\mathrm{d}{s}^2$$

$$+\frac{4\alpha \beta }{a}\rho \mathrm{d}t\mathrm{d}\rho +\frac{\mathrm{d}{\rho }^2}{a(2-{\rho }^2)}+$$

$$\displaystyle\qquad+\rho^{2}(2-\rho^{2})\bigg{[}\frac{1}{b^{2}}\Big{(}\mathrm{d}\psi+\cos^{2}(\tfrac{\theta}{2})\mathrm{d}\phi-\Omega\mathrm{d}t\Big{)}^{2}+$$

$$\displaystyle\qquad\qquad\qquad\qquad+\frac{b}{4}\Big{(}\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2}\Big{)}\bigg{]}\bigg{\}}\;,$$

$$\Pi$$

$${\partial }_{\rho }\phi +A_{\phi }[\rho ,a]\phi =g_{\phi }[\rho ,a,q]$$

$$\left[\begin{array}{c}{\partial }_t\phi \\ {\partial }_{t}q\\ {\partial }_{t}p\end{array}\right]=\left[\begin{array}{ccc}0& 0& 0\\ \gamma A_{\gamma }& A_d& A_q\\ 0& A_p& A_d\end{array}\right]\left[\begin{array}{c}{\partial }_{\rho }\phi \\ {\partial }_{\rho }q\\ {\partial }_{\rho }p\end{array}\right]+\left[\begin{array}{c}{f}_{\phi }\\ f_q+\gamma f_{\gamma }\\ f_p\end{array}\right],$$

$$A_u[\rho ]{\partial }_{\rho }u+B_u[\rho ,\phi ,p]u$$

$$A_v[\rho ]{\partial }_{\rho }v+B_v[\rho ,\phi ,q,p,u]v$$

$${\partial }_{t}u=f_u,$$

$$\Pi =(1-{\rho }^2{)}^4⟨\Pi ⟩+\mathcal{O}[(1-{\rho }^2{)}^5].$$

$${\Pi }_{\mathfrak{I}}{|}_{t=0}$$

$$\displaystyle\frac{\partial_{t}\Pi_{\mathfrak{R}}}{\alpha\sqrt{a}}\bigg{|}_{t=0}$$

$$\mathrm{d}{s}^2$$

$$\displaystyle\qquad\qquad\qquad\quad+\rho^{2}(2-\rho^{2})\bigg{[}\frac{1}{b^{2}}\Big{(}\mathrm{d}\psi+\cos^{2}(\tfrac{\theta}{2})\mathrm{d}\phi-\Omega\mathrm{d}t\Big{)}^{2}+\frac{b}{4}\Big{(}\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2}\Big{)}\bigg{]}\bigg{\}}\;,$$

$$\Pi$$

$$b$$

$${\Pi }_{\mathfrak{R}}$$

$${\Pi }_{\mathfrak{I}}$$

$$\beta$$

$$a$$

$$\alpha$$

$$\Omega$$

$$P^2$$

$$W^2$$

$$[PW{]}^2$$

$$Q^2$$

$$\displaystyle\qquad\qquad+\Big{[}\rho(2-\rho^{2})q_{r}+(1-\rho^{2})^{2}\varphi_{r}\Big{]}^{2}+\Big{[}\rho(2-\rho^{2})q_{i}+(1-\rho^{2})^{2}\varphi_{i}\Big{]}^{2}\;,$$

$$[QP]$$

$$\displaystyle\qquad\qquad+\Big{[}\rho(2-\rho^{2})q_{r}+(1-\rho^{2})^{2}\varphi_{r}\Big{]}p_{r}+\Big{[}\rho(2-\rho^{2})q_{i}+(1-\rho^{2})^{2}\varphi_{i}\Big{]}p_{i}\;,$$

$${\Phi }^2$$

$$[Q\Phi {]}^2$$

$$A^2$$

$$[B_{\delta \alpha }]$$

$$S_a$$

$$S_{\delta }$$

$$\displaystyle\qquad\qquad\qquad\qquad-\bigg{(}2-\rho^{2}(2-\rho^{2})(1-\rho^{2})^{4}\Big{[}3+\rho^{2}(2-\rho^{2})\Big{]}u_{a}\bigg{)}[Q\Phi]^{2}\bigg{\}}+$$

$$\displaystyle\qquad+\rho^{2}(2-\rho^{2})\bigg{\{}1-(1-\rho^{2})^{2}\bigg{(}3-\rho^{2}(2-\rho^{2})\Big{[}1-2\rho^{2}(2-\rho^{2})\Big{]}\bigg{)}u_{a}\bigg{\}}u_{a}\;.$$

$$(1-{\rho }^2){\partial }_{\rho }{\phi }_b-8\rho \left[1+\frac{(1-{\rho }^2{)}^6{u}_a}{2\sqrt{a}(1+\sqrt{a})}\right]{\phi }_b$$

$$(1-{\rho }^2){\partial }_{\rho }{\phi }_r-8\rho \left[1+\frac{(1-{\rho }^2{)}^6{u}_r}{4\sqrt{a}(1+\sqrt{a})}\right]{\phi }_r$$

$$(1-{\rho }^2){\partial }_{\rho }{\phi }_i-8\rho \left[1+\frac{(1-{\rho }^2{)}^6{u}_i}{4\sqrt{a}(1+\sqrt{a})}\right]{\phi }_i$$

$${\partial }_t{\phi }_b$$

$${\partial }_t{\phi }_r$$

$${\partial }_t{\phi }_i$$

$${\partial }_t{q}_b$$

Full text

Introduction – Spacetimes with anti-de Sitter (AdS) boundary conditions play a central role in our understanding of gauge/gravity duality Maldacena (1999); Gubser et al. (1998); Witten (1998); Aharony et al. (2000), where solutions to the Einstein equation with a negative cosmological constant are dual to states of strongly coupled field theories. This correspondence has inspired the study of gravitational physics in AdS over the past two decades.

It is perhaps surprising that the issue of the nonlinear stability of (global) AdS was only raised nine years after AdS/CFT was first formulated Dafermos and Holzegel (2006); Dafermos (2006). Dafermos and Holzegel conjectured a nonlinear instability where the reflecting boundary of AdS allows for small but finite energy perturbations to grow and eventually collapse into a black hole. This is in stark contrast with Minkowksi and de Sitter spacetimes, where nonlinear stability has long been established Friedrich (1986); Christodoulou and Klainerman (1993).

The first numerical evidence in favour of such an instability of AdS was reported in Bizon and Rostworowski (2011). This topic has since attracted much attention both from numerical and formal perspectives Dias et al. (2012a, b); Buchel et al. (2012, 2013); Maliborski and Rostworowski (2013a); Bizon and Jamuna (2013); Maliborski (2012); Maliborski and Rostworowski (2013b); Baier et al. (2014); Jamuna (2013); Basu et al. (2013); Gannot (2012); Fodor et al. (2014); Friedrich (2014); Piotr (2014); Maliborski and Rostworowski (2014); Abajo-Arrastia et al. (2014); Balasubramanian et al. (2014); Bizon and Rostworowski (2015); Balasubramanian et al. (2015); da Silva et al. (2015); Craps et al. (2014); Basu et al. (2015a); Deppe et al. (2015); Dimitrakopoulos et al. (2015); Horowitz and Santos (2015); Buchel et al. (2015); Craps et al. (2015a); Basu et al. (2015b); Yang (2015); Fodor et al. (2015); Okawa et al. (2015); Bizon et al. (2015); Dimitrakopoulos and Yang (2015); Green et al. (2015); Deppe and Frey (2015); Craps et al. (2015b, c); Evnin and Krishnan (2015); Menon and Suneeta (2015); Jalmuzna et al. (2015); Evnin and Nivesvivat (2016); Freivogel and Yang (2015); Dias and Santos (2016); Evnin and Jai-akson (2016); Deppe (2016); Dimitrakopoulos et al. (2016a, b); Rostworowski (2016); Jalmuzna and Gundlach (2017); Rostworowski (2017); Martinon et al. (2017); Moschidis (2017a, b); Dias and Santos (2017). Remarkably, this instability has recently been proved for the spherically symmetric and pressureless Einstein-massless Vlasov system Moschidis (2017a, b).

The collapse timescale is dual to the thermalisation time in the field theory, and is important for characterising and understanding this instability. For energies $E$ much smaller than the AdS length $L=1$, early evolution is well-described by perturbation theory. However, irremovable resonances generically cause secular terms to grow, leading to a breakdown of perturbation theory at a time $t\sim 1/E$. Numerical evidence suggests that horizon formation occurs shortly thereafter, i.e. at this same timescale. It is not fully understood why collapse seems to occur at the shortest timescale allowed by perturbation theory, though see Freivogel and Yang (2015) for some recent progress.

However, all numerical studies have been restricted to zero angular momentum. Though perturbation theory breaks down at $t\sim 1/E$ for systems with rotation as well Dias et al. (2012a); Dias and Santos (2016, 2017), this only places a lower bound on the timescale for gravitational collapse. It therefore remains unclear whether rotational forces could balance the gravitational attraction and delay the collapse time.

The inclusion of angular momentum also enriches the phase diagram of solutions. In addition to the Myers-Perry (and Kerr) family of rotating black holes, there are “black resonators” Dias et al. (2015) which can be described as black holes with gravitational hair, and “geons” Dias et al. (2012a); Horowitz and Santos (2015); Martinon et al. (2017) which are horizonless gravitational configurations held together by their own self-gravity. The nonlinear dynamics of these solutions remain largely unexplored.

Due to the lack of symmetries, the inclusion of angular momenta poses a numerical challenge (though see Bantilan et al. (2017) for recent progress away from spherical symmetry). For instance, the dynamical problem for pure gravity in four dimensions requires a full 3+1 simulation. To reduce numerical cost (see Bizon et al. (2005); Piotr (2014); Jalmuzna et al. (2015); Jalmuzna and Gundlach (2017) for the use of a similar strategy), we will rely on the fact that in odd dimensions $(d\geq 5)$, black holes have an enhanced symmetry when all of their angular momenta are equal. This simplification alone is insufficient for our purposes since gravitons that carry angular momenta break these symmetries. We therefore introduce a complex scalar doublet $\Pi$ given by the action

[TABLE]

As we shall see, this theory admits an ansatz with which one can study gravitational collapse with angular momentum using a 1+1 numerical simulation.

Moreover, this ansatz has a phase diagram of stationary solutions that is similar to that of pure gravity Dias et al. (2011). In particular, this theory contains hairy black holes and boson stars, which are somewhat analogous to black resonators and geons, respectively. Consequently, in addition to gravitational collapse, we are also able to investigate the dynamics of hairy black holes and boson stars.

Again, in this context, hairy black holes are much like black resonators, only with scalar hair instead of gravitational hair. Both the hairy black holes and boson stars exist for energies and angular momenta where Myers-Perry black holes are super-extremal and singular. For these conserved quantities, the weak cosmic censorship conjecture therefore implies that the final state following gravitational collapse cannot be a Myers-Perry black hole. For evolution respecting the symmetries of our ansatz, we wish to test cosmic censorship by identifying the endpoint of collapse.

Boson stars are horizonless solutions with a stationary metric and harmonically oscillating scalar field Dias et al. (2011); Liebling and Palenzuela (2012); Buchel et al. (2013). They are important objects for the study of the AdS instability since, like geons (and oscillons Maliborski and Rostworowski (2013a); Fodor et al. (2015) for a real scalar field), they can be generated as nonlinear extensions of normal modes of AdS. Such solutions can avoid the resonance phenomenon that leads to perturbative breakdown. Indeed, simulations of some of these solutions indicate stability well past $t\sim 1/E$. Initial data near these solutions therefore lie within an “island of stability”. We wish to investigate whether this stability applies for rotating boson stars.

However, boson stars far from AdS (i.e., past a turning point in their phase diagram) are expected to be unstable. We aim to determine the endpoint of this instability.

Method – We describe our ansatz and equations of motion schematically here; a full account is given in the Appendix. We take our metric and scalar to be

[TABLE]

where $\alpha$, $\beta$, $a$, $\Omega$, $b$, $\Pi_{\mathfrak{R}}$, and $\Pi_{\mathfrak{I}}$, are real functions of $t$ and $\rho$ only. This ansatz has $SU(2)$ rather than $U(1)\times U(1)$ symmetry due to the equal rotation in the $\psi$ and $\psi+\phi$ angles in orthogonal planes. This symmetry is preserved by $\Pi$. Without the scalar field, one finds that horizonless solutions have $\Omega=0$ and hence do not rotate.

Gauge freedom is fixed by maximal slicing, where the trace of the extrinsic curvature $K=0$ Choptuik et al. (1999). Unlike the choice $\beta=0$ in other studies, this gauge allows for evolution beyond the formation of an apparent horizon.

To bring the equations to first-order form, we introduce $w$ and $\delta$ as functions related directly to $\partial_{\rho}\Omega$ and $\partial_{\rho}\alpha$, respectively. Let the vector $\varphi=\{b,\Pi_{\mathfrak{R}},\Pi_{\mathfrak{I}}\}$, and introduce the vectors $p$, $q$ related to $\partial_{t}\varphi$ and $\partial_{\rho}\varphi$, respectively. Finally, define $u=\{w,\beta,a\}$ and $v=\{\delta,\alpha,\Omega\}$.

The vectors $\varphi$, $p$, $q$, $u$, $v$ represent our 15 dynamical functions. Their equations of motion can be summarized as follows. First, the definition of $q$ takes the form

[TABLE]

where $A_{\varphi}$ is a matrix and $g_{\varphi}$ is a vector.

Second, we have the evolution equations

[TABLE]

where the $A$’s are matrices that do not depend on $\varphi$, $q$, or $p$, and the $f$’s are vector-valued nonlinear expressions that can depend on all the dynamical functions. The first row of the equation above includes the definition of $p$. We have incorporated (3) into these evolution equations as damping terms with $\gamma$ acting as a damping coefficient.

Next, we have the slicing equations for $u$ and $v$

[TABLE]

where the $A$’s and $B$’s are matrices and the $g$’s are vectors. Included within these equations are the definitions of $w$ and $\delta$. These are nonlinear systems due to the dependence of the $g$’s on $u$ and $v$. However, given $\varphi$, $q$, and $p$, one can find $u$ and $v$ by solving a sequence of linear problems, with the solution to one linear problem entering nonlinearly as the source term for the next. (See Appendix for details.) Obtaining this property in the equations partially motivated the metric ansatz (2).

Finally, we have the evolution equations for $u$

[TABLE]

where $f_{u}$ is a vector that depends on all the functions.

Initial data is supplied as a choice of $\varphi$ and $p$. The remaining functions can be obtained by solving the nonlinear systems (3) and (5) with Newton-Raphson iteration. We evolve the system with a fourth order Runge-Kutta method that steps $\varphi$, $p$, and $q$ through (4), and then obtains $u$ and $v$ by solving (5) as a sequence of linear problems. We compute expansion coefficients from the metric to determine if a horizon has formed. We also monitor (6) and (3) as a check on numerics.

At infinity ($\rho=1$), we fix the boundary metric to be that of global AdS and require $\Pi=0$. The energy $E$ and angular momentum $J$ are read from the metric at infinity, and are conserved by (6). The response of the scalar field $\langle\Pi\rangle$ is obtained by

[TABLE]

Prior to horizon formation, we require regularity at the origin $\rho=0$. After horizon formation, our numerical grid will be excised, and boundary conditions at the excision surface are supplied for $u$ through (6), and the value of $\delta$ is held fixed 111A choice of $\delta$ at the excision surface fixes residual gauge freedom.. No other boundary conditions are required at the excision surface.

We use a spectral element mesh with Legendre-Gauss-Lobatto nodes, and inter-element coupling handled by a discontinuous Galerkin method with Lax-Friedrichs flux. Adaptive mesh refinement for splitting elements and increasing/decreasing polynomial order is decided by monitoring the Legendre spectrum within each element. Linear systems are solved via sparse LU decomposition. For all data presented here, relative energy and momentum violation, and violation of (3) are within or below $10^{-8}$. See the Appendix for more on numerical checks.

Gaussian Data – Consider Gaussian initial data

[TABLE]

with the other functions within $\varphi$ and $p$ vanishing. This data is parametrised by $\epsilon$ and $\lambda$. At fixed $\lambda$ and small varying $\epsilon$, we have $E\propto\epsilon^{2}$ and $J\propto E$. At larger $\epsilon$, deviation from this scaling may occur. This is a natural choice of parameters since individual normal modes of AdS that carry angular momentum also obey $J\propto E$ at small $E$. We take two families of initial data: one with fixed $\lambda=0$ where $J=0$, and another with $\lambda=1$, where $J\approx 0.155E$.

Let us describe the stationary black holes that can serve as final states of gravitational collapse. These black holes must fall within the symmetry class of our ansatz and have the same conserved quantities as our initial data. For $\lambda=0$, we have $J=0$ so the only stationary black holes are Schwarzschild-Tangherlini solutions.

For $\lambda=1$, there are two competing families of regular black hole solutions. Myers-Perry black holes have the most entropy where they exist, i.e. for $E>E_{\mathrm{extr}}\approx 0.0691$. For all energies $E<E_{\mathrm{extr}}$, hairy black holes have the most entropy (by being the only existing solution). We wish to verify that gravitational collapse for the $\lambda=1$ family will eventually settle into one of these black holes in accord with their respective energy ranges. This can be viewed as a test of cosmic censorship, since Myers-Perry would be superextremal for $E<E_{\mathrm{extr}}$.

A hairy black hole can be distinguished from the Myers-Perry case by the presence of the scalar field. In Fig. 1, we show the evolution of the normed square of the scalar response $|\langle\Pi\rangle|^{2}$ for two cases. In the first case, $E\approx 0.0873>E_{\mathrm{extr}}$, so the Myers-Perry configuration is the preferred solution. Indeed, we find that the scalar field vanishes at late times.

In the second case, $E\approx 0.0560<E_{\mathrm{extr}}$, so the Myers-Perry black hole is superextremal, and we find that the scalar field approaches a constant non-zero value at late times, indicative of a hairy black hole. We have also matched this value to that of the expected final hairy black hole solution which was first obtained in Dias et al. (2011).

In both cases, we have also matched the final entropy and angular frequency to their respective final stationary solutions Dias et al. (2011). In the subextremal case, we have matched quasinormal modes as well. (See Appendix.)

Now we compare the timescale for horizon formation between the $\lambda=0$ ($J=0$) and $\lambda=1$ ($J\approx 0.155E$) families of initial data. In Fig. 2, we show a log-log plot of the collapse time versus the energy. We see that at fixed energy, the initial data with nonzero angular momentum takes a longer time to collapse. However, the collapse times for both sets of initial data exhibit a power-law that is consistent with a $t\sim 1/E$ scaling. We conclude that in this case, angular momentum increases the collapse time but does not affect the timescale.

Boson Star Data – Boson stars within our ansatz have been constructed in Dias et al. (2011), and can be found by setting the metric to be time-independent and the scalar field to have a harmonically oscillating complex phase. They can be parametrised by their harmonic frequency $\omega$. For small energies, $\omega$ is close to a normal mode frequency of AdS. We focus on the lowest frequency mode with $\omega=5$ near AdS. For small energies, these particular boson stars have angular momentum $J\simeq 0.2E$.

As one increases the energy of the boson star, $\omega$ decreases until a turning point is reached around $\omega\approx 4.35$, where $E$ and $J$ are both maximal. Boson stars that lie on the AdS side of this turning point are expected to be nonlinearly stable (at least up to $t\sim 1/E$), and are otherwise expected to be unstable 222Morse theory typically implies that solutions on one side of a turning point are unstable. For completeness, we demonstrate linear instability in the Appendix..

We perturb boson stars near the turning point with frequencies $\omega=4.3$ (in the ‘unstable’ branch) and $\omega=4.4$ (in the ‘stable’ branch) with a Gaussian profile similar to (8a). Their scalar response $|\langle\Pi\rangle|^{2}$ is shown in Fig. 3. Note that though the scalar field oscillates with frequency $\omega$, these oscillations are canceled out in $|\Pi|^{2}$, and consequently are not seen in Fig. 3 nor in the metric.

Indeed, the $\omega=4.3$ boson star is unstable, but the endpoint of its evolution depends on the perturbation. For one perturbation (top panel of Fig. 3), evolution proceeds rapidly towards gravitational collapse, and eventually settles to a Myers-Perry black hole. While a competing hairy black hole also exists, it has less entropy than Myers-Perry in this region of parameter space.

Perturbing the same boson star with the opposite sign yields drastically different results. As one can see from the middle panel of Fig. 3, large $O(1)$ deformations develop in $|\langle\Pi\rangle|^{2}$ (and the metric) that oscillate for long times. The frequency of these oscillations is much smaller than the boson star frequency $\omega$. The metric and scalar both oscillate, so the final state (assuming continued stability) might be characterised as an oscillon. Since the frequency $\omega$ is still present in the scalar, this solution is, in a sense, a multi-frequency oscillon.

In contrast to the above, the lower panel of Fig. 3 shows that the perturbed boson with $\omega=4.4$ remains stable at long times, with no large deviations from the initial data.

We have repeated this study for different perturbations and boson stars, and also for oscillons (where $\Pi_{\mathcal{I}}=0$, see also Fodor et al. (2015)). We find no qualitative difference to the above.

Conclusions – Our numerical results suggest that much of our understanding of the instability of AdS carries over to situations with angular momenta as well. In particular, for generic data, the timescale for gravitational collapse $t\sim 1/E$ is preserved in the presence of rotation. Additionally, like oscillons in spherical symmetry, there are solutions that are nonlinear extensions of normal modes of AdS that are stable past $t\sim 1/E$.

We have also found that, depending on the perturbation, unstable boson stars will either collapse or oscillate. A comparison can be made to situations in flat space where the endpoint of unstable solutions can also depend upon the perturbation (see, e.g. Lee and Pang (1989); Seidel and Suen (1990); Lai and Choptuik (2007); Liebling and Palenzuela (2012); Sanchis-Gual et al. (2017); Carracedo et al. (2017)). In flat space, energy and angular momentum can be carried away, and the non-collapsing evolution is well-approximated by a perturbation of a stable boson star, presumably settling towards a stable boson star at asymptotically late time. By analogy, we suspect that the oscillating solution we find in AdS can be described as a nonlinear extension of a perturbed stable boson star. In AdS, however, there is a reflecting boundary which may cause the oscillations to persist indefinitely.

Finally, let us comment on interesting regions of parameter space that we have not studied. There is a range of energies and angular momenta where hairy black holes have more entropy than Myers-Perry black holes. This happens to be where Myers-Perry black holes are unstable to superradiance. In fact, hairy black holes branch off from Myers-Perry configurations precisely at the onset of this instability, for particular perturbations Dias et al. (2011).

This region is therefore a natural place to study the rotational superradiant instability Hawking and Reall (2000); Cardoso and Dias (2004); Kunduri et al. (2006); Dias and Santos (2013); Cardoso et al. (2014); Dias et al. (2015); Niehoff et al. (2016); Green et al. (2016) for which little is known fully dynamically. However, typical growth rates for this instability are around $10^{-5}$ Dias et al. (2011), which requires a longer simulation than we can feasibly perform with our methods. Furthermore, our ansatz implies that such a study will necessarily be incomplete. High angular wavenumbers are expected to play an important role in this instability Cardoso et al. (2014); Dias et al. (2015); Niehoff et al. (2016); Green et al. (2016), but our ansatz is restricted to only the $m=1$ azimuthal wavenumbers.

Acknowledgements – We would like to thank Mihalis Dafermos, Gary Horowitz, Luis Lehner and Harvey Reall for reading an earlier version of the manuscript. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by CENTRA/IST. Some computations were performed at the cluster “Baltasar-Sete-Sois” and supported by the H2020 ERC Consolidator Grant “Matter and strong field gravity: New frontiers in Einstein’s theory” grant agreement no. MaGRaTh-646597. Some computations were performed on the COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility and funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. O.J.C.D. is supported by the STFC Ernest Rutherford grants ST/K005391/1 and ST/M004147/1. BW is supported by NSERC.

Appendix

Equations of Motion

Here, we describe the equations of motion in full. The ansatz, as we have presented it, is reproduced here:

[TABLE]

We first perform a number of function redefinitions as follows:

[TABLE]

With these redefinitions, the finiteness of the new functions will ensure that the boundary conditions are satisfied. Vacuum global AdS is recovered when all of the new functions vanish.

From here, we introduce a number of functions that will put the equations of motion into first-order form. These are $q_{b}$, $q_{r}$, $q_{i}$, $p_{b}$, $p_{r}$, $p_{i}$, $u_{w}$ and $v_{\delta}$. We will give the definitions of these functions when presenting the equations of motion.

For ease of presentation, let us also define a number of auxiliary expressions:

[TABLE]

Now we present the equations of motion. We have a total of 15 functions which can expressed as $\varphi=\{\varphi_{b},\varphi_{r},\varphi_{i}\}$, $q=\{q_{b},q_{r},q_{i}\}$, $p=\{p_{b},p_{r},p_{i}\}$, $u=\{u_{w},u_{\beta},u_{a}\}$, and $v=\{v_{\delta},v_{\alpha},v_{\Omega}\}$. The full set of equations of motion are all generated from a combination of the Einstein equation, the Klein-Gordon equation for the scalar doublet, the maximal slicing gauge condition $K=0$, and definitions of first-order functions. There are a total of 21 equations.

The first three equations are linear slicing equations for $\varphi$ that define the functions in $q$:

[TABLE]

Note that the coefficient of the derivative term vanishes only at $\rho=1$. There is therefore a natural direction of integration for these equations, which is from the boundary $\rho=1$ towards the origin $\rho=0$. No external boundary conditions are required for these equations.

The next nine equations are evolution equations for $\varphi$, $q$, and $p$. The evolution equations for $\varphi$ are given by the definition of $p$:

[TABLE]

Note that finiteness of $\varphi$ requires that $p=0$ at $\rho=1$. We must enforce this condition in our numerical evolution.

The evolution equations for $q$ come from the commutation of time and spatial derivatives:

[TABLE]

where we have included damping terms with a non-negative constant $\gamma$ acting as a damping coefficient. Note that the $\gamma$ terms vanish when (A.4) is satisfied. One can show that deviations away from (A.4) will decay exponentially in time with a decay rate proportional to $\gamma$. As in the evolution equations for $\varphi$ (A.5), these equations require $p=0$ at $\rho=1$.

The evolution equations for $p$ are

[TABLE]

Finiteness of $p$ requires that $q=0$ at $\rho=0$. This condition is enforced during numerical evolution.

This collection of nine evolution equations (A.5), (A.6), and (A.7) take the form of an advection system

[TABLE]

for some vector of functions $\mathbf{u}$, advection operator $\mathbf{A}$, and nonlinear terms $\mathbf{f}$. The advection operator $\mathbf{A}$ has eigenvalues

[TABLE]

each with a three-fold degeneracy. These eigenvalues define the characteristics of this system. The two signs give velocities to the ingoing and outgoing characteristics. A computation of expansion coefficients indicates that an apparent horizon forms when

[TABLE]

which also corresponds to a situation where the outgoing characteristics switch sign. These eigenvalues are also used in our implementation of spectral element methods to determine flux across elements.

The next six equations are slicing equations for $u$ and $v$:

[TABLE]

where the last two equations (A.11e) and (A.11f) define $u_{\delta}$ and $u_{w}$, respectively. Each of these equations has a derivative term with a coefficient that vanishes only at $\rho=0$ or $\rho=1$. These determine the direction of integration for these equations. Namely, the first four are integrated from the interior out, and the last two are integrated from the boundary in. Prior to horizon formation, the origin is part of the numerical grid, and no external boundary conditions are required. After horizon formation, the numerical grid interior to the horizon will be excised, and boundary conditions will need to be supplied for the first four equations.

This set of slicing equations (A.11) is a nonlinear system, but can be solved as a series of linear systems. Suppose $\varphi$, $q$, and $p$ are known at a particular time slice. Then (A.11a) and (A.11b) can each be solved independently as a linear differential equation, obtaining $u_{w}$ and $u_{\beta}$. These functions enter into the source term $S_{a}$, and $u_{a}$ can then be obtained by solving (A.11c). Then $u_{w}$, $u_{\beta}$, and $u_{a}$ are placed in $[B_{\delta\alpha}]$ and $S_{\delta}$, and the coupled linear system (A.11a) and (A.11b) can be solved for $v_{\delta}$ and $v_{\alpha}$. Finally, $u_{w}$ and $v_{\alpha}$ enter into (A.11f) whose solution yields $v_{\Omega}$.

Note that instead of evolving $\varphi$ through (A.5), we could instead solve the slicing equation for $\varphi$ (A.4) along with the slicing equations (A.11). But, it is not possible to solve both (A.4) and (A.11) as a series of linear equations. We wish to avoid using nonlinear solvers during evolution, so we have opted to solve just (A.11) and handle (A.4) using damping terms in evolution. Note, however, that to obtain initial data, we do directly solve (A.4) and (A.11) using Newton-Raphson iteration.

Finally, the remaining three equations are evolution equations for $u$:

[TABLE]

These are degenerate evolution equations. That is, they do not contain spatial derivatives (or equivalently, their advection operator has only zero eigenvalues). Since the energy $E\propto u_{a}(t,\rho=1)$ and angular momentum $J\propto u_{w}(t,\rho=1)$, the vanishing of the equations (A.12a) and (A.12c) at $\rho=1$ imply that $E$ and $J$ are conserved.

After excision, these equations supply boundary conditions for $u_{w}$, $u_{\beta}$, and $u_{a}$ in the first three slicing equations of (A.11). The remaining boundary condition is determined by keeping $u_{\delta}$ fixed. One can show that this last condition just fixes an integration constant in the maximal slicing gauge condition $K=0$.

Lastly, for completeness, we express our Gaussian initial data in terms of the variables in this section:

[TABLE]

with the remaining functions in $\varphi$ and $p$ vanishing. The remaining functions can be obtained by solving (A.4) and (A.11) using Newton-Raphson iteration.

Numerical Validation

In this section, we present a number of additional numerical checks. Recall that we evolve the system by stepping $\varphi$, $q$, and $p$ in time using (A.5), (A.6), and (A.7), then obtain $u$ and $v$ by solving the slicing equations (A.11). There are two sets of equations that we do not solve directly: the slicing equations for $\varphi$ given in (A.4), and the evolution equations for $u$ in (A.12). We can use these equations for numerical checks.

We verify (A.4) by solving it for $\varphi$ and comparing the result to what we obtained through our main evolution. We call the difference of these quantities $\Delta_{\mathrm{Evol}}\varphi$. Similarly, given the numerical functions at a time slice, we take a single time step of the left hand side of (A.12) using (A.12) directly, and compare that result to what we obtain using our main evolution code. We call the difference of these $\Delta_{\mathrm{Evol}}u$. Finally, we monitor the relative violation of energy and angular momentum which are quantities conserved though (A.12).

In Fig. A.1, we show these various measures of error as a function of time. The differences $\Delta_{\mathrm{Evol}}\varphi$ and $\Delta_{\mathrm{Evol}}u$ fluctuate erratically, so we only show the running maximum. All of these errors are within or below $10^{-8}$.

Next, we perform an independent residual test on our solution. That is, we evaluate the residuals of all equations of motion using finite differences and demonstrate convergence with increasing resolution. By using spectral interpolation, we replace our spectral element mesh with uniform meshes at various resolutions, and we compute spatial derivatives via fourth-order finite differences. Time derivatives are also computed using fourth-order finite differences, but we fix the time step.

We take a time step near collapse and compute the residuals of all 21 equations as described. The results are shown as a log-log plot in the left panel of Fig. A.2. We have divided the equations into slicing equations and evolution equations, and have taken the infinity norm of these sets of equations. We find that the residuals converge with a fourth-order power law, as expected.

Finally, we test convergence of our code on a perturbed unstable boson star that leads to gravitational collapse. We carry out this computation from start to horizon formation on a series of fixed meshes, each with 16 elements and varying numbers of (Legendre-Gauss-Lobatto) nodes within each element. At the formation of an apparent horizon, we compute residuals as before. The right panel of Fig. A.2 shows that the convergence is exponential, as expected of spectral methods.

Matching of Final States

In Fig. 1 of the main paper, we presented the evolution of the scalar response $|\langle\Pi\rangle|^{2}$ from two sets of collapsing initial data. The scalar response $|\langle\Pi\rangle|^{2}$ in one of these vanishes at late times, which suggests that the final state is a Myers-Perry black hole. In the other, $|\langle\Pi\rangle|^{2}$ approaches a constant at late times, suggesting that the final state is a hairy black hole Dias et al. (2011). In this section, we present further evidence of the final states by matching other quantities.

One of these quantities is a pressure $\langle B\rangle$ extracted from the boundary stress tensor Balasubramanian and Kraus (1999); de Haro et al. (2001) defined as

[TABLE]

where $\varphi_{b}$ is defined in (A.2a), and we have chosen a conformal frame where the boundary metric is $\mathbb{R}^{(t)}\times S^{3}$. For solutions within the symmetry class of our ansatz, the boundary stress tensor is fully specified by $\langle B\rangle$, $E$, and $J$.

Once a horizon forms, we can also compute the entropy $S_{\mathrm{H}}$ and angular frequency $\Omega_{\mathrm{H}}$ of the apparent horizon. Though these are not gauge-invariant quantities during time evolution, they must approach that of a stationary solution at late times.

The evolution of these quantities is shown in Fig. A.3, expressed in units of a stationary black hole solution Dias et al. (2011) with the same $E$ and $J$ as the initial data. We see that, at late times, these quantities approach that of the stationary black hole solution. The final value for $\langle B\rangle$ in the superextremal case is the least conclusive of these since the remaining oscillations are still large compared to the hairy black hole value $\langle B\rangle_{\mathrm{HBH}}\approx 3.314\times 10^{-5}$. The smallness of this value in this case makes matching difficult.

For the simulation that shares conserved quantities with a Myers-Perry black hole, we have also compared quasinormal modes. The late-time behaviour of these simulations should be well-approximated by linear perturbation theory of the final stationary state. Using Prony analysis, we can extract quasinormal modes from the time evolution of the scalar response $\langle\Pi_{\mathfrak{R}}\rangle$, $\langle\Pi_{\mathfrak{I}}\rangle$, as well as the pressure $\langle B\rangle$. These can be compared to quasinormal modes computed directly from linear analysis of a stationary solution, which we have performed for Myers-Perry black holes. The comparison is made in Table 1, from which we see that the agreement is quite good.

Linear Perturbations of Boson Stars

In this section, we demonstrate the linear instability of boson stars that lie beyond the turning point (see Dias et al. (2011) for details on this turning point, and the top panel of Fig. A.4 for a phase diagram). One can typically prove that the existence of a turning point in the phase diagram implies an instability on at least one side of that point. Studies of boson stars and related solutions in situations different from ours suggests that the side connected to AdS should be stable, and the other side unstable. Nevertheless, for completeness, we perform a linear stability analysis to confirm this expectation.

For this purpose, we have decided to work in a different gauge than the rest of this paper. We choose a Schwarzschild-like gauge with $\beta=0$ in the ansatz . Boson stars are found in this ansatz by setting the metric to be time independent, and the scalar field to have a harmonic time dependence $\Pi(t,\rho)=e^{i\omega t}\tilde{\Pi}(\rho)$. The coupled gravity-matter equations for the boson stars are solved by Newton-Raphson iteration using pseudospectral methods on a Chebyshev grid (see Dias et al. (2016) for a review on these methods). A seed is provided by perturbation theory about AdS.

The functions are perturbed by choosing the form

[TABLE]

where $f$ stands for any metric or scalar field function, $f_{0}$ is the function on the boson star background solution, and the $\delta f$’s are their perturbations. The equations of motion are expanded to linear order, yielding an eigenvalue problem with eigenvalue $\tilde{\lambda}$. For the form we have chosen, real $\tilde{\lambda}$ is indicative of an instability. There is a symmetry $f_{s}\leftrightarrow-f_{s}$, $\tilde{\lambda}\leftrightarrow-\tilde{\lambda}$, so without loss of generality, we take $\tilde{\lambda}$ to be positive. Each eigenvalue will have a multiplicity of two. This degeneracy arises because one can shift time by a phase $t\to t+\phi_{0}$ and preserve the form of the functions (A.15) above.

We again employ pseudospectral methods to solve the eigenvalue problem Dias et al. (2016). We first solve the resulting matrix (generalized) eigenvalue problem by QZ factorization. We find no unstable modes in the branch of boson stars connected to AdS. Past the turning point, we identify an unstable mode. We track this mode in parameter space using Newton-Raphson iteration. To eliminate the degeneracy, we demand (in addition to a normalisation condition) that one of the perturbation functions takes a certain value at a spatial point. Since the resulting matrix problem is overconstrained, we solve it by linear least squares. We have checked that the solutions remain the same with different numerical resolutions.

In Fig. A.4, we show the phase diagram of boson stars and the eigenvalue corresponding to the unstable mode. We see that the eigenvalue approaches zero near the turning point (around $J\approx 0.083$), suggesting that the turning point is indeed the onset of this instability.

Interestingly, this eigenvalue also decreases rapidly near the edge of our numerics around $J\approx 0.045$. Near this region of parameter space, the phase diagram has a number of further turning points (and plots of $E$ and $J$ versus the boson star frequency $\omega$ produce spirals). This opens the possibility that boson stars may stabilize again past another turning point. Similar behaviour has been conjectured in Figueras et al. (2012) in the context of non-uniform black strings. However, it is difficult to proceed further since the solution approaches a singularity.