Well-posedness of cubic Horndeski theories

TL;DR

This paper investigates the local well-posedness of cubic Horndeski theories by extending three hyperbolic formulations of Einstein equations, aiming to facilitate numerical simulations in modified gravity models.

Contribution

It introduces three strongly hyperbolic formulations of the Einstein equations adapted for cubic Horndeski theories, extending existing methods to this class of modified gravity.

Findings

Extended elliptic-hyperbolic system formulation.

Adapted BSSN formulation with generalized slicing.

Extended CCZ4 formulation with shift conditions.

Abstract

We study the local well-posedness of the initial value problem for cubic Horndeski theories. Three different strongly hyperbolic modifications of the ADM formulation of the Einstein equations are extended to cubic Horndeski theories in the weak field regime. In the first one, the equations of motion are rewritten as a coupled elliptic-hyperbolic system of partial differential equations. The second one is based on the BSSN formulation with a generalised Bona-Mass\'o slicing (covering the 1+log slicing) and non-dynamical shift vector. The third one is an extension of the CCZ4 formulation with a generalised Bona-Mass\'o slicing (also covering the 1+log slicing) and a gamma driver shift condition. This final formulation may be particularly suitable for applications in non-linear numerical simulations.

Tables

Table 1. Table 1 : The list of variables used in the BSSN formulation

Quantity/Definition	Notation	Linearized version
Conformal factor	$\Omega$	$\omega$
Conformal metric	${\tilde{h}}_{ij}$	${\tilde{\gamma }}_{ij}$
Lapse function	$N$	$\alpha$
Scalar field	$\varphi$	$\psi$
$A\equiv {ℒ}_n\varphi$	$A$	$a$
Trace of the extrinsic curvature	$K$	$\kappa$
Conformal traceless extrinsic curvature	${\tilde{Q}}_{ij}$	${\tilde{q}}_{ij}$
${\tilde{V}}^i\equiv -{\mathring{D}}_j{\tilde{h}}^{ij}$	${\tilde{V}}^i$	${\tilde{v}}^i$

Table 2. Table 2 : The list of eigenvalues and eigenvectors of the principal symbol with ψ ^ = a ^ = 0 ^ 𝜓 ^ 𝑎 0 \hat{\psi}=\hat{a}=0 .

	$\hat{\alpha }$	$\hat{\omega }$	$\hat{\kappa }$	${\hat{\tilde{\gamma }}}_{ij}$	${\hat{\tilde{v}}}_i$	${\overline{\xi }}_0$
I.	$0$	$0$	$0$	${\tilde{\gamma }}_{ij}^{TT}$	$0$	$\pm 1$
II.	$0$	$0$	$0$	${\tilde{h}}_{ij}-\frac{{\xi }_i{\xi }_j}{{\left|\xi \right|}_{\tilde{h}}^2}$	$0$	$\pm 1$
III.	$0$	$0$	$0$	${\xi }_{(i}{e}_{j)}$	$-i\frac{m-1}{2}{\left|\xi \right|}_{\tilde{h}}^2{e}_i$	$\pm \sqrt{m}$
IV.	$0$	$\hat{\omega }$	$0$	$-\frac{2}{{\left|\xi \right|}_{\tilde{h}}^2}{\left(i{\xi }_{(i}{\hat{\tilde{v}}}_{j)}+2{\xi }_i{\xi }_j\hat{\omega }\right)}^{\text{TF}}$	${\hat{\tilde{v}}}_i$	$0$
V.	$N\sigma$	$\frac{1}{12}$	$-\frac{i{\xi }_0}{2N}$	${\left(\frac{{\xi }_i{\xi }_j}{{\left|\xi \right|}_{\tilde{h}}^2}\right)}^{\text{TF}}$	$i\frac{2}{3}{\xi }_i$	$\pm \sqrt{2\sigma }$
VI.	$0$	$0$	$0$	${\left(\frac{{\xi }_i{\xi }_j}{{\left|\xi \right|}_{\tilde{h}}^2}\right)}^{\text{TF}}$	$2i\left(m-1\right){\xi }_i$	$\pm \sqrt{\frac{4m-1}{3}}$

Table 3. (a)

	$\hat{\alpha }$	$\hat{\omega }$	${\hat{\beta }}^i$	${\hat{\tilde{\gamma }}}_{ij}$	${\overline{\xi }}_0$
I.	$0$	$0$	$0$	${\hat{\tilde{\gamma }}}_{ij}$	$\pm 1$
II.	$\mp i\frac{2\sigma \left(4f-3\right){\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{f\left(6\sigma +1\right)}$	$\mp i\frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }\left(f-1-2\sigma f\right)}{2fN\left(6\sigma +1\right)}$	$i{\xi }^i$	${\hat{\tilde{\gamma }}}_{ij}$	$\pm 1$
III.	$0$	$-\frac{1}{8{\left|\xi \right|}_{\tilde{h}}^2}i{\xi }_k{\hat{\tilde{v}}}^k$	$0$	$-\frac{2i}{{\left|\xi \right|}_{\tilde{h}}^2}{\left({\xi }_{(i}{v}_{j)}-\frac{{\xi }_i{\xi }_j}{4{\left|\xi \right|}_{\tilde{h}}^2}{\xi }_k{\hat{\tilde{v}}}^k\right)}^{\text{TF}}$	$0$
IV.	$\pm i\frac{\sqrt{2\sigma }{e}^{2\Omega }{\left|\xi \right|}_{\tilde{h}}\left(3\sigma -2f\right)}{2f}$	$\pm i\frac{\sqrt{2\sigma }{e}^{2\Omega }{\left|\xi \right|}_{\tilde{h}}}{8Nf}$	$i{\xi }^i$	$\pm i\frac{3\sqrt{2\sigma }{e}^{2\Omega }}{2fN{\left|\xi \right|}_{\tilde{h}}}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\pm \sqrt{2\sigma }$
V.	$0$	$0$	$i{e}^i$	$\pm i\frac{2e^{2\Omega }}{N\sqrt{f}{\left|\xi \right|}_{\tilde{h}}}{\xi }_{(i}{e}_{j)}$	$\pm \sqrt{f}$
VI.	$0$	$\pm i\frac{e^{2\Omega }{\left|\xi \right|}_{\tilde{h}}}{4\sqrt{3f}N}$	$i{\xi }^i$	$\pm i\frac{\sqrt{3}{e}^{2\Omega }}{N\sqrt{f}{\left|\xi \right|}_{\tilde{h}}}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\pm \sqrt{\frac{4f}{3}}$

Table 4. (a)

	$\hat{\alpha }$	$\hat{\omega }$	${\hat{\beta }}^i$	${\hat{\tilde{\gamma }}}_{ij}$	${\overline{\xi }}_0$
I.	$0$	$0$	$0$	${\hat{\tilde{\gamma }}}_{ij}$	$\pm 1$
II.	$\mp i\frac{2\sigma \left(4f-3\right){\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{f\left(6\sigma +1\right)}$	$\mp i\frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }\left(f-1-2\sigma f\right)}{2fN\left(6\sigma +1\right)}$	$i{\xi }^i$	${\hat{\tilde{\gamma }}}_{ij}$	$\pm 1$
III.	$0$	$-\frac{1}{8{\left|\xi \right|}_{\tilde{h}}^2}i{\xi }_k{\hat{\tilde{v}}}^k$	$0$	$-\frac{2i}{{\left|\xi \right|}_{\tilde{h}}^2}{\left({\xi }_{(i}{v}_{j)}-\frac{{\xi }_i{\xi }_j}{4{\left|\xi \right|}_{\tilde{h}}^2}{\xi }_k{\hat{\tilde{v}}}^k\right)}^{\text{TF}}$	$0$
IV.	$\pm i\frac{\sqrt{2\sigma }{e}^{2\Omega }{\left|\xi \right|}_{\tilde{h}}\left(3\sigma -2f\right)}{2f}$	$\pm i\frac{\sqrt{2\sigma }{e}^{2\Omega }{\left|\xi \right|}_{\tilde{h}}}{8Nf}$	$i{\xi }^i$	$\pm i\frac{3\sqrt{2\sigma }{e}^{2\Omega }}{2fN{\left|\xi \right|}_{\tilde{h}}}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\pm \sqrt{2\sigma }$
V.	$0$	$0$	$i{e}^i$	$\pm i\frac{2e^{2\Omega }}{N\sqrt{f}{\left|\xi \right|}_{\tilde{h}}}{\xi }_{(i}{e}_{j)}$	$\pm \sqrt{f}$
VI.	$0$	$\pm i\frac{e^{2\Omega }{\left|\xi \right|}_{\tilde{h}}}{4\sqrt{3f}N}$	$i{\xi }^i$	$\pm i\frac{\sqrt{3}{e}^{2\Omega }}{N\sqrt{f}{\left|\xi \right|}_{\tilde{h}}}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\pm \sqrt{\frac{4f}{3}}$

Table 5. (b)

	$\hat{\theta }$	$\hat{\kappa }$	${\hat{b}}^i$	${\hat{\tilde{q}}}_{ij}$	${\hat{\tilde{u}}}_i$
