Evolution equations for a wide range of Einstein-matter systems

We use an orthonormal frame approach to provide a general framework for the first order hyperbolic reduction of the Einstein equations coupled to a fairly generic class of matter models. Our analysis covers the special cases of dust and perfect fluid. We also provide a discussion of self-gravitating elastic matter. The frame is Fermi-Walker propagated and coordinates are chosen such as to satisfy the Lagrange condition. We show the propagation of the constraints of the Einstein-matter system.


Introduction
Einstein's theory of General Relativity provides us with the most appropriate tool for studying the dynamics of self gravitating objects. It is therefore of clear interest to study the structural properties of the Einstein field equations and to provide a framework for studying their solutions. The Cauchy problem (or, initial value problem) provides a setting for the analysis of generic solutions to the field equations parametrised in terms of the initial conditions -for details, see [7,16,23]. In particular, one is interested in showing that the Einstein equations admit a wellposed initial value formulation [25]. The standard strategy to address this issue is to show that the Einstein equations imply evolution equations that are on a hyperbolic form. Physical considerations associated to causality lead to the expectation of the Einstein equations admitting a hyperbolic formulation despite the fact that the immediate form of the equations is not manifestly hyperbolic due to general covariance. Thus, it is necessary to find a subset of the Einstein equations which indeed admits hyperbolicity. This procedure is called hyperbolic reduction -see [20] and [14] for details; for an overview of the different reduction methods, see [24].
The well-posedness of the vacuum Einstein equations was first shown in [13] -and later in the case for dust and the Einstein-Euler by the same author [12]. These results were obtained using a harmonic gauge to reduce the field equations to a form which is second order hyperbolic (Leray hyperbolicity). In [2] this method is extended to show existence of solutions locally for a selfgravitating, relativistic elastic body with compact support. Furthermore, in [9] well-posedness of a viscous fluid coupled to the Einstein equations is presented and in [11] a viable first order system is constructed. In [17] the concept of first order symmetric hyperbolic (FOSH) equations was developed. The same author showed later [18] that the Einstein-Euler system could be put on a FOSH form. In [15] a different approach, which makes use of a formulation in terms of frame fields, is employed to construct evolution equations for the Einstein-Euler system which also are on the form of a FOSH system. This method has the advantage that the reduced equations are symmetric hyperbolic while still maintaining a Lagrangian form -which is important in order to keep track of a boundary in the case of matter distributions with compact support.
As is made clear by the discussion above, it has been customary up until now to apply the hyperbolic reduction procedure to individual matter models separately -i.e. for every particular type of energy-momentum tensor. The motivation for our study is provided by the observation that the energy momentum tensor for a perfect fluid, elastic matter -see [4,19,3] for detailsand bulk viscosity -e.g. see [5,21,10,6,22] and references therein-may be put on a form consisting of a part involving the 4-velocity u and energy density ρ and a part involving a spatial symmetric tensor Π. Thus, by "hiding" the specific matter variables in the tensor Π one cannot differentiate between elastic matter, perfect or viscous fluid by considering the energy momentum tensor alone. By employing a hyperbolic reduction of the Einstein field equations coupled to an energy-momentum tensor on such a general form, we provide the necessary conditions for such a matter model to form FOSH evolution equations. We show that one can avoid the details of the specific matter models in the construction of a FOSH system by introducing an auxilliary field. Furthermore, we provide a constraint on the Π on the form of a wave equation, which must be satisfied for the system to be put on a FOSH form.
The procedure we employ to obtain these evolution equations is similar to that of [15] and may be described as follows: we introduce a frame field to replace the metric tensor as a variable and fix the gauge by choosing Lagrangian coordinates -i.e. one of the vectors of the frame field is chosen as to coincide with the 4-velocity of the particle trajectories; we also let the rest of the frame be Fermi propagated. By virtue of the Bianchi identity and assuming the connection to be Levi-Civita we show that the solution to a set of new field equations constructed with so called zero-quantities implies the existence of a metric solution to the Einstein field equations. A subset of these equations provides the symmetric hyperbolic evolution equations. As part of this construction, it turns out to be necessary to introduce an auxilliary field to remove derivatives of the energy-momentum tensor from the principal part of the evolution equation of some of the geometric fields. The evolution equation of Π -which encodes the matter fields-is given in terms of the electric decomposition of the auxilliary field. Finally, we make use of the evolution equations, Cartan's identity and the Bianchi identities to show the propagation of constraints. It is important to stress that due to the generality of the procedure, we do not provide an equation defining ρ. It is therefore necessary to provide an equation of state (or the equivalent) when using our equations for a specific matter model. We treat dust and perfect fluid as examples at the end and briefly discuss elastic matter.
A limitations of our procedure is in requirement of Π being a purely spatial tensor -indeed, without this requirement the energy momentum tensor would take its most general decomposition form. The difficulty of allowing Π to have timelike components resides in the procedure of keeping the hyperbolicity of the theory. We have used the spatial property of Π extensively in the process of eliminating problematic derivative terms from the principal part of the equations. We also assume that the equations of motion for a matter system may be entirely determined by the divergence-free condition of the energy-momentum tensor. Thus, any matter models which require additional equations to close the evolution of the matter variables, are not considered herein.
Lastly, we should mention that a very good discussion of the Einstein-Euler-entropy system is found in [8] where a complete discussion of the arguments of the framework put forward in [15] is given.

