{'fraction_non_alphanumeric': 0.041771252663982954, 'fraction_numerical': 0.013686952403504617, 'mean_word_length': 4.468211834779231, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equilibrium configurations in the dynamics of irrotational dust matter

Irrotational dust solutions of Einstein's equations are suitable models to describe the general-relativistic aspects of the gravitational instability mechanism for the formation of cosmic structures. In this paper we study their state space by considering the local initial-value problem formulated in the covariant fluid approach. We consider a wide range of models, from homogeneous and isotropic to highly inhomogeneous irrotational dust models, showing how they constitute equilibrium configurations (invariant sets) of the dynamics. Moreover, we give the characterization of such configurations, which provides an initial-data characterization of the models under consideration.


Introduction
Current theories of the formation of cosmic structure explain structure as the outcome of the amplification of small energy density fluctuations left after the inflationary epoch. The usually proposed mechanism for amplification is gravitational instability, which is essentially the tendency of the self-gravity of a density fluctuation to induce collapse. Within the framework of Newtonian Cosmology several successful techniques have been set up to describe this dynamical process (see [1,2] and references therein). Some remarkable and illustrative results are: (i) The Jeans instability [3]. For small density fluctuations there is a critical wavelength, the Jeans wavelength λ J , such that fluctuations with a greater wavelength will grow with time producing collapse, whereas fluctuations with a smaller wavelength will be dispersed as sound waves. (ii) The Lin-Mestel-Shu instability [4]. Given a self-gravitating, uniform, non-rotating, and pressurefree distribution of matter at rest, spherically-symmetric collapse takes place through a series of spherically-symmetric states. However, the slightest departure from sphericity is systematically magnified, giving rise to an instability. In particular, if the initial distribution is composed of oblate spheroids, the system will approach an infinitely thin disk (pancake formation), whereas if it is composed of prolate spheroids, the system will approach an infinitely thin cylinder (filament formation). Other similar studies can be found, e.g., in [5,6]. (iii) The Zel'dovich approximation [7]. This perturbative approach and its improvements (see, e.g., [1,2]) allow one to follow analytically the growth of density fluctuations into the non-linear regime and until caustics form. It predicts that the generic outcome of collapse is the formation of pancakes, although other structures like filaments or clumps (point-like singularities) can also form, but this is less probable.
However, despite the fact that Newtonian Cosmology provides a wide-ranging picture of the dynamical process of the formation of cosmic structures, there are situations in which it cannot make reliable predictions. Of particular importance is the evolution of superhorizon inhomogeneities, where the causal structure of the Newtonian theory (in which the gravitational interaction is transmitted instantaneously) is not suitable (see, e.g., [2,8]). Another interesting issue for which the relativistic approach is needed is the generation of gravitational waves by non-linear perturbations [9]. Therefore, the relativistic study of the gravitational instability mechanism is well motivated. On the other hand, the relativistic theory involves approaches that, from a technical point of view, are very different from the Newtonian ones. For instance, in order to use numerical methods with a fluid approach we would need to consider an initial-value problem formulation and to build a numerical code for the evolution equations. In the Newtonian theory we would need to solve Poisson's equation at each time whereas in General Relativity we only need to solve the constraints initially and this ensures that they will be satisfied for all later times.
Most of the relativistic approaches make assumptions on the geometric structure of the spacetime, like the imposition of symmetries (e.g., spherical symmetry), that make it impossible to obtain generic results like those described above within the Newtonian framework. In this sense, the covariant fluid approach [10,11] is an adequate starting point since it deals only with physical variables. Several studies have appeared within this framework. In [12], Irrotational Dust Models (IDMs) were studied assuming only that the magnetic part of the Weyl tensor, H ab , was negligible. The dynamical analysis favoured the spindle as the generic outcome of the collapse [13], in contrast to the Newtonian prediction, the pancake. However, the study of the dynamical consistency of the constraints led to a conjecture that only some known models were included [14] [the Friedmann-Lemaître-Robertson-Walker (FLRW), Bianchi I, and Szekeres dust models], leaving this question open. In [15] it was shown that if we relax the condition H ab = 0 , assuming instead H ab to be transverse, D b H ab = 0, the whole set of constraints is, in general, unstable (a chain of infinite new constraints appears). The main conclusion we can extract is that the dynamics of the gravitational instability mechanism is non-local (see also [16]), and hence we must be careful when imposing conditions on E ab and H ab [14,15,17] (other studies in this direction can be found in [18]).
This suggests that the use of exact analytic approaches to understand the gravitational instability mechanism is quite unrealistic. Therefore, we should consider perturbative and/or numerical methods. The aim of this paper is to contribute to the development of the latter. To that end, we will consider IDMs. Thus, we assume that pressure is negligible and that the matter fluid flow is irrotational. These assumptions are analogous to those made in the original Zel'dovich approximation [7]. Here, we will study equilibrium configurations of the dynamics (IDMs constitute an infinitedimensional dynamical system), as we will call the invariant sets of the dynamics. We will see that some exact IDMs, from homogeneous and isotropic to inhomogeneous models, constitute equilibrium configurations and we will find their characterization. On one hand, this study provides a deeper knowledge of the dynamics, especially of the role of the gradients in the evolution equations, responsible for the non-locality of the dynamics, and which were neglected in most of the particular models considerer until now in the literature. On the other hand, the characterization of particular IDMs as equilibrium configurations is at the same time an initial-data characterization. Then we know the initial data whose development corresponds to those IDMs, which provides interesting examples to check numerical codes and also information about what kind of data we must not prescribe if we want to study new behaviours. Moreover, it can be used to identify attractors, repellers or asymptotic states in the dynamics. Furthermore, an initial-data characterization provides also a solution of the constraint equations, the solving of which presents in general strong difficulties, especially if we are interested in prescribing generic initial data. This is even more important in cosmology, where we do not have physically reasonable boundary conditions in order to get a well-posed elliptic problem from the constraints. In that respect, the information that our study provides can be useful for solving the constraints in other situations.
The plan of the paper is as follows: In Sec. 2, using the covariant fluid approach and tetrad methods, and following [19,20], we introduce a local Initial-Value Problem (IVP) formulation for IDMs. In Sec. 3, using this formulation, we see that a wide range of exact IDMs constitute equilibrium configurations of the dynamics. Moreover, we find the characterization of those states, which will provide an initial-data characterization of the irrotational dust models considered. We end with some comments and remarks in Sec. 4. The notation we use in this paper follows that of previous works, especially [21]. We use units in which 8πG = c = 1, round brackets enclosing indices denote symmetrization and square brackets antisymmetrization. Throughout the paper we use coordinate charts as well as tetrads in order to express tensorial components. The conventions for indices are the following: we denote spacetime coordinate indices by the lower-case Latin letters a, . . . , l = 0, 1, 2, 3, and spacetime indices with respect to an arbitrary basis {e 0 , e 1 , e 2 , e 3 } by the remaining of lower-case Latin letters m, . . . , z = 0, 1, 2, 3. When we choose a basis adapted to the dust velocity field, i.e., e 0 = u, we will denote indices with respect to a spatial basis {e 1 , e 2 , e 3 } by lower-case Greek letters α, . . . , λ = 1, 2, 3 and spatial coordinate indices by the remaining of lower-case Greek letters µ, . . . , ω = 1, 2, 3.

An initial-value problem description of irrotational dust models
Irrotational dust solutions of Einstein's equations are models suitable for describing the gravitational instability mechanism and the dynamics of the Universe in the matterdominated era. Their energy-momentum distribution is completely described by the fluid velocity u and the energy density associated with it, ρ, the energy-momentum tensor then being The fluid flow is irrotational (u [a ∇ b u c] = 0), i.e. u generates orthogonal spacelike hypersurfaces, and geodesic (u b ∇ b u a = 0) by virtue of the vanishing of the pressure and the energy-momentum conservation equations. Hence, there is locally a function τ (x a ) such the fluid velocity is given by Then, τ is at the same time the proper time of the matter and the label of the hypersurfaces orthogonal to u, Σ(τ 1 ) ≡ {τ (x a ) = τ 1 : constant} . Choosing three independent first integrals of u, y µ (x a ) (u a ∂ a y µ = 0), {τ, y µ } is a set of comoving geodesic normal coordinates. The line element in such a coordinate system has the following form where h µν (τ, y σ ) are the non-zero components of the orthogonal projector to the fluid velocity, h ab ≡ g ab + u a u b , whose restriction to a hypersurface Σ(τ ) coincides with its first fundamental form. In order to establish an IVP for IDMs in terms of variables with a clear physical and geometrical meaning, it will be very convenient to describe them in terms of the covariant fluid approach introduced by Ehlers [10] (see [11] for more details) and using an orthonormal basis adapted to the fluid velocity, {e 0 = u, e 1 , e 2 , e 3 } [e m · e n = η mn ≡ diag(−1, 1, 1, 1)]. We can partially fix the freedom in the choice of the triad {e α } by requiring it to be parallely propagated along the fluid floẇ This choice makes the local angular velocity vanish Ω α ≡ 1 2 ε αβδ e β ·ė δ = 0 . The advantage of this choice is that we avoid having as a variable Ω α , which is not a dynamical quantity since there is no evolution equation for it.
The variables that we will use to describe IDMs can be divided into five groups (see [22,23] In this work we will also use the equivalent variables introduced by Schücking, Kundt and Behr (see [24] and references therein) The expansion Θ (≡ ∇ a u a ) and the shear tensor of the fluid worldlines. The shear is a symmetric and trace-free spatial [25] tensor It is worth noting that the restriction of the quantity to a hypersurface Σ(τ ) coincides with its second fundamental form. In this work we will consider as variables the expansion, Θ and the five independent components of the shear in a triad, σ αβ , or their following representation: (iv) Matter variables. Apart from the fluid velocity, the only matter variable in our case is the energy density ρ (= T 00 ). The Ricci tensor is then determined through Einstein's equations where Λ denotes the cosmological constant. (v) Weyl tensor variables. The Weyl tensor C abcd describes the spacetime curvature not determined locally by matter. Its ten independent components can be divided into two spatial, symmetric and trace-free tensors where η abcd is the spacetime volume 4-form. Using the analogy with electromagnetism, they are called the gravito-electric and -magnetic fields respectively. Whereas the gravito-electric field produces tidal forces having a Newtonian analogue (the tracefree part of the Hessian of the Newtonian potential), the gravito-magnetic field has no Newtonian analogue. The five independent components of the gravito-electric and -magnetic fields can be represented by the quantities (E + , , which are defined as in (4). The equations governing the behaviour of these quantities come from the Ricci identities applied to u, the second Bianchi identities, Einstein's equations (5) and the Gauss equations applied to the foliation {Σ(τ )} (see [11] for details). In the covariant fluid approach the covariant derivative ∇ a is decomposed into a time derivative along the fluid velocity,Ȧ a··· b··· ≡ u c ∇ c A a··· b··· , and a spatial covariant derivative tangent to the hypersurfaces Σ(τ ), where A a··· b··· is an arbitrary tensor field. It will be useful to introduce the spatial divergence and curl of an arbitrary 2-index tensor A ab [26] (see [21] for details on this notation) f ) is the volume 3-form of the hypersurfaces Σ(τ ). The projection of these definitions onto a triad {e α } gives ‡ where ∂ α ≡ e α µ ∂ y µ . Finally, we define the commutator of two spatial symmetric tensors, A ab and B ab , as The set of dynamical equations can be divided into two groups: (i) Evolution equations, which give the rate of change of our variables along the fluid world-lines. (ii) Constraint equations, which are relations only between spatial derivatives of our variables. The form of these equations, written in an orthonormal basis adapted to the fluid velocity, is the following

Constraint equations
‡ Angled brackets on indices denote the spatially projected, symmetric and tracefree part: In these equations, 3 R and 3 S αβ are the scalar curvature and the trace-free part of the Ricci tensor of the hypersurfaces Σ(τ ) respectively. They must be understood as given in terms of γ δ αβ and their derivatives, through the expression are the Ricci rotation coefficients associated with the triad {e α }, related to γ α βδ by δ αǫ Γ ǫ βλ = δ αǫ γ ǫ λβ + δ βǫ γ ǫ αλ + δ λǫ γ ǫ αβ . Now, we will study the evolution equations (6)(7)(8)(9)(10)(11)(12) to show that they form a firstorder symmetric hyperbolic system of partial differential equations (general accounts of this subject can be found in [27], and specializations in general relativity related with this work in [28,29]). Following [20] we can write the evolution equations (6)(7)(8)(9)(10)(11)(12) in matrix form as where we have grouped the variables (38 in all) into the vector U where the superscript T stands for the transpose of a row vector, and where U connection = (a 1 , a 2 , a 3 , n 11 , n 12 , n 13 , n 22 , n 23 , n 33 ) T , The right-hand side of (20) is the principal part of the system of equations, whose form is determined by the matrices M a , which in our case are (see also [20]) where Id N and 0 N denote the N × N identity and zero matrices respectively, and where the matrices W α are given by As we can see, A, B, and C are antisymmetric, and hence the matrices W α are symmetric. Moreover, N is a vector which depends analytically on the variables U . This formulation shows that the evolution system (20) satisfies the conditions to be a first-order symmetric hyperbolic system (FOSHS) of partial differential equations (see [19,20]). Indeed, the matrices M a are symmetric, and there is a 1-form t such that t a M a is a positive definite matrix; this is immediate upon taking t = −u (for more details see [27]). Moreover, if every 1-form t satisfying the last property is pastdirected, the system is said to be causal. We can check that (20) is causal (see also [20]). Furthermore, as we have said before, FOSHSs admit a well-posed IVP [27], which implies that given smooth initial data o U ≡ U | Σ(τo) , there exists a unique smooth solution * U . Now let us consider the constraints. First of all, we have to point out that not all of them are independent. We find the following relationships between them The last one was already given in [21,30]. This means that if constraints C 1 α , C 2 αβ , C 5 , and C 6 αβ are satisfied, constraints C 3 α and C 4 α will be satisfied automatically. Another important point concerning the constraints is their consistency with evolution. In [21,30], consistency was proved for analytical initial data. In [29], a general treatment for perfect fluids is given and consistency has also been shown for smooth initial data, but the constraints have not been split according to the variables used. In our treatment, we can construct a closed system of evolution equations for V = (U , C ≡ (C A )). It is worth noting that spatial derivatives only appear in the equations for C 3 α and C 4 α , and they have essentially the same principal part as Maxwell's equations. Therefore, the evolution equations found for V form a FOSHS. Then, given smooth initial data o U satisfying the constraints (C A [ o U ] = 0), and with * U being the only solution of (20), it is clear that ( * U , 0) is a solution, and is unique.

Equilibrium configurations and their characterization
We have just seen that the dynamical equations for IDMs constitute a symmetric hyperbolic system, which means we can formulate a well-posed IVP for irrotational dust matter. As a dynamical system, one important question is the study of equilibrium configurations in the state space, or in a more technical terminology (see, e.g., [31]), the study of invariant sets in the evolution. Here we will look for such configurations, associated with particular IDMs, finding the relations that characterize them, which also characterize the associated IDMs in terms of initial data.
Let us start by introducing some ideas, concepts and the procedure that we will use systematically. First, the state of the system will be described by the dynamical variables contained in the vector U [Eq. (21)], which we will call the state vector. Then, we will say that the state of the system at a time τ 1 is determined by U on Σ(τ 1 ), i.e. U | Σ(τ 1 ) . On the other hand, by configuration we will mean sets of states characterized by relations which do not involve explicitly time or time derivatives. That is, relations of the form For instance, the constraints [Equations (13)(14)(15)(16)(17)(18)(19)] have this form. Then a particular configuration whose characteristic relations (22) are preserved by the evolution will be called an equilibrium configuration (an invariant set). For these, R A | Σ(τ 1 ) = 0 implies that R A = 0 for τ > τ 1 . The constraints are again a good example. A consistent set of relations of the form (22), and defined in an open domain of the spacetime, determines a particular IDM in this domain. Therefore, this set of relations constitutes both a spacetime and an initial-data characterization of that IDM. An interesting example is the configuration defined by the vanishing of the gravitoelectric and -magnetic fields: As is well-known, the vanishing of these two tensors in an open domain of the spacetime implies that the metric there is a FLRW metric. However, assuming that these quantities vanish only on an open domain of a hypersurface Σ(τ 1 ), E ab | Σ(τ 1 ) = H ab | Σ(τ 1 ) = 0, this does not imply the vanishing of these quantities at a different time. This is a consequence of the evolution equations (11,12): from (11) we see that when ρ = 0 the shear acts a source generating gravito-electric field E ab , and from (12) we see that curl(E) ab acts as a source for the gravito-magnetic field H ab . Therefore, we conclude that E ab = H ab = 0, which is a spacetime characterization of the FLRW models, is not an equilibrium configuration. This reflects the obvious fact that a spacetime characterization is not in general an initial-data characterization. On the other hand, the vacuum configuration, defined by is an equilibrium configuration as we can see from the evolution equation (10).
To study whether a particular IDM determines an equilibrium configuration we will look for a set of relations (22), R ≡ (R A ), satisfied by that IDM, and such that they are preserved by the evolution. To that end, we will consider the relations R as new variables and will study their evolution. For all the particular IDMs we are going to consider in this paper, we will look for evolution equations having the following forṁ and with the property As is clear, this property implies that R = 0 is a solution. Then, if we can show that this is the only possible solution, it follows that the relations in R determine an equilibrium configuration. One way to tackle this question is to consider the variables V ≡ (U , R) and their evolution equations (20,23). In all the cases we will consider here, the coefficients appearing in these equations are analytic functions of the variables V . Therefore, the system (20,23) for V will satisfy the conditions of the Cauchy-Kowaleski theorem, which does not assume hyperbolicity at all (a statement of this theorem can be found, e.g., in [27]), and hence, the uniqueness of the solution R = 0 is ensured for analytical initial data (assuming, of course, this initial data o U is such that o R = 0).
In the cases where the complete system (20,23) is symmetric or strongly [32] hyperbolic (see, e.g., [27]) we will get a characterization also valid for smooth initial data. The constraints constitute again an example since we obtain a FOSHS for them.
In what follows we will apply this procedure to the following IDMs: In Sec. 3.1, to the FLRW models. In Sec. 3.2 to the spatially homogeneous IDMs and the Bianchi I subcase. In Sec. 3.3 to the Szekeres models. Finally, in Sec. 3.4, to the IDMs in which the hypersurfaces Σ(τ ) are flat. In all these cases we will give the dynamical system for V that shows explicitly that R = 0 is a solution of the evolution equations for the relations (22) which determine the equilibrium configuration. However, in some cases this system for the variables V will not provide a well-posed IVP for smooth initial data. In these cases we can use the fact that the relations R are defined in terms of U , and hence we can get several different systems of evolution equations for R. Moreover, we can also use the fact that equations (20) are not coupled with the equations for R, and that given smooth initial data o U , there exists a unique smooth solution * U (τ, y α ) such that * U | Σ(τo) = o U . Then given a system of equations for R , we can substitute U for * U in them so that we get a linear system of partial differential equationṡ where all the coefficients are smooth. If we can obtain a system like this such that it admits the formulation of a well-posed IVP for smooth initial data, and if we can show that R = 0 is a solution, then we would have shown that it is the solution. Therefore, our characterization would be also valid for smooth initial data. We have applied this procedure successfully to some of the cases enumerated above. For the sake of brevity, we will not write down explicitly these alternative systems, especially since in most cases they can be obtained by making a few changes in the systems (23) for R . Finally, it is important to remark that the initial-data o U in each case must satisfy, apart from the relations (22), the constraints. Taking this into account, the characterizations of the different equilibrium points that we are going to give will consist of the minimum conditions that, added to the constraints, imply the relations defining such configurations.

FLRW dust models
The FLRW models [33] are the standard cosmological models. They are motivated by the so-called Cosmological Principle in the sense that they are homogeneous and isotropic cosmological models (they have a six-dimensional group of motions). In the case of dust, the line element can be expressed in terms of elementary functions as follows (see, e.g., [34]): E is an arbitrary constant such that sign(E) = −sign(k), and the hypersurfaces Σ(τ ) have constant curvature ( 3 R = 6kR −2 ). There are several covariant spacetime characterizations of the FLRW models, and in particular of the dust subcase (see, e.g., [35]). Here, we will be concerned with the following two characterizations: (i) The kinematical characterization. In terms of the kinematical quantities a dust FLRW model is characterized simply by the vanishing of the shear σ ab = 0 .
(ii) The Weyl tensor characterization. The FLRW dust models can be characterized by the vanishing of the Weyl tensor, or equivalently, by the vanishing of the gravitoelectric and -magnetic fields From these spacetime characterizations we can see that the FLRW dust models constitute an equilibrium configuration in the dynamics of the IDMs, finding at the same time an initial data characterization for these models. To that end, we need to consider the following relations R 1 ab ≡ σ ab = 0 , R 2 ab ≡ E ab = 0 , R 3 ab ≡ H ab = 0 . Their evolution follows directly from (6-12) As is clear, R = 0 is a solution, and since the equations for (U , R) constitute a FOSHS, it is the only one. To sum up, we can say that if the initial state belongs to the configuration characterized by the vanishing of the shear and the gravito-electric field, the system will remain in this configuration. Moreover, the resulting spacetime belongs to the FLRW class. In this statement we have not taken into account the vanishing of the gravito-magnetic field because, as we have said before, we assume that the constraints hold. Apart from the vanishing of H ab , the constraints imply that is, Θ and ρ must be constant on the hypersurfaces Σ(τ ).
We can find another initial-data characterization of the FLRW models, equivalent to that given above, in terms of geometrical quantities only. It is based on the equations (3) and (18,19), and the relations that characterize it are ab , R 2 ab and the constraints imply R 3 ab . The remaining conditions imposed by the constraints, which are equivalent to (24), are All these relations together, tell us that the initial hypersurface Σ(τ o ) must be of constant curvature and the second fundamental form (which describes how it is immersed in the spacetime manifold) must be proportional to the first fundamental form, with the proportionality factor the trace K of K ab , which is a constant.

Spatially homogeneous IDMs
This well-known class of spacetimes is characterized by the existence of a threedimensional group of motions (G 3 ) simply transitive on the hypersurfaces orthogonal to the fluid velocity, {Σ(τ )}. These models were studied systematically in [24], where a classification was given. Following this work we can see that a spacetime characterization of these models is the following: "a perfect-fluid cosmological model contains a threedimensional group of motions simply transitive on three-surfaces orthogonal to the fluid velocity u if and only if there exits an orthonormal basis {e 0 = u, e α } such that the commutators satisfy γ n pq = γ n pq (τ ) ⇐⇒ ∂ α γ n pq = 0 ." The triad {e α } spans the surfaces of transitivity of the group.
Starting from the above characterization we will prove the following statement: If the initial state is such that the quantities o γ α βδ , o Θ, and o σ αβ are constant, the system will stay in this configuration. The cosmological models constructed from the development of initial data satisfying these conditions will correspond to spatially homogeneous (Bianchi) dust models. To show this, let us consider the following relations Whereas the relations in (25) will implement the characterization, the relations in (26) are auxiliary conditions needed for the proof, but they can be derived from (25) and the constraints (15)(16)(17)(18)(19). Using the evolution and constraint equations and the commutators of spatial derivatives (2) we obtain the following system of evolution equations for Ṙ As we can see, only the two last equations contain spatial derivatives, and their structure is essentially the same as that of equations (11,12). Then the extended system for (U , R) is a FOSHS, with R = 0 the only solution. This shows that Bianchi models constitute an equilibrium configuration characterized by (25). This result provides an initial-data characterization of the dust Bianchi models within the class of IDMs, but it says nothing about the particular classes of Bianchi models (see [24] for a classification). We can study which kind of equilibrium configurations particular Bianchi models may constitute by refining the characterization given in expressions (25,26). To that end, we will have to add more conditions depending on the class of Bianchi models under consideration. As an example, let us consider the Bianchi I dust models [36], in which the G 3 group is Abelian and the line-element is given by where o τ , p α are constants such that 3 α=1 p α = 3 α=1 p 2 α = 1. These models can be characterized by adding the following conditions: αβ is a consequence of R 3 αβδ [Eq. (25)], R 7 αβ , R 8 αβ and the constraint (15). We can obtain the following evolution equations for these quantities: |λ|β>δ . The whole system of evolution equations for the quantities U and R I (I = 1, . . . , 9) still constitutes a FOSHS, and R I = 0 is a solution, the only possible one. Therefore, we can say that conditions (25) and (R 7 αβ , R 8 αβ ) determine an equilibrium configuration in the dynamics of the IDMs, which corresponds to the Bianchi I dust models. In a similar way and using the information given in [24] we can study which equilibrium configurations come from the other classes of Bianchi models. The characterizations given above do not have a covariant translation. However, in some cases it is possible to get a covariant characterization, an example being Bianchi I IDMs. As in previous cases we start from covariant spacetime characterizations. For Bianchi I models two such characterizations follow from [37]. The first one is determined by the vanishing of the spatial covariant derivative of the shear D a σ bc = 0 .
The second one says that an IDM belongs to the class of the Bianchi I dust models if and only if the hypersurfaces Σ(τ ) are flat 3 R ab = 0 , and E ab (or equivalently σ ab ) is non-degenerate (algebraically general). Taking this into account, a covariant initial-data characterization of the Bianchi I dust models is given by The evolution equations for these relations form a closed system: where K ab must be understood as given in terms of the expansion and shear (3). Again, R = 0 is a solution, but now the system of evolution equations we have for (U , R) is not a FOSHS. However, using the procedure described above we can show that R = 0 is the only possible solution, showing that (28) constitutes a stationary configuration associated with the Bianchi I models (27).

Szekeres dust models
The Szekeres dust cosmological models [38] are inhomogeneous models. As was shown in [39], they do not have in general any Killing vector field. From this point of view they can be considered to be generic cosmological models, but from other points of view, as for instance the algebraic classification of spacetimes, they are not so generic since they are Petrov type D models, a special algebraic case. The line element can be written in the following diagonal form ds 2 = −dτ 2 + Σ −2 U(τ, y 1 )(dy 1 ) 2 + A 2 (τ, y α ) (dy 2 ) 2 + k 2 (y 2 , y 3 )(dy 3 where Σ coincides with the repeated shear eigenvalue σ 2 = σ 3 . From their spacetime characterization (see [35] and references therein) we can see that the evolution equations (6-12) become local equations [12], that is, the spatial derivatives disappear (there are no fields propagating in these spacetimes). Hence, the evolution equations are ordinary differential equations.
We can obtain an initial-data characterization of these models starting from the spacetime characterization (30). Taking into account the form of the third condition in (30), the use of an orthonormal basis adapted to the fluid velocity simplifies considerably the relations we will use to get the characterization. These relations are R 1 −,1,2,3 ≡ (σ − , σ 1 , σ 2 , σ 3 ) = 0 , R 2 11,22,33,23 ≡ (n 11 , n 22 , n 33 , n 23 ) = 0 , (31) where we have introduced the following definitions M ab ≡ curl(E) ab , and N ab ≡ curl(H) ab , which become new constraints For these quantities we did not give evolution equations, but we can find them using the evolution equations (6)(7)(8)(9)(10)(11)(12). The result can be expressed as follows: Moreover, the quantities CC ab and CC ′ ab represent combinations of the constraints, which vanish when the constraints hold. Here, we have used the well-known fact that adding combinations of the constraints to the evolution equations does not change the dynamics, but they can be useful for technical reasons and can modify the type of differential equations that we obtain. We will use this fact to find the initial-data characterization of the Szekeres models. Now we are ready to compute the evolution equations for the quantities R I . They can be expressed in the following forṁ The main point is to show that R = 0 is a solution for this system of equations. As we can see, this question depends only on the form of the different terms that make up the tensor F 1 αβ [see Eq. (37)]. The explicit expressions for its components are given in Appendix A, where the combination of constraints CC ab has been chosen in the following way: CC + = CC − = CC 1 = 0, and [see Eqs. (15) and (36)] With this choice we have shown in Appendix A that F 1 αβ = 0 when R I = 0. In the case of the evolution equation for R 6 αβ we have not added any constraint, i.e., we have chosen CC ′ αβ = 0. To sum up, we can find a system such that R I = 0 is a solution. The whole system of evolution equations for (U , R) is not a FOSHS system, but following the discussion above in this section we can modify the system for R to show that R I = 0 is the only possible solution. The conclusion is given in the following statement: If the initial state of the matter fluid belongs to the configuration defined by the following covariant conditions then it will belong to the configuration for all later time. It is worth noting that conditions (31)(32)(33)(34)(35) have a covariant translation. Apart from defining an equilibrium configuration, conditions (45) provide an initial-data characterization of the Szekeres models (the Szafron models [40] with constant pressure in the case of a non-zero cosmological constant). The remaining conditions are auxiliary, since they follow from (45) and the constraints, or they can be obtained just by fixing the remaining freedom in the choice of the triad {e α } .
On the other hand, we can have a covariant characterization equivalent to (45) expressed only in terms of geometrical quantities. It is given by the following relations where C ab is proportional to the conformally-invariant Cotton-York curvature tensor [41] associated with the hypersurfaces Σ(τ ),β ab . The proportionality factor is given bỹ One can check that C ab can be expressed in terms of the trace-free part of the Ricci tensor of the hypersurfaces Σ(τ ) as follows C ab has essentially the same properties as the Cotton-York tensor: i) It is spatial, symmetric and trace-free. ii) It is transverse (divergence-free): D a C ab = 0. iii) It vanishes if and only if the hypersurfaces Σ(τ ) are locally conformally-flat. The vanishing of this tensor has been sometimes related with the absence of gravitational waves [42]. For the Szekeres models this interpretation seems to be correct since C ab vanishes, and as we have said before, there is an absence of propagation signals in these models (the evolution equations are completely local). The first two conditions in (46) are clearly a direct translation of the first two conditions in (45). Therefore, to show this equivalence we only have to show that the third conditions are equivalent. To that end, we can use the constraint (19) to write C ab in the following form Then the problem reduces to checking that the terms in T ab vanish when the relations (31)(32)(33)(34)(35) hold. From (47), T ab contains only the shear and its first derivatives and hence we only have to check that its components do not contain any term of the type σ + ∂ α σ + or Γ σ 2 + , where Γ = a α , n 12 , or n 13 . Since the terms contained in T ab have the same structure as those in F 1 ab , we can use the calculations done for F 1 ab (see Appendix A). For the components T + , T − , and T 1 , it is a matter of algebra to check that the terms of the type mentioned above do not appear, and for the components T 2 and T 3 we need to use the components 2 and 3 of the constraint C 2 ab [Eq. (15)] to see that they cancel. These calculations finish the proof that the conditions (46) provide an initial-data characterization of the Szekeres models equivalent to that given in (45).

