Gravitational Collapse in Cubic Horndeski Theories

We study spherically symmetric gravitational collapse in cubic Horndeski theories of gravity. By varying the coupling constants and the initial amplitude of the scalar field, we determine the region in the space of couplings and amplitudes for which it is possible to construct global solutions to the Horndeski theories. Furthermore, we identify the regime of validity of effective field theory as the sub-region for which a certain weak field condition remains small at all times. We evolve the initial data using the CCZ4 formulation of the Einstein equations and horizon penetrating coordinates without assuming spherical symmetry.


Introduction and Summary
The detections of gravitational waves produced in mergers of compact objects [1,2] have revolutionised the field of gravitational physics, giving rise to the era of gravitational wave astronomy. Thanks to recent upgrades of the detectors, gravitational waves detections are made almost on a weekly basis. Therefore, we now have an unprecedented amount of data that gives us access to the strong field regime of gravity. The situation is only going to get better in the future, with new detectors gradually added to the network in the coming years and a forthcoming third generation of detectors such as the Einstein Telescope and ultimately Lisa, a space-based observatory. Therefore, very soon we will enter the era of precision gravitational wave astronomy. These advancements offer the opportunity (and carry the duty) to test Einstein's theory of general relativity (GR) using gravitational waves. One of the main challenges in doing these tests is to come up with templates of waveforms in alternative theories of gravity. One possibility is to focus on those phases of the binary that can be treated using perturbation theory, namely the inspiral [3] and the ringdown phases [4,5] respectively. However, the present data suggests that the corrections to GR are small. Therefore, one may hope that there is a better chance to detect some deviations from GR in the strong field regime, namely in the merger phase, where some effects may be enhanced. This would be the case for deviations from GR that are sourced by spacetime curvature, such as higher derivative corrections. So far the merger phase has been modelled phenomenologically [6,7], or by treating the deviations from GR perturbatively [8][9][10][11][12]; only the so called scalar-tensor and scalar-vector-tensor theories of gravity have been considered in their full non-linear glory in all phases of the binary [13][14][15][16].
Another difficulty is that there are many alternative theories of gravity, and each one of them modifies GR in a different way: adding new fields, breaking some symmetries, adding new terms to the action, etc.. At the moment there is no theoretical consensus nor any experimental evidence that favours a particular theory. Each modification of GR should be reflected in a unique way in the corresponding waveforms and hence the interest in analysing gravitational waves in alternative theories of gravity. However, in many of these theories it is not known whether the initial value problem is well-posed. Without a well-posed initial value problem, one cannot possibly simulate the non-linear regime of the theory on a computer and obtain the desired waveforms. There have been some recent efforts that have successfully managed to construct well-posed formulations of certain modified theories of gravity of physical interest [17][18][19]. 1 Alternatively, [22,23] have proposed to find well-posed formulations of alternative theories of gravity extending the Müller-Israel-Stewart formalism of viscous relativistic hydrodynamics [24][25][26][27] to those theories of gravity. Very recently [28] succeeded in applying this formalism to theories of gravity with higher curvature corrections assuming spherical symmetry.
Treating the modifications to Einstein's gravity perturbatively may seem justified given that the present data indicates that they are small. In this case, there are no issues with the well-posedness of the equations and this is the approach that has been adopted in a number of papers [5,[8][9][10][11]29]. However, it has some serious limitations: it is well-known that small effects can accumulate over time and eventually lead to a breakdown of perturbation theory in a regime where it should still be valid. Furthermore, this approach is completely insensitive to certain non-perturbative effects encoded in the full non-linear theory. For instance, the non-linear perturbation theory around anti-de Sitter space breaks down precisely before a black hole forms [30].
GR is a classical theory and, as such, it should be understood as low energy effective field theory (EFT) of gravity. Indeed, on general grounds, one expects that at sufficiently small distances, Einstein's theory will be modified by quantum corrections. From the point of view of EFT, these corrections can be organised in a series expansion involving increasing powers of the curvature tensor, and consequently higher derivatives of the spacetime metric. Since in current experiments we are only probing gravity at low energies, we should only be sensitive to a finite number of terms in the otherwise infinite series of corrections to GR. Moreover, the details of the UV completion of gravity should not be important at such low energies. Higher derivative corrections are just one example of the myriad of possible modifications to GR that have been considered. Any of these alternative theories of gravity should be understood as truncated low energy EFT and, as such, they only make sense if the corrections to GR are small.
One particular modified theory of gravity which is known to have a well-posed initial value problem is Horndeski theory [17][18][19]. This is the most general theory of a metric tensor coupled to a scalar field with second order equations of motion arising from a diffeomorphism invariant action in four spacetime dimensions. 2 The general action for this theory is (1.2) where κ := 16πG is related to the 4-dimensional Newton's constant; φ is a scalar field and 3,4,5) are freely specifiable functions, and R and G µν are the Ricci scalar and Einstein tensor of the spacetime metric g µν , respectively. Having only second-order equations is essential to avoid Ostrogradsky instabilities [35,36]. This theory has found numerous applications to cosmology; the literature on the subject is vast and we will not attempt to review it here. We refer the reader to the recent reviews [37][38][39]. In this work we study the non-linear regime of a subclass of Horndeski theories for which [17] found a well-posed CCZ4 formulation of the Einstein equations. In this paper, unlike [8][9][10][11], we consider the theory in its full non-linear baroque splendour, which allows us to explore its distinctive non-perturbative physics; our goal is to identify the weakly coupled regime of the theory so that it can be consistently treated as a valid EFT from which one can obtain meaningful predictions. Rather than studying a specific phenomenologically viable theory, our ultimate goal is to identify general features in the waveforms that do not depend on the details and that can be attributed to the higher derivatives and non-linearities in the action. Therefore, we treat it as a toy model that can give us a glimpse of the type of effects that one can expect in more complicated theories which involve higher derivatives of the spacetime metric tensor.
For clarity of the presentation, we have split our work in a series of two articles, of which this is the first one. In this paper we study gravitational collapse and black hole formation in Horndeski theory. Our goal is to identify the region in the space of couplings for which the Horndeski theories under consideration are weakly coupled throughout the evolution. Using these results, in a companion paper we study black hole binary mergers, treating the theory fully non-linearly while remaining the regime of validity of EFT in all phases of the binary. In the following subsection, we summarise the main results in the present article, and refer the reader to the companion paper [40] for the results on black hole binaries.

