The behavior of solutions to quasi–linear equations is not yet fully understood. Most of the solutions develop singularities in a finite time for most initial data, even if they are in
${M}_{0}$
. This is the case for convective systems, or more generally for genuinely nonlinear systems, see [
43,
53]
, for the definition and main results, a class which includes systems like perfect fluids and relativistic dissipative fluids –for they contain as part of the system the perfect fluid equations. This is also the case for general relativity, where singularity theorems (see [
38,
56]
) tell us about the development of singularities, although of a different type. Thus the concept of well posedness has to be modified to account for the fact that solutions only last for a finite time and this time depends on the initial data. Basically, the most we can pretend to show in the above generality is the same type of well posedness one requires from an ordinary system of equations. Which is quite a lot! The nonlinear aspect of the equations implies also that it is not possible to generalize their solutions to be distributions. The minimum differentiability needed to make sense of an equation depends on the particular equation. Furthermore, there are cases (e.g. convection) in which, for some function spaces of low differentiability, the equation makes sense and some solutions exist, but they are not unique^{
$\text{3}$
}
.
Definition 4
Let
${u}_{0}(t,x),t\in [0,{T}_{0})$
,
${T}_{0}<+\infty $
be a smooth solution of a quasi–linear evolution system. We shall say the system is well posed at the solution
${u}_{0}$
with respect to a norm
$\left\right\left\right$
if given any
$\delta >0$
there exists
$\varepsilon >0$
such that for any smooth initial data
$f\left(x\right)$
such that
$\left\rightf{f}_{0}\left\right<\varepsilon $
, with
${f}_{0}\left(x\right):={u}_{0}(0,x)$
there exists a smooth solution
$u(t,x)$
defined in a strip
$0\le t<T$
, with
$\leftu\right(t,\cdot ){u}_{0}(t,\cdot \left)\right<\delta $
,
$T{T}_{0}<\delta $
.
In order not to worry about the possibility that the smoothness of the solutions be too stringent a requirement, one can smooth out the equation using a one parameter family of mollifiers, and require that the relation
$\delta \left(\varepsilon \right)$
be independent of that parameter family.
To obtain results about well posedness, we just have to slightly modify the concepts of hyperbolicity already discussed in the constant coefficient case. Since in the constant coefficient case the matrices did not depend on the points of spacetime, nor on the solution itself, we had only two cases. In one case, the norm
$H$
did not depend on
${\omega}_{a}$
, and so in some base the matrix
${A}^{a}$
was symmetric. In the other case, the norm
$H$
did depend on
${\omega}_{a}$
, and we had a general strongly hyperbolic system.
In the latter case, it can be seen that
$H\left(\omega \right)$
is piece–wise continuous and so integrable, which is, in that case, all that is needed to proceed with the proof. In the general case with which we are now dealing,
$H$
would in general depend not only on
${\omega}_{a}$
, but also on the point of spacetime and on the solution,
$H=H({\omega}_{a},t,x,u)$
.
This difference has caused terminology to be not uniform in the literature, so I have taken advantage of this and establish terms in the way I consider best suited for the topic.
Certain authors call some systems symmetric hyperbolic and others symmetrizable. They call symmetric hyperbolic only those systems where the symmetrizer does not depend on the unknown variables nor on the spacetime variables (or at most depends only on the base space variables
$H:=H(t,x)$
); they call the other systems symmetrizable. This is a rather arbitrary distinction, since the methods of proof used are valid for both with no essential difference. Thus, if
$H$
does not depend on
${\omega}_{a}$
but depends smoothly on all other variables,
$H:=H(t,x,u)$
, then we shall still say the system is symmetric hyperbolic. In this case the nonsingular transformation which symmetrizes
${A}^{a}(t,x,u)$
is smooth in all its variables.
The existence and smoothness proof is based, as in the constant coefficient case, on energy norm estimates, but now supplemented by Sobolev inequalities. Since the norm is built out of
$H$
and it does not depend on
${\omega}_{a}$
, no passage to Fourier space is needed.
If
$H$
does also depend on
${\omega}_{a}$
, and is smooth on all variables,
$H:=H(t,x,u,{\omega}_{a})$
, we shall say the system is strongly hyperbolic. The existence and smoothness proof now requires the construction of a pseudodifferential norm out of
$H$
, and so pseudodifferential calculus is needed, which implies that
$H$
has to be smooth in all its entries, in particular in
${\omega}_{a}$
.
We shall not discuss weak hyperbolic systems, for they are generically unstable under perturbations, nor shall we discuss strictly hyperbolic systems, i.e. systems with strictly different eigenvalues of
${A}^{a}{\omega}_{a}$
, for they seldom appear in physical processes in more than one dimension.
With this concept of well posedness we have the following theorem [See for instance [
54]
pg.
123]:
Theorem 4
Let
$r>\frac{n}{2}+1$
, then a strongly hyperbolic system is well posed with respect to the Sobolev norm
$\left\right{}_{r}$
. The solution is in
$C\left(\right[0,T),{H}^{r})$
, the time of existence depends only on
$\left\rightf{}_{r}$
.
Remarks:

∙
In the generic case the value of
$r$
cannot be reduced from the above value, but of course it can for certain special types of systems. In general relativity, a slight improvement, (
$r>\frac{n}{2}$
), is obtained from the fact that the matrix
${A}^{a}$
only depends on a subset of variables (the metric).

