Black Hole Formation in AdS Einstein-Gauss-Bonnet Gravity

AdS spacetime has been shown numerically to be unstable against a large class of arbitrarily small perturbations. In arXiv:1410.1869, the authors presented a preliminary study of the effects on stability of changing the local dynamics by adding a Gauss-Bonnet term to the Einstein action. Here we provide further details as well as new results with improved numerical methods. In particular, we elucidate new structure in Choptuik scaling plots. We also provide evidence of chaotic behavior at the transition between immediate horizon formation and horizon formation after the matter pulse reflects from the AdS conformal boundary. Finally, we present data suggesting the formation of naked singularities in spacetimes with ADM mass below the algebraic bound for black hole formation.


Introduction
It is well known that, in Einstein gravity, Minkowski space is stable against arbitrarily small perturbations. The simplest way to understand this is by noting that the formation of microscopic black holes exhibits critical behaviour, usually referred to as Choptuik scaling [2]. Specifically, as the amplitude of a small initial perturbation gradually decrease, the dynamics undergoes a phase transition between black hole formation (for large amplitudes) and dispersion to infinity (for small amplitudes). Infinitesmally small perturbations invariably disperse. As with all critical phenomena, the end state is highly sensitive to small variations in the initial data near the transition region. The properties of the phase transition depend both qualitatively and quantitatively on the nature of the local dynamics. The transition is second order in the absence of a fundamental scale in the problem, but can be either first or second order when new scales are introduced. New scales can arise either in the matter action or due to modifications to the gravitational dynamics via, for example, the addition of higher curvature terms. It has been shown [3] that adding a Gauss-Bonnet (GB) curvature squared term in five and six spacetime dimensions radically affects the critical behaviour in microscopic black hole formation. 1 Naively one might expect Choptuik scaling, which was originally discovered as a local phenomenon, to be insensitive to boundary conditions at infinity. In particular, it was thought to be unaffected by the inclusion of a cosmological constant. To the contrary, [5] argued from numerical results that, in (global) anti-deSitter space (AdS), which exhibits reflecting boundary conditions at the conformal boundary, black holes form from arbitrarily small perturbations for massless scalar matter with a large class of initial data. The instability arises because subcritical matter that initially disperses is able to return from the boundary in finite time to form a horizon near the origin after additional gravitational focusing. Subsequent analysis by many authors  has demonstrated the existence of "islands of stability," i.e. non-negligible regions of the initial data parameter space for which black holes never form. In fact, some perturbative analysis suggests that stability against horizon formation may be generic in the parameter space of initial conditions, and it is still an open question whether stability, instability, or both are technically generic at arbitrarily small but finite amplitude. Other work has considered massive scalars [33][34][35], a gauge field and charged scalar [36], and holographic models of confining theories (related to the Poincaré patch rather than global AdS) [37,38].
The stability of AdS spacetime is an interesting question in mathematical physics in its own right, but the issue takes on particular significance in the context of the AdS/CFT correspondence, in which gravity in AdS spacetime is dual to a Yang-Mills theory on the conformal boundary. Since black hole formation in the bulk AdS spacetime corresponds to thermalization in the spatially compact boundary CFT, it is perhaps less surprising to think that generic initial conditions lead to black holes. Indeed, islands of stability are more surprising as they imply that some low-energy perturbations of Yang-Mills theories on S 3 need not thermalize. However, the high degree of symmetry in AdS (integrability of the boundary theory) can lead to quasiperiodic behavior. It is clear that AdS/CFT is a rich system with many lessons about nonperturbative dynamical behavior.
The end state of gravitational collapse in AdS spacetime results from the interplay of local (weak turbulence) and global (resonance) dynamics of the spacetime. Quantum theory generically suggests the need for higher-derivative terms in the gravitational and matter actions that necessarily alter the short distance, high curvature dynamics near the final stages of gravitational collapse, i.e. the local dynamics. In the AdS/CFT correspondence, highercurvature terms in the gravitational action correspond to finite N and 't Hooft coupling effects in the dual theory, including differing a and c central charges in 4D CFTs. Our focus in this paper is the gravitational collapse of a massless scalar field in Einstein-Gauss-Bonnet (EGB) gravity in the AdS 5 /CFT 4 correspondence. We are motivated in part by the possible relation of the boundary CFT to strong dynamics in QCD.
While one expects a tower of higher-derivative couplings suppressed by powers of the string scale, 5D EGB gravity has been an important model of higher-curvature effects in the AdS/CFT correspondence because it is the first example of Lovelock gravity [39,40] beyond the Einstein-Hilbert action. 2 The key feature of Lovelock terms in the gravitational action is that the equations of motion remain second order in derivatives of the metric despite the fact that the action is higher order. Not only does this imply that the theory is ghost-free when linearized around a flat background, but it also makes the study of AdS stability tractable. In 5D, only the lowest order Lovelock term (beyond Einstein), the Gauss-Bonnet term, is relevant. As a result, it is the unique higher-curvature theory of gravity with second-order equations of motion.
The present authors initiated a study of the stability of AdS in EGB gravity in [1]. The purpose of the current paper is to provide further details of our calculations as well as new results with improved numerical methods. In particular, we present an additional discussion of structure in critical behavior near transitions between collapse before and after reflection from the conformal boundary, evidence for self-similar (that is, chaotic) behaviour in the black hole formation time vs amplitude plots in transition regions, and data hinting at the formation of naked singularities in spacetimes with ADM mass below the algebraic bound for black hole formation.
The paper is organized as follows. In section 2, we review EGB gravity and derive via Hamiltonian techniques the relevant equations of motion in Schwarzschild coordinates. We also briefly describe our numerical methods there. We discuss our results on the above topics in section 3. We close with a summary and prospects for future work. An appendix contains the derivation of the equations of motion for the same system but in the AdS analogue of flat slice coordinates for future reference.

