{'fraction_non_alphanumeric': 0.03610959898385048, 'fraction_numerical': 0.010291106099489333, 'mean_word_length': 4.56199538638985, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalarized Black Hole dynamics in Einstein dilaton Gauss-Bonnet Gravity

We report on a numerical investigation of the stability of scalarized black holes in Einstein dilaton Gauss-Bonnet (EdGB) gravity in the full dynamical theory, though restricted to spherical symmetry. We find evidence that for sufficiently small curvature-couplings the resulting scalarized black hole solutions are nonlinearly stable. For such small couplings, we show that an elliptic region forms inside these EdGB black hole spacetimes (prior to any curvature singularity), and give evidence that this region remains censored from asymptotic view. However, for coupling values"superextremal"relative to a given black hole mass, an elliptic region forms exterior to the horizon, implying the exterior Cauchy problem is ill-posed in this regime.


I. INTRODUCTION
We present numerical results on the nonlinear evolution of spherically symmetric black hole solutions in a modified gravity theory: Einstein dilaton Gauss-Bonnet (EdGB) gravity. EdGB gravity is one of the few known scalar tensor theories that may admit sensible classical evolution (at least for some open subsets of initial data; see our earlier studies [1,2]), yet does not allow the Schwarzschild or Kerr stationary black hole solutions. Instead the expected solutions are conjectured to be "scalarized" black holes [3][4][5] (the detailed form of this statement depends on the functional form of the Gauss-Bonnet coupling and scalar field potential [6,7]). The variant of EdGB gravity we consider is poorly constrained by weak field gravity measurements (e.g. from binary pulsars [8]), though to be consistent with the speed of gravitational waves inferred by the binary neutron star merger GW170817 [9] requires a negligible cosmological background for the dilaton field [10]. Assuming the latter, the strongest constraints on the theory may then come from gravitational wave observation of the final moments of binary black hole inspiral. The theory thus provides an interesting alternative to general relativity (GR) to compare against when gravity is in the strong field dynamical regime (see e.g. [11] and references therein).
From the perspective of effective field theory, EdGB gravity can be motivated as the leading correction to GR in a low energy expansion of quantum gravity that incorporates mixing between a scalar degree of freedom and the tensor (metric) degrees of freedom of GR [12,13]. Then one would not expect significant modifications to GR away from the Planck scale, and in particular not for astrophysical black holes. Alternatively, as we do here, one could consider the coupling parameter of the theory to be arbitrary, and EdGB gravity taken verbatim as a classical theory of gravity with a scale dependent modification to GR. Such a theory may not be mathematically well-posed in some regimes (or in a generic sense not at all [14,15]), though a "healthy" sector of solutions could still be extracted by treating it as an effective field theory and limiting consideration to small perturbative corrections to GR. Several groups are pursuing such an "order reduction" approach to understanding EdGB gravity, and related theories where beyond-Ricci curvature scalars are added to the Einstein Hilbert gravitational action [11,[16][17][18][19][20][21][22]. Another approach, inspired by the Israel-Stewart "fix" of relativistic hydrodynamics [23], is to explicitly modify the GR-extensions to lead to well-posed equations [24,25].
Our approach is instead to attempt to solve the complete classical field theory, and discover which classes of initial data (if any) lead to well-posed evolution. Our motivation is two fold. The first is the desire to know how classical gravity can in principle differ from the predictions of GR in the dynamical strong field regime, as is applicable to the last stages of binary black hole coalescence. This could give more meaning to quantitative statements of consistency of observed waveforms with the predictions of GR, help constrain EdGB gravity, or discover modifications to GR of a class similar to that offered by EdGB gravity. Though an effective field theory approach as described above is likely "guaranteed", by construction, to give well-posed evolution schemes for small deviations to GR, it is still unknown if this approach could be pushed to solve for modifications large enough to provide waveforms distinguishable from GR in an observation, given the typical signal-to-noise ratios expected from the current generation of ground based detectors, and taking waveform degeneracies into account.
The second reason is that nonlinear modifications of gravity have been introduced in attempts to address the discovery of dark matter, dark energy, solve the flatness and horizon problems of early universe cosmology, and resolve the issues of the initial cosmological singularity and singularities formed during gravitational collapse. EdGB is an important and representative member of a class of modified gravity theories, Horndeski theories, that have been invoked to solve these various issues within classical GR (for a recent review on Horndeski theories see [26]), and it thus is interesting to see if the theory is mathematically sensible as a classical field theory.
Previous studies of EdGB black holes have focused on static solutions to the theory (e.g. [5][6][7]), the dynamics of the EdGB scalar in the decoupling limit (e.g. [16,19,27]), or linear perturbations of static EdGB black hole backgrounds (e.g. [28]). Additionally, a recent work [21] explored the dynamics of the scalar and metric fields of the theory from an effective-field theory framework. There it was shown that scalarized, rotating EdGB black hole solutions are stable for small enough couplings, up to leading order in metric and scalar field perturbations. Restricted to spherical symmetry, our results extend this to all orders in the Gauss-Bonnet coupling, showing consistency for small couplings, and showing where and how the theory breaks down for large couplings.
In a previous work, we studied the dynamics of EdGB gravity in spherically symmetric collapse using horizonavoiding coordinates [1,2]. There we considered collapse of a dense concentration of the dilaton field, and for sufficient mass could show a horizon was forming, though we could not evolve beyond that time. Here, we give results from a new code solving the EdGB equations in horizon penetrating coordinates, allowing us to investigate the long-term stability of scalarized black holes (for times up to of order t ∼ 10 3 m, where m is the mass of the black hole). Also, we begin with the Schwarzschild solution as initial data (with an optional exterior dilaton field perturbation).
Upon evolution of Schwarzschild initial data, we find that the scalar hair grows, and an elliptic region forms in the interior of the black hole. This indicates black hole physics in EdGB gravity has aspects of it governed by a mixed elliptic-hyperbolic equation (or simply mixed-type equation), and it is unclear how this could affect the Cauchy problem exterior to the horizon. I.e., there is no a-priori reason to expect this elliptic region to "obey" cosmic censorship, and leave the scalar in the exterior domain to be governed by a hyperbolic partial differential equation (PDE). Instead, we will simply assume that this is possible, and excise the elliptic region from the domain. If during subsequent evolution no new elliptic region forms, and the solution settles to a stationary state, we will claim this is a self-consistent application of excision, and the resulting hairy black hole is stable (to within limitations of numerical evolution) 1 . We do find this to be the case for Gauss-Bonnet couplings below an extremal limit for a given black hole mass. We compare these solutions to the scalarized decoupled black hole solutions of EdGB gravity, and find good agreement, the better the smaller the Gauss-Bonnet coupling is (for a fixed black hole mass). However, above the extremal limit, an elliptic region does form outside the horizon, indicating a break-down of the exterior Cauchy problem for small black holes (relative to the EdGB coupling scale).
An outline of the rest of the paper is as follows. In Sec. II we describe the equations of motion, variables, and metric ansatz we use. In Sec. III we describe aspects of the numerical code, including our excision strategy, as well as some diagnostic quantities we monitor. In Sec.IV we describe the results mentioned above in detail, and end in Sec. V with concluding remarks. We leave some convergence results, and a derivation of the decoupling limit about a Schwarzschild black hole in Painlevé-Gullstrand, to the appendices. We use geometric units (8πG = 1, c = 1) and follow the conventions of Misner, Thorne, and Wheeler [30].

II. EQUATIONS OF MOTION
The action for the EdGB gravity theory we consider is where f (φ) is a (so far unspecified) function, and G is the Gauss-Bonnet scalar where δ µναβ ρσγδ is the generalized Kronecker delta. In geometric units, the Gauss-Bonnet coupling constant λ has dimensions length squared. Varying (1) with respect to the metric and scalar fields, the EdGB equations of motion are In this work we will only consider the coupling function While other coupling functions are often considered in the literature on EdGB black holes, this is the simplest which is thought to give rise to stable scalarized black hole solution; see [3,5,16,19,31] and references therein. This coupling may additionally be motivated as the lowest order term in the effective field theory expansion of a metric theory coupled to a scalar field (e.g. [8]). From the symmetry λ → −λ, φ → −φ, we only consider λ ≥ 0. We evolve this system in Painlevé-Gullstrand (PG)-like coordinates (e.g. [32][33][34][35]) so-named since t = const. cross sections are spatially flat (the Schwarzschild black hole in these coordinates is given by α = 1, ζ = 2m/r). We define the variables and take algebraic combinations of Eq. (3b) and the tr, rr, and ϑϑ components of Eq. (3a) (c.f. [2]) to obtain the following evolution equation for the {P, Q} variables: The quantities A (P ) and F (P ) are lengthy expressions of {α, ζ, P, Q} and their radial derivatives. In the limit λ = 0 Eq. (7b) reduces to Interestingly, in PG coordinates the Hamiltonian and momentum constraints do not change their character as elliptic differential equations going from GR to EdGB gravity: where ρ ≡ n µ n ν T µν = 1 2 f ≡ df /dφ, and n µ ≡ (−α, 0, 0, 0). While Eqs. (9a) and (9b) hold for any f , as mentioned above we only consider f (φ) = φ in this article.

A. Diagnostics
As PG coordinates are spatially flat the ADM mass prescription always evaluates to zero, and does not capture the correct physical mass of the spacetime. Instead then we use the Misner-Sharp mass [36] evaluated at spatial infinity to define the spacetime mass (and this does give the correct mass for the Schwarzschild solution in GR) As described in the next section we compactify so that infinity is at a finite location on our computational grid. The asymptotic mass m is preserved up to truncation error in our simulations of both EdGB gravity and GR. Given a spacetime with mass m, we define the dimensionless curvature-coupling We will classify different solutions based on their curvature couplings C, with GR the limit C = 0; empirically (as we discuss in our results IV) we find strong EdGB corrections arising when C 0.1.
Following the procedure used in [1,2], we calculate the radial characteristics of the scalar degree of freedom via Eqs. (7a) and (7b), after having removed the spatial derivatives of α and ζ from these equations using the constraints (9a, 9b). The corresponding characteristic speeds c ± are where ξ a ≡ (ξ t , ξ r ) solves the characteristic equation In the limit λ = 0, these speeds reduce to the radial null characteristic speeds in PG coordinates c We see that ζ = 1 marks the location of a marginally outer trapped surface (MOTS) (e.g. [37] and references therein). We take the location of the MOTS to represent the size of the black hole on any given time slice. The characteristic equation, Eq. (17), takes the following form when expressed as an equation for the characteristic speeds c where Where the discriminant D ≡ B 2 − 4AC > 0 are regions of spacetime where the equations are hyperbolic, where D < 0 the equations are elliptic, and following the language of mixed-type PDEs (e.g. [38] and references therein), the co-dimension one surfaces where D = 0 separating elliptic and hyperbolic regions are called sonic lines. In the GR limit λ = 0 the scalar equations are always hyperbolic (D > 0), though as we found in [1,2], for sufficiently strong couplings C the discriminant D is not of definite sign, and the scalar equations are then of mixed-type (similar conclusions have been drawn for other member of the Horndeski class of theories; see e.g. [39][40][41]).

B. Spatial compactification
In vacuum when P = Q = 0, the general solution to Eqs. (9a) and (9b) is ζ ∝ r −1/2 and α = const. 2 . We found that this falloff in ζ made it difficult to impose stable outer boundary conditions at a fixed, finite r. To alleviate this problem, we spatially compactify through a stereographic projection where L is a constant, and discretize along a uniform grid in x, with spatial infinity x = L now the outer boundary of our computational domain. For all the simulations presented in this article we chose L = 5m, where m is the mass of the initial Schwarzschild black hole. At x = L we impose the conditions α| x=L = 1 , ζ| x=L = 0, P | x=L = 0, Q| x=L = 0, φ| x=L = 0. These conditions are consistent with our initial conditions and asymptotic fall off of the metric and scalar field. For the latter, generally φ → 1/r, though if we impose exact Schwarzschild initial data outside some radius r 1 (such that φ(r > r 1 , t = 0) = 0 and ∂ t φ(r > r 1 , t = 0) = 0), the Gauss-Bonnet curvature will source an asymptotic field that decays like 1/r 4 ; by causality (as long as the equations are hyperbolic) the 1/r component sourced by the BH, or any scalar radiation from a field we put in at r < r 1 , will never reach spatial infinity.

C. Initial data
The computational domain covers , where x exc (r exc ) is the excision radius, and can vary with time (described in the following section). We set initial data at t = 0 by specifying the values of P and Q, and then solve for α and ζ using the momentum and Hamiltonian constraints. These ordinary differential equations (ODEs) are discretized using the trapezoid rule and solved with a Newton relaxation method, integrating from x = x exc to x = L. At x = x exc (some distance inside the horizon, as discussed in the next section) we set α and ζ to their Schwarzschild values: If we begin with zero scalar field energy (Q| t=0 = P | t=0 = 0), solving the constraints recovers the Schwarzschild solution on t = 0 to within truncation error.

D. Excision
At every time step we solve for α, ζ, P , and Q by alternating between an iterative Crank-Nicolson solver for P and Q and the ODE solvers for α and ζ, until the discrete infinity norm of all the residuals are below a pre-defined tolerance (this is a similar strategy to that used in our earlier code based on Schwarzschild-like coordinates, and more details can be found in [2]). The excision strategy assumes all characteristics of hyperbolic equations are pointing out of the domain at the excision surface x = x exc . This implies that for P and Q we cannot set boundary conditions there, rather their evolution equations must be solved, with the finite difference stencils for the radial derivatives appropriately changed to one-sided differences. For α, as with the initial data, the inner boundary condition is arbitrary, and after each iteration we rescale it so that α(t, x = L) = 1. For ζ, to obtain a consistent solution to the full field equations requires that the boundary condition ζ(t, x = x exc ) be set by solving the corresponding evolution equation for ζ there (optionally ζ could be evolved over the entire domain using this equation). Taking algebraic combinations of the equations of motion (3a), an appropriate evolution equation for ζ can be obtained The expressions for A (P ) and F (P ) are lengthy expressions of λ, {α, ζ, P, Q}, and their radial derivatives. In the limit λ = 0, Eq. (23) reduces to FIG. 1: Schematic illustration of the solution of ODEs at synchronized levels with fixed mesh refinement in our setup: we integrate from left to right (the compactified radial coordinate x increases from left to right).
Eq. (23) provides the boundary condition for ζ at the excision surface. We then integrate outwards in r using the Hamiltonian constraint, Eq. (9a) as described above to solve for ζ. The formation of a MOTS is signaled by ζ = 1 (Eq. (16)), and we always place the excision surface inside the MOTS 3 . The location of x exc is updated every time step before solving for the scalar and metric fields. The location of the excision point is chosen so that it is always interior to the MOTS, but lies exterior to (or directly on) the sonic line (for further discussion see Sec. III A and Sec. IV C). The location for r exc on our initial data slice depended on the strength of the curvature coupling C, which we detail in our Results, Sec. IV.

E. Fixed mesh refinement with a hyperbolic-ODE system
To achieve the necessary long term accuracy over thousands of m in evolution using limited computational resources, we evolved some simulations using a Berger and Oliger (BO) style mesh refinement algorithm [42]. Due to the nature of our initial data and perturbations, a fixed hierarchy suffices, with the higher resolution meshes confined to smaller volumes centered about the origin. For those runs, we typically used 4 additional levels beyond the base (coarsest) level, with a 2 : 1 refinement ratio between levels; specifically, we set the inner boundary for all levels at x exc , and the outer boundary locations at x l /m = {5.00, 2.30, 2.00, 1.75, 1.55}, from coarsest to finest (the initial horizon location r = 2m, and we chose the excision radius x exc so that r exc (x exc ) = 0.8 × 2m; see Eq. (21)).
The original BO algorithm was designed for purely hyperbolic systems of equations; to include the ODE constraint equations, we employ the extrapolation and delayed solution modification developed for such coupled elliptic/hyperbolic systems [43]. Here, for the hyperbolic equations (governing P and Q), the solution is obtained on the mesh hierarchy with the usual BO time-stepping procedure : one time step is first taken on a coarse parent level before two 4 steps are taken on the next finer child level, and this is repeated recursively down the mesh hierarchy. During this phase the ODEs are not solved, and where the values of the corresponding constrained variables (α and ζ) are needed to evaluate terms in the hyperbolic equations, approximations for these variables are obtained via extrapolation from earlier time levels. Instead, the ODEs are solved after the fine-to-coarse level injection phase of the hyperbolic variables, when the advanced time of a given parent level is in sync with all overlapping child levels (thus, on the very finest level this scheme reduces to the unigrid algorithm described in the previous section). For more details see [43].
One difference with our system of equations compared to that described in [43], is there some form of global relaxation method was assumed for the elliptics, while here the ODE nature of our constraint equations requires integration from the inner to outer boundary. This might complicate things for a general hierarchy with disconnect grids on a given level. Here, since we only have one grid per level, and each always includes the physical inner boundary, it is reasonably straight forward to integrate the ODEs during the solution phase of the algorithm (see also section 6.7 of [44]): we begin on the finest level, setting the boundary conditions as required at x = x exc , then integrate outward, using the solution at the last point on a child level as an initial condition for continued integration on the parent level (with the solution at interior points on the child level injected to common points on the parent level). This is schematically shown for two levels in Fig. 1. For ζ, the boundary condition at x = x exc is obtained by its corresponding evolution equation as described in the previous section; α at x = x exc is set by extrapolation from prior time levels, and we only globally rescale α to satisfy our outer boundary condition α(t, x = L) = 1 at times when all levels of the hierarchy are in sync.

IV. SCALARIZED BLACK HOLES: NUMERICAL RESULTS
A. Growth of "hair" from Schwarzschild initial data For most of our simulations we begin with exact Schwarzschild black hole initial data, which is a solution to the initial value problem in EdGB gravity. Specifically, at t = 0, for x > x exc (the initial excision radius as described in Sec. III D), we set We then evolve this, performing a survey of outcomes varying the EdGB coupling parameter λ (in the GR case λ = 0, as expected, the resultant numerical solution is static to within truncation error). Previous studies of static scalarized black hole solutions in EdGB gravity have found that regularity of the scalar field at the horizon places an upper limit on the coupling parameter. For the linear coupling case we consider, this is (see e.g. [5]) where r h is the areal radius of the horizon. For black holes much larger than this, r h ≈ 2m (the space time is close to Schwarzschild), though approaching the extremal limit a non-negligible amount of the spacetime mass m can be contained in the scalar field, and simply replacing r h with 2m to express the above in terms of our curvature-coupling parameter C (13) gives a poor estimate of the corresponding extremal value C extr . From Figure 4 of [5] we can infer a more accurate translation : We find the extremal limit is not characterized by the appearance of a naked curvature singularity, but instead the formation of a sonic line (and elliptic region) outside of the horizon of the black hole. Our measured extremal limit of C extr ∼ 0.23, as shown in Sec. IV C, is quite close to the above limit from [5]. Note that the "extremal limit" for a scalarized EdGB black hole is of a different nature than the extremal limits of Kerr or Reissner-Nordstrom black holes. The spin or electric charge of a black hole is set by the black hole's formation history: black holes of the same mass can have different spins or charges depending on the initial configuration and net charge and angular momentum of the matter that fell in to form the black hole. By contrast for an EdGB black hole the final scalar charge is set by the Gauss-Bonnet coupling, and Gauss-Bonnet curvature at the horizon, independent of its formation history [3,5].
Given that the Schwarzschild solution, of any mass, is valid initial data in EdGB gravity, we can certainly begin with superextremal black holes in our evolution (and again to be clear, here we use the term "superextremal" to refer to C > C extr ; there is no spin or charge in our numerical solutions). As we show below however, these develop elliptic regions outside the horizon. Moreover, our results in [1,2] show that trying to form a superextremal black hole from gravitational collapse of the dilaton field (in spherical symmetry) will result in an elliptic region appearing before a horizon. This suggests superextremal black holes in EdGB gravity exist in the regime of the theory governed by mixed-type equations, and their presence or "formation" (however that could be interpreted in a mixed-type problem) would mark a breakdown of the Cauchy problem. Also note that failure of the Cauchy problem is not a priori connected to regions of strong curvature or black hole formation; as further shown in [1,2], strong coupling and mixed-type character can be present for arbitrarily small spacetime curvature.
We first present results from evolution of Schwarzschild black hole initial data, and curvature couplings below the extremal limit. In all cases, if we move our excision radius sufficiently far interior to the horizon, we find that at some time an elliptic region forms in the interior. However, for these cases we can choose an excision radius closer to the horizon so that the evolution settles to a stationary state without any elliptic region forming in this new domain. As discussed in the introduction, we view this as a consistent initial boundary value evolution of EdGB gravity where the elliptic region is "censored" from the exterior hyperbolic region. In Figs. 2 and 3 we show examples of scalar hair growth for these cases (with the elliptic region excised), and their difference from the static "decoupled" scalar field profiles for a Schwarzschild black hole background (see Appendix B), for various curvature-couplings. These runs employed the fixed mesh refinement algorithm described above, with the base level grid having N x = 2 8 + 1 grid points, Courant-Friedrichs-Lewy (CFL) number of 0.25, and an excision radius at fixed at r exc = 0.95 × 2m. We find the scalar field settles down to solutions that differ little from the static decoupled scalar field profile, although the difference grows as the curvature coupling approaches the extremal limit (c.f. Fig. 3 of [5]). Nevertheless, in agreement with the results of [5], the difference of the full solution from the static decoupled limit solution remains small outside the black hole horizon. From convergence studies we find we can resolve the difference of the scalar field profile from its decoupled value well within truncation error; see is at x/m ≈ 1.48, and spatial infinity is at x/m = 5. As expected, the decoupling limit approximation improves the further C is away from the extremal limit C extr ∼ 0.23. For other run parameters see Sec. IV A.

B. Perturbed Schwarzschild initial data
Schwarzschild initial data is not generic; in particular, the scalar field is (initially) only growing in response to the Gauss-Bonnet curvature source. To investigate a slightly broader class of initial conditions, here we add a small, mostly ingoing propagating perturbation to φ outside the horizon: This family of initial data (rescaled "bump functions" multiplied by a Gaussian) is smooth and compactly supported outside the initial black hole horizon for a > 2m. With a fixed curvature-coupling C, we find that we can stably evolve an initial black hole plus scalar field bump if the amplitude of the latter is sufficiently small; or equivalently if the metric curvature measured by the Ricci scalar R induced by the scalar field bump is sufficiently small. For our initial data for α and ζ, we set their values at the excision surface as in Eq. (25), and then integrate outwards in r. An example of such a case is shown in Fig. 7. When the induced curvature is large, an elliptic region forms outside the black hole horizon (and soon after that the code crashes). As a rough estimate, we find this occurs when |R × λ| ∞ 0.1. This result is consistent with our earlier findings of collapse of a scalar field pulse without any interior black hole [1,2].
C. Internal structure of an EdGB black hole, and the near extremal limit For small values of the curvature coupling we can consistently excise any elliptic region that forms interior to the horizon. We always excise the Schwarzschild curvature singularity at r = 0, and in none of the cases we have run (small or large curvature coupling) do we see any signs of a curvature or scalar field singularity forming away from r = 0 while the equations remain hyperbolic. We can typically evolve for some time after an elliptic region forms without excising it before the code crashes, however this as an artifact of finite resolution, and we can only expect FIG. 4: Convergence and truncation error estimate at t ∼ 100m for the scalar field φ. Shown are point-wise differences of the solution computed with different resolutions; the decrease going to successively higher resolutions is consistent with second order convergence, and the magnitude for a given pair is an estimate of the error in the scalar field profile at those resolutions. Comparing with Fig. 3, we see we can resolve the difference of the scalar field from the decoupled value well above truncation error for the range of curvature-couplings considered here. For other run parameters see Sec. IV A. The lack of convergence of the scalar field near our compactification point x/m = 5 is a transient effect from small reflections of the scalar field off the grid refinement boundaries as the scalar field grows around the black hole. simulation. The small "features" in the Ricci scalar at the t = 10m time slice are located at grid refinement boundaries, and converge away with higher base resolution (compare with Fig. 6). For other run parameters see Sec. IV A.  (28)) outside the horizon. The mass of the initial Schwarzschild black hole is m s = 20, while the scalar bump adds an additional ∼ 3 in mass, giving a net mass m = 23; thus with λ = 1 the curvature coupling is C = λ/m 2 ≈ 2 × 10 −3 . Some of the scalar field falls into the black hole, while some disperses to infinity, and at late times the scalar field approaches the decoupled scalar field profile. (This is from unigrid evolution, with N x = 2 12 + 1, and CFL number 0.5.) convergence to the continuum limit using our hyperbolic scheme until the first appearance of the sonic line. In other words, to say anything rigorous about what might occur within the elliptic region would require solving a mixed-type problem, and we do not have the tools for doing so.
With increasing values of the curvature coupling approaching the extremal limit (27), as the scalar field grows, the location at which the sonic line first appears moves closer to the MOTS. Prior to this, we excise some distance within the MOTS, though when the sonic line appears we increase the excision radius to be at the sonic line 5 . We then employ a "high water mark" strategy during subsequent evolution, increasing the excision radius to match the location of the sonic line if it grows, though do not reduce the excision radius if the sonic line shrinks (presumed to be happening if the characteristic discriminant on the excision boundary increases in magnitude away from zero).
For cases where the elliptic region remains censored, we typically find that initially the sonic line does grow, and then (presumably) shrinks within the excision radius as the solution settles to a stationary state. For interest, we estimate the location of the sonic line by extrapolation, as follows. Recall, the equation for the radial characteristics is (Eq. (17)) where A, B, C are functions of α, ζ, P, Q, and their radial derivatives. The characteristics thus satisfy  is sufficiently late to give a good estimate of their static end state values (see Fig. 9). Due to our high water mark excision strategy, the excision point represents the largest measured radius the sonic line attained during evolution. Extrapolating the curves for the MOTS and sonic lines, we estimate the "extremal" coupling for our initial data to be C extr ∼ 0.23. These results were obtained with unigrid evolution using N x = 2 13 + 1 (corresponding to the highest resolution curves shown in Fig. 9), and CFL number 0.2 (0.1) for C ≥ 0.17 (C < 0.17).
using a set of points adjacent to the excision boundary. In typical cases for Schwarzschild initial data, and subextremal curvature couplings, this estimate suggests the true location of the final stationary sonic line lies within ∼ 94% of its maximum value (the excision point); see Fig. 8 for a survey of the late time values of the excision radius, MOTS location and sonic line estimate, and Fig. 9, for the evolution of these quantities for one example (including several resolutions). In the latter figure, the shrinking of the MOTS after some initial growth coincides with violation of the null convergence condition (R µν k µ k ν ≥ 0 for all null vectors k µ ; see e.g. [45]), which is known to occur in EdGB gravity (for more details in the spherical collapse problem see [2]). A plot of R µν k µ k ν is shown in Fig. 10. We note that the stable violation of the NCC is thought to be a key ingredient in the construction of singularity free cosmological and black hole solutions (for a review, see e.g. [46]). This violation of the NCC appears to be transient: as the scalarized black hole settles to a stationary solution, we find the horizon stops shrinking and the region of NCC violation disappears. The slow increase in the horizon size for t/m 50 is due to numerical error; we find it converges to zero with increasing resolution. For curvature couplings above C extr , the sonic line can move outside the MOTS, or initially appear outside it. Linearly extrapolating the data shown in Fig.8 to the location where the late time MOTS will cross the sonic line, we estimate C extr ∼ 0.23, close to but slightly larger than the value C extr ∼ 0.22 coming from seeking exactly static EdGB black hole solutions with a non-singular φ field on the horizon [5] (though even beyond caveats with our extrapolations, we do not expect these two methods to give identical numerical values for an extremal coupling). In Fig. 11 we show an example evolution of Schwarzschild initial data with superextremal curvature coupling. We see the sonic line quickly overtakes the black hole horizon, leading to a "naked" elliptic region.

V. DISCUSSION
In this article we have presented numerical evidence that in spherical symmetry, and for sufficiently small curvature couplings (what we call subextremal), EdGB black holes are nonlinearly stable. For subextremal couplings even moderately close to the extremal limit, solving the decoupled scalar equation for the scalar field profile provides a good fit to the numerical solution obtained in the full theory. Beginning from Schwarzschild initial data, instead of developing a non-central curvature singularity in the interior as was found for static solutions [5], we find the formation of a sonic line and elliptic region in the interior. Our treatment of the EdGB equations as hyperbolic does not allow  Fig. 8), and from runs at three resolutions : the labels 0, 1, 2 refer to unigrid resolutions N x = 2 11 + 1, N x = 2 12 + 1, and N x = 2 13 + 1, respectively (a CFL factor of 0.2 was used in all cases). At early times as the sonic line grows, we increase the location of the excision surface to match; after reaching a maximum radius, the sonic line presumably starts to shrink again, and then the curves in the figure show an estimate of this location based on extrapolation of the characterstic speeds (see Sec. IV C). The resolution study demonstrates that at late times we are converging to a static solution (in the vicinity of the horizon).  11: Example evolution of Schwarzschild initial data with superextremal curvature coupling: C = 0.275. We excise along the sonic line, so the excised region can be thought of as the elliptic region. The sonic line very quickly reaches and overtakes the MOTS, and the elliptic region subsequently becomes "naked", shortly after which the code crashes. We note our code crashes if we do not excise the elliptic region, regardless if it is interior or exterior to the black hole horizon. (This is from a unigrid evolution, with N x = 2 13 + 1, and CFL number 0.2).
us to conclude anything about possible extensions of the spacetime into the elliptic region. For subextremal black holes, our statement about their stability relates to the region exterior to the horizon, and assumes that our excision strategy used to eliminate the interior elliptic region is self-consistent (which is supported by the stability/convergence of the corresponding numerical evolutions). For superextremal cases, the sonic line forms or evolves to be outside the horizon, meaning we cannot excise it, and we would need to treat the exterior equations as mixed-type to obtain sensible solutions (or said another way, then the exterior evolution ceases to satisfy a well-posed Cauchy initial value problem). The particular value of the curvature coupling we find for the extremal limit is similar to, but slightly different from that given for static EdGB black holes solutions [3,5]; this is not particularly surprising given we are dynamically forming them from Schwarzschild initial data.
There are various ways in which this work could be extended. One is to explore a wider class of initial conditions; for example, collapse to a black hole from a regular matter source, whether the pure EdGB scalar field as in [2], or coupled to another source of matter driving most of the collapse (e.g. extending the study of [27], which only considered the decoupled EdGB field on top of Oppenheimer-Snyder style collapse, to the full EdGB equations). Recent work suggests that whether or not scalarized black holes form in the theory depends on the functional form of f (φ) (see Eq. (1)) [6,47]; with the methods presented in this paper one could explore these questions with numerical solutions to the full theory in spherical symmetry. Another future direction is to study critical collapse in EdGB gravity using adaptive mesh refinement.
Finally, this work could be extended by considering numerical solutions of EdGB gravity in axisymmetry, or without any symmetry restrictions. This would couple in propagating metric degrees of freedom, and hence introduce a qualitatively different aspect of the theory not available in spherical symmetry. If, similar to the conclusions found here and in [1,2], there exist subsets of initial data that offer well-posed hyperbolic evolution, then EdGB gravity may still be viable as an interesting modified gravity theory to confront with gravitational wave binary merger data. On the other hand, if the linear analysis in [14,15] that EdGB gravity is generically ill-posed in a particular gauge applies to all gauges, then the well-posed cases we have found could be an artifact of spherical symmetry, and including any gravitational wave degrees of freedom would render the theory ill-posed. versus time computed on the base level of a fixed mesh refinement run (other levels also show second order convergence to zero). The value N x is the resolution of the base grid.