∙
The time of existence is the same for all
$r$
. That is, when a solution loses differentiability, it loses all of it at the same time. This is reflected in the following result: [54] . If a solution exists until
${T}_{0}$
and there
$\left\rightu(t,\cdot ){}_{{C}_{\u25c6}^{1}}$
is bounded, (where
${C}_{\u25c6}^{1}$
is a Zygmund space, see Appendix A in [54] ), then the solution can be extended further. Again in general relativity a slight improvement can be obtained, see [39] , and [17] .
2.4 Hyperbolicity and Numerical Simulations.
Knowing that the field equations for GR can be cast in a symmetric hyperbolic form, we can now ask how this fact can be of help for numerical calculations, (besides the extremely important fact that the problem would be well posed and so tractable!). There are at least two reasons why one should use variables in which the system is hyperbolic when performing numerical simulations.
The first reason is that having a strongly hyperbolic system allows for standard constructions of numerical codes which are stable. If the system is symmetric, hyperbolic codes with better properties can be constructed. In particular, schema for symmetric hyperbolic systems can be constructed with numerical dissipation of lower order than what is needed for generic strongly hyperbolic ones (see chapter VI of [
36]
). In variables where the system becomes diagonal, one can also use methods which take advantage of that structure.
The second reason for using a hyperbolic formulation for numerical analysis is that well posedness of the system gives bounds on the growth of the solution and its derivatives, as long as the solution is smooth
^{
$\text{4}$
}
. This property, when used in conjuction with stable algorithms, implies that one can bound the errors made on the simulation. That is, one knows not only that the error goes to zero as some power of the step size, but also the proportionality factor of that power law. In simulating phenomena whose observation would require hundreds of millions of dollars, a tight control of the accuracy reached should be required. Nevertheless it should be noted that raw hyperbolicity estimates alone usually give exponential bounds with very large growth coefficients, and that they are not of much value for numerical work.
Many of the systems we shall analyze can be cast as flux conservative equations with sources.
This is a direct consequence of the facts that the principal part of the equations depends only on the metric variable, and that the equation for the time derivative of it does not contain derivatives of any of the dynamical variables. This property is important when using codes with variable grid spacing, even more if one considers that there are many standard codes for fluids –which are truly flux conserving– with adaptive grid schema.
It has to be said that flux conservation is important when dealing with systems that develop shock waves, that is in convective or more precisely in genuinely non–linear systems (for a definition of this term and many of the results, see [
43,
53]
). One should be cautious about any expectation of improvement by using flux conservative properties in general relativity, since here the shocks would probably not develop –in particular the systems are not genuinely nonlinear. Rather, when singularities appear, they would be much worse than mere discontinuities of some of the dynamical fields.^{
$\text{5}$
}
Due to bad gauge choices, discontinuities resembling shocks have been observed in numerical simulations, see [
5]
. Perhaps, instead of trying to devise an algorithm which allows one to go through these discontinuities, one should concentrate on finding better gauges, where it could even happen that the system cannot be put in flux conserving form. Thus, it is not clear whether flux conservative forms are relevant for vacuum general relativity.^{
$\text{6}$
}
Going Further For the reader wishing to delve further into this precious theory of hyperbolic systems, while keeping a physicist's approach, I recommend [
34]
. For those wishing to see more of the machinery at work, I recommend the book of Kreiss and Lorentz [
48]
. Finally, for those who really want to get the latest on the technical aspects and the modern approach to the problem, I recommend Taylor's book, [
54]
. Considerations about numerical analysis and algorithms can be found in [
36]
.
In particular, that book contains general stable algorithms for strongly and symmetric hyperbolic systems and numerical error bounds in terms of analytic bounds of the exact solutions applicable to non–linear systems.
3 The Problem of Hyperbolicity in General Relativity
What follows are descriptions of the problem of hyperbolicity in general relativity and of the main approaches that have been proposed to deal with the problem.
In studying Einstein's field equations we are faced with a problem (see [
20,
34,
29]
): While the theory of partial differential equations has developed as a theory for tensor components in a given coordinate system, or at best for tensor fields in a given metric space, Einstein's equations acquire their full meaning –and the characteristics which distinguish them from all other physical theories– when they are viewed as equations for geometries, that is, for equivalent classes of metric tensors.
The object of the theory is not a metric tensor, but the whole equivalence class to which it belongs –all other metrics related to the first one by a smooth diffeomorphism. This fact is contained in the equations, for they are invariant under those diffeomorphisms. To clarify the concept, and see the problem, let us assume we have a solution to Einstein equations in a given region of a manifold.
Take a space–like hypersurface across it,
${\Sigma}_{0}$
, and a small “lens shaped” region which can be foliated by smooth spacelike surfaces
${\Sigma}_{t}$
starting at
${\Sigma}_{0}$
. If Einstein's equations were hyperbolic for the metric tensor, then uniqueness of the solution (the metric tensor) in the lens shaped region would follow from the standard theory once proper initial data at
${\Sigma}_{0}$
is given. But we know that if we apply a diffeomorphism to the original metric tensor solution, which is different from the identity only in a region inside the lens shaped one but which does not intersect the initial
${\Sigma}_{0}$
slice, the resulting metric would also satisfy Einstein's equations, thus contradicting uniqueness, and so the possibility that the system be hyperbolic.
Since, as shown in
§ 2 , hyperbolicity is equivalent to the existence of norms which are bounded under evolution, we see that for Einstein's equations there cannot be such norms on the space of metric tensors. Norms are not only important for well posedness, but also for other related issues which often appear in general relativity, in particular when one tries to see whether some approximation schema is indeed an approximation. Examples of this appear in very unrelated cases, for instance, in numerical algorithms and postNewtonian approximations. Thus, a method is needed to find relevant norms on metric tensors, that is, to break the diffeomorphism invariance.
The norms thus obtained are not natural, and so by themselves do not imply any physical closeness of metrics in numerical values. They have to be considered in their topologically equivalent class.
Physically relevant notions of closeness can still be obtained by building, out of the metric tensor and its derivatives, diffeomorphism invariant quantities and making the comparisons with then.
Can we avoid this detour into tensors and make a theory of diffeomorphism invariant objects?
It is not clear whether this can be done. Some attempts in this direction have been made by trying to build norms which have some partial diffeomorphism invariance. Here the norms are made out of scalars built out of curvature tensor components of the metric, in particular see [
20,
29]
. But I think a fully geometrical theory needs other types of mathematics than the theory of partial differential equations, a theory which might be emerging from parallel questions in quantum gravity.
3.1 The Standard Approach, or the 4D Covariant Approach.
The standard approach to overcome this problem has been to “fix the gauge”, that is, by imposing some extra condition on the metric coordinate components which would select one and only one representative from each equivalent class of Einstein's geometries. With a clever choice of gauge fixing, commonly the socalled harmonic gauge, Einstein's equations can be “reduced” to a hyperbolic system by removing from them the parts which the gauge condition would make vanish.
This reduced system is equivalent to the full Einstein equations if one can prove that, possibly after imposing further initial data set constraints, the solutions of the reduced system satisfy the gauge conditions previously imposed to the system. Alternatively, one can think of this method as fixing some components of the tensor obtained by taking the difference between two connections, the one associated with the metric tensor at which we are looking and the one associated with some other arbitrary background metric [
23,
38]
.
A convenient way to describe this scheme is by introducing a background metric,
${\stackrel{~}{g}}^{ab}$
, thus the gauge is not a coordinate condition, but rather a condition which links the physical metric with the background one. In this approach, see [
38]
, the basic variable is a densitized symmetric tensor,
${\Phi}^{ab}:=g\delta {g}^{ab}$
, where
$g$
is the metric determinant with respect to the background one, ^{
$\text{7}$
}
and
$\delta {g}^{ab}:={g}^{ab}{\stackrel{~}{g}}^{ab}$
.
In these variables Einstein's equations become,
$$\begin{array}{ccc}& & \frac{1}{2}{g}^{cd}{\stackrel{~}{\nabla}}_{c}{\stackrel{~}{\nabla}}_{d}{\Phi}^{ab}{g}^{c(a}{\stackrel{~}{\nabla}}_{c}{\Psi}^{b)}+\frac{1}{2}{g}^{ab}{\stackrel{~}{\nabla}}_{c}{\Psi}^{c}\end{array}$$ 
(2)

$$\begin{array}{ccc}& & +\left(termsin{\stackrel{~}{\nabla}}_{c}\delta {g}^{de}and\delta {g}^{de}\right)\end{array}$$ 
(3)

$$\begin{array}{ccc}& & =8\pi g{T}^{ab}g{\stackrel{~}{G}}^{ab},\end{array}$$ 
(4)