I.	$0$	$0$	$0$	$\mp \frac{1}{2}{\left|\xi \right|}_{\tilde{h}}{e}^{-2\Omega }{\hat{\tilde{\gamma }}}_{ij}$	$0$
II.	$\frac{{\left|\xi \right|}_{\tilde{h}}^2\left(4f-3\right)\left(2\sigma -1\right)}{2fN\left(6\sigma +1\right)}$	$\frac{3-4f}{Nf\left(6\sigma +1\right)}{\left|\xi \right|}_{\tilde{h}}^2$	$\mp \frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{Nf}{\xi }^i$	$\mp \frac{1}{2}{\left|\xi \right|}_{\tilde{h}}{e}^{-2\Omega }{\hat{\tilde{\gamma }}}_{ij}-\frac{1}{N}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\mp \frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{Nf}{\xi }_i$
III.	$0$	$0$	$0$	$0$	${\hat{\tilde{v}}}_i$
IV.	$0$	$\frac{3\sigma -2f}{2fN}{\left|\xi \right|}_{\tilde{h}}^2$	$\mp \sqrt{2\sigma }\frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{Nf}{\xi }^i$	$\pm i\frac{3\sigma -2f}{2fN{\left|\xi \right|}_{\tilde{h}}}{\left({\xi }_i{\xi }_j\right)}^{\text{TF}}$	$\mp \sqrt{2\sigma }\frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{Nf}{\xi }_i$
V.	$0$	$0$	$\mp \frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{N\sqrt{f}}{e}^i$	$0$	$\mp \frac{{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{N\sqrt{f}}{e}_i$
VI.	$0$	$0$	$\mp \frac{2{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{N\sqrt{3f}}{\xi }^i$	$0$	$\mp \frac{2{\left|\xi \right|}_{\tilde{h}}{e}^{2\Omega }}{N\sqrt{3f}}{\xi }_i$

Equations

$$S=\frac{1}{16\pi G}\int {\mathrm{d}}^{4}x\sqrt{-g}\left({\mathcal{L}}_1+{\mathcal{L}}_2+{\mathcal{L}}_3+{\mathcal{L}}_4+{\mathcal{L}}_5\right)$$

$${\mathcal{L}}_1$$

$${\mathcal{L}}_2$$

$${\mathcal{L}}_3$$

$${\mathcal{L}}_4$$

$${\mathcal{L}}_5$$

$$S=\frac{1}{16\pi G}\int {\mathrm{d}}^{4}x\sqrt{-g}\left(\mathcal{R}+X+G_2(\varphi ,X)+G_3(\varphi ,X)\square \varphi \right).$$

$$E_{ab}\equiv$$

$$+\frac{1}{2}{\partial }_X{G}_3\left(-\square \varphi {\nabla }_a\varphi {\nabla }_b\varphi +2{\nabla }_{(a}\varphi {\nabla }_{b)}{\nabla }_c\varphi {\nabla }^c\varphi -{\nabla }_c{\nabla }_d\varphi {\nabla }^c\varphi {\nabla }^d\varphi g_{ab}\right)=0.$$

$$E_{ab}-\frac{1}{2}E{g}_{ab}=$$

$$+\frac{1}{2}{\partial }_X{G}_3\left(-\square \varphi {\nabla }_a\varphi {\nabla }_b\varphi +2{\nabla }_{(a}\varphi {\nabla }_{b)}{\nabla }_c\varphi {\nabla }^c\varphi \right)=0.$$

$$E_{\varphi }\equiv$$

$$-\left({\partial }_X^2{G}_2+2{\partial }_{X\varphi }^2{G}_3\right)\left((\partial \varphi {)}^2\square \varphi -{\nabla }^a\varphi {\nabla }^b\varphi {\nabla }_a{\nabla }_b\varphi \right)-{\partial }_X{G}_3\left((\square \varphi {)}^2-{\nabla }_a{\nabla }_b\varphi {\nabla }^a{\nabla }^b\varphi \right)$$

$$+{\partial }_X^2{G}_3{\nabla }^a\varphi {\nabla }^b\varphi \left(\square \varphi {\nabla }_a{\nabla }_b\varphi -{\nabla }_a{\nabla }^c\varphi {\nabla }_c{\nabla }_b\varphi \right)+2X\left({\partial }_{\varphi }^2{G}_3+{\partial }_{X\varphi }^2{G}_2\right)-{\partial }_{\varphi }{G}_2=0$$

$${\tilde{E}}_{\varphi }\equiv$$

$$-\left({\partial }_X^2{G}_2+2{\partial }_{X\varphi }^2{G}_3\right)\left((\partial \varphi {)}^2\square \varphi -{\nabla }^a\varphi {\nabla }^b\varphi {\nabla }_a{\nabla }_b\varphi \right)-{\partial }_X{G}_3\left((\square \varphi {)}^2-{\nabla }_a{\nabla }_b\varphi {\nabla }^a{\nabla }^b\varphi \right)$$

$$+{\partial }_X^2{G}_3{\nabla }^a\varphi {\nabla }^b\varphi \left(\square \varphi {\nabla }_a{\nabla }_b\varphi -{\nabla }_a{\nabla }^c\varphi {\nabla }_c{\nabla }_b\varphi \right)+2X\left({\partial }_{\varphi }^2{G}_3+{\partial }_{X\varphi }^2{G}_2\right)-{\partial }_{\varphi }{G}_2$$

$$+{\partial }_X{G}_3\left(2X^2+X{G}_2+X^2{\partial }_X{G}_2+X^2{\partial }_X{G}_3\square \varphi +4X^2{\partial }_{\varphi }{G}_3+2X{\partial }_X{G}_3{\nabla }^a\varphi {\nabla }^b\varphi {\nabla }_a{\nabla }_b\varphi \right)=0$$

$${\left(\frac{\partial }{\partial t}\right)}^a=N{n}^a+N^a.$$

$$K_{ab}=-\frac{1}{2}{\mathcal{L}}_n{h}_{ab}=-\frac{1}{2N}\left({\partial }_t-{\mathcal{L}}_N\right)h_{ab}.$$

$$A\equiv n^a{\nabla }_a\varphi =\frac{1}{N}\left({\partial }_t-{\mathcal{L}}_N\right)\varphi$$

$$A_a\equiv h_a^b{\nabla }_b\varphi .$$

$$\displaystyle\left(\partial_{t}-\mathcal{L}_{N}\right)K_{ij}-\frac{1}{2}\partial_{X}G_{3}(Xh_{ij}+A_{i}A_{j})\left(\partial_{t}-\mathcal{L}_{N}\right)A=-D_{i}D_{j}N+N\biggl{\{}R_{ij}+KK_{ij}-2K_{ik}K^{k}_{j}$$

$$+\frac{1}{2}{h}_{ij}\left(G_2-X{\partial }_X{G}_2-X{\partial }_X{G}_3\left(D^k{A}_k+AK+A^k{D}_k\ln N\right)\right)-\frac{1}{2}\left(1+{\partial }_X{G}_2+2{\partial }_{\varphi }{G}_3\right)A_i{A}_j$$

$$\displaystyle+\frac{1}{2}\partial_{X}G_{3}\left(A^{k}D_{k}\left(A_{i}A_{j}\right)-A_{i}A_{j}D^{k}A_{k}-A_{i}A_{j}AK-A_{i}A_{j}A^{k}D_{k}\ln N-2AA_{(i}D_{j)}A\right)\biggr{\}}.$$

$$X=\frac{1}{2}\left(A^2-A^k{A}_k\right).$$

$$2\mathbf{H}\equiv 2E_{\mu \nu }{n}^{\mu }{n}^{\nu }=$$

$$-{\partial }_X{G}_3\left(A^{3}K+A^2{D}^i{A}_i-K_{ij}{A}^i{A}^{j}A-A^i{A}^j{D}_i{A}_j\right)=0,$$

$${\mathbf{M}}_i\equiv E_{\mu i}{n}^{\mu }=$$

$$-\frac{1}{2}{\partial }_X{G}_3\left(A^2{A}_{i}K+A{A}_i{D}^k{A}_k-K_{kl}{A}^k{A}^l{A}_i-A{A}^j{D}_j{A}_i-A_i{A}^k{D}_{k}A+A^2{D}_{i}A\right)=0.$$

$$E_{\varphi }\equiv$$

$$\displaystyle+\tfrac{1}{4}\bigl{(}\partial_{X}G_{3}\bigr{)}^{2}\Bigl{(}3A^{4}-2{A}_{i}{A}_{j}A^{2}{h}^{ij}-{A}_{i}{A}_{j}{A}_{k}{A}_{l}{h}^{ik}{h}^{jl}\Bigr{)}$$

$$\displaystyle+\partial_{{X}}^{2}G_{3}\Bigl{(}A^{2}{\Phi}_{ij}{h}^{ij}-{A}_{k}{A}_{l}{\Phi}_{ij}{h}^{ik}{h}^{jl}\Bigr{)}\biggr{)}$$

$$-{\partial }_{\varphi }{G}_2-{\Phi }_{ij}{h}^{ij}-2{\Phi }_{ij}{\partial }_{\varphi }{G}_3{h}^{ij}+{\partial }_{X\varphi }^2{G}_2(A^2-A_i{A}_j{h}^{ij})$$

$$+{\partial }_{\varphi }^2{G}_3(A^2-A_i{A}_j{h}^{ij})+\frac{1}{2}{G}_2{\partial }_X{G}_3(A^2-A_i{A}_j{h}^{ij})$$

$$+{\partial }_X^2{G}_2(-2A_{i}A{\Phi }_j{h}^{ij}+A_k{A}_l{\Phi }_{ij}{h}^{ik}{h}^{jl})$$

$$\displaystyle+\partial_{X}G_{2}\Bigl{(}-{\Phi}_{ij}{h}^{ij}+\tfrac{1}{4}\partial_{X}G_{3}\bigl{(}A^{4}-2{A}_{i}{A}_{j}A^{2}{h}^{ij}+{A}_{i}{A}_{j}{A}_{k}{A}_{l}{h}^{ik}{h}^{jl}\bigr{)}\Bigr{)}$$

$$\displaystyle+\partial_{X\phi}^{2}G_{3}\Bigl{(}-4{A}_{i}A{\Phi}_{j}{h}^{ij}+A^{2}{\Phi}_{ij}{h}^{ij}+2{A}_{k}{A}_{l}{\Phi}_{ij}{h}^{ik}{h}^{jl}-{A}_{i}{A}_{j}{\Phi}_{kl}{h}^{ij}{h}^{kl}\Bigr{)}$$

$$\displaystyle+\partial_{X}G_{3}\Bigl{(}\partial_{\phi}G_{3}\bigl{(}A^{4}-2{A}_{i}{A}_{j}A^{2}{h}^{ij}+{A}_{i}{A}_{j}{A}_{k}{A}_{l}{h}^{ik}{h}^{jl}\bigr{)}+A^{4}-2{\Phi}_{i}{\Phi}_{j}\ {h}^{ij}$$

$$\displaystyle-2{A}_{i}{A}_{j}A^{2}{h}^{ij}+{A}_{i}{A}_{j}{A}_{k}\ {A}_{l}{h}^{ik}{h}^{jl}-{\Phi}_{ij}{\Phi}_{kl}{h}^{ij}\ {h}^{kl}+{\Phi}_{ij}{\Phi}_{kl}{h}^{ik}\ {h}^{jl}\Bigr{)}$$

$$\displaystyle+\partial_{{X}}^{2}G_{3}\Bigl{(}-A^{2}{\Phi}_{i}{\Phi}_{j}{h}^{ij}-2{A}_{i}A{\Phi}_{jl}{\Phi}_{k}\ {h}^{ik}{h}^{jl}+{A}_{i}{A}_{j}{\Phi}_{k}\ {\Phi}_{l}{h}^{ik}{h}^{jl}$$

$$\displaystyle+2{A}_{k}A{\Phi}_{ij}\ {\Phi}_{l}{h}^{ik}{h}^{jl}-{A}_{l}{A}_{m}{\Phi}_{in}\ {\Phi}_{jk}{h}^{il}{h}^{jm}{h}^{kn}+{A}_{m}{A}_{n}{\Phi}_{ij}\ {\Phi}_{kl}{h}^{im}{h}^{jn}\ {h}^{kl}\Bigr{)}$$

$$\displaystyle+\tfrac{1}{4}\bigl{(}\partial_{X}G_{3}\bigr{)}^{2}\Bigl{(}-8{A}_{i}A^{3}{\Phi}_{j}{h}^{ij}+A^{4}{\Phi}_{ij}{h}^{ij}+4{A}_{k}{A}_{l}A^{2}{\Phi}_{ij}\ {h}^{ik}{h}^{jl}$$

$$-2A_i{A}_k{A}^2{\Phi }_{jl}{h}^{ik}{h}^{jl}+A_i{A}_j{A}_k{A}_l{\Phi }_{mn}{h}^{ik}{h}^{jl}{h}^{mn}$$

Full text

1 Introduction

In the past few decades there has been a growing interest in modifying general relativity (GR) and exploring various properties of these modified theories. The most common reason for this is that the cosmological constant problem led some people to believe that GR may not be the correct theory of gravity, even at large distance scales and low energies. Investigating the properties of extensions of GR is also relevant from an effective field theory (EFT) point of view: since GR is non-renormalizable, one expects that it is only an EFT valid up to some energy scale. When describing strong field phenomena, other operators (besides the Einstein-Hilbert term) might become relevant in the action of the underlying UV complete theory. Finally, testing the rigidity of the predictions and mathematical properties of GR to small deformations can also help us better understand GR itself and discover new techniques.

Trying to modify GR with the purpose to cure one of its shortcomings, however, typically introduces new unwanted pathologies elsewhere. For example, theories with equations of motion containing higher than second derivatives generically suffer from the so-called Ostrogradsky instabilities111However, in the EFT context it may be argued that such a “runaway” behaviour can be discarded as the mass of the ghost is usually around the cutoff and solutions of this type fall beyond the range of validity of the EFT (see e.g. [34]).. To avoid such instabilities, investigations are usually restricted to theories with second order equations of motion. Examples of such theories include Horndeski theories [1]. Horndeski theories are the most general diffeomorphism-invariant scalar-tensor theories with an action principle that has second order equations of motion. The action for general Horndeski theory is

[TABLE]

with

[TABLE]

and $X\equiv-\frac{1}{2}(\partial\phi)^{2}$, $G_{i}$ $(i=2,3,4,5)$ are freely specifiable functions, $\mathcal{R}$ and $\mathcal{G}_{ab}$ denote the Ricci scalar and the Einstein tensor, respectively (corresponding to the spacetime metric $g$) .

Even a theory with second order equations of motion cannot automatically be considered as a viable physical theory, unless it possesses a well-posed initial value formulation (Cauchy problem). This means that the theory can be formulated in terms of a system of time evolution equations that satisfies the following two physically reasonable properties. Given suitable initial data (that satisfies certain constraints), then i) a unique solution must exist to the evolution equations and ii) the solution must depend continuously on the initial data (in a suitable norm). Apart from being a crucial mathematical criterion, the problem of well-posedness is also interesting from a numerical and experimental point of view. Observations of gravitational wave signatures from coalescing black holes provide new possibilities to test general relativity and its alternatives. In order to do this, however, one needs a stable numerical scheme to solve the equations of motion which is not possible without a well-posed initial value formulation of the theory. Ill-posed systems are unsuitable for numerical computer simulations because initial numerical errors tend to grow drastically during the evolution.

In general, the equations of motion in these theories is a nonlinear system of PDEs. A generic feature of these types of equations is that there may not exist global in time solutions for all data, solutions may blow up in a finite time. The best one can hope for (at least for generic initial data) is to establish local well-posedness, that is to say, the above two criteria are only required to hold for a finite (but strictly non-zero) time. To establish local well-posedness of the non-linear equations, it is sufficient to study the properties of the highest derivative (principal) terms in the linearized equations of motion in a generic background.

The set of algebraic conditions on the principal terms in a system of PDEs that are relevant for well-posedness is generally called hyperbolicity. There exist multiple notions of hyperbolicity in the literature such as weak, strong, strict and symmetric hyperbolicity. It can be shown [3] that strongly hyperbolic equations possess a locally well-posed initial value formulation, whereas weaker notions of hyperbolicity do not guarantee that (see definitions and more precise statements in Appendix A) For this reason, strong hyperbolicity is considered to be a minimal requirement to perform numerical simulations and in this paper, we shall be primarily concerned with strong hyperbolicity.

The problem of well-posedness in theories of gravity is exacerbated by the fact that most of these theories are diffeomorphism-covariant. This implies that solutions to the equations of motion are never unique in a mathematical sense. This problem is usually solved by fixing the gauge. This involves imposing a condition on the components of the metric and/or its derivatives, and modifying the equations of motion by terms that vanish when the gauge condition holds. In general relativity the simplest gauge condition leading to a well-posed formulation is provided by the harmonic gauge [7, 16] which reduces the Einstein equations to a system of quasilinear wave equations on each component $g_{\mu\nu}$.

The question of well-posedness in modified theories is far from settled. In fact, for non-trivial extensions of GR, it had not been discussed until recently. In [5, 6] Papallo and Reall studied the linearized equations of motion of Lovelock [2] and Horndeski theories in a generic weak field background and in (generalized) harmonic gauge. The restriction to a weak field background is due to the fact that in these theories well-posedness is known to break down when the background fields are strong [5, 18, 19, 20, 21]. However, one may still hope that the theory does not lose its predictive power when it is used to describe small deviations from GR. It turns out that the only subclass of these theories that has well-posed Cauchy problem in a generalized harmonic gauge is the so-called $k$-essence-type theories ($G_{3}=\partial_{X}G_{4}=G_{5}=0$). It should be noted that the linearized equations of motion in theories with non-trivial $G_{3}$ were shown to be strongly hyperbolic in a specific choice of generalized harmonic gauge. However, this result does not extend to the non-linear case in a generic background.

Despite this result, one cannot immediately conclude that more general theories are useless. There exists several different well-posed formulations of GR so one might hope that some other formulation and a different choice of gauge could be suitably extended to more general Horndeski theories. In this paper, we focus on the cubic subclass $(G_{4}=G_{5}=0)$ of Horndeski theories222This is related to the theory given by $G_{2}\neq 0$, $G_{3}=b(\phi)X$, $G_{4}=G_{4}(\phi)$ and $G_{5}=0$ by a field redefinition. and study its initial value formulation in more detail. Recently, this class of theories has received some attention among cosmologists due to the fact that this model can describe a non-singular bouncing universe [29, 30, 31, 32]. Furthermore, cubic Horndeski theories naturally arise as certain low energy limits of massive gravity theories [38, 39].

The main result of this paper is that cubic Horndeski theories do possess a well-posed initial value formulation (at least in the weak field regime) and we provide three examples of strongly hyperbolic formulations. Since readers with different backgrounds may find different parts of this paper interesting, we intend to organise our results accordingly.

We begin with a general discussion of the ADM formulation of cubic Horndeski theories in Section 2. More specifically, we present the standard ADM evolution and constraint equations of cubic Horndeski theories and show that a suitable linear combination of these equations give a scalar evolution equation which contains no second derivatives of the spacetime metric. This observation has already been made in [29] but we emphasize this fact here again, since it is a key step to obtain well-posed formulations. The section is concluded by a preliminary discussion of constraint propagation.

Section 3 is mathematically more involved333Some additional information is provided on pseudodifferential calculus and its applications to hyperbolic and elliptic PDEs in Appendix A. so readers interested in numerical applications may jump straight to Section 4. In Section 3 we present an elliptic-hyperbolic formulation of cubic Horndeski theories, using ideas put forward by Andersson and Moncrief in [8] for vacuum GR. After briefly reviewing [8], we show how a suitable modification of the constant (or arbitrarily prescribed) mean curvature and spatial harmonic gauge conditions lead to second order elliptic equations for the lapse function and the shift vector. In the weak field regime and on slices with negative Ricci curvature, existence and uniqueness of solutions to these elliptic equations is guaranteed. It is finally shown that under these assumptions, the strong well-posedness result of Andersson and Moncrief for GR extends to cubic Horndeski theories.

In Section 4 we consider a version of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [23, 24] with a generalized Bona-Massó slicing condition [25] and non-dynamical (i.e. arbitrary but a priori fixed) shift vector. This formulation contains $2$ free parameters: one that describes the slicing condition and one that describes how we modify the evolution system by the momentum constraint. It is shown that when these parameters obey a lower bound then the system of equations is strongly hyperbolic in the weak field regime.

Finally, Section 5 is the most relevant to those with interests in non-linear numerical computer simulations. Here we review the so-called covariant conformal Z4 (CCZ4) formulation [26] which was constructed to enhance the accuracy of numerical simulations in GR. This was achieved by an appropriate modification of Einstein’s equation so that constraint violations are damped away during the evolution. Together with a $2$-parameter family of dynamical gauge conditions (generalized Bona-Massó slicing and “gamma driver” conditions), a straightforward generalization of the CCZ4 system to cubic Horndeski theories constitutes a strongly hyperbolic system of PDEs whenever a simple lower bound on these parameters is imposed and the fields are sufficiently weak. In particular, the slicing conditions selected by strong hyperbolicity include the 1+log slicing which is used in many numerical applications. We also comment on the issue of constraint damping in cubic Horndeski theories.

2 Setting up the problem

2.1 Equations of motion

In this section we provide the ideas that all three formulations (presented in the subsequent sections) share.

We adapt the following notation. We are going to use the Latin letters $(a,b,c,...)$ for abstract indices and Greek letters $(\mu,\nu,\rho,...)$ for coordinate indices. The Latin letters $(i,j,k,...)$ will be used for spatial indices. As mentioned before, we use calligraphic letters ($\mathcal{R}$, $\mathcal{R}_{ab}$, $\mathcal{G}_{ab}$, etc.) for spacetime curvature tensors, whereas curvature tensors defined on spatial slices are denoted by regular ($R$, $R_{ab}$, etc.) letters. For a metric $m$, $|\cdot|_{m}$ denotes the pointwise norm with respect to $m$ (e.g. for a vector field $v^{a}$ we have $|v|_{m}=m_{ab}v^{a}v^{b}$). Our convention on the metric signature is $(-,+,+,+)$.

As mentioned in the Introduction, the class of theories under consideration can be described by the action ($X\equiv-\frac{1}{2}(\partial\phi)^{2}$)

[TABLE]

The reason for separating out $X$ in this action is that we are going to view Horndeski theories as small deformations of Einstein’s theory with a minimally coupled scalar field (referred to as the Einstein-scalar-field theory later, similarly to [5]).

Varying the action (2) with respect to the metric yields the equation of motion [6]

[TABLE]

In the ADM-type formulations of general relativity, it is often beneficial to use the linear combination $E_{ab}-\frac{1}{2}Eg_{ab}=\mathcal{R}_{ab}=0$ as equation of motion, rather than $E_{ab}=\mathcal{G}_{ab}=0$. In fact, it turns out that it is useful to consider the same combination of the gravitational equations of motion in Horndeski theories:

[TABLE]

In the scalar evolution equation

[TABLE]

(obtained by varying the action (2) with respect to $\phi$) the only term involving second derivatives of the metric is $\mathcal{R}_{ab}\nabla^{a}\phi\nabla^{b}\phi$. We will see that it is useful to express this from $\left(E_{cd}-\frac{1}{2}Eg_{cd}\right)\nabla^{c}\phi\nabla^{d}\phi$. In other words, instead of the equation $E_{\phi}=0$, we are going to use

[TABLE]

as the scalar evolution equation. The reason for this is that this equation contains derivatives of the scalar field up to second order and derivatives of the metric only up to first order. (Some benefits of the use of this particular linear combination were also noticed in [29].)

Now we assume that the spacetime manifold $(M,g)$ is globally hyperbolic $M=\mathbb{R}\times\Sigma$ and $h_{ab}$ is the spatial metric induced on the spacelike Cauchy surfaces $\Sigma_{t}$. Let $n^{a}$ be the future directed unit normal to $\Sigma_{t}$. The lapse function $N$ and the shift vector $N^{a}$ are then defined by

[TABLE]

The convention on the extrinsic curvature used here is

[TABLE]

We also need ADM variables for the derivatives of the scalar field, let

[TABLE]

and

[TABLE]

For convenience, we also introduce a fixed, smooth background metric on the spatial slices $\mathring{h}$ and denote the corresponding covariant derivative and Christoffel symbol by $\mathring{D}$ and $\mathring{\Gamma}^{i}_{jk}$, respectively.

Now we are ready to provide the standard ADM-type equations of motion in cubic Horndeski theories. Taking the spatial projection of (2.1) in both indices yields the tensor evolution equation (some helpful formulas for carrying out ADM decompositions are provided in Appendix B)

[TABLE]

Note that

[TABLE]

Similarly to general relativity, the projections $E_{ab}n^{a}n^{b}$ and $E_{cb}h^{c}_{a}n^{b}$ in Horndeski theories yield constraint equations: the Hamiltonian constraint is

[TABLE]

while the momentum constraint reads as

[TABLE]

Even though we are not going to use the explicit form of the scalar evolution equation (2.1), only some of its properties, we rewrite it in terms of the ADM variables, for reference. One obtains444The Mathematica package xAct [42] was of great help in the derivation of the equations.

[TABLE]

where we used the following auxiliary variables

[TABLE]

In the modifications of the ADM formulation considered in this paper, the system of evolution equations takes the general form

[TABLE]

[TABLE]

where $u$ and $v$ are both column vectors of size $n$, corresponding to the dynamical variables; $\mathbf{a}$ and $\mathbf{b}$ are $n\times n$ matrices, depending on the fields $u$, $v$ (and their derivatives). We assume that the matrix $\mathbf{a}$ is invertible so that the hypersurfaces $\Sigma_{t}$ are non-characteristic and spacelike.

In Horndeski theories, for example, the variables $u$ include fields like the components of the induced metric $h_{ij}$ and the scalar field $\phi$, while the variables $v$ are auxiliary variables such as the extrinsic curvature $K_{ij}$ and the normal derivative of the scalar field $A$. Since the formulations to be discussed have different dynamical fields, this will be made explicit later, on a case-by-case basis.

The hyperbolicity of the system of equations (19a,19b) is found by analyzing its characteristic equation. This is obtained by linearizing the evolution equations, selecting the highest derivative (principal) terms and replacing all the derivatives with $\partial_{\mu}\to i\xi_{\mu}\equiv i(\xi_{0},\xi_{i})$. Note that in equation (19a) the principal terms are the terms proportional to $v$ and first derivatives of $u$; whereas in equation (19b) the principal terms are first derivatives of $v$ and second derivatives of $u$. The characteristic equations have the general form

[TABLE]

where $\mathbb{A}$ is a $2n\times 2n$ matrix depending on the background fields, $\mathbb{L}(\xi_{k})$ is also a $2n\times 2n$ matrix depending on the background fields and the spatial Fourier variable $\xi_{k}$, $U$ is a size $2n$ column vector associated with the dynamical fields (i.e. $u$ and $v$). The condition that $t=$constant surfaces are non-characteristic implies that $\mathbb{A}$ is invertible.

Equation (20) can be regarded as the eigenvalue problem for the matrix $\mathbb{M}(\xi_{k})\equiv\mathbb{A}^{-1}\mathbb{L}(\xi_{k})$. The system (19a,19b) is called weakly hyperbolic if and only if the eigenvalues of $\mathbb{M}(\xi_{k})$ are real for any $\xi_{k}$ of unit norm $\xi_{k}\xi_{k}=1$. When the matrix $\mathbb{M}(\xi_{k})$ i) has real eigenvalues, ii) is diagonalizable and iii) has a complete set of eigenvectors that depend smoothly on $\xi_{k}$, for any $\xi_{k}$ of unit norm, then the system (19a,19b) is said to be strongly hyperbolic. In order for the Cauchy problem of (19a,19b) to be well-posed, the system must be at least strongly hyperbolic. For more precise definitions and statements, see Appendix A and references cited therein.

