Hyperbolic extensions of constrained PDEs

Systems of partial differential equations (PDEs) comprising a combination of constraints and evolution equations are ubiquitous in physics. For both theoretical and practical reasons, such as numerical integration, it is desirable to have a systematic understanding of the well-posedness of the Cauchy problem for these systems. In this article, we first review the use of hyperbolic reductions, where the evolution equations are singled out for consideration. We then examine in greater detail the extensions, namely, systems in which constraints are evolved as auxiliary variables alongside the original variables, resulting in evolution systems with no constraints. Assuming a particular structure of the original system, we provide sufficient conditions for the strong hyperbolicity of an extension. Finally, this theory is applied to the examples of electromagnetism and a toy model of magnetohydrodynamics.

1 Introduction

In this work, we continue [1–4] the study of first-order systems of equations in which there are more equations than unknowns, but with a structure that permits, in principle, splitting suitable linear combinations of them into “evolution” and “constraint” equations. We restrict to the case of consistent systems, in which the number of equations is equal to the number of constraints plus the number of independent variables, and furthermore to the special case in which the number of independent variables matches the number of evolution equations. The latter means that we do not consider systems with gauge freedom remaining, which would imply the existence of variables with unspecified equations of motion. In this case, one can attempt a solution by carefully restricting the initial data and then directly solving the evolution equations. For an introductory review, see Hilditch [5]. One must then check that the constraint equations are satisfied in the time development. For this, integrability identities among the whole system of equations must be satisfied. These conditions will be assumed and spelled out in detail below. This “free evolution approach” requires us to establish the well-posedness of the Cauchy problem Gustafsson et al. [6]; Kreiss and Ortiz [7] (for a review of well-posedness applied to general relativity, see Sarbach and Tiglio [8]). We restrict ourselves to the concepts arising from the theory of strongly hyperbolic systems, in which well-posedness is determined by algebraic properties of the principal symbol of the equation system. For first-order systems, the principal symbol is simply the set of matrices multiplying the derivatives of the variables. The algebraic properties leading to well-posedness have several equivalent characterizations summarized in the Kreiss matrix theorem Kreiss [9]. To assert well-posedness for the systems under consideration, we need to find a suitable square system, that is, a system where the number of variables equals the number of equations. This can be achieved by taking a subset of the equation system, called a reduction, resulting in a pure evolution system. The use of reductions is customary, but another possibility, which is often employed in numerical schemes, consists of making an extension, that is, extending the system by adding more variables. These extensions are commonly referred to as divergence cleaning [10]; Munz et al. [11, 12], from their use in magnetohydrodynamics, or as $\lambda$ [13] or Z-systems [14] from their use in general relativity.

A paradigmatic example is given by the Maxwell equations,

${\nabla }_a{F}^{ab}=J^b,{\epsilon }^{\mathit{dabc}}{\nabla }_a{F}_{bc}=0,{\nabla }_a{J}^a=0,$

where the unknowns are the components of the Faraday tensor $F_{ab}$, an anti-symmetric tensor (so there are a total of six independent variables). $J^a$, the current vector, is a given vector fixed in space-time, which has vanishing divergence. This is necessary due to the integrability identity ${\nabla }_b({\nabla }_a{F}^{ab})=0$. We work here in four-dimensional space-time $(M,g_{ab})$ with the Levi–Civita derivative ${\nabla }_a$ associated with $g_{ab}$. There are thus a total of $8=4+4$ equations for $F^{ab}$, so six of them should be evolution equations, and the remaining two should be constraints. Introducing a time-like covector $n_a$, one finds that contraction with that vector on both equations gives constraint, that is, equations which have derivatives only in directions perpendicular to $n_a$; while projection on the space perpendicular to $n_a$ gives equations that have derivatives along $n^a$ for each of the independent components of $F^{ab}$. Thus, in the terminology introduced above, a reduction is obtained by taking only these projections as the evolution equations. The integrability identity and divergence property of $J^a$ together imply that constraints are satisfied in the time development if they are at an initial surface.

On the other hand, an extension is given by adding two auxiliary constraint variables $(Z_1,Z_2)$, one for each Maxwell constraint, and making a choice for their equations of motion. To accomplish this in a covariant fashion, we need to define two tensor fields $(g_1,g_2)$. The proposed extended system is

${\nabla }_a{F}^{ab}+g_1^{ba}{\nabla }_a{Z}_1=J^b,{\epsilon }^{\mathit{dabc}}{\nabla }_a{F}_{bc}+g_2^{ba}{\nabla }_a{Z}_2=0,{\nabla }_a{J}^a=0,(1)$

It turns out that if the symmetric parts of $(g_1,g_2)$ are Lorentzian metrics whose cones have non-zero intersections among each other and with the cone of $g$, then the extended system is well-posed. (We use the mathematical notion of a cone; when needed, we use the term light cone to refer to their boundaries). The equations that were constraints are now evolution equations for $(Z_1,Z_2)$, and the others acquire spatial derivatives of these fields. As mentioned above, such extensions have been employed with enormous success in numerical relativity [15–20] and computational astrophysics, with works introducing this approach for magnetohydrodynamics [11, 12]; Dedner et al. [10] is particularly influential. Here, we investigate the space of possible extensions that lead to well-posed Cauchy problems and how to construct them in a natural, covariant fashion.