Summary of the main results
In this paper we consider gravitational collapse in Horndeski theories using as initial data a spherically symmetric lump of scalar field (2.5). Even though the initial data is spherically symmetric, we evolve it using a 3+1 evolution code based on GRChombo [41], without symmetry assumptions. We have also considered gravitational collapse of some non-spherical scalar field configurations but we did not observe significant differences from the spherically symmetric case. However, a thorough study of gravitational collapse beyond spherical symmetry in Horndeski theories is beyond the scope of this paper.
Before we describe our results, we comment on previous works that are directly related to ours. Gravitational collapse and black hole dynamics in spherical symmetry in Einsteindilaton-Gauss-Bonnet (EdGB) theory has been studied before [42][43][44][45]. This theory can be considered to be a member of the Horndeski class, but the mapping between the two is highly non-trivial [46]. In these papers the authors study, among other things, the hyperbolicity of the equations of motion in various regions of the spacetime, including the interior of black holes, as a function of the coupling. They show that for large enough couplings the equations of motion can change character from hyperbolic to elliptic, even outside black holes, in which case one cannot solve them as an evolution problem. In a related work, [47] considers the conditions under which one may be able construct global solutions of Horndeski theories. In this paper, the authors study in detail the hyperbolicity of the equations of motion and the pathologies that may arise during the evolution in some specific examples. They also perform numerical simulations of spherically symmetric scalar field collapse to illustrate the breakdown of the hyperbolicity at strong coupling in different situations. Our work can be considered as an extension of these papers in different directions, as we now explain.
In this article we consider the so called cubic Horndeski theories (2.1), for which [17] showed that they have a well-posed initial value problem in the CCZ4 formulation of the Einstein equations and in puncture gauge. Because we are not particularly interested in a specific theory but rather in identifying general features of the non-linear dynamics of Horndeski theories, we consider two particularly simple and illustrative cases, see equation (2.4). In order for these theories to make sense as EFTs, the Horndeski terms have to be suitably small compared to the GR terms. Indeed, the well-posedness result of [17] only holds if a certain weak field condition is satisfied. For the class of theories that we consider, the relevant weak field conditions are given by (2.18). The main goal of this paper is to identify the region in the space of initial conditions and couplings for which the weak field conditions (2.18) are small at all times.
In our simulations of scalar field collapse we keep the radius r 0 and width ω of the initial Gaussian lump fixed, and vary both the amplitude A and Horndeski coupling (g 2 or g 3 depending on the theory under consideration). For every pair (A, g 2 ) or (A, g 3 ), we monitor both the character of the equations of motion of the scalar field 3 and the weak field conditions (2.18) everywhere in spacetime, except in a certain region of the interior of black holes when they form. It seems reasonable to accept the breakdown of EFT in a region sufficiently close to a singularity as long as this region is covered by a horizon. In this case, there is no loss of predictivity since this region is causally disconnected from the Universe outside the black hole, where EFT remains valid. The same criterion was adopted in [44].
Our main results for the G 2 = 0, G 3 = 0 theory are summarised in Fig. 1. The analogous figure for the G 3 = 0, G 2 = 0 theory is qualitatively similar and can be found in Section 3.2, Fig. 10. For the sake of definiteness, in the following we shall focus our discussion on the G 2 = 0, G 3 = 0 theory but essentially the same conclusions apply to the G 3 = 0, G 2 = 0 theory.
The dimensionless coupling constants η 2 and η 3 , see eqs. (2.8)-(2.9), control the future development for our initial data; in other cubic Horndeski theories one should be able to define analogous dimensionless couplings, and therefore the conclusions of this paper should apply to those theories as well. In Fig. 1 we show the various dynamical regimes of the G 2 theory as a function of η 2 and the initial scalar amplitude A. As one would expect, the weakly coupled regime of the theory corresponds to suitably small values of η 2 , but the boundary of this region depends non-trivially on the scalar amplitude.
The blue region in Fig. 1 denotes the values of (A, η 2 ) for which the scalar equation is hyperbolic at all times. The yellow region corresponds to the values of (A, η 2 ) for which the scalar equation is initially hyperbolic, and hence the initial value problem is well-posed, but it changes character during the evolution, signalling a breakdown of the theory. The green region corresponds to the values of (A, η 2 ) for which the scalar equation is not hyperbolic on the initial data slice and hence the initial value problem is not well-posed. The black band in Fig. 1 Figure 1: Dynamical regimes of the G 2 = g 2 X 2 theory as a function of the initial amplitude A and the dimensionless coupling constant η 2 , see eq. (2.8). The black band denotes the region near critical collapse; black holes form to the right of this band. The orange curve on the right marks the region where the initial data contains a trapped surface. The scalar equation is hyperbolic at all times in the blue region; EFT is valid in the interior of this region. In the yellow region, the scalar equation is initially hyperbolic but it changes character during the evolution. In the green region the initial value problem is not well-posed.
initial data gets close to Choptuik's critical solution [48], which is a naked singularity. This band splits the figure into two regions corresponding to the small and large data regimes: for initial data in the blue region to the left of the black band, the scalar field disperses to infinity. On the other hand, initial data in the blue region to the right of the black band collapses into a black hole. For initial data in any of blue regions in Fig. 1 it is possible to construct global solutions to the G 2 = g 2 X 2 Horndeski theory. Away from the boundary of this region, the deviations from GR are "small" everywhere on and outside black holes (if there are any) for all times. By "small" here we mean that the weak field conditions (2.18) are satisfied. Therefore, we identify the interior of the blue region as the regime of validity of EFT for the corresponding Horndeski theory. Of course, when black holes form during the evolution, EFT will break down near the singularity, just as GR does. In this case, we excise a portion of the interior of the black hole since it is causally disconnected from external observers. For values of (A, η 2 ) close to the boundary of the blue region, the weak field conditions (2.18) can become O(1) during the evolution while the scalar equations remain hyperbolic. In this case one may argue that even though the theory has a well-posed initial value problem, higher derivative corrections not included in the action (2.1) should become important and hence one should not trust the theory as it stands.
As mentioned above, the yellow region in Fig. 1 denotes the values of (A, η 2 ) for which the evolution breaks downs due to the change of character of the scalar equation and this breakdown cannot be hidden behind a horizon. The change of character of the scalar equation is typically associated to the weak field conditions becoming O(1) or larger but this is not always the case. Indeed, for certain values of (A, η 2 ), and in particular for η 2 < 0, the weak field conditions can be reasonably small during the evolution and yet the equations change character. Beyond this point it is no longer possible to solve the theory as an initial value problem. In Section 3 we study in detail how and where in spacetime the loss of hyperbolicity of the scalar equation happens depending on the Horndeski couplings and we correlate it to the weak field conditions (2.18). For 0 < A 0.05, the boundary between the blue and yellow regions is given by a constant value of η 2 ∼ (1.44 ± 0.06) × 10 −4 and η 2 ∼ (−1.94±0.07)×10 −5 respectively. This is non-trivial since the location of this boundary is obtained from the non-linear evolution of the initial data. As we will see in Section 3.1, for η 2 > 0 the breakdown of the evolution happens through a Tricomi-type-of transition while for η 2 < 0 the transition is of the Keldysh type. Fig. 1 indicates that η 2 → 0 as one approaches the critical regime from both sides. This is expected since for A near the critical amplitude A * = 0.13 ± 0.01, the gradients of both the metric tensor and the scalar field become very large as the solution approaches the critical solution, which leads to a change of character of the scalar equation unless g 2 → 0 as A → A * . Since the regime of validity of EFT is essentially the empty set at the critical solution, in the rest of the paper we will purposely avoid the region near criticality. For values of A > A * , a black hole forms during the evolution of the initial data. The larger the value of A, the larger the black hole that forms and the sooner it forms. Since larger black holes result in lower curvatures on the horizon scale, larger values of the couplings are allowed and yet the theory remains weakly coupled on and outside the black hole. This is the reason why η 2 increases for larger A. For sufficiently large A, the initial data already contains a trapped surface. Since we are interested in studying gravitational collapse, we do not consider those values of A.
It is clear from the previous discussion that our weak field conditions (2.18) bear some relation with the hyperbolicity condition of the scalar equation of motion (2.3) but such a relation is not a direct one. It is possible that one can come up with refined and sharp weak field conditions that also capture the change of character of the equations when they are violated but finding them is beyond the scope of this paper. It follows from our analysis that the regime of validity of EFT corresponds to the weak field conditions (2.18) being satisfied (to justify that higher derivative terms in (2.1) can be neglected) and that the initial value problem is well-posed, i.e., the scalar equations of motion are hyperbolic everywhere in spacetime, perhaps except in a small region inside black holes. These two conditions are satisfied in the interior of the blue region in Figs. 1 and 10. For initial data in this region, the Horndeski theories that we have considered are valid EFTs and global solutions can be constructed.
The rest of the paper is organised as follows. In Section 2 we present the theories that we consider and we analyse the corresponding hyperbolicity conditions. In Section 3 we present and analyse the results of our numerical simulations. Subsection 3.1 discusses in detail the dynamics of the G 2 = 0 theories, while the G 3 = 0 theories are dealt with in Subsection 3.2. We conclude with some final remarks in Section 4. We have relegated some technical details to the Appendices. In Appendix A we write down the equations of motion for scalar field and the effective scalar metric in a 3+1 form. We collect some technical results in Appendices B and C, and the convergence tests are presented in Appendix D. Appendix E contains the results of certain numerical simulations that are also relevant for the main text. In this paper we adopt the following notation. We use Greek letters (µ, ν, ρ, ...) to denote full spacetime indices and Latin letters (i, j, k, ...) for the spatial ones. We adopt the mostly plus metric signature, and we set G = c = 1.

