††thanks:
The two authors contributed equally. The order of authors is listed
alphabetically.

# Nonlinear Dynamics of Quadratic Gravity in Spherical Symmetry

Aaron Held<EMAIL_ADDRESS>Blackett Laboratory, Imperial College
London, London SW7 2AZ, United Kingdom  Hyun Lim<EMAIL_ADDRESS>Computational Physics and Methods (CCS-2), Los Alamos National Laboratory, Los
Alamos, NM 87545 USA  Center for Theoretical Astrophysics, Los Alamos National
Laboratory, Los Alamos, NM 87545 USA  Applied Computer Science (CCS-7), Los
Alamos National Laboratory, Los Alamos, NM 87545 USA

###### Abstract

We present the first numerically stable nonlinear evolution for the leading-
order gravitational effective field theory (Quadratic Gravity) in the
spherically-symmetric sector. The formulation relies on (i) harmonic gauge to
cast the evolution system into quasi-linear form (ii) the Cartoon method to
reduce to spherical symmetry in keeping with harmonic gauge, and (iii) order-
reduction to 1st-order (in time) by means of introducing auxiliary variables.
Well-posedness of the respective initial-value problem is numerically
confirmed by evolving randomly perturbed flat-space and black-hole initial
data. Our study serves as a proof-of-principle for the possibility of stable
numerical evolution in the presence of higher derivatives.

##  I Introduction

General Relativity (GR) is in excellent agreement with an ever-growing body of
experimental tests. At the same time, theoretical considerations strongly
suggest that the theory is incomplete: It is not known how to consistently
couple GR to quantum field theories of matter once back-reactions are non-
negligible. Moreover, even classical matter distributions are prone to
gravitational collapse and thus to the formation of curvature singularities,
typically accompanied by geodesic incompleteness. See Penrose (1965); Hawking
(1965); Geroch (1966); Hawking and Penrose (1970) for singularity theorems in
static and highly symmetric settings and, e.g., Dafermos (2005); Harada et al.
(2002) for numerical explorations in less symmetric and dynamical settings.

Said theoretical breakdowns strongly motivate to embed GR into the modern
framework of effective field theory (EFT) and treat the Einstein-Hilbert
action merely as the leading-order term in a local-curvature expansion of a
general diffeomorphism-invariant action of gravitational (and matter) degrees
of freedom. Inconsistencies in the coupling to matter and the formation of
singularities can then be interpreted as a consequence of extrapolating the
EFT beyond its regime of validity: As the curvature grows during gravitational
collapse, higher-order terms in the EFT will eventually become non-negligible
and alter the dynamics of GR at some, as of now untested, curvature scale. We
may encounter this scale anywhere between the largest currently probed
curvature scales and the Planck scale. Indeed, quantum fluctuations are widely
expected to induce such EFT curvature corrections, cf. ’t Hooft and Veltman
(1974); Stelle (1977); Goroff and Sagnotti (1985); Avramidi and Barvinsky
(1985); van de Ven (1992) for perturbative quantum gravity, Boulware and Deser
(1985); Zwiebach (1985) for string theory, and Percacci (2017); Eichhorn
(2019); Reuter and Saueressig (2019) for reviews in the context of asymptotic
safety.

Since GR tends to hide regimes of growing curvature behind event horizons
Penrose (1969), experimentally probing the horizon-scale physics of black
holes is presumably one of the most promising ways to push the limits of the
EFT of gravity. A rapidly increasing number of gravitational-wave (GW) events
from black-hole binary mergers Abbott et al. (2016, 2020) provide access to
this, previously uncharted, strong-gravity regime.

Utilizing this data to constrain new physics beyond GR requires obtaining
alternative predictions for the leading-order corrections in the above
gravitational EFT. As large-curvature regimes reveal the nonlinear character
of (beyond) GR dynamics, such predictions require numerical relativity
simulations, cf. Shibata and Nakamura (1995); Baumgarte and Shapiro (1998);
Pretorius (2005a, b) for pioneering work in numerical GR. It is, therefore,
crucial to find a well-posed numerical evolution for said leading-order EFT
corrections. The existence of a well-posed initial value problem (IVP), see
Sarbach and Tiglio (2012); Isenberg (2014) for reviews in GR, could pose a
restriction for any viable theory. Certainly, a well-posedness IVP (for
physically meaningful sets of initial data) is necessary to perform any stable
numerical evolution.

First numerical simulations have been achieved in specific beyond-GR theories
such as dynamical Chern-Simons Okounkova et al. (2019, 2020), Einstein-
dilaton-Gauss-Bonnet Ripley and Pretorius (2020a, b); Okounkova (2020); East
and Ripley (2021); Silva et al. (2020), Horndeski theories Figueras and França
(2020), K-essence Bezares et al. (2021), and a study of EFT terms at quadratic
order in curvature Cayuso and Lehner (2020). A stable numerical evolution is
guaranteed by either iterative treatment Okounkova et al. (2019, 2020);
Okounkova (2020), a dampening of high-frequency modes Cayuso and Lehner
(2020), or by an established well-posed evolution at weak coupling Kovács
(2019); Kovács and Reall (2020a, b); East and Ripley (2021). In contrast, at
strong coupling the onset of ill-posed regimes has been observed in Ripley and
Pretorius (2020a, b).

Here, we investigate Quadratic Gravity (QG), i.e., the gravitational EFT
including all independent terms up to quadratic order in curvature – sometimes
also referred to as Stelle-gravity Stelle (1977, 1978). At the formal level,
it has been shown – without any constraint to weak coupling – that QG admits a
well-posed IVP Noakes (1983), see also Morales and Santillán (2019).

Accounting for EFT corrections to GR naturally implies higher-order equations
of motion. The latter are theoretically disfavored by the Ostrogradski theorem
Ostrogradsky (1850) which states that non-degenerate higher-derivative
theories result in linearized degrees of freedom with opposite-sign kinetic
terms. Any non-vanishing coupling between these modes implies the onset of a
linear instability, cf. eg. Becker et al. (2017); Anselmi and Piva (2018);
Donoghue and Menezes (2019); Salvio (2019) for related recent developments in
the context of a unitary quantum evolution of QG.

We emphasize that Ostrogradski (in)stability and well-posedness are not
necessarily related. On the one hand, the Ostrogradski theorem is a physical
statement: theories that exhibit opposite-sign kinetic terms can develop
physical instabilities. On the other hand, well-posedness is a mathematical
property of partial-differential equations (PDEs), crucial for numerical
simulations but not necessarily related to physical implications.

In fact, the present work can be viewed as a proof of principle for a higher-
derivative gravitational theory with Ostrogradski ghosts which nevertheless
admits a well-posed IVP and in which we can thus simulate spacetime dynamics
numerically. The proof in Noakes (1983) establishes the existence of a well-
posed IVP at the level of the four-dimensional equations of motion. As in GR,
it remains non-trivial to translate it into (3+1) form to obtain a well-posed
IVP suitable for actual numerical evolution. Here, we do so in the reduced
sector of spherically symmetric dynamics. This allows us to present the first
nonlinear evolution of physical initial data in QG.

The rest of the paper is organized as follows: in Section II, we review the
well-posed IVP formulation of QG, cf. Noakes (1983); in Section III, we use
the Cartoon method Alcubierre et al. (2001); Pretorius (2005a) to reduce the
evolution equations to spherical symmetry and perform the order-reduction in
the symmetry-reduced setup; in Section IV, we present the resulting numerical
evolution which we find to be well-posed, both for perturbations of flat
spacetime as well as for perturbations of the Schwarzschild solution; in
Section V we conclude with a summary and discuss the implications of our
results for future works.