The article is organized as follows. In Section 2, we define the type of systems to be considered, including the necessary conditions they must satisfy in order to have a well-posed Cauchy problem. In Section 3, we introduce the Kronecker decomposition of matrix pencils and explain its implications to the study of strongly hyperbolic systems. In Section 4, we formalize the framework for extensions. Given the considerable freedom in choosing them, we use the Kronecker decomposition as a guide for making these choices. In Section 5, we demonstrate how this framework applies to two concrete examples: Maxwell’s electrodynamics and a toy model of magnetohydrodynamics (MHD). Finally, in Section 6, we conclude with discussions and provide comments on how this line of research is being further developed.

2 Preliminaries and notation

To fix notation, we specify the systems we consider, following the notation of Geroch [1]; Abalos and Reula [3]; Abalos [4]. We consider a manifold $M$ of dimension $n$, and the following system of the quasi-linear first-order partial differential equations on the fields $\varphi$,

$E^A≔{\mathfrak{N}}_{\alpha }^{Aa}\left(x,\varphi \right){\nabla }_a{\varphi }^{\alpha }-J^A\left(x,\varphi \right)=0,(2)$

where the indices $A$, $a$, $\alpha$ are abstract, grouping several tensorial indices into one and merely indicating where the contractions are. We use lower-case Latin indices to denote single vector indices, lower-case Greek indices to indicate variable fields, and upper-case Latin to label the equations space. The $|\cdot |$ function on indices indicates their total dimension.

We impose the following conditions on ${{\mathfrak{N}}^{Aa}}_{\alpha }\left(x,\varphi \right)$:

Condition 1: the generalized Kreiss condition.

We assume that the matrix ${\mathfrak{N}}_{\alpha }^{Aa}\left(x,\varphi \right)$ is smooth in all arguments and that there exists a hypersurface orthogonal covector $n_a$ such that for all values of $k_a$, not proportional to $n_a$, the matrix pencil

${\mathfrak{N}}_{\alpha }^{Aa}{l}_a\left(\lambda \right)={\mathfrak{N}}_{\alpha }^{Aa}\left(\lambda n_a+k_a\right),$

has a kernel only for a finite set of real values $\{{\lambda }_i(k)\}$ of $\lambda$ (the term matrix pencil refers here to the uni-parametric combination $\lambda \mathfrak{N}+\mathfrak{B}$, where $\mathfrak{N}$ and $\mathfrak{B}$ are matrices that do not depend on the parameter $\lambda$).

In addition, the corresponding singular values of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a(\lambda )$ approach zero in a linear way, that is, $\sigma (\lambda )\ge c_i|(\lambda -{\lambda }_i)|$, with $c_i>0$ in a neighborhood of ${\lambda }_i$. We recall that the singular values are the square roots of the eigenvalues of ${({\mathfrak{N}}_{\alpha }^{Aa}{l}_a)}^T({\mathfrak{N}}_{\beta }^{Ab}{l}_b)$. Because this is an $|\alpha |\times |\alpha |$ matrix, there are $|\alpha |$ singular values (see Abalos [2] for more details and for a more general definition).

These conditions imply two things: i) the rank of ${\mathfrak{N}}_{\alpha }^{Aa}\left(x,\varphi \right)n_a$ is maximal. Therefore, by defining any vector $t^a$ transversal to the surface flat defined by $n_a$ (i.e., $t^a{n}_a\ne 0$), we can obtain all field derivatives along $t^a$ from their values and their derivatives at that surface. This means that we have enough evolution equations for each field ${\varphi }^{\alpha }$. Observe that once we have a choice of $n_a$ satisfying Condition 1, then there is an open set of covectors satisfying the same condition. Thus, we can always form hypersurfaces in a neighborhood of any point, leading to a local initial value problem; ii) In the case that the number of equations equals the number of variables, these conditions imply there is a well-posed Cauchy problem, in the usual sense for strongly hyperbolic systems, off of the mentioned surface. This is the classic Kreiss condition.

In case there are more equations than variables, we need to make sure that there are no more linearly independent equations having derivatives along the transversal vector $t^a$; otherwise, we would have an inconsistency because two equations could give different values for the same transversal derivative. To accomplish that, we impose:

Condition 2: the Geroch constraint condition.

If the number of equations is larger than the number of variables $|A|>|\alpha |$, then we assume there exists a set of matrices $C_A^{\mathrm{\Gamma }a}$, which are labeled by upper-case Greek indices, with

$\begin{array}{l}\hfill {C_A^{\mathrm{\Gamma }(a}{\mathfrak{N}}^{\left|A\right|b)}}_{\alpha }=0,\end{array}$

and that $\text{rank}(C_A^{\mathrm{\Gamma }a}{n}_a)=|A|-|\alpha |=|\mathrm{\Gamma }|$. This condition ensures that the rest of the equations do not have derivatives off of the surface defined by $n_a$, so that the system is consistent. Indeed, the following linear combination of equations, called constraints,

${\psi }^{\mathrm{\Gamma }}≔C_A^{\mathrm{\Gamma }a}{n}_a\left({\mathfrak{N}}_{\alpha }^{Ab}{\nabla }_b{\varphi }^{\alpha }-J^A\right),$

have only derivatives on the flat defined by $n_a$.

There is a further consistency condition that would guarantee that if the initial data are such that constraint quantities vanish at the initial surface, then they would also vanish along evolution [4]. We require the following:

Condition 3: integrability.

${\nabla }_d\left(C_A^{\mathrm{\Gamma }d}{E}^A\right)=L_{1A}^{\mathrm{\Gamma }}\left(x,\varphi ,\nabla \varphi \right)E^A\left(x,\varphi ,\nabla \varphi \right),$

In other words, there is a particular on-shell identity among derivatives of our equation system. In most cases of physical interest, this identity is a consequence of gauge or diffeomorphism invariance.

3 Kronecker decomposition

