Phase Diagram of Stability for Massive Scalars in Anti-de Sitter Spacetime

We present the phase diagram of stability of 5-dimensional anti-de Sitter spacetime against horizon formation in the gravitational collapse of a scalar field, treating the scalar field mass and width of initial data as free parameters. We find that the stable phase becomes larger and shifts to smaller widths as the field mass increases. In addition to classifying initial data as stable or unstable, we identify two other phases based on nonperturbative behavior. The metastable phase forms a horizon over longer time scales than suggested by the lowest order perturbation theory, and the irregular phase can exhibit non-monotonic and chaotic behavior in the horizon formation times. Our results include evidence for chaotic behavior even in the collapse of a massless scalar field.


I. INTRODUCTION
Through the anti-de Sitter spacetime (AdS)/conformal field theory (CFT) correspondence, string theory on AdS 5 × X 5 is dual to a large N conformal field theory in four spacetime dimensions (R × S 3 when considering global AdS 5 ). The simplest time-dependent system to study in this context is the gravitational dynamics of a real scalar field with spherical symmetry, corresponding to the time dependence of the expectation value of the zero mode of a single trace operator in the gauge theory. Starting with the pioneering work of [1][2][3][4], numerical studies have suggested that these dynamics may in fact be generically unstable toward formation of (asymptotically) AdS d+1 black holes even for arbitrarily small amplitudes. While perhaps surprising compared to intuition from gravitational collapse in asymptotically flat spacetimes, the dual picture of thermalization of small energies in a compact space is more expected. In terms of the scalar eigenmodes on a fixed AdS background, the instability is a cascade of energy to higher frequency modes and shorter length scales (weak turbulence), which eventually concentrates energy within its Schwarzschild radius. In a naive perturbation theory, this is evident through secular growth terms.
However, some initial scalar field profiles lead to quasi-periodic evolution (at least on the time scales accessible via numerical studies) at small but finite amplitudes; even early work [1,5] noted that it is possible to remove the secular growth terms in the evolution of a single perturbative eigenmode. A more sophisticated perturbation theory [6][7][8][9][10][11][12][13][14][15][16][17] supports a broader class of quasi-periodic solutions that can contain non-negligible contributions from many modes, and other stable solutions orbit the basic quasi-periodic solutions [14]. Stable solutions exhibit inverse cascades of energy from higher frequency to lower frequency modes due to conservation laws following from the high symmetry of AdS (integrability of the dual CFT). Stable behavior also appears in the full non-perturbative dynamics for initial profiles with widths near the AdS length scale [18][19][20]; however, analyses of the perturbative and full dynamics in the literature have not always been in agreement at fixed small amplitudes. For example, some perturbatively stable evolutions at finite amplitude actually form black holes in numerical evaluation of the full dynamics [6,21,22]. Understanding the breakdown of the approximations used in the perturbative theory, as well as its region of validity, is an active and important area of research [23][24][25][26][27].
Ultimately, the main goal of this line of inquiry is to determine whether stability or instability to black hole formation (or both) is generic on the space of initial data, so the extent of the "islands of stability" around single-mode or other quasi-periodic solutions and how it varies with parameters of the physics on AdS are key questions of interest. The biggest changes occur in theories with a mass gap in the black hole spectrum, such as AdS 3 and Einstein-Gauss-Bonnet gravity in AdS 5 , which cannot form horizons at small amplitudes. While small-amplitude evolution in AdS 3 appears to be quasi-periodic [28,29], there is some evidence to point toward late-time formation of a naked singularity in AdS 5 Einstein-Gauss-Bonnet gravity [30,31] (along with a power law energy spectrum similar to that at horizon formation). Charged scalar and gauge field matter [32] also introduces a qualitative change in that initial data may lead to stable evolution or instability toward either Reissner-Nordström black holes or black holes with scalar hair.
In this paper, we extend the study of massive scalar matter initiated in [33,34]. Specifically, using numerical evolution of the full gravitational dynamics, we draw the phase diagram of gravitational collapse as a function of scalar field mass and initial scalar profile width. By considering the time to horizon formation as a function of the initial profile's amplitude, we identify several different classes of stable behavior and indicate them on the phase diagram. Finally, we analyze and characterize these different stable behaviors. Throughout, we work in AdS 5 , due to its relevance to strongly coupled gauge theories in four dimensions and because previous literature has indicated massless scalars lead to greater instability than in AdS 4 (the main other case considered), which makes the effects of the scalar field mass more visible.
We note briefly two caveats for the reader. First, horizon formation always takes an infinite amount of time on the AdS conformal boundary due to the usual time dilation effects associated with horizons; this agrees with the understanding of thermalization in the CFT as an asymptotic process. Horizon formation times discussed in this paper correspond to an approximate notion of horizon formation that we will describe below, but alternate measures of thermalization may be of interest. Second, the black holes we discuss are smeared on the compact X 5 dimensions of the gravitational side of the duality, as in most of the literature concerning stability of AdS, and we are particularly interested in small initial amplitudes that lead to black holes small compared to the AdS scale. As described in [35][36][37], small black holes in this situation suffer a Gregory-Laflamme-like instability toward localization on X 5 (which may in fact lead to formation of a naked singularity). At the same time, certain light stable solutions for charged scalars (boson stars) are stable against localization on X 5 [38]. We therefore provisionally assume that the onset of the Gregory-Laflamme-like instability occurs only at horizon formation, not at any point of the earlier horizon-free evolution.
The plan of this paper is as follows: in section II, we review the time scales associated with horizon formation with an emphasis on the behavior of massive scalars and briefly discuss our methods. Then, in section III, we present the phase diagram of different stability behaviors, and an attempt at heuristic analytic understanding appears in IV. We close with a discussion of our results.