where
${\stackrel{~}{\nabla}}_{c}$
is the covariant derivative associated with
${\stackrel{~}{g}}^{ab}$
,
${\stackrel{~}{G}}^{ab}$
its Einstein tensor, and
${\Psi}^{a}:={\stackrel{~}{\nabla}}_{b}{\Phi}^{ab}={\stackrel{~}{\nabla}}_{b}\sqrt{\frac{g}{\hat{g}}}{g}^{ab}$
.
In that gauge, Einstein's equations are “reduced” to a hyperbolic system by removing from them all terms containing
${\Psi}^{a}$
, for this quantity is assumed to vanish in this gauge. By doing this one gets a set of coupled wave equations, one for each metric component. Thus by prescribing at an initial (spacelike) hypersurface values for
${\Phi}^{ab}$
and of its normal derivatives one gets unique solutions to the reduced system. When are such solutions solutions to Einstein's equations? That is, under what conditions does
${\Psi}^{a}$
vanish everywhere? It turns out that the Bianchi identity grants that when
${\Phi}^{ab}$
satisfies the reduced equations, then
${\Psi}^{a}$
satisfies a linear homogeneous second order hyperbolic equation. Standard uniqueness results for such systems implies that if initially
${\Phi}^{ab}$
is chosen so that
${\Psi}^{a}$
and its normal derivative vanish at the initial surface, then they vanish everywhere on the domain of dependence of that surface. Thus the question is now posed on the initial data, that is, on whether it is possible to choose appropriate initial data for the reduced system,
$({\Phi}^{ab},{\dot{\Phi}}^{ab})$
, in such a way that
$({\Psi}^{a},{\dot{\Psi}}^{a})$
vanish initially. It turns out, using that the reduced equations are satisfied at the initial surface, that one can indeed express
${\Psi}^{a}$
and its normal derivative at the initial surface, in terms of
${\Phi}^{ab}$
and its derivatives (both normal and tangential to the initial surface). Thus one finds there are plenty of initial data sets for which solutions to the reduced system coincide with solutions to the full Einstein system. Are they all possible solutions to Einstein's equations, or are we loosing some of them by imposing this scheme? The answer to the first part of the question is affirmative (subject so some asymptotic and smoothness conditions), for one can prove that given “any” solution to the Einstein equations, there exists a diffeomorphism which makes it satisfy the above harmonic gauge condition.
It is important to realize that it is not necessary to set
${\Psi}^{a}$
to zero to render the Einstein equations hyperbolic; it just suffices to set it equal to some given vector field on the manifold, or any given vector function of the spacetime points and on the metric, but not its derivatives. So there are actually many ways to hyperbolize Einsteins's equations via the above scheme. We shall call all of them harmonic gauge conditions, and reserve the name full harmonic condition to the one where
${\Psi}^{a}\equiv 0$
.
An important advantage of this method is that some gauge conditions, like the full harmonic gauge, are fourdimensional covariant –although a background metric is fixed– a condition which can be very useful for some considerations.
One drawback of this method, at least in the simplest version of the harmonic gauge, i.e. the full harmonic gauge, was recognized early, [
14]
. The drawback is the fact that this gauge condition can be imposed only locally, and generically breaks down in a finite evolution time. A related problem has been discussed recently in [
5]
in the context of the hyperbolizations of the ADM variables with the harmonic gauge along the temporal direction. The above disadvantage can be considered just a manifestation of another: the lack of ductility of the method, that is the fact that one has been able do very little besides imposing the full harmonic gauge condition, and that for each new harmonic gauge condition one would like to use, a whole study of the properties of the reduced equations would have to be undertaken. Although there are many other gauge conditions besides the harmonic one, the issue of the possibility of their global validity, or the search for other properties of potential use, do not seem to have been considered. For a detailed discussion of this topic, see [
28,
29]
, and also [
38]
.
One can summarize the situation by noticing that in this setting one needs to prescribe a four vector as a harmonic gauge condition. Since the theory keeps its four dimensional covariance, then it is hard to choose any other vector but zero, that is the full harmonic gauge. Since recently there have been no advances in this area, I do not elaborate on it.
3.2 The Modification of the Field Equations Outside the Constraint Submanifold, or the 3+1 Decomposition Point of View.
Another approach to deal with the diffeomorphism freedom of Einstein's equations is by first removing the diffeomorphism invariance. This is done by prescribing the time foliation along evolution, that is, by prescribing a lapseshift pair along evolution. This removes the diffeomorphism invariance up to threedimensional diffeomorphisms at the initial surface. Sometimes four dimensional covariance is also broken by splitting Einstein's equations, and possibly other supplementary equations, with respect to that foliation, and then recombining the split pieces in a suitable way. The resulting equations are equivalent to the original ones –for it is just a linear combination of them– so they have the same solutions. Notice that here we are doing more than just a
$3+1$
decomposition, since in general one is recombining spacespace components of the Einstein tensor with timetime, and timespace components in a noncovariant way, and taking as equations this combination, or even transforming the equations to first order in derivatives of the variables by defining new variables and equations and modifying them.
After this procedure is done, one obtains a system which is symmetric hyperbolic for most choices of
given lapseshift functions, once they are suitably rescaled. Subsequent arguments go very much on adding equations for the lapseshift vector in order to make the whole system well posed, and presumably useful for some application. It is instructive to think of these modifications of the evolution equations from the point of view of the initial value formulation. There one starts by solving the constraint equations, the timetime and timespace components of the Einstein tensor, at the initial surface. With the initial data thus obtained, one finds the solution to the evolution equations which are taken to be the spacespace components of the Einstein tensor. Since that evolution preserves the constraints, (The vector field generating the flow in phase–space is tangent to the constraint submanifold.), one can forget about the constraint equations and think of the evolution equations as providing an evolution for the whole phase–space. In this sense, the modification one is making affects the evolution vector outside the constraint submanifold, leaving the vector intact at it. Uniqueness of solutions, which follows from the well posedness of the system, then implies that the solutions stay on the submanifold. Nevertheless, and we shall return to this point, as shown in [
30]
, there is no guarantee that the submanifold of constraint solutions is stable with respect to the evolution vector field as extended on the whole phase–space.
This is an important point for numerical simulations.
4 Recent Approaches to the Problem.
4.1 The ADM representation.
In the ADM formulation of Einstein's equations, the fields involved are detached from the underlaying spacetime and brought into an abstract three dimensional manifold product a segment of the real line. To obtain the relation between these abstract fields and the metric tensor defining a solution of Einstein's equations, we start by pretending we have such a solution, that is a spacetime
$(M,{g}_{ab})$
. In this spacetime we look at a Cauchy surface,
${\Sigma}_{0}$
, that is, an everywhere spacelike hypersurface such that any inextendible timelike piecewise smooth curve pierces it once and only once, and a time flow, that is, a smooth timelike vector field,
${t}^{a}$
. From the definition of the Cauchy surface,
${t}^{a}$
is never tangent to it and every point of
$M$
falls in an integral curve of
${t}^{a}$
.
Thus, in assuming the existence of
${\Sigma}_{0}$
we are restricting attention to manifolds of the type
$S\times \Re $
.
To make this structure more apparent we define a function
$t$
by setting to zero the parameter defining the integral curves of
${t}^{a}$
at
${\Sigma}_{0}$
, that is, the value of
$t$
at
$p\in M$
is defined as the value of the parameter the integral curve of
${t}^{a}$
takes at
$p$
if defined in such a way that at
${\Sigma}_{0}$
is takes the value
$0$
. We shall call the surfaces of constant
$t$
by
${\Sigma}_{t}$
. Thus
${t}^{a}{\nabla}_{a}t=1$
, notice that nevertheless, in general, they are not spacelike. When they are spacelike, we shall call
$t$
a time function. In this case we also say we have a spacelike foliation of spacetime
$(M,{g}^{ab})$
. We shall assume that this is the case, but one must take into account that when we are solving for a spacetime we do not know for how long this would continue to hold. Using this structure we can split tensor into “space” and “time” parts with respect to the surfaces
${\Sigma}_{t}$
. If
${n}^{a}$
is the normal to them, that is,
${n}^{a}:={g}^{ab}\frac{{\nabla}_{b}t}{\sqrt{{g}^{cd}{\nabla}_{c}t{\nabla}_{d}t}}$
then
${h}^{ab}:={g}^{ab}{n}^{a}{n}^{b}$
is the induced metric on each
${\Sigma}_{t}$
. We also define the lapse shift pair,
$(N,{N}^{a})$
, as the “timelike” and “spacelike” parts of
${t}^{a}$
with respect to
${\Sigma}_{t}$
, that is,
${t}^{a}:=N{n}^{a}+{N}^{a}$
.
Given the foliation, and the triad
$({h}^{ab},N,{N}^{a})$
we can reconstruct the metric as
${g}^{ab}={h}^{ab}{N}^{2}({t}^{a}{N}^{a})({t}^{b}{N}^{b})$
. Given another foliation,
$({\Sigma}_{0},{{t}^{\prime}}^{a})$
, but the same triad, we get another metric tensor, the relation between both is a diffeomorphism which leaves invariant
${\Sigma}_{0}$
and is generated by the integral curves of
${t}^{a}{{t}^{\prime}}^{a}$
. Alternatively we can choose another pair
$(N,{N}^{a})$
and so construct another metric tensor, the relation between both metric thus obtained is also a diffeomorphism.
Since one is interested in geometries, that is in metric tensors up to diffeomorphisms, both metric obtained are equivalent, but one needs to pick up an specific one in order to write down, (and solve), Einstein's equations.
Using the above splitting vacuum Einstein's equations become two
evolution equations:
$$\begin{array}{ccc}{\dot{h}}_{ab}=& & 2N{h}^{1/2}({P}_{ab}\frac{1}{2}{h}_{ab}P)+2{D}_{(a}{N}_{b)},\end{array}$$  
$$\begin{array}{ccc}{\dot{P}}^{ab}=& & N{h}^{1/2}\left({R}^{ab}\right(h)\frac{1}{2}R(h\left){h}^{ab}\right){h}^{1/2}({D}^{a}{D}^{b}N{h}^{ab}{D}^{c}{D}_{c}N)\end{array}$$  
$$\begin{array}{ccc}& +& \frac{1}{2}N{h}^{1/2}{h}^{ab}({P}_{cd}{P}^{cd}\frac{1}{2}{P}^{2})2N{h}^{1/2}({P}^{ac}{P}_{c}^{b}\frac{1}{2}P{P}^{ab})\end{array}$$  
$$\begin{array}{ccc}& +& {D}_{c}\left({N}^{c}{P}^{ab}\right)2{P}^{c(a}{D}_{c}{N}^{b)},\end{array}$$  
where the dot means a Lie derivative with respect to
${t}^{a}$
,
${P}^{ab}:=\sqrt{h}({K}^{ab}K{h}^{ab})$
and
${K}^{ab}$
is the extrinsic curvature of
${\Sigma}_{t}$
with respect to the fourgeometry, that is,
${K}^{ab}:={h}^{ac}{\nabla}_{c}{n}^{b}$
.
${D}_{a}$
is the covariant derivative associated with
${h}^{ab}$
on each
${\Sigma}_{t}$
, and
${R}^{ab}$
its Ricci tensor.
And two
constraint equations:
$$\begin{array}{ccc}& & hR\left(h\right){P}^{ab}{P}_{ab}+\frac{1}{2}{P}^{2}=0,\end{array}$$  
$$\begin{array}{ccc}& & {D}_{a}{P}^{ab}=0.\end{array}$$  
Note that there are two “dynamical” variables,
${h}_{ab}$
, and
${P}^{ab}$
, while the lapse–shift pair
$(N,{N}^{a})$
, although necessary to determine the evolution, is undetermined by the equations. Note also that they do not enter into the constraint equations, for as said above a change on the lapse–shift pair leave the fields at initial surface unchanged.
It is important step back now and see that these equations can be thought as “living” in a structure completely detached from space–time. To see this, identify all points laying in the same integral curve of
${t}^{a}$
, thus the equivalence class is a three dimensional manifold
$S$
, homeomorphic to any
${\Sigma}_{t}$
. On it, for each
$t$
we can induce spacecontravariant tensors, such as
${h}_{ab}\left(t\right)$
,
${P}_{ab}\left(t\right)$
,
${N}_{a}\left(t\right)$
, and scalars, as
$N\left(t\right)$
. As long as the surfaces are space–like, the induced metric is negative definite and we can invert it, thus we can perform all kinds of contractions and write the above equations as dynamical equations on the parameter
$t$
on fields on the same manifold,
$S$
. This is of course the setting in which one sets most of the schemes to solve the equations, and it is hard to keep control, even awareness, that the surfaces defining the foliation can become null or nearly so. Einstein's equations “feel” that effect since they are causal, and this is abruptly fed back via the development of singularities on the solutions . They do not have any thing to do with real singularities of space–time, but rather with foliations becoming null.
The first application of this idea to Einstein's equations appears to have been [
18]
. There a hyperbolic system consisting of wave equations for the time derivative of the variable
$h{h}^{ab}$
is obtained when the shift is taken to vanish and the lapse is chosen so as to impose the time component of the harmonic gauge. The shift is set to zero, but, as stated in the paper, this is an unnecessary condition. Thus, it is clear that one gains in flexibility compared with the standard method above. This seems to correspond with the fact that in one of the equations forming the system, the time derivative of the evolution equation for the momentum, the momentum constraint has been suitably added, thus modifying the evolution flow outside the constraint submanifold.
As stated in the paper, they could not use the other constraint, the Hamiltonian one, to modify that equation.
The condition for that system to be (symmetric) hyperbolic is that the term below should not have second derivatives of
${P}^{ab}$
.
$$({h}^{ab}\Delta {D}^{a}{D}^{b})(\frac{\dot{N}}{N}\frac{N}{2}P),$$
where
${h}^{ab}$
is the induced three metric on a hypersurface,
${D}_{a}$
is the covariant derivative at that hypersurface compatible with
${h}_{ab}$
,
$\Delta :={h}^{ab}{D}_{a}{D}_{b}$
,
$N$
is the lapse function, and
$P:={h}_{ab}{P}^{ab}$
, with
${P}^{ab}$
the momentum field conjugated to
${h}_{ab}$
.
The simplest condition to guarantee this is:
$$\dot{N}=\frac{{N}^{2}}{2}P,$$
which in view of the definition of
${P}^{ab}$
, which implies
$\dot{h}=hNP$
, has as a solution,
$$N=(\frac{h}{e}{)}^{1/2},$$
where
$\frac{h}{e}$
is the determinant of the metric
${h}^{ab}$
with respect to a constant in time background metric
${e}^{ab}$
. This is precisely the harmonic condition for the time component of
${\Psi}^{b}$
in notation introduced in the paper's introduction. That is,
${\Psi}^{b}{n}_{b}\equiv 0$
, where
${n}_{b}$
it the normal to the foliation. If the determinant of
${e}_{ab}$
is not taken to be constant in time, then one gets,
$$\frac{\dot{N}}{N}\frac{N}{2}P=\frac{\dot{e}}{e},$$
and so the system remains hyperbolic. Thus, we see that, up to the determinant of the metric, the lapse can be prescribed freely. This freedom is very important because it gives ductility to the approach, since this function can be specified according to the needs of applications. We shall call this a generalized harmonic time gauge.
Although in the introduction of [
18]
there is a remark dedicated to numerical relativists about the possible importance of having a stable system, the paper did not spark interest until recently, when applications required these results to proceed. In recent years, a number of papers have appeared which further elaborate on this system, [
2,
3,
4]
. In particular I would like to mention [
3]
, where the authors look at the system in detail, writing it as a first order system, and introduce all variables which are needed for that. In these recent papers, the generalized harmonic time gauge has been included, as well as arbitrarily prescribed shift vectors. If one attempts attempts to write down the system as a first order one, that is, to give new names to the derivatives of the basic fields until bringing the system to the form of equation 1 , the resulting system is rather big, it has fifty four variables, without counting the lapseshift pair. We shall see that there are first order hyperbolic systems with half that number of variables.
Two similar results are of interest: In [
9]
a system is introduced with basically the same properties, but of lower order, that is, only first derivatives of the basic variables are taken as new independent variables in making the system first order. In this paper, it is realized that the same trick of modifying the evolution equations using the constraints can be done by modifying, instead of the second time derivative of the momentum, the extra equation which appears when making the ADM equations a first order system, that is the equation which fixes the time evolution of the space derivatives of the metric, or alternatively the time evolution of the Christoffel symbols. When this equation is suitably modified by adding a term proportional to the momentum constraint, and when the harmonic gauge in the generalized sense used above is imposed, a symmetric hyperbolic system results. In [
10]
the generalized harmonic time gauge is included, as well as arbitrarily prescribed shift vectors. For the latest on this approach see, [
11]
. I shall comment more on this in next section, § 4 .
In [
31]
a similar system is presented. In this case, the focus is on establishing some rigorous results in the Newtonian limit. So a conformal rescaling of the metric is employed using the lapse function as conformal factor. The immediate consequence of this transformation is to eliminate from the evolution equation for
${P}^{ab}$
the term with second space derivatives of the lapse function, precisely the term giving rise to one of the terms in equation 4.1 . The end consequence is that the conformal metric is flatter to higher order. With this re–scaling, and using the same type of modification of the evolution equation for the space derivatives of the metric that the above two approaches use, a symmetric hyperbolic system is found, for arbitrary shift and lapse. ^{
$\text{8}$
}
This freedom of the lapse and shift was used to cancel several divergent terms of the energy integrals in the Newtonian limit by imposing an elliptic gauge condition on the shift, which also determined uniquely the lapse. This resulted in a mixed symmetric hyperbolic–elliptic system of equations. In [
32]
an attempt is made to explore what other possibilities there are of making symmetric hyperbolic systems for general relativity with arbitrarily prescribed lapse–shift pairs. A set of parametrized changes of field variables and of linear combinations of equations are made, and it is shown that there exists at least a one parameter family of symmetric hyperbolic systems. In these systems generalized harmonic time gauge is replaced by:
$$N:=(\frac{h}{e}{)}^{\delta}\delta >0.$$
So, the dependence of the lapse on the determinant of the metric can be modified, but never suppressed. Since the changes in the parameter imply changes in the dynamical variables, while the factor proportional to the momentum added to the evolution equation for the connection is unique, and so fixed, it is not clear whether this can be of help for improving numerical algorithms.
We shall see this type of dependence arise in one of the developments of one of the above mentioned approaches.
In [
29]
a similar system, in the sense of using variables from the
$3+1$
decomposition, is obtained by imposing also the same generalized harmonic time gauge. This system, as is the original system of [
18]
, is of higher order because it includes the electric and magnetic parts of the Weyl tensor in the
$3+1$
decomposition. As such, it contains more variables (fifty) than the two discussed above (thirty).
4.2 The Frame Representation.
Another way of dealing with Einstein equations is through the frame representation. There, instead of using the metric tensor as the basic building block of the theory, a set of frame fields is used. At first this representation seems to be even less economical than the metric tensor representation of geometries, since, on top of the diffeomorphism freedom, one has the freedom to choose the frame vector fields. In short, one has sixteen variables instead of ten. But the antisymmetry of the connection coefficients compared with the symmetry of the Christoffel symbols levels considerably the difference, ending, after adding the second fundamental form and others, with twentyeight variables, instead of the thirty of the conventional system. But this count is not entirely correct. As mentioned above, in order to close the system of equations one has to add the evolution equations for the electric and magnetic part of the Weyl tensor, thus ending with a total of thirty four variables. Actually one can close the system with the twenty four variables at the expense of making the equations into second order wave equations, (See for instance [
55]
.), so effectively adding more variables when re–expressing it as a first order system.
The more important application of the frame representation has been the conformal system obtained by Friedrich, [
27,
26]
, (see § 1 ), where in a fixed gauge he got a symmetric hyperbolic system which allowed him to study global solutions. Later, using similar techniques and spinors, he found a symmetric hyperbolic system with the remarkable property that lapse and shift appear in an undifferentiated form, allowing for greater freedom in relating them to the geometry without hampering hyperbolicity [
28]
. In [
29]
he introduces new symmetric systems for frame components where one can arbitrarily prescribe the gauge functions, which in this case does not only include the equivalent to the lapse–shift pair, but also a three by three matrix fixing the rotation of the frame.
In this case, these gauge functions enter up to first derivatives. This compares very favorably with the ADM representation schema where the lapseshift entered up to second order derivatives.
Friedrich also finds a symmetric hyperbolic system with the generalized harmonic time condition.
Contrary to the systems in the ADM formalism, where the issue is rather trivial, these systems do not seem to allow for a writing in flux conservative form. We do not consider that a serious drawback. The structure of Einstein's equations is very different than those of fluids, where the unavoidable presence of shocks makes it important to write them that way. Indeed the reason fluids have shocks can be attributed to their genuinely nonlinear character, [
43]
, a property not shared by Einstein's theory. (More about this in the next section, § 4 .)
4.3 Ashtekar's Representation.
The example of the one dimensional wave equation 2.2 can be slightly improved by the following construction: We define a two dimensional vector
$u:=({u}^{1},{u}^{2}):=(\phi ,{\phi}_{t}\alpha {\phi}_{x})$
. We then have
${u}_{t}^{1}=\alpha {u}_{x}^{1}+{u}^{2}$
, while,
${u}_{t}^{2}=\alpha {u}_{x}^{2}+(1{\alpha}^{2}){\phi}_{xx}$
. Thus we have a diagonal –and so symmetric hyperbolic– system if
$\alpha =\pm 1$
, namely,
$${u}_{t}=\left[\begin{array}{cc}\alpha & 0\\ 0& \alpha \end{array}\right]{u}_{x}+\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]u,$$
The two possible values
$\alpha $
can take correspond to the two characteristic directions the wave equation defines. This trick can not be extended two more dimensions, basically because the space derivative of
$\phi $
is a vector, and so can not be properly mixed with its time derivative. But in other dimensions one can implement similar schema if the fields are not scalars but appropriate tensors.
Einstein's equations in Ashtekar's' variables [
7,
8]
is one beautiful example of this, since they have the remarkable property of naturally constituting a first order evolution system. Because of this reason it is also a compact system with twentyseven unknowns, even before imposing the reality condition on the connection variable. Recently it has been proven [
40]
that such a system is symmetric hyperbolic if suitable combinations of the constraint equations are added to their evolution equations, thus effectively changing the flow outside the constraint submanifold of phase–space.
In Ashtekar's representation the basic variables are a densitized
$SU\left(2\right)$
soldering form,
${{\stackrel{~}{\sigma}}_{A}^{a}}^{B}$
and a
$SU\left(2\right)$
connection
${{{A}_{a}}_{A}}^{B}$
which are tangent to a spacelike foliation of space time determined by given “lapse”–shift pair
$\stackrel{~}{N}=N/det\sigma $
,
${N}^{a}$
.
The symmetric hyperbolic evolution equation system is:
$${\mathcal{\mathcal{L}}}_{t}{\stackrel{~}{\sigma}}^{b}=\frac{i}{\sqrt{2}}{\mathcal{D}}_{a}\left(\stackrel{~}{N}\right[{\stackrel{~}{\sigma}}^{a},{\stackrel{~}{\sigma}}^{b}\left]\right)+\frac{i}{\sqrt{2}}\stackrel{~}{N}[\stackrel{~}{C},{\stackrel{~}{\sigma}}^{b}]+2{\mathcal{D}}_{a}\left({N}^{[a}{\stackrel{~}{\sigma}}^{b]}\right)+{N}^{b}\stackrel{~}{C}+[{A}_{a}{N}^{a},{\stackrel{~}{\sigma}}^{b}]$$
$${\mathcal{\mathcal{L}}}_{t}{A}_{b}={\mathcal{D}}_{b}\left({A}_{a}{N}^{a}\right)+{N}^{a}{F}_{ab}+\frac{i}{\sqrt{2}}\stackrel{~}{N}[{\stackrel{~}{\sigma}}^{a},{F}_{ba}]+\frac{i}{{\sigma}^{2}\sqrt{2}}\stackrel{~}{N}{\stackrel{~}{\sigma}}_{b}C+\frac{i}{{\sigma}^{4}}\stackrel{~}{N}{{\epsilon}_{b}}^{dc}{\stackrel{~}{\sigma}}_{c}{C}_{d},$$
where
$(\mathcal{D})$
is the
$SU\left(2\right)$
derivative whose difference with respect to a flat connection is
${{{A}_{a}}_{A}}^{B}$
.
$C$
,
${C}^{a}$
, and
$\stackrel{~}{C}$
are the constraint equations,
$$\begin{array}{ccc}C(\stackrel{~}{\sigma},A)& :=& tr\left({\stackrel{~}{\sigma}}^{a}{\stackrel{~}{\sigma}}^{b}{F}_{ab}\right)\end{array}$$ 
(5)