When studying the well-posedness of the Cauchy problem, the relevant aspect is the behavior of the system in the limit of high frequencies. We can thus restrict our attention to a neighborhood of each point and work in the frequency domain, employing the Fourier–Laplace transform in space and time, respectively. Explicitly, we consider a time function $t$ and a foliation given by its level surfaces. We define $n_a={(dt)}_a$ and take a vector $t^a$ transversal to the foliation and adjust it such that $t^a{n}_a=1$. We choose covectors $k_a$ such that $t^a{k}_a=0$ and define $l_a=\lambda n_a+k_a$. We perform Fourier in $k_a$, and Laplace in $\lambda$. Thus, we replace space derivatives by $i{k}_a$ and time derivatives by $i\lambda$. Furthermore, in what follows, once any particular $k_a$ is chosen, we take a coordinate base so that $n_a={(d{x}^0)}_a$, and $k_a={(d{x}^1)}_a$, and so $l_a=(\lambda n_a+k_a)={(\lambda d{x}^0+d{x}^1)}_a$. Finally, in the high frequency limit, we obtain ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a{\tilde{\varphi }}^{\alpha }=0$.

The Kronecker decomposition of a matrix pencil is a canonical transformation that generalizes the Jordan decomposition of a square matrix pencil. Considering the (square or non-square) pencil $\mathfrak{N}\lambda +\mathfrak{B}$, the Kronecker decomposition is achieved by multiplying this pencil by specific matrices $W$ and $Q$, which are independent of $\lambda$ (as in the square Jordan decomposition case). This transformation results in a new pencil $(W\mathfrak{N}Q)\lambda +(W\mathfrak{B}Q)$ that has a block structure with particular canonical blocks (see Gantmakher [21, 22], for a detailed description and Equation 3 for an example).

It turns out that the Kronecker decomposition can be used naturally in the analysis of systems with constraints or gauge freedom. With the first two conditions assumed above, the Kronecker decomposition of the pencil ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a(\lambda )$ is given by

${\mathfrak{N}}_{\alpha }^{Aa}{l}_a=\left[\begin{array}{ccccccc}\hfill \lambda -{\lambda }_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \lambda -{\lambda }_d\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill \lambda \hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \lambda \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill \hfill & \hfill 0\hfill & \hfill \lambda \hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \dots \hfill & \hfill \hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right](3)$

Ultimately, this represents a change of basis of both the variable and equation spaces, which depends on $k_a$ but not on $\lambda$. The first block is a diagonal $d\times d$ block, this diagonal represents the true degrees of freedom of the entire system. It contains as many elements as the “zeros” of the singular value decomposition, counting their multiplicity. The $2\times 1$ blocks, called $L_1^T$ in the literature, are due to the constraints; there are a total of $r=|\alpha |-d$ blocks. Because each block occupies two rows, we see that the number of zero rows is $s=|A|-d-2r$. The zero rows are present in many systems; they represent differential constraints among the constraints themselves. The numbers defined above also satisfy:

$\begin{array}{ll}\hfill d& :=\dim \left(\text{right}\text{\_}\ker \left(C_A^{\mathrm{\Gamma }a}{n}_a{\mathfrak{N}}_{\alpha }^{Ai}{k}_i\right)\right),\hfill \\ \hfill r& :=\text{rank}\left(C_A^{\mathrm{\Gamma }a}{n}_a{\mathfrak{N}}_{\alpha }^{Ai}{k}_i\right),\hfill \\ \hfill s& =:\dim \left(\text{left}\text{\_}\ker \left(C_A^{\mathrm{\Gamma }a}{n}_a{\mathfrak{N}}_{\alpha }^{Ai}{k}_i\right)\right).\hfill \end{array}$

With this decomposition at hand, it is easy to see how to choose among them linear combinations that give evolution equations for all ${\varphi }^{\alpha }$. Observe that the equations (rows) with a $\lambda$ are certain to contain derivatives transversal to the $n_a$ flats. So, we must include them, but we can add any combination of the other rows to them. It turns out that, by simply adding to each of these rows the immediate row below, multiplied by any number ${\pi }_i,i=1,\dots ,r$, and discarding all the remaining rows, we obtain the evolution equations.

$\begin{array}{l}\hfill h_A^{\beta }{\mathfrak{N}}_{\alpha }^{Aa}{l}_a:=\left[\begin{array}{cccccc}\hfill \lambda -{\lambda }_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill 0\hfill & \hfill ...\hfill & \hfill 0\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \lambda -{\lambda }_d\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \lambda -{\pi }_1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ...\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \lambda -{\pi }_r\hfill \\ \hfill \hfill \end{array}\right].\end{array}$

Thus, we have constructed a map from the equation space to the variable space, which we refer to as a reduction and denote by $h_A^{\beta }$. Thus, $h_A^{\beta }{\mathfrak{N}}_{\alpha }^{Aa}{l}_a$ is a map from the variable space into itself consisting of a set of diagonal matrices satisfying the classic Kreiss conditions (see point ii. within Condition 1). Notice that we can choose the extra roots of $\lambda$ (i.e., the $\{{\pi }_i\}$) as we please. They are the propagation speed of extra constraint modes. This simple observation is the principle behind the results in Reula [23]; Abalos [4].

Thus, there is a reduction (a linear combination of the equations) such that the Cauchy problem of the system is well-posed. Furthermore, Condition 3 asserts that if the initial data satisfy all equations (including the vanishing of the constraints), then all the equations are satisfied for all times as long as the solution exists. See Abalos [4] for details.

4 Extensions

