{'fraction_non_alphanumeric': 0.03269098898638342, 'fraction_numerical': 0.010943644429517053, 'mean_word_length': 4.579604903568361, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mixed Hyperbolic - Second-Order Parabolic Formulations of General Relativity

Two new formulations of general relativity are introduced. The first one is a parabolization of the Arnowitt, Deser, Misner (ADM) formulation and is derived by addition of combinations of the constraints and their derivatives to the right-hand-side of the ADM evolution equations. The desirable property of this modification is that it turns the surface of constraints into a local attractor because the constraint propagation equations become second-order parabolic independently of the gauge conditions employed. This system may be classified as mixed hyperbolic - second-order parabolic. The second formulation is a parabolization of the Kidder, Scheel, Teukolsky formulation and is a manifestly mixed strongly hyperbolic - second-order parabolic set of equations, bearing thus resemblance to the compressible Navier-Stokes equations. As a first test, a stability analysis of flat space is carried out and it is shown that the first modification exponentially damps and smoothes all constraint violating modes. These systems provide a new basis for constructing schemes for long-term and stable numerical integration of the Einstein field equations.


I. INTRODUCTION
For years people have tried to obtain analytic solutions of the complex field equations of Einstein's general theory of relativity (GR). Apart from few cases where symmetry is invoked, it is almost impossible to analyze the complicated dynamics in the strong gravitational field regime as that is described by GR. Approximation methods have been developed over the course of time, but the most promising tool for tackling problems such as gravitational waves arising from binary black hole (BBH) or binary neutron star (BNS) mergers, gravitational collapse etc., is numerical relativity.
In the past few years remarkable progress has been made towards achieving long term and stable evolution of the Einstein field equations. Recently particular cases of the long standing problem of the evolution of BBH were solved [4,5,6]. Despite this exceptional achievement, the general problem of long term and stable evolution of the GR equations remains open and there is still much work left to be done. There is no theory or prescription to chose what formulation(s) and what gauge conditions are suitable for the numerical solution of a given problem. For example, there does not seem to exist a definitive explanation of why the approaches of [4,5,6] perform so well when contrasted to previous efforts. Furthermore, there are astrophysical and theoretical problems of great interest for which the formulations above have not been applied yet and it is not known whether they will prove successful in such cases. Such problems are the study of the internal structure of black holes and astrophysical phenomena, where except for black holes matter is also involved.
The numerical integration of the Einstein equations is not an easy task, because the computations can become unstable and an exponential blow up of the numerical error may occur, even when the formulation employed admits a well-posed initial value problem. If the numer-ical techniques, gauge and boundary conditions chosen do not suffer by pathologies, perhaps the most important cause of instabilities during a free evolution is the growth of the constraint violating modes which are excited by numerical errors.
Several methods have been proposed to deal with this last type of instability. A path followed by many relativists is referred to as constrained evolution [9,10,11,12,13,14,15], in which a subset of the dynamical variables is evolved in time by use of the evolution equations and the remaining variables are obtained by solving the constraints after each time-step. This method has proven to be very accurate and has been successfully applied in cases where some type of symmetry is present. Another method, which resembles the approach of constrained evolution, is called constraint projection [20,21]. In this method all variables are evolved in time using the evolution equations and they are periodically "projected" via an optimal algorithm to the closest configuration of dynamical fields which satisfies the constraints. This method is quite robust, but it can be computationally demanding because it requires the repeated solution of non-linear elliptic equations.
Another strategy takes advantage of the fact that in the ideal case where the constraint equations are satisfied, one has the freedom to add combinations of the constraints to the right-hand-side (RHS) of the evolution equations. This approach changes the mathematical properties of the equations, but not the physics these equations describe. A number of strongly and symmetric hyperbolic evolution equations have been constructed using this method, see [16] for a review, which resulted in improved stability of numerical simulations of GR.
By virtue of the freedom to add the constraint equations to the RHS of the evolution equations and/or by extension of the number of dynamical fields, it is also possible to add terms that act as "constraint drivers", i.e., terms which turn the constraint surface into a local attractor. This technique is nowadays usually termed as "constraint damping".
The constraint damping strategy naturally splits in two different categories: a) Additional, non-physical variables are appropriately introduced in the evolution system, so that their time evolution drives the physical variables toward the surface of constraints [18,19], and b) no additional (non-physical) variables are introduced in the system, but only combinations of the constraints, which act as constraint drivers, are added to the RHS of the evolution equations of the physical fields, e.g. [7]. In this case the constraint dissipation is inherent to the evolution equations of the physical variables themselves.
The goal of many formulations in numerical relativity has been to incorporate such drivers in symmetric or strongly hyperbolic systems without changing the principal part of the evolution equations [18,19]. The common property of all these formulations is that the damping rate of the short-wavelength constraint-violating modes is independent of the wavelength of these modes. As a result, the constraint violations are not smoothed during a free evolution, unless numerical dissipation is artificially added to a code. The existence of short-wavelength components may lead to spurious oscillations and potentially terminate the computations [4].
For this reason, it may be advantageous to consider formulations of GR which possess wavelength-dependent constraint damping. One way to achieve this goal is to construct formulations of GR that under free evolution force all or some of the constraints to evolve according to parabolic equations. Parabolic equations damp short wavelengths extremely efficiently and for this reason they posses desirable smoothing properties [29].
The purpose of this paper is to present two new parameterized 3+1 formulations which fall in category b) of the constraint damping strategy. These two formulations will be referred throughout this work as system 1 and system 2. The difference between these two systems and most formulations of GR, which belong to the same category, is that constraint dissipation does not only enter through low order terms, but also through higher order terms since the order of the of the evolution equations is increased. In system 1 and system 2 all or some of the constraint violating modes have been parabolized.
The desirable property of the well-posed system 1 is that the evolution equations of the constraints become second-order parabolic, independently of the gauge conditions and spacetime configuration, which in turn implies that the constraint surface becomes a local attractor. The evolution equations of system 1 may be classified as a set of mixed hyperbolic-second-order parabolic partial differential equations (PDEs).
The well-posed system 2, is a parabolization of what is usually referred to as the Kidder, Scheel, Teukolsky (KST) formulation which was introduced in [2] and later extended in [3]. A generalization of the gauge condition chosen in [2] is given and a certain combination of the terms of system 1 is added to the RHS of the KST formu-lation. Unlike system 1, system 2 forces only a subset of the constraint variables to evolve according to parabolic equation and the constraint propagation system is mixed strongly hyperbolic -second-order parabolic (MHSP).
It is here noted that hyperbolic and parabolic drivers have been constructed for the MHD equations, by introducing additional non-physical fields, whose time evolution drive the physical fields on the constraint surface, see for example [22] and references therein. A similar approach has also been applied in numerical relativity via a Hamiltonian constraint relaxation method [23], which introduces a parabolic equation to approximately solve the Hamiltonian constraint. The parabolic equation is solved not in real (coordinate) time, but in some fiducial time until the Hamiltonian constraint has been relaxed. The difference between the suggested approach in this work and the aforementioned two efforts is that here there is no need for extra variables and that the constraint damping is inherent to the evolution equations of the physical dynamical variables, which implies that the damping takes place in real time and not in some fiducial time. In addition to that, in contrast to the approach of [23], system 1 dissipates the perturbations which violate not only the Hamiltonian constraint, but all constraints of the formulation.
Note also that a general, robust constraint driver has been developed and successfully applied to two formulations of Maxwell's equations in [24], where the evolution equations were cast into a type of MHSP set of equations. In [25], the same approach has been successfully applied to a "1+1" case of the ADM evolution equations, linearized about flat space, where the resulting evolution system resembled the structure of a mixed hyperbolic -fourth-order parabolic set of PDEs. The approach of [24] can in general be applied to the full non-linear equations of the first-order ADM formulation, which is re-introduced in this paper (overcoming this way the obstacle of fourth-order evolution systems). However, the main disadvantage of this method is that the resulting evolution system would turn from quasi-linear into highly non-linear, see [25], because the constraints would be added in quadrature. On the other hand, the ideas suggested in this work result in quasi-linear PDEs and this constitutes one of the major differences between this paper and the work in [24].
Finally, it must be mentioned that the idea of adding multiples of the derivatives of the constraints to the RHS of the evolution equations has been independently suggested in [26], where the approach was applied successfully to different formulations of Maxwell's equations.
This paper is organized as follows. In section II, a model equation is presented to illustrate the parabolization approach. The first order ADM formulation and gauges are presented in section III. In section IV system 1 is presented and studied. A review of the KST formulation is presented, a generalization of the gauge of the original KST formulation is made and system 2 is given and studied in section V. In section VI a study of the stability properties of flat space with three different formulations is presented. In section VII a brief discussion of the choice of parameters is given. Finally, the work is concluded in section VIII.