$$\begin{array}{ccc}{C}_{a}(\stackrel{~}{\sigma},A)& :=& tr\left({\stackrel{~}{\sigma}}^{b}{F}_{ab}\right)\end{array}$$ 
(6)

$$\begin{array}{ccc}\stackrel{~}{C}(\stackrel{~}{\sigma},A)& :=& {\mathcal{D}}_{a}{\stackrel{~}{\sigma}}^{a},\end{array}$$ 
(7)

Note that here there is an extra vectorial constraint, a
$SU\left(2\right)$
valued scalar, which corresponds to the fact that the system has extra degrees of freedom, the
$SU\left(2\right)$
rotations. The extra constraint is just a strange way of asserting the symmetry of the second fundamental form, and is of the type of substitutions we made above to improve on the wave equation system. The constraints on themselves satisfy a symmetric hyperbolic system of equations.
Note that the principal part of the system is block diagonal and the eigenvectors–eigenvalues are very simple combinations of
$\stackrel{~}{\sigma}$
with the elements of an orthogonal basis
$\{{\omega}^{a},{m}^{a},{\overline{m}}^{a}\}$
, where
${\omega}^{a}$
is the wave vector.
In this new system, the “lapse”shift pair can be chosen arbitrarily. But in fact the “lapse” that appears here is a scalar density which has already incorporated the square root of
$h$
on it. So the freedom is actually the same as in the ADM representation. As in the frame representation, the lapse–shift pair appears only with derivatives up to first order. In this case it is relatively easy to see the freedom in making up evolution equations for the lapse–shift pair. As said above, the system is symmetric for Ashtekar's variables, since the lapse–shift pair enters as terms with up to first derivatives, one can take those terms from the non–principal part of the system and promote them into the principal part of a bigger system which incorporates the lapse–shift as extra variables. Thus, these terms constitute an off–diagonal block of the bigger principal part matrix. Imposing symmetry to the bigger matrix fixes the opposite off–diagonal part of the matrix.
The only freedom left is on the lapse–shift block–diagonal part, which can be chosen to be any symmetric matrix we like. The non–principal part of the equation system on the lapse–shift sector can also be chosen arbitrarily. Of course, in contrast with the ADM representation results, one can also choose a gauge condition via elliptic equations on the lapse–shift pair. In this case, the elliptic system can be of first order in the lapse–shift or related variables. For instance, one could use Witten's equation to evolve them.
As we have seen, the generalized harmonic time gauge seems to appear naturally in most attempts to get well posed evolution systems. Thus it seems to be really a key ingredient, perhaps with some physical content. One could argue in that direction from the circumstances in which it appears in [
31]
, namely effectively improving the estimates of the behavior of solutions admitting a Newtonian limit that is in a way related to the longitudinal modes of the theory. This longitudinal modes are part of the evolution, although they are not expected to behave in a hyperbolic manner.
See also the comments around equation (9) in [
11]
.
5 Beyond the Prescribed Gauge.
In this section we will look at several attempts to, once a hyperbolic system has been obtained for the evolution variables, extend it into a bigger well posed system where the lapse–shift pair, or more generally the gauge variables of a particular system, can be determined by some prescription which in general depends on the dynamical variables. The bigger system does not need to be hyperbolic, it only needs to be well posed, and so for instance it can be mixed elliptic–hyperbolic.
5.1 Trial and Error Method.
One alternative would be to start with some arbitrary prescription for the gauge variables, evolve the solution for a while, stop, look for troublesome regions, and modify the gauge prescription there. One would do that only a finite number of times and choose smooth prescriptions and smooth transitions between them, so no problem of well posedness or numerical stability would be created by that procedure. With some experience, and luck, this procedure could work.
Perhaps the mayor drawback of this approach is the fact that the strict harmonic gauge, the most conspicuous choice, does not behave well under evolution. This has been known for a long time [
14]
, and more recently new indications have been found in [
5]
. These findings are of course not valid for the generalized harmonic time gauge because one can trivially take any solution, draw a well behaved foliation on it, and identify the generalized gauge that works for that solution. Some attempts have been made in order to get loose from the generalized harmonic time condition, presumably with the intention of later imposing equations on the free variables, otherwise independent fields. In [
1]
a nonstrictly hyperbolic system is found by taking yet another time derivative of the evolution equation for the momentum variable. In that way, they are able to prescribe in a completely free way the lapse–shift pair. In doing this, they obtain a nonstrictly hyperbolic system, in the sense of LerayOhya, [
51]
, which I presume in the language of first order systems means that it is a weakly hyperbolic, but with certain other properties which imply that the system is well posed in Gevrey classes of functions ^{
$\text{9}$
}
.
The resulting system, once brought to first order, has a rather big number of variables. It is not clear that one can stablish numerical stability and convergence for these types of systems, for at least in the continuum estimates an infinite number of derivatives are involved.
In [
29]
, as said in the previous section, a system in the frame representation is made where the corresponding gauge variables can also be given arbitrarily. The importance of this freedom is that in this case one can prescribe directly the lapse. This is in contrast to the case of the generalized harmonic time condition, in which one prescribes the lapse up to the square root of the metric and so finds out what the lapse really was only after solving the problem.
5.2 Hyperbolic Extensions.
One would also like to have recipes which could be used automatically during evolution, that is, algebraic or differential equations, which would not only fix uniquely the evolution of the gauge variables, but which would also result in a well posed evolutionary problem.
One approach has been taken in [
10,
12]
, and [
11]
where the equation 4.1 has been modified to (in the notation of [
18]
):
$$\dot{N}=\frac{{N}^{2}}{2}f\left(N\right)P,$$
With this new evolution equation for the lapse function, they analyze the principal part of the equations and see that if
$f>0$
then it has imaginary eigenvalues and a complete set of eigenvectors.
Thus, up to the smoothness requirement on the eigenvalues with respect to the wave vector, the system seems to be at least strongly hyperbolic and so well posed. This prescription enlarges the system a bit; as one also has to correctly include in it the corresponding equations for the first space derivatives of the lapse function, which become now a dynamical variable.
Although the system is well posed in the sense of the theory of partial differential equations, it has some instabilities from the point of view of the ordinary differential equations. A quick look at the toy model in [
5]
shows that if we take constant initial data for (in that paper's notation)
$\alpha ={\alpha}_{0}$
,
$g={g}_{0}$
, and
${K}_{0}$
, and null data for
$A$
, and
$D$
, then the resulting system is just a coupled set of ordinary equations. One can see that
$g={g}_{0}(K/{K}_{0}{)}^{2}$
, and so
$$\dot{\frac{\alpha}{K}}=\left(\frac{\alpha}{K}{)}^{2}\frac{{K}_{0}^{2}}{{g}_{0}}\right(1f).$$
If
$f=1$
, the harmonic time gauge in this notation, nothing happens at first sight. See nevertheless [
5]
. If
$f>1$
and
$\frac{{\alpha}_{0}}{{K}_{0}}>0$
then we have a singularity in a finite time. The same happens if
$f<1$
and initial data is taken so that
$\frac{{\alpha}_{0}}{{K}_{0}}$
is negative. Thus we see that this gauge prescription can generate singularities which do not have much to do with the propagation modes, and so with the physics of the problem. In [
5]
and [
6]
numerical simulations have been carried out to study this problem.
Needless to say, these instabilities would initially manifest themselves in numerical calculations via the forming of large gradients on the various fields coupled to the above fields, and the time at which they appear depends on the size of the trace of the momentum variable. In [
6]
a proposal to deal with this problem is made which consists of smoothing out the lapse via a parabolic term.
In view of the fact that this problem already arises for constant data, it is doubtful that such a prescription can cure it.
Note that the above prescription for the evolution of the lapse for
$f={f}_{0}>0$
is identical to the one considered in [
32]
, namely equation 4.1 . It is most probably the case then that the same sort of instability would be present there, although the equations considered there are different, due to the inclusion of terms proportional to the scalar constraint in order to render the system symmetric hyperbolic.
In [
5]
there is also a study of another type of singularity which is not ruled out with the choice of the harmonic gauge,
$f=1$
. This singularity seems to be of a different nature, and is probably related to the instability of the harmonic gauge already mentioned. It clearly has to do with the non–linearities of the theory.
It should be mentioned that there are a wide variety of possibilities for making bigger hyperbolic systems out of those which are hyperbolic for a prescribed lapseshift pair, or for the generalized harmonic gauge variant. In that respect, perhaps the systems which are more amenable to a
methodological and direct study are the ones in the frame or in Ashtekar's representations, for there, as discussed in the previous section for the Ashtekar's representation systems, § 4.3 , the possibilities to enlarge the system and keep it symmetric hyperbolic are quite clear and limited.
5.3 Elliptic Extensions.
These are other types of approaches which take more into account the longitudinal modes of the theory, namely those which are related to the energy or matter content of the space time, and those which do not propagate as waves. The nature of these modes implies that these approaches seem to need a global knowledge of the solution, which in practice appears by imposing either elliptic or parabolic equations, the latter as a way to drive the solution close to satisfying an elliptic equation for larger times. Systems of that sort have already been used in applications: In [
20]
, a hyperbolic system with lapse given by an elliptic equation is used in the proof of global existence of small data. The elliptic equation is used to impose the maximal slice condition during evolution, that is
$P\equiv 0$
. In that work, the first order system is for the electric and magnetic parts of the Weyl tensor, while the metric, connection, and extrinsic curvature tensor are obtained by solving elliptic equations on each slice. For their aims, obtaining a priori estimates, this suffices. For numerical simulations of evolution, it is better to solve, as much as possible, evolution equations, and not elliptic ones. Thus for this aim, equations –hopefully hyperbolic or at least parabolic– should be added to evolve the above mentioned (lower order in derivatives) variables. This has improved recently in [
19]
, and [
1]
with a slight generalization to [
20]
in admitting arbitrarily prescribed
$P$
's. In particular, in [
19]
a complete proof of well posedness of mixed symmetric hyperbolic–elliptic systems is given. Such a proof must be implicit somewhere in [
20]
, and a general argument has been given in [
31]
. Surprisingly, such a result, which has a rather simple argument based on the standard elliptic and hyperbolic estimates, has not before had the clean proof it deserves. This gauge has been used to show existence of near Newtonian solutions by [
52]
.
In [
31]
a different elliptic condition is imposed in order to study near Newtonian solutions. An elliptic system is considered for both lapse and shift. It is similar, but not equal, to the above gauge, for in this work a much stronger condition is required on the order of approach of relativistic solutions to the Newtonian limit. This implies globally controlling not only the lapse, but also the shift.
The last two works mentioned hint at some interplay between the problems of finding well behaving gauges for near Newtonian solutions and for long term evolution. The argument has been that, since in this gauge the principal part of the equations is well behaved near the –singular– Newtonian limit, and since the rest of the terms of the hyperbolic system go to zero on that limit, one expects for the time the solution exists to go to infinity as one approaches the Newtonian limit. Thus the gauge should be well behaved until then. I cite [
42]
for recent work on this and [
41]
for a well behaved system in asymptotically null slices amenable to study slow solutions near null infinity.
In the frame and in Ashtekar's representations one could even consider first order elliptic (spinorial) equations to fix the gauge variables. In the frame representation one can even fix gauge variables via an algebraic condition.
6 The Role of the Constraints
In this section, I analyze the role the constraint equations play in these new formulations. I will start by analyzing how the equations were handled in the original proof of well posedness. Then I discuss what has been said recently and what has to be done in the future, both for theoretical considerations and for numerical work. In particular, I discuss the role of the constraints in the initial–boundary value problem, a problem which has still to be fully solved, but which nevertheless is applied in most of the numerical simulations.
6.1 The Constraints in the Harmonic Gauge.
In the harmonic gauge, the role of the constraint equations as such is hidden. This is because the constraint equations become evolution equations due to the gauge choice. Indeed, in the full harmonic gauge all components of the Einstein tensor become a big second order hyperbolic equation for the metric components. The need for the constraints to be satisfied at the initial surface enters the picture because it implies that the first time derivative of the gauge condition must vanish there in order to guarantee that the constraints vanish everywhere. Thus the question about the preservation of the constraint equations does not appear here. At most one can say it has been traded for the question of the gauge consistency.
6.2 The Constraints in the New Systems: Theoretical Considerations.
In the new hyperbolic systems, where covariance is lost, one solves only for the evolution equations. Thus, the question of whether the constraint equations hold during evolution if they hold at the initial surface arises again. If the problem is about the evolution of the whole space–time, or about evolution on the domain of dependence of some space–like surface, then there is a good argument showing that the constraint equations would be satisfied as a consequence of the uniqueness of the system under consideration:
Assume initial data is given at some space–like hypersurface which satisfies the constraints there. We use the new evolution system and get a solution in the domain of dependence of the system, (which, if gauges propagate at speeds greater than light, might be smaller than the domain given by the metric). But using the harmonic gauge I know that there is a solution to the Einstein equation on a maximally extended domain of dependence. If one can diffeomorphically transform metric corresponding to that solution into one satisfying the gauge used for the evolution with the new system, then, since it satisfies all the equations, including the constraints, it follows that it will also satisfy the equations of the new system. Uniqueness of solutions of the new system implies it must be the one found initially and so it also satisfies the constraints. Thus we see that no particular consideration for the constraint equations is needed.
6.3 The Constraints in the New Systems: Numerical Considerations.
For numerical simulations, the role of the constraint equations is delicate: Since there are always numerical errors, although the vector field defined by Einstein's equations in any of the above approaches is tangent to the constraint submanifold, we can only expect to be in the neighborhood of that submanifold. So, since one is effectively modifying the evolution equations outside the submanifold, that the vector field is tangent to it is not enough. For if that submanifold were unstable, it could very well be that a spurious numerical solutions could start growing during evolution and takes us completely away from it. This problem has been noticed by several people and has been considered in detail in [
30]
. There it is shown that, while in some evolution systems, the constraints themselves obey hyperbolic evolution equations, in others that is not the case so are presumably unstable.
It is not clear to me that the condition that the constraint system be well posed is the one needed for considering a system free of this problem. First because well posedness as such is not enough to guarantee the possibility of a numerical scheme: The system could be well posed but still depart exponentially from the constraint submanifold, thus making impossible any reliable calculation. So the non–principal part of the system must also be considered, and probably suitably modified in the neighborhood of the constraint submanifold. Second, since one is never solving, or simulating, the constraint evolution equations, that is, they play no role in the scheme, why should one consider them at all? I think the emphasis should be put on guaranteeing a numerical scheme without spurious solutions; because, as argued above, uniqueness should imply that the constraints are satisfied. Thus, what seems to be needed is a connection between well posedness, or rather no exponential departure from the constraint submanifold, and lack of spurious solutions on the numerical schema.
6.4 The Constraints in the Initial–Boundary Value Problem.
A completely different situation arises if one is considering an initial–boundary value problem for Einstein's equations. Although this problem has not been solved in its full generality for Einstein's equations, it is clear that in order for the constraints to be satisfied during evolution, some of the boundary values have to be chosen in a special way. It is here that the type of equations the constraints satisfy is most important. In particular, if they also form a hyperbolic system, then a study of its principal part at the boundary would tell which conditions are needed to guarantee uniqueness of the solutions, in particular the trivial solution, and so which are the boundary conditions we must force upon the evolution system for the dynamical fields. Since most numerical simulations are in fact initial–boundary value problems, the problem of well posedness and the problem of the propagation of the constraints are central.
References

