Gravitational Collapse in Einstein dilaton Gauss-Bonnet Gravity

TL;DR

This paper numerically investigates gravitational collapse in Einstein dilaton Gauss-Bonnet gravity, revealing how modifications to general relativity influence black hole formation and scalar field behavior, especially in strong coupling regimes.

Contribution

It provides the first detailed numerical analysis of spherical collapse in EdGB gravity, highlighting the transition from hyperbolic to elliptic PDE behavior and the emergence of scalar hair near black holes.

Findings

Weak EdGB coupling solutions approach GR results.

Strong coupling leads to elliptic regions where PDE character changes.

Negative energy densities and violation of null convergence condition occur near elliptic regions.

Abstract

We present results from a numerical study of spherical gravitational collapse in shift symmetric Einstein dilaton Gauss-Bonnet (EdGB) gravity. This modified gravity theory has a single coupling parameter that when zero reduces to general relativity (GR) minimally coupled to a massless scalar field. We first show results from the weak EdGB coupling limit, where we obtain solutions that smoothly approach those of the Einstein-Klein-Gordon system of GR. Here, in the strong field case, though our code does not utilize horizon penetrating coordinates, we nevertheless find tentative evidence that approaching black hole formation the EdGB modifications cause the growth of scalar field "hair", consistent with known static black hole solutions in EdGB gravity. For the strong EdGB coupling regime, in a companion paper we first showed results that even in the weak field (i.e. far from black hole formation), the EdGB equations are of mixed type: evolution of the initially hyperbolic system of partial differential equations lead to formation of a region where their character changes to elliptic. Here, we present more details about this regime. In particular, we show that an effective energy density based on the Misner-Sharp mass is negative near these elliptic regions, and similarly the null convergence condition is violated then.

Equations

$$S=\frac{1}{2}\int d^{4}x\sqrt{-g}\left(R-(\nabla \varphi {)}^2+2\lambda \varphi \mathcal{G}\right),$$

$$\mathcal{G}\equiv \frac{1}{4}{\delta }_{\rho \sigma \gamma \delta }^{\mu \nu \alpha \beta }{R}^{\rho \sigma }{}_{\mu \nu }{R}^{\gamma \delta }{}_{\alpha \beta },$$

$$E_{\mu \nu }^{(g)}$$

$$E^{(\varphi )}$$

$$d{s}^2=-e^{2A(t,r)}d{t}^2+e^{2B(t,r)}d{r}^2+r^2\left(d{\vartheta }^2+{\sin }^2\vartheta d{\phi }^2\right).$$

$$Q\equiv$$

$$P\equiv$$

$$E^{(A)}\equiv$$

$$E^{(B)}\equiv$$

$$E^{(Q)}\equiv$$

$$E^{(P)}\equiv$$

$$\mathcal{B}\equiv (1-e^{-2B})/r^2,\mathcal{I}\equiv 1-8\lambda e^{-2B}Q/r,\mathrm{and}$$

$$\mathcal{K}\equiv e^B\frac{\frac{PQ}{2}+4\lambda \mathcal{B}\left(-P{\partial }_{r}B+{\partial }_{r}P\right)}{e^{2B}+4\lambda \left(-3+e^{2B}\right)\frac{Q}{r}}.$$

$$\displaystyle\partial_{r}A(t,r)\Big{|}_{r=0}$$

$$\displaystyle B(t,r)\Big{|}_{r=0}=0,\ \ \ \partial_{r}B\Big{|}_{r=0}$$

$$\displaystyle Q(t,r)\Big{|}_{r=0}$$

$$\displaystyle\partial_{r}P(t,r)\Big{|}_{r=0}$$

$$\displaystyle\partial_{t}Q+\frac{1}{r}\partial_{r}\left(rQ\right)\Big{|}_{r=r_{\max}}$$

$$\displaystyle\partial_{t}P+\frac{1}{r}\partial_{r}\left(rP\right)\Big{|}_{r=r_{\max}}$$

$$\phi(t,r)\Big{|}_{t=0}=a_{0}\left(\frac{r}{w_{0}}\right)^{2}\mathrm{exp}\left(-\left(\frac{r-r_{0}}{w_{0}}\right)^{2}\right),$$

$$P(t,r)\Big{|}_{t=0}=-\frac{1}{r}\phi(t,r)-Q(t,r)\Big{|}_{t=0}.$$

$$E^I\left(v^J,{\partial }_a{v}^K\right)=0,$$

$${\mathfrak{p}}_J^I\left({\xi }_a\right)\equiv \frac{\delta E^I}{\delta ({\partial }_a{v}^J)}{\xi }_a,$$

$$\det \left({\mathfrak{p}}_J^I\left({\partial }_a\Sigma \right)\right)=0,$$

$${\tilde{E}}^I\left(v^J,{\partial }_a{v}^K\right)=0,$$

$$\mathfrak{p}(\xi )=\tilde{\mathfrak{a}}{\xi }_t+\tilde{\mathfrak{b}}{\xi }_r,$$

$$\tilde{\mathfrak{a}}\equiv$$

$$\tilde{\mathfrak{b}}\equiv$$

$$c_{\pm }=\frac{1}{2}\left(\mathrm{Tr}\left(\tilde{\mathfrak{c}}\right)\pm \sqrt{\mathrm{Tr}{\left(\tilde{\mathfrak{c}}\right)}^2-4\mathrm{Det}\left(\tilde{\mathfrak{c}}\right)}\right),$$

$$\tilde{\mathfrak{c}}\equiv {\tilde{\mathfrak{a}}}^{-1}\cdot \tilde{\mathfrak{b}},$$

$$\mathcal{D}\equiv \mathrm{Tr}{\left(\tilde{\mathfrak{c}}\right)}^2-4\mathrm{Det}\left(\tilde{\mathfrak{c}}\right)$$

$$\displaystyle\tilde{\mathfrak{a}}\big{|}_{\lambda=0}=$$

$$\displaystyle\tilde{\mathfrak{b}}\big{|}_{\lambda=0}=$$

$$\mathrm{Im}(\tilde{c})=\mathrm{Im}(c)\frac{\det [{\Lambda }^a{}_{\tilde{a}}]}{Z{Z}^{\ast }},$$

$$m_{MS}(t,r)\equiv \frac{r}{2}\left(1-(\nabla r{)}^2\right)=(r/2)\left(1-e^{-2B(t,r)}\right),$$

$${\mathcal{T}}_{\mu \nu }\equiv -2\lambda {\delta }_{\alpha \beta \rho \sigma }^{\gamma \delta \kappa \lambda }{R}^{\rho \sigma }{}_{\kappa \lambda }\left({\nabla }^{\alpha }{\nabla }_{\gamma }\varphi \right){\delta }_{(\mu }^{\beta }{g}_{\nu )\delta }+{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\left(\nabla \varphi \right)}^2.$$

Full text

• February 2019

Keywords: modified gravity, numerical relativity

1 Introduction

While General Relativity (GR) has passed all experimental and observational tests so far (caveats with dark energy and dark matter aside), there are well known reasons to suspect that GR is not a complete theory of gravity. One reason is that at the level of the classical equations of motion, black hole (BH) and most cosmological solutions are geodesically incomplete [1], with the expectation that these spacetimes also generically contain curvature singularities. Another is that as matter is quantum in nature, the Einstein field equations relating a classical description of geometry to a classical stress energy tensor of matter can only be an approximate theory; though the predominant opinion today seems to be that the resolution of this issue is that geometry is also fundamentally quantum (as opposed to, for example, a completely novel theoretical construct that reduces to the modern, tested physical theories in appropriate limits), there is no evidence for this at present. Though regardless of theoretical reasons to think GR is incomplete, from a purely empirical point of view we have only recently entered the era where we can begin to verify the dynamical strong field predictions of GR, through gravitational wave (GW) observation of compact object mergers. Though initial tests are consistent with the GR description of these events [2], we are still in the early days of GW astronomy, and the data cannot yet provide high precision tests of this regime.