A generic extension would imply the addition of an extra matrix, ${\tilde{\mathfrak{N}}}^{\mathrm{\Delta }Aa}\left(x,\varphi \right)$ (and extra variables $Z_{\mathrm{\Gamma }}$), to obtain a square system

${\mathfrak{N}}_{\alpha }^{Aa}\left(x,\varphi \right){\nabla }_a{\varphi }^{\alpha }+{\tilde{\mathfrak{N}}}^{\mathrm{\Gamma }Aa}\left(x,\varphi \right){\nabla }_a{Z}_{\mathrm{\Gamma }}-J^A\left(x,\varphi \right)+B^A\left(x,\varphi ,Z\right)=0.(4)$

Here, $B^A\left(x,\varphi ,Z\right)$ is a term we can also freely choose that does not include derivatives of $\varphi$ or $Z$ and that goes to 0 when Z goes to 0. In general, $B^A$ represents damping terms [13]; [10]; [24], which are important in numerical applications. For simplicity in our discussion, however, we omit it.

Because we are interested in solving Equation 2 for $\varphi$, our extension proposal only makes sense if we can show that for suitable initial data (for $(\varphi ,Z)$), the solution of Equation 4 has $Z=0$ throughout the development, thereby ensuring that $\varphi$ is a solution of Equation 2.

As we explained before, if we assume Conditions 1, 2, and 3 and take any initial data for $\varphi$ satisfying the constraints, we know that the initial value problem for Equation 2 is “well-posed” and has a unique solution ${\varphi }_{\text{sol}}$. (Here, by well-posed, we mean that the map from Cauchy data to solutions is continuous. To establish this, one finds a hyperbolic reduction from which we may assert that the reduced system is well-posed for arbitrary initial data. Then, one shows that if the initial data satisfy the constraints, then the solutions of the reduced system also satisfy them. Thus, they are solutions to the whole system, and we call the whole system well-posed). Therefore, if we choose ${\tilde{\mathfrak{N}}}^{\mathrm{\Gamma }Aa}$ such that the extended system, Equation 4, is well-posed, then for any initial data, there will be a unique solution. If we choose as initial data $({\left.{\varphi }_{\text{sol}}\right|}_{t=0},{\left.Z\right|}_{t=0}=0)$, then $({\varphi }_{\text{sol}},Z=0)$ will be a solution, and by uniqueness is the solution. Therefore, we only need to show that system Equation 4 satisfies Kreiss’s condition.

4.1 Strong hyperbolicity of the extensions

A particularly interesting set of extensions is obtained by noticing the symmetry between the Kronecker decomposition of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a(\lambda )$ and ${(C_B^{\mathrm{\Delta }a}{l}_a(\lambda ))}^T$. So, we start by computing it:

${\left(C_A^{\mathrm{\Gamma }b}{l}_b\right)}^T=\left[\begin{array}{cccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \lambda \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \dots \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill -1\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \lambda \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \lambda -{\rho }_1\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \dots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \lambda -{\rho }_s\hfill \end{array}\right]$

Recalling that the matrices $C_B^{\mathrm{\Delta }a}{l}_a$ can be thought of as a basis, labeled by $\mathrm{\Delta }$, for the kernel of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$, it is easy to understand its structure. Here, the rows with zeros are $d$ in number. This is so because the diagonal part of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$ cannot contribute to the kernel. We then have $r$ blocks ${\left[\begin{array}{cc}\hfill -1\hfill & \hfill \lambda \hfill \end{array}\right]}^T$, observing that they have a minus sign on them. This is because they are kernels for the corresponding $L_1^T$ blocks of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$. Finally, there is a block that is a kernel of the zero rows of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$. This part is completely undetermined, so we have simply added a diagonal matrix.

To make more apparent the extension we proposed, we reorganize the rows of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$ and ${(C_A^{\mathrm{\Gamma }b}{l}_b)}^T$ such that

${\mathfrak{N}}_{\alpha }^{Aa}{l}_a=\left[\begin{array}{cc}\hfill J\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \lambda I_r\hfill \\ \hfill 0\hfill & \hfill I_r\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right],{\left(C_A^{\mathrm{\Gamma }b}{l}_b\right)}^T=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill -I_r\hfill & \hfill 0\hfill \\ \hfill \lambda I_r\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill J_c\hfill \end{array}\right].(5)$

Here, all the matrices are blocks matrices where $J=(\lambda -{\lambda }_1,\dots ,\lambda -{\lambda }_d)$ of size $d\times d$, $J_c=(\lambda -{\rho }_1,\dots ,\lambda -{\rho }_s)$ of size $s\times s$, and $I_r$ is the identity matrix of size $r\times r$. The zero rows of ${\mathfrak{N}}_{\alpha }^{Aa}{l}_a$ are of size $s\times \left|\alpha \right|$, and the zero rows of ${\left(C_A^{\mathrm{\Gamma }b}{l}_b\right)}^T$ are of $d\times \left|\mathrm{\Gamma }\right|$.

From this reorganization, it is apparent that a natural choice of ${\tilde{\mathfrak{N}}}^{\mathrm{\Gamma }Aa}$ is given by

${\tilde{\mathfrak{N}}}^{\mathrm{\Gamma }Aa}=G^{AB}{C}_B^{\mathrm{\Gamma }a},$

where $G^{AB}$ now must be chosen to render the system diagonalizable. This is, of course, not the most general extension but is a natural and fully covariant proposal for ${\tilde{\mathfrak{N}}}^{\mathrm{\Gamma }Aa}$. The principal symbol of Equation 4 becomes then

$M_D^{Aa}{l}_a=\left[\begin{array}{cc}\hfill {\mathfrak{N}}_{\alpha }^{Aa}\hfill & \hfill G^{AB}{C}_B^{\mathrm{\Delta }a}\hfill \end{array}\right]l_a,$