EGB Gravity and EoMs
In this section, we briefly review features of Einstein-Gauss-Bonnet gravity in AdS 5 and the Hamiltonian derivation of both the mass function and scalar equations of motion.

Einstein-Gauss-Bonnet Gravity In AdS
5-dimensional Einstein-Gauss-Bonnet gravity is a special case of Lanczos-Lovelock gravity [39,40]. The action is with λ > 0 in AdS (we use R for the Riemann tensor and its contractions). The covariant equations of motion are [41,42]:

4)
A key feature of spherically symmetric EGB is the existence of a generalized mass function where R is the areal radius and |DR| 2 = g µν R ,µ R ,ν [41]. In vacuum the mass function is constant on shell: ∂ µ M = 0 → M = M = constant. The most general vacuum solution with compact (positive curvature) horizon, given here in Schwarzschild-like coordinates, is where M is the on-shell value of the mass function and F 2 (R; M ) = |DR| 2 is obtained by inverting (2.5). The minus sign in front of the square root corresponds to the physical sector because it yields the Schwarzschild-Tangherlini-AdS solution in the limit that λ 3 → 0. Note that the GB term yields a modified effective cosmological constant as can be seen by taking either the M → 0 or R → ∞ limit in (2.6). The vacuum solution describes a single horizon black hole spacetime. In terms of the mass function, the horizon condition |DR| 2 R H = 0 is which implies that R H → 0 as M(R H ) → M crit ≡ λ 3 /2 even in the dynamical context. This suggests that it is impossible to form a black hole when the ADM mass is less than this critical value. This is a special feature of 5D EGB, as it depends critically on the exponent of R H in the third term of the mass function. It is similar to the existence of a critical mass for black holes in AdS 3 with Einstein gravity.