As for conventions, we use the $(-,+,+,+)$ signature, we work in geometrized
units where ($c=1$, $G=1$), and use Latin letters as Lorentzian spacetime
indices. Round (square) brackets denote (anti-)symmetrization of the enclosed
pair of indices.

##  II Quadratic Gravity and the Noakes equations

Quadratic Gravity (QG) incorporates the leading-order, i.e., curvature-
squared, corrections to GR and can be parameterized by the action

$\displaystyle S_{\text{QG}}=\int d^{4}x\sqrt{|g|}\left[\frac{1}{16\pi
G}R+\alpha R_{ab}R^{ab}-\beta R^{2}\right]\;.$ (1)

Note that we have chosen the minus sign in the last term to agree with
conventions in the equations of motion in Noakes (1983). The first term is the
common Einstein-Hilbert term where we have, for the present section, re-
instated Newton’s constant $G$. In 4D, the most general corrections of
quadratic order in curvature can be parameterized by the two dimensionless
constants $\alpha$ and $\beta$. A potential Riemann-squared term can be
rewritten into the former through the Gauss-Bonnet identity.

The associated equations of motion for quadratic gravity (QG) contain (up to)
4th-order derivative terms Stelle (1977) and read

$\displaystyle 2\,T_{ab}=$
$\displaystyle(\alpha-2\beta)\nabla_{a}\nabla_{b}R-\alpha\Box
R_{ab}-\left(\frac{1}{2}\alpha-2\beta\right)g_{ab}\Box R$
$\displaystyle+2\alpha R^{cd}R_{acbd}-2\beta RR_{ab}-\frac{1}{2}(\alpha
R^{cd}R_{cd}-\beta R^{2})$ $\displaystyle+\frac{1}{16\pi
G}\left(R_{ab}-\frac{1}{2}g_{ab}R\right)\;.$ (2)

Here, we have included an energy-momentum tensor $T_{ab}$ for potential matter
sources.

In the linearized theory, it has been shown Stelle (1977) that, in addition to
the massless spin-2 mode of GR, the linearized Ricci scalar and linearized
traceless Ricci tensor propagate a massive scalar and massive spin-2 modes
with respective masses

$\displaystyle m_{0}^{2}=\frac{1}{32\pi G(3\beta-\alpha)}\;,\quad\quad\quad
m_{2}^{2}=\frac{1}{16\pi G\alpha}\;.$ (3)

The massive spin-2 mode is an Ostrogradski ghost, i.e., in the linearized
theory, its kinetic term has the opposite sign in comparison to the massless
graviton. In Noakes (1983), Noakes finds that the nonlinear evolution can be
formulated with the same degrees of freedom. Following this insight, the Ricci
scalar $\mathcal{R}\equiv g^{ab}R_{ab}$ and the traceless Ricci tensor
$\widetilde{\mathcal{R}}_{ab}\equiv R_{ab}-1/4\,g_{ab}\,R$ can be elevated to
independent variables, as indicated by the curly notation. From here on, these
should no longer be evaluated on the metric. Rather, they are treated as
independent evolution variables. Re-expressing the parameters $\alpha$ and
$\beta$ by the two masses $m_{0}^{2}$ and $m_{2}^{2}$, the equations of motion
can be separated into a trace (2nd equation below) and a traceless (3rd
equation below) part, i.e.,

$\displaystyle R_{ab}(g)=$
$\displaystyle\;\widetilde{\mathcal{R}}_{ab}+\frac{1}{4}g_{ab}\mathcal{R}\;,$
(4) $\displaystyle\Box\,\mathcal{R}=$
$\displaystyle\;m_{0}^{2}\mathcal{R}+2T^{c}_{\phantom{c}c}\;,$ (5)
$\displaystyle\Box\,\widetilde{\mathcal{R}}_{ab}=$
$\displaystyle\;m_{2}^{2}\widetilde{\mathcal{R}}_{ab}+2T^{\text{(TL)}}_{ab}$
$\displaystyle-\frac{1}{3}\left(\frac{m_{2}^{2}}{m_{0}^{2}}-1\right)\left(\nabla_{a}\nabla_{b}\mathcal{R}-\frac{1}{4}g_{ab}m_{0}^{2}\mathcal{R}\right)$
$\displaystyle+2\widetilde{\mathcal{R}}^{cd}C_{acbd}-\frac{1}{3}\left(\frac{m_{2}^{2}}{m_{0}^{2}}+1\right)\mathcal{R}\widetilde{\mathcal{R}}_{ab}$
$\displaystyle-2\widetilde{\mathcal{R}}_{a}^{\phantom{a}c}\widetilde{\mathcal{R}}_{bc}+w\frac{1}{2}g_{ab}\widetilde{\mathcal{R}}^{cd}\widetilde{\mathcal{R}}_{cd}\;.$
(6)

Here, we have supplemented the trace and traceless equations by the definition
of the Ricci curvature in terms of the metric (on the left-hand side) and in
terms of the fiducial variables (on the right-hand side), cf. Eq. (4). The
latter provides a 2nd-order evolution equation for the metric in which the
fiducial variables appear as sources, i.e., without derivatives. The matter
sources are also split into trace and traceless part (indicated by a (TL)
superscript) and, in turn, only source the fiducial variables. Furthermore, we
have introduced the Weyl-tensor $C_{abcd}$ for brevity. The latter can
equivalently be expressed in terms of $R_{abcd}$,
$\widetilde{\mathcal{R}}_{ab}$, and $\mathcal{R}$ as

$\displaystyle
C_{abcd}=R_{abcd}+g_{b[c}\widetilde{\mathcal{R}}_{a]d}+g_{d[a}\widetilde{\mathcal{R}}_{c]b}+\frac{1}{6}g_{a[d}g_{b]c}\mathcal{R}\;.$
(7)

Therefore, the metric appears in the trace and traceless equation, cf. Eqs.
(5)-(II), only as part of the covariant derivatives as well as in $R_{abcd}$.

The trace and traceless Ricci variables $\mathcal{R}$ and
$\widetilde{\mathcal{R}}_{ab}$ appear merely as ‘fiducial sources’, i.e., as
terms without derivatives. Thus, the metric equation, cf. Eq. (4), can be
treated as in GR where (generalized) harmonic coordinates allow to express the
Ricci tensor as a strongly hyperbolic quasi-linear 2nd-order differential
operator. More explicitly, the Ricci tensor on the left-hand side of Eq. (4)
can be expressed as

$\displaystyle
R_{ab}(g)=-\frac{1}{2}g^{cd}g_{ab,cd}+g_{c(a}F^{c}_{\phantom{c},b)}+\mathcal{O}^{1}_{ab}(\partial
g)\;,$ (8)

where we use the usual comma-notation for partial derivatives and
$F^{a}=-g^{cd}\Gamma^{a}_{cd}$ is a gauge potential which, provided one works
in (generalized) harmonic gauge, i.e., $F^{a}=0$ ($F^{a}=\text{const}$),
reveals the strongly hyperbolic quasi-linear character. Additional lower-order
derivative terms are denoted by $\mathcal{O}^{1}_{ab}(\partial g)$.