In particular, for cubic Horndeski theories555For different formulations of the theory, the number of dynamical variables may be different so we continue the discussion with general $n$., we can write $U=(U_{g},U_{\phi},U_{A})^{T}$ where $U_{g}$ is a size $2n-2$ vector that corresponds to the gravitational variables (e.g. $h_{ij}$ and $K_{ij}$); the components $U_{\phi}$ and $U_{A}$ correspond to the variables $\phi$ and $A$. Since the equations (9) and (2.1) contain no principal terms associated with the spacetime metric $g$, the matrices $\mathbb{A}$ and $\mathbb{M}$ have the upper triangular form

[TABLE]

[TABLE]

The matrix blocks labelled by subscripts $gg$, $g\phi$ and $\phi\phi$ have sizes $(2n-2)\times(2n-2)$, $2\times(2n-2)$ and $2\times 2$, respectively. It is also worth noting that the matrices $\mathbb{A}_{gg}$, $\mathbb{A}_{g\phi}$, $\mathbb{L}_{gg}$, $\mathbb{L}_{g\phi}$ depend only on the fields $u$, $Du$ and $v$, that is to say, the tensor (gravitational) evolution equations are quasilinear (see e.g. equation (2.1)).

It will be useful (especially in Section 5) to separate the Einstein-scalar-field theory and the Horndeski (i.e. $G_{2}$ and $G_{3}$-dependent) parts in $\mathbb{A}$ and $\mathbb{L}$:

[TABLE]

where $\mathbb{L}_{0}$ and $\mathbb{M}_{0}$ correspond to the Einstein-scalar-field theory, $\delta\mathbb{A}$ and $\delta\mathbb{L}$ are the Horndeski terms. The specific forms of these terms will also be given on a case-by-case basis.

Let us consider in more detail the characteristic equation corresponding to (2.1) (or (2.1)). Selecting the second derivatives of $\phi$ in the linearized version of (2.1) (recall that there are no second derivatives of $g$ in (2.1)) and substituting the derivatives $\partial_{\mu}\to i\xi_{\mu}\equiv i(\xi_{0},\xi_{i})$, the characterstic equation for the scalar mode is given by

[TABLE]

with

[TABLE]

and

[TABLE]

The notations $P_{g\phi}$, $P_{\phi\phi}$ and $P^{\prime}_{\phi\phi}$ refer to the coefficients of the second derivatives of $\phi$ in the linearized versions of equations (2.1), (2.1) and (2.1), respectively.

The same characteristic equation corresponding to the scalar degree of freedom has been previously found in [5] and [29]. More precisely, to make a comparison with [5], we note that with the preferred gauge choice ${\mathcal{H}}_{ab}=-\partial_{X}G_{3}\nabla_{a}\phi\nabla_{b}\phi$ (and $G_{4}=0$) made therein, equation (230) of [5] agrees with (24).

In a regime in which the fields are sufficiently weak, $P_{\phi\phi}^{\prime}$ is close to the spacetime metric $g$ and therefore, it is a Lorentzian metric. By sufficiently weak fields we mean field configurations such that the Horndeski terms are small compared to the Einstein-scalar-field terms. More precisely, let $E=\text{max}\left\{|\mathcal{R}_{\mu\nu\rho\sigma}|^{1/2},|\nabla_{\mu}\phi|,|\nabla_{\mu}\nabla_{\nu}\phi|^{1/2}\right\}$ in all orthonormal bases. Then the weak field condition is equivalent to

[TABLE]

Regarding (24) as an equation for the “characteristic speeds” $\xi_{0}$ for given $\xi_{i}\neq 0$, this equation has two distinct real solutions $\xi_{0}^{\phi,\pm}$. Furthermore, the weak field assumptions (2.1) also ensure that the spacelike $t=$constant hypersurfaces are non-characteristic.

2.2 Constraint propagation

In addition to studying the hyperbolicity of the equations of motion in different formulations, we need to address the issue of constraint propagation. That is to say, we need to check whether solutions to the equations used in these formulations remain solutions of the original Horndeski equations of motion during the evolution. Here we present a fairly detailed derivation of the equations governing the propagation of gauge conditions and constraints, even though the individual steps are quite standard. The purpose of this is to demonstrate that the Bianchi identity (and its generalization) leads straightforwardly to a homogeneous system of PDE for the constraint variables, without making any reference to a specific form of the equations of motion. The only assumption we make is that the equations of motion are second order PDEs obtained by varying a diffeomorphism invariant action. This guarantees that the normal projection of the equations of motion is a constraint equation. Note that in this section we derive the equations without gauge fixing, the effect of the gauge fixing terms on the constraint propagation system will be discussed later. Let $E_{ab}=0$ be the equations of motion obtained by varying the action with respect to the spacetime metric $g_{ab}$. Let us decompose it as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These variables denote the spatial evolution equation, the Hamiltonian constraint and the momentum constraint, respectively.

First, we consider

[TABLE]

using $\mathbf{E}_{ab}n^{b}=0$, $\mathbf{M}_{a}n^{a}=0$, $n_{a}n^{a}=-1$ and $\nabla^{a}n_{a}=-K$. Furthermore, the following identities hold:

[TABLE]

We would like to use the spatial projection of the trace reversed version of $\mathbf{E}_{ab}$ as evolution equation666The reason why we prefer to use $\mathcal{E}_{ab}=0$ as evolution equation rather than $\mathbf{E}_{ab}=0$ is that in general relativity the latter approach yields only a weakly hyperbolic system of equations for the constraint variables [4]., i.e.

[TABLE]

Hence, we set

[TABLE]

With these identities we have

[TABLE]

Next, we consider the spatial projection

[TABLE]

In this case we use

[TABLE]

to obtain

[TABLE]

(cf. eqn. (103,104) in [4]).

Now we consider the generalized version of the Bianchi identity, that is,

[TABLE]

which is a consequence of the diffeomorphism invariance of the Horndeski action. Recall that we wish to use $\tilde{E}_{\phi}=0$ as the scalar evolution equation. For this reason we set

[TABLE]

Putting together equations (33), (36), (37) and (38) gives the equations governing the evolution of the momentum and Hamiltonian constraints

[TABLE]

[TABLE]

In both of these equations, the terms in the second line arise due to the fact that we use equation (2.1) instead of (2.1) as the scalar equation of motion. Note that in each of the formulations studied in this paper, the tensor evolution equation is modified with gauge fixing terms. In other words, the equation $\mathcal{E}_{ij}=0$ is replaced with a different equation which introduces additional terms into equations (2.2-2.2). This will be analysed on a case-by-case basis.

3 Elliptic-Hyperbolic formulation

3.1 Review of Andersson and Moncrief’s results

In this subsection we briefly summarize the work done by Andersson and Moncrief in [8] on the vacuum Einstein’s equations. Firstly, we describe how they derived a coupled elliptic-hyperbolic system equivalent to the vacuum Einstein’s equations. Then we sketch their arguments establishing local well-posedness.

We start from the ADM formulation in which the vacuum Einstein equations

[TABLE]

are rewritten as two sets of first order in time evolution equations

[TABLE]

[TABLE]

[TABLE]

complemented by the Hamiltonian constraint

[TABLE]

and the momentum constraint

[TABLE]

Andersson and Moncrief consider a modified version of the system (42) by imposing constant mean curvature (CMC) slicing777We see that the mean curvature $K$ is constant over the slices $\Sigma_{t}$, but not necessarily in time. As it is mentioned in [8], they could also have considered a prescribed mean curvature slicing, i.e. $K=s(t,x)$.

[TABLE]

and a spatial harmonic (SH) gauge condition

[TABLE]

The evolution equations are the defining equation of the extrinsic curvature (42a) and (42b) modified by adding $-D_{(i}V_{j)}$ to the RHS of (42b):

[TABLE]

[TABLE]

[TABLE]

Equations (42c,42d) are replaced by the modified constraints which can be regarded as the equations which determine the lapse function and the shift vector:

[TABLE]

[TABLE]

Equation (45c) can be obtained by taking the trace of (45b), using the Hamiltonian constraint to trade in the Ricci curvature $R$ for lower order terms and using the CMC condition to set $\left(\partial_{t}-\mathcal{L}_{N}\right)K=1$.

Equation (45d) can be derived as follows. Taking the time derivative of $V^{i}$ and commuting $\partial_{t}$ with $h^{kl}$ and the spatial derivatives, one easily obtains

[TABLE]

Using the momentum constraint, the CMC slicing condition ($D_{i}K=0$) and the spatial harmonic condition $(\partial_{t}-\mathcal{L}_{N})V^{i}=0$ to eliminate second derivatives of the spatial metric and first derivatives of the extrinsic curvature then yields (45d). It is worth emphasizing that the CMC slicing condition was used to arrive at both equation (45c) and (45d).

Now we move on to the question of well-posedness of the system (45a-45d) in Sobolev spaces and consider initial data $h_{ij},N,N^{i}\in H^{s}$ and $K_{ij}\in H^{s-1}$ ($s>\tfrac{5}{2}$ in 4-dimensional spacetime) that satisfies the Hamiltonian and momentum constraints.

The modified constraints (45c,45d) are equations relating derivatives of $N,N^{i}$ up to second order to derivatives of $h_{ij}$ up to first order (including the extrinsic curvature), when written in coordinates. This statement is obvious for (45c) but one can check that the second derivatives of $h_{ij}$ cancel each other out on the LHS of (45d). More precisely, the modified constraints have the form

[TABLE]

with $u=(N,N^{i})^{T}$ and $\mathbf{A}$ being a second order, linear elliptic differential operator, with coefficients depending only on the spatial metric, its first spatial derivatives and the extrinsic curvature. Moreover, the elliptic operator $\mathbf{A}$ is lower triangular:

[TABLE]

with

[TABLE]

Standard results in the theory of elliptic PDEs (see Appendix A.2) show that the scalar elliptic operator $-D^{i}D_{i}+K_{ij}K^{ij}$ is an isomorphism $H^{s}\to H^{s-2}$. Furthermore, it is proved in [8] that the elliptic operator $\mathbf{C}^{i}_{~{}j}$ is an isomorphism $H^{s}\to H^{s-2}$, if $(\Sigma,h)$ is a compact manifold with negative Ricci curvature888Note that this formulation was ultimately used to prove a global existence theorem for non-linear perturbations of spatially compact versions of FRW spacetimes with $k=-1$ [9].. These results and the lower triangular structure of $\mathbf{A}$ then implies that $\mathbf{A}$ is also an isomorphism $H^{s}\to H^{s-2}$. Therefore, if we a priori assume that $h_{ij}\in H^{s}$ and $K_{ij}\in H^{s-1}$, then the unique solutions $N,N^{i}$ to the modified constraints are in $H^{s+1}$, i.e. they have an extra regularity compared to $h_{ij}$. This is necessary because this way the terms involving first derivatives of $N^{i}$ in (45a) and second derivatives of $N$ in (45b) are non-principal.

It follows that by solving the modified constraints to determine $N$ and $N^{i}$, the evolution equations become a first order quasilinear system of pseudodifferential equations. It is easy to see now that this system of evolution equations is strongly hyperbolic. Linearizing the evolution equations and substituting the derivatives with the Fourier variables $(\partial_{t},\partial_{k})\to(i\xi_{0},i\xi_{k})$, the coefficients of the highest derivative (principal) terms define the principal symbol. Based on the definitions and general arguments presented in Appendix A.3, strong hyperbolicity means that the principal symbol has real eigenvalues with a complete set of eigenvectors with smooth dependence on $\xi_{i}$. Linearizing (45a) and (45b) around a generic solution $h_{ij}\to h_{ij}+\gamma_{ij}$, $K_{ij}\to K_{ij}+\kappa_{ij}$, the eigenvalue problem of the principal symbol reads as

[TABLE]

where $\mathbb{M}_{0}$ is the $2\times 2$ block matrix (recall that the terms involving second derivatives of $N$ and first derivatives of $N^{i}$ are non-principal)

[TABLE]

It is easy to see that the principal symbol has eigenvalues $\xi_{0}^{\pm}=N^{k}\xi_{k}\pm N|\xi|_{h}$ with a complete set of eigenvectors: for any symmetric matrix $u_{ij}$,

[TABLE]

is an eigenvector with $\xi_{\mu}=(\xi_{0},\xi_{i})$. Note that this means that $\xi_{\mu}$ is a null vector. Since all eigenvalues are real and the eigenvectors can be chosen to be independent of $\xi_{k}$, the system of evolution equations is strongly hyperbolic when the modified constraints are solved. (In fact, it is symmetric hyperbolic so one can demonstrate well-posedness by standard energy methods in physical space, as was done in [8].)

3.2 Equations of motion and gauge fixing in cubic Horndeski theories

We will now show that the above formalism can be extended to cubic Horndeski theories. For this, we first discuss the generalization of the SH-CMC gauge condition. Recall that the CMC condition in general relativity was used to set $(\partial_{t}-\mathcal{L}_{N})K$ to $1$ in the trace of the evolution equation. When taking the trace of (2.1), it is possible to get a constraint equation by a choice of an appropriate slicing condition which sets the terms involving $(\partial_{t}-\mathcal{L}_{N})K$ and $(\partial_{t}-\mathcal{L}_{N})A$ to an a priori fixed function. For this reason, we take an approach very similar to the one above: in the trace of (2.1) we trade in $R$ using the Hamiltonian constraint to get

[TABLE]

We are seeking a gauge condition of the form

[TABLE]

to eliminate the time derivatives in equation (3.2). Taking the normal derivative of (54) gives

[TABLE]

Therefore, the desired choice is an $f$ satisfying

[TABLE]

Note that this slicing condition has an interesting relationship with the canonical momentum $\pi^{ij}$ conjugate to $h_{ij}$. If we switch to a Hamiltonian description, equation (3.2) (which is the trace of (2.1)) is equivalent to the trace of

[TABLE]

where $\mathcal{H}$ is the Hamiltonian. Hence, it is clear that the time differentiated terms in (3.2) come from $\partial_{t}(h_{ij}\pi^{ij})$, that is, the preferred slicing condition is equivalent to $\pi=s(x,t)$.