a $\left|A\right|\times \left|A\right|$ square matrix.

We now propose a particular expression for $G^{AB}$, namely,

$G^{AB}=\left[\begin{array}{cccc}\hfill I_d\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -D^2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill I_r\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_s\hfill \end{array}\right],(6)$

where $D=\text{diag}\left({\pi }_1,\dots ,{\pi }_r\right)$ is of size $r\times r$, and $I_s$ is the identity matrix of size $s\times s$.

Using expressions Equations 5, 6, we conclude

$M_D^{Aa}{l}_a=\left[\begin{array}{cccc}\hfill J\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \lambda I\hfill & \hfill D^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill I\hfill & \hfill \lambda I\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill J_c\hfill \end{array}\right],$

It is easy to verify that this matrix is pencil-similar to the following diagonal matrix:

$M_D^{Aa}{l}_a\sim \text{diag}\left(\dots ,\lambda -{\lambda }_i,\dots ,\lambda +{\pi }_j,\lambda -{\pi }_j,\dots ,\lambda -{\rho }_k,\dots \right)$

and so it satisfies Kreiss’s condition. The extra $2r$ eigenvalues $\{{\pi }_i,-{\pi }_i\}$, introduced by $G^{AB}$, come in pairs, which means that there are $r$ new null cones as characteristic. We shall see this in the examples below, where Lorentzian metrics are used to realize these null cones.

5 Examples

In this section, we present two implementation examples of our proposal, showing that they produce well-posed systems while largely preserving the covariance of the original theories. In all cases, extra Lorentzian metrics are introduced to avoid light cone intersections.

5.1 Maxwell’s equations

We start with the example given in the introduction Equation 1. For them, we have a space of variables $F^{ab}$ (anti-symmetric tensors), which is $|\alpha |=6$ dimensional in a four-dimensional space-time of metric $g_{ab}$. The space of equations is $|A|=8$, namely, two space-time vectors. We have (see Geroch [1])

$\begin{array}{l}\hfill {\mathfrak{N}}_{\alpha }^{Aa}=\left(\begin{array}{c}\hfill {\delta }_{\left[c\right.}^a{\delta }_{\left.d\right]}^q\hfill \\ \hfill {{\epsilon }^{pa}}_{bc}\hfill \end{array}\right)C_A^{b\mathrm{\Gamma }}=\left(\begin{array}{c}\hfill {{\delta }^b}_q\hfill \\ \hfill {{\delta }^b}_p\hfill \end{array}\right)C_A^{b\mathrm{\Gamma }}{l}_b=\left(\begin{array}{c}\hfill l_q\hfill \\ \hfill l_p\hfill \end{array}\right)\end{array}$

Given a time-like $n_a$, we have

$\begin{array}{l}\hfill {\mathfrak{N}}_{\alpha }^{Aa}{n}_a=\left(\begin{array}{c}\hfill n_{\left[c\right.}{\delta }_{\left.d\right]}^q\hfill \\ \hfill {{\epsilon }^{pa}}_{bc}{n}_a\hfill \end{array}\right).\end{array}$

So, it is the map $F_{ab}\to (E_a,B_a)$, which is of the maximal rank. This system satisfies Condition 1; see Abalos and Reula [3] for more details.

The tensor $C_A^{b\mathrm{\Gamma }}{l}_b$ is also of maximal rank for any $l$b¹. Since the dimension of the image is 2-dimensional, we have $|A|=|\alpha |+|\mathrm{\Gamma }|$, and the system is consistent, satisfying Condition 2.

We also have

$\begin{array}{l}\hfill {\nabla }_b(C_A^{b\mathrm{\Gamma }}{\mathfrak{N}}_{\alpha }^{Aa}{\nabla }_a{\varphi }^{\alpha })={\nabla }_b\left(\begin{array}{c}\hfill {\delta }_{\left[c\right.}^a{\delta }_{\left.d\right]}^b{\nabla }_a{F}^{cd}\hfill \\ \hfill {{\epsilon }^{ba}}_{cd}{\nabla }_a{F}^{cd}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\nabla }_b{J}^b\hfill \\ \hfill 0\hfill \end{array}\right)=0\end{array}$

and so Condition 3 is also satisfied.

A suitable reduction is

$\begin{array}{l}\hfill h_{\beta B}=(g_{q\left[r\right.}{t}_{\left.s\right]},-{\frac{3}{2}\epsilon }_{pars}{t}^a).\end{array}$

This renders the evolution equations symmetric hyperbolic. As we saw above, a simple extension is obtained introducing two tensors $(g_1^{pq},g_2^{pq})$ and defining

$G^{AB}=\left(\begin{array}{cc}\hfill g_1^{pq}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill g_2^{pq}\hfill \end{array}\right)$

If we take their symmetric parts to be any two Lorentzian metrics, each one of them sharing a common time-like covector $n_a$ with $g_{ab}$, but not touching their null cones (for brevity, we do not consider here such degenerate cases), then the system is strongly hyperbolic and so has a well-posed Cauchy problem. To check this, we compute the characteristics of the system and the corresponding eigenvectors and see when we get a complete set, that is, a total of eight eigenvectors.

The characteristic equations are

$l_b\delta F^{ab}+g_1^{ab}{l}_b\delta Z_1=0$

${\epsilon }^{\mathit{abcd}}{l}_b\delta F_{cd}+g_2^{ab}{l}_b\delta Z_2=0,$