Overview of the article
In Section 2 we introduce the geometric tools necessary for the subsequent discussion; in particular the frame formalism is introduced. The Einstein equations together with the energy momentum tensor is presented in Section 3; we also give a brief review of the projection formalism. In Section 4 we outline the gauge choices and in Section 5 we present the zero-quantities used in the propagation of constraints. The evolution equations for our system are derived in Section 6. In Section 7 we show propagation of constraints and in Section 8 we present the reduced equations for the special cases of dust and perfect fluid. A brief discussion of self-gravitating elastic matter is also presented. Final remarks are given in Section 9. An appendix provides an extended discussion of a framework for relativistic elasticity -this model provided the main motivation for the present analysis.

Geometric background
In what follows, let (M, g) denote a spacetime represented by a 4-dimensional manifold, M, with a Lorentzian metric g. The motion of particles of some matter filling spacetime give rise to a natural splitting by constructing frames comoving with the flow lines of the particles. One advantage with such a view is that it does not require a foliation. We shall denote the tangent vector to the flow lines as u satisfying g(u, u) = −1.
At each point p ∈ M the frame field {e a } is such that g(e a , e b ) = η ab .
The frames {e a } give rise to a co-frame, {ω a } satisfying In the following all indices will be given in terms of the frame and co-frame unless otherwise stated. The metric tensor give rise to a natural connection ∇ such that ∇g = 0, which is the metric compatibility condition. In terms of the frames, this condition takes the form where the frame connection coefficients are defined by the directional derivative along the direction of the frame indices Furthermore, if the connection ∇ is torsion-free, we have that where the frame components of the torsion tensor are defined by The commutation of the connection may be expressed in terms of the Riemann curvature tensor and the torsion tensor The frame components of the Riemann curvature tensor is given by -see [20] for details. The Riemann tensor has all the usual symmetries, and it satisfies the Bianchi identity for a general connection Furthermore, we recall that the Riemann tensor admits the irreducible decomposition with C c dab the components of the Weyl tensor and Rη ab (7) denotes the components of the Schouten tensor. The connection ∇ is called the Levi-Civita connection of g if it satisfies (1) and (2). In what follows we will assume the connection to be Levi-Civita.

The Einstein equations
In this work we consider the Einstein equations with energy-momentum tensor on the form where ρ is a positive function of the matter fields. We require Π ab to be a symmetric and purely spatial tensor -i.e.
We do not put any further restrictions on Π ab other than it satisfies the divergence-free condition of (9) ∇ a T ab = 0.
Remark 1. An energy momentum tensor of the form given in (9) is of a very general form and the conditions (10a), (10b) are not stringent restrictions. Thus, the power of the formalism developed herein lies in its generality: given an equation for ρ in terms of the matter fields, one can ignore the matter specific equations of motion and instead solve equations for Π ab . The equations obtained will then be symmetric hyperbolic. This assumes that one can extract the complete set of equations of motion for the matter fields from (11).

A projection formalism
At each point in the spacetime manifold M the flow lines give rise to a tangent space which can be split into parts in the direction of u and those orthogonal. This means that without implying a foliation, we may decompose every tensor defined at each point p ∈ M into its orthogonal and timelike part. This may be done by contracting with u and the projector defined as Thus, a tensor T ab may be split into its time-like, mixed and space-like parts given, respectively, by where ′ denotes that the free indices left are spatial -e.g. T ′ a0 u a = 0. Decomposing ∇u we obtain where χ a b and a b are the components of the Weingarten tensor and 4-acceleration, respectively, defined by In the literature (e.g. see [25] p.217) the trace, trace-free and antisymmetric part of (12) is called, respectively, the expansion, shear and the twist of the fluid. By decomposing (11) we obtain an equivalent system of equations in terms of the above quantities The decomposition of the four volume is Given a tensor T abc which is antisymmetric in its two last indices, we may construct the electric and magnetic parts with respect to u. In frame indices this is, respectively, defined by where the Hodge dual operator, denoted by * , is defined by and has the property that T * * abc = −T abc . Depending on the symmetries and rank of the tensor, the above definition for electric and magnetic decomposition may vary slightly. Central for our discussion is that E ab and B ab are spatial and symmetric.