With (generalized) harmonic coordinates at hand and by expanding the Riemann
tensor as well as all covariant derivatives in terms of the metric, the
evolution of the variables $g_{ab}$, $\mathcal{R}$, and
$\widetilde{\mathcal{R}}_{ab}$ can be written as

$\displaystyle g^{cd}g_{ab,cd}$
$\displaystyle=-2\widetilde{\mathcal{R}}_{ab}-\frac{1}{2}g_{ab}\mathcal{R}+\mathcal{O}^{1}_{ab}(\partial
g)\;,$ (9) $\displaystyle g^{cd}\mathcal{R}_{,cd}$
$\displaystyle=m_{0}^{2}\mathcal{R}\;,$ (10) $\displaystyle
g^{cd}\widetilde{\mathcal{R}}_{ab,cd}$
$\displaystyle=\mathcal{O}^{2}_{ab}(\partial\partial\mathcal{R},\partial\widetilde{\mathcal{R}},\partial\partial
g)\;.$ (11)

Again, we summarize lower-order terms with $\mathcal{O}^{1}_{ab}(\partial g)$
and
$\mathcal{O}^{2}_{ab}(\partial\partial\mathcal{R},\partial\widetilde{\mathcal{R}},\partial\partial
g)$ with the notation in brackets indicating the highest order of derivatives
of these terms. As will become clear below, their explicit form is irrelevant
regarding well-posedness. Nevertheless, we explicitly provide their form in
App. A, where we also correct some typos in comparison to Noakes (1983);
Morales and Santillán (2019).

The above system is not yet of _diagonal_ quasi-linear form due to both
$\partial\partial\mathcal{R}$\- and $\partial\partial g$-terms appearing in
the traceless equation, cf. Eq. (11). However, the system is amenable to
diagonalization because of the lack of
$\partial\widetilde{\mathcal{R}}$-contributions in the (already diagonal-form)
equations for $\mathcal{R}$ (Eq. (10)) and $g$ (Eq. (9)), see Noakes (1983)
for more formal and general statements. More explicitly, the given PDEs can be
diagonalized by introducing extra variables $V_{a}\equiv\mathcal{R}_{,a}$ and
$h_{abc}\equiv g_{ab,c}$ and adding derivatives of the former two equations to
the system, i.e.,

$\displaystyle g^{mn}V_{a,mn}$ $\displaystyle=\mathcal{O}_{a}(\partial
V,h)\;,$ (12) $\displaystyle g^{mn}h_{abc,mn}$
$\displaystyle=\mathcal{O}_{abc}(\partial h)\;,$ (13) $\displaystyle
g^{mn}\widetilde{\mathcal{R}}_{ab,mn}$
$\displaystyle=\mathcal{O}^{2}_{ab}(\partial V,\partial
h,\partial\widetilde{\mathcal{R}})\;.$ (14)

This extended system is now of diagonal quasi-linear form and the standard
theorems for hyperbolicity Leray and Institute for Advanced Study
(Princeton(1953) apply, cf. Noakes (1983) for a more detailed proof and our
separate publication Lim and Held (2021) for a full (3+1)-formulation of QG.

Treating $\widetilde{\mathcal{R}}_{ab}$ and $\mathcal{R}$ as independent
variables, Eqs. (4)-(II) only describe the physical evolution of QG if
additional constraints guarantee that $\widetilde{\mathcal{R}}_{ab}$ and
$\mathcal{R}$ equate to their metric counterparts. This can be captured by
introducing a constraint variable
$\mathcal{C}_{ab}=G_{ab}-\widetilde{\mathcal{R}}_{ab}+1/4\,g_{ab}\mathcal{R}$,
describing the deviation of the fiducial variables from the physical Einstein
tensor. Demanding that $\mathcal{C}_{ab}$ and its first time derivative vanish
(Bianchi constraints), ensures that the initial data is physical. Similarly,
since the above formulation as a quasi-linear system requires harmonic
coordinates, we need to ensure that the initial data does so too. Overall, the
initial-data constraints read

$\displaystyle F^{a}$ $\displaystyle=0\;,$ (15)
$\displaystyle\pounds_{n}F^{a}$ $\displaystyle=0\;,$ (16)
$\displaystyle\mathcal{C}_{a\phantom{b};b}^{\phantom{a}b}$
$\displaystyle=0\;,$ (17)
$\displaystyle(\pounds_{n}\mathcal{C}\indices{{}_{a}^{b}})$
$\displaystyle=0\;,$ (18)

where $\pounds_{n}$ denotes Lie derivatives along a timelike normal vector
$n^{a}$ orthogonal to the initial-data hypersurface. A similar order-reduction
as for the evolution equations above, cf. Noakes (1983), ensures that also the
propagation of the constraints can be written as a quasilinear diagonal 2nd-
order system.

In addition, initial data has to obey the usual Gauss-Codazzi (or Hamiltonian
and shift) constraints of GR which read

$\displaystyle n^{a}n^{b}\mathcal{C}_{ab}$ $\displaystyle=0\;,$ (19)
$\displaystyle n^{a}h_{c}^{b}\mathcal{C}_{ab}$ $\displaystyle=0\;.$ (20)

Finally, also the physical constraints are preserved in the temporal direction
(secondary constraints).

Counting the number of pieces of independent initial data reveals that the
nonlinear theory retains the same number of degrees of freedom as the
linearized one. Eqs. (4)-(5) constitute a set of 20 2nd-order evolution
equations (equivalent to the 10 4th-order equations in Eq. (II)). At the same
time, the harmonic constraints (Eqs. (15)-(16)), the Bianchi constraints (Eqs.
(17)-(18)), and the Gauss-Codazzi constraints (Eqs. (19)-(20)) require $8$
relations among the initial data (including the respective secondary
constraints), each. Overall this amounts to $40-24=16$ pieces of independent
initial data, corresponding to eight degrees of freedom, i.e., a massless
spin-2 (two degrees of freedom), a massive scalar (one degree of freedom), and
a massive spin-2 mode (five degrees of freedom). We will keep track of this
counting in the order reduction and the reduction to spherical symmetry below.

##  III Reduction to spherical symmetry via the Cartoon method

Unfortunately, choosing coordinates in which spherical symmetry is explicit
makes it impossible to maintain the harmonic gauge condition $\Box x^{a}=0$
Sorkin and Choptuik (2010). However, as we have seen in the last section, the
latter is a crucial to achieve a well-posed formulation, cf. discussion below
Eq. (8). We, therefore, remain in Cartesian coordinates $(t,x,y,z)$ and follow
the Cartoon method Alcubierre et al. (2001); Pretorius (2005a) to make use of
the Killing vector fields associated with spherical symmetry, i.e.,

$\displaystyle\xi_{1}^{\mu}$
$\displaystyle=x(\boldsymbol{\partial}_{y})^{\mu}-y(\boldsymbol{\partial}_{x})^{\mu}\;,$
$\displaystyle\xi_{1}^{\mu}$
$\displaystyle=y(\boldsymbol{\partial}_{z})^{\mu}-z(\boldsymbol{\partial}_{y})^{\mu}\;,$
$\displaystyle\xi_{1}^{\mu}$
$\displaystyle=z(\boldsymbol{\partial}_{x})^{\mu}-x(\boldsymbol{\partial}_{z})^{\mu}\;.$

Note that the bold font $\boldsymbol{\partial}_{x,y,z}$ indicate the basis
one-forms and _not_ partial derivatives. Expanding the respective vanishing
Lie-derivatives $\mathcal{L}_{\xi_{i}}\boldsymbol{X}=0$ acting on tensorial
objects $\boldsymbol{X}$, one can re-express partial derivatives in two of the
spatial directions, for instance $(y,z)$, in terms of the third, for instance
$x$. For scalars $\Phi$, vectors $\Psi_{a}$ and tensors $\Pi_{ab}$ which are
preserved along the Killing vector field, one finds

$\displaystyle\partial_{y}\Phi$ $\displaystyle=\frac{y}{x}\partial_{x}\Phi\;,$
$\displaystyle\partial_{y}\Psi_{a}$
$\displaystyle=\frac{1}{x}\Big{(}y\,\partial_{x}\Psi_{a}+\Psi_{x}\,\delta_{a}^{y}-\Psi_{y}\,\delta_{c}^{x}\Big{)}\;,$
$\displaystyle\partial_{y}\Pi_{ab}$
$\displaystyle=\frac{1}{x}\Big{(}y\,\partial_{x}\Pi_{ab}-2x\,\delta^{x}_{(a}\Pi_{b)y}+2y\,\delta^{y}_{(a}\Pi_{b)x}\Big{)}\;,$

and equivalently for $(y\leftrightarrow z)$. These relations allow us to
reduce all spatial derivatives in Eqs. (9)-(11) to those with respect to a
single coordinate, e.g., $x$.

Furthermore, a coordinate transformation from a coordinate system in which the
spherical symmetry is manifest, i.e., transforming
$\overline{X}=(t,r,\theta,\phi)$ back to Cartesian coordinates $X=(t,x,y,z)$,
gives symmetry relations between tensor components via

$\displaystyle\Pi_{\overline{a}\overline{b}}=\frac{\partial
X^{a}}{\partial\overline{X}^{\overline{a}}}\frac{\partial
X^{b}}{\partial\overline{X}^{\overline{b}}}\Pi_{ab}\;.$ (21)

For the present case of spherical symmetry, transformation of the symmetry
identities $\Pi_{t\theta}=0$, $\Pi_{t\phi}=0$, $\Pi_{r\theta}=0$,
$\Pi_{r\phi}=0$, $\Pi_{\theta\phi}=0$, and
$\Pi_{\theta\theta}\sin^{2}\theta=\Pi_{\phi\phi}$ back to Cartesian
coordinates implies the relations

$\displaystyle\Pi_{ty}$
$\displaystyle=\frac{y\,\Pi_{tx}}{x}\;,\quad\quad\quad\quad\;\,\Pi_{tz}=\frac{z\,\Pi_{tx}}{x}\;,$
(22) $\displaystyle\Pi_{xy}$
$\displaystyle=\frac{xy\left(\Pi_{xx}-\Pi_{yy}\right)}{x^{2}-y^{2}}\;,\;\Pi_{xz}=\frac{xz\left(\Pi_{xx}-\Pi_{yy}\right)}{x^{2}-y^{2}}\;,$
(23) $\displaystyle\Pi_{yz}$
$\displaystyle=\frac{yz\left(\Pi_{xx}-\Pi_{yy}\right)}{x^{2}-y^{2}}\;,$ (24)
$\displaystyle\Pi_{zz}$
$\displaystyle=\frac{(x^{2}-z^{2})\Pi_{yy}-(y^{2}-z^{2})\Pi_{xx}}{x^{2}-y^{2}}\;.$
(25)

Naturally, all of the above also holds for raised indices. For
$\widetilde{\mathcal{R}}_{ab}$, we can additionally make use of tracelessness
to remove one further component.

With the above relations, the evolution equations in Eqs. (9)-(11) (with the
explicit form of lower-order terms provided in App. A) can be expanded into
evolution equations for only eight independent variables which we group into

$\displaystyle\mathbf{u}=(\mathcal{R},g_{tt},g_{tx},g_{xx},g_{yy})\;\text{and}\;\mathbf{v}=(\widetilde{\mathcal{R}}_{tt},\widetilde{\mathcal{R}}_{tx},\widetilde{\mathcal{R}}_{xx}),$
(26)

according to whether their associated evolution equations are already quasi-
linear or not. The set of eight 2nd-order equations takes the form

$\displaystyle\partial_{t}^{2}\mathbf{u}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\partial_{t}{\mathbf{u}}\right)\;,$
(27) $\displaystyle\partial_{t}^{2}\mathbf{v}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\partial_{t}{\mathbf{u}},\,\partial_{t}{\mathbf{v}},\,\partial_{t}^{2}{\mathbf{u}}\right)\;.$
(28)