Equations of motion
In this paper we consider the special subset of Horndeski theories for which [17] proved well-posedness of the initial value problem in both the BSSN and CCZ4 formulations of the Einstein equations in the usual gauges used in numerical relativity. This class of theories is given by setting G 4 = G 5 = 0 in the general Horndeski action (1.1). This results in the so called cubic Horndeski theories described by the action Here, X = − 1 2 (∇ µ φ)(∇ µ φ) and V (φ) are the usual kinetic and potential terms respectively in the standard action for a minimally coupled scalar field, and G 2 (φ, X) and G 3 (φ, X) are arbitrary functions of their arguments. In this paper, we have explicitly separated the canonical kinetic and potential terms from G 2 so that G 2 and G 3 parametrise the higher derivative terms and non-minimal couplings of the scalar field to gravity. The resulting Einstein equations are: where G µν is the Einstein tensor. The equation of motion for the scalar field is: 4 We write down equations (2.2) and (2.3) in the usual 3+1 conformal decomposition and implement the CCZ4 form of the Einstein equations that is suitable for the numerical simulations. The equations that we have implemented in our code as well as the details of the numerical simulations are given in Appendix A. In the remainder of this Section, we describe the specific cubic Horndeski theories that we have studied, our initial data, the analysis of the hyperbolicity of the scalar equations and the weak field regime.