Gauge considerations
The gauge to be considered in our hyperbolic reduction procedure for the Einstein field equations follows the same considerations as in [15]. In particular, we make the following choices: i. Orientation of the frame. We align the time-leg of the frame with the the flow vector u tangent to the worldlines of the particle -that is, we set ii. Basis in a coordinate system. Given a coordinate system x = (x µ ) we expand the basis vectors as e a = e a µ ∂ µ .
Given an initial hypersurface, S ⋆ , then the coordinates (x j ) defined on S ⋆ remain constant along the flow and, thus, specify the frame.
iii. Lagrangian condition. The implementation of a Lagrangian gauge is equivalent to requiring that e 0 = ∂ t where t is a suitable parameter along the world-lines of the material -e.g. the proper time. In terms of the components of the frame, this condition is equivalent to requiring that iv. Fermi Propagation of the frame. We require the vector fields e a to be Fermi propagated along the direction of e 0 -i.e.
This implies the following conditions on the connection coefficients: for i, j = 1, 2, 3. The second condition is a consequence of the metric compatibility condition. A frame satisfying the above equation is a frame where e 0 = u and {e i } is orthonormal at every point along the trajectory for which u is the tangent vector.

Zero-quantities
In the subsequent discussion it will prove convenient to introduce, as a book-keeping device, the zero-quantities where L ce denotes the components of the Schouten tensor as defined by equation (7) and Π ≡ Π a a . Moreover, byR d abc it is understood the expression for the Riemann tensor in terms of the connection coefficients Γ a b c and its frame derivatives. We have also defined whereĈ d abc is defined as having the same symmetries as the components of the Weyl tensor C d abc . Remark 2. The components ρ d abc are known as the algebraic curvature and encode the decomposition of the Riemann curvature tensor in terms of the Weyl and Schouten tensors while F c abc are the components of the Friedrich tensor. The latter provides a convenient way to encode the second Bianchi identity for the curvature.
Remark 3. The tensor Z cab , hereafter to be referred to as the Z-tensor, is introduced in order for the evolution equations of the electric and magnetic part of the Weyl tensor to be expressed in terms of lower order terms -i.e. preventing any derivatives of Π ab to appear in the equations and hence keeping their hyperbolicity.
In terms of the objects introduced in the previous paragraphs, the Einstein field equations (8) can be encoded in the conditions More precisely, one has the following result: implies the existence of a metric g solution to the Einstein field equations (8) with energy-momentum tensor defined by the components T ab . Moreover, the fieldsĈ d abc are, in fact, the components of the Weyl tensor of g.
Remark 4. Note that equations (19a)-(19d) do not provide a closed system of evolution equations for the unknowns of our system. They are only the necessary equations for giving Lemma 1.
Proof. The frame {e a } obtained from the solution to equation (19b) implies, in turn, by the condition ω b , e a = δ a b the existence of a coframe {ω b } from which one can construct a metric tensor g via the relation g = η ab ω a ⊗ ω a .
Since the coefficients Γ a c b satisfy the no-torsion and metric compatibility conditions (19b) and (1), then they must coincide with the connection coefficients of the metric g with respect to the frame {e a }. Moreover, by equation (3) we have that where R d abc denotes the frame components of the Riemann curvature tensor. Using the Riemann decomposition as defined by equation (6) together with equation (19c) we obtain Taking the trace of equation (20) with respect to the indices b and d and using the trace-free property of the Weyl tensor andĈ d abc we obtain Finally, taking the trace of equation (21) and using equations (11) and (5), we get the identity, The latter shows thatL ab are, in fact, the components of the Schouten tensor of the metric g.
Using the definition of the Schouten tensor in terms of the Ricci tensor, equation (7), it follows readily that the metric g satisfies the Einstein field equations with an energy-momentum tensor defined by the components T ab . Returning to equation (21) we conclude by the uniqueness of the decomposition of the Riemann tensor that the fieldsĈ d abc are, in fact, the components of the Weyl tensor of g.
Remark 5. In the following to ease the notation, and in a slight abuse of notation we simply write C d abc instead ofĈ d abc .

Evolution equations
Given the gauge conditions introduced in Section 4, the next step in our analysis involves the extraction of a suitable (symmetric hyperbolic) evolution system from equations (19a)-(19d). We do this in a number of steps.