II. A SIMPLE MODEL EQUATION
To illustrate the strategy that will be employed in order to modify the Einstein equations a scalar wave equation for the variable g(x, t) is introduced. A linear onedimensional scalar wave equation is given by where ∂ x denotes a spatial derivative and ∂ t a time derivative. Equation (1) is equivalent to the following first-order system Equations (2) have the same totality of solutions as equation (1) only if the following constraint is satisfied at all times Equation (3) is essentially the definition of the D variable. In addition to the evolution equations of the dynamical variables g, K, D, one can also obtain the evolution equation of the constraint C. By taking a time derivative of the constraint one finds The initial value problem in GR is similar to that of equations (2) and (3), in the sense that there is a set of evolution equations for the dynamical variables and a set of constraints that have to be satisfied at all times. In GR however, the equations are so much more complex that during their numerical integration the constraints may grow exponentially with time leading the computations to a termination. This is not the case with equation (2), because (2) is a set of strongly hyperbolic PDEs (definitions and theorems pertaining to hyperbolic PDEs can be found in [29]) and hence a perfectly well-behaved one. It can be solved numerically, without any difficulty in satisfying the constraint (3) at all times. Therefore, the purpose of the model equation is to demonstrate a strategy in order to suppress very efficiently the growth of constraint violating modes.
This strategy is based on the fact that not only the constraint C has to be zero at all times, but also its spatial derivative. This implies that one can add combinations of those two to the RHS of the evolution equations (2), without changing the solutions of the equations. With that in mind, the following modification is suggested in order to damp any constraint violating modes and still retain the well-posedness of the system.
In order to prove that (5) admits a well-posed Cauchy problem, it is instructive to present from [29] some definitions and theorems concerning parabolic and mixed parabolic-hyperbolic systems of quasi-linear partial differential equations .
Consider a quasi-linear even-order operator of the form where and where i=1,d ν i = ν, ∂ s νs denotes the partial derivative with respect to x s of order ν s , s = 1, . . . , d with d the number of spatial dimensions, u is the column vector of the n dynamical variables and A ν are n × n matrix functions of the dynamical variables, but not their derivatives.
b. Theorem 1. The Cauchy problem for a parabolic system ∂ t u = P (∂ x )u is well-posed.
Assume now that a system of partial differential equations can be written in the following form where u and v are dynamical variables, P 2 = P 2 (∂ x ) is a second-order differential operator, P 1 = P 1 (∂ x ) is a firstorder operator and the coupling terms R 12 = R 12 (∂ x ) and R 21 = R 21 (∂ x ) are general first-order differential operators. Equations (10) may or may not have zerothorder couplings. Further assume that the uncoupled systems are second-order parabolic and first-order strongly hyperbolic, respectively.
c. Definition 2. The coupled system (10) is called mixed strongly hyperbolic -second-order parabolic.
The coupling terms R 12 and R 21 or any low order terms do not destroy the well-posedness of the system and the following theorem can be proved.
d. Theorem 2. The Cauchy problem for mixed second-order parabolic-strongly hyperbolic systems (10) is well-posed.
Equation (5) is a system of mixed strongly hyperbolic and second-order parabolic PDEs, provided that λ > 0. This can be most easily seen, if the explicit form of (5) is written It can be readily realized that equations (12)- (14) have the form of (10) with u = g and v = {K, D}. The coupling terms are R 12 = −λ∂ x D and R 21 = ζ∂ x g which are first order. Without these first-order couplings, equation (12) is second-order parabolic when λ > 0, whereas the system of (13) and (14) is strongly hyperbolic. Therefore, the total set of equations is MHSP and admits a well-posed Cauchy problem. The most interesting aspect of the modification introduced in (5) is the structure of the evolution of the constraint C. A straightforward calculation yields Equation (15) is parabolic, provided that λ > 0. Therefore, not only the evolution equations of the dynamical variables are well-posed, but also the evolution of the constraint. By Fourier transforming (15), one can show that damping of the constraint C takes place, if λ > 0 and ζ > 0. This is because the damping rate of C is of the form exp(−(λκ 2 + ζ)t), where κ is the magnitude of the wavenumber.
By the form of the damping rate of C it is evident that the value of ζ can have a major impact on the long-term stability of the numerical integration of (5), even though theorem 2 dictates that equation (5) admits a well-posed Cauchy problem irrespectively of the value of ζ.
At this point it is instructive to try to understand why this is so. Note that the constraint and its derivative were added to the RHS of equation (2). Note also that the derivative of the constraint C may not be violated, while the constraint itself may be violated, as long as the violation occurs by a constant. Constant violations correspond to infinite wavelength and hence κ = 0. This is the reason why there is no damping, if ζ = 0 and κ = 0. To be more specific, at the RHS of the evolution of variable g the derivative of the constraint has been added. So, a violation of the constraint C by a constant is not "felt" by the λ term. The λ term damps only the finite wavelength constraint violating perturbations, because κ = 0 perturbations violate not only the constraint, but also the derivatives of the constraint. If ζ < 0 then in order for damping to take place, the following inequality has to be satisfied for all allowed κ In a computer there are limits on the wavenumber of the perturbations, because of the finite length of the computational domain and the mesh size. The maximum κ corresponds to the Nyquist frequency and the minimum to κ min = 2π/L, with L the length of the computational domain. Hence, condition (16) is satisfied, if λ > −ζ/κ min 2 . However, because of the efficiency of the damping very soon all allowed finite wavelength constraint violating perturbations go to zero. As a result the λ term in (5) is switched off, since there are no further violations of the derivatives of the constraints. Therefore, the evolution of the constraint C is effectively given by equation (15) with λ = 0. If in that case there is any residual κ = 0 violation of the constraints (which is inevitable since an exact zero initial value of a constraint is never achieved in a computer), this violation will obtain explosive behavior (because ζ < 0 was assumed), which will occur at a time-scale ζ −1 . Therefore, for long-term and stable numerical integration both λ > 0 and ζ > 0 are required.
A final contradiction has to be resolved here. It was stated earlier that system (5) admits a well-posed Cauchy problem irrespectively of the value of ζ, as long as λ > 0. However, in the previous paragraph it was argued that if ζ < 0 an exponential blow up of the numerical solution occurs. The resolution of this contradiction takes place, if one realizes that well-posedness does not in general guarantee global in time existence of solutions, but only short time existence.
The difference between the suggested modification of this section and previous constraint damping strategies is the λ term. The advantage of the new approach is that the shortest the wavelength of the constraint-violating perturbations the more efficiently those perturbations are getting damped. And it is the short wavelength perturbations, which produce the most spurious oscillations that may lead to blow up of the numerical error. Of course ζ is important to damp the long wavelength (including κ = 0) perturbations, which do not get damped very efficiently by the λ modification.
An important result of the wavelength-dependent damping timescale is that the constraint violating noise is smoothed with time. Another important property of the suggested modification is that it assigns non-zero "speed" of propagation to the static mode present in (2). This is may be an advantage because, according to [17], the presence of zero eigenvalues can be responsible for numerical instabilities caused by perturbations by low order terms.
The same parabolization strategy will be applied to the Einstein equations in section IV and it will be demonstrated that with appropriate modifications all constraints of the theory can be damped.