Hamiltonian Analysis
The total action describing the gravitational collapse of a massless scalar field in EGB gravity is The Hamiltonian analysis of spherically symmetric EGB without cosmological constant was performed in [43,44] and extended to generic Einstein-Lanczos-Lovelock gravity in [45]. Following [43,44] we use the ADM parametrization (2.10) and integrate out the angular coordinates in (2.9) to obtain a two dimensional dimensionally reduced action. With this metric, we define for future convenience. Here and in the following, a dot is the partial derivative with respect to t, and a prime is the partial with respect to the radial coordinate x. The dimensionally reduced Hamiltonian is a linear combination of the Hamiltonian constraint G and diffeomorphism constraint F where we have dropped an overal factor equal to the integral over the unit three sphere. As shown in [43], 14) The momentum conjugate to Λ is given by which determines y = y(R, Λ, P Λ ) implicitly in terms of the other gravitational phase space variables. Note that y is independent of P R and that we do not require the expression for P R in the following. By taking suitable linear combinations of the Hamiltonian and diffeomorphism constraints, the total Hamiltonian (2.12) can, up to boundary terms, be written [45] as where ρ m = 1 2 and M is the mass function (2.5) expressed in terms of phase space variables. It is important for the following that P R appears only in the diffeomorphism constraint F.
We choose as spatial coordinate R = R(x) with consistency conditionṘ = 0, which requires We can now set the diffeomorphism constraint, gauge fixing condition, and consistency condition strongly to zero to obtain the partially reduced Lagrangian where The remaining coordinate freedom can be fixed in two distinct ways. The first is to set the metric function Λ = Λ(x) to be a specific function of x. We outline this procedure in an appendix.
The more common choice, namely Schwarzschild-like coordinates, is used for numerical studies in much of the current literature. This class is obtained by the choice y(R, Λ, P Λ ) = 0, which in turn implies that P Λ = 0. When y = 0, (2.5) reduces to The consistency condition,ẏ = 0, for this gauge choice is Using the Hamiltonian constraint, We note that in vacuum (N Λ/R ) = 0 and the constraint M = 0 can be solved algebraically for Λ and N to give (2.6). The dynamical equations can be obtained by varying the following fully reduced Hamiltonian with respect to ψ and P ψ (Π and Φ are not canonical variables): In the above, N and Λ are implicitly defined by the consistency condition (2.22) and Hamiltonian constraint (2.24), respectively. They do not need to be varied, however, since the corresponding variations of H red are proportional to the Hamiltonian constraint and consistency condition. The resulting evolution equations arė The above, along with (2.21), (2.22), and (2.24) are the complete set of equations to solve. We now put these equations into the form used in [1] by making the substitutions We choose a compactified spatial coordinate R(x) = l tan(x/l) with l = 1/λ ef f . The metric in terms of dimensionless coordinates x → x/l, t → t/l is while the Hamiltonian constraint in terms of the new metric functions becomes where λ 2 = 1 − 2λ 3 /l 2 and we have used the identity We make the scalar field and its conjugate dimensionless by rescaling Φ → κ 5 Φ/ √ 3 and Π → κ 5 Π/ √ 3. Finally, we absorb l 2 into the mass function and λ 3 to make them dimensionless as well.
In the end, we solve the following equations: Since nonlinear self-gravitation effects drop off sufficiently quickly at large radius due to the dilution of energy density, the scalar field satisfies the same asymptotic expansion as in the linearized theory, Φ = ρ 3 Φ 0 + Φ 2 ρ 2 + · · · and Π = ρ 2 Π 0 + Π 2 ρ 2 + · · · , where ρ = π/2 − x. These are the boundary conditions of the normalizable linear eigenmodes e k (x), which can be defined in terms of Jacobi polynomials; the leading terms in these expansions correspond to expectation values of operators in the dual field theory. 3 At the origin, we require that Π be an even function of x and Φ be odd for smoothness.

Numerical Methods
We briefly outline our numerical methods and how these have changed since our previous work [1]. A detailed description is provided in an appendix of [35]. The key improvement to our code is that we now solve the spatial ordinary differential equations using an adaptive fifth-order Dormand-Prince stepper. We set the desired relative and absolute tolerances and the stepper will adjust the step size over the spatial mesh so that the desired tolerances are met locally. The adaptive method requires scalar field data in between grid points which we supply using a cubic spline. We find that the stepper takes many small steps near the origin and much larger steps further out.

Results
As in [1], we consider initial data of the form ie, Gaussian in Π, and a GB parameter of λ = 0.002. Figure 1 provides an overview of our results for the horizon formation time t H , which cover an amplitude range = 27 − 48. In the figure, blue circles represent formation of a horizon, while red triangles represent lower limits on t H for amplitudes which do not form a horizon for t < 100. For Einstein gravity (λ = 0), t H would be approximately piecewise constant appearing as "steps" with t H with decreasing amplitude. Physically speaking, at large amplitude, gravitational collapse can proceed immediately, but lower amplitude initial data disperses, reflects from the conformal boundary one or more times, and finally collapses after more gravitational focusing. As in the earlier results of [1], gravitational collapse in EGB gravity exhibits the same as well as additional features. First, there is a transition from immediate collapse to collapse after one or more reflections. There is critical behavior at these transitions, which has been studied in some detail in [46,47] for Einstein gravity in AdS. In the following subsection, we study the first critical point at large amplitude, when horizon formation stops occurring immediately, extending the analysis of this region in [1].
Another key feature of figure 1 is that t H appears to demonstrate sensitivity to initial conditions in certain amplitude ranges. That is, while there are some steps in horizon formation time where t H remains approximately constant with , t H varies wildly in transition regions between the steps. At low enough amplitude, the steps are apparently so narrow that they dissolve into the transition regions. In subsection 3.2, we explore in more detail whether the transition regions exhibit chaotic behavior such as self-similarity.
Because horizon formation is apparently sensitive to initial conditions in some regions of the amplitude, we have opted to keep all the simulations of figure 1 at a fixed resolution of 2 12 + 1 grid points, even when they begin to lose convergence (as illustrated by a loss of conserved ADM mass). Otherwise, an increase in resolution could act as a small shift in amplitude. At this resolution, simulations lose up to 2.5% of the conserved ADM mass by t = 100, so simulations that do not form a horizon by this time are shown only as lower limits on t H . We have tested several amplitudes with 2 13 + 1 grid points and found that subcritical simulations remain subcritical while horizon formation times in the step regions (which are stable vs change of initial conditions) have a relative difference of 5 × 10 −7 .
Finally, at the lowest amplitudes shown, none of the simulations form a horizon. As noted earlier, horizons cannot form below a critical conserved mass M crit in EGB gravity. In other words, all initial data must be stable at low amplitudes, in apparent contrast to the case of Einstein gravity. For our choice of initial data, the critical mass corresponds to an amplitude of approximately crit ∼ 21.86; figure 1 hints that higher amplitudes may also be dynamically stable against horizon formation. It is also an interesting question whether evolution of initial data below the critical amplitude is quasi-periodic or evolves toward a naked singularity. In section 3.3, we study the evolution for two amplitudes, one just larger than and one just smaller than the critical amplitude, and present evidence suggestive of naked singularity formation at finite origin time below the critical amplitude.