IDMs with flat hypersurfaces
The spacetime characterization of these models is obvious; they are defined by the vanishing of the Ricci tensor of the hypersurfaces Σ(τ ), i.e., From a geometrical point of view these models can be considered as inhomogeneous generalizations of the standard Einstein-de Sitter cosmological model [43]. Their relevance is enhanced by recent observations [44], which suggest that the spatial geometry of the local universe is close to flat. For the case of a vanishing cosmological constant (Λ = 0), the IDMs satisfying (48) were determined in [37]. There are two classes of solutions, both having a vanishing gravito-magnetic field, H ab = 0. The first one consists of the Bianchi I dust models (27), which are Petrov type I models. The second class, which has a degenerate gravito-electric tensor (Petrov type D), is the subclass of the Szekeres models [38] satisfying (48), which were given in [45]. Since we have already studied the first class, we will focus here on the second class, which is composed of two families of solutions. The line element of the first family is given by [(y α ) = (x, y, z)] where A, B and C are arbitrary functions of x. The line element of the second family is where U = V W , W = a y 2 + z 2 + 2by + 2cz + d Now we will show that the IDMs (49,50) constitute an equilibrium configuration in the dynamics of IDMs, and we will also find its characterization. To that end, we need to consider the following relations: (31), (33) and the alternative to (32), The evolution equations for them are given by the equations (39)(40)(41)(42), which do not change, and by the equationṡ where T αβ is defined by Eq. (47). Therefore, we have a closed system for which R I = 0 is a solution, which can be shown to be unique. Then we can say that if initially the spatial geometry is flat, the gravito-magnetic field vanishes and the shear is degenerate, the system will stay in this configuration. Again, we have been able to extract a covariant initial-data characterization from the relations that define the equilibrium configuration.