A. 3+1 decomposition of spacetime
The four-dimensional spacetime of the theory of general relativity can be cast into a 3+1 decomposition [1]. This is achieved by assuming that the spacetime can be foliated by a one-parameter family of spacelike hypersurfaces. The four-dimensional metric is then written in the following form where γ ij is the positive definite 3-metric on the t = const. hypersurfaces, α is the lapse function, β i the shift vector, and x i are the spatial coordinates, i = 1, 2, 3.
To specify the spacetime foliation the gauge equations for the lapse function and the shift vector have to be given. These can generally be written as follows [28]. In this work only algebraic gauges of the following form are considered where γ is the determinant of the three-metric γ ij , and the shift vector β i may be either constant or a fixed function of the spacetime coordinates (t, x i ). In [27,28] it was shown that (19) makes the evolution well-posed on the surface of constraints, if Working with (19), one has the flexibility to use either "1+log" slicing or a densitized lapse where Q is constant or a fixed function of t, x i and σ the densitization parameter. If (21) is used, A = 1/α, whereas if (22) is employed, A = σ. Both gauge conditions have been used rather successfully in numerical relativity. However, note that "1+log" slicing has better singularity avoidance properties than the densitized lapse [30].

B. Evolution equations and constraints
The Einstein equations in 3+1 split form are referred to as the ADM formulation and consist of two subsets of equations. The first subset is that of the evolution equations, which describe how the dynamical variables evolve in time. The second subset is that of the constraint equations, which have to be satisfied for all times. The standard second-order ADM formulation [1] has as dynamical variables the 3-metric γ ij and the extrinsic curvature K ij of the 3D spacelike hypersurfaces.
The simplest first-order form of the ADM formulation is derived from the second-order one by introducing additional dynamical variables and then deriving the evolution equations for D kij . The evolution equations of the first-order ADM arê where K = γ mn K mn is the trace of the extrinsic curvature, ∇ i is the covariant derivative operator associated with the three-metric, R ij is the Ricci tensor associated with the three-metric, and∂ o = ∂ t − £ β , with £ β the Lie derivative along the shift vector β i . The Lie derivatives of the dynamical variables are and The subset of constraint equations is where R is the trace of the three-Ricci tensor. H is the Hamiltonian constraint, M i the momentum constraints and C kij are new constraints due to the introduction of D kij . The equations as presented above do not include matter terms, because this work focuses on vacuum solutions of the Einstein equations.
A final set of twelve constraints can be derived by virtue of (23), via taking a spatial derivative and keeping in mind that partial derivatives commute. Then D kij must also satisfy Constraints (33) relate to (32) via the following equation It is important to note here that in equations (24)-(31), all spatial partial derivatives of the three-metric are replaced by the corresponding D kij variable. However, in doing so an ambiguity arises due to the fact that there is not a single way to replace the mixed second-order derivatives of the three-metric by first order derivatives of the D kij variables. For example the derivative ∂ 1 ∂ 2 γ ij can be replaced by either ∂ 1 D 2ij or ∂ 2 D 1ij . There is no rule or recipe on how to map those mixed derivatives of the standard ADM formulation to the first-order derivatives of the D kij variables, when deriving the 1st order ADM formulation. Hence, one has to make a choice to lift this ambiguity. No matter what that choice is, in the continuum limit where C kij = 0, there is no difference between the numerous final forms the first-order ADM formulation may obtain. However, in the discrete limit where the satisfaction of the C kij constraints is at best approximate, there may be some difference.
In this work the aforementioned ambiguity is lifted by choosing the most natural way in replacing the mixed derivatives of the three-metric, i.e., they are written in as most symmetric way as possible. This is done by following the suggestion of [2], according to which This last modification is equivalent to arbitrarily replacing partial derivatives of γ ij by D kij in the spatial Ricci tensor R ij , and then recalculating the Ricci tensor as followsR whereR ij stands for the recalculated Ricci tensor. It is noted here, that it seems that the last term in the parenthesis makes the Ricci tensor asymmetric under the interchange of its indices, because that term satisfies C ijkl = −C jikl . However, a straightforward calculation shows that the addition of this term is the only way that the Ricci tensor can be cast into a manifestly symmetric object, in the context of the 1st order ADM formulation. Note that a more general Ricci tensorR ij can be calculated by adding to the RHS of (36) extra terms as followsR This modification is redundant for system 1 because the system is well-posed and has constraint damping when ψ = 0, see section IV. However, ψ = 0 is essential for the well-posedness of system 2, see section V.
C. Evolution of the constraints of 1st-order ADM By taking the time-derivative,∂ o , of the constraint equations of the first order ADM formulation one can derive the evolution equations for those constraints by replacing the time derivatives of the dynamical variables through the evolution equations and then eliminate the appropriate combinations of spatial derivatives of the dynamical variables via the constraint equations. After some tedious algebra, one finds that the evolution of the constraints, without using (36), is given bŷ are functions of the dynamical variables and their spatial derivatives, as well as of the lapse function and its spatial derivatives. Their exact forms are long and they are not presented here for brevity, but the full equations are available by the author upon request. The Lie derivatives of the constraints along the shift vector are From equations (38)-(41) it is expected that if the constraints are satisfied on the initial time-slice, they should be satisfied throughout the entire evolution. It is noted here, that if equation (36) is used for the Ricci tensor, it is straightforward to show that the Ricci scalar and hence the Hamiltonian constraint remain unaffected. However, since the evolution equations of the extrinsic curvature are modified, the resulting expressions of the evolution of the momentum and Hamiltonian constraints change. Their exact expressions are not presented here, but their derivation is straightforward, if one uses equations (A1)-(A6) of appendix A.
The principal part of (38)-(41), has only two non-zero eigenvalues, which propagate with the fundamental speed ±α, when β i = 0. All remaining eigenvalues are zero (static modes). Static modes are present in the ADM constraint evolution system, even if β i = 0 because the time derivative of C kij is zero independently of the shift vector. This is where the argument in [17] fits in, to explain why several formulations of GR blow up when integrated numerically. Just like the ADM constraint propagation equations, the evolution of the constraints with other formulations, such as KST, have a set of 0 eigenvalues. This implies that when constraint violating perturbations arise during numerical integration of the evolution equations, these perturbations remain in the system because they cannot propagate off the computational domain. As a result, continuous perturbations by numerical errors lead to accumulation of constraint violating noise which eventually terminates the computations.
Note, however, that the reason why the ADM formulation has not been successful in multi-dimensional simulations is because the evolution system is ill-posed in more than one dimensions. This section is concluded by noting that since the ADM constraint equations have to be satisfied at every single point of a spacelike hypersurface, not only the constraints themselves must be zero, but also their spatial partial derivatives of any order have to be zero, i.e., This means that one is free to add to the RHS of the ADM equations not only combinations of the constraints themselves, but also combinations of their spatial derivatives. This property, will be used in subsequent sections to obtain parabolized formulations of general relativity.