Critical Phenomenon
Critical phenomena in the gravitational collapse of a spherically symmetric scalar field in Einstein gravity (for 4-dimensional asymptotically flat spacetime) was first observed by Choptuik [2]. Choptuik found that geometrical quantities such as the mass of the black hole obey the scaling law where p is a parameter in the initial data profile, p is the critical value of p, 4 , and γ is the critical exponent. A detailed semi-analytic study by Gundlach in four dimensions [48] found that γ = 0.374 ± 0.001. For small amplitude initial data in asymptotically AdS spacetime, any horizon that forms will be small compared to the AdS scale, so the critical behavior at any transition (ie, collapse after n reflections transitioning to collapse after n + 1 reflections) should have the same critical exponent as the asymptotically flat case. In the case of Einstein gravity, [47] confirms the expectation, finding a critical exponent consistent with the Gundlach value independent of the width of initial data or the number of reflections before collapse. The critical behavior of EGB gravity differs from Einstein gravity even in asymptotically flat spacetime. For one, the Gauss-Bonnet term contributes to the equations of motion only in 5 dimensions or more; in 5D Einstein gravity, the critical exponent is γ ≈ 0.416 [49,50]. Critical phenomena in 5D EGB gravity has been studied in [3,4], which found that the new, short distance length scale alters the near-critical behavior such that no black hole forms below a minimum horizon radius [3]. This is similar to the case of a massive scalar field in asymptotically flat spacetime, which also has a dynamically determined minimum horizon radius [51].
Again, it is natural to ask which features of the critical behavior persist or differ in asymptotically AdS 5 spacetime. As in 4D, we expect critical behavior near each transition (n to n + 1 reflections) to match that in asymptotically flat spacetime because the black holes formed are initially much smaller than the AdS curvature scale. In figure 2, we show log(r AH ) as a function of log(( − )/ ), where is the amplitude above which scalar field configurations collapse immediately (the 0 to 1 reflection transition). We find a critical exponent γ ≈ 0.42 in agreement with results in asymptotically flat spacetime.
For EGB gravity, we expect a minimum horizon radius at each critical point, as in asymptotically flat spacetime. Figure 3 shows the scaling of the initial apparent horizon radius with amplitude near the critical point for immediate collapse, ≈ 45.3315. It is initially apparent from figure 3a, which shows values of far from where the black holes form very quickly, that there is in fact a radius gap as seen in [1]. Continuing to amplitudes with − 10 −5 in figure 3b, we observe persistence of the radius gap R min ∼ 10 −1.9 along with sudden jumps, or steps, in the horizon radius.  An explanation for this behavior is apparent in animations of our simulations. As the scalar field collapses, the initial profile fragments, with the main portion of the mass remaining near the origin and driving horizon formation while several pulses of mass disperse toward the boundary. For − small, one or even two of these subsidiary pulses can reflect from the boundary and return to the neighborhood of the origin (possibly multiple times) before A(t, x) reaches the threshold for approximate horizon formation. These subsidiary pulses are responsible for the multiple local minima of A(t, x) noted in [1]. Animations showing subpulses reflecting from the boundary once and twice are available at http://ion.uwinnipeg. ca/~afrey/AdSGB16.html. Although the horizon formation times t H for these amplitudes are small, it is important to remember that t is the proper time at the origin. As it turns out, there is a significant redshift factor between this time and the proper time outside x r AH . To quantify the time dilation factor, in figure 4 we plot the proper time of an observer at the AdS boundary, given by as a function of the proper time at the origin (in one particular collapse). While an insufficient amount of time apparently passes for the scalar field to reflect off of the AdS boundary according to observers at the origin, the relevant time is actually better approximated by the proper time at the AdS boundary since δ(t, x) is roughly spatially constant outside the main portion of matter. Specifically, while only a time ∆t ≈ 0.37 passes, the corresponding boundary time elapsed is ∆τ ≈ 4.4, enough for the subsidiary pulses to reflect off the boundary and interact with the forming black hole. Interestingly, this effect should be observable in Einstein gravity close enough to the critical amplitude (since infinite boundary time passes before the metric function A(t, x) → 0), but it appears to be much more challenging to resolve. Some progress on the subject has been made [52], but a different gauge choice and black hole excision may be needed to fully explore this behavior. The GB term seems to enhance time dilation effects significantly.