Rewriting $G_{3}(\phi,X)$ as a function depending on $\phi$, $A$, $A_{i}$ and $h^{ij}$, the condition (56) can be integrated in $A$ and so $f$ can be determined up to the addition of an arbitrary function of $\phi$, $h_{ij}$ and $A_{i}$. Hence, the elliptic equation for $N$ reads as

[TABLE]

Based on standard results in the theory of elliptic PDEs (see Appendix A.2 and the references cited therein), this equation has a unique solution for $N$ if the coefficient of the zeroth order term (i.e. terms enclosed by the curly brackets $\{\}$) is non-negative. This condition is satisfied if the Horndeski terms are small corrections to GR, i.e. in the weak field regime. Note that the weak field requirements on $f$ are

[TABLE]

with $E=\text{max}\left\{|\mathcal{R}_{\mu\nu\rho\sigma}|^{1/2},|\nabla_{\mu}\phi|,|\nabla_{\mu}\nabla_{\nu}\phi|^{1/2}\right\}$ in an orthonormal basis.

The generalization of the spatial harmonic gauge is more straightforward: we require

[TABLE]

for some $H^{i}$ ($V^{i}$ is as in equation (44)). Once again, we can derive an elliptic equation for the shift vector by requiring $(\partial_{t}-\mathcal{L}_{N})J^{i}=0$ and eliminating derivatives of the extrinsic curvature by using the momentum constraint and the generalized CMC condition. The result is

[TABLE]

Furthermore, using

[TABLE]

and

[TABLE]

gives

[TABLE]

The operator on the LHS acting on $N^{i}$ is exactly the same as in GR (c.f. (49)) which means that (3.2) has a unique solution for $N^{i}$ on spatial slices with negative spatial Ricci curvature.

Finally, we need to decide how to use the generalized spatial harmonic condition in the evolution equations. As mentioned before, we do not modify the scalar equation $\tilde{E}_{\phi}$, only the tensor equation in the most natural way, i.e. by replacing $\mathcal{E}_{ij}$ by $\tilde{\mathcal{E}}_{ij}\equiv\mathcal{E}_{ij}-D_{(i}J_{j)}$:

[TABLE]

Therefore, the Cauchy problem for cubic Horndeski theories can be formulated as follows. Consider initial data $h_{ij},N,N^{i},\phi\in H^{s}$ and $K_{ij},A\in H^{s-1}$ ($s>\tfrac{7}{2}$)999The lower bound on $s$ is stronger than in vacuum GR due to the fact that the scalar evolution equation (2.1) is not quasilinear (see Appendix A.3). that satisfies the Hamiltonian and momentum constraints. Then the system of equations to be solved consists of the evolution equations (8), (9), (2.1) and (3.2), together with the elliptic equations (3.2), (3.2).

3.3 Constraint propagation

Before moving on to the question of well-posedness, in this section we explain how to get the equations describing the propgation of the gauge conditions and the original constraints using the gauge-fixed equations of motion. As described in the previous section, we use $\tilde{\mathcal{E}}_{ij}\equiv\mathcal{E}_{ij}-D_{(i}J_{j)}=0$ as tensor evolution equation. Following the argument started in section 2.2, we set $\tilde{\mathcal{E}}_{ij}=\tilde{{E}}_{\phi}=0$ and switch to the new variables

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

we have the following system of homogeneous linear evolution equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The first two equations follow easily by recalling the steps we used to get the elliptic equations (45c,45d) from the evolution equations and the constraints. To show that the quantities $(\mathbf{F},\mathbf{J}_{i},\mathbf{\tilde{H}},\mathbf{\tilde{M}}_{i})$ remain zero during the evolution, we first note that it follows from equations (67a,67b) that if $(\mathbf{F},\mathbf{J}_{i},\mathbf{\tilde{H}},\mathbf{\tilde{M}}_{i})$ vanish initially then $\partial_{t}\mathbf{F}=\partial_{t}\mathbf{J}_{i}=0$ on the initial surface. It turns out that one can obtain a simple energy estimate for the system (67) without solving the eigenvalue problem of the principal symbol. Consider the energy 101010The steps of the proof of constraint propagation, as well as the expression for the energy, are the same as in [8] (Lemma 4.1.), since the extra terms entering to the equations due to the presence of a scalar field are non-principal. Nevertheless, we provide a sketch of the proof for completeness.

[TABLE]

Specifically, we want to show that $|\partial_{t}E_{\text{constraint}}|\leq CE_{\text{constraint}}$ for some constant $C(h,K,N)$. Clearly, the action of $\partial_{t}$ on the volume form can be bounded by a constant. When the time derivative acts on the gauge and constraint quantities, we use (67) to exchange the time derivatives. Since the energy (68) is invariant under spatial diffeomorphisms, the terms involving Lie derivatives will vanish.

The non-principal terms can be estimated by the energy itself. For example,

[TABLE]

[TABLE]

The potentially problematic (principal) terms are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\simeq$ denotes equivalence up to principal terms and $\sim$ denotes equivalence of the integrands up to integration by parts. We see that the terms containing higher derivatives cancel each other out, giving us the desired result. Therefore, if $E_{\text{constraint}}$ vanishes initially, then it remains zero during the evolution as well, implying $(\mathbf{F},\mathbf{J}_{i},\mathbf{\tilde{H}},\mathbf{\tilde{M}}_{i})=0$.

3.4 Proof of strong hyperbolicity

Now we linearize the equations of motion (8), (9), (2.1), (3.2), (3.2), (3.2) around a generic weak field background with negative spatial Ricci curvature. The linearized quantities are as follows:

[TABLE]

Now we assume that the elliptic equations (3.2) and (3.2) have unique solution for any $(\gamma_{ij},\kappa_{ij},\psi,a)$ satisfying the gauge conditions and constraint equations. Recall that this is true for weak fields and when the background spacetime has negative spatial Ricci curvature.

The important difference compared to GR is that in this case the first derivatives of $\beta^{i}$ appearing in the defining equation of the extrinsic curvature (8) and the second derivatives of $\alpha$ appearing in the tensor evolution equation (3.2) are principal terms. To see how these terms affect the hyperbolicity of the evolution equations, we have to look at the principal terms in the linearized versions of the elliptic equations (3.2), (3.2):

[TABLE]

[TABLE]

Solving the linearized elliptic equations for given $\psi$, $a$, $\gamma_{ij}$ and $\kappa_{ij}$, one obtains a (non-local) map $(\hat{\alpha},\hat{\beta}):(\psi,a,\gamma,\kappa)\to(\alpha,\beta)$. This solution map is a pseudodifferential operator of class $\mathcal{OP}_{\text{cl}}^{-2}$ with principal symbol (see Appendix A)

[TABLE]

[TABLE]

where the $~{}\hat{}~{}$ on the fields $\psi$ and $a$ denotes the Fourier transform. Note that the principal symbol of the solution map does not depend on $\hat{\gamma}$ and $\hat{\kappa}$.

By solving the elliptic equations, the evolution equations take the form of a first order pseudodifferential system of evolution equations. To determine the hyperbolicity of this system, we consider the eigenvalue problem of the principal symbol, c.f. (20). The vector $U$ is now a column vector of size $14$

[TABLE]

The characteristic equation can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the specific form of the linearized version of (2.1) is quite long and unessential for our purposes. Nevertheless, we note that substituting (79b) into (79d) gives (recall equation (24))

[TABLE]

It follows from the upper triangular structure of the characteristic equation (see Section 2) that

[TABLE]

is an eigenvector of (79a-79d) with eigenvalues $\xi_{0}^{\pm}=N^{k}\xi_{k}\pm N|\xi|_{h}$ for any symmetric $\hat{u}_{ij}$. (Recall that $\xi_{\mu}=(N^{k}\xi_{k}\pm N|\xi|_{h},\xi_{i})$ is a null vector with respect to the spacetime metric.) One can easily find $6$ linearly independent vectors $\hat{u}_{ij}$. Taking into account the two sign choices in $\xi_{0}^{\pm}$ and in (81), this gives $12$ eigenvectors: a $6$-dimensional eigenspace for both eigenvalues $\xi_{0}^{\pm}=N^{k}\xi_{k}\pm N|\xi|_{h}$. Note that these include the $2$ pairs of transverse-traceless modes corresponding to the gravitational degrees of freedom, i.e.

[TABLE]

with $\hat{u}_{ij}^{TT}$ satisfying $h^{ij}\hat{u}_{ij}^{TT}=0$ and $\xi^{i}\hat{u}_{ij}^{TT}=0$.

The remaining eigenvalues $\xi_{0}^{\phi,\pm}$ are found by solving

[TABLE]

(see Section 2 for notation). Recall that for weak fields, $\xi_{0}^{\phi,\pm}$ are distinct and real. The corresponding eigenvectors have the form

[TABLE]

with

[TABLE]

[TABLE]

These expressions are clearly smooth functions of $\xi_{i}$ for any choice of $H_{i}$, since $h_{ij}$ is a positive definite metric and $\xi_{i}\neq 0$ by assumption ($\xi_{i}$ has unit norm).

These two “physical” eigenvectors corresponding to the scalar degree of freedom satisfy the high frequency limit of the Hamiltonian and momentum constraints

[TABLE]

[TABLE]

To summarize, we have found that the principal symbol (see equations (79a-79d)) has $14$ real eigenvalues and the corresponding eigenvectors are linearly independent and have smooth dependence on $\xi\in\mathbb{S}^{2}$. This implies that the evolution equations are strongly hyperbolic when the modified constraint equations have a unique solution for arbitrary $(\phi,A,h,K)$. In particular, this is the case in a weak field regime and in spacetimes that can be foliated with generalized prescribed mean curvature slices with negative Ricci curvature.

4 BSSN-type formulation with non-dynamical shift vector

4.1 Equations of motion

The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation (together with its modifications) is a popular method used to numerically integrate Einstein’s equations. Several versions of this approach give rise to a strongly hyperbolic reformulation of the vacuum Einstein equations [11, 12, 13]. Here, we extend this approach to cubic Horndeski theories.

The equations of motion are obtained from the standard ADM equations given in Section 2 as follows. First, we introduce the conformal metric $\tilde{h}_{ij}$ as a new variable, defined by

[TABLE]

where the conformal factor $\Omega$ is

[TABLE]

for an arbitrary smooth background metric $\mathring{h}_{ij}$. Note that this implies that $\det\mathring{h}=\det\tilde{h}$. The inverse conformal metric, denoted by $\tilde{h}^{ij}$ is then

[TABLE]

Next, we define the quantity

[TABLE]

where $\mathring{\Gamma}$ and $\mathring{D}$ denote the Christoffel symbol and the covariant derivative corresponding to $\mathring{h}_{ij}$. Similarly, let $\tilde{D}$ be the covariant derivative corresponding to the metric $\tilde{h}$ and let

[TABLE]

and

[TABLE]

Clearly, the definition of $\tilde{A}_{i}$ does not depend on $\tilde{h}$ but the index of the vector field $\tilde{A}$ is raised and lowered with $\tilde{h}$, whereas the index of $A$ is raised and lowered with $h$. We continue to use a similar convention in the further discussion: indices of tensor fields denoted by letters with a tilde are raised and lowered with $\tilde{h}$, whereas indices of tensor fields without a tilde are raised and lowered with $h$.

The extrinsic curvature is decomposed to its trace and conformal traceless parts

[TABLE]

or alternatively,

[TABLE]

The evolution equations for the variables $\tilde{h}_{ij}$ and $\Omega$ are

[TABLE]

[TABLE]

where $\partial_{0}$ is given by

[TABLE]

To write the remaining equations in a more compact way, we use the conformal versions of the auxiliary variables introduced in (16):

[TABLE]

The evolution equation for $K$ is the trace of the tensor evolution equation, i.e. the same as (3.2), except that the variables $h_{ij}$ and $K_{ij}$ (and the covariant derivatives) are now replaced by the corresponding expressions (86-92).

[TABLE]

The equation describing the evolution of $\tilde{Q}_{ij}$ is obtained by taking the trace free part of (2.1)

[TABLE]

where $\tilde{T}_{ij}^{TF}$ denotes the trace free part of a symmetric tensor $T_{ij}$,

[TABLE]

and the conformal decomposition of the spatial Ricci tensor is

[TABLE]

Finally, the equation for $\tilde{V}^{i}$ is obtained by commuting $\partial_{0}$ with $\mathring{D}$ to get

[TABLE]

and then adding the momentum constraint times $2mN$ to (104) gives

[TABLE]

where $m$ is arbitrary. We will show that for a range of values for the parameter $m$, the system is strongly hyperbolic. Finally, once again we use $\tilde{E}_{\phi}=0$ as the scalar evolution equation (defined in (2.1)) which can also be rewritten in terms of the variables introduced here:

[TABLE]

[TABLE]