IV. PARABOLIZATION OF THE 1ST-ORDER ADM FORMULATION
In this section equations (24) to (26) are modified in such a way as to have as many parabolic modes as possible. Parabolic modes are modes which decay exponentially with time at a rate exp (−ǫκ 2 t), with ǫ > 0 a damping parameter and κ 2 the magnitude squared of the wavevector. With that in mind, the following modification is proposed The added derivatives of the Hamiltonian and momentum constraints are given by rather complicated and long expressions and for this reason they are not presented here. However, their full expressions are available by the author upon request. It is only noted that in those expressions whenever derivatives of the 3-metric, ∂ k γ ij , are present they are replaced by the variables D kij . This ways one ensures that the resulting system of PDEs is quasi-linear instead of non-linear.
A. Structure and Well-Posedness of system 1 Any new formulation of GR intended for numerical simulations is required to admit a well-posed initial value problem. This is because an ill-posed formulation cannot lead to long-term and stable numerical simulations of the Einstein Equations. This sub-section establishes the wellposedness of system 1 and discusses the structure of its evolution equations.
In the analysis that follows ζ is considered to be zero because the results with ζ = 0 are rather complicated. It has been found that, just like in the model equation of section II, ζ has no impact on the well-posedness of system 1, but it must satisfy ζ > 0 for long-term and stable numerical simulations.
For the study of well-posedness of a given set of evolution equations the zeroth-order terms are omitted and only the higher-order terms have to be considered. For the family of gauges (19) the highest order operator of system 1 is of second-order and according to [29] it has to be analyzed first.
The second-order part is given by the following equations where the symbol "≈" means equal to the second-order part. It can be readily seen that the evolution equations of γ ij , K ij and D kij decouple to second order. This implies that one can study each set of equations separately. The evolution equations for the three-metric (50) are quasi-linear second-order parabolic, if λ > 0. All six eigenvalues of P (iκ) are equal to and for λ > 0 condition (9) is satisfied.
(54) Therefore, K ij does not evolve according to parabolic equations because condition (9) cannot be satisfied due to the existence of zero eigenvalues. An analogous analysis of equations (52) yields the following eigenvalues Thus, condition (9) cannot be satisfied and D kij do not evolve according to parabolic equations either.
Therefore, system 1 is not parabolic and thus theorem 1 cannot be used to prove well-posedness. However, the eigenvalues of the Fourier transformed second-order part of the evolution equations constitute a diagnostic tool for well-posedness and in order for system 1 to be well-posed, the following conditions are necessary Before a sufficient proof of well-posedness is given, note that in appendix B it is shown that if the Fourier transform P (iκ) of the principal part of a quasi-linear spatial operator P (∂ x ) a) satisfies the Petrovskii condition [33] For general second order systems which are not parabolic, the first order terms have to be taken into consideration in the definition of the principal part [29]. Hence, one has to consider both second and first order terms in the principal part of system 1.
For the family of gauges (19) with zero shift vector, the Fourier transformed principal part of system 1 has 28 non-zero eigenvalues, 22 of which have a non-zero real part and the remaining 6 are purely imaginary. Due to the Petrovskii condition, only the real part of the eigenvalues is important for establishing well-posedness. It turns out that the Petrovski condition can be satisfied only when If (57) is true, the real parts of the 22 complex eigenvalues are equal to the 22 non-zero eigenvalues (53)-(55) of the second-order part of system 1. Hence, if both (56) and (57) hold true, the Petrovskii condition is satisfied with c = 0 because the real part of all eigenvalues is strictly non-positive. The Fourier transform of the principal part of system 1 is found to have a complete set of eigenvectors, if one of the following conditions holds true Finally, a calculation of the eigenvectors of P (iκ) of system 1 shows that the eigenvectors can be chosen as analytic functions for κ = 0. In particular, their components can be chosen as sums and products of polynomials of κ i and polynomials of κ.
Therefore, if (56) and (57) are satisfied and one of equations (58) or (59) is true, then conditions a), b) and c) above are met and system 1 admits a well-posed initial value problem. This holds true, if a non-zero fixed shift vector is assumed.
Note here, that the presence of the lapse function as a multiplicative factor in the added terms of system 1 is not necessary for well-posedness and omitting it may prove useful. This is because the damping rate of the constraints is of the form exp(−αλκ 2 t), and hence when the lapse function is small, the damping time-scale may become large. Examples where the lapse function attains small values are black hole spacetimes. It is well-known that α becomes small near the singularities. A way to overcome the difficulty of inefficient damping, when the lapse function is small, is to allow large values for λ. However, this could result in numerical complications (see section VII). Therefore, it may eventually be better to remove the lapse function from the damping terms because this results in a damping rate of the form exp(−λκ 2 t).
Another important subject is the structure of the evolution equations of system 1. In addition to the complex eigevalues, the principal part of system 1 has 8 eigenvalues with zero real parts. These are a subset of the eigenvalues of the principal part of the ADM equations and are: 0 with multiplicity two, ±iακ each with multiplicity two and ±iα √ 2Aκ. Note that, if a non-zero shift vector is considered, the two zero eigenvalues become β i κ i . These 8 modes correspond (as is illustrated in section VI) to constraint-satisfying degrees of freedom, which are the four gauge modes and the four physical modes. In [27] it was shown that for gauge (19) the constraint satisfying modes are described by strongly hyperbolic equations.
Thus system 1 has a set of hyperbolic modes and a set of parabolic modes and for this reason it may be classified as a mixed hyperbolic -second-order parabolic (MHSP) set of equations. However, notice that it does not really fulfill the requirements of definition 2.
The goal of maximazing the number of non-zero eigenvalues has been achieved by system 1. All degrees of freedom now evolve according to a principal part. This, according to [17], may be an advantage for numerical applications.
In addition, something even more important has been achieved: The parabolic modes correspond to constraint violating perturbations, which are now getting damped. This can be more easily realized, if one considers a solution to the equations and apply small amplitude constraint satisfying and small amplitude constraint violating perturbations separately. If constraint satisfying perturbations are imposed on the system, then an exact cancellation of the added terms of system 1 occurs and the system reduces to that of ADM. A thorough study of the evolution of constraint satisfying perturbations has been carried out in [27] and hence this is not the subject of this work. If however, one imposes constraint violating perturbations on system 1 then the parabolic terms are turned on. As a result, the constraint violating perturbations and in particular the short-wavelength ones, have no other choice, but to decay exponentially with time. This will also be realized by the structure of the evolution equations of the constraints in the following section.
B. Evolution of the Constraints with System 1 By virtue of system 1 and equations (A1)-(A6), it can be shown that the constraint evolution equations with system 1 arê where the exact evolution equations of C kij and C klij have been presented, but only the principal parts of the evolution equation of H and M i have been given ("≃" means equal to the principal part).
The eigenvalues of the Fourier transform of the principal part of system (60)-(63) are all non-zero and they equal the non-zero eigenvalues of the second-order part of system 1 (53)-(55). Thus, equations (60)-(63) are second-order parabolic, when condition (56) is satisfied and hence they admit a well-posed Cauchy problem.
The parabolic structure of the evolution equations of the constraints implies that the constraints are damped with time. In particular, the shortest the wavelength of perturbations which violate the constraints, the faster those perturbations are damped because the damping rate is of the form exp(−αλκ 2 t). The new evolution equations of the constraints suggest that the surface of constraints is a local attractor.
Another important property of system 1 is that the constraint violatiions are smoothed with time. This is due to the wavelength-dependent nature of the damping timescale.
Furthermore, the structure and hence the smoothing and damping properties of equations (60)-(63) is independent of the gauge condition or the spacetime solution, since no specific gauge or spacetime geometry was assumed in the derivation of the new constraint propagation equations.
Finally, equations (60), (61) and (63) are completely decoupled to highest-order from the rest of the system and are manifestly second-order parabolic, when (56) is satisfied. However, (62) is coupled to second-order to (61).
At this point, one may wonder if the aforementioned discussion on damping and well-posedness is relevant when equation (34) is enforced at the RHS of (62). This is so, because this enforcement turns (62) into thirdorder. Even in this case, it is straightforward to show that the eigenvalues remain the same and that the constraint propagation system remains well-posed. The proof of well-posedness of the constraint evolution system, with (34) enforced at the RHS of (62), is exactly the same as the proof of well-posedness for parabolic systems presented in [29]. For this reason it is not repeated here. Therefore, the properties of the evolution equations of the constraints are independent of whether (34) is enforced or not. This is to be anticipated because all forms the constraint propagation equations (60)-(63) may obtain by use of (34) are equivalent with each other and physical conclusions do not depend on the form of the equations.