Self-Similarity
As we noted above, the initial horizon formation time t H exhibits a much richer structure in EGB gravity than in pure Einstein gravity. By now it is well-known that t H typically increases piecewise monotonically with decreasing for massless scalar matter in Einstein gravity, forming the ubiquitous step structure seen in many references. 5 In contrast, while figure 1 also exhibits some steps in t H , the transitions from step to step exhibit a significant non-monotonic scatter in t H . For example, while initial data with 45.3 collapses immediately and initial data with 39. 6 44.0 collapse after one reflection from the boundary, amplitudes 44.0 45.3 vary wildly. The appearance of smaller steps and apparent sensitivity to initial conditions in the transition regions led [1] to speculate that the t H vs curve may have a fractal structure. Here, we investigate the transition region 44.0 ≤ ≤ 45.3 in more detail with the aim of uncovering signatures of chaotic behavior. Since changing resolution amounts to a change in initial conditions, all simulations discussed in this subsection are carried out at a fixed resolution of 2 14 + 1 grid points, following the same reasoning explained above. Figure 5 shows the transition region in detail for three ranges; figure 5a shows the entire region, figure 5b shows a small area surrounding the t H ≈ 132 point at = 44.4, and 5c shows a smaller area to the right of that point. Blue circles represent horizon formation, while red triangles represent simulations that do not form a horizon within t = 305, which can be taken as a lower limit on t H for those amplitudes. These simulations lose several percent of the ADM mass by that time, however, so a conservative reader may prefer to read these as lower limits of t H 170, just greater than the largest values of t H for collapsing simulations. Regardless, the plots for the three amplitude ranges show a similar structure of rapidly varying horizon formation times with amplitude. This remains suggestive of fractal-like, self-similar behavior, at least on the scales shown.
To test the self-similarity of the t H vs curve quantitatively, we use a variation on the box-counting-dimension estimate. Specifically, we draw grid lines at each of the tick marks on figure 5a and count the number of boxes so created that are occupied by data points. For a first estimate, we include subcritical simulations as if they have t H = 300. Data points lying on a grid line are counted as occupying the box above or to the right as appropriate. In this case, a curve of dimension D should occupy N = W/s D boxes, where W is the total horizontal range (W = 1.3 in figure 5a) and s is the length between grid lines (s = 0.26 in figure 5a); the box-counting dimension is defined as D = − lim s→0 ln(N/W )/ ln s. To take an approximate limit, we repeat the procedure for figures 5b and 5c, keeping vertical grid lines at the tick marks shown but scaling the vertical distance between horizontal grid lines with s. We find N = 9, 10, 11 and consequently D s = 1.44, 1.22, 1.14 respectively for the three subfigures (s = 0.26, 0.04, 0.004). However, it is reasonable to argue that the amplitudes that do not form a horizon may have different values of t H from each other (or be truly stable), so we should not count them. If we repeat the box-counting test while ignoring the apparently  stable points and also subtracting from W the width of any boxes that contain no collapsing data points, we find D s = 1.35, 1.15, 1.13 for the three subfigures. This provides weak but suggestive evidence that the t H vs curve has a fractal dimension of around 1.14, somewhat greater than unity.
Another characteristic of chaotic behavior is exponential growth of some measure of distance between two systems with similar initial conditions, |∆| ∼ exp(λt) for Lyapunov exponent λ. We consider three neighboring amplitudes in figure 5c, 1,2,3 = 44.413, 44.412, 44.411 and take as a measure of the distance between them the difference in the upper envelope of the Ricci scalar at the origin ∆ 12 =R 1 (t) −R 2 (t), etc. The bar indicates the maximum (at the origin) over one full reflection from the conformal boundary ∆t = π. Figure 6a shows thē R for the three amplitudes (in green dashed, blue solid, and red dotted curves); note that the three amplitudes lead to different values of t H , so theR 1,2 curves do not extend across the  entire plot. Figure 6b shows ∆ 12 (blue points) and ∆ 23 (magenta squares) versus time along with best fit exponential functions (blue solid and magenta dashed lines respectively). Both fits are consistent with Lyapunov exponents λ ∼ 0.31, or mild chaotic behavior. The actual differences appear to have oscillation superimposed on the exponential growth. While [1] first suggested that AdS gravitational collapse in EGB gravity exhibits chaotic behavior, [53] have also found evidence for chaos in the gravitational collapse of two thin shells of matter in AdS with Einstein gravity. As in our figure 5a, [53] finds hints for selfsimilarity in the t H curve as a function of initial conditions (in their case, the common initial radius of the two shells of matter). In this system, energy transfers gravitationally between the shells as they pass through each other; in the self-similar region, the transfer back and forth leads to chaotic behavior in the horizon formation time. The two shells are not both near the origin when the horizon forms; instead, the horizon forms when one of the shells happens to have accumulated a large enough density to form a horizon on its own. In addition, [53] also finds a small but positive Lyapunov exponent for the deviations between nearby initial conditions in the chaotic region of parameter space. Clearly this is similar behavior, and there may be a deeper analogy between scalar collapse in EGB gravity and the two-shell system. Specifically, at least for some amplitudes, the GB term causes the initial scalar field pulse to break into multiple pulses, each of which behaves as an independent shell of matter. For shells with large radii, the GB term is negligible, so we are in fact also studying the collapse of multiple transparent shells in Einstein gravity. Examining one of our evolutions as an animation is instructive; an animation of M for = 44.412 is available at http://ion.uwinnipeg.ca/~afrey/AdSGB16.html. We see that the initial pulse slowly separates into two (groups of) pulses of matter, which are approximately completely out of phase by t ∼ 15 and each of which contains one tall, thin shell of matter. Eventually, one of the pulses forms a horizon while the other is far away. So, once the GB term separates the initial matter distribution into two pulses, it seems that energy transfer between pulses may be responsible for the chaotic behavior, as in the two-shell system. We have also examined our simulations in the piecewise-constant regions of figure 1 for comparison; while collapses that reflect from the boundary multiple times do exhibit some pulse fragmentation, only the main pulse ever has a high, thin shell of matter.