Cases explored
The action (2.1) comprises several well-known particular cases that have been extensively studied in other contexts, mostly cosmology (see [37,38,49]). For instance, quintessence, which consists of a simple scalar field minimally coupled to GR; this model is obtained by setting G 2 = G 3 = 0 in (2.1). On the other hand, models of k-essence are obtained by setting G 3 = 0 in (2.1), with the common choice of G 2 (φ, X) = f (φ)g(X) for arbitrary functions f and g of their arguments. Finally, kinetic gravity braiding [50], also referred as Cubic Galileons [51,52], are obtained from (2.1) by choosing G 3 = 0; this class of models is often simplified to the shift symmetric case, corresponding to G 3 (φ, X) = g(X), for an arbitrary function g. Therefore, the subclass of Horndeski theories that we consider is very rich and has multiple applications to gravitational physics and cosmology.
In our work we are not interested in a particular model but rather in exploring general features of the non-linear physics encoded in cubic Horndeski theories. From the point of view of EFT, one would expect (2.1) to be valid when the G 2 and G 3 terms are suitably small, which corresponds to X being small. Therefore, one can consider Taylor-expanding some general (smooth) functions G 2 and G 3 for small X and keep only the leading order terms. With this in mind, we therefore focus on the simplest non-trivial functions G 2 and G 3 : where g 2 and g 3 are arbitrary coupling constants with dimensions of Length 2 that we can tune. These or similar choices have been considered in the literature before, namely in models of dark energy [50,[53][54][55][56][57], and in studies of the fate of the Universe in cosmological bounces or inflationary models [56,[58][59][60], among others [61,62].

Initial data
For the present analysis, motivated by the objective of studying gravitational collapse, we choose a family of initial data for the scalar field (φ, Π) modelling a spherically symmetric bubble centred at c: where r = x − c and r = || r|| 2 with the Euclidean 2-norm. Notice that the class of theories in (2.4) have a reflection symmetry φ → −φ , g 3 → −g 3 and hence, we can choose A > 0 without loss of generality. Regarding the scalar momentum, assuming an approximately Minkowski initial background, we choose an ingoing wave pulse: To explore the relevant phenomenology of these theories, we have studied many different scenarios. Using a full 3D code, we were able to verify that all the features hereafter described are not a peculiarity of spherical symmetry, and also occur when the symmetry is broken, without any seemingly interesting new features emerging. However, we have not attempted to carry out a thorough analysis of non-spherically symmetric scalar field collapse. Hence, in the following we only present the results for the spherically symmetric case.
With the choices (2.5) and (2.6) for the initial scalar profile and momentum, we obtain the initial data for the metric by solving the Einstein constraints using the conformal transverse-traceless decomposition [63,64]. We choose a conformally flat initial metric and vanishing trace and transverse-traceless part of the extrinsic curvature. Hence, we solve for the conformal factor of the spatial metric and three leftover degrees of freedom of the traceless part of the extrinsic curvature (which reduce to one in spherical symmetry).
To get some intuition about how the modifications of GR affect our initial data, we can expand the initial ADM mass for small amplitudes and couplings around a Minkowski background. We find, and we have included the contribution of a mass term in the scalar potential V (φ) = 1 2 m 2 φ 2 . From (2.7) we see that for our initial data, the strength of the modifications of GR due to the Horndeski terms is measured by the dimensionless couplings: These dimensionless couplings play an important role in the future development of the initial data and determine the weakly coupled regime of these theories.
In Fig. 2 we show the initial conformal factor χ and scalar profile φ for some representative cases. From this figure we see that even for relatively large amplitudes within the range that we have considered, the conformal factor has a very small dependence on the Horndeski couplings. For the specific case of A = 0.22, the difference between g 2 = 1.5 and GR at r = 0 is 0.2%, which is in accordance with the fact that for this case the dimensionless coupling is small (η 2 ∼ 3 × 10 −3 ). One can also notice that for sufficiently small amplitude, the conformal factor is almost 1 for any reasonable value of g 2 .

Effective metric and characteristic speeds
To identify the regime of validity of EFT, we need to first determine the character of the equations of motion for the scalar field (2.3) and the conditions under which they are hyperbolic. To do so, we consider the principal part of the scalar equation (2.3), which is a wave equation governed by an effective metric (2.10) The eigenvalues of h µν determine the character of the equation: if the product of the eigenvalues is negative then equation is hyperbolic; if the product is positive then the equation is elliptic, and if it is zero the equation is parabolic.
Having a well-posed initial value problem is the minimum requirement that we should demand on any classical theory; therefore, the breakdown of hyperbolicity of the scalar equation in this case can be associated to the breakdown of the theory itself. As [47] noted, the fact that the effective metric (2.10) depends on the scalar field itself and its gradients implies that shocks can generically form from smooth initial data, at which point the wellposedness is lost. Therefore, the local character of the scalar equation is a useful proxy to establish the regime of validity of the theory [42][43][44]47] and to measure the size of the non-linearities and deviations from GR. We will come back to this point below.
When considering spacetimes containing black holes, the evolution of the spatial slices in puncture gauge is such that the determinant of the inverse spacetime metric goes to zero near the puncture, i.e., det(g µν ) = − χ 3 α 2 → 0 (see Appendix A). Consequently the same happens for the effective metric (2.10). To distinguish this gauge effect from an actual breakdown of the hyperbolicity of the scalar equation, we note that h µν = g µρ h ν ρ and therefore: with det (h µ ν ) = 1 in GR. Clearly, deviations of this quantity from 1 encode the dynamics of the Horndeski theories and hence we will focus our attention on det (h µ ν ). The characteristic speeds, also called front velocities, are important since they correspond to the local speed of propagation of the scalar modes and hence they tell us about the effective causal cone that the scalar field "sees". 5 The characteristics are given by the zeros of the characteristic polynomial which, for the scalar field equation, is for some covector ξ µ that defines the characteristic surface. Physically this corresponds to considering the high frequency and small amplitude limit of a wave with wave vector ξ µ . To calculate the propagation speeds without symmetry assumptions, we specify a direction of propagation, n i suitably normalised n i n j δ ij = 1, where δ ij is the Euclidean 3D metric (since the space is locally flat). Then, the speed of propagation in the n i direction is: In spherical symmetry one can naturally use a radial vector for the direction of propagation, (2.14) In our conventions, v + and v − correspond to the ingoing and outgoing modes of the scalar field respectively and they are normalised so that they tend to +1 and −1 at infinity. When v − ≥ 0 and v + ≥ 0 in a certain region, scalar modes cannot reach asymptotic observers; the boundary v − = 0 of this region is the sound horizon [65].
As discussed in [43,47], the equations can change character from hyperbolic to parabolic and elliptic in a manner which is qualitatively similar to what happens in the two standard equations of mixed type, namely the Tricomi equation, 15) and the Keldysh equation, Both equations are hyperbolic for y < 0 and they change character at the transition line y = 0. Related to this change of character are the appearence of ghosts, gradient instabilities and formation of caustics [66][67][68].
In the case of the Tricomi equation, the characteristic speeds go to zero at y = 0 where the equation becomes parabolic, while in the Keldysh equation the characteristic speeds diverge at y = 0. If the characteristic speeds of both the ingoing and outgoing modes vanish then the evolution freezes. This can happen because of the choice of gauge; for instance, in coordinates that are not horizon penetrating, the lapse asymptotically goes to zero at the horizon, effectively resulting in zero characteristic speeds. However, in this case the freezing of the evolution is a consequence of the gauge choice and it does not correspond to a breakdown of EFT. Therefore, in the case of a Tricomi-type-of transition, we also need to check that the deviations from GR are suitably large to conclude that the loss of hyperbolicity corresponds to a breakdown of the theory. On the other hand, a Keldysh-type-of transition involves diverging characteristic speeds 6 , which will typically signal a breakdown of EFT. This case is more difficult to handle numerically since one is forced to take prohibitively small time steps. 7 Note from (2.10) that h 00 has a factor of −1/α 2 coming from g 00 , and hence the deviations from GR are measured by −α 2 h 00 . Therefore, the Keldysh-type-of transition without symmetry assumptions is signalled by −α 2 h 00 → 0, which implies that the t = const. hypersurface being evolved is no longer spacelike with respect to the scalar effective metric [69]. We associate this breakdown of the evolution to a Keldysh-type-of transition since the characteristic speeds diverge. However, strictly speaking, at this point the equation may not have changed character yet but the two effects go hand in hand. 8 We discuss in detail the different types of transitions in the G 2 = 0 and G 3 = 0 cases in the next subsection.
The previous discussion only relates to the existence of a well-posed initial value problem but it does not fully address the issue of whether the theory under consideration makes sense as a truncated EFT [70]. We now turn to this point. As mentioned in [17], local wellposedness is only guaranteed in the weak field regime, meaning that the Horndeski terms are small compared to GR ones. One possible weak field condition that compares the size of the Horndeski terms versus GR is: where L is a length scale estimate for the system: L −1 = max{|R αβµν | 1 2 , |∇ µ φ| , |∇ µ ∇ ν φ| 1 2 } in all orthonormal bases. For the cases (2.4), this is explicitly: In order for the Horndeski theories under consideration (2.4) to be in the regime of validiy of EFT, in this paper we require that the evolution equation of the scalar field is hyperbolic and that (2.18) is satisfied. These two conditions ought to be imposed on and outside black hole horizons, should there be any in the spacetime.

Case of
To monitor the character of the scalar equation, we compute the determinant of the scalar effective metric. Even though it is possible to find an analytic expression for the full determinant (using Cayley-Hamilton's theorem and Newton's identities), for simplicity we consider the G 2 = 0, G 3 = 0 and the G 3 = 0, G 2 = 0 cases separately. As explained in the discussion surrounding eq. (2.11), we only need to consider the determinant of the effective metric with one index up and one index down, which is significantly simpler. For the G 2 = 0 case, we have Realising that, up to scalars, this metric is the identity plus the tensor product of two vectors, one can use the Weinstein-Aronszajn identity to calculate the determinant of the full 4D 8 We would like to thank Luis Lehner for discussions on these issues. metric without assuming any symmetries. We find: , (2.20) where in the last line we have used that G 2 = g 2 X 2 . We can compute the eigenvalues and eigenvectors v ν by noting that so we conclude that ∇ µ φ is an eigenvector with eigenvalue (1 + ∂ X G 2 + 2 X∂ 2 XX G 2 ). The other three eigenvectors are orthogonal to the 4-vector ∇ µ φ and have degenerate eigenvalues equal to (1 + ∂ X G 2 ), in accordance to (2.20).
To monitor a Keldysh-type-of transition, we have to compute −α 2 h 00 . For the G 2 = 0, G 3 = 0 case, this is given by, where Π = n µ ∇ µ φ is the scalar momentum, and in the last line we have used that . All in all, for G 2 as in (2.4), the two quantities that inform us about the breakdown of the initial value problem for the scalar equation are (2.20) and (2.22). The scalar equation is hyperbolic as long as these two quantities are non-negative. 9 These are the same conditions found in [69], and it is evident that if the weak field conditions (2.18) are satisfied then the scalar equation is hyperbolic. For a non-constant scalar profile φ, Π i Π i > 0 but X can be either positive or negative, depending on the balance between scalar gradients and momentum. In a dynamical evolution, both can become large. As this happens, g 2 X can decrease to make either (2.20) or (2.22) zero, see Fig. 3. If g 2 > 0, the fact that Π i Π i > 0 implies that det (h µ ν ) will reach zero before −α 2 h 00 , and the equation will become parabolic on a co-dimension one surface, where at least one of the characteristic speeds goes to zero while the others remain bounded. This will correspond to a Tricomi-type-of transition. On the other hand, if g 2 < 0 the opposite is true and −α 2 h 00 may become zero before det (h µ ν ) does, leading to infinite speeds of propagation and a very abrupt termination of the evolution associated to a Keldysh-type-of transition. Both behaviours were identified in [47]. 10 The changes of character described in the previous paragraph can only occur if |g 2 X| is suitably large and hence outside the weak field regime. While generically one can expect that weak data eventually enters the strong field regime, one question that we need to address is whether or not the region where EFT breaks down can be hidden inside a black hole. If the answer is positive, then one can hope that classical observers at infinity will be protected from any potential pathologies that arise in the scalar equations and EFT will retain its 9 Note that we have pulled out a minus sign in (2.11) so det (h µ ν ) > 0 corresponds to h µν having one negative eigenvalue and three positive ones, as it should for a hyperbolic equation. 10 Reference [47] uses a coupling g with the opposite sign as our g 2 .