One problem with achieving the tightest possible constraints on deviations from GR in the strong field (or discovering them), is at present we have no interesting, viable alternatives to GR that can give quantitative predictions to the analogue of the merger regime of BH inspiral in GR. This is where the predominant share of SNR (signal-to-noise ratio) is coming from with current detections (in particular GW150914), and where one might expect to see the first hints of corrections to GR. By interesting we mean theories that when restricted to regimes that are consistent with existing non-GW tests nevertheless still offer significant differences for BH mergers; by viable we mean theories that possess a well-posed initial value problem (IVP) that can be solved to make predictions of mergers to confront with data. For example, a class of viable but uninteresting theories in this regard (and we emphasize we certainly do not mean “uninteresting” for any other reason) are the typical scalar tensor theories, such as Brans-Dicke, as they have the same vacuum sector of solutions, hence BH mergers, as GR. Another example that is viable but likely uninteresting is Einstein-Maxwell-Dilaton gravity; here the problem is to obtain mergers with noticeable deviations in the GW emission requires what is expected to be astrophysically unrealistic amounts of electric charge for the BHs [3].

Over the past several years two likely interesting modified gravity theories have attracted the attention of researchers attempting to study the full non-linear BH merger problem [4, 5, 6, 7]: dynamical Chern-Simons (dCS) gravity (see e.g. [8]), and Einstein-dilation-Gauss-Bonnet (EdGB) gravity(see e.g. [9, 10, 11]), the latter being the focus of this paper. One the main motivations of these studies, and ours here, is to understand how in principle strong field merger dynamics could differ from GR (as opposed to any observational or theoretical impetus arguing for such modifications on the scale of astrophysical BHs). In terms of testing GR though BH mergers, what is interesting about the particular variant of EdGB gravity we consider here is it does not admit the Schwarzschild or Kerr BH solutions of GR. Instead, the analogue BH solutions only exist above a minimum length scale related to the coupling constant $\lambda$ in the theory, and feature scalar “hair” [9, 12, 13]. Moreover, for values of $\lambda$ that would produce significant changes in stellar mass BHs, the corresponding effect on material compact objects such as neutron stars is insignificant [11], implying this theory could be consistent with current GR tests, yet give different results for stellar mass BH mergers.

One problem with EdGB gravity relates to whether it is viable in the above sense of the word. This paper is a follow up to a first study [14] of gravitational collapse in EdGB gravity in spherical symmetry to begin to address this issue from the level of fully non-linear dynamical solutions. Earlier work on the well-posedness of EdGB gravity [15, 16] considered the linearized equations in the small coupling parameter limit, and found that these equations are at best weakly hyperbolic about generic backgrounds, within a class of “generalized harmonic” gauges. This is certainly a sign for concern, however demanding that a theory be well-posed in all possible situations might be unnecessarily restrictive if problems do not occur in scenarios of interest, here in particular for binary BH mergers. Considering the recent results of [17], we also mention there may exist other gauges for which EdGB may have well posed initial value problem for generic small field initial data.

Given how challenging solving for BH merger spacetimes is in GR alone, it makes sense to tackle this problem in EdGB gravity beginning with simpler scenarios that capture some aspects of the final problem, uncover any issues that might arise, and if there are no show-stoppers, move forward. One approach along these lines follows an effective field theory interpretation of EdGB gravity, beginning with the GR solution, then examining perturbative corrections. The advantage to this approach is one can begin with fully non-linear GR BH merger solutions. The first step here is the so-called decoupling limit (see D), where the EdGB scalar is not yet allowed to back-react on the geometry; this has successfully been carried out in  [7] (a similar approach has been taken in dCS gravity in [4], and even recently extended to first order metric perturbations [5]).

Another approach, that we follow here, is to begin with the fully non-linear, non-perturbative EdGB equations, but in a symmetry reduced setting. The benefit of this is we can immediately begin looking for non-perturbative deviations from the predictions of GR. The natural, simplest symmetry to consider for our purposes is spherically symmetry, as this allows us to study black hole formation in asymptotically flat, 4-dimensional spacetimes. One key result of our initial study, described in [14], is within this symmetry class we do identify a regime of EdGB gravity that is “pathological” from the perspective of having a well-posed IVP : specifically, in the strong coupling regime, we find scenarios where evolution of initial data leads to the EdGB dilaton equation changing character from hyperbolic to elliptic within a region of the spacetime (or said another way, this equation is then actually of mixed type). However, as discussed more in [14], given how this phenomenon scales with the magnitude of the coupling parameter, there are regimes of EdGB gravity that may yet offer a viable, interesting modified gravity scenario for application to GW astronomy (the main limitation of this first study, aside from symmetry considerations, is since we do not use horizon penetrating coordinates, we cannot address the long time, non-linear stability of BH’s regardless of the magnitude of the EdGB coupling parameter).

We should also mention that variants of EdGB gravity have been extensively studied in the cosmological context. There is a vast literature on this, and since this topic is outside of our scope, we will not attempt to cite the relevant papers, but instead refer the reader to the recent review articles [18, 19] for further reference. In cosmology the equations are usually analyzed in the form of Horndeski theories (roughly Einstein gravity coupled to a scalar field with all possible non-standard kinetic terms that yield second order equations of motion), with EdGB gravity being a sub-class of the most general Horndeski theory. The motivations there are more often driven by the need to explain dark energy, or to come up with models of the early universe, including inflation or bouncing models. What is particularly interesting with regard to the last-named problem is some Horndeski theories can violate the null convergence condition (NCC, or equivalently the null energy condition (NEC)), leading to concrete realizations of non-singular classical bounces. Here we also find that EdGB gravity can violate the NCC in spherical gravitational collapse.

Also as we find here, for some Horndeski theories in cosmological settings the equations appear to have regimes where mixed type character is present. However, as far as we are aware, all these analysis have been carried out at the level of linear perturbations about a cosmological background solution, and often the corresponding elliptic regions are ascribed to be subject to a gradient, or Laplace instability. This is a misnomer in a sense, as the “instability” is an artifact of analyzing an elliptic region of a partial differential equation (PDE) assuming it were hyperbolic (as opposed to a physical instability in a system described by hyperbolic PDEs where, for example, exponentially growing modes can be excited). As in our case, this means those Horndeski scenarios do not admit a well-posed hyperbolic IVP, but does not imply that a sensible interpretation as a mixed type problem is impossible.

Note that when many of these modified gravity theories are directly applied to address questions on a cosmological scale, things can “break” on smaller scales such as the solar system, compact objects etc., and vice versa (hence the need to invent screening or “chameleon” mechanisms). This is certainly the case with EdGB gravity, and so for it to have a chance of still being interesting and viable for BH mergers (and barring invention of a screening mechanism) we must assume the EdGB scalar field is irrelevant on cosmological scales. If it did have a large cosmological value (where large means its contribution to the normalized energy density of the universe is $\Omega_{EdGB}\sim O(1)$), to avoid formation of elliptic regions on smaller scales and subsequent breakdown of the IVP would require a coupling parameter so small it would be completely uninteresting for GW tests of GR [14]. Moreover, even if one assumed such mixed type character was benign and would not lead to unexplained phenomena on smaller scales, measurement of the speed of GWs implied by the binary neutron star merger GW170817 together with counterpart electromagnetic signals [20] rules out large couplings if there are cosmologically relevant scalars [21].