Abrahams, A., Anderson, A., ChoquetBruhat, Y., and York Jr, J.W., “A nonstrictly hyperbolic system for the Einstein equations with arbitrary lapse and shift”, C. R. Acad. Sci. Ser. II, 323, 835–841, (1996). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Abrahams, A.M., Anderson, A., ChoquetBruhat, Y., and York Jr, J.W., “Einstein and Yang–Mills theories in hyperbolic form without gauge fixing”, Phys. Rev. Lett., 75, 3377–3381, (1996). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Abrahams, A.M., Anderson, A., ChoquetBruhat, Y., and York Jr, J.W., “Geometrical hyperbolic systems for general relativity and gauge”, Class. Quantum Grav., 14, A9–A22, (1997). Related online version (cited on 16 April 1997):
.
☻ open access ✓

Abrahams, A.M., Anderson, A., ChoquetBruhat, Y., and York Jr, J.W., “Hyperbolic Formulation of General Relativity”, (March, 1997). URL (cited on 18 January 1998):
.
☻ open access ✓

Alcubierre, M., “Appearance of coordinate shocks in hyperbolic formalisms of general relativity”, Phys. Rev. D, 55, 5981–5991, (1997).

Alcubierre, M., and Massó, J., “Pathologies of hyperbolic gauges in general relativity and other field theories”, Phys. Rev. D, 57, R4511–R4515, (1998). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Ashtekar, A., “New Hamiltonian formulation of general relativity”, Phys. Rev. D, 36(6), 1587–1602, (1987).