V. PARABOLIC MODIFICATIONS TO THE KST FORMULATION
In this section system 2 is introduced. This system is based on the KST formulation as presented in [2]. Thus, before going into the details of system 2, first the KST formulation is reviewed and the original gauge choice of the KST system is extended to the family (19).

A. The KST formulation
The KST modification to the first-order ADM formulation, as originally introduced, is a parameterized family of strongly hyperbolic formulations.
The KST modification is In addition to these modifications to the RHS of the ADM system, Kidder, Scheel and Teukolsky used a densitized lapse (22) and a fixed shift vector β i .
In this work the lapse condition is extended to the algebraic family (19). The resulting KST system is essentially the same as the original one, but with one difference: σ must be replaced by the quantity A defined in equation (20). For illustration, the principal part of the second-order covariant derivative of the lapse is The RHS of (65) is exactly the same as the RHS it has in the original KST system [2], but with σ replaced by A.
The 30 characteristic speeds of the principal part of the KST formulation in the generalized gauge (19) are equal to {0, ±α, ±c 1 α, ±c 2 α, ±c 3 α}. Eighteen of the speeds are 0. Speeds c 1 , c 2 , c 3 are different than those of the original KST formulation and are given by where the squares of the speeds is shown. For the system to be well-posed the RHS of equation (66) must be positive. The analysis of well-posedness of the KST system in the extended gauge is the same as that of the original one. Hence, if the RHS of equation (66) is positive, the new evolution system is strongly hyperbolic unless one of the following conditions is satisfied Ensuring globally that conditions (67) will not be satisfied is not as straightforward anymore due to the spatial dependence of A. However, there is a way around this difficulty. Notice that if speeds c 2 , c 3 are set to unity, the strong hyperbolicity conditions are met. This can be done irrespectively of whether densitized lapse or "1 + log" slicing is chosen, by requiring that the following conditions be satisfied Solving these equations parametrically one obtains The choice χ = −8/5 in (69) renders the system weakly hyperbolic and hence ill-posed. The solutions above coincide for the following values The generalization of the gauge constitutes the first extension of the original KST formulation.

B. Evolution of the constraints of the KST modification to ADM
By virtue of equations (A1)-(A6) , one can show that the evolution of the constraints with the KST formulation becomes∂ where the exact evolution equations of the C kij , C klij constraints are presented, but only the principal part of the evolution of H, M i is given. The apparent discrepancy in these equations with the equations presented in [2], arises because of the difference of the definitions of H and C klij . In particular, H = 2Ĥ and C klij = 2Ĉ klij , whereĤ and C klij denote the definition of those constraints in [2].
The characteristic speeds of the Fourier transform of the principal part of equations (72)-(75) are {0, ±c 2 α, ±c 3 α}, i.e., they are a subset of the characteristic speeds of the evolution equations, as expected. It was shown in [2] that if the evolution equations of the dynamical variables are strongly hyperbolic, then the constraint propagation equations are also strongly hyperbolic.
Although this property is really desirable, the KST formulation has a serious disadvantage. Equation (72), shows that C kij evolve according to lower order terms. This implies that "static" modes correspond to the principal of (72). These low order terms of (72) are functions of the momentum constraints. This in turn means that any violations of the momentum constraints will directly lead to violations of C kij . This is not desirable and therefore damping of C kij should be introduced in the KST system.

C. Mixed Hyperbolic-Parabolic KST Formulation
In this section the generalized KST system of section V A is modified so that the resulting formulation has the desirable properties of strong damping and smoothing of C kij . This is achieved by the same parabolization strategy that was applied to the scalar wave equation of section II. In particular, the following modification is suggested∂ where (KST ) stands for the RHS of the extended KST formulation as that is given in section V A. Equations (76) together with the evolution equations of K ij from equations (64) are called system 2 throughout this work. For λ > 0 the evolution equation of the three-metric is a parabolic equation (see section V C 1). Notice also the term ζαC kij which has been added to the evolution of D kij . This is important to couple the C kij constraint violating modes which arise from the evolution equations of the three-metric, to those which arise from the evolution equations of the D kij variables. It is straightforward to check that if λ = 0, then there still are eighteen zero eigenvalues in the principal part of the formulation, whereas if ζ = 0 then there are twelve zero eigenvalues. As in equation (5), if both λ > 0 and ζ > 0, strong damping of C kij occurs. The λ term takes care of the short wavelength perturbations, whereas the ζ term takes care of the long wavelength ones.
Note that the proposed modification can be easily integrated with the system of redefined variables, which was introduced in [2]. This is so, because in that system only K ij and D kij were transformed. This may be advantageous, because the redefined system has shown better stability properties in numerical simulations than those of the KST formulation presented in section V A.