Unfortunately, the explicit expressions are too large to meaningfully be
displayed here. Instead, we provide them, along with all the subsequent
reductions, in ancillary files 111See the GitHub repository
(https://github.com/aaron-hd/QG-sphSymm-ancillary) for the Mathematica Inc.
script and ancillary files. Parts of the derivation make use of the xAct
package Martin-Garcia et al. (2007) (http://www.xact.es/). and restrict the
following discussion to a schematic form. We note that, in distinction to Sec.
II, the above notation only keeps track of the order of time derivatives such
that all instances of $\mathcal{O}$ can potentially contain up to 2nd-order
spatial derivatives.

The above 2nd-order evolution system naively propagates 16 free initial data
functions, i.e., 8 degrees of freedom. These are subject to _physical
constraints_ , which we will come back to in Sec. III.3. For now, we only want
to keep track of additional _auxiliary constraints_ which appear due to the
order-reduction below.

###  III.1 Reduction to quasi-linear 2nd-order form

In analogy to the 4D-diagonalization procedure in Noakes (1983), we introduce
additional auxiliary variables

$\displaystyle\dot{\mathbf{u}}=(\dot{\mathcal{R}},\,\dot{g}_{tt},\,\dot{g}_{tx},\,\dot{g}_{xx},\,\dot{g}_{yy})\equiv\partial_{t}\mathbf{u}$
(29)

and differentiate the first two equations by time. Adding the resulting
equations to the evolution system results in

$\displaystyle\partial_{t}^{2}\dot{\mathbf{u}}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}},\,\partial_{t}\dot{\mathbf{u}},\,\partial_{t}\mathbf{v}\right)\;,$
(30) $\displaystyle\partial_{t}^{2}\mathbf{v}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}},\,\partial_{t}\dot{\mathbf{u}},\,\partial_{t}\mathbf{v}\right)\;,$
(31) $\displaystyle\partial_{t}\mathbf{u}$
$\displaystyle\equiv\dot{\mathbf{u}}\;,$ (32)
$\displaystyle\partial_{t}\dot{\mathbf{u}}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}}\right)\;.$
(33)

In the above, we have added the definitions
$\dot{\mathbf{u}}\equiv\partial_{t}\mathbf{u}$ to the evolution equations.
Indeed, these defining equations for the auxiliary variables act as 1st-order
evolution equations for $\mathbf{u}$ while the original evolution equations
become 1st-order _auxiliary constraints_ on $\dot{\mathbf{u}}$ since there are
other 2nd-order evolution equations for the latter. Thereby, the evolution
system remains consistent with $21-5=16$ free functions of initial data.

###  III.2 Reduction to 1st-order form

By introducing eight further auxiliary variables,

$\displaystyle\ddot{\mathbf{u}}\equiv\partial_{t}\dot{\mathbf{u}}\quad\text{and}\quad\dot{\mathbf{v}}\equiv\partial_{t}\mathbf{v}\;,$
(34)