Mixed type behavior and elliptic region formation has been observed and discussed in the context of collapse simulations of other modified gravity theories [22, 23, 24, 25]. Reference [26] discusses the appearance of mixed type PDEs in loop quantum gravity models of the early universe, and in the Hartle-Hawking no boundary proposal. Interestingly, there the signature change is interpreted as a property of the model, rather than signalling a pathology, and proposals are made to solve the corresponding mixed type equations. Reference [27] provides a similar discussion of mixed type problems, and proposes methods to solve this class of PDE in the context of numerical relativity. A more complete account of the appearance of mixed type PDEs in physics and applied mathematics, and some of the attempts to systematically understand them, may be found in [28].

1.1

Layout of the remainder of the paper, and conventions

An outline for the remainder of the paper is as follows. In Section 2, we describe the particular variant of EdGB gravity we study, write out the form of the equations within the spherically symmetric ansatz, discuss relevant initial, boundary and regularity conditions, and briefly mention the numerical methods we use to solve these equations (more details on the numerics are given in E). We then describe the main analysis tools we employ to understand properties of the solutions : the characteristics of the theory in Section 3, a quasi-local mass measure in Section 4, and the NCC in Section 5. Following that we give results from numerical solutions of several representative members of our initial data family : Section 6 contains a case in the weak field, weak coupling regime, Section 7 contains a case from the strong field, weak coupling regime, and Section 8 discusses several cases from the (moderately) weak field, strong coupling regime (this was the regime initially presented in [14], where we also give results scaling to the truly weak field, strong coupling limit). We discuss potential future directions in the conclusion; in particular to study long term BH stability in the weak coupling regime, or early time behavior in the strong field, strong coupling regime will require the use of horizon penetrating coordinates.

We give some details of the derivation of the EdGB equations in A, the specific form of the components of the tensor equations of motion within our spherically symmetric ansatz in B, a second method to compute the characteristics in C (largely equivalent to the method described in Section 3) a derivation of the ‘decoupled’ EdGB scalar profile about a Schwarzschild black hole background in D, and a description of all the numerical methods we employed to solve the EdGB PDEs in E.

We used geometrized units where $G=c=1$, and use MTW [29] sign conventions for the metric tensor, etc.

2 Basic equations

2.1

Shift-symmetric dilaton Gauss-Bonnet gravity

The action for the EdGB gravity theory we consider is

[TABLE]

where $R$ is the Ricci scalar, $g$ is the determinant of the metric tensor $g_{\mu\nu}$, $\phi$ is the dilaton field, and $\mathcal{G}$ is the Gauss-Bonnet scalar that can be written in terms of the Riemann tensor $R_{\rho\sigma\mu\nu}$ as

[TABLE]

with $\delta^{\mu\nu\alpha\beta}_{\rho\sigma\gamma\delta}$ the generalized Kronecker delta. In our units the Gauss-Bonnet coupling constant $\lambda$ has dimension $[L]^{2}$. Varying (1) with respect to $g^{\mu\nu}$ and $\phi$ (see A) we obtain

[TABLE]

There are several theories which go under the name of EdGB gravity, each of which differs by the functional form of the coupling between the dilaton and the Gauss-Bonnet scalar, or the presence of a potential $V(\phi)$ for the dilaton in the action (we consider $V(\phi)=0$). For example, one variant of dilaton Gauss-Bonnet gravity appears as a leading order term in the low-energy effective action to certain string theories (e.g. [30, 31]): there the coupling goes as $\alpha e^{-\gamma\phi}$, where $\alpha$ and $\gamma$ are constants that are set by the string theory in question. The theory $(\ref{eq:EdGBAction})$ we consider is equivalent to this to leading order in the dilaton coupling (with $\alpha\gamma\propto\lambda$, and recalling that any constant times $\mathcal{G}$ in the action in 4-dimensional spacetime can be replaced by a boundary term that does not affect the equations of motion), and goes by several names: ‘shift symmetric Truncated Einstein dilaton Gauss-Bonnet gravity’ (e.g. [10, 11]) or ‘shift symmetric Einstein dilaton Gauss-Bonnet gravity’ (e.g. [7]); for brevity we will refer to it simply as EdGB gravity.

EdGB gravity is invariant under $\phi\to-\phi$ and $\lambda\to-\lambda$. In this work we present simulations with initial data for $\phi$ that is everywhere positive, and consider $\lambda$ with positive and negative values.

2.2 Spherical symmetry and coordinate ansatz

We work in spherical symmetry in polar coordinates, and use the following ansatz for the line element:

[TABLE]

To reduce the dilaton equation of motion (3b) to a set of first order PDEs, we define the variables

[TABLE]

The equations for the metric (3a) and scalar (3b) retain a similar structure to the Einstein-Klein-Gordon equations ($\lambda\rightarrow 0$) in spherical symmetry; namely, there are no gravitational degrees of freedom, and all dynamics are driven by the scalar field. Hence, we can consider a fully constrained evolution scheme to solve for the coupled system of equations. In such a scheme, the metric variables are solved for using what are essentially the elliptic Hamiltonian and momentum constraint equations in GR, and the scalar equation is treated as hyperbolic. The Klein-Gordon equation is always hyperbolic in GR (away from coordinate or geometric singularities), though of course one of the main results of our study (first described in [14]) is this property fails to generically hold in the strong-coupling limit of EdGB gravity. Also akin to the Einstein equations, the EdGB system of PDEs permit a partially constrained evolution scheme, where a hyperbolic evolution equation is used for the metric variable $B$ instead of a constraint equation; as mentioned in the appendix, we have also implemented such a scheme, and verified we obtain consistent results to within truncation error, though for brevity all the characteristic analysis and simulation results presented here pertain to the fully constrained scheme.

Following the fully constrained evolution strategy, we take particular algebraic combinations of the nontrivial components of the EdGB equations of motion (see B) to give the following closed system of PDEs that we solve using numerical methods:

[TABLE]

where

[TABLE]

In particular (3efc) and (3efd) define the PDE evolution equations for $Q$ and $P$ respectively; (3efa) and (3efb) contain no time derivatives and are the ODE (ordinary differential equation) constraint equations for $A$ and $B$ respectively.

2.3 Boundary and regularity conditions

We discretize the above equations over a domain $r\in[0..r_{\rm max}]$. At the origin $r=0$ we require regularity of the fields, leading to

[TABLE]

Equations (3efa,3efb) are first order ODEs for $A$ and $B$, so strictly speaking we can only impose one boundary condition at one of the boundaries for each. In practice we integrate from $r=0$ to $r_{\rm max}$, setting $A(t,r=0)=0$ and $B(r,r=0)=0$; with the scalar field variables $P$ and $Q$ appropriately regular as above, the structure of the field equations guarantees that $A$ and $B$ also satisfy the above regularity conditions. Our coordinate system (3d) has residual gauge freedom in that we can rescale $t$ by an arbitrary function of itself, and we use this to rescale $A(t,r)$ after each ODE integration step so that $A(t,r)|_{r=r_{max}}=0$. In that way our time coordinate $t$ measures proper time of static observers at the outer boundary.

For $Q$ and $P$ at the outer boundary we impose the following approximate outgoing radiation boundary conditions:

[TABLE]

2.4

Initial data

For initial data, we are free to choose $P(t=0,r)$ and $Q(t=0,r)$ (subject to the regularity conditions described in the previous subsection). For the simulation results presented here, we begin with the following family of initial data for $\phi(t=0,r$):

[TABLE]

where $a_{0}$, $w_{0}$, and $r_{0}$ are constants. This then gives $Q(t,r)|_{t=0}=\partial_{r}\phi|_{t=0}$, and we choose $P$ so that the scalar pulse is initially approximately ingoing:

[TABLE]

Because of spherical symmetry and our constrained evolution scheme, the only “free” data for the metric variables $A$ and $B$ is the overall scale of $A$, which as discussed in the previous subsection we set so that $t$ measures proper time for static observers at the outer boundary of our domain.

2.5 Numerical solution methods in brief