Ashtekar, A., New Perspectives in Canonical Gravity, (Bibliopolis, Naples, Italy, 1988).

Bona, C., and Massó, J., “Hyperbolic evolution system for numerical relativity”, Phys. Rev. Lett., 68, 1097–1099, (1992).

Bona, C., Massó, J., Seidel, E., and Stela, J., “New Formalism for Numerical Relativity”, Phys. Rev. Lett., 75, 600–603, (1995). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Bona, C., Massó, J., Seidel, E., and Stela, J., “First order hyperbolic formalism for numerical relativity”, Phys. Rev. D, 56, 3405–3415, (1997). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Bona, C., Massó, J., and Stela, J., “Numerical black holes: a moving grid approach”, Phys. Rev. D, 51, 1639–1645, (1995). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Bruhat, Y., “Theoreme d' existence pour certain systemes d'equations aux deriveés partielles nonlinaires”, Acta Math., 88, 141–225, (1952).

Bruhat, Y., “Un theoreme d'inestabilite pour certain equations hyperboliques nonlinaires”, C. R. Acad. Sci., 276A, 281, (1973).

ChoquetBruhat, Y., “Espaces temps eineiniens généraux, chocs gravitationnels”, Ann. Inst. Henri Poincare, 8, 327–338, (1968).