Equations for the components of the frame
The evolution equations for the components of the frame e a µ are obtained from the no-torsion condition (19b). In order to do so we exploit the freedom available in the choice of the frame and require it to be adapted to the world-lines of the material particles and the gauge conditions outlined above.
Making use of the expansion (14) in equation (19b) one readily finds that Setting a = 0 in the above expression and making use of the Lagrangian gauge condition (iii) we obtain This last equation will be read as an evolution equation for the frame coefficients e b ν with b = 1, 2, 3. As it only contains derivatives along the flow lines of the matter, it is, in fact, a transport equation along the world-lines. Observe that for b = 0 the equation is satisfied automatically -recall that as a consequence of the Lagrangian condition (15) the coefficients e 0 µ are already fixed.
Remark 6. Assuming that the gauge conditions (i), (ii) and (iii) above hold, equation (22) can be succinctly written as Σ 0 c b = 0. This observation will be of use in the discussion of the propagation of the constraints.

Evolution equations for the connection coefficients
The evolution equations for the frame components are given in terms of the frame connection coefficients. Due to the Fermi propagation and the metric compatibility, equation (1), the independent, non-zero components of the connection coefficients are Γ i k j , Γ 0 0 j and Γ i 0 j . Evolution equations for Γ i k j may be extracted from the equation for the algebraic curvature (19c). More precisely, we consider the condition The Riemann tensor can be expressed in terms of the connection coefficients via equation (3). Furthermore, using equation (18a) and the gauge condition (iv), we obtain where i, j, k, . . . = 1, 2, 3. In the above calculation we have used that Π a0 = 0 and η i 0 = 0.

Remark 7.
Assuming that the gauge condition (iv) holds, equation (23) is equivalent to Observe again, that the resulting equations are, in fact, transport equations along the world-line of the material particles.
The evolution equations for the remaining connection coefficients will be obtained by splitting Π ab into its trace and trace-free part, where Π {ab} denotes the trace-free part of Π ab . Plugging this into (13a) and (13b), we obtain The last term may be written, where we have used the symmetry of the Riemann tensor in the last step. A straight forward calculation shows that, Using that, Substituting this back into (25), we may now write the {0, i} components of J db as In the above expression we have used the Lagrangian gauge condition to set u i = 0. Finally, using the definition for a i , we readily obtain To obtain the evolution equation for the remaining connection coefficient, we consider the {i, j} components of J ab : From equation (3) we have, Substituting this back into (29) gives the final evolution equation: Equations (28) and (30) are on a form which is known to be symmetric hyperbolic -we refer again to [1] for details.
Remark 8. The presence of a spatial derivative of ρ in equations (28) and (30) means that an equation for ρ is necessary to ensure the hyperbolicity of the equations.

Evolution equations for the decomposed Z-tensor
It is well known that in vacuum the Bianchi equation leads to a symmetric hyperbolic equation for the independent components of the Weyl tensor. By contrast, an inspection of the definition of the Friedrich tensor F abcd , equation (17b), reveals that the condition F abc = 0 involves both derivatives of C abcd and the matter variables. This potentially destroys the symmetric hyperbolicity of the equation for the components of the Weyl tensor. In the following we will show that it is possible to deal with this difficulty by providing two auxiliary fields -the Z-tensor and σ-tensor as defined by equations (18c) and (6.3), respectively.
We first define some important quantities and identities used in the following discussion. The Z-tensor has the symmetries The symmetry of the Z-tensor thus allows for a decomposition in terms of its electric and magnetic parts defined respectively as where, Z * ebd , is the dual Z-tensor defined in the customary way. The electric and magnetic part of the Z-tensor are symmetric tensors defined on the orthogonal space of u -i.e. one has that Ψ ac = Ψ (ac) , Ψ ac u a = 0, Φ ac = Φ (ac) , Φ ac u a = 0.
As such, the Z-tensor and its dual may be expressed in terms of the spatial fields By plugging the definition for the Z-tensor into the definitions of Ψ ab and Φ ab , respectively, we obtain an evolution equation for the matter tensor Π ab in terms of Ψ ab together with a constraint equation. Namely, one has that where D b denotes the Sen connection defined as, It is worth noting that due to the 1 + 3 split of space time, we do not have a spatial metric on the 3-dimensional hypersurfaces. Hence, we cannot define a spatial derivative -i.e. a spatial metric satisfying the metric compatibility condition does not exist on the three surfaces. Equation (32b) is regarded as a constraint equation.
Remark 9. Note that equation (32b) will always hold as long as the definition of the Z-tensor (i.e. equation (18c)) propagates. This will be showed in Section 7.
In order to close the system and to ensure hyperbolicity a set of evolution equations for the fields Ψ ab and Φ ab are needed. In the rest of this section we shall develop these equations and show they form a first order symmetric hyperbolic system. The evolution equations for Ψ ab is obtained by taking the divergence of the Z-tensor -i.e we have the equation Expanding the above equation and using the decomposition of the Z-tensor -i.e. equation (31a)-we obtain after a number of steps the equation where W ab denotes the lower order terms and is explicitly given by In the above, we have defined Note that the derivatives of χ ab and a a may be expressed in terms of the connection coefficients and thus dealt with by (23) and the definition for the Riemann tensor -i.e. we have In obtaining equation (33) we have used standard tensor manipulations involving the commutation of derivatives using the Riemann tensor and frequently making use of the spatial property of Π ab to get rid of derivatives. In particular, we have used the identities where L.O.T is a shortening for "lower order terms". The first is obtained from considering ∇ b Z ′ c0b and the second is a trivial result of the definition for the Z-tensor.