the set of evolution equations is cast into a from in which it becomes purely
1st-order in time, i.e.,

$\displaystyle\partial_{t}\ddot{\mathbf{u}}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}},\,\ddot{\mathbf{u}},\,\dot{\mathbf{v}}\right)\;,$
(35) $\displaystyle\partial_{t}\dot{\mathbf{v}}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}},\,\ddot{\mathbf{u}},\,\dot{\mathbf{v}}\right)\;,$
(36) $\displaystyle\partial_{t}\dot{\mathbf{u}}$
$\displaystyle\equiv\ddot{\mathbf{u}}\;,$ (37)
$\displaystyle\partial_{t}\mathbf{u}$ $\displaystyle\equiv\dot{\mathbf{u}}\;,$
(38) $\displaystyle\partial_{t}\mathbf{v}$
$\displaystyle\equiv\dot{\mathbf{v}}\;,$ (39) $\displaystyle\ddot{\mathbf{u}}$
$\displaystyle=\mathcal{O}\left(\mathbf{u},\,\mathbf{v},\,\dot{\mathbf{u}}\right)\;.$
(40)

These evolution equations, which we have supplemented by the defining
equations for the auxiliary variables $\dot{\mathbf{u}}$, $\dot{\mathbf{v}}$,
and $\ddot{\mathbf{u}}$, now denote $21$ 1st-order evolution equations.

It is now apparent that the original evolution equations for $\mathbf{u}$ take
the role of five _auxiliary constraints_ , cf. Eq. (40), such that the overall
number of free functions to be specified as initial data remains $21-5=16$.

We emphasize, once more, that we have hidden up-to-2nd-order spatial
derivatives in the notation such that the auxiliary constraint in Eq. (40) is
by no means algebraic but rather a 2nd-order spatial ODE.

A formal proof of well-posedness of the above evolution equations would either
require a further reduction of the remaining higher-order spatial derivatives,
or to analyse strong hyperbolicity as a first-order in time and arbitrary-
order in space (FTNS) system along the lines of Gundlach and Martin-Garcia
(2006); Richter and Hilditch (2011); Hilditch and Richter (2015). Instead of
pursuing such a formal proof, we focus on demonstrating a stable numerical
evolution in Sec. IV.

###  III.3 Physical Constraints

In spherical symmetry, the 24 physical constraints reduce to 12 since in each
four-component constraint equation two Cartoon relations (arising from a
relation equivalent to Eq. (21) but for vectors) can be used. We are therefore
left with $21-5(\text{auxiliary})-12(\text{physical})=4$ independent pieces of
initial data, i.e., two degrees of freedom.

More explicitly, the Hamiltonian and shift constraint reduce to

$\displaystyle\mathcal{C}_{tt}$ $\displaystyle\equiv
G_{tt}-\widetilde{\mathcal{R}}_{tt}+\frac{1}{4}\,g_{tt}\mathcal{R}=0\;,$ (41)
$\displaystyle\mathcal{C}_{tx}$ $\displaystyle\equiv
G_{tx}-\widetilde{\mathcal{R}}_{tx}+\frac{1}{4}\,g_{tx}\mathcal{R}\;.$ (42)

These physical constraints reproduce (some components of) the auxiliary
constraints, cf., Eq. (40), upon using $R(g)=\mathcal{R}$. Since
$R(g)=\mathcal{R}$ and $\widetilde{R}_{ab}(g)=\widetilde{\mathcal{R}}_{ab}$
these are, of course, the very relations which are supposed to be enforced by
the auxiliary constraints. This implies that Hamiltonian and momentum
constraints are automatically fulfilled once the auxiliary constraints are
fulfilled.

In the subsequent numerical analysis, we will monitor the Hamiltonian
constraint to confirm that the numerical evolution is well-posed.

Solving the constraints, both physical and auxiliary, is a non-trivial task
and will be addressed in future work. Below, we will restrict to perturbations
of exact solutions of the theory for which all constraints are fulfilled.
Moreover, we restrict to vacuum solutions of GR which, since they are Ricci-
flat, are exact solutions of QG as well.

##  IV Well-posed numerical evolution

Having derived a set of 1st-order (in time) evolution equations for the
spherically-symmetric sector, cf. Eqs. (35)-(39), we now proceed to solve
these numerically. Looking at perturbations of (i) flat spacetime, cf. Sec.
IV.2, and (ii) Schwarzschild spacetime, cf. Sec. IV.3, we monitor the
Hamiltonian constraint $\mathcal{C}_{tt}$, cf. Eq. 41, in order to confirm
that we are indeed evolving a well-posed IVP. Since the other constraints (and
the evolution equations) are coupled, it can be expected that violations of
any other constraint will percolate into the Hamiltonian constraint.

###  IV.1 Numerical method

We use a fourth-order finite difference method to evaluate spatial derivatives
and a fourth-order Runge-Kutta method to evolve in time. The computational
domain is chosen as $x\in(0,10M]$, working in units of the mass $M$. We evolve
all our equations in a unigrid with $N_{x}=1025$ points. Hence the grid
resolution is $\Delta x\simeq 0.01\,M$ with the Courant–Friedrichs–Lewy
condition Courant et al. (1967) set to 0.25.

We perform a robust stability test in order to test whether the evolution is
well-posed. Given the respective initial data
$(\mathbf{u}_{0},\mathbf{v}_{0},\dot{\mathbf{u}}_{0},\dot{\mathbf{v}}_{0},\ddot{\mathbf{u}}_{0})$
for flat (Sec. IV.2) and Schwarzschild (Sec. IV.3) spacetime, we add random
noise to all components of $\mathbf{u}$ and $\mathbf{v}$, i.e,

$\displaystyle\mathbf{u}^{i}(t=0)$
$\displaystyle=\mathbf{u}^{i}_{0}+A_{\textrm{noise}}\textrm{RAND}(x)\;\;\forall\;i=1,\dots,5\;,$
$\displaystyle\mathbf{v}^{i}(t=0)$
$\displaystyle=\mathbf{v}^{i}_{0}+A_{\textrm{noise}}\textrm{RAND}(x)\;\;\forall\;i=1,\dots,3\;.$
(43)

Here, $A_{\textrm{noise}}$ is a noise amplitude which we vary from $10^{-10}$
to $10^{-5}$ and $\textrm{RAND}(x)$ generates random values between -1 and 1.
We do not add random noise to $\dot{\mathbf{u}}_{0}$, $\dot{\mathbf{v}}_{0}$,
and $\ddot{\mathbf{u}}_{0}$. As the random noise violates the constraints, the
above constitutes a robust stability test for whether the Hamiltonian
constraint exhibits growth modes that would signal an ill-posed IVP.

###  IV.2 Well-posed evolution for flat spacetime

Figure 1:  Constraint plot ($l2$-norm of the Hamiltonian constraint in Eq.
(41)) for flat-space initial data, performing a noise (robust stability) test
with different noise amplitudes, ranging from $10^{-5}$ to $10^{-10}$, top to
bottom. Each curve represents an increase by a factor of ten in the initial
amplitude over the curve below.

In order to investigate numerical stability (i.e., well-posedness) of the
evolution in flat spacetime, we initialize the evolution at

