{'fraction_non_alphanumeric': 0.03402688155119771, 'fraction_numerical': 0.010576348012213164, 'mean_word_length': 4.622013802866749, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fixing the dynamical evolution in scalar-Gauss-Bonnet gravity

One of the major obstacles to testing alternative theories of gravity with gravitational-wave data from merging binaries of compact objects is the formulation of their field equations, which is often mathematically ill-suited for time evolutions. A possible way to address these delicate shortcomings is the fixing-the-equations approach, which was developed to control the behaviour of the high-frequency modes of the solutions and the potentially significant flow towards ultra-violet modes. This is particularly worrisome in gravitational collapse, where even black hole formation might be insufficient to shield regions of the spacetime where these pathologies might arise. Here, we focus (as a representative example) on scalar-Gauss-Bonnet gravity, a theory which can lead to ill-posed dynamical evolutions, but with intriguing stationary black hole physics. We study the spherical collapse of a scalar pulse to a black hole in the fixing-the-equations approach, comparing the early stages of the evolution with the unfixed theory, and the later stages with its stationary limit. With this approach, we are able to evolve past problematic regions in the original theory, resolve black hole collapse and connect with the static black hole solutions. Our method can thus be regarded as providing a weak completion of the original theory, and the observed behaviour lends support for considering previously found black hole solutions as a natural outcome of collapse scenarios.


I. INTRODUCTION
The observation of gravitational waves emitted by binaries of compact objects opened a new possible channel to confront general relativity (GR) with alternative theories of gravity. So far, in the dynamical and strong-field regime of binary mergers, consistency with the expectations of GR has been confirmed mostly through null and consistency tests [1][2][3][4]. These tests are however insufficient to rule out specific alternative theories, as predictions for the merger waveforms of compact objects beyond GR have only been obtained in a handful of cases, if no approximations are made (see e.g. [5][6][7][8][9][10]). Such waveforms require numerically solving the field equations of a specific theory, i.e. a full system of partial differential equations (PDEs). However, a major obstacle in this enterprise is the generic lack of well-posedness of the initial data (Cauchy) problem in most alternative theories of gravity [11][12][13][14][15][16]. This problem, arising from the mathematical structure of the underlying field equations, stands in the way of obtaining predictions of beyond-GR effects, unless a suitable approach is introduced to address this issue.
As a representative example of this problem, we focus here on Horndeski gravity [17], a generalized class of scalar-tensor theories, in which a new dynamical degree of freedom, a scalar field, modifies GR through a nonminimal coupling. The Horndeski action is constructed so as to encompass all the possible terms leading to second-order field equations, thus avoiding Ostrogradski ghosts [18]. If the scalar field is invariant under constant displacements (shift symmetry), a no-hair theorem states that stationary spherically symmetric and slowly rotating black holes (BHs) are equivalent to their GR counterpart [19,20]. The only exception to this theo-rem arises when the scalar field couples linearly 1 to the Gauss-Bonnet (GB) invariant [20,21], defined in terms of the Riemann tensor R abcd , the Ricci tensor R ab and the Ricci scalar R as When shift symmetry is broken, the phenomenology of BH solutions is even richer. A notable example is the occurrence of spontaneous scalarization in a large sub-class of Horndeski theories [22][23][24][25][26][27][28][29][30][31]. This phenomenon, like the analogous spontaneous scalarization of neutron stars in Brans-Dicke-like scalar-tensor theories [32], amounts to a growth of the scalar field around BHs, prompted by a linear tachyonic instability of the scalar and eventually quenched by the non-linearities of the problem. A tachyonic instability arises naturally in many sub-classes of Horndeski gravity [33,34], including scalar-GB (sGB) gravity, which represents the case study of this paper. Scalarized compact objects acquire an additional (scalar) charge, which not only affects the gravitational interaction but which can also yield monopole/dipole emission, which would impact the gravitational wave signal [35][36][37].
The appearance of intriguing physics in sGB gravity is not limited to stationary BHs, as dynamical scalarization and de-scalarization might occur in binaries [38][39][40][41] or in gravitational collapse [42][43][44][45][46]. However, as already stated, the well-posedness of the Cauchy problem is a necessary condition to study these timedependent systems without approximations. The criteria of Hadamard [47] states that the Cauchy problem for a system of PDEs is well-posed if there exists a unique solution that depends continuously on its initial data. Since proving this statement is a difficult task [48], one usually attempts to show that the PDE system is strongly hyperbolic, which is a sufficient condition for well-posedness under suitable initial and/or boundary conditions [49,50]. This requires writing the system of equations as a quasilinear first-order system and showing that the principal part has real eigenvalues and a complete set of eigenvectors.
For sufficiently weak couplings and smooth initial data, local well-posedness of sGB gravity can be established [51,52] and successful studies have been presented in this regime [12][13][14]. However, for stronger couplings and/or scenarios where higher frequency modes develop, the hyperbolic structure of the system can break down. Specifically, it can be shown that the underlying equations can switch their character from hyperbolic to elliptic [12,13,39]. More specifically, during the evolution one of the characteristic speeds (i.e. one of the eigenvalues of the characteristic matrix) becomes imaginary. This breakdown of the character of the system resembles the Tricomi equation [12][13][14].
Faced with this possibility, a potential way through is to modify the system of equations so that short wavelength features (with respect to some chosen scale) in the solution are somewhat controlled, while longer wavelength modes are faithfully tracked [53]. This idea, applied successfully in many studies [54][55][56][57][58], is rooted in the expectations that the theory at hand is the truncation of a putative theory which has a sensible behavior, and that the problems encountered are an artifact of truncating it at a given order. As an analogy, consider the pinch-off of a soap bubble. Some time prior to and after the pinch-off, its surface can be faithfully described with hydrodynamics, but around the pinch-off itself one needs to use a different theory (molecular dynamics) or a suitable effective regularization of the hydrodynamics equations to bridge the earlier and later stages. If the solution obtained is largely insensitive to the details of the modification introduced and, stronger yet, such solution recovers the characteristics expected in the early and late stages of the dynamics, strong support for the approach yielding a faithful representation of the UV complete underlying theory can be argued.
In the current work, we focus on a problem that illustrates this approach and behavior. Namely, we study the collapse within sGB and show that the latter gives rise to a scalarized BH after passing through a regime where the original system of field equations develops a pathological behavior. We implement a technique that controls analytical obstacles and related numerically induced problems in the time evolution, and show that this approach allows for bridging regimes where the theory is sensible. This latter property -not necessarily expected a priori -is undoubtedly a welcome sight that indicates the robustness of scalarized BHs as fixed point solutions of the theory.
This work is organized as follows. In section II we present the theory of sGB, its equations of motion, and how we replace those equations with the fixed system. Then, in section III, we specialize to spherical symmetry, expressing the equations of motion in the full and in the fixed theory, and presenting the characteristic speeds. Moreover, we discuss the numerical set up of the simulations. The results of the latter are shown in section IV. Finally, we draw our conclusions and possible implications of our results in section V. By convention, we use units such that c = G N = 1.

II. THE THEORY AND THE FIXING
The action for sGB gravity in vacuum (also referred to as full theory henceforth) is given by [59] where the GB invariant G given in Eq. (1) couples to the scalar field through the function f (φ). A variation of the action gives the equations of motion where the sources of the Einstein and Klein-Gordon equations are with As discussed in the Introduction, the system (3)-(4) is not strongly hyperbolic for a large class of choices of the coupling function f (φ). In order to circumvent this problem we introduce here a fixing-the-equations technique [53], which was inspired by dissipative relativistic hydrodynamics [60][61][62]. The main idea is to modify in an ad hoc way the higher-order contributions to the equations of motion, by replacing them with some auxiliary fields, and let the latter relax towards their correct value through a driver equation. We replace the system (3)-(4) with the following fixed system of equations.
Here, we gathered the auxiliary fields in the vector u = (Γ ab , Σ), while {ξ, τ } are constant timescales (although in principle they can differ for each component of u) controlling how the auxiliary variables approach to the full theory solution S = T (G ) ab , S (G ) . Notice also that equation (10) is one of many options that can force the auxiliary variables to approach their corresponding correct values. Crucial in this respect is the presence of the first time derivatives ∂ t , which introduce a "preferred" time direction, i.e., dissipation. The choice, made above, of a wave-like driver equation for the auxiliary fields is in a sense natural given the underlying hyperbolic structure of the Einstein equations. We will discuss some issues related to this choice later. In the following, we focus our analysis and numerical evolution on spherical symmetry.

III. EVOLUTION IN SPHERICAL SYMMETRY
In order to study the non-linear dynamics we write down the equations (3)-(4) or (8)-(10) as a first-order (in time) evolution system, by applying the techniques used in Refs. [12,15].

A. Evolution equations
We adopt the following ansatz for the metric in polar coordinates and we introduce two new functions to express derivatives of the scalar field-which is a generic function of time and radial coordinate φ = φ(t, r)-by defining In spherical symmetry, the tensor T (G ) ab has four nontrivial components. This means that in the fixed theory we make use of five auxiliary variables, Σ(t, r), Γ 11 (t, r), Γ 12 (t, r), Γ 22 (t, r) and Γ 33 (t, r). Henceforth, we place a bar on top of any function evaluated in the full theory, in order to distinguish it from the fixed theory.
Schematically, the full system (3)-(4), written in spherical symmetry with the ansatz considered and after some manipulation of the equations, reduces to three evolution equations and two constraint equations whereĒ GB P ,C GB A ,C GB B have a lengthy and uninformative expression that can be found in the appendix of reference [13]. The first two equations come from the consistency of the partial derivatives of the scalar field, the Full theory Fixed theory Table I. List of the variables used and characteristic speeds computed in the full and in the fixed theories.
third comes from the Klein-Gordon equation (4) and the last two come from the constraints of the theory (3).
On the other hand, the evolution equations for the fixed theory in spherical symmetry are explicitly given by where H ≡ e −A+B ∂ t u and J ≡ ∂ r u. The explicit form for S is given in Appendix B.
The two systems we are considering are described by two constraint equations and M evolution equations, where M = 3 for the full theory and M = 3 + 3 · 5 = 18 for the fixed theory. If we label with V the variables corresponding to the evolution equations (with dimension M ) and with W = (A, B) the two metric variables, we can write the system in the following compact form An overview of the variables and of the equations can be found in Table I. A useful tool to check the hyperbolicity of the system during the evolution is to evaluate the characteristic speeds. The latter can be calculated from the characteristic matrix of the full system where we introduced the covector χ a and defined with the indices I and J running from 1 to either 2 or M according to the function and the equation considered. The characteristic speeds are then defined as where, if S is invertible, The expression for the characteristic speeds for the full theoryv GB ± is very long and we do not show it here. On the other hand, the speeds in the fixed theory are equivalent to those of GR For the rest of the paper, we choose the coupling function where η and ζ are parameters of the theory. The quadratic term proportional to η is the one responsible for spontaneous scalarization provided that η > 0. The addition of a quartic term proportional to ζ ensures the existence of radially stable spherically symmetric BH solutions, provided that ζ −0.8η [25]. In all the cases considered we choose ζ = −6η, which largely ensures this requirement. The hyperbolicity of these class of theories was largely studied in previous works. The ansatz we chose for spherically symmetric evolution in polar coordinates leads to a breaking of the hyperbolicity in the linear theory [12,13], and we confirmed that it also happens for this choice of the coupling (see e.g., Figure 1 in the Results section). Attempts to evolve these theories were made with the help of excision [14] and/or generalised harmonic gauge [39,51,52]. Both help to ensure a (local) well-posed formulation of the problem by adopting a convenient gauge as well as excising a significant region inside the BH, where short wavelength modes are excited. However, even with these techniques, there exist regions inside whose the system dynamically develops elliptic character, thus breaking well-posedness. The approach taken with the fixing-the-equations method aims to provide a robust way to explore a potentially significantly larger region of the parameter space.

B. Diagnostic tools
For diagnostic of the solutions obtained we find it useful to define the following norm, valid for any quantity F with two different definitions in the full and fixed theory where ε = 10 −5 is to ensure that the denominator never goes to zero. We stress that since the dynamical evolution cannot be followed in the full theory past the formation of the elliptic region, in our analysis we often compute the quantities F GB using the variables of the fixed theory.
From the metric at infinity we can define the Misner-Sharp mass [63] as During the evolution, we rescale all the dimensionful quantities such as M MS (0) = M → 1. For example, we redefine r → r/M and η → ηM 2 . Finally, we compute the null convergence condition (NCC) as R ab a b ≥ 0, where the contraction of the Ricci tensor is taken with null vectors defined as a = {e −A , e −B , 0, 0}. A satisfied NCC is expected to hold in GR for BHs and cosmological settings, and it is a powerful tool in many theorems [64,65]. We expect anyway that the NCC is violated in sGB gravity [13]. It is worth noting that since R ab contains time derivatives, one can get rid of them using the equations of motion. This means that the NCC has a different expression when evaluated in the full or in the fixed system.

C. Initial data
As initial data (ID) for φ, P and Q, we take two different families of function. The ID labelled as type I is a static Gaussian pulse where a 0 and r 0 are the amplitude and grid location of the pulse respectively, and w 0 is the root-mean-square width of the Gaussian pulse. The ID labelled as type II is an (approximately) ingoing pulse, defined as For the metric variables A and B we do not impose any initial data, as one can just solve the constraints equations given the initial configuration for the scalar, as explained below. Finally, we set all the auxiliary variables to 0 at t = 0.

D. Numerical scheme
We solve the two systems of equations with the same methodology, by using a fully constrained evolution scheme which works as follows.
We discretize all the equations over the domain r ∈ [0, r max ]. By choosing a resolution N we determine the radial step ∆r = r max /N . To define the time step ∆t we adopt a Courant parameter λ r = 0.25 such that ∆t = λ r ∆r satisfies the Courant-Friedrichs-Levy condition. To ensure the regularity of the equations at the origin, we must impose These conditions apply both to the full and the fixed system. Moreover, for the latter, we also set the following conditions for the auxiliary variables At the outer boundary of integration, we set approximately outgoing boundary conditions for the variables V, while the metric variables are automatically determined by the solution of the equations. We perform the evolution as follows. At each time step, we solve the two first order equations C A and C B from 0 to r max in space by using a fourth-order Runge-Kutta (RK) scheme to obtain a solution for A and B. As to carry the integration one needs to know the value of the functions V at intermediate virtual points, we evaluated this with fifth-order Lagrangian interpolator. At r = 0, it is sufficient to set A(t, 0) = 0 and B(t, 0) = 0 to automatically satisfy conditions (45). Once a solution is found, we exploit the remaining gauge freedom to set A(t, r) → A(t, r)−A(t, r max ). This ensures that the proper time of an observer located at r = r max is identically t.
Once the metric functions are obtained, we integrate the other variables V in time through the method of lines by using a fourth-order accurate strong stability preserving RK integrator. For convenience, we also added sixthorder Kreiss-Olliger dissipation. In every equation, we discretized the radial derivatives using fourth-order finite differences operators satisfying summation by parts [66].

IV. RESULTS
In this section, we show how the fixing-the-equation method allows for a stable evolution in the regime where the full theory would break down due to the appearance  Table II. Summary of the runs made, where we specify an identifying label, in which theory we are solving the equations, the kind of initial data, the parameters of the theory, the parameters of the fixing and in which figure they are used. The label sGB indicates the full theory, while the label GR indicates a theory with minimally coupled scalar field. of an elliptic region. A summary of the numbered evolutions (EV#) studied in this Section is shown in Table II.
First of all, we show that gravitational collapse in the full and fixed theories produces approximately the same dynamics. We evolve the same pulse in the two theories with the choice of parameters summarized in Table II, with labels EV1 and EV2. In the comparison, we pay attention to the characteristic speeds of the theories. In Figure 1 we plot, as a function of time, the maximum and minimum ofv GB ± and v GB ± . We stress that the functional form of the two quantities is the same, and a bar (no bar) denotes quantities evaluated in the full (fixed) theory. During the first stage of the gravitational collapse, around t ≤ 15M , the characteristic speeds are almost indistinguishable and constant during the evolution, showing that the character of the full theory is strongly hyperbolic and well reproduced by the fixed one. Nevertheless, at t ∼ 22M, one of the characteristic speeds crosses zero becoming imaginary shortly after, and therefore the system changes character from hyperbolic to elliptic. It is worth noting that even if the equations in the fixed system are different, their variables reproduce quite well the breaking of hyperbolicity of the full theory. We recall that the true characteristic speeds of the fixed problem are v F ± , and that they never change sign. This means that for every case, for t 22M the full theory cannot be evolved, but we can always follow the evolution of the fixed system.
With this in mind, we evolve a pulse in the fixed theory, for a choice of parameters that would break hyperbolicity in the full theory. In Figure 2, we give an estimate of how the full theory would behave, by computing the quantities of this theory with the variables of the fixed one. Here, we plot time snapshots of the evolution EV3, summarized in Table (II). In the top panels, we show the characteristic speeds, v GB ± and v F ± , and the norm ||C A ||. In the bottom panels we show the scalar profile φ and g tt . The first panel shows results at t = 2.3M . One can see that the full system would develop an elliptic character, confined only in a limited region of space where the characteristic speeds would become imaginary and D < 0. This region is denoted by an orange vertical strip. Moreover, we plot a blue region defined by the conditions v GB ± < 0 or v GB ± > 1. This region roughly corresponds to that where the deviation between C F A and C GB A becomes noticeable. This pattern remains in the second and third panels, the former being a snapshot at BH formation at t 14.6M , while the latter is a snapshot of the later evolution at t = 35.4M , where we expect scalarization of the BH.
It is worth noticing that the region where ||v ± || is of order 1 also corresponds to that where the equations of the two theories differ significantly. Outside this region, the equations and speeds of the two theories are very close. Thus, we conjecture that the physical evolution of the scalar field outside the region where ||v ± || ∼ 1 is trustworthy. We checked that this behavior persists at late times, when the scalar field grows considerably, as a result of spontaneous scalarization.
Another way to understand how the fixing works consists of evaluating the NCC. In Fig. 3 we compare the condition for the fixed system, for the full one evaluated with quantities in the fixed system, and for an evolution with equivalent initial data, but with a minimally coupled scalar field (labelled as EV4 in Table II). We can see that at the formation of the elliptic region (upper panel), the NCC of the fixed system resembles that of GR at small radii, while it almost overlaps with the condition evaluated for the full system at large radii. This latter trend persists throughout the whole evolution, as one can notice in the lower panel of the figure, evaluated at t = 35.4M . The absolute value of the NCC is here expressed in logarithmic scale, we denote the regions where it becomes negative by a dashed line style. These results agree with expectations based on the analysis carried in [13].
The previous analysis heuristically shows that the numerical evolution of the fixed theory is reliable, especially outside the apparent horizon (when the latter forms). Thus, we evolved a few cases with different combinations of η and ζ to see whether spontaneous scalarization of the final BH that forms from the collapse is achieved. We consider one representative case, summarized in Table II as EV5. The results for the gravitational collapse of this case are shown in Fig. 4. The blue lines indicate time snapshots of the radial profile of the scalar field φ after the formation of an apparent horizon. After the collapse, the scalar field surrounding the BH settles down to a static solution. We compare this equilibrium configuration with the one obtained by solving the static spherical problem in the full theory (c.f. [25]). This static solution can be obtained simply by requiring all the functions to be time-independent, or equivalently by replacing ∂ t φ = ∂ t P = ∂ t Q = 0, in Eqs. (3)-(4). By doing so, the equations are reduced to a set of ordinary differential equations, which can be easily integrated by imposing proper boundary conditions at the horizon of the BH and asymptotic flatness at the outer boundary. We clearly see that the evolution drives the scalar field to match the static profile. We stress that we obtained similar results with different values of the coupling and different choices of initial data.
A few caveats to discuss are the following. The dynamical evolution of the scalar field does not seem to asymptote well to the static solution in a region very close to the apparent horizon. We conjecture that this behaviour is due to a small portion of the elliptic region of the full theory leaking out of the horizon. Therefore, even if we can "heal" the character of the field equations in the fixed theory, we still expect some discrepancy in this region of large gradients and high-frequency energy flows. On the other hand, the scalar field does not fully match the static case at large radii either, but this is simply because we could not run the evolution any longer. At the level of implementation, possible improvements include using horizon penetrating coordinates rather than polar coordinates and/or excising the region inside the apparent horizon.
In order to understand how the fixing behaves, we studied the effect of different values of ξ and τ . We are interested in the limit τ → 0 and ξ → 0. However, note that for values of ξ lower than the coupling constant η, the evolution breaks down at early stages and cannot be followed, as expected from [53]. This is an issue that can be addressed by coupling both parameters in a suitable manner. Otherwise, since the decay rate is ≈ τ /ξ, adopting too small a ξ value would yield a stiff system, which needs to be handled with care at the numerical level. In our simulations we could set τ identically to zero as the results do not change with respect to cases where τ /η 1. Notice that a static, radially dependent, solution would induce a discrepancy between u and S given by ≈ ξ∇ 2 u. One can recover the static behavior as ξ → 0 for the par- ticular fixing employed. 2 Here, we employ a single value for ξ. For completeness, we also checked cases where we adopt five different constants-one for each auxiliary field-and found that we can lower a subset of these constants, while still getting a proper evolution. Namely, by lowering by a factor 10 the constants controlling the equations for Σ, Γ 22 and Γ 33 , we observed that the results obtained do not qualitatively change.
On the other hand, in Fig. 5 we show the effects of increasing ξ. The three panels are analogous to Fig. 4, but with η = 6 and ζ = −6η, and each of them assumes an increasing value of ξ. The runs correspond to EV3, EV6 and EV7 in Tab. II, respectively. It is clear that the scalar evolution matches well the static prediction only for the lowest value of ξ. This is not surprising, as large values of ξ imply that the fixed system's faithfulness to the original one degrades.
In the current work we studied the dynamical collapse and related hair BH formation in sGB gravity. Such a theory has been the subject of scrutiny in recent years, which has revealed interesting behaviors for BH physics [22][23][24][25][26][27][28][29]. However, as it has been discussed elsewhere [12-14, 67, 68], it also presents important mathematical challenges in regimes of high couplings and/or short wavelengths. Such challenges stem from the equations of motion breaking their hyperbolic character. This makes it impossible to explore important questions such as the formation and non-linear stability of compact objects, unless further steps are taken to address the the theory's shortcomings. Further, even in the weak coupling regime, physical effects abound that induce short wavelength behavior from long wavelength, e.g., focusing and collapse. For instance, a gravitational (or scalar) wave passing by a BH could be focused and excite significantly short wavelengths, potentially breaking (local) hyperbolicity. Unless a BH generically shields such regions away from physical observers, one could argue that extensions of GR that present these shortcomings would  have limited physical relevance. It behooves theorists to try and devise ways to explore such question.
Here, we adopted a strategy dubbed as fixing-theequations, which introduces further auxiliary fields, constrained within some given scale, to force a suitably modified system to agree with the original equations of mo- tion, while at the same time controlling the behavior at shorter wavelengths [53]. This approach, to a certain extent, recovers some of the higher energy degrees of freedom, which in an effective field theory sense are integrated out and whose influence results in corrections to the Einstein-Hilbert action. Controlling the high frequency modes of the solution, this approach eliminates (or significantly reduces) the mathematical pathologies of the original system, thus allowing one to push evolutions further. Note that the fixing approach breaks Lorentz symmetry, in particular through the τ -term in Eq. (10) which sets a scale-in time and wavelength-below which such breaking is especially evident. The presence of this term is in agreement with the fact that the UV completion of the theory has to be non-standard [69].
In our particular application, the fixing allowed us to resolve the collapse to a scalarized BH. Remarkably, this dynamically formed BH agrees with static solutions of the original unfixed theory. We stress that this was not necessarily expected, as our method does introduce modifications in the theory. As such, the agreement of the solutions found is highly non-trivial and indicates that the fixed system successfully bridges between stages where/when the original system is well behaved. We also stress that we find regimes where hyperbolicity breaks outside the apparent horizon for sGB as already noted in previous works [12,14,39,44]. Two comments are in order here: (i) that sGB shows fails to be hyperbolic outside the horizon implies one can not assume that a horizon will hide such failure from outside observers. (ii) this failure should not happen in a proper UV completion of the theory, and the fixed version allows to by pass such issue. Last, we note that the fixed system yields a solution that approaches the one expected in the static regime. This is analogous to the behavior observed in an effective-fieldtheory-derived gravitational theory where BH accretion was studied and the late time solution agreed with the static one [55]. If this observed robustness could be extrapolated to more general solutions and other extensions to GR, it would lend strong support to the physical relevance of such solutions and conclusions drawn from perturbative studies. Furthermore, this agreement supports the idea that the approach employed can be regarded as furnishing a weak completion of the original theory.  gence factor c p and the order of convergence as follows We use the Misner-Sharp mass (38) to compute c p and p, as shown in Fig. 6. In the top panel, we show the convergence for EV8, finding a fourth order convergence as expected. In the lower panel the convergence for EV3 is displayed. In this case, until the collapse we find fourthorder convergence, but after the formation of the apparent horizon the convergence is somewhat degraded. This is expected due to the large gradients in the collapse front.