ChoquetBruhat, Y., and Christodoulou, D., “Elliptic systems in
${H}_{s,\delta}$
spaces on manifolds which are Euclidean at infinity”, Acta Math., 145, 129–150, (1981).

ChoquetBruhat, Y., Christodoulou, D., and Francaviglia, M., “Cauchy data on a manifold”, Ann. Inst. Henri Poincare A, 29, 241–255, (1978).

ChoquetBruhat, Y., and Ruggeri, T., “Hyperbolicity of the 3+1 System of Einstein Equations”, Commun. Math. Phys., 89, 269–275, (1983).

ChoquetBruhat, Y., and York Jr, J.W., “Mixed elliptic and hyperbolic system for the Einstein equations”, in Ferrarese, G., ed., Gravitation, Electromagnetism and Geometric Structures, International Conference in honour of A. Lichnerowicz, Villa Tuscolana, 19–23 October 1995, (Pitagora, Bologna, Italy, 1996). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Christodoulou, D., and Klainerman, S., The Global Nonlinear Stability of the Minkowski Space, vol. 41 of Princeton Mathematical Series, (Princeton University Press, Princeton, U.S.A., 1993).

Christodoulou, D., and Ó Murchadha, N., “The Boost Problem in General Relativity”, Commun. Math. Phys., 80, 271–300, (1981).