We numerically solve Equations (3efa)-(3efd) using (overall) second order accurate finite difference techniques. We implemented several different solution methods as detailed in  E. In particular, to gain confidence that the late time convergence problems for strong coupling cases are due to the character of the $(P,Q)$ subsystem changing from hyperbolic to elliptic in a certain region, and not due to an unstable numerical evolution method, we explored two completely different discretization and evolution schemes for $(P,Q)$ : (i) a Crank-Nicolson method (with both a Newton-Gauss-Seidel iterative solver, and Newton iteration together with a fully implicit matrix inversion for each linear step of the Newton iteration), and (ii) a fourth order in time Runge-Kutta (method of lines) solver. We integrate the constraint ODEs (3efa) and (3efb) with a trapezoidal method, which is ‘A stable’ [32] (we also experimented with a few of variants to deal with the non-linearities in the equation for $B$, as outlined in the appendix). All schemes give solutions consistent with each other to within truncation error. Here then, all results we show were obtained using the iterative Crank-Nicholson scheme. We further experimented with a range of Courant-Friedrichs-Lewy (CFL) factors $0.01-0.5$, confirming results do not qualitatively depend on this; the simulations discussed here used a CFL factor of $0.25$.

3 Hyperbolicity analysis

We briefly summarize the theory of characteristics; standard references include [33, 34, 35]. Consider a system of first order PDEs 111Through field redefinitions essentially any system of PDEs may be written in this form.

[TABLE]

where $I,J,K$ index the $N$ equations of motion and dynamical fields $v^{J}$, and $a$ indexes the $n$ coordinates $\{x^{a}\}$ of the underlying (spacetime) manifold $M$ (and here because of our restriction to spherical symmetry $a$ only runs over the $(t,r)$ coordinates). The principal symbol is defined to be

[TABLE]

where $\xi_{a}$ is an $n$ dimensional covector. A characteristic surface $\Sigma\subset M$ by definition satisfies the characteristic equation

[TABLE]

and (3efhin) is the eikonal equation for the characteristic surface.

For a physical interpretation of characteristics, we consider a system of $N$ first order PDEs for $N$ fields (i.e. of the form (3efhil)) that is totally hyperbolic: i.e it has $N$ real (possibly degenerate) characteristic surfaces222For a system that is not totally hyperbolic we could instead consider a totally hyperbolic subsystem; see e.g. Section 3.1. Our treatment of characteristics roughly follows that of [36].. Consider the solution to small amplitude high frequency wave solutions: $v_{0}^{I}e^{ik_{a}x^{a}/\epsilon}$, with $0<\epsilon\ll 1$. Solutions of this form to leading order in $\epsilon$ satisfy $\mathfrak{p}^{I}_{J}(k_{a})v_{0}^{J}=0$. Nontrivial solutions to this equation exist if and only if $\mathrm{det}\left(\mathfrak{p}^{I}_{J}(k_{a})\right)=0$; i.e. if and only if the wave vector satisfies (3efhin). Thus the wavefronts of high frequency wave solutions propagate on the characteristic surfaces. The characteristic surfaces locally delimit the causal region of influence for hyperbolic PDE (e.g. [37]).

In local coordinates, letting $t$ index the timelike coordinate and $i$ index the spacelike coordinates of the background geometry, the speed of these perturbations for the $n^{th}$ characteristic is given by $c^{(n)}=(v^{(n)})^{i}/v^{t}$, where $(v^{(n)})^{\mu}$ is a vector parallel to the $n^{th}$ characteristic surface. We may relate $c^{(n)}$ to the characteristic covector by noting that since locally the $n^{th}$ characteristic covector is equal to the gradient of the $n^{th}$ characteristic surface, $\xi^{(n)}_{\mu}=\partial_{\mu}\Sigma^{(n)}$, then $(v^{(n)})^{\mu}\xi^{(n)}_{\mu}=0$, from which we find $c^{(n)}=-\xi^{(n)}_{t}/\xi^{(n)}_{i}$. For a simple example of this procedure, consider the $1+1$ dimensional scalar transport equation $\partial_{t}\psi+v\partial_{x}\psi=0$. The symbol is $\mathfrak{p}=\xi_{t}+v\xi_{x}$, the characteristic equation is $\xi_{t}+v\xi_{x}=0$, and the speed of propagation along the characteristic is $-\xi_{t}/\xi_{x}=v$.

We compute the characteristic vectors and speeds for the system of PDEs (3efa), (3efb), (3efc), and (3efd) in two different ways. In the first, discussed below, we only consider the $P,Q$ evolution subsystem, eliminating all $A$ and $B$ gradients from these equations using the constraints. In the second (discussed in C), which is more for a consistency check than anything else, we apply the characteristic analysis verbatim to the full system of equations, obtaining the same results for $P,Q$ as the first, and confirming that $A$ and $B$ are elliptic.

3.1

Characteristics calculation

Eliminating $\partial_{r}A$ and $\partial_{r}B$ from Equations (3efc) and (3efd) using Equations (3efa) and (3efb), we write the scalar field system in the same form as before,

[TABLE]

but now $I,J,K$ only index the fields $Q$ and $P$. The principal symbol then reads

[TABLE]

where

[TABLE]

Solving the characteristic equation for the characteristic speeds $c\equiv-\xi_{t}/\xi_{r}$, we obtain

[TABLE]

where

[TABLE]

From standard PDE theory, the sign of the discriminant

[TABLE]

of (3efhiqr) at any point of the spacetime determines the character of the PDE there: when $\mathcal{D}>0$ it is hyperbolic, when $\mathcal{D}=0$ it is parabolic, and when $\mathcal{D}<0$ it is elliptic.

In general when $\lambda\neq 0$, $\mathrm{Tr}\left(\mathfrak{c}\right)\neq 0$, so that $c_{+}\neq-c_{-}$. In GR ($\lambda=0$), we have

[TABLE]

so that $\mathrm{Tr}\left(\mathfrak{c}\right)=0$, and the characteristic speeds are c_{\pm}\big{|}_{\lambda=0}=\pm e^{A-B}. The general expressions for the components of the matrices $\tilde{\mathfrak{a}}$, $\tilde{\mathfrak{b}}$, and $\tilde{\mathfrak{c}}$ can be obtained through straightforward algebraic manipulation of Equations (3efa)-(3efd); the resultant expressions are long and not particularly insightful, so we do not write out their full forms here.

3.2

Invariance of the characteristics under coordinate transformations

As mentioned, one of the main results of our study is that the EdGB equations in spherical symmetry can be of mixed elliptic/hyperbolic type in certain scenarios. Specifically then, evolution, beginning with initial data where the scalar equation is everywhere hyperbolic, leads to formation of a region where the characteristic structure switches to elliptic (separated by a parabolic so-called sonic line, though generically is a co-dimension one surface and not a “line”). The elliptic region is particularly problematic for the validity of EdGB gravity as a classically well-posed, predictive modified theory of gravity (see the discussion in [14]), and so we would like to be certain that our identification of the elliptic region is not somehow a coordinate artifact. It is well-known that the characteristic structure of a PDE is invariant under so-called point transformations (essentially coordinate transformations treating all the dependent variables as scalars), though it is unclear that this must hold when solving the PDEs of EdGB gravity as a Cauchy IVP problem in an arbitrary gauge. The problem is that the structure and even rank of the principal symbol is unknown until the gauge equations have been chosen (in the ADM [Arnowitt-Deser-Misner] language, that would be the equations governing the lapse $\alpha$ and shift vector $\beta^{i}$: we have effectively chosen the equation for the shift to be the algebraic condition $\beta^{i}=0$, and our choice of polar-areal coordinates in spherically symmetry fixes $\alpha$ to within an overall scale).