Case of
In equation (B.1) of Appendix B we present the full analytic form of the determinant of the scalar effective metric in the G 3 = 0, G 2 = 0 case. For clarity, in this subsection we analyse (B.1) for small g 3 , which is the relevant limit in the weak field regime.
To obtain the expansion of (B.1) for small g 3 , we use the scalar equation of motion (several times if necessary) to replace φ in (B.1) by V (φ) and terms which are higher order in g 3 , in the spirit of order reducing schemes. We then obtain, up to second order: Let us focus on the case of zero scalar potential, V (φ) = 0, which is the relevant one for this paper. In this case, the correction to GR in det(h µ ν ) comes at O(g 2 3 ), while in −α 2 h 00 it comes at leading order. Because τ does not have a definite sign, then regardless of the sign of g 3 , there will be regions in spacetime where −α 2 h 00 will vanish before det (h µ ν ) Figure 4: Evolution of the profile of a (massless) scalar field in GR for the Cases 1-4.
does, resulting in a Keldysh-type-of transition. This should be the generic behaviour in the V (φ) = 0 case for the G 3 = 0, G 2 = 0 theory, and it is indeed what we observe in our numerical simulations, see Section 3.2. The picture changes for V (φ) = 0; then det (h µ ν ) receives a contribution to leading order in g 3 and the type of transition will depend on the details of the scalar potential and the initial data.

Numerical results
In this section we present the results of our numerical simulations of the gravitational collapse of a single massless scalar bubble with intial data as in Section 2.3. In all our simulations we keep the radius r 0 and the width ω of the initial Gaussian profile (2.5) fixed, and we vary both the amplitude A and Horndeski coupling g 2 or g 3 . The reason is that varying r 0 and ω leads to similar results and varying A alone makes the analysis simpler. We choose r 0 = 5 and ω = √ 0.5, which set the length scale in our simulations. Since we consider spherically symmetric scalar field collapse (even though we do not assume spherical symmetry in our simulations), there are essentially two relevant regimes depending on whether the initial data disperses to infinity (small data) or it collapses into a black hole (large data). We consider four representative values of the initial amplitude A, so that we can probe the regimes far and close to critical collapse for both small and large data: 11 • Case 1: A = 0.05 -dispersion far from the critical regime, with initial ADM mass of ≈ 0.022.
• Case 3: A = 0.22 -collapse into small black hole with initial ADM mass of ≈ 0.5.
• Case 4: A = 0.33 -collapse into a larger black hole with initial ADM mass of ≈ 1.6.
For each of these cases, we vary the Horndesky couplings (g 2 or g 3 ) while ensuring that the initial value problem is well-posed. We then evolve the initial data by solving the coupled equations of motion (2.2)-(2.3) numerically, and we monitor both the hyperbolicity of the scalar equation and the weak field conditions (2.18). In this way we can identify the regime of validity of the EFT for both small and large data. We shall refer to the different cases as "weakly" or "strongly" coupled depending on the whether the hyperbolicity of the scalar equation breaks down at some point during the evolution; this breakdown is typically associated to the weak field conditions (2.18) becoming large but, as we shall see, this need not always be the case. The evolution of the scalar field in GR (i.e., g 2 = g 3 = 0) for the Cases 1-4 is shown in Fig. 4. It is worth emphasising that Cases 2 and 3 above do not exhibit Choptuik's critical behaviour, as the amplitude A is purposely chosen to be sufficiently 'far' from the critical amplitude A * ≈ 0.13 ± 0.01. The reason is that Choptuik's critical solution is a naked singularity and EFT will necessarily break down close to it. Indeed, it is apparent from Figs. 1 and 10 that the coupling constants have to be tuned down to maintain the hyperbolicity of the scalar equation as we approach the critical regime from both sides. In addition, the weak field conditions (2.18) become large the closer we get to the critical solution, as expected.
Since the G 2 = 0 and G 3 = 0 theories do not exhibit significant qualitative differences in terms of the dynamics of collapse of the scalar field, in the next subsection we will focus the discussion on the G 2 = 0 theory considering different values and signs of the coupling constant g 2 . In subsection 3.2 will only highlight the main differences in the G 3 = 0 case.

G 2 theories
In the following subsections we will discuss gravitational collapse in Horndeski theories with G 2 = g 2 X 2 for different values of the coupling constant g 2 . For our scalar field initial data, during collapse a positive and negative peak in X form; these peaks grow as the evolution progresses and the scalar shell approaches the origin. After reaching the origin, they bounce back and smaller peaks of opposite signs form, eventually resulting in the formation of a black hole or dispersion to infinity. See Fig. 4 for the evolution of the scalar field profile in GR; in the Horndeski theories it is qualitatively similar. With an initial ingoing momentum, as in our initial data, momentum dominates over spatial gradients and the positive peak will be much larger in amplitude than the negative or the subsequent peaks that form after the bounce. Considering the expressions for det(h µ ν ) and −α 2 h 00 in (2.20) and (2.22) for the G 2 = 0 theory, this implies that generically a negative g 2 will lead to a breakdown of the hyperbolicity of the equations for a significantly smaller |g 2 | and it will happen sooner than for a positive g 2 . Furthermore, as described in Section 2.4.1, for g 2 < 0 the change of character will be of the Keldysh type while for g 2 > 0 it will be of the Tricomi type.