$\displaystyle\mathbf{u}_{0}$
$\displaystyle=(\mathcal{R},\,g_{tt},\,g_{tx},\,g_{xx},\,g_{yy})=(0,-1,0,1,1)\;,$
$\displaystyle\mathbf{v}_{0}$
$\displaystyle=\dot{\mathbf{v}}_{0}=\dot{\mathbf{u}}_{0}=\ddot{\mathbf{u}}_{0}=0\;,$
(44)

adding random noise as specified in Eq. (43). We evolve all tests until
$t/M\simeq 1000$ which provides a reasonably long evolution to determine
whether any growth modes are developed.

We monitor the behavior of the Hamiltonian constraint, cf. Eq. (41), during
evolution, cf. Fig. 1. This clearly shows that the constraint decays and
stabilizes. There appear to be no growth modes in the constraint for flat-
space initial data, at least during the monitored time.

We conclude that the performed noise test finds no indications of numerically
unstable behavior which strongly suggests that the evolution equations with
flat-space initial data constitute a well-posed IVP.

###  IV.3 Well-posed evolution for Schwarzschild spacetime

Figure 2:  Constraint plot ($l2$-norm of the Hamiltonian constraint in Eq.
(41)) for Schwarzschild initial data with added noise. Different noise
amplitudes were chosen as described in Fig. 1

Going beyond flat spacetime, we also consider initial data for a Schwarzschild
black hole in Cartesian Kerr-Schild coordinates, i.e.,

$\displaystyle g_{tt,0}$ $\displaystyle=\sqrt{\frac{r}{r+2M}}\;,\quad
g_{tx,0}=\frac{2M}{r}\frac{x}{r+2M}\;,$ $\displaystyle g_{xx,0}$
$\displaystyle=1+\frac{2Mx^{2}}{r^{3}}\;,\quad
g_{yy,0}=1+\frac{2My^{2}}{r^{3}}\;,$ $\displaystyle\mathcal{R}_{0}$
$\displaystyle=0\;,$ $\displaystyle\mathbf{v}_{0}$
$\displaystyle=\dot{\mathbf{v}}_{0}=\dot{\mathbf{u}}_{0}=\ddot{\mathbf{u}}_{0}=0\;,$
(45)

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$ and $M$ is the mass of the Schwarzschild
black hole. We evolve simulations until $t/M\simeq 200$ which suffices for the
constraints to settle into near-stable state if no growth modes develop.

The noise test, cf. Fig. 2, confirms that, also for Schwarzschild initial
data, the evolution constitutes a well-posed IVP, cf. Fig. 2.

Figure 3:  Constraint plot for Schwarzschild initial data for different QG
mass parameters $m_{0}$ and $m_{2}$ (in units of the Schwarzschild mass $M$).

In addition, we vary the masses $m_{0}$ and $m_{2}$ of the QG spin-0 and
spin-2 modes, cf. Eq. (3), in order to confirm that the well-posed constraint
evolution persists for a range of values for $m_{0}$ and $m_{2}$ in the
vicinity of the Schwarzschild mass $M$, cf. Fig. 3. We observe _numerically_
stable evolution for all tested values of $m_{0}$ and $m_{2}$. As for the
noise test, the constraint decays and stabilizes over the investigated time.

We emphasize that this does not test physical stability. Even in mass-ranges
(for the QG masses $m_{0}$ and $m_{2}$) in which the Schwarzschild solution
could become _physically_ unstable, we expect the constraint evolution to be
_numerically_ stable, i.e., well-posed. In particular, even during a potential
physical decay of Schwarzschild spacetime, possibly to some other solution of
QG Pravda et al. (2017); Podolský et al. (2020), the constraints should be
preserved. We will investigate such physical (in)stability in future work.

Overall, we conclude that, also for Schwarzschild initial data, no indications
of numerically unstable behavior are found which strongly suggests that the
evolution is well-posed.

##  V Conclusion

As a proof of principle for well-posed numerical evolution in beyond-GR
theories with higher-curvature operators, we successfully obtain the first
fully nonlinear time evolution in Quadratic Gravity, i.e., for a gravitational
theory including the lowest-order EFT corrections to GR. The system of
evolution equations is obtained by (i) use of harmonic gauge to treat the
Ricci scalar and traceless Ricci tensor as independent variables, cf. Noakes
(1983), (ii) reduction to spherical symmetry, which preserves harmonic gauge
by use of the Cartoon method Alcubierre et al. (2001); Pretorius (2005a), and
(iii) order-reduction to a set of manifestly 1st-order (in time) evolution
equations.

We perform fully non-linear numerical simulations that confirm that there
exists physically significant initial data for which the system of equations
is well-posed. In particular, we observe numerically stable dynamics in the
spherically-symmetric sector for random perturbations about flat spacetime and
about the Schwarzschild solution. (Being Ricci-flat, the latter is also a
solution of QG.)

This opens up several opportunities for future work. As a direct application,
we will investigate physically (un)stable branches of Schwarzschild BHs and
other exotic BHs in QG, cf. Lu et al. (2015); Lü et al. (2015, 2017); Holdom
and Ren (2017); Kokkotas et al. (2017); Bonanno and Silveravalle (2019);
Podolský et al. (2020), to determine the final state of spherical
gravitational collapse.

As a long-term goal, the present study strongly motivates that one can also
establish a well-posed IVP in full 3+1. We are currently working on a well-
posed formulation within a (3+1)-dimensional computational infrastructures
such as Fernando et al. (2019); Haas and et al. (2020). We emphasize that
while the proof in Noakes (1983) guarantees the existence of a well-posed IVP
also in (3+1) dimensions, the construction of such a formulation remains non-
trivial. Establishing such a well-posed IVP will eventually enable us to
perform binary-black-hole mergers to extract GW signals in this theory.

Finally, the presented methodology is, in principle, applicable also to other
theories with higher-order equations of motion. As long as the order-reduced
equations of motion are amenable to diagonalization to quasi-linear form,
other higher-derivative theories may also admit a well-posed IVP.

#### Acknowledgements.

We are grateful for hospitality at Perimeter Institute for Theoretical
Physics, where the original idea for this manuscript was developed. Research
at Perimeter Institute is supported in part by the Government of Canada
through the Department of Innovation, Science and Economic Development Canada
and by the Province of Ontario through the Ministry of Colleges and
Universities. We thank Frans Pretorius and Justin Ripley for correspondence
regarding the Cartoon method as well as Eric W. Hirschmann and Astrid Eichhorn
for discussions. HL is supported by the LANL ASC Program and LDRD grants
20190021DR. This work used resources provided by the LANL Institutional
Computing Program. LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S.DOE (Contract No.
89233218CNA000001). AH is supported by a Royal Society Newton International
fellowship [NIF/R1/191008]. This work is authorized for unlimited release
under LA-UR-21-22739

##  Appendix A Explicit form of lower-order terms

The explicit form of the lower-order terms in the metric Eq. (9) is given by

$\displaystyle\mathcal{O}^{1}_{ab}(\partial g)=$
$\displaystyle\;F^{c}\Gamma_{(ab)c}+2g^{ed}\Gamma^{c}_{e(a}\Gamma_{b)cd}+g^{cd}\Gamma^{e}_{ad}\Gamma_{ecb}\;,$
(46)

where the first term vanishes upon use of harmonic coordinates. This comes
about from the well-known expansion of the Ricci-tensor in harmonic
coordinates.