where we need to solve these equations for $\lambda$ with $l_a=\lambda n_a+k_a$ and $n_a,k_a$ given and for the eigenvectors $\delta F^{ab}$ and $\delta Z_{\mathrm{1,2}}$. The solutions split into three cases: first, when $l_a$ is null with respect to $g^{ab}$ (physical case), then when it is null with respect to $g_1^{ab}$ or $g_2^{ab}$ (extended cases), as we explain below.

We already know four of the eigenvectors, namely, the physical ones arising from the original system. To recover these, we set $\delta Z_1=\delta Z_2=0$ and search for the value of $\delta F_{ab}$. The second equation then implies that $\delta F_{cd}=2l{[}_c{A}_d]$ for some vector $A_d$, while the first implies that $(l_a{l}^a)A^b-(l_a{A}^a)l^b=0$ where indices are raised with the space-time metric. Because $A_a$ cannot be proportional to $l_a$ (otherwise $\delta F_{cd}$ would vanish), both terms must vanish and so we conclude

$g^{ab}{l}_a{l}_b=0,$

which admits two real solutions for $\lambda$. Hence, $A^a$ is orthogonal to $l_a$, which leaves two options remaining for $A^a$ for each of the two values of $\lambda$.

Now, we want to find the rest of the eigenvectors. For that, we first choose $\delta Z_1=1$, $\delta Z_2=0$. Contracting the first equation with $l_b$, and using the anti-symmetry of $\delta F$, we get a condition for $l_a$,

$g_1^{ab}{l}_a{l}_b=0,(7)$

which again admits two real values of $\lambda$. Repeating the argument above, the first equation becomes

$\left(l_a{l}^a\right)A^b-\left(l_a{A}^a\right)l^b+g_1^{ab}{l}_b=0(8)$

Because the null cones of $g^{ab}$ and $g_1^{ab}$ are by assumption not touching, we have $g^{ab}{l}_a{l}_b\ne 0$. It follows that $A^a=-g_1^{ab}{l}_b/(l_c{l}^c)$ satisfies Equation 8 provided that Equation 7 holds. Observe furthermore that $A^a+\alpha l^a$ satisfies the same equations and results in the same Faraday tensor $\delta F_{ab}$ for any $\alpha$. Thus, Equation 7 gives two additional eigenvectors.

If we drop the assumption that the null cones of $g^{ab}$ and $g_1^{ab}$ are non-touching and assume that they touch at $l_a$, then to have a solution, we need that $g^{ab}{l}_b$ must be proportional to $g_1^{ab}{l}_b$.

The final case is similar to the second. We choose $\delta Z_1=0$, $\delta Z_2=1$ and obtain

$g_2^{ab}{l}_a{l}_b=0$

and the same equations for the dual of $\delta F_{ab}$, so we need not discuss it separately.

In summary, we have obtained the eight eigenvectors we require to satisfy the Kreiss condition and conclude that the system is strongly hyperbolic.

5.2 Toy MHD

Here we look at the evolution of a magnetic field $b^a$ driven by a given velocity field $u^a$ in a space-time $(M,g_{ab})$. The system is

$\begin{array}{l}\hfill {\nabla }_a(b^{\left[a\right.}{u}^{\left.b\right]})=0\end{array}(9)$

Here, we take $u^a$ to be time-like and of norm one, $u^a{u}^b{g}_{ab}=-1$. We also take $u^a{b}^b{g}_{ab}=0$. This last is a gauge condition to make the solutions unique for the whole system because otherwise, if $(u^a,b^a)$ is a solution, then $(u^a,b^a+\eta u^a)$ also is a solution, with $\eta$ an arbitrary function.

We observe that there are four equations for three variables. Three of them are evolution equations for the three components of $b^c$. We shall see below that the other is a constraint. Thus, Condition 2 is also satisfied.

The principal part of system Equation 9 is

$\begin{array}{l}\hfill {\mathfrak{N}}_c^{ba}{\nabla }_a{b}^c=u^{\left[a\right.}{\nabla }_a{b}^{\left.b\right]}={\delta }_c^{\left[a\right.}{u}^{\left.b\right]}{\nabla }_a{b}^c.\end{array}$

It is easy to check that Condition 1 is satisfied. The Geroch matrices are also easy to obtain as $C_b^d{l}_d≔{\delta }_b^d{l}_d$. They form a basis of the left kernel of ${\mathfrak{N}}_c^{ba}{l}_a$ and, as we explained before, this means that when Equation 9 is contracted with $C_b^d{u}_d=u_b$, a constraint is generated; this is

${\nabla }_a{b}^a-b^a{a}_a=0,$

where $a^a\equiv u^b{\nabla }_b{u}^a$. We notice that this is the spatial divergence of $b^a$ in disguise.

On the other hand, the following integrability condition $C_b^d{\nabla }_d{\nabla }_a(b^{\left[a\right.}{u}^{\left.b\right]})={\nabla }_b{\nabla }_a(b^{\left[a\right.}{u}^{\left.b\right]})=0$ holds trivially; thus, the system satisfies Condition 3.

The extended system consists of adding a term $g_1^{ba}{\nabla }_{a}Z$ to Equation 9, with $g_1^{ba}$ as in the previous example and with the extra variable $Z$. Its principal part is $u^{\left[a\right.}{\nabla }_a{b}^{\left.b\right]}+g_1^{bc}{C}_c^d{\nabla }_{d}Z=0$, with $C_b^a={\delta }_b^a$. The characteristic equation is

$\frac{1}{2}\left(u^a{l}_a\delta b^b-u^b{l}_a\delta b^a\right)+g_1^{bd}{l}_d\delta Z=0(10)$