Weak coupling
We first consider the case of a small and positive coupling constant g 2 ; we choose g 2 = 0.005 as a representative example. This is a case of a theory that remains in the regime of validity of EFT throughout the whole evolution, both for small and large initial data. For this choice of parameters, the maximum of the weak field conditions (2.18) is small everywhere on and outside horizons (if they form) at all times.
In Fig. 5 we display det(h µ ν ) (top) and −α 2 h µ ν (bottom) for Cases 1-4. The white lines in Figure 6: Outgoing (top) and ingoing (bottom) scalar characteristic speeds for g 2 = 0.005. The evolution freezes inside black holes as a consequence of the 1 + log slicing condition that we use. Sounds horizons form for large initial data.
these plots indicate the trajectories of the initial scalar field peak and serve to guide the eye. In Cases 1 and 2, the scalar field bounces at the origin and eventually disperses to infinity; as the amplitude increases from Case 1 to Case 2, the scalar field spends more time near the origin where gravitational focusing is stronger. For sufficiently large amplitudes (Cases 3 and 4) it collapses into a black hole. In all cases, both det(h µ ν ) > 0 and −α 2 h 00 > 0 throughout the evolution so the scalar equations are hyperbolic at all times. The long dashed line in Fig.  5 indicates the contour where the maximum of weak field conditions (2.18) is equal to one; as we can see, for Cases 1 and 2 the weak field conditions are always less than one everywhere in spacetime, while in Cases 3 and 4 and the weak field conditions are greater than one only inside the apparent horizon (solid black line). Only in Case 3 there is a small region near the origin where the weak field condition is greater than one and for a short period of time is not covered by an apparent horizon. Note however that this region is already cloaked by the sound horizon (dotted black line), so the scalar modes emanating from this region cannot reach asymptotic observers. For Cases 2-4, det(h µ ν ) can significantly deviate from 1 (its GR value) when the scalar field is most contracted at the origin. Likewise, the bottom plots in Fig. 5 show that −α 2 h 00 also exhibits some deviation from its GR value near the origin but it never gets anywhere close to 0. Therefore, despite the weak field conditions being small at all times, the Horndeski terms can have a significant impact on the dynamics of the system, especially near the origin where the gravitational focusing is strongest.
In Fig. 6 we display the characteristic speeds for both the outgoing (top) and the ingoing (bottom) modes. Notice that in Case 2, both speeds approach zero at the origin when the scalar field collapses but their sign does not change. This is indicative of strong gravitational dynamics, as one would expect since Case 2 is "close" to the critical regime. Also, note that there are no scalar horizons in this case and all the scalar field eventually disperses to infinity. The dynamics changes in Cases 3 and 4, where a black hole forms. First, notice that v − changes sign inside the black hole, from negative to positive; this implies that inside the black hole, outgoing modes travel inwards, as expected. Eventually both speeds become close to zero in the region near the singularity. This is just a consequence of using 1 + log slicing in our simulations, which effectively freezes the evolution inside black holes. Second, we do observe the formation of scalar horizons, where v − = 0 and v + > 0. In both Cases 3 and 4, the characteristic speed of the outgoing modes is small in the vicinity of the sound horizon; consequently, even though the scalar field can eventually reach infinity, it will remain near the black hole for a long time, thereby interacting with itself and with the black hole.
It is apparent from the results shown here that even though the weak field conditions (2.18) are small everywhere, the scalar field still exhibits strong dynamics, such as the dynamical formation of scalar horizons. The latter is a non-perturbative effect and it that can only be seen if one treats the Horndeski theory fully non-linearly. Evidently, if the couplings are small then the scalar horizon will be close to the metric horizon. In the case of a black hole binary in a Horndeski theory of gravity, even if the effects of the strong scalar dynamics are locally small, over a sufficiently long time they can lead to significant deviations from GR that may be observable [40].

Strong coupling
In this subsection we analyse the case for which the g 2 coupling is large and positive. We choose g 2 = 1.5 as a representative example. For this value of the coupling constant, the weak field conditions (2.18) can be O(1) for large initial data, see Fig. 7 Cases 3 and 4. Therefore, strictly speaking, in these cases the theory is already outside the regime of validity of EFT even though the initial value problem is well-posed. Nevertheless, we choose this value of the coupling constant as an illustrative example of the dynamics of Horndeski theories for large and positive g 2 .
In Fig. 7 we display det(h µ ν ) and −α 2 h 00 during the evolution for our four cases. Unsurprisingly, this figure shows that in all cases the evolution breaks down at some point.  is because det(h µ ν ) → 0 in a certain region at some instant of time and hence the scalar equation changes character, becoming parabolic. Beyond this point and it is not possible to solve the equations as an initial value problem. For this value of the Horndeski coupling, for all Cases 1-4 the weak field conditions (2.18) have become large before the equations change character. Also, note that for large initial data (Cases 3 and 4), the evolution breaks down before an apparent horizon has had time to form and hence the pathology in the scalar equations of motion cannot be hidden behind the horizon. Fig. 7 (bottom) shows that in all cases −α 2 h 00 deviates significantly from its GR value and but remains well above zero up until the breakdown of the evolution. Likewise, we observe that in these simulations the characteristic speeds of both the ingoing and outgoing modes remain bounded at all times; in fact, both v − → 0 and v + → 0 as det(h µ ν ) → 0. Therefore, the loss of hyperbolicity for the g 2 > 0 theories is due to a Tricomi-type-of transition, in accordance with the discussion in Section 2.4.1. By lowering the coupling constant a bit, it is possible to hide the strong scalar field dynamics that causes the breakdown of the hyperbolicity of the equations inside a large enough black hole. This is illustrated in Appendix E, Fig. 13. For such "intermediate" couplings, the evolution still breaks down in Cases 2 and 3, while in Case 4 the pathologies that develop in the scalar equation can be hidden behind the horizon. In this case, one can continue the evolution without encountering any issues. Moreover, the weak field conditions in Case 4 remain small on and outside the black hole horizon despite the fact that g 2 is large. Clearly, from the expression for the dimensionless coupling η 2 , eq. (2.8), one can achieve the same results by increasing the initial amplitude A instead of decreasing g 2 .

Negative coupling
In this subsection we discuss the case of a strong and negative coupling constant g 2 . As an illustrative example, we consider g 2 = −0.2.
As anticipated in Section 2.4.1, the dynamics of the scalar field changes quite significantly for negative couplings. First, a smaller absolute value of g 2 is enough to cause a breakdown of the hyperbolicity of the scalar equations for both small and large data. The results are shown in Figs. 8 and 9. In all cases we find that −α 2 h 00 → 0 before det(h µ ν ) → 0, even though this is not easily seen from Fig. 8. This implies that, in our gauge, the t = const. surfaces are no longer spacelike with respect to the scalar effective metric before the scalar equation changes character. The fact that for g 2 < 0, −α 2 h 00 → 0 first results in infinite characteristic speeds of propagation for both the ingoing and outgoing modes, see Fig. 9. Therefore, we associate the breakdown of the hyperbolicity of the scalar equation to a Keldysh-type-of transition, in accordance to the discussion in Section 2.4.1 and [47]. The diverging characteristic speeds near the transition point imply that the dynamics of the scalar field becomes increasingly fast right before it breaks down; to adequately resolve it, in our simulations we had to significantly reduce the Courant factor. However, at some point it is no longer feasible in practice to keep reducing it, and numerical errors eventually build up until the simulation inevitably crashes. A possible way out would be to change our slicing conditions to ensure that the t = const. hypersurfaces remain spacelike with respect to both h µν and g µν , but we have not attempted to do so here.
Notice that for this value of the coupling constant, the weak field conditions (2.18) are always less than one everywhere in spacetime, including the region near the origin where gravitational focusing is strongest, except immediately before the breakdown. This is simply a consequence of the fact that |η 2 | is small in all Cases 1-4. Related to this last observation, the breakdown occurs before either sound horizons or apparent horizons have had time to form, so the pathologies cannot be hidden from asymptotic observers. We expect that a refined weak field condition should be able to capture that this case also becomes strongly coupled before the breakdown of the evolution.
Needless to say, for sufficiently small absolute values of |g 2 | the scalar equations remain hyperbolic at all times for Cases 1-4. In this situations the evolution is qualitatively similar to the small and positive g 2 case that we have already discussed in Subsection 3.1.1. Likewise, for a given g 2 < 0 and a sufficiently large A, the pathologies in the scalar equation can be hidden inside the black hole horizon.