Remark 10.
To preserve the hyperbolicity of equation (33) it is understood that the tensor σ ab is given in terms of lower order terms. It can readily be shown that σ 00 = L.O.T from the temporal part of equation (33); with the symmetry of σ ab , it is thus 9 components which need to be specified: To obtain the evolution equation for the field Φ ab , we proceed in a similar way as above by considering the dual equation, By applying the decomposition (31b) and expanding we obtain after a few manipulations the evolution equation on the desired form, with U ac denoting the lower order terms -explicitly given by Equations (33) and (34) are on a form known to be symmetric hyperbolic -see [1] for a more detailed discussion. We note that we have made ample use of the suite xAct 1 to obtain W and U.
Remark 11. In the expressions for W ac and U ac it is understood that wherever R d abc appears, it is to be evaluated using the decomposition in terms of the C d abc andL ab .

Evolution equations for the decomposed Weyl tensor
The construction of suitable evolution equations for the components of the Weyl tensor follows a similar approach as in the previous discussion. Again, the strategy is to decompose the Weyl tensor into parts orthogonal to the 4-velocity -i.e. one needs to understand the form the Weyl tensor takes on the orthogonal space to the 4-velocity.
Due to the symmetries of the Weyl tensor, the essential components are encoded in what are called the electric and magnetic parts of the Weyl tensor defined, respectively, as where C * abcd denotes components the Hodge dual of the Weyl tensor. In terms of these spatial tensors, the frame components of the Weyl tensor and its dual admit the decomposition -see e.g. [20] for details. For convenience we have written In order to obtain evolution equations for E ab and B ab , we make use of the following decomposition of the Bianchi identity (19d) and its dual: The term F a0b = −F ab0 in equation (35a) gives the evolution equations for E ab . More precisely, after a long computation one finds that Similarly, the symmetric part of the term F * ab0 in equation (35b) gives the evolution equations for B ab -namely, We note that equations (36) and (37) are on the same form as the one given in [15] and constitutes a symmetric hyperbolic system of equations. We refer once again to [1] for an explicit discussion.
Remark 12. The standard approach to show that equations (36) and (37) constitute a symmetric hyperbolic system ignores the tracefreeness of the fields E ab and B ab and list 12 components in a vector. Thus, a posteriori it is necessary to show that the fields are tracefree if they were so initially. This is discussed in Section 8.

Summary
We summarise the results of the long computations of this section by the following proposition: Proposition 1. The evolution equations for the matter fields as expressed by Π ab is given by (32a) and satisfy the constraint (32b). Furthermore, the evolution equations for the geometric fields e a µ , Γ a b c , E ab , B ab and the auxiliary fields Ψ ab , Φ ab are given, respectively, by The above evolution equations constitutes a symmetric hyperbolic system. The remaining equations from (19a)-(19d) are considered constraint equations.

Propagation equations
In order to complete our analysis of the evolution system, we need to show that the equations that have been discarded in the process of hyperbolic reduction (i.e. the constraints) propagate. In this section we will, therefore, construct a subsidiary system for the zero-quantities Σ e ab , ∆ d cab and F abc . The task is then to show that either the Lie derivative of the constraints vanish, or that it may be written in terms of zero-quantities. A key observation in this strategy is the fact that several of the zero-quantities can be regarded as differential forms with respect to a certain subset of their indices -thus, Cartan's identity can be readily be used to compute the Lie derivative in a very convenient way.