Equations (94), (95), (4.1), (4.1), (4.1) and (4.1) must be complemented with the evolution equations for the lapse function and the shift vector. This can be done by choosing an appropriate slicing condition and a spatial coordinate condition. A popular choice for the slicing condition is harmonic slicing and its generalizations. Harmonic slicing means that the harmonic coordinate condition is imposed only on the time coordinate

[TABLE]

Writing this out in terms of the ADM variables gives an evolution equation for the lapse function:

[TABLE]

Sometimes it is more convenient to consider a generalization of this condition, called the Bona-Massó slicing condition

[TABLE]

for a suitable function $F$. The choice $F(N)=\frac{2}{N}$ called the 1+log slicing is the most widely used in numerical applications. We will generalize this condition even further for cubic Horndeski theories

[TABLE]

To simplify the discussion of the linearized equations, we introduce

[TABLE]

There are a number of ways to impose a dynamical gauge condition on the shift vector in general relativity. However, it has also been demonstrated that it is possible to write Einstein’s equations in a strongly (even symmetric) hyperbolic form by choosing an arbitrary (but a priori fixed) shift vector [12]. In these latter formulations, the shift vector is only a source term. In this section we take this latter approach and show that the equations of motion of the subclass of Horndeski theories under consideration can be written in a strongly hyperbolic form with arbitrarily fixed shift vector.

4.2 Proof of strong hyperbolicity

In this section, we are going to show that the BSSN equations of motion (consisting of equations (94), (95), (4.1), (4.1), (4.1), (110), (106), (4.1)) is a strongly hyperbolic system for the dynamical variables. Our strategy will be as follows. First, we linearize the equations and select the highest derivative (principal) terms in the equations. (The list of the variables and their linearized versions is summarized in Table 1.) These are the terms that are at the first derivative level in equations (94), (95), (110), (106) and the ones at the second derivative level in equations (4.1), (4.1), (4.1), (4.1)). The second step is to convert the principal terms in the equations to the eigenvalue problem of the principal symbol. This involves switching to Fourier variables:

[TABLE]

[TABLE]

Note that in this section our choice of basis is slightly different from the one used in Section 3: here we use $\partial_{0}\equiv\partial_{t}-N^{k}\partial_{k}$, rather than $\partial_{0}\equiv\partial_{t}$. Clearly, this only amounts to a shift in the variable $\xi_{0}\to\xi_{0}-N^{k}\xi_{k}$. Nevertheless, from now on, we will denote the eigenvalues of the principal symbol by $\xi_{0}$ in this basis. This includes the solutions $\xi_{0}^{\pm,\phi}$ to the scalar characteristic equation (24).

Next, we solve the eigenvalue problem, by determining the eigenvalues $\xi_{0}$ and the corresponding eigenvectors explicitly. Finally, we show that the conditions of strong hyperbolicity are met for an appropriate choice of the parameters $m$ and $\sigma$: the principal symbol has real eigenvalues and a complete set of eigenvectors with smooth dependence on $\xi_{k}$.

Since the linearization is straightforward, we simply just state the eigenvalue problem for the $20$ variables $(\hat{\alpha},\hat{\omega},\hat{\kappa},\hat{\tilde{\gamma}}_{ij},\hat{q}_{ij},\hat{\tilde{v}}_{i},\hat{\psi},\hat{a})$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the coefficients $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ depend on the background fields and $\xi_{i}$ and their explicit form is not important for our purposes (although it is straightforward to obtain them by linearizing (4.1)). Substituting (79b) into (79d) gives

[TABLE]

that is to say, equation (24). In the discussion of the eigenvalues and eigenvectors it is simpler to use the variable

[TABLE]

instead of $\xi_{0}$.

We can identify a subset of eigenvalues and eigenvectors of (112) by setting $\hat{\psi}=\hat{a}=0$. These $18$ eigenvalues and eigenvectors are the same as in the formulation described in section IV of [12] for vacuum GR:

• I.

Transverse-traceless (physical) modes with null characteristics $\bar{\xi}_{0}=\pm 1$. Since there are two linearly independent transverse-traceless symmetric tensors $\hat{\tilde{\gamma}}_{ij}^{TT}$, the eigenvectors span a $4$-dimensional space.

• II.

Transverse modes with null characteristics, spanning a $2$-dimensional space.

• III.

A $4$-dimensional space of modes with $\hat{\tilde{\gamma}}_{ij}=\xi_{(i}e_{j)}$ for any $e_{i}$ orthogonal to $\xi_{i}$ w.r.t. $\tilde{h}_{ij}$ and $\bar{\xi}_{0}=\pm\sqrt{m}$.

• IV.

Zero speed modes ($\bar{\xi}_{0}=0$) with $\hat{\tilde{\gamma}}_{ij}=-\frac{2}{|\xi|_{\tilde{h}}^{2}}\left(i\xi_{(i}\hat{\tilde{v}}_{j)}+2\xi_{i}\xi_{j}\hat{\omega}\right)^{\text{TF}}$ where $\hat{\omega}$ and $\hat{\tilde{v}}_{i}$ are arbitrary, spanning a space with dimension $4$.

• V.

A two-dimensional space of modes with $\bar{\xi}_{0}=\pm\sqrt{2\sigma}$.

• VI.

A two-dimensional space of modes with $\bar{\xi}_{0}=\pm\sqrt{\frac{4m-1}{3}}$.

The expressions for the eigenvectors are listed in Table 2. These eigenvalues are real and the eigenvectors are smooth functions of their arguments when $\sigma>0$ and $m>\frac{1}{4}$.

For strong hyperbolicity, we need $20$ linearly independent eigenvectors with real eigenvalues and smooth dependence on $\xi_{i}$. Since we have already found $18$, this amounts to finding two additional eigenvalues and eigenvectors with non-trivial $\hat{\psi}$ and $\hat{a}$. The eigenvalues corresponding to these eigenvectors are found by solving the system consisting of equations (112d) and (112h), or equivalently, (24). As mentioned before, this equation has two distinct real solutions $\xi_{0}^{\phi,\pm}$ in a weakly curved background (see Section 2). Therefore, for a weakly curved background, there must indeed be 2 additional eigenvectors. The only thing that needs to be shown is the condition on smooth dependence.

The easiest way to obtain the corresponding eigenvectors is to derive a closed equation containing $\hat{\psi}$, $\hat{a}^{\phi}\equiv\tfrac{1}{N}i\xi_{0}^{\phi,\pm}\hat{\psi}$ and the Fourier transform of the linearized extrinsic curvature

[TABLE]

Equations (112) imply

[TABLE]

with $\hat{\kappa}_{i\xi}\equiv\hat{\kappa}_{ij}\xi^{j}$ and

[TABLE]

Note that in (4.2) we raise and lower indices with $h$.

It follows from (115) and (112) that if the solutions $\hat{\kappa}_{ij}^{\phi,\pm}$ of (4.2) are smooth functions of $\xi_{i}$, then the same is true for the auxiliary variables $\hat{\tilde{q}}_{ij}$, $\hat{\kappa}$, $\hat{\tilde{\gamma}}_{ij}$, $\hat{\tilde{v}}_{i}$, $\hat{\omega}$ and $\hat{\alpha}$. Based on the tensorial structure of (4.2), we look for an eigenvector of the form

[TABLE]

Note that

[TABLE]

and

[TABLE]

Plugging this into (4.2) we get the following system of linear equations for the coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The solution is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In order for $c_{1}$ to be smooth for any $\xi_{i}\in\mathbb{S}^{2}$, $-(\xi_{0}^{\phi,\pm})^{2}+2\sigma N^{2}|\xi|_{h}^{2}$ must be non-zero for any $\xi_{i}$. This can be achieved by choosing $\sigma$ to be large enough. To see why this is true, we first note that the zeros of the function

[TABLE]

define a cone for any $\sigma$, $\xi_{0}$. In the weak field regime, the null cone of $P^{\prime}_{\phi\phi}$ is close to the null cone of the spacetime metric $g$ for any $\xi_{i}\in\mathbb{S}^{2}$, they might even intersect for special values of $\xi_{i}$. Recall that

[TABLE]

Since $P^{\prime}_{\phi\phi}(\xi_{0};\xi_{i})$ and $\mathcal{F}_{\sigma=1/2}(\xi_{0};\xi_{i})$ may intersect for some $\xi_{i}$, and $\xi^{\phi,\pm}_{0}$ are the solutions of $P^{\prime}_{\phi\phi}(\xi_{0};\xi_{i})=0$, this means that the expression on the RHS of (127) could blow up for some $\xi_{i}$ if $\sigma$ is close to (or equal to) $\frac{1}{2}$. To avoid this, we can just choose $\sigma$ to be large enough (i.e. larger than $\frac{1}{2}$), so that the cones given by $\mathcal{F}_{\sigma}(\xi_{0};\xi_{i})$ lie entirely “inside” the null cones of $g$ and $P\textquoteright_{\phi\phi}$. In other words, for an appropriate $\sigma>\tfrac{1}{2}$, $\mathcal{F}_{\sigma}(\xi_{0}^{\phi,\pm};\xi_{i})$ vanishes for no choice of $\xi_{i}$. Choosing larger values of $\sigma$ makes it possible to deal with stronger background fields.

To summarize, we have shown that the equations of motion for cubic Horndeski theories form a strongly hyperbolic system in a version of the BSSN formulation. The system was obtained using a generalization of the harmonic slicing condition and an arbitrary but (non-dynamically) fixed shift vector. The system is strongly hyperbolic for any $m>\frac{1}{4}$ and for suitable $\sigma\equiv\partial_{K}F>\frac{1}{2}$, regardless of how the source function $F$ depends on the scalar field and its derivatives. (For weakly curved backgrounds, choosing a large enough constant, e.g. $\sigma=1$ is enough.) This means that the original harmonic slicing $\sigma=\frac{1}{2}$ does not work for cubic Horndeski theories.111111One might wonder if $c_{1}$ could be set to [math] in (127) by an appropriate choice of $\rho$ and $\rho_{k}$, even for $\sigma=\frac{1}{2}$. However, this particular choice of $\rho$ and $\rho_{k}$ fails to satisfy the integrability condition $\partial_{A}\rho^{k}=\partial_{A_{k}}\rho$, a similar issue to the one encountered in [5] (see equations (237-241) and the corresponding discussion). On the other hand, the so-called 1+log slicing often used in numerical general relativity, corresponds to the choice $\sigma=\tfrac{1}{N}$ and hence remains a good slicing condition as long as $N<2$. Note that for GR the condition for strong hyperbolicity is $m>\frac{1}{4}$ and $\sigma>0$, whereas the condition for symmetric hyperbolicity is $6\sigma=4m-1>0$.

4.3 Propagation of constraints

To show that the solutions of the BSSN system are also solutions of the original Horndeski equations of motion, we derive a system of evolution equations for the Hamiltonian constraint, the momentum constraint and the variable

[TABLE]

and show that the system of equations is strongly hyperbolic. By uniqueness of the solutions to strongly hyperbolic systems, it follows that if the constraints are satisfied initially then they continue to hold throughout the evolution.

Starting from equations (2.2) and (2.2) and setting

[TABLE]

the constraint evolution equations become

[TABLE]

[TABLE]

[TABLE]

We remind the reader that indices of tensors denoted with letters with tilde are raised and lowered with the conformal metric $\tilde{h}$, indices of tensors without a tilde are raised and lowered with the original induced metric $h$.

Note that we have two additional constraints: the tracelessness of $\tilde{Q}_{ij}$ and $\det\mathring{h}=\det\tilde{h}$. Introducing the constraint variables

[TABLE]

and

[TABLE]

it follows easily that

[TABLE]

[TABLE]

Substituting (135) into the time derivative of (136) implies a wave equation for $\mathbf{D}$, decoupled from the rest of the constraint propagation system (132):

[TABLE]

Therefore, it is clear that starting from initial data that satisfies $\mathbf{D}=0$ and $\mathbf{T}=0$, these conditions will continue to hold throughout the evolution. For this reason, now we only need to deal with the system (132).

The eigenvalue problem of the principal symbol of (132) can be written as

[TABLE]

[TABLE]

[TABLE]

The eigenvalues and eigenvectors of the principal symbol are as follows

[TABLE]

[TABLE]

and finally,

[TABLE]

for any vector $e_{i}$ orthogonal to $\xi_{i}$ with respect to $h_{ij}$. Therefore, the principal symbol of the system describing the evolution of the constraints possesses a complete set of smooth eigenvectors with real eigenvalues, provided that $m>\frac{1}{4}$. Hence, the system (132) is strongly hyperbolic.

5 CCZ4-type formulation

5.1 Constraints

In this section we discuss how the so-called covariant conformal Z4 (CCZ4) [26, 27] formulation extends to the class of theories under consideration.

The idea behind the CCZ4 formulation is to introduce a 4-vector field $\mathcal{Z}^{a}$ that measures the deviations from the actual tensor equations of motion $E_{ab}$ (see (2.1)) and add terms containing $\mathcal{Z}^{a}$ and its first derivatives to the equations so that $\mathcal{Z}_{a}=0$ is an attractor of the modified equations. The modification is carried out in the following way:

[TABLE]

or in the trace reversed version

[TABLE]

The parameters $k_{1}$ and $k_{2}$ here are real constants. Splitting the four vector $\mathcal{Z}^{a}$ as $\mathcal{Z}^{a}\equiv Z^{a}+n^{a}\Theta$ with $Z^{a}\equiv h^{a}_{b}\mathcal{Z}^{b}$ and $\Theta\equiv-n^{a}\mathcal{Z}_{a}$, one can write the normal-normal and normal-spatial projections of (142) as

[TABLE]

[TABLE]

where the expressions for the Hamiltonian and momentum constraints ($\mathbf{H}$ and $\mathbf{M}_{i}$) are given by equations (2.1) and (2.1).

When the generalized Bianchi-identity (37) holds, the evolution equations for the Hamiltonian and momentum constraints are

[TABLE]

[TABLE]

Once again, the system (144-5.1) describing the propagation of constraint violations has the same principal symbol as in general relativity, cf. equations (7), (8), (11), (12) of [28]. Therefore, the hyperbolicity and the high frequency behaviour of that system is not altered by the Horndeski terms. This has the following implications for constraint damping. Similarly to [28], one can carry out a preliminary mode analysis by linearizing around a generic weak field configuration and studying the high frequency limit of (144-5.1). For large frequencies, the Horndeski terms become insignificant. Making a plane wave ansatz for the constraint variables then reduces the high frequency limit of (144-5.1) to the same eigenvalue problem as in [28] (see equation (19)). Hence, we come to the same conclusion as in vacuum GR: the real parts of all eigenfrequencies are negative if $k_{1}>0$ and $k_{2}>-1$. This suggests that with such choice of the parameters $k_{1}$ and $k_{2}$, large frequency constraint violating modes will be damped away in cubic Horndeski theories.

5.2 Equations of motion

Next, we provide the full system of evolution equations, in the conformal decomposition. We introduce $\tilde{Z}_{i}\equiv Z_{i}$, $\tilde{Z}^{i}\equiv\tilde{h}^{ij}\tilde{Z}_{j}$ and

[TABLE]

Similarly to the BSSN case, we use auxiliary variables

[TABLE]

to write the equations more compactly. (For a more complete list of formulas for the ADM and conformal decompositions, see Appendix B.)

We use a natural generalization of the harmonic slicing condition (recall $\partial_{0}\equiv\partial_{t}-N^{k}\mathring{D}_{k}$)

[TABLE]

the shift vector is evolved using the standard "Gamma driver" condition [22]

[TABLE]

[TABLE]

The parameters $f$ and $\eta$ are to be chosen and $B^{i}$ is an auxiliary variable. The evolution equations for the variables $\tilde{h}_{ij}$ and $\Omega$ are just defining equations and hence left unaltered compared to (94) and (95):

[TABLE]

[TABLE]

The evolution equation for $\Theta$ is the same as (144) (using the expression of the Hamiltonian constraint)

[TABLE]

Equation (5.1) is no longer kept as a separate equation, instead, it is added to the evolution equation for $\tilde{V}^{i}$:

[TABLE]

This is equivalent to adding the momentum constraint times $2N$ to (104).

The evolution equations for $K$ and $\tilde{Q}_{ij}$ are the same as before, except for the constraint damping terms.

[TABLE]

[TABLE]

Finally, we also have a pair of scalar evolution equations: (106), (4.1).

We conclude this section with a technical remark. For general relativity, it has been noted that the CCZ4 equations of motion can be derived from an action principle [26] (at least if we ignore the lower order terms with $k_{1}$ and $k_{2}$). If we insisted on a similar action principle for cubic Horndeski theories, then upon taking the linear combination (2.1) of the gravitational and scalar equations of motion, the resulting equation would contain principal terms from the ADM decomposition of $\nabla Z$. However, for cubic Horndeski theories, it is more beneficial to keep equation (4.1) as the scalar evolution equation. It appears that for these theories, it is more useful to introduce the $Z$-terms at the level of the equations, rather than at the level of the action.

5.3 Strong hyperbolicity

Now we linearize the system of equations (153),(106), (4.1) and study its hyperbolicity. The linearly small quantities corresponding to the new variables are

[TABLE]

The eigenvalue problem of the principal symbol for the 27 variables $U=(\hat{\alpha},\hat{\omega},\hat{\theta},\hat{\kappa},\hat{\beta}^{i},\hat{b}^{i},\hat{\tilde{u}}^{i},\hat{\tilde{\gamma}}_{ij},\hat{\tilde{q}}_{ij},\hat{\psi},\hat{a})$ can be written as follows.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, again, we do not need to deal with the precise expressions of the coefficients $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$, we only need to keep in mind that substituting (79b) into (79d) yields

[TABLE]

i.e., equation (24).

Our strategy is analogous to the one described in section 4: it is easy to find 25 eigenvalues and eigenvectors of the system (154) with $\hat{\psi}=\hat{a}=0$. For a more compact notation, we use $\bar{\xi}_{0}$ (see equation (114)) instead of $\xi_{0}$ to list the eigenvalues:

• I.

A 12-dimensional space of modes for arbitrary $\hat{\tilde{\gamma}}_{ij}$ and $\hat{\alpha}=\hat{\omega}=\hat{\beta}^{i}=\hat{\tilde{u}}^{i}=0$, with eigenvalues $\bar{\xi}_{0}=\pm 1$.

• II.

A 2-dimensional space of modes for arbitrary $\hat{\tilde{\gamma}}_{ij}$ and non-trivial $\hat{\alpha}$, $\hat{\omega}$, $\hat{\beta}^{i}$, $\hat{\tilde{u}}^{i}$, with eigenvalues $\bar{\xi}_{0}=\pm 1$.

• III.

A 3-dimensional space of zero speed modes ($\bar{\xi}_{0}=0$) for arbitrary $\hat{\tilde{v}}^{i}$.

• IV.

A 2-dimensional space of modes with eigenvalues $\bar{\xi_{0}}=\pm\sqrt{2\sigma}$.

• V.

A 4-dimensional space of modes with $\bar{\xi}_{0}=\pm\sqrt{f}$, for arbitrary $e^{i}$ orthogonal to $\xi^{i}$ (w.r.t. $\tilde{h}_{ij}$).

• VI.

A 2-dimensional space of modes with $\bar{\xi}_{0}=\pm\sqrt{\frac{4f}{3}}$.

The full expressions of the corresponding eigenvectors are given in Table 3. Clearly, these expressions depend smoothly on $\xi_{i}$ if $f>0$ and $\sigma>0$.

To show that the system consisting of equations (153), (106), (4.1) is strongly hyperbolic, it remains to be shown that (154) has two eigenvectors corresponding to the eigenvalues $\xi_{0}^{\phi,\pm}$ (obtained by solving (24)), with smooth dependence on $\xi_{i}$ and with non-zero $\hat{\psi}$ and $\hat{a}^{\phi}\equiv\frac{i}{N}\xi_{0}^{\phi,\pm}\hat{\psi}$ (see equation (154g)). These eigenvectors can be found as follows. Recall the general form of the characteristic equation (20) from Section 2. When grouping the terms in the characteristic equations (154) to Einstein-scalar-field and Horndeski parts (terms containing a factor $\partial_{X}G_{3}$) as in (2.1), we see that the Horndeski terms only act on the $\hat{\psi}$, $\hat{a}$ components of $U$. In other words, the matrices $\delta\mathbb{A}$ and $\delta\mathbb{L}$ in (2.1) are projections to the subspace associated with the scalar variables:

[TABLE]

The eigenvectors $U^{\phi,\pm}$ corresponding to the eigenvalues $\xi^{\pm,\phi}$ then satisfy

[TABLE]

Thus we have:

[TABLE]

A straightforward but lengthy calculation gives the following result:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the notation $\tilde{T}_{.\xi}$ stands for the contraction $\tilde{T}_{.i}\xi_{j}\tilde{h}^{ij}$ and $I$ is an index for the subspace orthogonal to $\xi^{i}$.

Interestingly, the $\theta$ components of these two eigenvectors are [math]. However, this is not surprising at all: this variable measures the constraint violations but these two eigenvectors correspond to a physical degree of freedom and as such, they must satisfy the high frequency version of the constraints.

In order to avoid singularities in the expressions (158) we need to choose the parameters $\sigma$ and $f$ in such a way that the expressions appearing in the denominators

[TABLE]

[TABLE]

[TABLE]

are non-zero for all $\xi_{i}\in\mathbb{S}^{2}$. In a generic weak field background, the null cone of $P^{\prime}_{\phi\phi}$ is a slightly distorted version of the null cone of the spacetime metric $g$ which means that for specific $\xi_{i}$ they may intersect. Comparing the critical expressions to

[TABLE]

and by similar considerations as in the end of Section 4.2, we see that smooth dependence on $\xi_{i}$ might be violated for some $\xi_{i}\in\mathbb{S}^{2}$ if $\sigma=\frac{1}{2}$, $f=\frac{3}{4}$ and/or $f=1$. This can be easily avoided: choosing suitable $\sigma>\frac{1}{2}$ and $f>1$ resolves this issue, larger values of $\sigma$ and $f$ allows stronger background fields. Therefore, we arrive at a similar conclusion as in the previous section: choosing the gauge parameters $\sigma$ and $f$ large enough ensures strong hyperbolicity of the CCZ4 system, as long as the fields are sufficiently weak. In particular, the combination of the 1+log slicing ($\sigma=\frac{1}{N}$, $N<2$) and a Gamma driver shift condition with $f>1$ appears to be a good candidate for numerical applications.

6 Discussion

In this paper we have provided three well-posed formulations of weak field, cubic Horndeski theories by generalizing an elliptic-hyperbolic, a BSSN-type and the CCZ4 formulation of general relativity. The elliptic-hyperbolic formulation was obtained by generalizing a combination of the constant mean curvature slicing and spatial harmonic gauge condition. In the weak field regime, on spatial slices with negative Ricci curvature, the elliptic equations can be uniquely solved for the lapse function and the shift vector. Under these assumptions, the evolution equations are strongly hyperbolic in any generalized harmonic gauge. The $2$-parameter family of BSSN-type and CCZ4-type formulations are also strongly hyperbolic when the parameters describing the slicing and shift conditions obey suitable bounds. Two important steps in the proofs of these results were noticing that fistly, the gravitational equation of motion (2.1) is quasilinear and secondly, that the second derivatives of the metric disappear in a linear combination of the gravitational and scalar equations of motion.

Analyzing the hyperbolicity of the evolution equations in the elliptic-hyperbolic, BSSN-type or the CCZ4-type formulations of more general Horndeski theories is more difficult. The root of the problem lies in the degeneracy of the principal symbol. In any formulation of GR the principal symbol of the equations of motion has one or more degenerate eigenspaces. This may continue to hold even when small perturbing operators are added to the equations of motion such as Horndeski terms in a weak field regime, especially if the perturbing terms have some special structure. This may lead to the failure of diagonalizability of the principal symbol and hence the failure of strong (or even weak) hyperbolicity. An upcoming paper on this issue is currently under preparation.

Of course, it is possible that more general Horndeski theories do not admit a well-posed initial value formulation at all. However, these theories may still be valid as effective field theories121212In fact, it was shown recently [34], that scalar-tensor EFTs can be written in a Horndeski form up to a fairly high order in the mass dimension expansion since higher derivative operators can be removed by appropriate field redefinitions.. A common feature of the classical equations of motion of EFTs is that they admit runaway solutions. These solutions are typically artifacts of the truncation process and cannot be considered physical. Several techniques have been developed to deal with such pathologies by modifying the equations of motion, some of them are listed in [34]. An example of such methods is the so-called reduction of order process that can be explained as follows. Given an EFT with higher order equations of motion such that the higher derivative terms are subleading in a mass dimension expansion. To any perturbative order, it is possible to derive an alternative equation of motion which is equally accurate up to the given order and contains time derivatives of the fields only up to second order131313In general, it is not possible to get rid of the higher order spatial derivatives.. The advantage of this method is that under certain assumptions (see e.g. [41]) the modified equation no longer has runaway solutions. On the other hand, the procedure breaks Lorentz covariance. Although Horndeski theories are not higher order, the above described process may still be carried out and the “reduced order” equations of motion may have a well-posed initial value problem. Another method is based on the Israel-Stewart process [37] which is widely used in relativistic hydrodynamics to deal with the ill-posedness of the Navier-Stokes equation. Roughly speaking, consider equations of motion of the form $L(\phi)=\epsilon S(\phi)$ for a collection of fields denoted by $\phi$, [35, 36]. Here $L$ is a differential operator such that the zeroth order equation of motion is $L(\phi)=0$; $S(\phi)$ is another differential operator playing the role of a correction term, suppressed with a small parameter $\epsilon$. In this procedure, one introduces an auxiliary variable $\Pi$ for $L(\phi)$ which is forced to satisfy a time evolution equation of the form $\tau\partial_{t}\Pi=-\Pi+\epsilon S$ where $\tau$ is a suitable timescale. The solution to this modified equation asymptotically approaches the solution to the original equation in a timescale $\tau$. This procedure also restricts solutions to the infrared and it may also fix the hyperbolicity of the equations of motion [35, 36].