G 3 theories
In this subsection we will briefly comment the dynamics in Horndeski theories with G 3 = g 3 X and G 2 = 0. In all cases that we have explored, either for g 3 > 0 or g 3 < 0, the dynamics is qualitatively similar to the G 2 = g 2 X 2 theories with g 2 < 0 so we will not go into much detail.
As discussed in Section 2.4.2, we expect that for sufficiently small absolute values of g 3 , the breakdown of the scalar evolution equations would be due to a Keldysh-type-of transition. Figure 9: Characteristic speeds of the outgoing (top) and ingoing (bottom) scalar modes for g 2 = −0.2. Both speeds diverge, even though v + does so slower, when the evolution breaks down, in accordance with a Keldysh-type-of transition.
Our numerical simulations confirm that this is indeed the case for either signs of g 3 . In Figs. 14 and 15 of Appendix E we show the results for a representative case with g 3 = 0.4. In Fig. 14 we see that −α 2 h 00 → 0 before det(h µ ν ) does, resulting in infinite characteristic speeds, as expected in a Keldsyh-type-of transition. In this case we observe that v − diverges as −α 2 h 00 → 0 while v + remains finite, see Fig. 15. Note that in this particular example, for large data (Cases 3 and 4) the evolution breaks down before the first apparent horizon appears. However, just as in the G 2 = 0 theories, either by increasing the initial scalar amplitude so that a sufficiently large black hole forms or by lowering |g 3 |, it is possible to hide the pathologies that may arise in the scalar evolution inside a black hole so that the theory remains in the regime of validity of EFT on and outside black holes. This is precisely what happens in the interior of the blue region in Fig. 10 in the large data regime.  In the yellow region, the scalar equation is initially hyperbolic but it changes character during the evolution. In the green region the initial value problem is not well-posed.
as in Fig. 1 and the qualitative features are also the same. The black band corresponds to the range of A for which the future development of initial data becomes close to Choptuik's critical solution. Black holes form for A to the right of the black band while for A's to left, the scalar field disperses. As before, global solutions to this particular Horndeski theory can be constructed for values of (A, η 3 ) in the blue region. The regime of validity of EFT corresponds to the interior of the blue region, away from its boundaries. For 0 < A 0.05, the boundary between the blue and yellow regions is at a constant value of η 3 given by η 3 ∼ (9.20 ± 0.09) × 10 −6 and η 3 ∼ (−1.04 ± 0.03) × 10 −5 respectively.

Final remarks
In this paper we have studied the regime of validity of certain cubic Horndeski theories of gravity that have a well-posed initial value problem. We have chosen two particularly simple cases, namely (2.4), but we expect that our results should extend to other models as well, at least in the weak coupling regime which is where these theories should be valid EFTs. For instance, for a single massive scalar field the results are qualitatively unchanged during gravitational collapse. Nevertheless, one expects that a massive scalar field will stay trapped around the black hole for a much longer time, forming scalar clouds [71]. Over long periods of time, such as in a black hole binary inspiral, the locally small deviations from GR introduced by Horndeski theories can (and will!) accumulate, giving rise to significant deviations. For the particular class of models that we have studied, the reason why the evolution breaks down is because the scalar equation changes character. For the G 2 = g 2 X 2 theory the transition can be of the Tricomi type for g 2 > 0, while for g 2 < 0 the transition is of the Keldysh type. On the other hand, for the G 3 = g 3 X theory, we have only observed a breakdownà la Keldysh. However, this is not generic for the G 3 theories; other choices such as G 3 = g 3 X 2 can exhibit both behaviours. Furthermore, we have provided some level of analytic justification for the types of pathologies that may arise in each of the models that we have considered.
In order for the initial value problem be well-posed and the theory be a consistent (truncated) EFT, we need to impose that a certain weak field condition (2.18) is suitably small. For certain choices of initial conditions and couplings (no fine-tunning required) the conditions in (2.18) can be O(1) and yet the scalar equation of motion is perfectly hyperbolic. Conversely, the conditions in (2.18) can be small and yet the scalar equation changes character. Therefore, it would be very interesting to obtain a sharp condition that identifies the truly weakly coupled regime of the theory and provides some analytic understanding of it, at least for certain classes of initial data.
Having identified the regime regime of validity of the Hordneski theories that we have considered, we can proceed to study black hole binaries for initial data in this regime. These studies will be presented in the companion paper [40].

A.1 Equations of Motion
To carry out the numerical simulations presented in this paper, we used a code based on GRChombo, a multipurpose numerical relativity code [41] 12 that implements the BBSNOK [72][73][74] or CCZ4 [75][76][77][78][79] formulations of the Einstein equations. In this appendix we present the conformal 3+1 form of the stress tensor and the scalar equation (2.3) as we have implemented in our code.
Consider the usual timelike vector n µ normal to the spatial hypersurfaces; the projector γ µν = g µν + n µ n ν defines the spatial 3-metric γ ij with the corresponding covariant derivative D i . From these, we obtain the following decomposition for the first derivatives of the scalar field: where L n denotes the Lie derivative along n µ . It follows that ∇ µ φ = Π µ − n µ Π and X = 1 2 (Π 2 − Π i Π i ). We also decompose the second derivatives of the scalar field, defining the auxiliary variables: and hence τ := τ i i = KΠ + D i Π i . Therefore, we get with n µ τ µ = 0 and n µ τ µν = 0. In terms of the usual conformal spatial metricγ ij := χγ ij (with det(γ ij ) = 1) and its associated covariant derivativeD i , we define the conformal variables for the scalar field as, Note that the indices ofτ ij are raised with the conformal metricγ ij so thatτ :=τ i i = τ , and similarly for all other conformal variables. For example,Π i = 1 χ Π i , which implies X = 1 2 (Π 2 − χΠ iΠ i ). With these definitions in place, the 3+1 conformal decomposition of the scalar energy-momentum tensor is: Similarly, the scalar field evolution equation (2.3) in first order form is given by (A.1) and: (A.10) Note that one can obtain the standard 3+1 evolution equations without a conformal transformation by setting χ = 1 and dropping any '˜' superscripts.
Regarding gauge and numerical evolution parameters, we choose 1 + log slicing and hyperbolic gamma-driver condition with the standard parameters. We use CCZ4 paramaters {κ 1 = 0.1 α , κ 2 = 0, κ 3 = 1} and Kreiss-Oliger numerical dissipation with σ = 0.3. Typical simulations used a Courant factor of 0.2 (reduced for Keldysh-type-of transitions), a coarse grid resolution of ∆x = 1 and up to 7 additional refinement levels, and a box size of L = 96 with Sommerfeld boundary conditions. We use the gradients φ and χ as well as contours of χ to tag cells for regridding. Last but not least, we use the symmetry of the system to only simulate one octant of the full domain, which reduces the computational cost of the problem.