However, at least we can show that the characteristic structure is invariant under coordinate transformations in the following sense (this is effectively the point-transformation calculation). In our evolution scheme we compute the characteristic surfaces $\Sigma$ as outlined above, and find the corresponding co-vectors $\xi_{a}=\partial_{a}\Sigma$. If at a given point $c\equiv-\xi_{t}/\xi_{r}$ is purely real, we know the PDE is hyperbolic at that point, and information will propagate along the characteristic surface. This could be superluminal, luminal, or sub-luminal relative to the metric light-cone depending on whether $\xi^{a}\xi_{a}$ is negative, zero, or positive respectively, but the scalar equation is still hyperbolic and will have its own causal-cone of influence. If $c$ has an imaginary component this is no longer true, and the PDE is elliptic. The question then is whether this property of the characteristic surface is invariant under coordinate transformations, and the answer is yes. For consider a general coordinate transformation respecting the spherical symmetry of the spacetime : let $x^{a}=x^{a}(\tilde{x}^{\tilde{a}})$ where $x^{a}$ and $\tilde{x}^{\tilde{a}}$ denote the $(t,r)$ and $(\tilde{t},\tilde{r})$ coordinates respectively, and the Jacobian of the transformation is $\Lambda^{a}{}_{\tilde{a}}\equiv\partial x^{a}/\partial\tilde{x}^{\tilde{a}}$ (with all coordinates and metrics real). Then $\xi_{\tilde{a}}=\xi_{a}\Lambda^{a}{}_{\tilde{a}}$, the new coordinate speed is $\tilde{c}\equiv-\xi_{\tilde{t}}/\xi_{\tilde{r}}$, and it is straight-forward to calculate that

[TABLE]

where $Z\equiv-c\Lambda^{t}{}_{\tilde{r}}+\Lambda^{r}{}_{\tilde{r}}$, and $Z^{\ast}$ its complex conjugate. In other words, as long as the transformation is non-singular an imaginary piece to $c$ in one coordinate system implies one in all.

3.3

Horizons

As mentioned above, when hyperbolic, the causal cones of the scalar degree of freedom in EdGB gravity ($\lambda\neq 0$) do not generally coincide with those of the spacetime. The latter would govern the speed of propagation of fields minimally coupled to the metric, such as a massless scalar field or a Maxwell field, and gravitational waves (which are not present in spherical symmetry). Regarding metric horizons, our coordinate system does not allow evolution through formation of a black hole, as the geometric light speeds are c_{g\pm}\big{|}=\pm e^{A-B}; i.e. the metric is necessarily singular along horizons. In strong-field evolutions we estimate that gravitational collapse occurs when $c_{g+}$ starts to evolve to zero at a finite radius; evolution beyond horizon formation will require the use of horizon penetrating coordinates, which we leave to a future study.

4 Quasi-local mass

In spherical symmetry in GR, a standard definition of quasi-local mass is the Misner-Sharp mass [38, 39] (sometimes also referred to as the Hawking-Israel or Hernandez-Misner mass, e.g. [40])

[TABLE]

where $r$ is the areal radius, and the last term on the right is the specific form it takes in our coordinate system. The Misner-Sharp mass satisfies several useful criteria for a quasi-local mass (e.g. [41]). For example, in asymptotically flat spacetimes it reduces to the ADM mass at spatial infinity and the Bondi-Sachs mass at future null infinity. A further useful property of the Misner-Sharp mass is in spherical symmetry one can relate it to the charge associated with the Kodama current [42], which satisfies a conservation law purely from properties of the Einstein tensor in spherical symmetry: one does need to a priori connect the Einstein tensor to the matter stress energy tensor to prove this (see e.g. [40]). Therefore it is reasonable to use the Misner-Sharp mass in spherically symmetric EdGB gravity as a measure of geometric mass. Then (akin to GR), if desired we can use the EdGB equations of motion to relate it to an integral of an effective matter energy density. Specifically, we write the EdGB equations (3a) as $G_{\mu\nu}=\mathcal{T}_{\mu\nu}$, with

[TABLE]

Then, replacing $G_{\mu\nu}$ with $\mathcal{T}_{\mu\nu}$ in the Kodama current (see e.g.[40]), a short calculation gives the following integral for the Misner-Sharp mass in our coordinate system (assuming regularity at $r=0$):

[TABLE]

The effective stress tensor $\mathcal{T}_{\mu\nu}$ does not always satisfy the usual energy conditions, hence $m_{MS}(t,r)$ is not necessarily a monotonically increasing function of $r$, as it is in GR coupled to “ordinary” matter. We will show some examples below illustrating the non-monotonicity of $m_{MS}$.

5

Null convergence condition

The null convergence condition (NCC) is

[TABLE]

for all null vectors $l^{\alpha}$. The NCC plays a role in, for example the classical black hole and cosmological singularity theorems [1], the laws of black hole mechanics and dynamical horizons [43, 44, 45], and in the “topological censorship” theorems [46, 47]. It is often stated that these theorems and properties rely on the null energy condition (NEC), $T_{\mu\nu}l^{\mu}l^{\nu}\geq 0$, however that comes from replacing the Ricci tensor in the above with an equivalent function of the stress energy tensor using the Einstein equations. We could like wise recast our analysis in terms of a NEC using the effective stress energy tensor introduced in the previous section, though we prefer the geometric interpretation of the NCC.

Related to the fact that the Misner-Sharp mass does not always monotonically increase with radius as discussed in the previous section, EdGB gravity does not generically satisfy the NCC (coupling to matter than satisfies the NEC). This can be seen by contracting (3a) with $l^{\mu}l^{\nu}$ to compute the explicit form of $R_{\mu\nu}l^{\mu}l^{\nu}$:

[TABLE]

Here, the only term that is manifestly positive definite is the kinetic term of the scalar in the small coupling ($\lambda\rightarrow 0$) limit. We will show examples below of scenarios where the NCC is violated during evolution in EdGB gravity (and the regions where it does roughly coincide with negative effective energy density in the Misner-Sharp mass, and is present where the equations become elliptic); specifically, we numerically evaluate $R_{\mu\nu}l^{\mu}l^{\nu}$ for outgoing null vectors $l^{\mu}\equiv(e^{-A},e^{-B},0,0)$.

6

Numerical results: weak field, weak coupling

Figure 1 shows results from the evolution of a representative member of the initial data family (3efhij) corresponding to a weak field and weak coupling case: the compaction of the scalar field $m_{ADM}/w_{0}\ll 1$, and $\lambda/w_{0}^{2}\ll 1$ respectively. The scalar pulse bounces smoothly off the origin and disperses to infinity. Throughout the evolution, the characteristics remain real and close to the geometric null characteristics. To within truncation error, the Misner-Sharp mass is monotonically increasing in $r$ and the NCC is preserved. Qualitatively the evolution matches that of Einstein massless scalar field evolution ($\lambda=0$).

7

Numerical results: strong field, weak coupling

We ran simulations with initial data that in GR would form a geometric horizon at $r\sim r_{gh}$, but weakly perturbed by EdGB modifications with $\lambda/r_{gh}^{2}\ll 1$. As our coordinates are not horizon penetrating, we cannot evolve the spacetime to or beyond (geometric) apparent horizon formation, which in our coordinates is signaled by $e^{A-B}\to 0$ (see Section 3.3). At all the finite resolutions we have used there is some small value of $e^{A-B}$ below which we loose convergence (and both metric fields diverge here, with $A\to-\infty$ and $B\to\infty$). Though before this, we do observe that the EdGB scalar field begins to grow near the nascent horizon. The growth is in qualitative agreement with the conclusions of [9, 13, 6, 48, 7], where the value of the field at the horizon is expected to asymptote to a unique value $\phi_{gh}\propto\lambda/r_{gh}^{2}$ depending upon the mass of the black hole (though of course questions about the ultimate stability of such “hairy” black holes, and under what conditions no elliptic regions form outside the horizon will require numerical solutions using horizon penetrating coordinates).