The explicit form of
$\mathcal{O}^{2}_{ab}(\partial\mathcal{R},\partial\widetilde{\mathcal{R}},\partial\partial
g)$ in the traceless equation, cf. Eq. (11), is given by

$\displaystyle\mathcal{O}^{2}_{ab}(\partial\partial\mathcal{R},\partial\widetilde{\mathcal{R}},\partial\partial
g)=$
$\displaystyle-\frac{1}{3}\left(\frac{m_{2}^{2}}{m_{0}^{2}}-1\right)\mathcal{R}_{,ab}\underline{-\widetilde{\mathcal{R}}_{cd}g^{cf}g^{ed}\left(g_{ef,ab}+g_{ab,ef}-g_{ae,fb}-g_{fb,ae}\right)}$
$\displaystyle\dashuline{+\frac{1}{6}\left(\widetilde{\mathcal{R}}_{ab}g^{ef}g^{cd}g_{ef,cd}\right)+g^{ef}g^{cd}\widetilde{\mathcal{R}}_{e(a}\left(g_{df,b)c}+g_{b)f,cd}-g_{b)d,cf}\right)}$
$\displaystyle\underline{+2\widetilde{\mathcal{R}}_{cd}\,g^{ce}\Big{[}2\Gamma^{d}_{f[e}\,\Gamma^{f}_{a]b}+g^{df}\,g^{mn}g_{fn,[a}\left(g_{e]m,b}+g_{bm,e]}-g_{e]b,m}\right)\Big{]}}$
$\displaystyle\dashuline{+g^{cd}\Big{[}2\Gamma^{n}_{c(a}\left(\Gamma^{m}_{dn}\widetilde{\mathcal{R}}_{b)m}+\widetilde{\mathcal{R}}_{b)n,d}\right)+2\Gamma^{m}_{d(b}\left(\Gamma^{n}_{a)c}\widetilde{\mathcal{R}}_{mn}-\widetilde{\mathcal{R}}_{a)m,c}\right)}$
$\displaystyle\quad\quad\quad\dashuline{+g^{ef}g^{mn}g_{fn,c}\widetilde{\mathcal{R}}_{e(b}\left(g_{dm,a)}+g_{a)m,d}-g_{a)d,m}\right)\Big{]}}$
$\displaystyle+m_{2}^{2}\widetilde{\mathcal{R}}_{ab}+2T^{\text{(TL)}}_{ab}+\frac{1}{12}\left(\frac{m_{2}^{2}}{m_{0}^{2}}-1\right)g_{ab}m_{0}^{2}\mathcal{R}$
$\displaystyle-\frac{1}{3}\left(\frac{m_{2}^{2}}{m_{0}^{2}}+1\right)\mathcal{R}\widetilde{\mathcal{R}}_{ab}-2g^{cd}\widetilde{\mathcal{R}}_{ac}\widetilde{\mathcal{R}}_{bd}+\frac{1}{2}g_{ab}\widetilde{\mathcal{R}}^{cd}\widetilde{\mathcal{R}}_{cd}\;.$
(47)

We have sorted the contributions in terms of their order in time derivatives:
the first two lines collect all 2nd-order terms; lines three to five all 1st-
order terms; the last two lines collect the 0th-order contributions.

Furthermore, we have underlined terms arising from
${\widetilde{\mathcal{R}}^{cd}R_{abcd}}$ (continuous) and from
${\Box\widetilde{\mathcal{R}}_{ab}}$ (dashed). The former agree with Noakes
(1983) (apart from minor sign typos) but the latter do not. In any case, as
discussed in the main text, a modification of these terms (as long as their
derivative order is preserved) does not impact Noakes’ proof of well-
posedness.

##  Appendix B Convergence Tests

We perform a standard convergence test to confirm the validity of our
implementation. We choose three different numerical grid sizes $N_{x}=1025$,
$N_{x}=2049$, and $N_{x}=4097$ as low, medium, and high resolutions,
respectively. We use the Schwarzschild black hole which is described in Sec IV
as initial data.

Figure 4:  Convergence test with Schwarzschild BH. Mass of BH is chosen as one
and it was evolved until $t/M\simeq 200$. We choose $m_{0}=0.1,m_{2}=0.1$ as
QG mass parameters. We plot differences of $\mathcal{C}_{tt}$ with different
resolutions. Differences remain small during the evolution. We rescale the
smaller resolution difference by 16, which is the expected convergence rate of
our implementation since we use fourth-order finite difference stencils.

Fig. 4 summarizes the result of this convergence test. We plot differences of
the constraint value, $\mathcal{C}_{tt}$, with different resolutions as a
function of time. We re-scale the smaller resolution difference by a factor of
$16$, which corresponds to the expected convergence rate given the applied
fourth-order finite difference stencils that we use. The differences remain
small and both lines agree which confirms the fourth-order convergence.