Cutler, C., and Wald, R.M., “Existence of radiating Einstein–Maxwell solutions which are
${C}^{\infty}$
on all of
${I}^{}$
and
${I}^{+}$
”, Class. Quantum Grav., 6, 453–466, (1989).

DeTurk, D., “The Cauchy problem for Lorentz metrics with prescribed Ricci curvature”, Comp. Math., 48, 327–349, (1983).

Fischer, A., and Marsden, J., “The Einstein Evolution Equations as a FirstOrder Symmetric Hyperbolic Quasilinear System”, Commun. Math. Phys., 28, 1–38, (1972).

Fischer, A., and Marsden, J., “General relativity, partial differential equations, and dynamical systems”, in Spencer, D.C., ed., Partial Differential Equations, vol. 23 of Proceedings of Symposia in Pure Mathematics, 309–327, (AMS, Providence, U.S.A., 1973).

Friedrich, H., “The asymptotic characteristic initial value problem for Einstein's vacuum field equations as an initial value problem for a firstorder quasilinear symmetric hyperbolic system”, Proc. R. Soc. London, Ser. A, 378, 401–421, (1981).

Friedrich, H., “On the regular and the asymptotic characteristic initial value problem for Einstein's vacuum field equations”, Proc. R. Soc. London, Ser. A, 375, 169–184, (1981).

Friedrich, H., “On the hyperbolicity of Einstein's and other gauge field equations”, Commun. Math. Phys., 100, 525–543, (1985).

Friedrich, H., “Hyperbolic reductions for Einstein's field equations”, Class. Quantum Grav., 13, 1451–1469, (1996).

Frittelli, S., “Note on the propagation of the constraints in standard 3+1 general relativity”, Phys. Rev. D, 55, 5992–5996, (1997).

Frittelli, S., and Reula, O.A., “On the Newtonian limit of general relativity”, Commun. Math. Phys., 166, 221–235, (1994). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Frittelli, S., and Reula, O.A., “Firstorder symmetrichyperbolic Einstein equations with arbitrary fixed gauge”, Phys. Rev. Lett., 76, 4667–4670, (1996). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Geroch, R., “The Local Nonsingularity Theorem”, J. Math. Phys., 24(7), 1851–1858, (1983).

Geroch, R., “Partial Differential Equations of Physics”, (February, 1996). URL (cited on 19 January 1998):
. Scottish Summer School in Theoretical Physics.
☻ open access ✓

Geroch, R., and Xanthopolous, B.C., “Asymptotic Simplicity Is Stable”, J. Math. Phys., 19, 714–719, (1978).

Gustafsson, B., Kreiss, H.O., and Oliger, J., Time Dependent Problems and Difference Methods, Pure and Applied Mathematics, (Wiley, New York, U.S.A., 1995).

Hadamard, J., Lectures on Cauchy's Problem in Linear Partial Differential Equations, (Yale University Press, New Haven, U.S.A., 1923).

Hawking, S.W., and Ellis, G.F.R., The Large Scale Structure of SpaceTime, Cambridge Monographs on Mathematical Physics, (Cambridge University Press, Cambridge, U.K., 1973).

Hughes, T., Kato, T., and Marsden, J., “Wellposed quasilinear second order hyperbolic systems with applications to nonlinear elastodynamics and general relativity”, Arch. Ration. Mech. Anal., 63, 273–294, (1976).

Iriondo, M.S., Leguizamón, E.O., and Reula, O.A., “Einstein's equations in Ashtekar variables constitute a symmetric hyperbolic system”, Phys. Rev. Lett., 79, 4732–4735, (1997). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Iriondo, M.S., Leguizamón, E.O., and Reula, O.A., “The Newtonian Limit on Asymptotically Null Foliations”, (October, 1997). URL (cited on 17 January 1998):
.
☻ open access ✓

Iriondo, M.S., Leguizamón, E.O., and Reula, O.A., “Fast and slow solutions in General Relativity: The initialization procedure”, J. Math. Phys., 39, 1555–1565, (1998). Related online version (cited on 17 January 1998):
.
☻ open access ✓

John, F., “Formation of Singularities in Elastic Waves”, in Ciarlet, P.G., and Roseau, M., eds., Trends and Applications of Pure Mathematics to Mechanics, Invited and Contributed Papers presented at a Symposium at École Polytechnique, Palaiseau, France, November 28 December 2, 1983, vol. 195 of Lecture Notes in Physics, 194–210, (Springer, Berlin, Germany, 1984).

Klainerman, S., “Uniform decay estimates and the Lorentz invariance of the classical wave equation”, Commun. Pure Appl. Math., 38, 321–332, (1985).

Klainerman, S., “The null condition and global existence to nonlinear wave equations”, Lect. Appl. Math., 23, 293–326, (1986).

Klainerman, S., “Remarks on the global Sobolev inequalities in Minkowski Space”, Commun. Pure Appl. Math., 40, 111–117, (1987).

Kreiss, H.O., “Über sachgemässe Cauchyprobleme”, Math. Scand., 7, 71–80, (1959).

Kreiss, H.O., and Lorentz, J., InitialBoundary Value Problems and the Navier–Stokes Equations, vol. 136 of Pure and Applied Mathematics, (Academy Press, Boston, U.S.A., 1989).

Kreiss, H.O., Nagy, G.B., Ortiz, O.E., and Reula, O.A., “Global existence and exponential decay for hyperbolic dissipative relativistic fluid theories”, J. Math. Phys., 38, 5272–5279, (1997). Related online version (cited on 19 January 1998):
.
☻ open access ✓

Leray, J., Hyperbolic Differential Equations, (Institute for Advanced Studies, Princeton, U.S.A., 1953).

Leray, J., and Ohya, Y., “Equations et Systèmes NonLinéaires hyperboliques nonstricts”, Math. Ann., 170, 167–205, (1967).

Rendall, A.D., “The Newtonian Limit for Asymptotically Flat Solutions of the Vlasov–Einstein System”, Commun. Math. Phys., 163, 89–112, (1994). Related online version (cited on 7 January 1998):
.
☻ open access ✓

Sideris, T., “Formation of Singularities in 3d Compressible Fluids”, Commun. Math. Phys., 101, 475–485, (1985).

Taylor, M.E., Pseudodifferential Operators and Nonlinear PDE, vol. 100 of Progress in Mathematics, (Birkhäuser, Boston, U.S.A., 1991).

van Putten, M.H.P.M., and Eardley, D.M., “Nonlinear wave equations for relativity”, Phys. Rev. D, 53, 3056–3063, (1996). Related online version (cited on 17 January 1998):
.
☻ open access ✓

Wald, R.M., General Relativity, (University of Chicago Press, Chicago, U.S.A., 1984).
Note: The reference version of this article is published by Living Reviews in Relativity