Figure 2 shows results of an example from such a strong field, weak coupling case. The scalar pulse approaches the origin, then “freezes” interior to what will be the eventual horizon (since the lapse function $\alpha=e^{A}\to 0$ there). Outside, the scalar field begins to grow, and for a while we can follow its evolution before the code fails. In Figure 3 we show a zoom-in of the scalar field at such a late time, together with the expected solution in the decoupling limit for a regular EdGB scalar field on a static black hole background (see D). That this solution qualitatively matches well with the full nonlinear evolution is another indication that this is in the weak coupling regime of EdGB gravity.

8

Numerical results: strong coupling, weak field

We first presented results from the weak field, strong coupling regime in [14]; here we describe two additional examples, and give more details. Specifically, we consider initial data (3efhij) with $a_{0}=0.02$, $r_{0}=20$, $w_{0}=8$, and $\lambda=\pm 50$; $m\sim 0.93$ for both cases (so this is fairly compact initial data, but is “weak” in the sense that we are still a factor of a few in mass away from initial data that would form a black hole; in [14] further data was given showing scaling to the truly weak field (low compaction), strong coupling regime).

8.1

Characteristics and formation of elliptic regions

For both cases (i.e. independent of the sign of $\lambda$), the solutions develop an elliptic region—see Figures 4 and 5. Interestingly, even though the sign of the Gauss-Bonnet coupling $\lambda$ has little effect on the ADM mass of the spacetimeit significantly affects when and where the elliptic region forms, as is evident in these figures. Preceding formation of this elliptic region, the outgoing scalar field characteristic speeds near it become negative, akin to trapped surface formation in GR gravitational collapse. However, the spacetime outgoing null characteristic speeds $e^{A-B}$ remain positive and well away from zero throughout the integration domain. Hence, this elliptic region is not “censored” by spacetime causal structure (the ADM mass of the spacetimes are below the smallest known static black hole solutions in EdGB gravity [9, 12, 13], and even so, the elliptic regions form well outside $r=2m$, so it does not seem plausible that some spacetime trapped region could eventually form to hide the elliptic region from asymptotic view). At the sonic line bounding the hyperbolic from elliptic region, all field variables are smooth and finite. In particular, there is no geometric or scalar field singularity that might otherwise have suggested the classical theory has already ceased to give sensible predictions prior to this; see Figures 6 and 7 that show the Ricci scalar as an example.

That the character of the $(P,Q)$ subsystem is hyperbolic in some regions of the spacetime, and elliptic in others (separated by the parabolic sonic line), means the EdGB equations can be of mixed type (note of course that this is different from coupled elliptic/hyperbolic systems often encountered in GR evolution, where some equations are elliptic, others hyperbolic, but each equation maintains its definite character throughout the domain). Mixed-type equations are not as common in the literature, but do arise in several situations, such as steady transonic flow (see for example [28], which also discusses other areas where mixed type equations appear). There are two canonical mixed type equations that at least locally (near the sonic line) are expected to capture the nature of most mixed type equations : the Tricomi equation

[TABLE]

and the Keldysh equation

[TABLE]

These equations are hyperbolic/parabolic/elliptic for y<0\/\ y=0\/$\ y>0$. The main qualitative differences between these two equations are how the characteristics in the hyperbolic region meet the parabolic sonic line, and how the characteristic speeds become imaginary. For the Tricomi equation, the characteristics intersect the sonic line orthogonally, with the corresponding speeds going imaginary passing though zero there. For the Keldysh equation, the characteristics intersect the sonic line tangentially, with the characteristic speeds diverging there before becoming imaginary. This affects the degree of smoothness one can generally expect for solutions to these equations, with the Keldysh equation having weaker regularity of solutions on the sonic line (see e.g.  [28]). Though the EdGB equations are vastly more complicated than these simple prototypes, at least based on the way the characteristics intersect the sonic line, as is apparent in Figures 4 and 5, and that the characteristic speeds go to zero there, it appears that the EdGB equations are of Tricomi type. This is typical for all cases we have run where an elliptic region forms (though interestingly, for a certain class of $P(X)$ Horndeski theories in similar collapse scenarios,  [25] find either Tricomi or Keldysh behavior approaching the sonic line, depending upon the initial data in the hyperbolic region).

That the mixed type behavior here appears to be of Tricomi type is a somewhat promising sign for EdGB gravity in terms of regularity on the sonic line (as we explicitly find in our solutions); however, that an elliptic region forms regardless of its type is problematic for the theory being capable of serving as a viable, physical model that can make predictions in the sense of possessing a well-posed initial value problem (for further discussion on this see [14]).

8.2

Misner-Sharp mass and the null convergence condition

Figure 8 is a plot of the initial Misner-Sharp mass profiles for the two strong coupling ($\lambda=\pm 50$) weak field cases, together with initial data with equivalent parameters for the GR ($\lambda=0$) case. As discussed in Section 4, we may interpret $\partial_{r}m_{MS}(t,r)/4\pi r^{2}$ as an effective local energy density at $(t,r)$. As is apparent in the figure, for EdGB gravity there are cleary regions where this energy density is negative (this phenomenon has been notice before in static solutions, see e.g. [9]). Despite large variations in $m_{MS}$ in the interior as the Gauss-Bonnet coupling $\lambda$ is varied, we find that the ADM mass (estimated by evaluating the Misner-Sharp mass at $r=r_{max}$) depends much more weakly on $\lambda$. With fixed initial data ($a_{0}=0.02$, $w_{0}=8$, $r_{0}=20$), the ADM mass changes by at most 1 part in $10^{3}$ as $\lambda$ varies from -75 to 75, where we estimate the numerical error in this quantity to be less than $1$ part in $10^{4}$ (from truncation error and finite radius effects).

Related to the negative effective energy densities, we find that the NCC (3efhiqaa) is violated around these regions for the non-zero $\lambda$ cases : see Figures 9 and 10. We note that we find no correlation between the existence of negative energy density regions or regions of NCC violation and the formation of elliptic regions. While we always observe negative energy density regions and regions of NCC violation at the formation of an elliptic region, we also observe those regions in simulations where the evolution remains hyperbolic.

8.3 Convergence of simulations

In addition to convergence data we have already shown in Figures 3,7 and 10, in Figure 11 we show convergence plots from the two strong coupling ($\lambda=\pm 50$) cases for the independent residual of the $\vartheta\vartheta$ component of the EdGB equations (3a). That this converges to zero (at second order) prior to formation of an elliptic region is a rather non-trivial check of the correctness of our solution, as $E_{\vartheta\vartheta}$ depends on temporal and spatial gradients of all variables $(P,Q,A,B)$ in the problem (the EdGB equations, as GR, are over-determined, allowing for such non-trivial checks of a solution obtained from a complete subsystem of PDEs). That we loose convergence after formation of the elliptic region is consistent with the fact that we are attempting to solve a mixed type equation using hyperbolic methods, which are not well-posed in the elliptic region (for more discussion on this see [14]).

We report that in addition to the convergence tests we have discussed and presented in this paper, we achieved second order convergence before the formation of elliptic regions for all of the fields and diagnostics we implemented in our simulations, including the EdGB and null characteristics (as shown for exmaple in Figure 4), and the mass aspect, (Figure 8). Interestingly, as with the regions of NCC violation, with the resolutions reported in this paper we resolve the regions of negative energy density ($\partial_{r}m_{MS}<0$) seen in Figure 8.

9 Conclusion