Naked Singularity Formation
In EGB gravity in AdS 5 , horizons must contain at least a minimal mass (even at zero radius); since the asymptotic value of the mass function is conserved, this implies that horizons cannot form below a critical value of the amplitude. For our initial data and choice of GB parameter, crit ∼ 21.86. We have already noted that we have failed to find horizon formation for t < 100 for any amplitude ≤ 32, leaving several important questions. One, is there a dynamical mechanism preventing gravitational collapse for small amplitudes that are nontheless greater than crit ? We can attempt to answer this by studying these amplitudes with high-resolution simulations for long times. For another, do amplitudes that do not form horizons lead to a stable, quasiperiodic evolution, or can they form naked singularities? Is the behavior the same or different for amplitudes above and below crit ? To address these questions, we have carried out simulations at = 20 and = 22, increasing resolution as necessary to carry the simulations to as long a time as possible.
We have been unable to find horizon formation in either case to times of t ∼ 325, 295 and resolutions up to n = 18, 19 respectively for = 20, 22. The need for the high resolutions is clear when we considerR, the upper envelope of the Ricci scalar at the origin, which we show in figures 7,8. In both cases, we find strong growth of the Ricci scalar to very large values, eventually reaching values of orderR ∼ 10 7 while avoiding formation of a horizon. From visual inspection of the simulations, the key physics seems to be dispersal of the original matter pulse into two pulses, which individually narrow, leading to very high curvatures, but which are nonetheless not massive enough to form a horizon. Nonetheless, the extreme growth ofR and pulse narrowing (which also drives the need for increasingly higher resolution at late times) suggests the possibility that these amplitudes will eventually form a naked singularity, rather than behaving in a quasi-periodic fashion.
As further suggestive evidence of singularity formation, we have studied the late-time energy spectra of both evolutions. Figure 9 shows the energy spectra (to the j = 1024 eigenmode) for both amplitudes at the latest time we were able to simulate in each case. These show a slow power-law decay at large mode number, which is usually characteristic of horizon formation. 6 However, it is impossible for a horizon to form for = 20, so the distribution of energy through the higher modes suggests the possible development of a naked singularity. Another point suggestive of singularity formation is that we find over 1% of the total energy lies in higher modes (j > 1024) for times greater than t ∼ 322.4. Similarly, the = 22 evolution seems to be moving rapidly toward either horizon or naked singularity    formation for t ∼ 295. Over 1% of the total energy is in higher modes for t 287, and close to 3% is in higher modes by the end of our simulation at t ∼ 295. This degree of energy in high eigenmodes allows an extreme concentration of energy density near the origin, which could drive the formation of a singularity. It is important to note that these spectra differ from that at earlier times (see for example the supplemental information of [1]), which have an apparent exponential cut-off for mode numbers less than 10 3 , indicating that the power has continued to cascade into higher and higher modes as the evolution progressed. This difference is one reason for a potentially different conclusion about evolution at amplitudes near crit in comparison to [1].  Assuming that a naked singularity does form at a finite time t, which is the proper time at the origin, it is important to ask whether redshift effects push the singularity formation to an infinite conformal time at the boundary, which controls the physics in the dual CFT. Unlike the case of horizon formation, however, time dilation effects seem to be unimportant here; at the latest times probed by our simulations, the metric function δ takes values of −0.007 and −0.035 at the boundary for = 20, 22 respectively, which are the minima over their evolution.
It is worth considering at this point what the formation of a naked singularity at finite time would mean in the context of the AdS/CFT correspondence. If the gravity side of the correspondence is described by the pure EGB gravity with no additional higher-curvature terms, the translation of a naked singularity to the dual field theory is unclear. It is tempting, therefore, to postulate that a naked singularity is a sign of a pathology in the theory; in fact, [55] has already argued that theories dual to pure EGB gravity in AdS suffer from acausalities. On the other hand, if the Gauss-Bonnet term is just the first in a tower of higher-curvature corrections, the extreme curvatures of figures 7,8 suggest that the additional corrections will become significant (as seems to be the case in the self-similar transition regions of section 3.2 as well). That would be a signal that the effective gravity theory is breaking down and should be replaced by the full string theory, and it is possible that the end state of collapse is a gas of strings near the origin. 7 Due to the very high resolutions necessary, it was not computationally feasible to perform convergence testing for the entire simulations shown in figures 7,8. However, convergence tests for part of the = 20, n = 17 calculations showed the expected 4th-order convergence for the Ricci scalar at the origin, and the clear overlap of the different resolutions in much of the figures argues for the reliability of our results.