Well-posedness of system 2
To prove well-posedness for system 2, note that the equations can symbolically be written aŝ It can be readily realized that equations (77)-(79) have obtained the form of a MHSP system (10) with u = γ ij and v = {K ij , D kij }. The coupling terms are R 12 = −αλγ ab ∂ a D bij and R 21 = ζα∂ k γ ij which are first order. The uncoupled (77) is a second-order parabolic equation provided that λ > 0 and the uncoupled system of (78) and (79) is strongly hyperbolic. Therefore, the total set of equations is mixed second-order parabolic -strongly hyperbolic, bearing thus resemblance to the compressible Navier-Stokes equations. According to theorem 2, this system admits a well-posed Cauchy problem. The value of ζ does not play any role in the establishment well-posedness, but it is really important for the stability of numerical simulations. If ζ < 0, then at some point exponential explosion of the error of C kij will take place. Thus, as it was explained in section II, ζ > 0 is required for long-term and stable numerical simulations.
Finally, note is that the lapse function multiplying the λ and ζ terms in equation (76) may be dropped and the system will still be well-posed. As is discussed in section VII, this may be advantageous in spacetimes, where black holes are involved.

D. Evolution equations of the constraints of the parabolized KST formulation
The modifications introduced to KST in equation (76), result into modifications of the evolution equations of the constraints. These modifications can easily be derived by virtue of equations (A1)-(A6) and the resulting constraint propagation system iŝ where (. . .) stands for the RHS of equations (73)-(75) and where (1) H abij and (4) H ij are given in appendix A. In the equations above, the exact RHS has been written for the evolution of C kij and only the principal parts for the remaining constraint variables. Equation (80) resembles the form of (5) and the evolution of C kij is not left at the mercy of the low order terms as in the KST formulation because the principal part is not zero anymore.
Equations (80) and (81) have obtained the form of MHSP just like the evolution equations of system 2. In particular, equation (80) is parabolic, the uncoupled system (81) is strongly hyperbolic and the coupling terms are of first order. Therefore, the constraint propagation equations with system 2 admit a well-posed initial value problem. In addition, due to the parabolic nature of equation (80), the C kij constraints are damped and smoothed with time.