where we need to solve this equation for $l_a=-\lambda u_a+k_a$ with $k_a$ given, and for the eigenvectors $\delta Z$ and $\delta b^a$ (with $u_a\delta b^a=0$).

Without loss of generality, we choose $k^a$ such that $u^a{k}_a=0$, and we rewrite the characteristic equations projecting on to $u_a$ and perpendicular to it (with the projector $h_{ab}\equiv g_{ab}+u_a{u}_b$). We obtain

$\begin{array}{ll}\hfill \frac{1}{2}{k}_a\delta b^a+u_a{g}_1^{ab}{l}_b\delta Z& =0\hfill \\ \hfill \frac{1}{2}\lambda \delta b^a+h_c^a{g}_1^{cb}{l}_b\delta Z& =0\hfill \end{array}$

The physical solution comes from choosing $\lambda =0$, and the eigenvectors $\delta Z=0$ and $\delta b^a$ orthogonal to $k_a$. Because $\delta b^a$ has two possible directions, we obtain two eigenvectors.

The remaining eigenvectors come from choosing $\lambda$ such that

$l_a{g}_1^{ab}{l}_b=0,(11)$

and $\delta Z=\frac{1}{2}\lambda$, $\delta b^a=h_c^a{g}_1^{cb}{l}_b$. This expression satisfies the second characteristic equation trivially, and it is easy to verify that the first one reduces to

$\frac{1}{2}{k}_a\delta b^a+u_a{g}_1^{ab}{l}_b\delta Z=\frac{1}{2}{l}_a{g}_1^{ab}{l}_b=0.$

Because, as before, there are two solutions for $\lambda$ from Equation 11, we obtain two more eigenvectors. In summary, we have obtained the four eigenvectors we require to satisfy the Kreiss condition and conclude that the extended system is strongly hyperbolic. Finally, we notice that Equation 11 can also be rederived from the integrability condition, i.e., by multiplying Equation 10 by $C_b^d{l}_d=l_b$.

6 Conclusion

Similar extensions to those proposed here were previously known, starting with the divergence cleaning used in magnetohydrodynamics and later generalized as $\lambda$-systems for generic symmetric hyperbolic systems. To implement them, it was necessary to break the covariance of the system in the usual sense of performing a $3+1$ decomposition. For symmetric hyperbolic systems, such extensions can be obtained in our framework by committing to a frame and a reduction and adding an extra term that annihilates the time component of the constraint basis. This results in an extended symmetric hyperbolic system.

In this article, we have presented an extension scheme for first-order PDEs. With appropriate adaptation, however, these results can be applied to systems of two or even more orders. We will show in future articles how to apply these ideas to gravity theories to extend the system and to fix the gauge, allowing us to reinterpret and generalize known results such as those of Bona et al. [25]; Hilditch and Richter [26]; Kovács and Reall [27].

Although the existence of a strongly hyperbolic extension is performed in Fourier space and results in a system of pseudodifferential equations, our examples show that in cases of physical interest, one may obtain differential extensions. These extensions furthermore retain covariance of the theory in the sense that, contrary to earlier $\lambda$-system extensions, at least in the principal part, they do not rely on a preferred time direction but instead the addition of other Lorentzian metric tensors. Further details and a complete proof will be provided in a longer version of this work.

In our analysis, we resorted to previous work to argue that the constraints, if initially satisfied, are satisfied at later times. This helped us conclude that $Z_{\mathrm{\Gamma }}$ remains zero throughout the evolution. There are, however, more elegant ways to show this when the constraints do not have any kernel from the left, that is, no set of zero rows in their Kronecker decomposition (see Equation 3). In such cases, it can be shown that the $Z_{\mathrm{\Gamma }}$ fields satisfy a second-order evolution system that is decoupled from ${\varphi }^{\alpha }$ and has a well-posed initial value problem. Choosing these fields to vanish at the initial surface and the ${\varphi }^{\alpha }$ fields satisfy the original constraints of the system, all derivatives of $Z_{\mathrm{\Gamma }}$ vanish on the initial surface, in particular any transversal derivative, so the unique solution to the second-order system is 0, and the constraints are satisfied for all times. Unfortunately, the presence of zeros may prevent the second-order system from being well-posed, so more care is needed. This will be further considered in the aforementioned longer article.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

FA: conceptualization, formal analysis, investigation, methodology, project administration, validation, writing–original draft, and writing–review and editing. OR: conceptualization, formal analysis, methodology, validation, writing–original draft, and writing–review and editing. DH: investigation, validation, writing–review and editing, and formal analysis.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Grant PID2022-138963NB-I00 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE, by CONICET, SeCyT-UNC, and MinCyT-Argentina and by FCT Project No. UIDB/00099/2020.

Acknowledgments

We thank Carlos Palenzuela for several helpful discussions associated with the article. We also thank the “Programa de projectes de recerca amb investigadors convidats de prestigi reconegut” of the Universitat de les Illes Balears through which Oscar Reula visited the UIB, making it possible to finalize this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Here the target space is two copies of $R^4$, and the image is 1-dimensional on each one of them.

References

1. Geroch R. Partial differential equations of physics. Scotland: General Relativity, Aberdeen (1996). p. 19–60.

Google Scholar

2. Abalos F. A necessary condition ensuring the strong hyperbolicity of first-order systems. J Hyperbolic Differential Equations (2019) 16:193–221. doi:10.1142/s0219891619500073

CrossRef Full Text | Google Scholar

3. Abalos F, Reula O. On necessary and sufficient conditions for strong hyperbolicity in systems with constraints. Class Quant Grav (2020) 37:185012. doi:10.1088/1361-6382/ab954c