Remark 13.
In what follows one should be careful when evaluating the covariant derivative of a tensor fields in frame coordinates. The following order should be employed: first, evaluate the tonsorial expression for the derivative of the tensor, then write the expression in a frame basis, and lastly do any contractions if necessary -e.g contracting with the four velocity.

Propagation of Divergence-free condition
The divergence free condition gives rise to two equations: on the one hand, equation (62b) is an evolution equation for ρ and the other hand equation (62a) a constraint for Π. As such, the latter needs to be shown to hold on the whole space time if satisfied on an initial hypersurface.
We first define the zero quantity, As is obvious from the above, Q b = 0 must hold for the Einstein equations to be satisfied. By contracting with u it is readily shown that Q 0 = Q b u b = 0. Thus it is sufficient to only consider Q i in what follows. A simple calculation shows that, where we have used the commutation property of the connection followed by the definition of the Z-tensor, as well as, 2R e [ad]a Π e a = 0.
By using equation (26) and multiplying through with u d , followed by applying equations (17a) and (27), we obtain the propagation equation, Note that ρ j k 0k = 0 due to the divergence free property of Weyl and T 0i = 0 as a consequence of the gauge. We have also made use of the evolution equation Σ 0 c a = 0.

Propagation equations for the torsion
For fixed value of the index e, the torsion Σ a e b can be regarded as the components of a 2-form -namely, one has that Σ e ≡ Σ [b e c] ω b ⊗ ω c . Using Cartan's identity to compute its Lie derivative along the vector u one finds that The second term in the right-hand side of the above equation can be seen to vanish as a consequence of the evolution equations (cf. Remark 6) while the first one involves the exterior derivative of the torsion which can be manipulated using the general form of the Bianchi identity. For clarity, these computations are done explicitly using frame index notation.
Following the general discussion given above consider the expression ∇ [0 Σ a c b] which roughly corresponds to the first term in the right-hand side of equation (40). Expanding the expression one readily finds that in a different way using the general expression for the first Bianchi identity (i.e. the form this identity takes in the presence of torsion): Setting a = 0 and making use of the zero-quantity defined in (17a) to eliminate the components of the Riemann curvature tensor one finds that where we have used the fact that ρ d [cab] = 0 and the evolution equations (38a) and (38b) in the last step. From the above discussion it follows the propagation equation Remark 14. The main structural feature of equation (41) is the fact that it is homogeneous in the zero-quantities Σ a c b and ∆ d abc .

Propagation equations for the geometric curvature
Next we turn to equation (19c). For this we observe that the zero-quantity ∆ d cab for fixed values of d and c can be regarded as the components of a 2-form on the indices a and b. Using again Cartan's identity one finds that Now, the last term in the right-hand side vanishes due to the evolution equation for the connection coefficients (see Remark 7), while the first term takes the form As in the case of the torsion, the strategy is to rewrite this expression in terms of zero-quantities only. For convenience in the following calculations we set From equations (5) and (18a) it readily follows that To simplify the calculations, we multiply by ǫ l abc . The first term yields while the second term gives where we have made use of the definition (42) and the symmetry of the Schouten tensor. Now, from equations (17b), (18b) and (19d) we readily obtain that Plugging this result back into (44), we obtain Putting the result for calculation (45) and (43) together, we obtain where we made use of equation (4). Multiplying by ǫ l mnp , we recover the equation in its original form. Thus, the right-hand side of the propagation equation for ∆ d c[ab] is given by But we also have that, Plugging the above result back into equation (46), we obtain the final propagation equation Remark 15. As in the case of the propagation equation for the torsion the main conclusion of the previous discussion is that the propagation equation for the zero-quantitity ∆ d abc is homogeneous on zero-quantities.

Propagation of the N-tensor
It is also necessary to show that N abc -see equation (17c)-propagates. The strategy will be different than what has been employed in the above discussions; rather, we will follow the strategy employed for the propagation of the Friedrich tensor in [14].
In the subsequent discussion we shall make use of the observation that, respectively are equivalent to the evolution equations for Ψ ab , Φ ab and Π ab and the constraint equation as given in (32b). Furthermore, we define the fields, By decomposing N * abc in terms of the fields λ ab and ξ ab , we have Using the symmetry relation we obtain the expression λ ab = N ′ ba0 − N ′ ab0 . But from the evolution equation for Π ab , we have that N ′ ba0 = 0, thus λ ab = 0.
Applying the above result, and the divergence in equation (48) we obtain, To obtain the above we have used the evolution equation for Φ ab and multiplied through with the projector h d a to get rid of a divergence. This is permitted as the field ξ ab is spatial. Thus, we have established the following lemma: Lemma 2. If the constraint ξ ab = 0 -equivalently equation (32b) -holds initialy, and under the assumption that the evolution equations (32a) and (38g) holds everywhere on M, then the relation, also holds everywhere on M.