As a final remark, we point out another (mathematical) issue which may also lead to an interesting research direction. In general, for genuinely higher order EFTs, one might have to rethink what a “well-posed initial value formulation” means. For example, classically, well-posedness is defined in a suitable function space which is typically a Sobolev space. However, for theories that are only meant to be valid up to some finite energy scale, finding a suitable function space for solutions is far from obvious.

Acknowledgements

I would like to thank my supervisor, Harvey Reall, for suggesting the problem, for the many helpful discussions and for the comments on this paper. I am also grateful to Tiago França for carefully reading and commenting on a draft version of this paper, and to Pau Figueras, Joe Keir, Jorge Santos and Ulrich Sperhake for their valuable insights. This research was supported by a George and Lilian Schiff Studentship.

Appendix A Background information

In this section we provide some general information regarding pseudodifferential calculus, elliptic PDEs and hyperbolic PDEs in a self-contained manner.

First, we define the norms and function spaces used in these notes. Sobolev spaces on a manifold $(\Sigma,h)$ will be denoted by $W^{s,p}(\Sigma)$ and $H^{s}(\Sigma)\equiv W^{s,2}(\Sigma)$. The notation $F(\mathbb{R},W^{s,p})$ stands for the space of curves of class $F$ with values in $W^{s,p}$, i.e. if $u(t,x)\in L^{\infty}(\mathbb{R},H^{s})$ then $u(t,\cdot)\in H^{s}$ and $u(\cdot,x)\in L^{\infty}$. When discussing pseudodifferential operators, we will briefly mention Schwartz spaces, denoted by $\mathcal{S}(\mathbb{R}^{n})$. This is the space of smooth functions $f$ satisfying

[TABLE]

for any multi-indices $\alpha$, $\beta$. The Fourier transform of a function or tensor field $A$ will be denoted by $\hat{A}$.

A.1 Basics of pseudodifferential calculus

We briefly summarize the basic results from pseudodifferential calculus, based on [3].

For $m\in\mathbb{R}$, let $S^{m}$ be the space of symbols, that is, the space of smooth functions $p(x,\xi)$ satisfying

[TABLE]

for multi-indices $\alpha$ and $\beta$, some constant $C$ depending only on $\alpha$; $\langle\xi\rangle\equiv(1+|\xi|^{2})^{1/2}$. Given a symbol $p(x,\xi)\in S^{m}$, one can define a corresponding pseudodifferential operator $p(x,\partial_{x})$ as an operator acting on $\mathcal{S}(\mathbb{R}^{n})$ in the following way: for $u\in\mathcal{S}(\mathbb{R}^{n})$

[TABLE]

It is easy to verify that for any $u\in\mathcal{S}(\mathbb{R}^{n})$, $p(x,\partial_{x})u\in\mathcal{S}(\mathbb{R}^{n})$, i.e. $p(x,\partial_{x}):\mathcal{S}(\mathbb{R}^{n})\to\mathcal{S}(\mathbb{R}^{n})$. The space of such operators $p(x,\partial_{x})$ is denoted by $\mathcal{OP}^{m}$. In fact, it can be shown that if $p(x,\partial_{x})\in\mathcal{OP}^{m}$, then $p(x,\partial_{x}):H^{s+m}\to H^{s}$.

It is also possible to define the spaces $S^{\infty}$ and $S^{-\infty}$ by $S^{\infty}\equiv\cup_{m\in\mathbb{R}}S^{m}$ and $S^{-\infty}\equiv\cap_{m\in\mathbb{R}}S^{m}$. It can be shown that pseudodifferential operators with symbols in the class $S^{-\infty}$ are smoothing operators.

An important subspace of $S^{m}$ called classical symbols and denoted by $S_{\text{cl}}^{m}$ is the one consisting of symbols $p(x,\xi)\in S^{m}$ which admit an asymptotic expansion in the following sense. We say that $p(x,\xi)\in S_{\text{cl}}^{m}$ (and $p(x,\partial_{x})\in\mathcal{OP}_{\text{cl}}^{m}$) if there exists a sequence of symbols $p_{m-j}(x,\xi)\in S^{m-j}$, $j=0,1,2,...,\infty$ that are homogeneous in $\xi$ of degree $m-j$ for $|\xi|\geq 1$ such that

[TABLE]

for all $N\geq 1$. The asymptotic expansion is denoted by

[TABLE]

where the notation $\sim$ means equality up to a smoothing operator. Indeed, by definition, the operator

[TABLE]

is a smoothing operator. The first term in the expansion (164) (i.e. the one with $j=0$) is called the principal symbol of $p$ and will be denoted by $p^{(0)}\equiv p_{m}$.

In this paper, we deal with classical pseudodifferential operators. A useful property of operators in this class is that for $p_{1}\in\mathcal{OP}^{m_{1}}_{\text{cl}}$ and $p_{2}\in\mathcal{OP}^{m_{2}}_{\text{cl}}$, the product $q(x,\partial_{x})\equiv p_{1}(x,\partial_{x})p_{2}(x,\partial_{x})$ is also a pseudodifferential operator in $\mathcal{OP}^{m_{1}+m_{2}}_{\text{cl}}$ such that $q^{(0)}(x,\xi)=p_{1}^{(0)}(x,\xi)p_{2}^{(0)}(x,\xi)$, i.e. the principal symbol of the product is the product of the principal symbols.

A.2 Elliptic equations

Now we consider some results regarding elliptic equations based on [3, 16, 17]. An elliptic operator $p(x,\partial_{x})\in\mathcal{OP}^{m}_{\text{cl}}$ is such that $|p(x,\xi)|\geq C\langle\xi\rangle^{m}$ for $|\xi|\geq 1$. Let us first consider a basic result for second order, linear elliptic operators defined by

[TABLE]

with coefficients $a\in H^{s}(\mathbb{R}^{n})$, $b\in H^{s-1}(\mathbb{R}^{n})$ and $c\in H^{s-2}(\mathbb{R}^{n})$, $s>\frac{n}{2}$ and $a^{ij}$ is a positive definite metric. Then

• i)

$P$ is a map $H^{s}\to H^{s-2}$ with finite dimensional kernel ker $P\subset C^{\infty}$.

• ii)

Furthermore, if $c\leq 0$ then ker $P$ is empty and hence $P$ is an isomorphism $H^{s}\to H^{s-2}$.

Part i) of the theorem holds for higher order elliptic operators, whereas ii) generally fails. The operator $q(x,\partial_{x})\in\mathcal{OP}^{-m}_{\text{cl}}$ is said to be a parametrix for $p$ if it satisfies $q(x,\partial_{x})p(x,\partial_{x})\sim p(x,\partial_{x})q(x,\partial_{x})\sim I$ (see notation in previous subsection). Given an equation $p(x,\partial_{x})u=f$ for $u\in H^{s}(\mathbb{R}^{n})$ and $f\in H^{s-m}(\mathbb{R}^{n})$, the formal solution to this equation is $u=q(x,\partial_{x})f~{}\text{mod}~{}C^{\infty}$. In other words, one can say that the solution map $q:f\to u$ is a pseudodifferential operator of class $\mathcal{OP}^{-2}$. For the present purposes, it will suffice to find the principal symbol of the parametrix (or solution map) corresponding to a linear elliptic differential operator of the form (166). Using the product identity, it is easy to see that the principal symbol of the parametrix of $P$ is $Q^{(0)}(x,\xi)=(a^{ij}\xi_{i}\xi_{j})^{-1}$.

A.3 Hyperbolic equations

We conclude this section by discussing the initial value problem for hyperbolic equations [3]. Consider, firstly, a first order linear equation of the form

[TABLE]

where $u$, $g\in C(\mathbb{R},H^{s}(\mathbb{R}^{n}))$ and $f\in H^{s}(\mathbb{R}^{n})$ ($s>\frac{n}{2}+1$) are $N$-component column vectors; $L(t,x,\partial_{x})\in\mathcal{OP}^{1}_{\text{cl}}$ is an $N\times N$ matrix valued function with smooth dependence on $t$, satisfying

[TABLE]

and

[TABLE]

for $\xi\in\mathbb{S}^{n-1}$, with some constant $C>0$ and a positive definite, Hermitian matrix $K(t,x,\xi)$ that depends smoothly on its arguments. Such equations are called strongly hyperbolic and the $N\times N$ matrix $K$ is called a symmetrizer. It can be proved that this problem is locally well-posed in Sobolev spaces $H^{s}$, i.e. there exists a unique solution

[TABLE]

to (167) with $T>0$ depending only on $||f||_{H_{s}}$.

The proof of this result is based on an inequality of the form

[TABLE]

where $\left(\cdot,\cdot\right)_{K}$ is a scalar product constructed in terms of the symmetrizer $K$, that is equivalent to the $L^{2}$-product. This leads to an energy estimate (after an application of the Gronwall inequality) of the form

[TABLE]

The above result can be extended to first order pseudodifferential equations, obtained by a reduction of second (or higher) order equations. This is usually done by introducing $v=\partial_{t}u$ (or the new variable $v$ could also be a linear combination $v\equiv\partial_{t}u-X^{i}\partial_{i}u$ for some $X^{i}$). A second order equation can be rewritten in a form similar to (167) but the operator $L(t,x,\partial_{x},\partial_{x}^{2})$ is now a $2\times 2$ block matrix with $N\times N$ matrix blocks, acting on column vectors $U=(u,v)^{T}\in\mathcal{H}^{s}\equiv H^{s}\times H^{s-1}$. Then the above result holds with initial data $U(0)=F\in\mathcal{H}^{s}$ and the solution will be in

[TABLE]

It can be shown that the existence of a smooth symmetrizer implies that $L(t,x,\xi_{i})$ is diagonalizable with real eigenvalues but it is not always true the other way around. However, if $S$ is the matrix whose columns are the eigenvectors of $L$, i.e. $S^{-1}LS$ is diagonal, then a positive definite symmetrizer is given by $K=\left(S^{-1}\right)^{\dagger}S^{-1}$. One therefore needs to check if the smoothness and boundedness conditions are also met. If the eigenvectors have a smooth dependence on these variables then so does $S$. Since the entries of $S^{-1}$ are rational functions of the entries of $S$, $S^{-1}$ and hence $K$ is also smooth in $(t,x,\xi_{i})$. The boundedness condition follows straightforwardly from smoothness in a compact spacetime. ($\xi_{i}$ takes values in a compact set so no additional assumptions are required for the $\xi_{i}$ dependence). In practice, to demonstrate strong hyperbolicity, it is easiest to show that the matrix $L$ is diagonalizable, has real eigenvalues and the eigenvectors depend smoothly on $(t,x,\xi_{i})$.

The notion of strong hyperbolicity and the above well-posedness result can be extended to quasilinear or even fully non-linear PDEs. More precisely, the following theorems are proved in Chapter 5 of [3].

The initial value problem

[TABLE]

for the quasilinear equation with $u$, $g\in C(\mathbb{R},H^{s}(\mathbb{R}^{n}))$ and $f\in H^{s}(\mathbb{R}^{n})$ ($s>\frac{n}{2}+1$) size-$N$ column vectors; $L(t,x,\partial_{x},u)\in\mathcal{OP}^{1}_{\text{cl}}$ an $N\times N$ matrix valued function that satisfies the above properties of strong (symmetrizable) hyperbolicity for any $(t,x,u,\xi_{i})$, then a unique local solution $u\in C\left([0,T),H^{s}\right)$ exists with $T>0$.

For the nonlinear system

[TABLE]

with data $f\in H^{s}(\mathbb{R}^{n})$ ($s>\frac{n}{2}+2$), [3] shows that if the matrix

[TABLE]

possesses a symmetrizer $K(t,x,u,\partial_{x}u)$ with the above properties, then the Cauchy problem (173) is locally well-posed with a unique solution $u\in C\left([0,T),H^{s}\right)$, $T>0$. These results suggest that to check whether the conditions of strong hyperbolicity are met for quasilinear and nonlinear equations, it is sufficient to study the linearized equations in a generic background.

Appendix B ADM and conformal decomposition

Here we provide a list of formulas used during the ADM decomposition of the equations of motion. The decomposition rules for the derivatives of the scalar field and the curvature tensors:

[TABLE]

where the definitions of the auxiliary variables are

[TABLE]

The conformal versions of these auxiliary variables are given by the following formulas

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, we also provide the conversion rules between curvature tensors, used in the conformal decomposition:

[TABLE]

[TABLE]

[TABLE]