In this paper we presented studies of numerical solutions of EdGB gravity in spherical symmetry in gravitational collapse-like scenarios, focusing on how properties of the solutions differ from similar situations in Einstein gravity minimally coupled to a massless scalar field. For sufficiently weak EdGB coupling we find results similar to GR : a weak field limit where the scalar field pulse disperses beyond the integration domain, and a strong field were a geometric horizon begins to form. In the latter scenarios, the EdGB scalar begins to grow outsize the nascent horizon in a manner consistent with known static “hairy” BH solutions. In the strong EdGB coupling regime, we find markedly different behavior from GR : (1) the equations of motion can be of mixed type, where an initially hyperbolic system shows development of a parabolic sonic line in a localized region of the domain, beyond which the character of the PDEs switches to elliptic (2) there are regions of negative effective energy density, and (3) there are regions where the NCC is violated. In the cases we have studied these three properties occur together within roughly the same region of spacetime. While the potential physical consequences for negative energy density and NCC violation have been extensively discussed in the modified gravity literature, the physical interpretation of mixed type equations remains largely unexplored. At the very least, mixed type behavior signals loss of predictability in the theory in the sense of it ceasing to possess a well-posed IVP.

One of our main motivations for studying EdGB gravity is to discover a viable, interesting modified gravity theory to confront with LIGO/Virgo binary BH merger data, in particular the part of the signals attributable to coalescence. In that regard, our results reported here and in a companion paper [14] do not yet rule out a coupling parameter that gives a smallest possible static BH solution of around a few solar masses (which would give the most significant differences from GR for stellar mass BH mergers), if we assume there is no cosmological background for the EdGB scalar (i.e. it is only present as sourced by curvature produced by other matter/BHs in the universe, though even then we need to ignore problems that might arise in the very early, pre-Big-Bang-Nucleosynthesis universe). So we have a tentative green light to continue this line of exploration of EdGB gravity. The next step is to solve the EdGB equations in spherical symmetry in a horizon penetrating coordinate system. This will allow us to begin addressing issues of long term, non-linear stability of hairy BHs, and perform a more thorough investigation of the strong field, strong coupling regime. Considering the qualitatively different behavior for GR we see in the EdGB simulations in the strong coupling regime, it would also be interesting to understand the nature of critical collapse in EdGB gravity, where (at least in GR) one can dynamically evolve from smooth initial data to regions of potentially unbounded curvature. We are presently working on a code to study this phenomena as well.

We thank J. Anderson, L. Lehner, V. Paschalidis, I. Rodnianski, L. Stein, and M. Taylor for useful conversations on aspects of this project. We thank the organizers of the workshop ‘Numerical Relativity beyond General Relativity’ and the Centro de Ciencias de Benasque Pedro Pascual, 2018, where we completed some of the work presented here. F.P. acknowledges support from NSF grant PHY-1607449, the Simons Foundation, and the Canadian Institute For Advanced Research (CIFAR). Computational resources were provided by the Feynman cluster at Princeton University.

Appendix A Derivation of dilaton Gauss-Bonnet tensor

Here we derive equations of motion for the dilaton Gauss-Bonnet term

[TABLE]

Varying the Gauss-Bonnet term with respect to the metric, we have

[TABLE]

We focus on the variation of the Riemann tensor term:

[TABLE]

In four dimensions, a five index antisymmetric tensor is zero, so we may write (c.f. Appendix A and B of [49])

[TABLE]

which implies

[TABLE]

We conclude that in four spacetime dimensions, the variation of the Gauss-Bonnet term with respect to the metric is

[TABLE]

Relabeling indices, and noting that from the Bianchi identities $\delta^{\rho\kappa\alpha\beta}_{\lambda\sigma\gamma\delta}R^{\lambda\sigma}{}_{\rho\kappa}$ is divergenceless on all its indices (e.g. [29] §13.5) the variation of the dilaton Gauss-Bonnet term is

[TABLE]

plus surface terms. Using similar manipulations as presented above, we note that taking the divergence of the Gauss-Bonnet tensor is

[TABLE]

so that assuming $\nabla_{\nu}\phi\neq 0$, taking the divergence of (3a) gives us (3b) (the ‘generalized Bianchi identity’ [50]).

Appendix B EdGB equations of motion

In the coordinates (3d), the nontrivial components of the EdGB equations of motion (3a) are

[TABLE]

In the language of the $3+1$ ADM formalism (3efhiqam) is the Hamiltonian constraint, (3efhiqan) is the momentum constraint, and (3efhiqao) and (3efhiqap) are part of the evolution equations for the extrinsic curvature of spacelike slices with normal vector $n^{\mu}=\left(e^{-A},0,0,0\right)$.

Equation (3a) for the EdGB scalar is

[TABLE]

and the evolution equation for the constraint $\partial_{r}\phi=Q$ is

[TABLE]

When $\lambda=0$, it is clear from (3efhiqam) and (3efhiqan) that the gravity degrees of freedom, $A$ and $B$, are fully constrained. All the dynamics are driven by (3efhiqaq). The addition of the EdGB tensor terms introduces $\partial_{t}P$ and $\partial_{t}B$ terms into the constraint equations. The Gauss-Bonnet scalar introduces second derivative terms as well as $\partial_{t}B$, $\partial_{t}A$ terms to (3efhiqaq). These new $\partial_{t}A$ and $\partial_{t}B$ terms appear to change the PDE character of the EdGB field equations versus the GR field equations. As it turns out though, we can use algebraic combinations of (3efhiqap) to remove second derivative and $\partial_{t}A$ terms, and (3efhiqan) to remove $\partial_{t}B$ terms from Equations (3efhiqam), (3efhiqan), and (3efhiqaq). Doing so leads us to Equations (3efa)-(3efd).

Appendix C

A second procedure to compute the characteristics

Here we present another procedure we used to calculate the characteristics of the propagating degree of freedom for the EdGB system in spherical symmetry. Instead of substituting in for $\partial_{r}A$ and $\partial_{r}B$ at the level of the full equations of motion, we first compute the full principal symbol and then substitute them in from the constraints. We find that this method is more numerically stable near the origin at high resolutions. This is most likely because the length the equations to be evaluated in each component of the principal symbol in this method are much shorter than they are in the other, which makes them less susceptible to floating point roundoff errors. Both methods produce equivalent results to within truncation error.

This procedure to compute the characteristics goes as follows: we consider the full system of equations (3efa)-(3efd); which take the following form

[TABLE]

where now $I,J,K$ index the fields $(A,B,Q,P)$: the equations $E^{(Q)}$ and $E^{(P)}$ retain the terms $\partial_{r}A$ and $\partial_{r}B$. The characteristic matrix for the full system is

[TABLE]

where

[TABLE]

Provided $\mathfrak{s}$ is invertible333Note that when $\lambda=0$, $\mathfrak{s}$ is the identity matrix. In practice, we have never encountered a situation where $\mathfrak{s}$ is not invertible., we can use Gaussian elimination to write the characteristic equation as

[TABLE]

where $c\equiv-\xi_{t}/\xi_{r}$, $\mathfrak{i}$ is the identity matrix, and

[TABLE]

The two characteristics given by $\xi_{r}=0$ define the characteristic surfaces for the constrained degrees of freedom. The characteristics for the dynamical degree of freedom are determined by solving the nontrivial determinant; we then find that the characteristic speeds for this degree of freedom are given by the eigenvalues of $\mathfrak{c}$,

[TABLE]

Appendix D Static decoupled EdGB solution about a

Schwarzschild black hole background

In Section 7 we compared the profile of our scalar field to that of the ‘decoupled’ EdGB scalar profile about a Schwarzschild background. For completeness we present the calculation of the profile of $\phi$. In the decoupling limit of EdGB (e.g. [6, 48]) the geometry is determined by the Einstein equations coupled to matter fields but not the EdGB scalar field, and the equation of motion for the EdGB scalar is given by

[TABLE]

We consider static solutions to this equation with a fixed Schwarzschild black hole background

[TABLE]

With this, (3efhiqbc) reduces to

[TABLE]