Propagation equations for the Bianchi identity
Lastly, we need to show propagation of the Bianchi identity, equation (19d). Again the strategy is to use the decomposition of the Friedrich tensor and its dual and use the divergence to obtain propagation equations for the constraints. First, we shall express the divergence of F abc in terms of known zero quantities. We have, where we have used the antisymmetry property of the Weyl tensor about the indices a and b and the commutator as defined in Substituting from the definition (17a) we obtain after a lengthy and somewhat messy calculation that where q a is the zero quantity defined in section 7.1.
Following the strategy outlined in [14], we define the fields, which encodes the information of the constraint equations of F abc and F * abc , respectively. Thus, the aim is to find evolution equations for p a and q a .
In terms of the above fields the decomposition (35a) takes the form, where we have used, To obtain the above, we used the evolution equations -i.e. that F bc0 = 0 and F * (bc)0 = 0 as well as the identity, which is a direct result of the symmetry properties of F abc . By taking the divergence of the first index of (50) and equating with u d h a c times (49), we obtain a propagation equation for the p a field, namely In the above, we have used that Σ 0 d b = 0 everywhere on M. Applying the same procedure to (35b) we obtain the propagation equations for the q a field, Remark 16. Again, the main observation to be extracted from the previous analysis is that equations (51) and (52) are homogeneous in the various zero-quantities. Moroever, their form is analogous to that of the evolution equations (36) and (37). Thus, it can be verified they imply a symmetric hyperbolic system. Note also that equation (51) is different in the principle part compared to equation (52). This is due to thhat the evolution equation F ab0 = 0 is not symmetrized. It is also understood in equation (51) that one apply equation (39) to eliminate the time derivative of Q a .

Main theorem
The homogeneity of the propagation equations for the various zero-quantities implies, from the uniqueness of symmetric hyperbolic systems that if the zero-quantities vanish on some initial hypersurface S ⋆ then they will also vanish at later times. We summarise the analysis of the previous subsections in the following statement: to the system of evolution equations given, respectively, by equations

Matter models
We will in the following exemplify the previous discussion with a number of particular matter models. We shall also note that although the equations given in the following resembles those found in [15], the treatment of the propagation of constraints for dust or perfect fluid was not treated therein. We fill this gap in this paper.

Dust
The simplest case is of course that of dust. In this case Π ab = 0 and the expression for the energy-momentum tensor, equation (9), reduces to Furthermore, as there are no internal interactions, each dust particle follows a geodesic -i.e the following hold a a = 0, and equation (23) takes the form Also, we have that Z abc = 0. Thus, the discussion of the Z-tensor and its evolution equations are irrelevant -i.e. there is no need for the construction of an auxiliary field. The evolution equations for the Weyl tensor reduce to and, Thus, equations (22), (54), (55), (56) and (53) provide the symmetric hyperbolic evolution equations for the fields e µ a , Γ i a b , E ab , B ab and ρ, respectively.

Perfect fluid
Before we discuss the details of a perfect fluid, we shall briefly review some important quantities in relativistic thermodynamics. Given a material with N different particle species, n A denote the number density of a particular species, where A = {1, 2, ..., N }. Furthermore, we denote by s the entropy density. The energy density of the system is a function of these quantities -i.e. we have ρ = f (s, n 1 , n 2 , ..., n N ) .
The function f is called the equation of state of the system. Finally, the first law of Thermodynamics is given by, where, denotes the temperature and chemical potential, respectively. In what follows we shall consider a simple perfect fluid -i.e a fluid of only one type of particles (A = 1) and with an energy momentum tensor with Π ab = ph ab .
Consequently, we have where p denote the pressure and is defined by Throughout we shall assume an equation of state of the form given by (57) with A = 1 and the law of particle conservation -i.e. u a ∇ a n = −nχ.
With these assumptions, equations (24a) and (24b) reduce to the well known Einstein-Euler equations, given by It follows from equations (60), (62b), (61) and (58) that the fluid is adiabatic -i.e we have u a ∇ a s = 0.
From the above discussion it follows that equation (28) takes the form Similarly, equation (30) takes the form where we have defined s a ≡ ∇ a s, n a ≡ ∇ a n.
The corresponding evolution equations are obtained by using equations (63) and (61), Now, writing equations (32a) and (32b) in terms of the above definitions we obtain where we have used equation (60) and the definition of µ to obtain,