A.2 Effective metric
As discussed in Section 2.4, the quantities −α 2 h 00 and det (h µ ν ) are useful to monitor the hyperbolicity of the scalar equation of motion and determine whether its change of character is of the Tricomi or Keldysh type. Here we present −α 2 h 00 and h µ ν in terms of the 3+1 conformal variables, which is how we have calculated them in our code: From (A.11) one can readily compute det(h µ ν ). If necessary, the effective metric with both indices up can also be obtained by raising the lower index in (A.11) with the spacetime metric.

B Determinant of the effective metric
We can compute det (h µ ν ) in full generality using Cayley-Hamilton's theorem and Newton's identities. The general case, with both G 2 = 0 and G 3 = 0, is not particularly insightful and in practice it is preferable to directly compute the determinant of the metric with a lowered index numerically. For clarity, in this Appendix we provide the explicit expression for the determinant in the case G 2 = 0 and G 3 = g 3 X: C Dealing with strong field regime inside black holes As described in Section 1.1, inevitably the evolution will exit the regime of validity of Horndeski theories inside black holes. To deal with this situation, in practice we excise a portion of the interior of the black hole. In this appendix we provide the details of our implementation. Rather than performing proper excision, i.e., cutting out a region of the domain, we found that it was easier to modify the evolution equations inside the black hole. The result should be the same as information cannot escape from this region. Note that in certain Horndeski theories, depending on the sign of the couplings, the scalar field can propagate faster than light and consequently the associated scalar apparent horizon will be inside the black hole horizon [65]. Therefore, to avoid unphysical effects leaking out of the black hole, any modification of the equations of motion should be done in a region contained within all apparent horizons.
Since puncture gauge can handle singularities very well in GR, in practice we turn off all Horndeski terms in a certain region inside the black hole and evolve the standard GR equations there. To do so, we first define a smooth transition function, valued between 0 and 1, as the sigmoid-like function: wherex represents the transition point, and w represents the transition width, relative tō x, such that wx is the actual width of the transition. 13 The metric apparent horizon can be accurately tracked during simulations, but for a sense of what is "well within the black hole", contours of the conformal factor χ are, in puncture gauge, an excellent measure. For example, for a Schwarzschild black hole, after puncture gauge settles, the apparent horizon corresponds to a contour of χ around 0.25, reducing to lower values as spin increases along the Kerr family of solutions. For reasonable choices of Horndeski couplings in the regime of validity of the theory, the scalar apparent horizon is close to the metric horizon. Therefore, the region inside a certain sufficiently small contour of χ should be contained in all apparent horizons. Denoting by W the maximum of all the weak field conditions (2.17), we define the excision function e(χ, W ) as: e(χ, W ) = σ(χ;χ, −w χ ) σ(W ;W , w W ) , (C.2) where w χ and w W are two adjustable parameters. In our simulations we typically used χ = 0.1, w χ = 0.2,W = 1, w W = 0.1. This choice is robust, in the sense that changing these barely affects the evolution across resolutions as long asχ is well within the black hole, which is the case for this choice. It follows from the definition (C.2) that e → 1 when χ <χ and W >W , and e → 0 otherwise. We then modify the right hand side of the evolution equations, collectively denoted by RHS, as: with e given by (C.2). In practice, we are only modifying the equations of motion in a region where the weak field condition is large and where the theory should not be trusted anyway.

D Convergence
In this appendix we provide details of some of the convergence tests that we have carried out. As an illustrative example, we consider the weak coupling g 2 = 0.005 case presented in 3.1.1. To carry out the convergence tests, we used simulations with coarsest level resolutions ∆x = 1 (low resolution, LR), ∆x = 0.75 (medium resolution, M R) and ∆x = 0.6 (high resolution, HR) respectively, all with the same 7 additional levels of refinement. The results of the simulations for the 4 cases analysed are shown in Fig. 11. The bottom panel shows the error estimates |M R − LR| (solid green curve) and |HR − M R| (solid purple curve), and compares them to the expected errors for 2 nd (dashed blue) and 4 th (dashed red) order convergence. The latter were obtained from the |HR − M R| error using the continuum limit of the convergence factor: (∆x LR ) n −(∆x M R ) n (∆x M R ) n −(∆x HR ) n . We see that our numerical results are consistent with convergence order between 2 and 4. Notice that it appears that the evolution has not reached a stationary state, but this should not be a concern since the outcome in terms of well-posedness and possible pathologies has already been determined after collapse occured.
We also monitor the behaviour of the Hamiltonian and Momentum constraints for the simulation with g 2 = 0.005 presented in 3.1.1. We measure the L 2 norm of a quantity Q by the volume average: where V is the volume of the box except the region excised inside black holes (if there are any present). We normalise the constraints by the norm of the sum of the absolute value of each term in the constraints. We show in Fig. 12  Hamiltonian Constraint Momentum Constraint Figure 12: L 2 norm of constraints for the g 2 = 0.005 run with coarsest level resolution of ∆x = 1 and 7 additional levels of refinement. We normalise the constraints by the norm of the sum of the absolute value of each term that composes it.
to infinity; as the scalar field propagates towards the boundaries, it moves away from the center of the grid into coarser refinement levels, and thus resolution is lost. Secondly, as the scalar field disperses or is absorbed by the black hole, matter terms in the constraints become increasingly smaller and, as a consequence, the normalisation factors used also significantly decrease. Therefore, we can conclude that we have a good numerical control over our simulations.

E Other cases of interest
In Figs. 13, 14 and 15 of this Appendix we collect the results of some simulations that are relevant for the discussion in the main text. For small enough initial data (Case 1) the evolution is perfectly consistent, while it breaks down in a Tricomi-type of transition in Cases 2 and 3. For large enough initial data (Case 4), the pathologies that may develop during the evolution are hidden behind the black hole horizon. In this case, the weak field conditions are small on and outside the black hole horizon. Figure 14: det(h µ ν ) (top) and −α 2 h 00 (bottom) for G 3 = g 3 X with g 3 = 0.4. The corresponding values of the dimensionless coupling η 3 , from left to right, are: 1.6 × 10 −5 , 3.2 × 10 −5 , 7 × 10 −5 , 1 × 10 −4 . In all cases the evolution breaks down because −α 2 h 00 → 0, signalling a Keldysh-type-of transition. Figure 15: Characteristic speeds of the outgoing (top) and ingoing (bottom) scalar modes for G 3 = g 3 X with g 3 = 0.4. v − diverges at the transition, but v + remains finite.