II. REVIEW
In this section, we review results on the stability of scalar field initial data as well as our methods (following the discussion of [34]).
A. Massive scalars, stability, and time scales As in most of the literature, we work in Schwarzschild-like coordinates, which have the line element (in asymptotic AdS d+1 ) in units of the AdS scale. In these coordinates, a horizon appears at A(x, t) = 0, but reaching zero takes an infinite amount of time (measured either in proper time at the origin or in conformal boundary time); following the standard approach, we define a horizon as having formed at the earliest spacetime point (as measured by t) where A drops below a specified threshold defined in §II B below. Of course, horizon formation represents a coarsegrained description since the pure initial state of the dual CFT cannot actually thermalize; a more precise indicator of approximate thermalization may be the appearance of a power law energy spectrum (exponentially cut off) in the perturbative scalar eigenmodes. This indicator is tightly associated with horizon formation (though see [30,31] for some counterexamples).
A key feature of any perturbative formulation of the gravitational collapse is that deviations from A = 1, δ = 0 appear at order 2 , where is the amplitude of initial data. As a result, horizons can form only after a time t ∼ −2 ; in the multiscale perturbation theory of [6,7,[9][10][11][13][14][15][16][17], there is in fact a scaling symmetry → , t → t( / ) 2 that enforces the proportionality t H ∝ −2 , where t H is the (approximate) horizon formation time for unstable initial data at small amplitude.
As a result, initial data can be divided into several classes with respect to behavior at low amplitudes, as illustrated in figure 1 for massless scalars. Stable initial data evolves indefinitely without forming a horizon. In practice, we identify this type of behavior in numerical evolutions by noting rapid horizon formation at high amplitude with a vertical asymptote in t H just above some critical amplitude. In our numerical results, we see a sudden jump at the critical amplitude to evolutions with no horizon formation to a large time t lim , possibly with a small window of amplitudes with large t H just above the critical amplitude. In a few cases, we have captured a greater portion of the asymptotic region. See figure 1a. Unstable initial data, in contrast, forms a horizon at all amplitudes following the perturbative scaling relation t H ∝ −2 as → 0. In our analysis, we will verify this scaling by fitting t H to a power law as described in section II B below; if we limit the fit to smaller values of , the scaling becomes more accurate. Figure 1b shows unstable data. The red curve is of the form t H = a −2 + b with a, b determined by matching the curve to the data for the largest two amplitudes with t H ≥ 60 (not a best fit); note that the data roughly follows this curve. The categorization of different initial data profiles with similar characteristic widths into stable and unstable is robust for massless and massive scalars [34]; small and large width initial data are unstable, while intermediate widths are stable. One of the major results of this paper is determining how the widths of initial data in these "islands of stability" vary with scalar mass.
A priori, there are other possible types of behavior, at least beyond the first subleading order in perturbation theory. Metastable initial data collapses with t H ∝ −p with p > 2 at small amplitudes (or another more rapid growth of t H as → 0). We will find this type of behavior common on the "shoreline" of islands of stability where stable behavior transitions to unstable. As we will discuss further below, metastable behavior may or may not continue as → 0; in principle, as higher order terms in perturbation theory become less important, the behavior may shift to either stable or unstable as described above. We in fact find circumstantial evidence in favor of the different possibilities. Figure  1c shows metastable initial data that continues to collapse to times t H ∼ 0.6t lim but more slowly than −2 ; note that t H for collapsed evolutions at small amplitudes lies significantly above the curve t H = a −2 + b (which is determined as in figure 1b). There was one additional type of behavior identified by [34], which was called "quasi-stable" initial data at the time since the low-amplitude behavior was not yet clear. We find here that these initial data are typically stable at small amplitude but exhibit irregular, often strongly non-monotonic or even chaotic, behavior in t H as a function of , so we will denote them as irregular initial data. Figure 1d shows an example of irregular initial data. Later, we will see more striking examples of this behavior for massive scalars.

B. Methods
For spherically symmetric motion, the Klein-Gordon equation for scalar mass µ can be written in first order form as where Π is the canonical momentum and Φ = φ ,x is an auxiliary variable. The Einstein equation reduces to constraints, which can be written as where the mass function M asymptotes to the conserved ADM mass at the boundary x = π/2. We will restrict to d = 4 spatial dimensions. Since results are robust against changes in the type of initial data [34], we can take the initial data to be a Gaussian of the areal radius in the canonical momentum and trivial in the field. Specifically, The width σ and field mass µ constitute the parameter space for our phase diagram. We solve the Klein-Gordon evolution equations (2,3) and Einstein constraint equations (4,5) numerically using methods similar to those of [20] on a spatial grid of 2 n + 1 grid points; we discuss the convergence properties of our code in the appendix. We denote the approximate horizon position x H and formation time t H by the first point such that A(x H , t H ) ≤ 2 7−n . In detail, we evolve the system in time using a 4th-order Runge-Kutta stepper and initially use a 4th-order Runge-Kutta spatial integrator at resolution n = 14. If necessary, we switch to a 5th-order Dormand-Prince spatial integrator and increase resolution near horizon formation. Due to time constraints, we do not increase the resolution beyond n = 21 for any particular calculation; if a higher resolution would be required to track horizon formation for a given amplitude, we exclude that amplitude.
To determine the stability class of initial data with a given width σ, we allow evolutions to run to a maximum time of t lim = 500 in AdS units, so t lim is a lower limit for t H for amplitudes that do not form a horizon within that time. Normally, however, if the initial data appears unstable, we only evolve amplitudes with t H 0.6t lim ; this is partly to save computational resources and partly to distinguish stable evolutions from collapsing ones. For unstable or metastable initial data, we find the best fit of the form t H = a −p + b to evolutions with t H > t f it , where t f it is a constant time chosen such that amplitudes with evolutions that last longer are usually roughly perturbative; in practice, t f it = 60 gives results close to the perturbative result p = 2 for evolutions expected to be unstable from the literature, but we will also consider t f it = 80, 100 as described below.

III. PHASES
Here we give our main result, the phase diagram of stability classes as a function of initial profile width and scalar mass, along with a more detailed discussion of the scaling of horizon formation time with amplitude for varying initial data.
The stability phase diagram for spherically symmetric scalar field collapse in AdS 5 , treating the width σ of initial data and scalar field mass µ as tunable parameters, appears in figure 2. Each (µ, σ) combination that we evolved numerically is indicated by a circle, with filled and empty circles representing the unstable and stable phases respectively. The metastable phase is represented by circles filled in the top half, while the irregular phase are filled in the right half. At a glance, two features of the phase diagram are apparent: as µ increases, the island of stability moves toward smaller values of σ and takes up a gradually larger range of σ. To be specific, the stable phase is centered at σ =σ ∼ 1.4 and has a width of ∆σ ∼ 0.7 for µ = 0, 0.5, withσ ∼ 1.2 for µ = 1. ∆σ increases to ∼ 1.1, and the island of stability is centered atσ ∼ 0.9 for µ = 5, 10, while ∆σ ∼ 1.2 for µ = 15, 20 with the stable phase centered atσ ∼ 0.8. Note that the transition between "light field" and "heavy field" behavior occurs for µ > 1 in AdS units.
The metastable and irregular phases appear at the shorelines of the island of stability, the boundary between unstable and stable phases. In particular, the slope of the power law t H ∼ −p as → 0 increases as the width moves toward the island of stability, leading to a metastable phase. We find metastability at the large σ shoreline for all µ values considered and also at the small σ shoreline for several scalar masses. It seems likely that metastable behavior appears in only a narrow range of σ for larger µ, which makes it harder to detect in a numerical search, leading to its absence in some parts of the phase diagram. We find the irregular phase at the small σ shoreline for every mass and at the large σ boundary for large µ, closer to stable values of σ than the metastable phase. This phase includes a variety of irregular and non-monotonic behavior, as detailed below. Truly chaotic behavior especially becomes more prominent at larger values of µ, as we will discuss below. (blue circles). The orange line is the best power law fit.

A. Metastable versus unstable phases
While the stable and irregular phases are typically apparent by eye in a plot of t H vs , distinguishing the unstable from metastable phase is a quantitative task. As we described in section II B, we find the least squares fit of t H = a −p +b to all evolutions with t H > t f it for the given (µ, σ), running over values t f it = 60, 80, 100. Using the covariance matrix of the fit, we also find the standard error for each fit parameter. We classify a width as having unstable evolution if the best fit value of p is within two standard errors of p = 2 for t f it = 60, 80 or one standard error for t f it = 100 (due to a smaller number of data points, the standard errors for t f it = 100 tend to be considerably larger). 1 Considering larger values of t f it helps to ensure that the particular initial profile does not reach the perturbative regime at smaller amplitude values.
The fits t H = a −p + b allow us to explore the time scale of horizon formation across the phase diagram, for example through a contour plot of one of the coefficients vs σ and µ. In most cases, this has not been informative, but an intriguing feature emerges if we plot the normalization coefficient a vs σ for unstable initial data at small σ, as shown in figure 3 for t f it = 60. By eye, the coefficient is reasonably well described by the fit a = 32.0(3)σ −2.01 (2) (values in parentheses are standard errors of the best fit values) independent of scalar field mass. This is not born out very well quantitatively; the reduced χ 2 for the fit is χ 2 /d.o.f.= 180, indicating a poor fit. However, the large χ 2 seems largely driven by a few outlier points with large scalar mass, so it is tempting to speculate that the gravitational collapse in this region of parameter space is driven by gradient energy, making all fields effectively massless at narrow enough initial σ. The picture is qualitatively similar if we consider the parameter a for t f it = 80, 100 instead. Several examples of metastable behavior appear in figure 4. These figures show both data from the numerical evolutions (blue dots and red triangles) and fits of the form t H = a −p + b for points with t H > t f it = 60 (magenta curves). The best fit parameters are given in table I. Figures 4a,4b demonstrate behavior typical of most of the instances of metastable initial data we have found; specifically, the initial data continue to collapse through horizon formation times of t H ∼ 0.6t lim but with p significantly greater than the perturbative value of p = 2. Note that the evolutions of figure 4b have been extended to larger values of t H to demonstrate that the evolutions continue to collapse to somewhat smaller amplitude values. Figure 4b is also of interest because its best fit value p ≈ 2.08 is approximately as close to the perturbative value as several stable sets of initial data but has a smaller standard error for the fit, so the difference from the perturbative value is more significant. Figure 4c shows metastable evolution to t H 0.6t lim but then a sudden jump to stability until t = t lim . In the  figure, the fit has been extended to the largest non-collapsing amplitude, which demonstrates that there is no collapse over a time period significantly longer than the fit predicts. This example argues that metastable data may in fact become stable at the smallest amplitudes. On the other hand, figure 4d shows a similar jump in t H to values t H < t lim ; evolution at lower amplitudes shows metastable scaling with p ≈ 6 for 360 < t H < t lim . The figure also shows a metastable fit with larger reduced χ 2 at larger amplitudes corresponding to t f it < t H < 0.4t lim . So this is another option: metastable behavior may transition abruptly to metastable behavior with different scaling (or possibly even perturbatively unstable behavior) at sufficiently small amplitudes. It is also reasonable to classify this case as irregular due to the sudden jump in t H ; we choose metastable due to the clean metastable behavior at low amplitudes.

B. Behaviors of the irregular phase
We have found a variety of irregular behaviors at the transition between the metastable and stable phases which we have classified together as the irregular phase; however, it may be better to describe them as separate phases. The phase diagram 2 indicates that the irregular phase extends along the "inland" side of the small σ shoreline and at least part of the large σ shoreline of the island of stability. What is not clear from our evolutions up to now is whether each type of behavior appears along the entire shoreline or if they appear in pockets at different scalar field masses. Examples of each type of behavior that we have found appear in figure 5.
The first type of irregular behavior, shown in figure 5a, is monotonic (t H increases with decreasing as usual), but it is not well fit by a power law. In fact, this behavior would classify as metastable by the criterion of section III A in that the power law of the best fit t H = a −p +b is significantly different from p = 2, except for the fact that the reduced χ 2 value for the fit is very large (greater than 10) and also that different fitting algorithms can return significantly different fits, even though the data may appear to the eye like a smooth power law. In any case, this type of behavior apparently indicates a breakdown of metastable behavior and hints at the appearance of non-monotonicity. So far, our evolutions have not demonstrated sudden jumps in t H typical of stability at low amplitudes, however. Figure 5b exemplifies non-monotonic behavior in the irregular phase. This type of behavior, which was noted already by [18], involves one or more sudden jumps in t H as decreases, which may be followed by a sudden decrease in t H and then resumed smooth monotonic increase in t H . There are suggestions that this type of initial data is stable at low amplitudes due to the usual appearance of non-collapsing evolutions, but it is worth noting that these amplitudes could instead experience another jump and decrease in t H , just at t H > t lim . Finally, [34] studied this type of behavior in some detail, denoting it as "quasi-stable." The last type of irregular behavior is apparently chaotic, in that t H appears to be sensitive to initial conditions (ie, value of amplitude) over some range of amplitudes. This type of behavior appears over the range of masses (see figure 1d for a mild case for massless scalars), but it is more common and more dramatic at larger µ. Figures 5c,5d represent the most extreme chaotic behavior among the initial data that we studied with collapse at t H < 50 not very far separated from amplitudes that do not collapse for t < t lim along with an unpredictable pattern of variation in t H . This type of chaotic behavior has been seen previously in the collapse of transparent but gravitationally interacting thin shells in AdS [39] as well as in the collapse of massless scalars in AdS 5 Einstein-Gauss-Bonnet gravity [30,31]. In both cases, the chaotic behavior is hypothesized to be due to the transfer of energy between two infalling shells, with horizon formation only proceeding when one shell is sufficiently energetic. In the latter case, the extra scale of the theory (given by the coefficient of the Gauss-Bonnet term in the action) leads the single initial pulse of scalar matter to break into two pulses.
We should therefore ask two questions: does this irregular behavior show evidence of true chaos, and is a similar mechanism at work here? To quantify the presence of chaos, we examine the difference in time evolution between similar initial conditions (nearby amplitudes), which diverge exponentially in chaotic systems. Specifically, any quantity ∆ should satisfy |∆| ∝ exp(λt) for Lyapunov coefficient λ. Our characteristic will be the upper envelope of the Ricci scalar at the origin per light crossing time,R(t). We consider the chaotic behavior exhibited by three states: a massless scalar of width σ = 1.1 with amplitudes = 1.02, 1.01, 1.00 (see figure 1d), a µ = 5 massive scalar of width σ = 0.34 and = 3.52, 3.51, 3.50, and a µ = 20 scalar of width σ = 0.19 and = 6.98, 6.95, 6.92 (figure 5d). Figure 6 details evidence for chaotic evolution in the µ = 5, σ = 0.34 case; figure 6a shows our characteristic function R(t) for the amplitudes 1 = 3.50, 2 = 3.51, and 3 = 3.52. By eye,R shows noticeable differences after a long period of evolution. These are more apparent in figure 6b, which shows the log of the differences ∆ ab ≡R a −R b , along with the best fits. Although there is considerable noise -or oscillation around exponential growth -in the differences (leading to R 2 values ∼ 0.2, 0.26 for the fits), the average slope gives Lyapunov coefficient λ = 0.007 (within the error bar of each slope), and each slope differs from zero by more than 3 standard errors. One interesting point is that the t H vs curve in figure 5b does not appear chaotic to the eye, even though it shows some of the mathematical signatures of chaos at least for 1 < < 3 , the amplitude values near the visible spike in t H .
The story is similar for the massless and µ = 20 cases we studied, which exhibit λ values that differ from zero by at least 1.9 standard deviations; see table II. This is a milder version of the chaotic behavior noted by [30,31,39], especially for the µ = 5 case studied. To our knowledge, this is the first evidence of chaos in the gravitational collapse of a massless scalar in AdS. One thing to note is that the strength of oscillation in log(|∆|) around the linear fit increases with increasing mass, so that the two best fit Lyapunov exponents for µ = 20 are no longer consistent with  each other at the 1-standard deviation level. The mechanism underlying the chaotic behavior seems somewhat different or at least weaker than the two-shell or Einstein-Gauss-Bonnet systems. When examining the time evolution of the mass distributions of these data, we see a single large pulse of mass energy that oscillates between the origin and boundary without developing a pronounced peak. However, there is also apparently a smaller wave that travels across the large peak. In the massless case examined, this wave deforms the pulse, leading at times to the double-shoulder appearance of figure 7a. In the µ = 5, σ = 0.34 case, the secondary wave is more like a ripple, usually smaller in amplitude but more sharply localized, as toward the right side of the main pulse in figure 7b. So the chaotic behavior may be caused by the relative motion of the two waves, rather than energy transfer between two shells. In this hypothesis, a horizon would form when both waves reach the neighborhood of the origin at the same time.
As a note, we have run convergence tests on several sets of irregular initial data and find that our calculations are convergent overall, as expected (even at lower resolution than we used). In particular, the massless scalar evolutions studied in table II are convergent already at resolution given by n = 12 (note that we typically start at n = 14). The only caveat may be for some of the strongly chaotic initial data with scalar mass µ = 20, which nonetheless appear well-behaved according to other indicators. The reader may or may not wish to take them at face value but should recall that we have presented other chaotic initial data with rigorously convergent evolutions. See the appendix for a more detailed discussion.

IV. SPECTRAL ANALYSIS
As we discussed in the introduction, instability toward horizon formation proceeds through a turbulent cascade of energy to shorter wavelengths or, more quantitatively, to 1st-order scalar eigenmodes with more nodes. Inverse cascades are typical of stable evolutions. Therefore, understanding the energy spectrum of our evolutions, both initially and over time, sheds light on the behavior of the self-gravitating scalar field in asymptotically AdS spacetime, providing a heuristic analytic understanding of the phase diagram.

A. Dependence on mass
The most visibly apparent feature of the phase diagram of figure 2 is that the island of stability both expands and shifts to smaller widths as the scalar mass increases. As it turns out, the energy spectrum of the Gaussian initial data (7) provides a simple heuristic explanation.
It is well established both in perturbation theory and numerical studies that initial data given by a single scalar linear-order eigenmode is in fact nonlinearly stable, and the spectra of many quasi-periodic solutions are also dominated by a single eigenmode. As a result, we should expect Gaussian initial data that approximates a single eigenmode (which must be j = 0 due to lack of nodes) to be stable. To explore how this depends on mass, we find the best fit values of , σ for the j = 0 eigenmode for each mass that we consider (defined by the least-square error from the Gaussian to a discretized eigenmode); this is the "best approximation" Gaussian to the eigenmode. Then we find the energy spectrum of that best-fit Gaussian; these are shown in figure 8a. From the figure, it is clear that the j = 0 eigenmode is closer to a Gaussian at larger masses. That is, other eigenmodes contribute less to the Gaussian's spectrum at higher masses (by several orders of magnitude over the range from µ = 0 to 20). Simply put, the shape  of the j = 0 eigenmode is closer to Gaussian at higher masses, which suggests that the island of stability should be larger at larger scalar field mass. Figure 8b compares the j = 0 eigenmode and best fit Gaussian for µ = 0 and 20; on inspection, there is more deviation between the eigenmode and Gaussian for the massless scalar.
In addition, the best-fit Gaussian width decreases from σ ∼ 0.8 for a massless scalar as the mass increases. At µ = 20, the best-fit width is σ ∼ 0.31. This suggests that Gaussians that approximate the j = 0 mode well enough are narrower in width at higher masses. An interesting point to note is that the island of stability for µ = 0, 0.5 is actually centered at considerably larger widths than the best-fit Gaussian. This may not be surprising, since the best-fit Gaussians at low masses actually receive non-negligible contributions from higher mode numbers; moving away from the best-fit Gaussian can actually reduce the power in higher modes. For example, the stable initial data shown in figure 9a below has considerably less power in the j = 2 mode.

B. Spectra of different phases
A key question that one might hope to answer is whether the stability phase of a given (µ, σ) can be determined easily by direct inspection of the initial data without requiring many evolutions at varying amplitudes. The initial energy spectra for examples of each phase, including monotonic, non-monotonic, and chaotic irregular behaviors, are shown in figure 9. These spectra are taken from among the smallest amplitudes we evolved in order to minimize backreaction effects.
Unfortunately, the initial energy spectra do not seem to provide such a method for determining the stability phase. Very broadly speaking, stable and metastable (µ, σ) correspond to initial spectra that drop off fairly quickly from the j = 0 mode as j increases, while unstable and irregular phases tend to have roughly constant or even slightly increasing spectra up to j = 5 or 10. However, figure 9d shows that some irregular initial data have spectra that decrease rapidly after a small increase from j = 1 to j = 2. Kinks in the spectrum are more prevalent for widths of the AdS scale or larger, while spectra for smaller widths tend to be smoother.

C. Evolution of spectra
While the initial spectrum for a given (µ, σ) pair does not have predictive value regarding the future behavior as far as we can tell, the time dependence of the spectrum throughout the evolution of the system is informative. Figure 10 shows the time-dependence of spectra for examples of the stable, unstable, metastable, and chaotic irregular phases. In each figure, the lower panel shows the fraction E j /M ADM in each mode up to j = 6, while the upper panel shows the cumulative fraction j E j /M ADM to the mode 2 k with k = 0 to 5.
The difference between stable evolution in figure 10a and unstable evolution in figure 10b is readily apparent. As the evolution proceeds, we expect a cascade of energy into higher mode numbers, but inverse cascades to lower modes can also occur. The stable evolution shows a slow pattern of cascades and inverse cascades, in fact. On the other hand,    : Initial (t = 0) energy spectra for the indicated evolutions. In order, these represent stable, unstable, metastable, monotonic irregular, non-monotonic irregular, and chaotic irregular initial data.
the unstable evolution shows a nearly monotonic cascade of energy into the highest modes along with a simultaneous cascade of energy into the lowest modes (therefore depleting intermediate modes). These are common observations in the literature and are included here for completeness. The metastable evolution shown in figure 10c is interesting in light of the stable and unstable spectra. The amplitude shown is from the "unstable" portion of figure 4d, the part consistent with the perturbative scaling t H ∼ −2 . However, the spectrum shows a similar pattern of slow cascades and inverse cascades to the stable phase example, though on a somewhat faster time scale in this case. While perhaps surprising, this is in keeping with the similarities noted between the initial spectra in figures 9a and 9c. We have also checked that the time-dependent spectrum at a higher amplitude with t H ∼ 100 follows the same pattern as 10c; in fact, it looks essentially the same but simply ends at an earlier time. This lends some credence to the idea that the metastable phase is stable at lowest nontrivial order in perturbation theory, with instability triggered by higher-order corrections. Alternately, the instability could be caused by an oscillatory singularity in the perturbative theory, as discussed in [15,[23][24][25] in the case of two-mode initial data. These divergences do not appear in the energy spectrum. Figure 10d shows the time-dependence of the spectrum in a chaotic irregular evolution, specifically µ = 20, σ = 0.19 at = 6.95, which is in the chaotic region of the t H vs plot in figure 5d. There is rapid energy transfer among modes, including cascades out of and inverse cascades into mode numbers j ≤ 32 over approximately a light-crossing time. It is easy to imagine that horizon formation might occur at any of the cascades of energy into higher modes, leading to seemingly random jumps in t H as a function of amplitude.
Finally, the time-evolved energy provide another possible measure of approximate thermalization in the dual CFT; namely, the spectrum should approach an (exponentially cut-off) power law at thermalization. In most cases, this occurs shortly before horizon formation, but there are exceptions, such as the late time behavior of initial data below the critical mass for black hole formation in Einstein-Gauss-Bonnet gravity [30]. In the case of chaotic behavior, it is particularly interesting to know if the spectra for similar amplitudes approach a power law at similar times even if horizons form at very different times. Figure 11 shows the energy spectra for two amplitudes in the chaotic region of the t H vs plot for µ = 0, σ = 1.1. Figure 11a is the spectrum just before horizon formation for = 1.01, while figure  11b is the spectrum at approximately the same time for = 1.02, which is very long before horizon formation. In this example, we see that the spectrum does approach a power law for the evolution that is forming a horizon, while the other evolution demonstrates apparently exponential decay. Therefore, this example suggests that a power law spectrum may yield similar results to horizon formation as a measure of thermalization in the dual CFT.

V. DISCUSSION
We have presented the phase diagram of stability of AdS 5 against horizon formation, treating the scalar field mass µ and width σ of initial data as free parameters. In addition to mapping the location of the so-called "island of stability," we have gathered evidence for two non-perturbative phases on the "shorelines" of the island, the metastable and irregular phases. While these must either exhibit stability (no collapse below some critical amplitude) or instability (collapse at arbitrarily small but finite amplitude) as the amplitude → 0, they are distinguished by their behavior at computationally accessible amplitudes. While perturbatively unstable evolutions obey t H ∝ −2 as → 0, metastable initial data follows t H ∝ −p for p > 2 over a range of amplitudes. The irregular phase is characterized by horizon formation times t H that are not well described by a power law and sometimes exhibit non-monotonicity or even chaos. Both of these phases appear across the range of µ values that we study and at both small-and large-width boundaries of the stable phase.
At this time, it is impossible to say whether metastable initial data is stable or unstable as → 0 (or if all metastable data behaves in the same way in that limit). Our numerical evolutions include cases in which the lowest amplitudes jump either to metastable scaling with smaller p or to evolutions that do not collapse over the timescales we study. In many cases, too, the power law t H ∝ −p with p some fixed value > 2 is robust as we exclude larger amplitudes from our fit. It is also possible that some of the metastable phase is stable in the perturbative theory (ie, to 3 order in a perturbative expansion) but not at higher orders.
The irregular phase seems likely to be (mostly) stable at arbitrarily small amplitudes based on our numerical evolutions, though we have not found a critical amplitude for monotonic irregular initial data. The irregular phase includes the "quasi-stable" initial data described in [18,34], which has a sudden increase then decrease in t H as decreases as well as truly chaotic behavior. In fact, we have found evidence for weakly chaotic behavior at the jump in t H for non-monotonic initial data in the form of a small but nonzero Lyapunov coefficient. Both non-monotonicity and chaos become stronger and more common at larger scalar masses; however, we have also found chaotic behavior for the massless scalar. To our knowledge, this is the first evidence of chaos in spherically symmetric massless scalar collapse in AdS, which is particularly interesting because there is only one physically meaningful ratio of scales, σ as measured in AdS units.
Aside from the ultimate stability or instability of the metastable and irregular phases, several questions remain. For one, black holes formed in massive scalar collapse in asymptotically flat spacetime exhibit a mass gap for initial profiles wider than the Compton wavelength 1/µ [42]. Whether this mass gap exists in AdS is not clear, and it may disappear through repeated gravitational focusing as the field oscillates many times across AdS; investigating this type of critical behavior will likely require techniques similar to those of [43]. Returning to our phase diagram, the physical mechanism responsible for chaos in the irregular phase is not yet clear. Is it some generalization of the same mechanism as found in the two-shell system? Also, would an alternate definition of approximate thermalization in the dual CFT, such as development of a power-law spectrum, lead to a different picture of the phase diagram? Finally, the big question is whether there is some test that could be performed on initial data alone that would predict in advance its phase? So far, no test is entirely successful, so new ideas are necessary.  initial data is defined analytically, so Q can appear poor at t = 0 since the errors are controlled by round off; in some cases, Q is therefore undefined and not plotted. First, we carried out convergence tests for mass µ = 0.5, width σ = 1, and amplitude = 1.12, which is monotonic irregular initial data presented in figure 5a. This amplitude collapses with t H ∼ 88. Figure 12a shows the (L 2 norm) order of convergence for the field variable φ, the mass function M , and the metric functions A, δ. While the order of convergence is initially poor and even negative, all these variables show approximately fourth order convergence for times t 70. The reason for the initially poor convergence is that the error between successive resolutions is already given by (machine limited) round off. As a demonstration, we tested the order of convergence with base resolution n = 11, as shown in figure 12b. The variables show order of convergece Q 3 already at this resolution for most of the evolution, losing convergence only for t > 80, where we see approximately 4th-order convergence in the n = 14 resolution computations.
Two of the authors have discussed the convergence properties of evolution for the nonmonotonic irregular initial data with µ = 20, σ = 0.1, = 11.74, which is in an amplitude region of increased t H surrounded by smaller values, in detail in [34]. In short, the variables φ, M, A, δ all exhibit fourth order convergence, as does Π 2 (t, 0), and the conserved mass actually has 6th order convergence.
Initial data for µ = 15, σ = 0.2 is also nonmonotonic, as shown in figure 13a. While we have not analyzed all aspects of the convergence, we see from the remainder of figure 13 that φ, M, A, δ exhibit convergent behavior at better than second order for = 7. amplitude in figure 13a). It is important to note that the larger amplitude also has the larger horizon formation time, contrary to the usual monotonic behavior. It is most crucial to validate the convergence of chaotic evolutions. In table II, we noted that the Ricci scalar at the origin has nonzero Lyapunov exponent at almost the 2 sigma level for amplitudes = 1.02, 1.01, 1.00 for µ = 0, σ = 1.1. We show the results of convergence tests for these amplitudes in figure 14; because these are longer evolutions, we consider the convergence at the lower resolutions n = 12, 13, 14. After a transient start-up period, these are all convergent with Q > 2.5 for all variables considered at all times; for most of the time, the order of convergence is Q > 3.5. It is worth noting that one of the amplitudes does not form a horizon through t = 500.
Initial data with µ = 1, σ = 1 is chaotic over a narrow range of amplitudes. We have carried out convergence testing for amplitudes = 1.15, 1.14, which are the two amplitudes with t H < 100 between amplitudes with t H 150 in figure 15a. The order of convergence was poor for these amplitudes in our initial tests with base resolution n = 14 because the error between resolutions was dominated by round-off, similar to the convergence tests we discussed above for µ = 0.5, σ = 1. In subsequent tests with lower resolutions n = 11, 12, 13, we find an order of convergence Q ∼ 4 for most of the evolutions (and always Q > 3). It is important to note again that our evolutions exhibit convergence while showing horizon formation at a later time for a larger amplitude in this case.
Finally, we ran convergence tests for the chaotic initial data with µ = 20, σ = 0.19 for = 6.95, 6.92, with t H ≈ 65.5, 40.8 respectively. As shown in figure 16, the simulations are close to fourth order convergence for most of the evolution, but there are periods where the order of convergence for evolution and constraint variabables becomes negative. This of course leads to the concern that the evolutions should have collapsed during those periods and extend into an "afterlife" evolution. We have therefore evolved these amplitudes through these regions (approximately t = 30 − 40 for = 6.95 and t = 18 − 30 for = 6.92) at high resolution (n = 18). If the evolutions are truly in an afterlife, this higher resolution calculation may include horizon formation. We do not observe this. Another tell-tale of would-be horizon formation is a decrease in the timestep size by an order of magnitude or more followed by an increase. We monitor the timestep size every 500 timesteps through this evolution but do not observe a decrease in timestep size by more than a factor of 2. As a result, we believe the values of t H found are reliable, though the reader may wish to consider them with some caution. Nonetheless, we emphasize that we have found convergent behavior for chaotic initial data at other scalar masses.