Imposing regularity of $\partial_{r}\phi$ at the geometric horizon $r=2M$, setting $\lim_{r\to\infty}\phi=0$, and changing variables to $x\equiv r/M$, we obtain

[TABLE]

which is what we compare against our numerical solutions in Figure 3.

Appendix E Numerical methods

We implemented three different finite difference PDE solution methods to solve equations (3efa)-(3efd), in order gain confidence that the code crashes occurring some time after formation of sonic lines are due to a property of the underlying continuum equations, rather than a numerical instability associated with a particular discretization scheme. The first two methods, described here, are fully constrained, the third is a partially constrained scheme, described below in E.5. All methods we implemented treat the ($P,Q$) subsystem as hyperbolic, and are (globally) second order accurate with fixed time and spatial steps. The two hyperbolic methods for $(P,Q$) we developed are an iterative Crank-Nicolson scheme (CN), and a fourth order in time Runge-Kutta (method of lines) scheme (RK4). We ran simulations with CFL numbers that varied from $10^{-2}$ to $0.5$. The different methods all give the same results to within truncation error, and once the elliptic region forms all crash in a qualitatively similar manner (growth of short wavelength solution components within the elliptic region at a rate proportional to their wave number; note though that the since our initial data is smooth, these short wavelength components are sourced by truncation error for the most part, and their “initial” amplitudes on the sonic line therefore decrease with resolution). This gives us confidence that the crashes are due to trying to solve a mixed type equation using hyperbolic methods, which are not well-posed in elliptic regions.

We use the notation $f^{n}_{j}$ for a discretized field, where $n$ stands for the time step and $j\in 0..N_{r}-1$ is the index within the spatial grid with $N_{r}$ points. The basic iteration loop we use for both the CN and RK4 evolution schemes, solving for the unknowns at time step $n+1$ given data at time step $n$, is as follows:

Initialize time step $n+1$ values for the fields $A$, $B$, $Q$, and $P$ with their values at time step $n$ (this step is unnecessary for the RK4 scheme). 2. 2.

For the CN scheme (E.1) perform one step of a Newton iteration to correct the unknown values of $Q^{n+1}_{j},P^{n+1}_{j}$; for the RK4 integration (E.2) take the next substep of the RK4 scheme, saving the results in temporary arrays, or $Q^{n+1}_{j},P^{n+1}_{j}$ for the final step. 3. 3.

Integrate the constraints for $A^{n+1}_{j}$ and $B^{n+1}_{j}$ given the current values of $Q^{n+1}_{j}$ $P^{n+1}_{j}$ (or the appropriate substep arrays when using RK4). Since equation (3efb) for $B$ does not depend on $A$, we first integrate this for $B$ (E.3), then substitute the result into (3efa) before integrating it for $A$ (E.4). 4. 4.

Repeat steps (ii) and (iii) until (a) for the CN iterative scheme the residuals for the full nonlinear set of equations are below a tolerance set to be a few orders of magnitude smaller than truncation error; (b) for RK4, we have completed all the RK substeps. 5. 5.

Apply a Kreiss-Oliger filter (e.g. [51]) to the now known variables $Q^{n+1}_{j}$ and $P^{n+1}_{j}$.

E.1 CN Hyperbolic PDE solver for $Q$ and

$P$

For the iterative methods we employ a Crank-Nicolson discretization in time (see e.g.[52]), where the equations (3efa) and (3efa) are discretized at a time half way between time steps $n$ and $n+1$, which we denote as time step $n+1/2$. Explicitly, we replace each field $f$ and its gradients with the following stencils

[TABLE]

We define the residual and field vectors $\mathcal{R}_{k}$ and $v_{k}$ respectively via

[TABLE]

where $0<k<2(N_{r}-1)$. For the iteration step (ii) above we compute the linear correction $\delta v_{j}$ by solving the following matrix equation

[TABLE]

for $\delta v_{j}$, where

[TABLE]

We invert the matrix $\mathcal{J}_{ij}$ in two different ways. For the first method we directly solve (3efhiqbn) with a banded matrix solver (the LAPACK routine dgbsv [53]). For the second method we solve (3efhiqbn) with Gauss-Seidel iteration (e.g. [54]).

E.2 RK4 PDE solver for $Q$ and

$P$

We use a standard fourth order in time Runge-Kutta algorithm (see e.g. [32]), so will not describe it here, but note that we still only employ a second order accurate discretization for spatial gradients; i.e. for each field $f$ we use the stencils

[TABLE]

For this study we are able to achieve the requisite accuracy with second order methods and reasonable computer power, though do not use a second order Runge-Kutta method, with the radial differences (3efhiqbp) and (3efhiqbq), as it is unconditionally unstable for the linear wave equation, as may be verified by a von-Neumann stability analysis.

E.3 ODE integrator for $B$

Equation (3efb) for the $B$ field takes the schematic form

[TABLE]

where both $c_{(B)}$ and $d_{(B)}$ are nonlinear functions of $B$, $P$, $Q$, and the radial derivatives of $P$, and $Q$. We solved this equation in two different ways. The first involves Newton’s method: we define the vectors $\mathcal{R}_{j}$ and $v_{j}$, with $0\geq j\geq N_{r}-1$ and

[TABLE]

where $(E^{(B)})^{n+1}_{j+1/2}$ is (3efb) with the fields finite differenced using the trapezoid stencil:

[TABLE]

Equation (3efb) is nonlinear in $B$, so we iteratively solve for $B_{j}$ by solving for the linear correction $\delta v_{j}$ in

[TABLE]

for $\delta v_{j}$, where

[TABLE]

As in E.1 we inverted $\mathcal{J}_{ij}$ two different ways: one using a banded matrix solver, and another iteratively using a Gauss-Seidel method. The Newton iteration was then repeated until the residual $\mathcal{R}_{j}$ was below some tolerance well below truncation error.

We also directly solved Equation (3efhiqbr) using a second order Runge-Kutta method, by writing the equation as $\partial_{r}B=-d/c$.

E.4 ODE integrator for $A$

The ODE for the $A$ field, (3efa) is of the form

[TABLE]

where $c_{(A)}$ and $d_{(A)}$ are functions of $B$, $P$, $Q$, and their radial derivatives. We discretize the fields and their derivatives using the trapezoidal rule as above (3efhiqbu,3efhiqbv). Since the ODE for $A$ is linear it is trivial to directly integrate it from the origin $j=0$ outward; specifically we directly solve for $A^{n+1}_{j+1}$ knowing $A^{n+1}_{j}$ and the other field values via

[TABLE]

E.5 Partially Constrained Evolution

In a partially constrained evolution, one (or more) variables are typically solved for using an evolution instead of constraint equation. Here, one can do that for $B$, with (3efhiqap) the corresponding second-order-in-time evolution equation for it. However, in Schwarzschild-like coordinates the momentum constraint (3efhiqan) is effectively a “first integral” for this equation, and instead then we consider this as our evolution equation for $B$ (recall for our constrained evolution we do not use the plain form of the momentum constraint, but first eliminate the time derivative of $B$ using the other equations). For initial data, we solve for $B$ using (3efb) with either an RK2 or a relaxation method. Once we begin evolving in time, we then use (3efb) as an independent residual to monitor the constraint.

We solved a discretized version of (3efhiqan) for $B$ using an iterative Crank-Nicolson method. On any given time step, we follow a similar procedure as above for the iterative constrained scheme, but now iterate over the evolution equations for $Q^{n+1}$, $P^{n+1}$, and $B^{n+1}$ a fixed number of times, then solve for $A^{n+1}$ using the constraint equation, (3efa). We repeat this process until the residuals of the evolutions equations for $Q$, $P$, and $B$ are below a tolerance set to be a few orders of magnitude below the truncation error. Afterward we apply a Kreiss-Oliger filter on the variables $Q^{n+1}$, $P^{n+1}$, and $B^{n+1}$, before advancing to the next time step.