Remarks and discussion
In this paper we have considered the general dynamics of IDMs from an Initial-Value Problem perspective. The equations governing these models constitute an infinitedimensional dynamical system due to the presence of spatial gradients. We have seen that a variety of well-known models, from the homogeneous and isotropic FLRW spacetimes to the inhomogeneous Szekeres cosmological models, constitute equilibrium configurations of the dynamics (invariant sets). It remains to find the characterizations of some important IDMs, as the remaining Bianchi classes of IDMs and the Lemaître-Tolman-Bondi dust models (the initial-value problem approach made in [46] could lead to such characterization). Moreover, we have found the relations characterizing these configurations, which in most cases can be expressed in a covariant way. At the same time, these relations provide an initial-data characterization of the IDMs considered. The information obtained together with the study of the evolution of the different physical quantities clarifies what the dynamical content of these IDMs is, as well as the role of the different variables and the spatial gradients. This work must be considered as a first step towards the implementation of a more general programme to understand the full dynamics of the gravitational instability mechanism. In that respect, two interesting issues to be elucidated are the question of how a triaxial configuration will evolve in collapse (among the models treated in this paper, only homogeneous models have an algebraically general shear and gravito-electric field), and the role of non-local effects due to the appearance of the curls of the gravitoelectric and -magnetic fields. Such a programme would entail a method to prescribe generic initial data. This means finding solutions of the constraints for any possible initial distribution of matter. The programme would also include a numerical code to evolve this initial data, i.e. to construct its development, and capable of extracting the physical information. The development of these two points, which does not seem an easy task, could be complemented by the information coming from studies based on perturbative methods (perturbations of a background or iterative approaches such as the long wave-length approximation scheme).
The results of this paper provide useful information for the development of such a programme. First, they can help us in the prescription of initial data: we know the role of the variables and what kind of initial data we have to avoid in order to study new behaviours. Moreover, since the characterizations given are solutions to the constraints, this information can be helpful in solving them for other different cases. On the other hand, we can use the characterizations to check numerical codes and to identify attractors, repellers and asymptotic states of the evolution.