CrossRef Full Text | Google Scholar

4. Abalos JF. On constraint preservation and strong hyperbolicity. Classical Quan Gravity (2022) 39:215004. doi:10.1088/1361-6382/ac88af

CrossRef Full Text | Google Scholar

5. Hilditch D. An introduction to well-posedness and free-evolution. Int J Mod Phys A (2013) 28:1340015. doi:10.1142/S0217751X13400150

CrossRef Full Text | Google Scholar

6. Gustafsson B, Kreiss H-O, Oliger J. Time dependent problems and difference methods, 24. John Wiley and Sons (1995). doi:10.1002/9781118548448

CrossRef Full Text | Google Scholar

7. Kreiss H-O, Ortiz OE. Introduction to numerical methods for time dependent differential equations. John Wiley and Sons (2014).

Google Scholar

8. Sarbach O, Tiglio M. Continuum and discrete initial-boundary value problems and einstein’s field equations. Living Rev relativity (2012) 15:9. doi:10.12942/lrr-2012-9

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kreiss H-O. Über die stabilitätsdefinition für differenzengleichungen die partielle differentialgleichungen approximieren. BIT Numer Mathematics (1962) 2:153–81. doi:10.1007/BF01957330

CrossRef Full Text | Google Scholar

10. Dedner A, Kemm F, Kröner D, Munz C-D, Schnitzer T, Wesenberg M. Hyperbolic divergence cleaning for the mhd equations. J Comput Phys (2002) 175:645–73. doi:10.1006/jcph.2001.6961

CrossRef Full Text | Google Scholar

11. Munz C-D, Omnes P, Schneider R, Sonnendrücker E, Voß U. Divergence correction techniques for maxwell solvers based on a hyperbolic model. J Comput Phys (2000) 161:484–511. doi:10.1006/jcph.2000.6507

CrossRef Full Text | Google Scholar

12. Munz C-D, Ommes P, Schneider R. A three-dimensional finite-volume solver for the maxwell equations with divergence cleaning on unstructured meshes. Computer Phys Commun (2000) 130:83–117. doi:10.1016/s0010-4655(00)00045-x

CrossRef Full Text | Google Scholar

13. Brodbeck O, Frittelli S, Hübner P, Reula OA. Einstein’s equations with asymptotically stable constraint propagation. J Math Phys (1999) 40:909–23. doi:10.1063/1.532694

CrossRef Full Text | Google Scholar

14. Bona C, Ledvinka T, Palenzuela C, Zacek M. General covariant constraint free evolution system for numerical relativity. arXiv (2002). doi:10.48550/arXiv.gr-qc/0209082

CrossRef Full Text | Google Scholar

15. Baumgarte TW, Shapiro SL. Numerical integration of Einstein’s field equations. Phys Rev D (1998) 59:024007. doi:10.1103/PhysRevD.59.024007

CrossRef Full Text | Google Scholar

16. Shibata M, Nakamura T. Evolution of three-dimensional gravitational waves: harmonic slicing case. Phys Rev D (1995) 52:5428–44. doi:10.1103/PhysRevD.52.5428

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Nakamura T, Oohara K, Kojima Y. General relativistic collapse to black holes and gravitational waves from black holes. Prog Theor Phys Suppl (1987) 90:1–218. doi:10.1143/ptps.90.1

CrossRef Full Text | Google Scholar

18. Pretorius F. Numerical relativity using a generalized harmonic decomposition. Class Quant Grav (2005) 22:425–51. doi:10.1088/0264-9381/22/2/014

CrossRef Full Text | Google Scholar

19. Bernuzzi S, Hilditch D. Constraint violation in free evolution schemes: comparing the BSSNOK formulation with a conformal decomposition of the Z4 formulation. Phys Rev D (2010) 81:084003. doi:10.1103/PhysRevD.81.084003

CrossRef Full Text | Google Scholar

20. Alic D, Bona-Casas C, Bona C, Rezzolla L, Palenzuela C. Conformal and covariant formulation of the z4 system with constraint-violation damping. Phys Rev D (2012) 85:064040. doi:10.1103/physrevd.85.064040

CrossRef Full Text | Google Scholar

21. Gantmakher FR. The theory of matrices, 1. New York: Chelsea Publishing Company (1998).

Google Scholar

22. Gantmakher FR. The theory of matrices, 2. Providence, RI: American Mathematical Soc. (1998). p. 131.

Google Scholar

23. Reula OA. Strongly hyperbolic systems in general relativity. J Hyperbolic Differential Equations (2004) 1:251–69. doi:10.1142/s0219891604000111

CrossRef Full Text | Google Scholar

24. Gundlach C, Calabrese G, Hinder I, Martín-García JM. Constraint damping in the z4 formulation and harmonic gauge. Classical Quan Gravity (2005) 22:3767–73. doi:10.1088/0264-9381/22/17/025

CrossRef Full Text | Google Scholar

25. Bona C, Ledvinka T, Palenzuela-Luque C, Zacek M. Constraint-preserving boundary conditions in the Z4 numerical relativity formalism. Class Quant Grav (2005) 22:2615–33. doi:10.1088/0264-9381/22/13/007

CrossRef Full Text | Google Scholar

26. Hilditch D, Richter R. Hyperbolicity of physical theories with application to general relativity. Phys Rev D (2016) 94:044028. doi:10.1103/PhysRevD.94.044028

CrossRef Full Text | Google Scholar

27. Kovács AD, Reall HS. Well-posed formulation of lovelock and horndeski theories. Phys Rev D (2020) 101:124003. doi:10.1103/PhysRevD.101.124003

CrossRef Full Text | Google Scholar