Classes of Stable Initial Data for Massless and Massive Scalars in Anti-de Sitter Spacetime

Since horizon formation in global anti-de Sitter spacetime is dual to thermalization of a conformal field theory on a compact space, whether generic initial data is stable or unstable against gravitational collapse is of great interest. We argue that all the known stable initial data for massless scalars are dominated by single scalar eigenmodes, specifically providing strong numerical evidence consistent with the interpretation that initial data with equal energies in two modes collapse on time scales of order the inverse square of the amplitude. We further scan the parameter space for massive scalar field initial data and present evidence for a novel class of stable or quasi-stable solutions for massive scalars with energy spread through several eigenmodes.


Introduction
Through the AdS/CFT correspondence, gravitational physics on (global) d + 1-dimensional anti-de Sitter spacetime (AdS d+1 ) is dual to a conformal field theory (CFT) on R × S d−1 (the boundary of AdS), with classical gravity valid in the strong coupling regime of the CFT. Thermal states of the CFT are dual to black holes in the bulk AdS spacetime, 1 so thermalization of an initial distribution of energy in the CFT is dual to horizon formation, ie, gravitational collapse, in the bulk AdS. It seems natural to expect that even small amounts of added energy confined to a compact space generically eventually thermalize, so we are led to a perhaps surprising conclusion that generic initial data for matter in AdS should lead to black hole formation. On the gravitational side of the correspondence, the reason is that massless fields can reach spatial infinity and reflect off the boundary in finite time, so energy cannot disperse as it would in asymptotically flat spacetime (similarly, massive fields are confined by the effective gravitational potential of AdS).
From this point of view, there are two natural questions. First, are there classes of initial data that are stable against gravitational collapse to a black hole in the bulk AdS (especially with measure greater than zero), and what is the dual CFT interpretation of this initial data? Second, in the case of gravitational collapse, how long does given initial data take to form a horizon -or, in CFT terms, how long does given initial data take to thermalize? Strictly speaking, both bulk horizon formation and CFT thermalization take infinite boundary time, so we can formulate the question more precisely by asking how long it takes for energy to spread through a large number of frequency modes, which typically happens just before an approximate horizon forms in the bulk. 2 We therefore take bulk horizon formation, as we will define approximately in section 2 below, as an indication of boundary CFT thermalization as well. This second question is actually easier to answer in the unstable case: at low perturbation amplitudes, self-gravitation can only become important on time scales of order −2 , where is the amplitude. In fact, the first question is often phrased in terms of stability over times of order −2 .
In 2011, [1] pioneered numerical work on these questions (see also [2][3][4]). Specifically, [1] presented intriguing numerical evidence that perturbations of AdS are generically unstable to black hole formation, at the expected time scale for low amplitudes. The original studies of massless scalar field matter have since been followed by studies of complex scalar fields [5] and modified theories of gravity [6]. Except when there is a mass gap in the black hole spectrum (such as AdS 3 [7,8] or AdS 5 Einstein-Gauss-Bonnet gravity [6]), these numerical studies have pointed to instability to horizon formation at arbitrarily small amplitudes for generic initial data, which occurs along with a cascade of energy into high frequency modes. On the other hand, certain initial data appear stable at low amplitudes, and additional types of initial data are the subject of some disagreement (see [9][10][11]).
The last several years have seen a simultaneous effort to understand AdS gravitational collapse in perturbation theory, particularly with the development of an expansion in terms of the scalar eigenmodes on a fixed AdS background [9,[12][13][14]. Among other symmetries, the perturbation theory obeys a scaling law A(t) → A(t/ 2 ) (where A is the coefficient of an eigenmode in the expansion), which is also observed in numerical calculations at low amplitude [1] and, as suggested above, leads to a universal prediction that horizon formation takes a time of order −2 in the perturbative regime regardless of field mass or higher curvature terms in the gravitational action [15]. The same scaling symmetry is also apparent in a perturbative calculation for a thin-shell scalar field profile [16]. However, the implications of the perturbation theory for stability of generic initial data at arbitrarily small amplitude are as yet unclear. Conservation laws of the perturbation theory imply the possibility of inverse cascades of energy returning to low-frequency modes, which have been verified in numerical simulations (see [9], for example, and our results) and which may be associated with stability at low amplitudes. Further, [16,17] have argued that (quasi-)periodic solutions of the perturbation theory should persist as stable solutions of the full theory at arbitrarily small amplitudes. On the other hand, [18] have argued for a singularity in the time derivative of the phases B n in generic solutions of the perturbation theory, 3 which may be related to horizon formation in the full nonlinear theory. As a result, stability, instability, or both may be generic in the space of initial data (in the sense of being open sets).
Given the existing tension, it seems useful to take stock of the characteristics of initial data that are agreed to be stable at low amplitudes. From the initial studies, both numerical and perturbative analysis agree that (nonlinearly modified) scalar eigenmodes are stable [1,20]; these are known as oscillons or, in the complex scalar case, boson stars. 4 In fact, perturbations around the oscillon solutions are expected to be stable as well [22], in agreement with numerical studies. We note here that all subsequently developed solutions for which there is numerical evidence of stability on long time scales are perturbed oscillons in the sense that their energy spectrum is dominated by the contributions of a single perturbative eigenmode; more precisely, one mode has at least twice the energy of any other individual mode. For example, initial data of Gaussian shape and width near the AdS scale is apparently stable [20,23,24] because it is actually dominated by a single scalar eigenmode; somewhat wider initial data requires a stronger admixture of other modes and is not quasistable. Likewise, the periodic solutions of [25,26] are dominated by a single mode as shown by their spectra, and, though [9,27] present a more general ansatz for quasi-periodic solutions to the perturbation theory, the physical solutions presented are all dominated by one mode. The highest temperature quasi-periodic solution shown in figure 1 of [27] still has the j = 0 mode having approximately twice as much energy as the j = 1 mode. It is also worth noting that, in most cases, a single mode dominates the spectrum even more strongly in the sense that it has more than 60% of the total energy. Those initial data for which apparently stable numerical solutions of the full theory have been presented in the literature are all of the latter type. While high temperature quasi-periodic solutions of the perturbation theory do exist with considerably less energy in the dominant mode, the behavior of these has not to our knowledge been evaluated in the full theory. As the reader may suspect, there are also controversial cases: [9,11] hinted at stability of initial data superposing two eigenmodes based on a perturbative calculation, which has been criticized by [10].
In this paper, we take a numerical approach to the stability of scalar fields in AdS at small amplitude, specifically looking at two questions. First, in section 3, we consider massless scalars in AdS 4 and provide an independent numerical analysis of the two-mode initial data of [9] as well as the superposition of three Gaussians used in [28]. We find that, when the field is stable against gravitational collapse over times that grow faster than −2 with decreasing amplitude, the energy spectrum is dominated by a single scalar eigenmode as described above.
3 [19] appeared while this work was in the final stages of preparation and argues the singularity is not present in a different gauge, specifically where the time coordinate is boundary time. 4 In the case of metric perturbations, geons as described in [21] are the equivalent stable eigenmodes.
(This condition is similar to the condition of [27] that stable solutions are "close" to a quasiperiodic solution of the perturbation theory.) Our calculations also point out difficulties in determining the reliability of the perturbative expansion even in the low amplitude regime when analytical arguments are not available.
In section 4, we then turn our attention to the case of scalars with mass µ = 0 (corresponding to backgrounds for irrelevant operators in the CFT), which have so far only been studied numerically in [28]. (We do not consider tachyonic scalars which are allowed for mass squared above the Breitenlohner-Freedman bound [29].) Given that massive fields have different stability properties in asymptotically flat spacetime, this is an important test case for stability in AdS.
As a review, consider the gravitational collapse of massive scalars in asymptotically flat spacetime, where an initial configuration can either form a black hole or disperse to infinity (ignoring self-interactions that, for example, lead to star birth in astrophysics). Given initial data of characteristic width λ and λµ 1, [30] found that collapse proceeds as for a massless scalar, with power-law scaling behavior for the horizon radius near the critical point between black hole formation and dispersion [31] (known as a type II transition, see [32]). However, for λµ 1, black holes form only above a finite mass gap (a type I transition). In essence, the difference in these two types of transition is whether the scalar potential or gradient energy dominates the energy of the initial data. For a massive scalar field in asymptotically AdS spacetime with curvature radius , the different ratios to be considered are λµ, µ and λ/ . In principle, the stability properties of the system can change whenever these ratios pass through unity.
We present an overview of the behavior of horizon formation time as a function of these ratios in section 4, followed by a more detailed analysis of the behavior at low amplitudes. For widths λ between the Compton wavelength and AdS scale, we find two types of particularly interesting behavior, depending on whether the scalar is light or heavy compared to the AdS scale. For light fields, we find a discontinuity in the horizon radius at formation as a function of the amplitude of initial data, possibly corresponding to a change in efficiency of thermalization in the boundary CFT. We analyze this behavior in more detail in section 5. Of particular interest, for heavy scalars, we find a class of stable (over time scales of order −2 ) initial data; we present evidence in section 6 that these solutions are qualitatively distinct from the single-eigenmode-dominated solutions discussed above in that the energy is distributed roughly evenly through several modes.
We conclude with a discussion of our results and begin here with an overview of our methods.

Methods
The evolution of a massive scalar field in AdS is governed by the Einstein equations and the wave equation, where the mass of the scalar field φ is µ. Following [1], we choose a spherically symmetric metric ansatz in Schwarzschild-like coordinates where d is the number of spatial dimensions. The areal radius is R(x) = tan(x/ ), and we henceforth work in units of the AdS scale (ie, = 1). The equations of motion are calculated using a Hamiltonian analysis, as in [33]. The evolution equations governing the nonlinear system are where Π = A −1 e δ φ ,t is the canonical momentum and Φ = φ ,x is an auxiliary variable. The metric functions are determined by where M is the mass function and M (t, x = π/2) ,t = 0. We choose δ(t, x = 0) = 0.
The linearized system is given by (see [34] for the mathematical properties of this system). The eigenfunctions of the operator L are given by Jacobi polynomials (see [29] for a review) (2.7) The normalization and eigenvalues are given by In this work we choose λ + , corresponding to the normalizable modes according to the inner product (f, As in [1], we define an "energy per mode" by projecting the field onto eigenmodes. Specifically, we define Π j := √ AΠ, e j , φ j := (φ, e j ), and , e j , so the energy spectrum is The total ADM energy is M ADM = ∞ j=0 E j . We use these projections to study how much energy is captured in modes up to some j max and how energy is transferred between modes throughout the evolution.
As far as we are aware, only initial data (ID) that is either one or more Gaussians in Π [1, 3, 5-7, 20, 24, 28], a superposition of two eigenmodes [9,10], (fake) boson stars [20], or a specially constructed time-periodic solution [25] has been studied. We consider four types of initial data. We primarily study Gaussian data in Π, similar to that originally studied in [1]. For a Gaussian pulse of width σ, like eq. (2.10), we define the wavelength to be λ = 2σ. As discussed in the introduction, one goal of this work is to explore the interplay among the length scales λ, , and 1/µ. To this end, we have considered several different values of σ and µ.
To explore the universality of our results, we also consider initial data in the form of an ingoing pulse This data is more difficult to evolve numerically, specifically as collapse ensues, so we do not perform simulations for all the values of σ and µ used for the Gaussian initial data in Π. Nonetheless, as we present in section 4.3 below, our results appear robust against changes in initial data. For comparison to existing literature, we also consider two other classes of initial data. First, [9] studied two-eigenmode initial data in AdS 4 , specifically, where κ is chosen freely. Finally, [28] considered a superposition of three Gaussian wavepackets in Π for AdS 4 . We comment on both these classes of initial conditions as they appear in the previous literature in addition to studying them ourselves.
Our numerical methods are similar to those used in [6] and are described in greater detail in [24]. We use a fixed grid of 2 n grid points (with the option to restart at higher resolution), mainly with a 4th order Runge-Kutta (RK4) spatial integrator. For some initial data we find it necessary to use a fifth-order Dormand-Prince (DoPr) method for the spatial integration. The DoPr method is typically only needed when the solution approaches collapse and the finer spatial scales need to be resolved more accurately. We terminate the simulation and determine that a horizon forms at A(t H , x H ) ≤ 2 k−n with k = 7 for µ ≤ 20, while for larger masses k = 8. We have found that requiring smaller values of k in some cases requires the formation of very narrow troughs in the horizon function A(t, x) which are not resolved by the spatial grid, so we have chosen these values to reliably report horizon formation times. A discussion of the convergence properties of our code is in appendix A.

Stability of massless scalars in AdS 4
In this section, we present results of simulations in AdS 4 as a comparison to the previous literature. We are particularly concerned with two sets of initial data which have been claimed to lead to stable solutions at low amplitude: two-mode initial data and three-Gaussian initial data.

Two-mode initial data
First, we consider massless scalars with two-mode initial data as in eq. (2.12) and κ = 3/5, following [9,11,27] and [10], which is a point of disagreement between the two sets of publications. Specifically, [9] suggests stability of sufficiently small amplitude two-mode data for the length of their simulations, including at = 0.09, while [10] found horizon formation at a long but finite time. We undergo an independent study of this initial data for amplitudes = 0.109, 0.09 and 0.08 and find that all three simulations end in collapse. Namely, comparing the bottom panel of figure 1a to figures 3 and 4 of [9] and figure 1 of [10], shows agreement with Bizoń and Rostworowski [10]. Below we provide evidence that the discrepancy is due to insufficient resolution in [9,11].
Our simulations also show the same scaling properties first noticed in [1], namely that −2 Π 2 ( 2 t, x = 0) is an approximate universal function for a given initial field profile as long is the simulation is still in the perturbative regime. This scaling symmetry is manifest in the improved perturbation theory of [9,12,13] (also known as "two-time formalism" or TTF), as discussed in the introduction. The close agreement after rescaling, as seen in figure  1a, provides evidence that the numerical solutions are still perturbative until shortly before horizon formation. Nonetheless, the perturbative TTF solution diverges from the nonlinear numerical solution with the same initial data (including our calculation and those of [9][10][11]27]). Specifically, in figure 3 of [9], Π 2 (t, x = 0) in the TTF evolution including modes up to j = 47 falls below the numerical result by approximately an order of magnitude at t ≈ 200, while the numerical solution still agrees with [10] and our results. Even the TTF solution including modes up to j = 200 [27] diverges strongly from the numerical solution;  To gain a better understanding of why the numerical and TTF solutions diverge, we perform a spectral decomposition of the energy for = 0.09. The top panel of figure 1b shows the total energy in modes up to j max , while the bottom panel shows the energy in the lowest six modes. 5 We first see a cascade of energy into higher modes, followed by an inverse cascade of energy back into the lower modes, in agreement with [9]. There is a strong cascade into higher modes just before horizon formation. Figure 1b shows that at this point about 4% of the energy is in modes greater than j = 47 and almost 2% in j > 256. Because higher modes are more sharply peaked around the origin, even a small amount of energy in higher modes can cause large differences in the magnitude of Π 2 (x = 0, t) and in determining whether or not a collapse occurs. The lesson is that perturbative solutions may only be reliable if a large number of eigenmodes are considered, even when the full solution is well within the perturbative regime. The correct claim of [11] that the TTF evolution with j max = 47 accurately captures the dynamics of the lowest modes is immaterial to the question of horizon formation because this is a local question near x = 0 in the small amplitude limit. It appears that a small amount of energy transferred to higher modes is sufficient to cause gravitational collapse/thermalization.
We have, however, not yet explained the difference in the non-perturbative numerical solutions of [9] and [10] at long times. In principle, one possibility is that these two references use different time coordinates (ours matches that of [10]   numerical coordinate transformation to the time coordinatet of [9,11], which is chosen such that δ(t, x = π/2) = 0 (this is conformal time on the boundary of AdS). The two coordinate times are related byt we show the relationship between t,t for late times in figure 2a for the = 0.09 initial data. Because the system is weakly gravitating until collapse, ∆t/∆t ≈ 1.02 for most of the simulation; however, as the horizon starts to form, time dilation becomes more significant, and the boundary time passes more quickly. We find that collapse is delayed byt H − t H ≈ 20 for the = 0.09 evolution for resolution and horizon formation threshold given by n = 19. This is insufficient to account for the lifetimet ≈ 1500 of the simulations of [9]. 7 Furthermore, time dilation only stretches the time axis of figures 1 and 2b and would not add extra oscillations as in figures 3 and 4 of [9]. Instead, the difference between [9,11] and [10] appears entirely due to insufficient resolution in the simulations of [9]. Specifically, the consistency and convergence tests presented in [9] were stopped much earlier than the times in the evolution that are in question, the longest ending att ≈ 600. As an illustration, the simulation depicted in red in figure 2b has grid spacing ∆ ≈ 9.6 × 10 −5 , a higher resolution than used in any of the tests in [11]. The rapid decrease in the ADM mass shows that late in the evolution our code also suffers from loss of convergence if insufficient resolution is used, and the evolution experiences a similar 7 Technically speaking, a lower threshold for horizon formation results in greater time dilation. However, even though we have been unable to locate a discussion of the threshold used in [9], we think this is an unlikely source of the error simply because our threshold is quite low.
"afterlife" to those of [9,11]. We have determined that this is due to the largest amplitude pulse sharpening and eventually being squeezed between grid points. Figure 2b also shows higher resolution simulations where the relative change in ADM mass is ∆M/M ≈ 7 × 10 −8 for the entirety of the evolution.
Unfortunately, the convergence tests presented in [11] are insufficient to evaluate whether or not the numerical solution is converging to the continuum solution. Typically when performing convergence tests either the number of grid points is increased by a constant amount or the resolution is increased by a constant factor. The convergence test presented in figure 1 of [11] does not do either of these. Rather, what the figure shows is that their code converges to the solution with ∆ ≈ 1.2 × 10 −4 as ∆ asymptotically approaches this value. Meanwhile, we have also compared numerical values of Π 2 (t, x = 0) with [10] and obtain the same results to several significant figures; an exact match is not expected due to slight differences in the algorithms used. The results of our study suggest that the disagreement between [9] and [10] would be resolved if [9] used sufficient resolution to be in the convergent regime late in the evolution.
At this point, it is worth assessing the importance of analyzing such apparently welltrodden ground (as [11] seemingly wonders). As we stated in the introduction above, one of the key questions about gravitational collapse in AdS is what classes of initial data can lead to stable evolution at small amplitude. While numerical studies cannot answer this question definitively, it is still critical to be careful about whether numerics have actually provided evidence of stability or not. In this case, equal-energy two-mode initial data, if stable, would present a qualitatively distinct class of stable initial data -all the apparently stable initial data for massless scalars is dominated by a single eigenmode. Indeed, some of the authors of [9] postulated in [14] that the two-mode initial data is close to a quasi-periodic solution (of a type formulated in [9]), but even the quasi-periodic solution closest to their two-mode initial data has 70% of its energy in the j = 0 eigenmode. It is also important to realize that, while perturbative methods are very powerful, they also suffer from stringent resolution requirements (ie, many scalar eigenmodes must be included) because gravitational collapse is an essentially local question at low amplitude. In other words, as [27] points out, it is the behavior of the very high j tail of the energy spectrum that determines the ultimate fate of the system. In summary, while it is impossible at this time to rule out stability of equal-energy two-mode initial data as → 0, all current numerical evidence is consistent with instability to horizon formation over times of order −2 precisely due to the contributions of high j modes.

Multiple-Gaussian initial data
We also study the multiple Gaussian initial data of eq. (2.13) for massless scalars in AdS 4 , which was first considered in [28] with parameters We have simulated the collapse of this initial data and find, as in [28], that it is apparently stable for small amplitudes. To understand why this data is stable, we performed a spectral decomposition as a function of time. The top   panel of figure 3a shows the total energy up to mode j max , while the bottom panel shows the energy in each of the lowest four modes, both for = 1. We find that approximately 82 per cent of the energy is in the lowest mode, essentially independent of time. This suggests that this initial data is a perturbation about a single mode, which is known to be stable [1,20,25].
For the same simulation we plot the upper envelope of Π 2 (t, x = 0) in the top left panel of figure 3b. At least for the duration of our simulation, there is no increase in the Ricci scalar at the origin, as is naively expected for stable solutions. Our initial data for = 1 is plotted in the top right panel of figure 3b and appears to match the bottom right panel of figure 9 in [28]. We observe qualitatively similar behavior to [28], but our quantitative results differ significantly. For large amplitudes we observe rapid horizon formation, but, with decreasing , we find a small region where the collapse time quickly increases and then decreases again. Finally, near = 5.756 we observe another rapid increase in formation time. Yet another small decrease in collapse time is observed followed by another increase. It is unclear if this behavior will recur indefinitely. In contrast, [28] observes a sudden increase in collapse time at ≈ 1.75 and no earlier increase (see their figure 8).

Collapse of massive fields in AdS 5
We now turn to an overview of the behavior of massive scalar fields in AdS. For specificity, we work in AdS 5 , which corresponds to a dual 4D gauge theory.
Our primary motivation is to consider different values of the dimensionless parameters λµ, µ, and λ/ and to determine when the behavior of massive fields diverges from that of massless fields. For µ < 1 (light fields), we consider initial conditions with width satisfying λµ < λ/ < 1, λµ < 1 < λ/ , and 1 < λµ < λ/ . Similarly, for µ > 1 (heavy fields), we consider the λ/ < λµ < 1, λ/ < 1 < λµ, and 1 < λ/ < λµ cases. In most cases, we find that decreasing initial amplitude leads to monotonically increasing horizon formation times, though the behavior may differ in detail from massless scalars. However, as for massless scalars, we also find evidence of stability against horizon formation at small amplitudes for λ ∼ , when a single eigenmode dominates the spectral distribution. We have also uncovered an apparently distinct class of (quasi-)stable initial conditions, which are not dominated by any single scalar eigenmode. We discuss those in more detail in section 6.
We consider multiple classes of initial conditions, demonstrating that the behavior displayed is robust.

Overview
For narrow pulses, λ less than both 1/µ and , we expect the scalar fields to act effectively massless, whereas wider pulses may exhibit different behavior. More precisely, since [30,35] found a phase transition at λµ ∼ 1, we also may expect a transition in behavior as the initial pulse width increases past either the AdS radius or Compton wavelength. However, because the scalar field can disperse and then reflect back to the origin (recall that massive fields are confined by a gravitational potential in AdS, even though they cannot travel to the AdS boundary), the initial distinction in pulse width may be erased after sufficient reflections, leading to common behavior at low amplitudes (indeed, this is suggested by the perturbative analysis of [15]). The AdS radius also sets the scale of the gravitational potential, so the cross-over between light and heavy fields and narrow initial pulses compared to may also change the behavior of the gravitational collapse. Heuristically, we expect the behavior of the collapse to be controlled by which type of energy dominates: gradient energy, scalar potential, or gravitational potential.
We consider first light scalars, picking µ = 0.5/ for specificity. To start, we choose Gaussian initial data in Π as given in equation (2.10) for a variety of widths σ. As expected, narrow initial pulses (λ = 2σ = 0.6 ) lead to gravitational collapse behavior similar to that of massless fields, as shown in figure 4a. In particular, the horizon formation time forms slightly sloped "steps" as a function of pulse amplitude, with each step separated by a jump of ∆t H increasing from approximately 2.8 at large amplitudes to π at small amplitudes. The initial horizon radius x H follows the characteristic arc pattern familiar from [1] with one arc per step in t H (though the steps occur rapidly enough at small amplitude that the arcs are difficult to resolve). The behavior for pulses wide compared to the AdS scale, shown in figure 4d for λ = 8 , is similar. The arcs in x H are less distinct at low amplitudes, and more investigation would be required to determine the behavior of the horizon radius in this case. Additionally, the steps in t H are more steeply sloped, leading to a noticeably larger collapse time even for amplitudes that collapse promptly (without reflections from the boundary). At intermediate widths, < λ < 1/µ, however, a different structure emerges. Specifically, in the amplitude regime where the scalar field undergoes prompt collapse, the initial horizon radius undergoes a sudden decrease (not associated with a step in t H ). This behavior is apparent in figure 4b    We now turn to the case of heavy scalars; figure 5 shows our results for a field of mass µ = 5/ . For narrow pulses λ < 1/µ, we find again, as expected, behavior similar to massless fields for the amplitudes we study; results are shown in figure 5a for σ = 0.05 . However, the case of heavy scalars and narrow initial data is numerically challenging, likely due to the narrow field pulse propagating close the boundary, where the φ/ cos 2 (x) term in the scalar equation of motion (2.3) becomes subject to large numerical error (since both numerator and denominator are small). Wider initial pulses are less numerically challenging. Very wide initial pulses (λ > ) as in figure 5d show a lengthened collapse time for prompt collapse surpassing t H = π/2, which corresponds to the length of time needed for the broadly distributed field to collapse toward the center of AdS. The subsequent steps in t H appear more compressed as a function of , leading to a rapid increase in horizon formation time with decreasing amplitude. This is apparent to some extent also in figure 6a for mass µ = 10/ . We also consider intermediate widths 1/µ < λ ≤ . At σ = 0.5 as in figure 5c, the step pattern in t H has mostly disappeared, replaced by a continuously varying behavior down to a critical amplitude ≈ 2.75. (There appear to be steps in t H just above this amplitude, but, since the step is ∆t H < π, it seems that the continuous behavior has simply become steep.) We have not been able to form a horizon below this critical amplitude, even while allowing the simulations to run to t = 500 (see section 4.2 below for further information). As it turns out, we should expect this initial data to be stable; over 91% of the energy of the initial data is carried by the j = 0 eigenmode, so this is a single-mode-dominated solution.
At a somewhat lower widths still in the intermediate range, we find more interesting behavior that becomes more apparent at larger masses. For µ = 5/ , σ = 0.3 , we find behavior very similar to the massless scalar, at least for solutions with t H 85 . However, at the larger mass µ = 20/ and intermediate width σ = 0.1 , the horizon formation time t H has a sudden narrow jump before decreasing again; see figure 6b. At lower amplitudes, t H increases rapidly again (faster than −2 ), leading to apparently stable behavior. Note that this behavior is distinct from single-mode-dominated solutions as exemplified in figure 5c -in fact, as discussed in section 6 below, the energy is distributed democratically throughout several eigenmodes. In distinction to the single-mode-dominated oscillon and their quasi-periodic generalizations, our results suggest the presence of stable oscillaton solutions for heavy scalar fields. For appropriate values of width, then, our initial data can approach the oscillaton solutions, leading to an extended collapse time. This effect may also be a manifestation of the dynamical mass gap found by [30] in a regime where the mass term dominates over gradient energy but the pulse is too narrow to be deformed by the AdS curvature. We revisit these issues in more detail in section 6.

Low amplitudes
As mentioned in the introduction, one of the most important questions is the dependence of t H on . At small amplitudes, self-gravitation is suppressed by two powers of , so gravitational collapse should not occur until a time of order −2 has past. In fact, the scaling symmetry of the perturbation theory [9,12] makes this argument precise; if a solution remains perturbative until a time t 0 , then a solution with lower amplitude will also remain perturbative until at least time t 0 / 2 . Since gravitational collapse should require non-perturbative gravitational physics, then t H 1/ 2 . While first presented for massless scalars, this argument generalizes to massive scalars and Einstein-Gauss-Bonnet gravity [15]. Since it would be noteworthy if gravitational collapse happened on a faster time scale (in the perturbative regime), we investigate the functional dependence of t H on in some of our simulations. Similarly, if t H grows faster than −2 with decreasing amplitude, that is an indication of quasi-stability (if not absolute stability).
Specifically, we fit t H = a p + b to low amplitude data points, both with and without the constraint b = 0, and take rough agreement between the values of p in the two fits as an indicator that we are in the perturbative regime. We have also checked that the results of the fit are robust against removal of some of the low amplitude data points. Results appear in table 1 and are generally consistent with t H ∝ −2 scaling for the initial data shown in the table. Due to time constraints, we have not carried out calculations for all our initial data into   that initial data that collapses with t H 100 is typically not yet in the perturbative regime. We also noted that scalars with mass µ = 5/ and initial data width σ = 0.5 lead to apparently stable evolution below a critical amplitude ≈ 2.75, and we have been unable to induce gravitational collapse below this amplitude. To support this interpretation, we have run simulations at amplitudes = 2.74, 2.67, 2.49 to time t = 500 at a base resolution of n = 14. The ADM mass is conserved to 2 parts in 10 12 over this time for all these simulations, and we have also run convergence tests (to t = 235, 500, 275 respectively) to verify agreement at higher resolution. As we noted above, we expect this initial data to be stable at low amplitudes, since it is single-mode dominated, and it remains dominated by the j = 0 mode to late times. Specifically, we find that over 80% of the energy remains in the j = 0 mode to t = 500 for all three amplitudes, and the high j spectrum decays exponentially. Interestingly, the initial data appears to become stable outside of the perturbative regime -the Ricci scalar at the origin does not respect the perturbative scaling symmetry at these three amplitudes, and the fraction of the energy held in the j = 0 mode averaged over long times increases with decreasing amplitude.

Robustness against change of profile shape
We have also examined the stability behavior of other initial data to determine the robustness of results. Specifically, we have considered an ingoing pulse as in eq. (2.11) and two-mode initial data (2.12), both in AdS 5 .
We considered two values of (µ, σ) for comparison of the ingoing pulse to the Gaussian initial data in Π, µ = 0.5 with σ = 0.3 and µ = 20 with σ = 0.1. The initial horizon radii and formation times are shown in figure 7 and share the basic characteristics of the corresponding mass and width in figures 4a and 6b. In particular, the µ = 20, σ = 0.1 case has the same striking sudden increase and reduction in t H as the amplitude decreases. This provides clear numerical evidence that our results are robust against change in the shape of initial data.
Since single mode oscillons and perturbations around them are stable [1,9], the simplest possible data that could produce a horizon is two-mode initial data. We choose the modes to have equal energy so they are as far away from a single-mode-dominated solution as possible and for comparison with previous work [1,9] in AdS 4 and [18] in AdS 5 . With this choice, the characteristic width λ of the initial data is fixed for a given mass; changing the ratio λµ of    Figure 7: x H and t H vs for ingoing pulse initial data the width to the Compton wavelength requires considering several masses. We consider mass values µ = 0, 0.5, 1.0 and √ 10 and estimate the width by fitting a Gaussian to the profile (the massless case was previously studied in [18]). For µ = 0.5, λµ < 1; for µ = 1, λµ ≈ 1; and for µ = √ 10, λµ > 1. The initial horizon radius x H and horizon formation time t H in these cases are qualitatively similar to the results shown in section 4.1 above.
In particular, we find no evidence that these solutions are close to a stable region. For µ = 0, 1/2, 1 and √ 10, we find only a direct cascade of energy to higher modes for the duration of our simulations, in agreement with the results of [18]. The top panel of figure 8 shows the turbulent energy cascade for the µ = √ 10 and = 0.02. While we cannot rule out that inverse cascades appear at lower amplitudes, we note that this initial data already appears to be in the perturbative regime, as Π 2 (t, x = 0) obeys the perturbative scaling law, 8 as we see in the bottom panel of figure 8. This suggests that AdS 5 may have smaller stable regions around single mode data than AdS 4 , at least for the values of µ we considered. While most of the behavior is similar for massive and massless scalars with two-mode initial data, Π 2 (t, x = 0) displays one interesting difference: it increases smoothly with small oscillations for massive scalars, while for massless scalars it increases in a piecewise linear fashion.
We also consider the late-time energy spectrum for all the types of initial data we study. In collapsing solutions, [24] provides numerical evidence that the spectrum takes the form is E j ∝ j −α where α = 6/5 + 4(d − 3)/5. We show the late-time spectrum in figure 9 for several masses of two-mode initial data as well as ingoing wave and Gaussian in Π initial data. These data have been evolved to long times (t ≈ 1100 for two-mode initial data and t ≈ 350 for the others). Remarkably the exponent α appears to be independent of both the type of initial data and the scalar mass. This suggests that it is the gravitational physics that sets the decay rate, rather than the matter content.
fit=0.57j Figure 9: Energy spectra late in evolution. Black line shows E j ∝ j −2 . Initial data are two-mode except Π ID as eq. (2.10) and φ ID as eq. (2.11).

Initial horizon radius decrease
In this section we present a more detailed discussion of the rapid decrease in the initial horizon radius for length hierarchy λ < 1/µ, as seen in figures 4b and 4c. As shown in figure 10a, we have also found similar behavior for a mass µ = 0.1/ and width σ = . While the decrease in x H is not as dramatic in this case, the left-most data point does show a clear drop, indicating that this behavior may be common. The astute reader will also notice several sudden increases in the initial horizon radius as decreases in figures 4b, 4c, and 10a. These increases, though again less striking, are familiar from the original work of [1] on massless scalar collapse in AdS 4 .
Both of these phenomena are related to the existence of multiple local minima in the metric function A, which we use as a determinant of horizon formation; a typical-looking example of this function just prior to t H for µ = 0.5/ , σ = 0.8 , = 1.98 is shown in figure 10b (solid red curve). 9 Recall that we terminate our calculations and declare horizon formation when A decreases below a resolution-dependent threshold; the minimum in A first reaches the threshold at (t H , x H ). When A has multiple local minima (which are associated with shells of infalling matter), the initial horizon radius x H is the position of the first local minimum to reach the threshold. It is easy to see how decreasing the initial amplitude might lead to a jump in x H , then: at some value of , the innermost shell of matter is dense enough to push the inner minimum of A below the threshold, but, at a slightly lower value of , the inner minimum does not contain quite enough mass to decrease past the threshold. Horizon formation must wait until another shell of mass approaches the origin, but the jump in mass corresponds to a jump in the Schwarzschild radius. It is worth noting that t H increases negligibly in this process; the formation of the inner minimum is associated with growth in the metric function −δ, leading to time dilation that allows the outer mass shells to fall inward while very little proper time t at the origin passes. Previous literature [6] has confirmed that this behavior remains under increasing resolution, which corresponds to decreasing the threshold, but the location of the jump may not be robust as decreasing the threshold can cause the jump to shift to higher amplitudes. The reverse behavior -a sudden decrease in the horizon radius -is more curious. For a fixed resolution (n = 14 in this case), decreasing the amplitude from = 1.98 to 1.95 changes the metric function A from the multi-minimum red solid curve to the single-minimum blue dashed curve in figure 10b. Counter-intuitively, the minimum for the lower amplitude is deeper than the inner minimum of the higher-amplitude curve. This behavior at fixed amplitude is not robust against decreasing the threshold for horizon formation (which, in our analysis, occurs with increased resolution); we have verified for = 1.95 that the initial minimum in A (pictured in figure 10b for resolution n = 14) does not pass the threshold for resolution n = 16 and that a second local minimum of A appears; the sudden decrease in horizon radius should shift to lower amplitudes. However, to be clear, this is a real feature of the solutions and not a convergence issue -the n = 14 and higher resolution n = 15, 16 solutions agree to roundoff, at least until the horizon threshold is reached for n = 14. We have also verified convergence of the solution at lower resolution (base resolutions n = 10, 12). This is yet another way in which gravitational collapse in AdS is sensitive to initial conditions and also numerical algorithms. These results also stress that it is crucial to use the same parameters when comparing simulations to each other -while most comparisons will not evince this type of sensitivity, we have not been able to predict when such sensitivity will appear. In vernacular, it is important to compare apples to apples. It is at first glance unclear whether we should interpret this phenomenon in the dual gauge theory; the initial horizon radius depends on the spatial slicing used to describe the gravitational collapse, and nothing physical in the boundary CFT should depend on a gauge choice in that way. In addition, the sharp decrease in x H apparently depends sensitively on the threshold for horizon formation. However, let us momentarily give in to temptation and give our results an interpretation in the boundary theory. Since we take approximate horizon formation as corresponding to partial thermalization, we can think of the mass inside the initial horizon as corresponding to the fraction of the boundary energy that is approximately thermalized. If we take these results at face value, they tell us that decreasing the initial amplitude past a critical value leads to approximate thermalization with less energy transferred into higher modes. In figure 11, we display energy spectra for our samples with µ = 0.5/ , σ = 0.8 , and amplitudes = 1.95 and 1.98 which support this view to some degree. Specifically, for = 1.95, a significantly greater fraction of the energy remains in the two lowest modes, indicating that a smaller amount of energy has approximately thermalized. Conversely, the energy fraction in higher modes (j ≥ 100, for example) is significantly larger for the = 1.98 evolution.
The boundary time of horizon formation provides a possible explanation for this phenomenon. Although both amplitudes lead to similar values of t H as measured in terms of proper time at the origin, they have quite different conformal timet H at the boundary due to time dilation between the inner minimum of A and the boundary. At resolution n = 14, we have t H ≈ 0.85 andt H ≈ 2.28 for = 1.95 as compared to t H ≈ 0.82 andt H ≈ 4.04 for = 1.98. In essence, our measure of thermalization (A passing a fixed threshold) has allowed the higher amplitude solution a longer time to thermalize in the boundary theory. Indeed, if we increase resolution to n = 16, the = 1.95 solution develops more minima in A, and the energy spectrum 11a evolves to resemble figure 11b more closely. For future work, it would be interesting to determine what other simple indicators beside a threshold for A make accurate measures of partial thermalization in the boundary theory.

Quasi-stable solutions
In section 4.1, we noted that massive scalars with intermediate pulse widths 1/µ < λ < exhibit a curious behavior of t H as the amplitude decreases, which we specifically observed in figure 6b for µ = 20, σ = 0.1 with initial data given by eq. (2.10). At large amplitudes, t H increases in the typical piecewise (roughly) constant fashion. The increase then becomes quite rapid, followed by a sudden jump in t H by more than a factor of 3. Very quickly, t H decreases again, followed by a steady but rapid increase. This type of behavior has appeared in the previous literature for massless [20] and massive [28] scalars. In the massless scalar case, the corresponding initial data leads to stable, single-mode-dominated oscillon solutions at low amplitudes; in this section, we provide preliminary evidence that the corresponding ID for massive scalars may represent a novel class of stable solutions -oscillatons -with energy spread democratically through multiple modes. While this behavior also appears for ingoing wave initial data as in figure 7b, in this section we will restrict to ID in the form (2.10) for simplicity.
At the lowest amplitudes we have studied for µ = 20, σ = 0.1, we find that t H rapidly increases as decreases, which provides some support for possible stability at arbitrarily small amplitude. These calculations require very high resolution (n ≥ 17) to remain convergent until collapse, and we have carried out convergence testing to be confident in our calculations. In table 2 below, we show fits t H = a p + b for the indicated amplitude ranges (which are entirely below the large jump in t H noted in section 4); due to computational complexity, these fits all have ≥ 9.62. Except for the smallest amplitude range 9.62 ≤ ≤ 9.8, we find a disagreement in p values for fits constrained to have b = 0 or not, and we also find that the smaller ranges have increasingly negative values for the power p when b = 0. This indicates that t H = a p + b does not provide a good fit to t H ; we have tried also fits of the form t H = a( − * ) p + b, but those also fail to provide a good fit (in fact, common fitting algorithms do not find sensible fits at all). The lesson of table 2 is not that the horizon formation time follows a specific power law in amplitude. What is evident, though, is that t H increases much more rapidly than −2 with decreasing amplitude, so this initial data appears   to be quasi-stable, that is, stable over longer time scales than expected from generic arguments about perturbation theory. Given constraints on computing time, we have not yet ascertained whether t H undergoes repeated jumps such as the one seen near = 11.74 in figure 6b or reaches a critical amplitude below which horizons do not (apparently) form. However, it seems reasonable to conclude that the low amplitude behavior indicates (quasi-)stability over longer than normal time scales at the least. We plot the energy in the lowest six modes alongside the upper envelope of Π 2 (t, x = 0) in figure 12 for = 9.5. An inverse energy cascade occurs at t ≈ 35 and then again at t ≈ 130 with a rapid increase in Π 2 occurring at t ≈ 307. There is a pattern of increases and decreases in Π 2 and an approximate recurrence in the spectrum before collapse, both consistent with direct and indirect energy cascades. This is similar to the orbits around quasi-periodic solutions studied by [27] for massless scalars; however, the spectrum in figure  12a is distinct in that the the j = 0, 1, 2 modes all carry the greatest share of the energy at some point in time (compare to the spectrum of equal-energy two-mode ID evolution in figure 1b). Interestingly, during the final growth in Π 2 (t ≈ 311), ∼ 99% of the energy is in the lowest 16 modest.
The presence of inverse cascades tend to coincide with rapid increase in formation time and has been used to argue indirectly for stability in certain solutions (for example, [9]). As figure 12b shows, it is possible for collapse to occur even long after an inverse cascade with no obvious way of predicting whether or not a rapid cascade of energy to higher modes will occur later in the evolution. As we found for two-mode data in section 3.1, the evolution may be apparently quasi-periodic for long times and then abruptly change, possibly leading to horizon formation. Another example of this type of behavior is found in the region near ≈ 11.74, where t H jumps from approximately 37.7 to 126 and back to 75.6. We have carried out convergence testing on numerical evolution of three amplitudes = 11.742 (t H ≈ 37.73), = 11.739 (t H ≈ 126.02), and = 11.736 (t H ≈ 78.71) to verify the expected 4th-order convergence. The data for the = 11.739 simulation is presented in appendix A as an example of our convergence testing. Furthermore, it is worth noting that we find the long-lived (t H ≈ 126.02) behavior for more than one amplitude value. As we show in figure 13, evolution at these three amplitudes is in remarkable agreement until immediately before horizon formation. While that is to some extent to be expected due to the very small changes in amplitude (about 0.06% over the whole region), it is striking how unpredictable horizon formation seems in a comparison of the three amplitudes. In fact, the curvature for the intermediate amplitude, which has the largest t H , starts growing but then turns over, allowing a lower amplitude evolution to form a horizon first. Figure 13 also shows the low j part of the spectrum, which evinces a similar pattern of recurrences as the lower amplitude ( = 9.5) evolution shown in figure 12a, though the higher-amplitude evolutions have a number of higher-frequency oscillations on top of the long-term recurrence pattern. Once again, the spectra for the three amplitudes near = 11.74 are difficult to distinguish except just before one of the evolutions forms a horizon.
To investigate how common this behavior is for massive scalars, we also study extremely massive scalars with µ = 100/ . This data is numerically much more difficult to evolve than the less massive fields. Therefore, all simulations discussed here were typically run with of 2 14 or more grid points. We then ran simulations at higher resolutions to test if the results agreed, as well as did global convergence tests of the longer simulations. We find that the solutions that collapse at "late" times (after t ≈ 15) require 2 17 or more grid points to be in the convergent regime. Running simulations to several hundred time units at this resolution is a daunting task, so we present only preliminary results here. See Appendix A for a detailed discussion of the convergence of our methods. We are interested in data that is not a perturbation about a single mode, so we choose a narrow pulse of σ = 1/16; the width of the Gaussian that best fits e 0 (x) for µ = 100 is σ ≈ 7/50 (for µ = 20, it is σ ≈ 3/10). We plot the energy spectrum at t = 0 in the top panel of figure 14a. The two lowest modes have nearly equal energy, and an exponential decay in the higher modes, so this is not single-mode dominated, at least initially. We also show the evolution of the energy in the lowest four modes in the bottom panel of figure 14a. Compared to other solutions (such as for µ = 20 discussed above), we find a strikingly periodic and rapid transfer of energy into and out of the lowest mode. Based on a fit p 0 cos(Ωt + p 1 ) + p 2 to E 0 , 10 we find a frequency approaching Ω ≈ 2 as the amplitude is decreased. Interestingly, we also find that the amplitude of the oscillations decrease for smaller perturbations.
As we have noted, the stable solutions found so far are perturbations about single mode solutions, typically with exponentially decaying spectra at higher modes. In figure 14b, we provide some suggestive evidence that this initial data may be stable despite initially having nearly equal energies in the j = 0, 1 modes (and significant energy in the next two modes). We plot the horizon formation time in the bottom panel, and Π 2 (t, x = 0) and R(t, x = 0) for = 15 in the top panel. Tracking the rapid increase in formation time has proven surprisingly difficult, which means that whether or not this data truly is stable is an open question. Nevertheless it is clear that increasing the mass of the scalar results in surprisingly different dynamics from what is seen in the massless case.

Discussion
The question of stability of perturbations of AdS against gravitational collapse (either absolutely or on time scales of order −2 for amplitude ) is intriguing both as a problem in gravity and through its connection to thermalization in strongly-coupled gauge theories via the AdS/CFT correspondence. We have investigated what classes of initial data might lead to stable evolution for both massless and massive scalar fields. As reviewed in the introduction, it is well-established that perturbations around single eigenmode oscillons (or boson stars) are stable, at least on long time scales. We propose, therefore, that one criterion for stability is that the energy spectrum is initially dominated by a single scalar eigenmode; this is compatible with the conjecture of [27] that stable solutions orbit quasi-periodic solutions in configuration space. It is simple to check that most of the stable solutions for massless scalars present in extant literature are, in fact, single-modedominated. We have provided a check that the stable three-Gaussian initial data of [28] are single-mode-dominated, as well. In addition, we studied the somewhat controversial equal-energy two-mode initial data discussed also in [9][10][11]27] and and find that claims of stability based on amplitudes studied in the literature so far are premature, though stability at smaller amplitudes remains an open question. Specifically, for the amplitudes of initial data studied in the extant literature, according to a fixed approximate measure of horizon formation, gravitational collapse occurs at a time proportional to −2 . We conclude that the disagreement in the literature over stability of this initial data at small amplitudes is due to the use of insufficient resolution in the numerical solutions of [9]. It is also important to note that the perturbative TTF solutions require high resolution (in the sense of requiring a large number of eigenmodes); even a small fraction of energy in high j modes can trigger gravitational collapse since horizon formation is an essentially local process. Furthermore, other than comparison to a numerical solution of the fully nonlinear problem, there is no simple diagnostic in the literature for the failure of a perturbative solution, so it is not clear that perturbative methods are a reliable gauge of stability even over a fixed time scale.
Starting in section 4, we present a thorough overview of gravitational collapse for massive fields in AdS 5 . We consider evolution for the six regions of parameter space determined by whether each dimensionless ratio λµ, µ , λ/ is larger or smaller than unity and confirm in some examples that our qualitative results for each parameter range are robust against changes in the initial data. We also confirm that initial data in most of these parameter ranges collapse with t H ∝ −2 , as expected from perturbation theory. In the process, we further discovered two types of novel behavior related to thermalization and stability for initial data with intermediate widths.
First, for scalars with mass less than the AdS scale, we noted, in addition to the more familiar sudden increases in the initial horizon radius x H , also the appearance of sudden decreases in x H with decreasing amplitude. Both phenomena are related to the threshold for the metric function A(t, x) used to determine approximate horizon formation and the appearance of multiple local minima in A at a fixed time. In some cases, as the amplitude decreases, the first-formed minimum of A never reaches the threshold, but subsequently formed local minima created by additional infalling matter at larger x can reach the threshold. This leads to sudden increases in x H . However, we have shown that the opposite case can occur; sometimes decreasing the amplitude of the initial data actually allows the first-formed local minimum of A to decrease beyond a threshold that it cannot reach at a slightly higher amplitude. This puts a spotlight on how approximate horizon formation is defined due to the sensitivity of results to the threshold. It also motivates the development of other criteria either for approximate horizon formation or thermalization in the dual field theory.
Our most important result is to provide evidence of a new type of (quasi-)stable field profile for massive fields and intermediate width 1/µ < λ < . For µ = 100, σ = 1/16, there is an apparent critical amplitude, below which we have not been able to find horizon formation. For µ = 20, σ = 1/10, we find with decreasing amplitude first the usual increase of t H , then a sudden jump from t H ≈ 37.7 to t H ≈ 126 followed by a decrease to t H ≈ 75.6, then a rapid increase in t H . While we have not yet ascertained whether there might be a critical amplitude, our results do argue for stability over time scales of order −2 . In both cases, the energy spectrum is not dominated by a single mode, unlike for other known stable solutions, but rather is initially spread democratically through several modes, possibly indicating the existence of a quasi-stable multi-mode oscillaton solution. Following the hypothesis [27] that stable evolutions are orbits of quasi-periodic solutions, it would be interesting to determine if quasi-periodic solutions of the perturbation theory for a massive scalar allow energy to be distributed more evenly between modes. An alternate possibility is that the "island of stability" for perturbations around single-mode oscillons or quasi-periodic solutions is much wider for heavy scalars in an AdS manifestation of the mass gap found in asymptotically flat spacetime [30]. Another important point is that our µ = 20, σ = 1/10 solutions in the region when t H jumps from 37.7 → 126 → 75.6 show a significant agreement up until the time of collapse, even though the final horizon formation time differs greatly. There is no clear indication in advance whether a given amplitude will collapse sooner or later. 11 We leave a more careful study of these potentially novel stable solutions for future work.
In summary, there are two major lessons from our work. First is that, while there is by now a well-developed leading-order perturbation theory for scalar gravitational collapse in AdS, comparison to numerical solutions of the full nonlinear theory shows that it is very difficult to predict when the perturbative theory will break down. We have given several explicit examples of this difficulty; as a result, it is difficult to know how long a perturbative solution remains reliable (for example in the sense of [17]). Second, all the known stable or quasi-stable evolutions of massless scalars in AdS are dominated by a single eigenmode of the scalar on a fixed AdS background, whereas quasi-stable evolutions of the massive scalar can have energy democratically spread through several eigenmodes. These massive scalar solutions call for future investigation. There are two quantities we monitor for consistency of the solutions. The system (2.3) has a constraint, C := φ ,x − Φ = 0, and the ADM mass should be conserved by the time evolution. Since explicit convergence tests for each simulation are too demanding, we typically monitor the ADM mass and C for consistency. For small and zero µ, this is usually sufficient because we find that the mass rapidly decreases as convergence is lost (see figure 2b). This happens because the pulses slip between the mesh if insufficient resolution is used. However, as µ is increased the ADM mass is conserved well and C will remain small even when convergence is lost. A careful study of the evolution when this occurs shows that it is the formation of several small minima in A that are under resolved. These minima each reach a different value and the narrower, farther out (radially) ones are ultimately those that trigger horizon formation. We find that in many cases at least 2 17 grid points are required to resolve these minima, making the computations incredibly expensive. Nevertheless for particularly long or interesting simulations, we have run explicit convergence tests to provide credibility for our results.
We perform an array of convergence tests to determine the reliability of our results, which are presented here for µ = 20, σ = 0.1, and amplitude = 11.74. This solution has a large t H compared to slightly larger or smaller amplitudes, so it is an interesting test case for convergence and consistency of the code. Figure 16a demonstrates that the scalar field, mass function, and metric functions A and δ are fourth-order convergent (in L 2 norm) through multiple reflections off the outer boundary, and figure 16b shows that the constraint residual C decreases with increasing resolution. Figure 16c shows that the ADM mass is well-conserved for the duration of our simulation. We  plot the ADM mass at all three resolutions as well as the Richardson extrapolated value. 12 The order of convergence of the ADM mass is plotted in figure 16d and exhibits approximately sixth order convergence. Because of the variable step size method used with the DoPr integration, the ADM mass has a higher order of convergence than expected solely from the time evolution. Since Π 2 (t, x = 0) has been used to detect the onset of the turbulent instability. We show convergence of this quantity in figure 16e. Finally, we plot the convergence of the constraint residual in figure 16f, which demonstrates that the during the evolution we stay on the constraint sub-manifold. In post processing we used a second order method, unlike in the actual time evolution where we use a fourth order method. The results of our testing demonstrate that our code is convergent and consistent even for long simulations with several reflections off the outer boundary, as well as for large scalar field mass values.