VI. STABILITY OF FLAT SPACE
In this section the stability of flat space is studied in conjunction with system 1, the first-order ADM and the KST formulations. System 2 will not be considered any further in this work so that the theoretical predictions of this section are made more relevant to the numerical simulations shown in [31]. A flat space stability analysis along with numerical simulations with system 2 will be considered elsewhere.
The stability analysis is carried out by linearizing the evolution equations and the constraints about flat space. In particular small perturbations are applied on all dynamical variables as follows γ ij = δ ij + δγ ij , K ij = 0 + δK ij , D kij = 0 + δD kij . The linearized evolution equations show how these perturbations evolve with time, while the linearized constraint equations serve as a tool for the classification of the individual modes in constraint satisfying and constraint violating. Both the evolution and the constraint equations are studied in the Fourier space in order to reduce the analysis to an eigenvalue problem.
In what follows to make the comparison among the three formulations appropriate the same gauge condition is employed. This gauge condition has zero shift vector and lapse function of the form α = α(γ) with A = ∂ ln α/∂ ln γ > 0. When linearization is carried out about flat space, "1+log" slicing cannot be distinguished from densitized lapse, because A = 1 for both "1+log" slicing and a densitized lapse with σ = 1.
In addition to a specific gauge condition, only perturbations along the x 1 -direction will be considered. The results for any direction are qualitatively the same, but because they are rather complicated, here only a one dimensional analysis is presented.
The linearized evolution equations are (84) Equations (82)-(84) contain all three formulations (namely ADM, KST, system 1) for brevity, because they all have the ADM one as the underlying formulation. Just by choosing the free parameters one can go from one formulation to another. So, the ADM equations are obtained by setting all free parameters to zero, the KST equations are extracted from (82)-(84) by setting {λ = ǫ = φ = θ = ξ = 0} and finally system 1 is obtained from (82)-(84), by setting ρ = ψ = χ = η = 0. Note that one could study the system of equations (82)-(84) as a whole, but this analysis is complicated and will be the subject of a future work.
The linearized constraint equations are The effect of equation (86) is the same as replacing a partial spatial derivative, ∂ x , with iκ and a partial timederivative, ∂ t , with iω. The Fourier transform of equations (82)-(85) is presented in appendix C.
To study flat space stability in conjunction with the ADM, KST and parabolized-ADM (system 1) equations, all one has to do is solve the eigenvalue problem which results by the Fourier transform of (82) and which can be written in the form where the transpose of the column vector v is v T = {γ 11 ,γ 12 ,γ 13 ,γ 22 ,γ 23 ,γ 33 ,K 11 ,K 12 ,K 13 ,K 22 ,K 23 ,   I: Classification of the ADM flat space modes into constraint satisfying and constraint violating. The first column gives the number of modes. The second column gives the characteristic speed, which corresponds to the modes of the first column. The third, fourth and fifth columns have either a "x" or a " " as values. A "x" means that the constraint is not violated by the specific modes, whereas a checkmark " " implies that at least one of the components of the constraint is violated by the specific modes.
One can check whether these modes violate or satisfy the constraints, via pure inspection of equations (C4). If the components of a given eigenvector of M are such that it makes all constraints vanish, the eigenvector is said to be a constraint-satisfying mode. Otherwise, it is called a constraint-violating mode.
The 30 × 30 matrix M and its eigenvectors will not be presented in this work for brevity, but they are available by the author upon request. What will be presented are its eigenvalues ω i or characteristic speeds v i = ω i /κ, where appropriate.

A. Flat space stability with the first-order ADM formulation
When the ADM equations are employed, M has ten non-zero eigenvalues and twenty-six eigenvectors, if A = 1/2, whereas it has twenty-four eigenvectors, if A = 1/2. A classification of these modes into constraint satisfying and constraint violating is presented in table I. No matter what the value of A is, there always are eight constraint satisfying modes (first three rows of table I). In addition to these eight modes, there always are two constraint violating wave modes (fifth row of table I).
Note here that the results shown in table I are valid only for A = 1/2. The case A = 1/2 has to be considered separately. The results with A = 1/2, are qualitatively the same, but as already mentioned matrix M has twenty-four eigenvectors. This 4-fold or 6-fold degeneracy is disadvantageous for the ADM formulation, because these degenerate subspaces lead to linear growth of the constraint violating modes.
To see this one has to understand the structure of the 4-dimensional degenerate subspace (A = 1/2). A particularly useful tool for doing this is a Jordan decomposition (JD) of matrix M . If A = 1/2, one can perform the exact same procedure to study the six-dimensional degenerate subspace.
The JD is a similarity transformation M = XJX −1 , where X is a non-singular transformation matrix and J is called the Jordan matrix, which has the eigenvalues of M on its diagonal. If J is diagonal then M has a complete set of eigenvectors, if however J is not diagonal, then the off diagonal elements (which can be placed either above or below the main diagonal) are equal to unity. These off-diagonal elements can give the information of the combination of dynamical variables, which are responsible for the fact that matrix M does not have a complete set of eigenvectors. This can be seen, if the eigenvalue problem (87) is written as follows Equation (88) is an eigenvalue problem for matrix J. One can find the Jordan matrix J and the non-singular transformation matrix X by standard methods and hence one can also find which components of the transformed vector X −1 v correspond to the off-diagonal elements of J.
By carrying out the JD in the ADM case one finds that the four zero modes, which are thought to be "nonevolving" and which belong to the four-dimensional degenerate subspace, are described by the following equations Plugging in the components of the zero-speed eigen- TABLE II: Classification of the KST flat space modes into constraint satisfying and constraint violating. The first column gives the number of modes. The second column gives the characteristic speed, which corresponds to the modes of the first column. The third, fourth and fifth columns have either a "x" or a " " as values. A "x" means that the constraint is not violated by the specific modes, whereas a checkmark " " implies that at least one of the components of the constraint is violated by the specific modes..
vectors shows that the RHS of (90) and (92) is 0. Hence, components γ 12 , K 13 , γ 13 , K 13 should not evolve with time. This is true only if the RHS is exactly zero. However, small numerical perturbations can result in an nonzero RHS. As a result, K 12 and K 13 will begin to grow linearly in time. Non-zero values of K 12 and K 13 , in turn lead to violations of M 2 and M 3 respectively. This linear growth in the non-linear regime, where the pertur-bations are strong, is likely to obtain explosive behavior and eventually terminate simulations. This is precisely the root of the ill-posed nature of the ADM formulation. So far, even the analysis of flat space-time shows that one has to seek for other formulations, which at least a complete set of eigenvectors. This is exactly what was achieved by the KST formulation.

B. Flat space stability with the KST formulation
If the KST formulation is employed, matrix M of equation (87) has twelve non-zero eigenvalues and a complete set of eigenvectors, irrespectively of the value of A. Nevertheless, it still has sixteen zero eigenvalues which correspond to constraint violating solutions.
A classification of the KST modes into constraint satisfying and constraint violating is given in table II. Note, that the results presented in table II are for the values of parameters given by equation (69).
A few remarks are worth making at this point: a) The components of the constraint satisfying modes of the KST formulation (first 3 rows of table II) are exactly the same as those of the ADM formulation independently of whether A = 1/2 or A = 1/2, as expected.
b) The KST formulation has cured the pathology of the incompleteness of eigenvectors of the ADM formulation. Note further that if χ = −2/5, the number of modes that violate the Hamiltonian constraint are minimized. It is not known yet, whether this is true for arbitrary spacetimes. For this, a detailed study is required, which is outside the scope of this work. However, it is interesting to notice the coincidence between this property and equations (71). c) Despite the desired property of strong hyperbolicity, the KST formulation has eighteen static modes sixteen of which violate the C kij constraints. According to [17], the existence of static modes can lead to exponential explosion caused by low-order term perturbations.
The disadvantages of the KST formulation were demonstrated in [31], where it was shown that numerical simulations with the KST system in conjunction with finite difference methods are unexpectedly terminated in certain cases. In particular, the KST simulations of gauge waves and expanding Gowdy spacetimes were terminated even earlier than the ADM simulations. In addition, linear growth of C kij was observed in the evolution of initially small random noise. The unexpected behavior of the KST formulation was found to be caused by the growth of "static modes", which were excited by the unavoidable small violations of the momentum constraints. The flat space stability analysis with system 1, shows that system 1 is the most stable of the 3 formulations. For the choice of parameters that make system 1 wellposed, i.e. when both conditions (56) and (57) are satisfied, and one of equations (58) or (59) holds true, M has twenty-eight non-zero eigenvalues and a complete set of eigenvectors.
The results presented in this section are specifically for  III: Classification of the flat space modes of system 1 into constraint satisfying and constraint violating. The first column gives the number modes. The second column gives the eigenvalues, which correspond to the modes of the first column. If the eigenvalue is purely real, then the corresponding mode propagates as a wave, whereas if the imaginary part of the eigenvalue is negative, the corresponding mode is getting damped. The third, fourth and fifth columns have either a "x" or a " " as values. A "x" means that the constraint is not violated by the specific modes, whereas a checkmark " " implies that at least one of the components of the constraint is violated by the specific modes.
parameters which satisfy (58) and ζ = 0. The analysis with parameters satisfying (59) and ζ = 0 is qualitatively the same. The reason why ζ is set to zero, is because the results with non-zero ζ are rather complicated. The effect of the ζ term has already been discussed both in section II and in [7] and that is the wavelength-independent damping of all modes which violate C kij . In addition to the choice of parameters described above, a harmonic lapse function (A = 1/2) has been used here for simplicity. The results remain qualitatively the same for any value of A. A classification of the system 1 flat space modes into constraint satisfying and constraint violating follows in table III. Notice that in this case the eigenvalues have been given, instead of the characteristic speeds.
A few remarks can be readily made: a) The eigenvalues of all constraint violating modes are non-zero and in particular they have negative imaginary part, which implies exponential damping. The damping timescale is inversely proportional to κ 2 , which means that the short wavelength perturbations are damped extremely efficiently. b) Another comment regards the nature of the system 1 modes. The first eight constraint-satisfying modes propagate as waves (the eigenvalues are purely real), except for two which are static. These modes are the same among all formulations, which are derived by addition of combinations of constraints to the RHS of the ADM formulation.
All remaining modes are only damped, except for the two of the last row of table III. These two modes represent damped traveling waves. The origin of these modes can be understood, if one takes a closer look at the ADM modes. In table I, it is seen that apart from the six traveling constraint satisfying modes, there are another two traveling wave modes which violate the constraints. The parabolic terms of system 1 are second order (see equations (82)-(84)) and therefore they cannot interact with the first order terms of the ADM formulation and lead to an exact cancellation of the two ADM traveling constraint-violating modes. However, because these ADM modes violate the constraints, system 1 damps them while they propagate. c) Finally, the zero eigenvalues (first row of table (III)) correspond to constraint satisfying modes and thus they do not pose any threat. Nevertheless, a fixed shift β i = constant assigns non-zero speed to these two modes. In particuarly they become equal to β i κ i , where κ i is the one-form which corresponds to the wavevector.
All aforementioned predictions about the damping properties of system 1 have been numerically confirmed in [31].