## References

  * Penrose (1965) R. Penrose, Phys. Rev. Lett. 14, 57 (1965).
  * Hawking (1965) S. Hawking, Phys. Rev. Lett. 15, 689 (1965).
  * Geroch (1966) R. P. Geroch, Phys. Rev. Lett. 17, 445 (1966).
  * Hawking and Penrose (1970) S. W. Hawking and R. Penrose, Proc. Roy. Soc. Lond. A314, 529 (1970).
  * Dafermos (2005) M. Dafermos, Commun. Pure Appl. Math. 58, 0445 (2005), eprint gr-qc/0307013.
  * Harada et al. (2002) T. Harada, H. Iguchi, and K.-i. Nakao, Prog. Theor. Phys. 107, 449 (2002), eprint gr-qc/0204008.
  * ’t Hooft and Veltman (1974) G. ’t Hooft and M. J. G. Veltman, Ann. Inst. H. Poincare Phys. Theor. A20, 69 (1974).
  * Stelle (1977) K. S. Stelle, Phys. Rev. D 16, 953 (1977).
  * Goroff and Sagnotti (1985) M. H. Goroff and A. Sagnotti, Phys. Lett. 160B, 81 (1985).
  * Avramidi and Barvinsky (1985) I. G. Avramidi and A. O. Barvinsky, Phys. Lett. B 159, 269 (1985).
  * van de Ven (1992) A. E. M. van de Ven, Nucl. Phys. B378, 309 (1992).
  * Boulware and Deser (1985) D. G. Boulware and S. Deser, Phys. Rev. Lett. 55, 2656 (1985).
  * Zwiebach (1985) B. Zwiebach, Phys. Lett. 156B, 315 (1985).
  * Percacci (2017) R. Percacci, _An Introduction to Covariant Quantum Gravity and Asymptotic Safety_ , vol. 3 of _100 Years of General Relativity_ (World Scientific, 2017), ISBN 9789813207172, 9789813207196, 9789813207172, 9789813207196.
  * Eichhorn (2019) A. Eichhorn, Frontiers in Astronomy and Space Sciences 5, 47 (2019), eprint 1810.07615.
  * Reuter and Saueressig (2019) M. Reuter and F. Saueressig, _Quantum Gravity and the Functional Renormalization Group_ (Cambridge University Press, 2019), ISBN 9781107107328.
  * Penrose (1969) R. Penrose, Riv. Nuovo Cim. 1, 252 (1969).
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], eprint 1602.03841.
  * Abbott et al. (2020) R. Abbott et al. (LIGO Scientific, Virgo) (2020), eprint 2010.14527, URL https://dcc.ligo.org/LIGO-P2000223/public.
  * Shibata and Nakamura (1995) M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
  * Baumgarte and Shapiro (1998) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 024007 (1998), eprint gr-qc/9810065.
  * Pretorius (2005a) F. Pretorius, Classical and Quantum Gravity 22, 425 (2005a), URL https://doi.org/10.1088/0264-9381/22/2/014.
  * Pretorius (2005b) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005b), eprint gr-qc/0507014.
  * Sarbach and Tiglio (2012) O. Sarbach and M. Tiglio, Living Reviews in Relativity 15, 9 (2012), ISSN 1433-8351, URL https://doi.org/10.12942/lrr-2012-9.
  * Isenberg (2014) J. Isenberg, _The Initial Value Problem in General Relativity_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2014), pp. 303–321, ISBN 978-3-642-41992-8, URL https://doi.org/10.1007/978-3-642-41992-8_16.
  * Okounkova et al. (2019) M. Okounkova, L. C. Stein, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 100, 104026 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.104026.
  * Okounkova et al. (2020) M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 101, 104016 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.104016.
  * Ripley and Pretorius (2020a) J. L. Ripley and F. Pretorius, Classical and Quantum Gravity 37, 155003 (2020a), URL https://doi.org/10.1088/1361-6382/ab9bbb.
  * Ripley and Pretorius (2020b) J. L. Ripley and F. Pretorius, Phys. Rev. D 101, 044015 (2020b), URL https://link.aps.org/doi/10.1103/PhysRevD.101.044015.
  * Okounkova (2020) M. Okounkova, Phys. Rev. D 102, 084046 (2020), eprint 2001.03571.
  * East and Ripley (2021) W. E. East and J. L. Ripley, Phys. Rev. D 103, 044040 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.044040.
  * Silva et al. (2020) H. O. Silva, H. Witek, M. Elley, and N. Yunes (2020), eprint 2012.10436.
  * Figueras and França (2020) P. Figueras and T. França, Class. Quant. Grav. 37, 225009 (2020), eprint 2006.09414.
  * Bezares et al. (2021) M. Bezares, M. Crisostomi, C. Palenzuela, and E. Barausse, JCAP 03, 072 (2021), eprint 2008.07546.
  * Cayuso and Lehner (2020) R. Cayuso and L. Lehner, Phys. Rev. D 102, 084008 (2020), eprint 2005.13720.
  * Kovács (2019) A. D. Kovács, Phys. Rev. D 100, 024005 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.024005.
  * Kovács and Reall (2020a) A. D. Kovács and H. S. Reall, Phys. Rev. Lett. 124, 221101 (2020a), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.221101.
  * Kovács and Reall (2020b) A. D. Kovács and H. S. Reall, Phys. Rev. D 101, 124003 (2020b), URL https://link.aps.org/doi/10.1103/PhysRevD.101.124003.
  * Stelle (1978) K. S. Stelle, General Relativity and Gravitation 9, 353 (1978), ISSN 1572-9532, URL https://doi.org/10.1007/BF00760427.
  * Noakes (1983) D. R. Noakes, Journal of Mathematical Physics 24, 1846 (1983).
  * Morales and Santillán (2019) J. O. Morales and O. P. Santillán, JCAP 03, 026 (2019), eprint 1811.07869.
  * Ostrogradsky (1850) M. Ostrogradsky, Mem. Ac. St. Petersbourg VI 4, 385 (1850).
  * Becker et al. (2017) D. Becker, C. Ripken, and F. Saueressig, JHEP 12, 121 (2017), eprint 1709.09098.
  * Anselmi and Piva (2018) D. Anselmi and M. Piva, JHEP 05, 027 (2018), eprint 1803.07777.
  * Donoghue and Menezes (2019) J. F. Donoghue and G. Menezes, Phys. Rev. D 100, 105006 (2019), eprint 1908.02416.
  * Salvio (2019) A. Salvio, Phys. Rev. D 99, 103507 (2019), eprint 1902.09557.
  * Alcubierre et al. (2001) M. Alcubierre, S. Brandt, B. Bruegmann, D. Holz, E. Seidel, R. Takahashi, and J. Thornburg, Int. J. Mod. Phys. D 10, 273 (2001).
  * Leray and Institute for Advanced Study (Princeton(1953) J. Leray and N. Institute for Advanced Study (Princeton, _Hyperbolic differential equations_ (Princeton Institute for Advanced Study, 1953), URL https://books.google.com/books?id=q_w-AAAAIAAJ.
  * Lim and Held (2021) H. Lim and A. Held (2021), to appear.
  * Sorkin and Choptuik (2010) E. Sorkin and M. W. Choptuik, Gen. Rel. Grav. 42, 1239 (2010), eprint 0908.2500.
  * Gundlach and Martin-Garcia (2006) C. Gundlach and J. M. Martin-Garcia, Class. Quant. Grav. 23, S387 (2006), eprint gr-qc/0506037.
  * Richter and Hilditch (2011) R. Richter and D. Hilditch, J. Phys. Conf. Ser. 314, 012102 (2011), eprint 1303.4433.
  * Hilditch and Richter (2015) D. Hilditch and R. Richter, J. Hyperbol. Diff. Equat. 12, 1 (2015), eprint 1412.6034.
  * Courant et al. (1967) R. Courant, K. Friedrichs, and H. Lewy, IBM Journal of Research and Development 11, 215 (1967).
  * Pravda et al. (2017) V. Pravda, A. Pravdová, J. Podolský, and R. Švarc, Phys. Rev. D 95, 084025 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.084025.
  * Podolský et al. (2020) J. Podolský, R. Švarc, V. Pravda, and A. Pravdova, Phys. Rev. D 101, 024027 (2020), eprint 1907.00046.
  * Lu et al. (2015) H. Lu, A. Perkins, C. N. Pope, and K. S. Stelle, Phys. Rev. Lett. 114, 171601 (2015), eprint 1502.01028.
  * Lü et al. (2015) H. Lü, A. Perkins, C. N. Pope, and K. S. Stelle, Phys. Rev. D 92, 124019 (2015), eprint 1508.00010.
  * Lü et al. (2017) H. Lü, A. Perkins, C. N. Pope, and K. S. Stelle, Phys. Rev. D 96, 046006 (2017), eprint 1704.05493.
  * Holdom and Ren (2017) B. Holdom and J. Ren, Phys. Rev. D 95, 084034 (2017), eprint 1612.04889.
  * Kokkotas et al. (2017) K. Kokkotas, R. A. Konoplya, and A. Zhidenko, Phys. Rev. D 96, 064007 (2017), eprint 1705.09875.
  * Bonanno and Silveravalle (2019) A. Bonanno and S. Silveravalle, Phys. Rev. D 99, 101501 (2019), eprint 1903.08759.
  * Fernando et al. (2019) M. Fernando, D. Neilsen, H. Lim, E. Hirschmann, and H. Sundar, SIAM J. Sci. Comput. 41, C97 (2019).
  * Haas and et al. (2020) R. Haas and et al., _The einstein toolkit_ (2020), to find out more, visit http://einsteintoolkit.org, URL https://doi.org/10.5281/zenodo.4298887.
  * (65) W. R. Inc., _Mathematica, Version 12.2_ , champaign, IL, 2020, URL https://www.wolfram.com/mathematica.
  * Martin-Garcia et al. (2007) J. M. Martin-Garcia, R. Portugal, and L. R. U. Manssur, Comput. Phys. Commun. 177, 640 (2007), eprint 0704.1756.