Discussion
In this manuscript, we have expanded on the analysis of black hole formation in AdS Einstein-Gauss-Bonnet gravity first presented in [1]. In our examination of the horizon formation time t H as a function of amplitude (figure 1), we have considered three particular physical phenomena in detail: critical behavior at a transition in t H , possibly chaotic behavior below the transition amplitude, and long-time evolution at low amplitudes possibly hinting at formation of naked singularities.
We first examined critical behavior at the transition from immediate collapse, i.e. when a black hole forms without the matter first dispersing, to collapse after one or more reflections from the conformal boundary and confirmed the existence of a dynamical radius gap of approximately R min ∼ 10 −1.9 (see figure 3). A number of questions remain about this critical behavior. For example, do the transitions from other step regions of figure 1 exhibit the same radius gap as in figure 3, or does the value change (or even vanish)? Is there a universal scaling law for the horizon radius for amplitudes below any transitions in figure 1, as was demonstrated in Einstein gravity in [46,47]? We expect that the radius gap is due to the existence of a massive critical solution, in which case the local features of the critical behaviour should be independent of the number of reflections from the boundary before horizon formation. Thus, the radius gap should persist for higher numbers of bounces. Some evidence in this direction was provided in figure 3 of [1], which compares the scaling plots after one bounce to those with no bounces.
Our simulations also shed light on a novel dynamical feature of the critical behavior, first presented in [1]. The step-like behaviour in the scaling plot (figure 3b) a time dilation effect; part of the initial lump of matter disperses from the origin rather than falling into the forming horizon. For amplitudes close enough to the transition value, one or more of these sub-pulses have enough time to reflect from the AdS boundary and return to the origin before the simulation reaches our criterion for horizon formation. Although not previously observed, this could in principle occur in any gravity theory in AdS that exhibits critical behaviour, including Einstein gravity. It would be interesting to check since this effect depends both on global features of AdS and local dynamics.
The transitions from one piecewise-constant region of figure 1 to another also demonstrate significant scatter. We have presented evidence that the t H vs curve is self-similar in the region from ∼ 44.0 to 45.3. Furthermore, a positive Lyapunov exponent between the evolutions of nearby amplitudes is a hallmark of truly chaotic behavior. The chaotic behaviour appears to have as its source the separation of the initial pulse into two (or possibly more) pulses and the subsequent transfer of energy between them as they pass through each other between reflections off infinity and the origin. Here, too, questions remain: Is chaos present for any gravitational system in which the matter forms multiple pulses, as is the case in the two-shell system of [53]? Does the matter distribution fragment whenever the physics has a second scale (other than the AdS scale), such as a mass for the scalar field?
Perhaps the most intriguing aspect of our analysis is the potential evidence for naked singularity formation in the model. Below the algebraic mass gap of M crit = λ 3 /2 ( crit = 21.86 for our GB parameter and initial data), there are only two possible end states: a quasi-periodic steady state or naked singularity formation. While it is impossible to provide definitive proof numerically, the observed dramatic increase in curvature and concentration of energy into higher modes in the absence of horizon formation suggest that the end state will be a naked singularity for our two long evolutions near the critical value (one just below, the other just above). In pure EGB gravity, naked singularity formation may indicate a pathology of any dual field theory. On the other hand, the extreme curvatures found may instead indicate the excitation of other higher curvature terms and string degrees of freedom, leading to the eventual production of a gas of strings rather than the actual development of a singularity. It is curious that horizon formation seems strongly suppressed (takes a lot longer, or does not form at all) even above the algebraic threshold crit . Our results allow us to speculate as to the cause of this suppression, which seems to be a highly nonlinear effect. In the case of our evolution just above threshold, = 22, the ADM mass is M ∼ 0.00101, which is just barely above the critical value (see figure 10a for the conserved mass as a function of amplitude). Assuming the existence of a dynamical radius gap of about R min ∼ 10 −1.9 , equation (2.8) implies that the minimum amount of mass dynamically required to form a horizon is actually close to M(R min ) ∼ 0.00108, so both evolutions discussed in section 3.3 were below this dynamical limit. The fact that we do not see a black hole form slightly above the critical value is perhaps not a surprise. More surprising is the apparent suppression of black hole formation for amplitudes near = 32 and below (see figure 1). The mass at this amplitude is close to 0.002, double M crit . However, we have seen in chaotic transition regions that the mass tends to split into at least two pulses. In this case, the splitting can produce smaller shells of matter that individually do not have enough energy to form a horizon, so that horizon formation depends on the subsequent energy transfer between pulses/shells. Indeed, this splitting occurs at long times for = 30.2, as indicated in figure 10b, so it seems that the lowest amplitudes shown in figure 1 are in a chaotic region. Since the two pulses appear to carry a substantial fraction of the mass, horizon formation will require significant energy transfer between pulses (more than in chaotic regions at higher amplitude). It seems likely that this significant transfer is unlikely and will occur only rarely, leading to very long horizon formation times. It is also worth speculating if similar physics occurs for gravitational collapse in AdS 3 , which also has a critical black hole mass and apparent suppression of horizon formation for amplitudes just above the critical mass [56,57].
It is clear that gravitational collapse in EGB gravity in AdS is an intricate system, showing first-order transitions, chaotic behavior, and possible formation of naked singularities.  A Appendix: Generalized Flat Slice (PG) Coordinates Working as before with R = R(x), we choose: Flat slice coordinates would correspond to the choice Λ = 1, but this is not appropriate when the cosmological constant is non-zero [58]. Instead we first write the mass function as and choose a function Λ(x) that yields a diagonal metric in vacuum. This requires the first three terms in the square brackets above to vanish, which in turn implies that The sign in front of the square root is chosen to give the usual answer Λ = R when λ = 0. When λ 3 = 0 and R = x, this gives: which is correct for AdS in Schwarzschild coordinates for M = 0. With this choice the mass function reduces to: The Hamiltonian constraint is which determines y = −(N r /N )R in terms of data on a slice. In these coordinates the consistency conditionΛ = 0 is found to be where using (A.5) The dynamical equations in terms of Π = P ψ /R n−2 and Φ = ψ arė (A.10)