Concluding remarks
As stressed previously, we have developed first order symmetric hyperbolic evolution equations for a wide range of matter models which solves the Einstein equations. Our formalism should thus be applicable to the development of a theory of Neutron stars as an relativistic elastic system. In this case one would proceed as with the perfect fluid case: one need to write the tensor Π in terms of its trace and trace-free part and provide equations for n and ǫ to close the system. The latter is likely obtained from thermodynamical considerations. In addition it is necessary to provide an equation for σ in terms of lower order terms. The treatment given in this paper is sufficiently general that showing symmetric hyperbolicity for a given matter model coupled to the Einstein equations, is reduced to the simple task of showing that the system admits an energy momentum tensor on the form (9) satisfying (10b) and (10a). It is understood that an equation for ρ and σ is provided.
A The frame components of elastic energy-momentum tensor In the following let B be a 3-dimentional manifold representing the ensemble of particles making up the elastic body. The key object in relativistic elasticity is a map such that if x = (x µ ) and X = (X M ) are, respectively, coordinates on the spacetime and body manifold we then have φ M (x µ ) = X M .
As the manifolds M and B have, respectively, dimension 4 and 3, the map φ is non-injective (one-to-one). In the following it will be assumed that the inverse image φ −1 (X) of a point on B with coordinates X = (X M ) is a timelike curve on M. We denote the tangent vector to the curve γ : R → M with γ ≡ φ −1 (X) by u. If we assume γ to be parametrised by its proper time, then we have that g(u, u) = −1.
The curve γ describes the worldline of the particle of the point on B with coordinates X.
The map φ represents the configuration of the elastic body. This means that φ associates to each spacetime event a material particle. In other words, φ relates the physical state of a material body with the notion of an event in spacetime. The deformation of the elastic body is represented by the deformation gradient, defined by in terms of the coordinates at M and B by For a fixed value of the coordinate indices on the body manifold, the components φ M µ give rise to a covector field φ M a on M satisfying the condition φ M a u a = 0.
We introduce the strain of the material by applying the push-forward to the inverse metric tensor g µν on M. Its the components are given by The body manifold is equipped with a volume form V ABC which allows us to construct a scalar function n interpreted as the particle density number of the material via the relation This interpretation of n is found reasonable by the observation that the equation for particle conservation hold -that is, one has that ∇ µ (nu µ ) = 0.
In order to formulate a frame version of the energy momentum tensor of a relativistic elastic material, we begin by consider a frame {E A } on B with associated coframe {Ω B }. As we have not introduced a metric on B, we do not assume any orthonormality condition on the frame and coframe.
The map φ, defines the pullback φ * which can be used to pull-back the coframe to M. More precisely, one has that Λ B ≡ φ * Ω B , Λ B a = Λ B , e a . As the map φ is surjective and has maximal rank, the set of covectors {Λ B } is linearly independent. The fields {Λ B a } will be used, in the sequel, to describe the configuration of the material body. The coefficients Λ A a are orthogonal with respect to u a -that is Λ A a u a = 0.
We denote the determinant of the frame field as e. It is related to the determinant of the metric tensor by e = √ −g. Furthermore, we have In the above δ is understood as an infinitesimal variation. More precicely, for a function f we have, δf ≡ ∂f ∂x µ δx µ . Equation (73a) can be obtained by using Jacobi's formula for a square matrix given by,
In terms of the above fields we construct a Lagrangian on the form L = L Λ A b , e a µ . The action thus reads where we have defined the Lagrangian density The variation of the action yields δS = ∂e ∂e a µ δe a µ L + e ∂L ∂e a µ δe a µ + e ∂L where we have made use of equations (73a) and (73b) and defined T a µ ≡ ω a µ L + ∂L ∂e a µ + ∂L ∂Λ A a Λ A µ .
By multiplying with e c µ η ac , and applying the chain rule to the second term, we obtain Assuming that the Lagrangian may be written on the form (see [26] for details) Substituting the above expressions back into equation (74) we find an expression for the components of the energy-momentum tensor of the form where, in the following, Π ab will be known as the components of the Cauchy stress tensor and is given by We further make the reasonable assumption that where h ab as usual denotes the frame components of the projector metric. To show that this is reasonable, we note the following: the equation holds identically both under multiplication of u a and η ca Λ C a -in the latter case, one has to invoke the definition of h AB on the right hand side of the equation to show this. Secondly, on a spatial hypersurface S ∈ M the map φ is a diffeomorphism which implies that the object h ab defined on S is physically equivalent to the corresponding object defined on B via φ. Using this assumption in (76) we obtain the desired form of the energy momentum tensor. Namely, one has that T ab = ρu a u b + Π ab , with