VII. ON THE CHOICE OF PARAMETERS
This section discusses some ideas, which are related to the choice of parameters of system 1 and 2. It should be noted that this is a complicated discussion and only preliminary arguments will be made which are by no means complete. An analysis of this delicate subject has also been given in [31].
Because the equations have been parabolized, the values of the free parameters depend mainly on two factors: a) The effectiveness of the damping rate and b) the numerical scheme employed to integrate the evolution equations.
As already mentioned in sections II and IV B, the damping rate is of the form exp(−αλκ 2 t). Therefore, what one would desire is for the minimum damping rate to be efficient enough. The minimum damping rate is exp(−α min λκ 2 min t), where α min is the minimum value of the lapse function and κ min = 2π/L with L the length of the computational domain. Therefore, if α min is very small then the damping rate is very small, unless the system 1 parameters are assigned large values.
In principle, setting the damping parameters to very large values results in a very efficient drive towards the surface of constraints. Nevertheless, there are numerical issues involved with parabolic equations which have to be addressed. Before the validity of the following discussion is taken for granted a Von Neuman stability analysis [34] of the evolution equations is required. However, general statements can still be made even if such an analysis is not available.
The numerical stability of an explicit numerical scheme for parabolic equations requires that the mesh size and damping parameter satisfy a Courant condition of the form λ∆t/(∆x) 2 ≤ c, where ∆t is the time step, ∆x is the size of the mesh and c is a constant. For example, c = 1/2 for a one dimensional linear diffusion equation with a four point explicit scheme with forward time difference [34].
The above Courant condition implies that the time step has to be really small, if λ is large. It is this tradeoff between time marching versus strong damping, which is the limiting factor on the values of the damping parameters with explicit numerical schemes.
Similar constraints on the maximum time step exist when explicit schemes are chosen for the ADM evolution equations. In this case, it is the speed of propagation of the information, which sets the limit on the time step. For most applications, numerical stability with the ADM formulation requires a Courant condition of the form µ = ∆t/∆x ≤c, where µ is the usual Courant number. The same is true for other hyperbolic systems. For example, the numerical integration of a one-dimensional wave equation with characteristic speed unity and a standard explicit scheme [34] requires µ = ∆t/∆x ≤ 1.
Systems 1 and 2 contain both parabolic and hyperbolic terms, and hence for numerical stability both the parabolic and the hyperbolic condition have to be satisfied [34]. Combining the two conditions, one finds that for system 1 in conjunction with explicit schemes the overall numerical stability is of the form λµ ≤c∆x. Thus, λ should have a value, such that µ and ∆x should be large enough in order to finish computations within reasonable time-scales, while maintaining satisfactory damping properties. This has been shown to be possible in [31].
Note, however, that the above limitation on the timestep is not present when the numerical integration is carried out with implicit schemes. For example, a Crank Nickolson numerical scheme for the one-dimensional diffusion equation [34] is always stable. An implicit scheme would allow one to set the damping parameters to arbitrarily large values without having to worry about limitations on the grid structure, but at the cost of cumbersome programming.
Taking all these facts into account, it is concluded that until more detailed numerical simulations become available the safe choice for the parameters of systems 1 and 2 are those which satisfy the conditions for well-posedness presented in sections IV A and V C 1 respectively.

VIII. CONCLUSIONS
Two new parameterized 3+1 formulations have been presented. These two formulations are referred to as system 1 and system 2.
The well-posed system 1 is derived by adding combinations of the constraint equations and their derivatives to the RHS of the first-order ADM formulation. As a result the evolution equations of the dynamical variables have become second order in spatial derivatives and it has been shown that system 1 has obtained mixed hyperbolic -second-order parabolic nature.
The evolution of the constraints with system 1 is governed by 2nd-order parabolic PDEs irrespectively of the gauge condition employed. This in turn means that the constraints are damped and smoothed in time, and hence the surface of constraints is a local attractor.
The well-posed system 2, is a parameterized parabolization of what is usually called the KST formulation. The original KST gauge condition was generalized and a certain combination of the terms of system 1 was added to the RHS of the KST formulation. System 2 is a manifestly mixed hyperbolic -second-order parabolic system of equations.
The effect of the terms introduced in system 2 is the efficient damping and smoothing of perturbations which violate a subset of the constraints of the formulation, while the remaining constraints evolve according to strongly hyperbolic equations.
As a first test for system 1, a study of the stability of flat space against small amplitude perturbations was carried out and contrasted to the stability of the ADM and the KST formulations. The results indicate that in this case system 1 is more stable than both the ADM and the KST formulations. In particular, it was shown that system 1 damps exponentially and smooths all finite wavelength constraint-violating modes. Numerical support for this prediction has been given in [31].
These two new formulations provide a new basis for obtaining well-posed formulations of GR which possess constraint damping and in turn a new basis for constructing numerical schemes for long-term and stable numerical simulations of the Einstein equations. the possibility of using higher order systems of partial differential equations in numerical relativity. P can be considered to coincide with its linearized and localized principal part [29], which is denoted by P .
Before the proposition of this appendix and its proof are presented it is important to introduce a definition and a theorem from [29]. a. Definition 3: An operator, P , is called semibounded if there exists a Hermitian operator such that (u, P u) H + (P u, u) H ≤ 2c(u, u) H , for all sufficiently smooth functions u(x) and c ∈ R is a constant.
The quantity (u, w) H above is called an energy inner product and its exact definition can be found in [29].
b. Theorem 3: If the operator P is semi-bounded the Cauchy problem for (B1) is well-posed. Now the following proposition can be proven. Proof: Proving well-posedness in the case κ = 0 is always trivial because then P (iκ) = 0 [29]. For κ = 0, if P (iκ) satisfies the Petrovskii condition then the eigenvalues ω i of P (iκ) satisfy where c ∈ R and is independent of κ.
If P (iκ) has a complete set of eigenvectors, then the matrix S(κ)P (iκ)S(κ) −1 is diagonal with the eigenvalues of P (iκ) on the diagonal, provided that S(κ) −1 contains these eigenvectors as columns. Now, the hermitian matrix H(κ) = S * S, where S * is the complex conjugate transpose of S, satisfies whereω i denotes the complex conjugate of ω i . The definition of inequalities between Hermitian matrices can be found in [29]. The energy inner product (u, w) H exists because of assumption c) [29] of proposition 1. Then, the following energy estimate follows (u, P u) H + (P u, u) H = = (P (iκ)ũ, H(κ)ũ) + (ũ, H(κ)P (iκ)ũ), = ũ, (H(κ)P (iκ) + P (iκ) * H(κ))ũ , whereũ is the Fourier transform of u. Parseval's relation, (u, w) H = (ũ,w) H , has been used in lines 1 and 3 and inequality (B3) has been used in going from line 2 to line 3 above. Inequality (B4) proves that P is a semi-bounded operator. Then according to theorem 3, (B1) admits a wellposed initial value problem. This completes the proof of proposition 1.
APPENDIX C: FOURIER TRANSFORM OF LINEARIZED EVOLUTION AND CONSTRAINT EQUATIONS OF ADM, KST AND SYSTEM 1 apparent discrepancy arises, because the first-order ADM evolution equations are not reduced to equations for the three-metric only. This is exactly what was done in [17], and hence the analysis presented in [17] was for a different set of equations than those of the first-order ADM formulation.
The reason why this approach is not followed in this work is that the totality of the modes present in the first-order ADM formulation is larger than that of the second-order ADM formulation. The first-order ADM formulation has opened up an eighteen-dimensional sub-space, where more instabilities than the standard ADM formulation may arise. The totality of solutions of the first order ADM is greater than that of the standard second-order ADM and the two coincide, if and only if the C kij constraints are satisfied [32]. The proposition in this work is to simulate the Einstein equations using variants of the first-order ADM formulation. Therefore, a study of the totality of the modes present in the firstorder ADM formulation and its variants is required and the evolution equations will not be reduced to equations for the three metric only.