A dynamical system approach to inhomogeneous dust solutions

A dynamical system approac h to inhomogeneous dust solutions. Rob erto A. Sussman ‡ ∗ Instituto de F ´ ısic a, Universidad de Guanajuato, L oma del Bosque 103, L e on, Guanajuato, 37150, M ´ exic o. ‡ On sabb atic al le ave fr om Instituto de Ciencias Nuclear es, Universidad Nacional A ut´ onoma de M ´ exico (ICN-UNAM), A. P. 70–543, 04510 M´ exic o D. F., M´ exic o. W e examine numerically and qualitatively the Lema ˆ ıtre–T olman–Bondi (L TB) inhomogeneous dust solutions as a 3–dimensional dynamical system c haracterized by six critical p oints. One of the co ordinates of the phase space is an av erage density parameter, h Ω i , which b ehav es as the ordinary Ω in F riedman-Lema ˆ ıtre–Rob ertson–W alker (FLR W) dust spacetimes. The other t wo coordinates, a shear parameter and a density con trast function, con vey the effects of inhomogeneit y . As long as shell crossing singularities are absen t, this phase space is bounded or it can be trivially compactified. This space con tains sev eral in v ariant subspaces whic h define relev ant particular cases, suc h as: “parab olic” ev olution, FLR W dust and the Sc hw arzschild–Krusk al v acuum limit. W e examine in detail the phase space ev olution of sev eral dust configurations: a lo w densit y v oid formation scenario, high densit y re–collapsing univ erses with open, closed and w ormhole topologies, a structure formation scenario with a black hole surrounded by an expanding background, and the Sch warzsc hild–Krusk al v acuum case. Solution curv es (except regular centers) start expanding from a past attractor (source) in the plane h Ω i = 1, asso ciated with self similar regime at an initial singularity . Depending on the initial conditions and sp ecific configurations, the curves approach several saddle p oints as they evolv e b et ween this past attractor and other tw o p ossible future attractors: perp etually expanding curves terminate at a line of sinks at h Ω i = 0, while collapsing curves reach maximal expansion as h Ω i div erges and end up in sink that coincides with the past attractor and is also asso ciated with self similar b ehavior. P ACS num b ers: 12.60.Jv, 14.80.Ly , 95.30.Cq, 95.30.Tg, 95.35.+d, 98.35.Gi I. INTR ODUCTION Inhomogeneous dust solutions with spherical symme- try are among the oldest, simplest and most useful ex- act solutions of Einstein’s equations. These solutions w ere initially deriv ed indep enden tly by Lema ˆ ıtre (1933) and T olman (1934) and then re-derived b y Bondi (1947), hence they are known as the Lema ˆ ıtre–T olman–Bondi (L TB) solutions (see [1] and references quoted therein for a comprehensive review on these solutions and their applications). In practically all applications of these solutions the standard original v ariables in which they w ere derived are used, though alternative v ariables amenable to an “initial v alue” treatment ha ve been prop osed (see [2] and references quoted therein). While not the only type of parametrization for these solutions, we find these v ari- ables particularly useful for a numerical treatment. W e prop ose in this pap er to examine these mo dels within the framework of numeric and qualitativ e tec h- niques kno wn generically as “dynamical systems” [3]. This approach requires as a first necessary step express- ing the ev olution equations as a proper system of first or- der autonomous ODE’s (ordinary differential equations). This necessary requirement w ould, apparen tly , exclude inhomogeneous L TB solutions because their evolution equations are necessarily PDE’s (partial differential equa- tions). How ever, under a “3+1” cov ariant decomposition ∗ Electronic address: sussman@nucleares.unam.mx based on the 4–velocity field [4], the ev olution equations for L TB dust solutions b ecome first order autonomous PDE’s con taining only time deriv ativ es [3]. W e pro ve in this case ho w the spacelik e gradient equations (con- strain ts) are automatically satisfied for all times once initial conditions are selected so that these constrain ts hold in an arbitrary initial h yp ersurface 3 T i of constan t time. Rigorously (see Appendix B), the solution curv es of these evolution equations are then equiv alent to a prop er dynamical system built with ODE’s, under the strong re- striction of complying with the sp ecial set of initial con- ditions that is compatible with the spacelike constraints of the PDE system. In order to describe L TB dust solutions as a dynamical system w e re–write the first order evolution equations of the 3+1 decomposition [4] in terms of initial v alue v ariables defined in reference [2]. These v ariables are in tro duced and generalized in section I I I and in sections IV and V w e show ho w they lead in a v ery natural w ay to a system of three autonomous evolution equations that is similar to those obtained with “expansion normalized” v ariables for other known space–times (FLR W, Bianc hi mo dels, Kanto wsky–Sac hs, see [3]). Initial conditions (section VI) are selected from tw o “primitiv e” initial v alue functions (densit y , scalar 3– curv ature) defined along an initial and regular h yp er- surface, 3 T i mark ed b y constan t comoving time. These functions determine the initial v alues for the phase space v ariables in suc h a w ay that the t wo spacelik e constrain ts c haracterizing the 3+1 evolution equations for the L TB solutions are satisfied at the 3 T i and at subsequent 3 T . A third initial function determines the num b er of Regular 2 Symmetry Centers (RSC) at the 3 T i (see App endix A), and thus it determines the top ological class of the 3 T , leading to h yp ersurfaces with “op en” top ology , (home- omorphic to R 3 , one RSC); “closed” top ology (homeo- morphic to S 3 , tw o RSC) “wormhole” (homeomorphic to S 2 × R or to S 2 × S 1 , zero RSC). The resulting 3–dimensional phase space is discussed in section VI I, with its critical p oints, in v ariant subsets and particular solutions given in section VI I I. This space is parametrized b y three functions: a dimensionless shear scalar, S , a densit y a v erage parameter, h Ω i , and a den- sit y contrast function ∆ ( m ) . As long as initial condi- tions are selected so that unphysical shell crossing sin- gularities do not arise [2, 5, 7], these v ariables are either b ounded or can b e trivially compactified for all configu- rations. These functions hav e a straightforw ard physical in terpretation: the av erage h Ω i b ehav es exactly as the standard Ω in a FLR W dust spacetime: its initial v al- ues h Ω i i determine the dynamical ev olution of eac h so- lution curve: re–collapse (if h Ω i i − 1 > 0) or perp etual expansion (if h Ω i i − 1 ≤ 0). Th us, h Ω i defines inv ari- an t subspaces h Ω i = 0, h Ω i = 1, 0 < h Ω i < 1 and h Ω i > 1. The densit y contrast function ∆ ( m ) pro vides a non–lo cal c haracterization of the t yp e of inhomogeneit y along the 3 T (rest frames of fundamental observ ers): a “clump” if − 1 < ∆ ( m ) ≤ 0 or a “v oid” if ∆ ( m ) ≥ 0, with ∆ ( m ) = 0 at the RSC. This v ariable also defines an in v ariant subset: ∆ ( m ) = − 1, corresp onding to v acuum Sc hw arzschild–Krusk al solutions. Another inv arian t sub- space is given by the line S = ∆ ( m ) = 0 marking the homogeneous and isotropic FLR W subcase, containing as sub cases the Einstein–de Sitter and Minko wski–Milne univ erses. W e examine in detail in section IX the phase space ev olution of several representativ e dust configurations: lo w densit y perp etually expanding (“h yp erb olic” dynam- ics) with h Ω i i − 1 ≤ 0, as w ell as a case of v oid forma- tion scenario in which an initial density clump b ecomes a void (a n umeric realization of configurations prop osed in [10]). W e also examine configurations with high den- sit y re–collapsing (“elliptic”) dynamics with op en, closed and wormhole top ologies ( h Ω i i − 1 > 0); a “structure formation” scenario: solution curv es near the RSC re– collapse into a blac k hole h Ω i i − 1 > 0 while “exter- nal” curv es perp etually expand in to a cosmic bac kground h Ω i i− 1 ≤ 0. W e also examine the v acuum sub case, whic h is the Sc hw arzschild–Krusk al spacetime giv en in comov- ing non–static co ordinates made up b y radial timelike geo desics. All solution curves (except the RSC) for all configura- tions start their evolution at a past attractor (source) at h Ω i = 1 , S = 1 / 2 , ∆ ( m ) = − 1 that represents an initial big bang lik e singularity and is associated with a self sim- ilar regime. The evolution of eac h solution curv e depends on the sign of h Ω i i − 1 and is fully con tained in the in v ari- an t subsets asso ciated with h Ω i . Curves with h Ω i i − 1 = 0 (“parab olic” evolution) remain in the plane h Ω i = 1 and terminate in a future attractor (sink) associated with the Einstein–de Sitter zero spacial curv ature FLR W universe. Curv es with h Ω i i − 1 < 0 (“hyperb olic” evolution) ev olve to wards a future attractor given b y a line of sinks at h Ω i = S = 0, while curves with h Ω i i − 1 > 0 (“elliptic” re–collapsing ev olution) evolv e into a maximal expansion state h Ω i → ∞ and then terminate at the sam e critical p oin t from which they started, though this point is no w a sink or a future attractor. The RSC’s of the configurations evolv e separately from the rest of the curv es along a line with ∆ ( m ) = S = 0, starting alwa ys from a source at h Ω i = 1 and going to- w ards a sink at h Ω i = 0 or tow ards h Ω i → ∞ , resp ec- tiv ely , for expanding and re–collapsing configurations. In all cases some of the curves pass near a saddle close to the RSC asso ciated with homogeneity ( S = ∆ ( m ) = 0), while in the structure formation scenario some curves also approach another saddle p oint that splits collapsing curv es h Ω i i > 1 from expanding ones h Ω i i ≤ 1. The solu- tion curves of the v acuum case are confined to the plane ∆ ( m ) = − 1, ev olving tow ards h Ω i = 0 or h Ω i → ∞ , de- p ending on whether the radial geodesics are bound or not. W e summarize and discuss the results from the article in section X. W e also provide three appendices dealing with important issues: RSC’s and geometric prop erties of h yp ersurfaces 3 T (App endix A), treatmen t of ev olu- tion equations given as PDE as a restricted dynamical system (App endix B) and analytic solutions in terms of the v ariables used in the article (App endix C). I I. L TB DUST MODELS IN THEIR ORIGINAL V ARIABLES. The Lema ˆ ıtre–T olman–Bondi (L TB) metric [1, 2, 5, 6] is the spherically symmetric line element ds 2 = − c 2 dt 2 + Y 0 2 1 − K dr 2 + Y 2 ( dθ 2 + sin 2 θ dφ 2 ) , (1) where Y = Y ( t, r ), Y 0 = ∂ Y /∂ r and K = K ( r ). The momen tum-energy tensor usually associated with (1) is that of a dust source: T ab = ρc 2 u a u b , (2) where ρ = ρ ( t, r ) is the rest–mass density (in gm cm − 3 ). F or the metric (1) and source (2) with a comoving 4– v elo city u a = δ a 0 , Einstein’s field equations G a b = κT a b , with κ = 8 π G/c 4 , yield: ˙ Y 2 = 2 M Y − K, (3) 2 M 0 = κρc 2 Y 2 Y 0 , (4) where M = M ( r ) has units of length and ˙ Y = ∂ Y /∂ x 0 . The only nonzero kinematic parameters are the ex- pansion scalar Θ = u a ; a and the shear tensor: σ ab = 3 ∇ ( a ; b ) u − (Θ / 3) h ab , given b y Θ = 2 ˙ Y Y + ˙ Y 0 Y 0 , (5) σ a b = diag [0 , − 2Σ , Σ , Σ] , Σ = 1 3 " ˙ Y Y − ˙ Y 0 Y 0 # . (6) Other imp ortant quantities are the “Electric” W eyl ten- sor: E ab = C abcd u c u d and the 3–dimensional Ricci ten- sor, 3 R , of the hypersurfaces orthogonal to u a (mark ed b y constant v alues of x 0 = ct ): E a b = diag [0 , − 2 E , E , E ] , E = M 0 3 Y 2 Y 0 − M Y 3 = κ 6 ρc 2 − M Y 3 . (7) 2( K Y ) 0 = 3 R Y 2 Y 0 . (8) The usual treatmen t of L TB dust solutions consists in solving analytically the evolution equation (3). The so- lutions are usually classified in terms of the sign of K , whic h determines the t ype of evolution of dust lay ers: p erp etually expanding/collapsing “parab olic” ( K = 0) and “hyperb olic” ( K < 0), and re-collapsing “elliptic” ( K > 0). These analytic solutions are given in App endix C. How ever, instead of using these solutions, w e will fol- lo w in this article a different approach based on defining new v ariables more suitable for a qualitative and n umer- ical analysis. I I I. V OLUME A VERA GES AND SCALING LA WS The form of ρ and 3 R in terms of 2 M 0 and ( K Y ) 0 in (4) and (8) suggests considering the following v olume a verages along an arbitrary hypersurface 3 T i orthogonal to u a and marked b y constant x 0 i = ct i h A i i ≡ R A i d V i R d V i , d V i = Y 2 i Y 0 i sin θ dr dθ dφ, (9) where the subindex i denotes ev aluation at t = t i and we ha ve assumed that 3 T i is fully regular in the integration range. Applying (9) to (4) and (8) w e get κc 2 3 h ρ i = κc 2 3 R ρ Y 2 Y 0 dr R Y 2 Y 0 dr = 2 M Y 3 , (10) 1 6 h 3 Ri = 1 6 R 3 R Y 2 Y 0 dr R Y 2 Y 0 dr = K Y 2 , (11) where w e hav e dropped the subindex i and ha ve tak en the lo wer bound of the integration range of r to b e deter- mined b y a suitable b oundary condition on ρ i and 3 R i , for example a RSC. See App endix A and section IX–E. The a v erages h ρ i and h 3 Ri are non–lo cal quan tities de- p ending on the form of ρ and 3 R along the integration range. W e define the following “con trast functions” com- paring these quan tities with their lo cal counterparts: ∆ ( m ) ≡ ρ − h ρ i h ρ i , ⇒ ρ = h ρ i [1 + ∆ ( m ) ] , (12) ∆ ( k ) ≡ 3 R − h 3 Ri h 3 Ri , ⇒ 3 R = h 3 Ri [1 + ∆ ( k ) ] , (13) The interpretation of these these con trast functions [2] follo ws by using the definitions (4), (10) and (12) and in tegrating by parts along the 3 T , leading to ∆ ( m ) = κc 2 6 M Z ρ 0 Y 3 dr , (14) where this integral is ev aluated from a RSC. Since M ≥ 0 and Y ≥ 0 and assuming ρ ≥ 0, the only quantit y that can change sign is ρ 0 , therefore, lo oking at the radial v ariation of a densit y profile from the symmetry cen ter along an arbitrary 3 T we hav e ∆ ( m )  ≤ 0 ⇔ ρ 0 ≤ 0 , ρ ≤ h ρ i , density clump ≥ 0 ⇔ ρ 0 ≥ 0 , ρ ≥ h ρ i , density void (15) Lik ewise, for the contrast function ∆ ( k ) , using the def- initions (8), (11) and (13) and integrating by parts along the 3 T we obtain an analogous expression to (14): ∆ ( k ) = 1 6 K Z 3 R 0 Y 3 dr , (16) Ho wev er, while a negative ρ is not ph ysically interesting, there is no physical ob jection for 3 R being p ositive or negativ e, or changing sign in the allow ed range of r . If 3 R > 0 for all the allow ed range of r , we hav e the sam e results as with ∆ ( m ) in (15), but if 3 R < 0, then curv a- ture clumps or v oids are defined b y the opposite signs [2]. W e define now the follo wing “scale factors” related to the metric functions: ` ≡ Y Y i , (17) Γ ≡ Y 0 / Y Y 0 i / Y i = 1 + ` 0 /` Y 0 i / Y i , (18) relating Y and Y 0 ev aluated at an arbitrary and a fiducial (or “initial”) hypersurface 3 T i . Since (10) and (11) are v alid at any 3 T , we obtain the scaling laws h ρ i = h ρ i i ` 3 , (19) h 3 Ri = h 3 R i i ` 2 . (20) while from (4), (8) and (17)–(18) we get ρ = ρ i ` 3 Γ , (21) 3 R = 1 ` 2 Γ  3 R i + 1 3 h 3 R i i (1 − Γ)  , (22) 4 Comparing (12)–(13) with (19)–(20) and (21)–(22) w e get the scaling laws for ∆ ( m ) and ∆ ( k ) 1 + ∆ ( m ) = 1 + ∆ ( m ) i Γ , (23) 2 3 + ∆ ( k ) = 2 / 3 + ∆ ( k ) i Γ , (24) The quantities Θ, Σ and E given in (5), (6) and (7) tak e the forms Θ = 3 ˙ ` ` + ˙ Γ Γ , (25) Σ = − ˙ Γ 3Γ , (26) E = κc 2 6 [ ρ − h ρ i ] = κc 2 6 h ρ i ∆ ( m ) , (27) the evolution equation (3) b ecomes ˙ ` 2 ` 2 = κc 2 3 h ρ i − 1 6 h 3 Ri , (28) while the L TB metric (1) takes the F riedmanian form ds 2 = − c 2 dt 2 + ` 2  Γ 2 Y 0 i 2 dr 2 1 − 1 6 h 3 R i i Y 2 i + Y 2 i  dθ 2 + sin 2 θ dφ 2   , (29) It is also helpful to express the spacial ( i.e. r adial ) gra- dien ts of h ρ i and h 3 Ri in terms of the contrast functions and Γ. With the help of (4), (8), (10), (11), (12) and (13) w e obtain the follo wing useful relations: h ρ i 0 h ρ i = 3 Y 0 Y ∆ ( m ) = 3 Y 0 i Y i ∆ ( m ) Γ , (30) h 3 Ri 0 h 3 Ri = 3 Y 0 Y ∆ ( k ) = 3 Y 0 Y ∆ ( k ) Γ , (31) Notice that these equations are strictly v alid along all h yp ersurfaces 3 T . IV. 3+1 DECOMPOSITION. F rom a co v ariant 3+1 decomposition of spacetime based on the 4–velocity field u a presen ted in [3, 4], the field and energy balance equations for a dust model char- acterized by (1) and (2) are equiv alent to the following set of first order evolution equations: ˙ Θ 3 = −  Θ 3  2 − 1 3 σ ab σ ab − κc 2 6 ρ = 0 , (32a) ˙ σ h ab i = − 2 3 Θ σ ab − σ h a c σ b i c − E ab = 0 , (32b) ˙ E h ab i = − κc 2 2 ρ σ ab − Θ E ab + 3 σ h a c E b i c , (32c) ˙ ρ = − ρ Θ , (32d) together with the constraints: ˜ ∇ b σ b a − 2 3 ˜ ∇ a Θ = 0 , (33a) ˜ ∇ b E b a − κc 2 3 ˜ ∇ a ρ = 0 , (33b) where ˙ A ab ≡ u c ∇ c A ab is the conv ectiv e deriv ativ e along the 4–velocity , ˜ ∇ a A bc ≡ h a d ∇ d A bc is the spacelik e gra- dien t, tangen t to the 3 T hypersurfaces and orthogonal to u a , and A h ab i ≡ A ( ab ) − (1 / 3) A c c h ab is the spacelike, symmetric, trace–free part of a tensor A ab . Equations (32a)–(32d) are 4 tensorial equations for 4 tensorial quantities Θ , ρ, σ a b , E a b , which together with the constraints (33a)–(33b) form a completely deter- mined system of partial differen tial equations (equiv a- len t to Einstein’s field equations). Bearing in mind that u a = δ a 0 and considering the forms of the trace–free ten- sors σ a b and E a b in (6) and (7), the system (32) reduces to the follo wing set of scalar equations: ˙ Θ 3 = −  Θ 3  2 − 2Σ 2 − κc 2 6 ρ = 0 , , (34a) ˙ Σ = − 2 3 Θ Σ + Σ 2 − E , (34b) ˙ E = −E − Σ h κ 2 ρc 2 + 3 E i , (34c) ˙ ρ = − ρ Θ , (34d) 5 and (33a)–(33b) b ecome Σ 0 + Θ 0 3 + 3 Σ Y 0 Y = 0 , (35a) E 0 + κρ 0 c 2 6 + 3 E Y 0 Y = 0 . (35b) It is straightforw ard to express the system (34) and the constrain ts (35) in therms of the new v ariables introduced in the previous sections. Considering (12), (23) and (25)– (27) and and defining H ≡ ˙ ` ` = Θ 3 + Σ , (36) the system (34) becomes the follo wing set of equiv alent ev olution equations: ˙ H = − H 2 − κ 6 h ρ i , (37a) ˙ Σ = 3Σ 2 − 2 H Σ + κ 6 h ρ i ∆ ( m ) , (37b) h ρ i ˙ = − 3 h ρ i H, (37c) ˙ ∆ ( m ) = 3[1 + ∆ ( m ) ] Σ , (37d) while the constraints (35b) and (35a) become resp ectively h ρ i 0 − 3 h ρ i ∆ ( m ) Γ Y 0 i Y i = 0 , (38a) H 0 − 3 Σ Y 0 Y = H 0 − 3 Σ Γ Y 0 i Y i = 0 . (38b) The constraint (33b) leads to (38a) which is iden tical to (30), so this constrain t is automatically satisfied at all t (all 3 T ) b y a system of ev olution equations based on the v ariables h ρ i and ∆ ( m ) . Bearing in mind (17), (18) and (26), the constrain t (33a) in its form (38b) can b e written as the in tegrability condition for ` : Γ ˙ ` ` ! 0 − ˙ Γ  Y 0 i Y i + ` 0 `  = Γ " ˙ ` ` ! 0 −  ` 0 `  ˙ # = 0 . (39) F urther, since (39) holds and (36) implies H = ˙ ` ` =  κc 2 3 h ρ i − 1 6 h 3 Ri  1 / 2 , (40) the constraint (38b) implies (with the help of (23), (24), (30), and (31)): Σ = − 1 2 1 3 κc 2 h ρ i ∆ ( m ) − 1 6 h 3 Ri ∆ ( k )  1 3 κc 2 h ρ i − 1 6 h 3 Ri  1 / 2 . (41) If the system (37) is self–consistent, then (41) must hold for every 3 T and so it must comply with (37b). In- serting (41) in to (37b) and eliminating time deriv atives of h ρ i , ∆ ( m ) , H , h 3 Ri and ∆ ( k ) with the help of (37a), (37c), (37d) and the scaling laws (23)–(24), it is straight- forw ard to see that (37b) is identically satisfied at every 3 T . V. DIMENSIONLESS V ARIABLES: A D YNAMICAL SYSTEM In order to transform system (37a)–(37d) into a dy- namical system, we need to recast H , Σ and h ρ i in terms of dimensionless “expansion normalized” v ariables. [3]. Ho wev er, instead of using the expansion parameter Θ in (25), it turns out to be easier to use H defined in (36). This leads to: H ≡ H H 0 , (42) h Ω i ≡ κc 2 h ρ i 3 H 2 , (43) S ≡ Σ H , (44) where H 0 is a characteristic length scale. Notice that we can also define a dimensionless Ω parameter for the lo cal densit y , which would b e related to that of (43) by: Ω = κc 2 ρ 3 H 2 = h Ω i [1 + ∆ ( m ) ] . (45) W e also need to mak e the “dot” deriv ative ∂ /∂ ct a di- mensionless op erator, so we define: ∂ ∂ τ ≡ 1 H 0 ∂ c∂ t . (46) Inserting (42)–(46) in to (37a)–(37d) we obtain the fol- lo wing dimensionless system: H ,τ = −  1 + 1 2 h Ω i  H 2 , (47a) S ,τ = H  S (3 S − 1) + 1 2  ∆ ( m ) + S  h Ω i  , (47b) h Ω i ,τ = H h Ω i [ h Ω i − 1 ] , (47c) ∆ ( m ) ,τ = 3 H S h 1 + ∆ ( m ) i . (47d) where ,τ = ∂ /∂ τ . The form of these equations suggests that further sim- plification follows if the three equations (47b)–(47d) can b e decoupled from (47a). Defining the following co ordi- nate transformation: τ = τ ( ξ , ¯ r ) , r = ¯ r (48) so that for an y function A = A ( τ , r ) we hav e A ( τ , r ) = A ( ξ ( τ , ¯ r ) , ¯ r ) and all partial time driv ativ es in (47) can b e expressed as  ∂ A ∂ τ  r = ∂ A ∂ ξ  ∂ ξ ∂ τ  r = H ∂ A ∂ ξ , (49) where ξ is selected so that:  ∂ ξ ∂ τ  r = H , ⇒ ξ = ln `, (50) 6 The in tro duction of the new v ariable ξ defined by (49) and (50) remov es the dependence on H of the ev olution equations (47b), (47c) and (47d): ∂ S ∂ ξ = S (3 S − 1) + 1 2  ∆ ( m ) + S  h Ω i , (51a) ∂ ∆ ( m ) ∂ ξ = 3 S h 1 + ∆ ( m ) i , (51b) ∂ h Ω i ∂ ξ = h Ω i [ h Ω i − 1] , (51c) while (47a) b ecomes: ∂ ∂ ξ (ln H ) = −  1 + 1 2 h Ω i  , (52) The system (51) is formally a dynamical system con- structed with “expansion normalized” v ariables [3]. How- ev er, it is still a system of PDE’s, even if it only con tains deriv atives with resp ect to ξ and r is basically a param- eter. As we show in App endix B, this system of PDE’s equiv alent authonomous ODE system with restricted ini- tial conditions (restricted in order to fulfill the spacelik e constrain ts). VI. INITIAL CONDITIONS Bearing in mind (4), (11), (12)–(13) and (40)–(41), in tial conditions for the system (51) can b e constructed through the follo wing steps: 1. Cho ose the top ological class of the initial h yp er- surface 3 T i in terms of the num b er of RSC (see App endix A). A conv enien t choice is: Y i = H − 1 0 f ( r ) , (53) where f ( r ) is a (at least a C 2 ) function whose zeros corresp ond to RSC’s and H 0 is the same character- istic length scale as in (42) and (46). Strictly sp eak- ing, the choice of f is a choice of radial co ordinate (see App endix A). 2. Construct dimensionless quan tities out of ρ i and 3 R i : m i ≡ κc 2 ρ i 3 H 2 0 , k i = 3 R i 6 H 2 0 , (54) 3. Obtain the orbit volume av erages and contrast functions h m i i = κc 2 3 H 2 0 h ρ i i , (55a) h k i i = 1 6 H 2 0 h 3 R i i , (55b) ∆ ( m ) i = m i h m i i − 1 , (55c) ∆ ( k ) i = k i h k i i − 1 , (55d) 4. The initial forms of the remaining v ariables are then: H i = [ h m i i − h k i i ] 1 / 2 , (56a) h Ω i i = h m i i h m i i − h k i i , (56b) S i = − 1 2 h m i i ∆ ( m ) i − h k i i ∆ ( k ) i h m i i − h k i i , (56c) Just as w e pro ceeded with (37), the fulfillment of con- strain ts (38a) and (38b) imply that the functional forms of H , h Ω i and S are the generalization of (56) for all times, that is: H = [ h m i − h k i ] 1 / 2 , (57a) h Ω i = h m i h m i − h k i , (57b) S = − 1 2 h m i ∆ ( m ) − h k i ∆ ( k ) h m i − h k i , (57c) where h m i = κc 2 h ρ i 3 H 2 0 = h m i i ` 3 , (58a) h k i = h 3 Ri 6 H 2 0 = h k i i ` 2 . (58b) Inserting (57) and (58) into (51) and eliminating deriv a- tiv es of h m i and h k i by means of (24), (26), (37a), (37c) and (37d), w e can see that (57) and (58) fully solv e (47). Th us (47) propagates the initial data giv en by (53)–(56c) with the constrain ts (38a) and (38b) satisfied for all 3 T . The standard approach to L TB dust solutions is based on solving the evolution equation (3). Therefore, it is il- lustrativ e to re–write suc h equation in terms of the initial v alue v ariables that we ha ve defined here: ` 2 ,τ = H 2 i  h Ω i i ` − ( h Ω i i − 1)  . (59) This form of (3) shows ho w the sign of h Ω i i − 1 determines the t yp e of ev olution of dust la y ers, just as the sign of K do es it in (3), leading to “parab olic” ( h Ω i i − 1 = 0), h yp erb olic” ( h Ω i i − 1 < 0) and “elliptic” ( h Ω i i − 1 > 0) ev olution like setting K = 0 , K < 0 , K > 0 in (3). See [5, 6]. Analytic solutions of (59) are given in App endix C. In the numerical examination of dust configurations we will need to solv e numerically system (37), b esides system (51). In terms of the v ariables and initial conditions (53)– (58) the system (37) takes the dimensionless form: H ,τ = −H 2 − h m i 2 , (60a) s ,τ = 3 s 2 − 2 H s + h m i 2 ∆ ( m ) , (60b) h m i ,τ = − 3 h m i H , (60c) ∆ ( m ) ,τ = 3(1 + ∆ ( m ) ) s, (60d) 7 where τ , H and h m i ha ve been defined in (42), (46), (58) and s = Σ /H 0 . This system will b e specially useful for lo oking at the collapsing stage of re–collapsing configu- rations, a feature that cannot b e studied by (51) b ecause maximal expansion takes place as H → 0, and so h Ω i di- v erges and the solution curves cannot b e extended further to study their collapsing stage. Also, (60) yields directly quan tities lik e h m i and H , while the metric functions ` and Γ in (29) follow from the v ariables in (60) as: ` =  h m i i h m i  1 / 3 , Γ = 1 + ∆ ( m ) i 1 + ∆ ( m ) , (61) and the so–called “curv ature radius” Y = √ g θθ can be easily found as Y = Y i ` . VI I. SOME QUALIT A TIVE AND ANAL YTIC RESUL TS. The solution curv es [ S ( ξ , r ) , ∆ ( m ) ( ξ , r ) , h Ω i ( ξ , r )] of (51) evolv e in a phase space which is a region of R 3 parametrized by the co ordinates [ S, ∆ ( m ) , h Ω i ]. The evo- lution v ariable is ξ , while each curve represents a funda- men tal comoving observer labelled by a constant v alue r . It is highly desirable that the phase space co ordinates remain b ounded as these curves evolv e along their maxi- mal range of ξ for all r . W e discuss in this section several useful analytic results and the relation b etw een the phase space v ariables and their initial v alues. Considering that ξ = ln ` and h Ω i i − 1 = h k i i h m i i − h k i i , (62a) h Ω i − 1 = h k i h m i − h k i = h k i i e ξ h m i i − h k i i e ξ , (62b) the functions H , h Ω i and S in (57)–(58) can b e given as h Ω i = h Ω i i h Ω i i − [ h Ω i i − 1]e ξ , (63a) H 2 = H 2 i  h Ω i i e − 3 ξ − ( h Ω i i − 1) e − 2 ξ  , (63b) S = − h Ω i i ∆ ( m ) − [ h Ω i i − 1] ∆ ( k ) e ξ 2 [ h Ω i i − [ h Ω i i − 1] e ξ ] . (63c) where we are emphasizing that these are now functions of ( ξ , r ). Notice that even if, in general, the surfaces of constan t ξ do not coincide with the 3 T (surfaces of con- stan t t or τ ), the initial hypersurface ξ = 0 do es coincide with the initial 3 T i . A. Constrain ts on h Ω i F rom(57), w e hav e the following imp ortant constraint: h Ω i − 1 = h k i H 2 = h k i i ` 2 H 2 = h Ω i i − 1 ` 2 H 2 . (64) Since h Ω i ≥ 0, this constrain t (together with (63a)) has imp ortan t qualitative consequences: it defines in v ariant subsets asso ciated with h Ω i . If initial conditions for a range of r are selected so that h Ω i i = 0 , 0 < h Ω i i < 1 , h Ω i i = 1, or h Ω i i > 1, then all solution curv es for suc h range resp ectively comply with h Ω i = 0 , 0 < h Ω i < 1 , h Ω i = 1, or h Ω i > 1 for all their maximal extensibility range of of the evolution parameter ξ . Notice from (63) that, irresp ectiv ely of the sign of h Ω i ( r ) i − 1, we ha ve h Ω i → 1 as ξ → −∞ (or, equivlen tly , as ` → 0). Thus, irrespectively of the type of dynamics in(40) given by the sign of h Ω i i − 1, all solution curves of (51) will approach the plane h Ω i = 1 near the initial big bang singularity (in the asymptotic range: ξ → −∞ or ` → 0). Since h Ω i ≥ 0, (64) implies that all solution curv es with initial conditions 0 ≤ h Ω i i ≤ 1 will b e constrained to the inv ariant subset given by the region 0 ≤ h Ω i ≤ 1 for all v alues of ξ , while (63) implies that the such curv es will ev olve from the plane h Ω i = 1 in ξ → −∞ ( ` → 0) to wards the plane h Ω i = 0 in the asymptotic range ξ → ∞ (or ` → ∞ ). Equations (63) place no restriction in the range of ξ of solution curv es when h Ω i i ≤ 1, but curv es with initial conditions given by h Ω i i > 1 the evolution parameter ξ will ha v e a bounded range of maximal extendibilit y giv en b y h Ω i → ∞ , as ξ → ln h Ω i i h Ω i i − 1 . (65) Suc h solution curves start from the plane h Ω i = 1 as ξ → −∞ ( ` → 0) and diverge for the limiting v alue of ξ = ln ` given by (65). Ho w ever, this is the same finite v alue of ξ that makes H in (63b) v anish, so the blo wing up h Ω i → ∞ of all solution curves in this case simply marks the v alues of ` asso ciated with the “maximal” ex- pansion of dust lay ers (maximal v alue of ` ) in the “turn around” b efore the collapsing stage. Since h Ω i blows up the collapsing stage cannot b e describ ed by (51) (though it can b e describ ed by (60)). B. Constrain ts on S and ∆ ( m ) In order to appreciate the relation b etw een the co ordi- nates of the phase space, [ S, ∆ ( m ) , h Ω i ], it is very useful to express S given by (57c) as: S = − 1 2 n h Ω i ∆ ( m ) − [ h Ω i − 1] ∆ ( k ) o , (66) where, from the scaling laws (23) and (24), we hav e ∆ ( m ) = 1 + ∆ ( m ) i Γ − 1 , (67a) ∆ ( k ) = 2 / 3 + ∆ ( k ) i Γ − 2 3 . (67b) 8 Th us, for finite ∆ ( m ) i and ∆ ( k ) i , it is eviden t that a suf- ficien t condition for ∆ ( m ) and ∆ ( k ) to remain bounded for all ξ is that initial conditions are selected so that a shell crossing singularity do es not emerge (see App endix A and references [2, 5, 6, 7]). That is, we must hav e for all solution curv es mark ed by ( ξ , r ) Γ( ξ , r ) > 0 , no shell crossing singularity , (68) where Γ has b een defined in (18). Notice, from (63c), (66) and (67), that even if (68) holds, S (lik e h Ω i ) div erges as H → 0 (maximal expansion) for curv es with h Ω i i > 1, though it remains b ounded if h Ω i i ≤ 1. In order to o vercome the blowing up of S and h Ω i in the cases when h Ω i i > 1 we hav e used the system (60) to examine the collapsing stage of the solution curves. It is imp ortan t to remark that for initial conditions for which b oth (68) and 0 < h Ω i i ≤ 1 hold, the three co ordinates of phase space, [ S, ∆ ( m ) , h Ω i ], of all solution curv es remain bounded and restricted to a finite region of R 3 . Initial conditions that guarantee the fulfillment of (68) w ere derived in terms of original v ariables in reference [7] (see also [5, 6]), in terms of initial v alue v ariables dis- cussed here in [2] (see App endix A). VI I I. CRITICAL POINTS AND P AR TICULAR CASES In the general case in which none of the v ariables [ S, ∆ ( m ) , h Ω i ] is restricted to tak e any special constant v alue associated with particular cases, the coordinates and nature of critical p oints of (51) are: • h Ω i = 1: C 1 S = 1 / 2 , ∆ ( m ) = − 1 , source , (69a) C 2 S = − 1 / 3 , ∆ ( m ) = − 1 , saddle , (69b) C 3 S = 0 , ∆ ( m ) = 0 , saddle , (69c) • h Ω i = 0: C 4 S = 1 / 3 , ∆ ( m ) = − 1 , saddle , (70a) C 5 S = 0 , ∆ ( m ) arbitrary , sink , (70b) C 6 S = 0 , ∆ ( m ) = 0 , sink , (70c) The nature of some of the critical points of (51) changes when [ S, ∆ ( m ) , h Ω i ] take certain specific particular v alues that corresp ond to v arious space–times that are particu- lar cases of dust L TB solutions. W e examine these par- ticular cases and their critical p oints in the remaining of this section. The phase space and the lo cation of critical p oin ts and lo cation of all particular cases is depicted in figure 1, FIG. 1: Phase space The figure shows all critical p oints and the in v arian t subspaces given by the surfaces h Ω i = 1 and ∆ ( m ) = − 1 (Sc hw arzschild v acuum case). Particular cases, suc h as FLR W, parab olic evolution and the v acuum case in Lema ˆ ıtre co ordinates, are also inv ariant subspaces shown as lines going from a source to a sink. The FLR W sub–case con tains the Einstein–de Sitter (EdS) univ erse h Ω i = 1 (the critical p oint C 3 ) and the Minko wski–Milne space–time (M) h Ω i = 0 (the critical p oint C 6 ). A. P arab olic (or “marginally b ound”) evolution F or initial conditions h Ω i i = 1 equations (63) imply h Ω i = 1 for all ξ . This is an inv ariant set associated with the so–called “parab olic” solutions of (59) and cor- resp onding to h k i = h 3 Ri = K = 0 (see [5, 6]). A full closed analytic solution (compatible with (57)) is giv en b y: ∆ ( m ) ( ξ , r ) = ∆ ( m ) i [1 + ∆ ( m ) i ] e 3 ξ/ 2 − ∆ ( m ) i , (71a) H ( ξ , r ) = H i e − 3 ξ/ 2 , (71b) S ( ξ , r ) = − 1 2 ∆ ( m ) ( ξ , r ) . (71c) so that solution curves remain for all ξ in the line S = − ∆ ( m ) / 2 that lies in the plane h Ω i = 1 (see figure 1). W e ha ve all v ariables fully determined. It is straightforw ard to verify that (71) are fully compatible and equiv alen t to analytic solutions given in terms of ` ( ct, r ) (see [2] and App endix C). Critical p oints on this case are: C 1 : source , C 3 : sink , (72) The phase space ev olution go es along the line [ S, − 2 S, 1], from C 1 to C 3 (resp ectiv ely , past and future global at- tractors). Since the critical p oint C 3 is characterized b y conditions of homogeneity ∆ ( m ) = S = 0, it corresp onds to the Einstein–de Sitter univ erse (dust FLR W model with zero spacial curv ature). 9 B. FLR W dust mo dels If ∆ ( m ) = S = 0 with unrestricted h Ω i , we ha ve an in v ariant set asso ciated with the particular homogeneous and isotropic sub–case of dust FLR W mo dels. The criti- cal p oints are in this case: C 3 : source , C 6 : sink , (73) Dep ending on the sign of h Ω i i − 1, the phase space evolu- tion of these mo dels will take place in the line [0 , 0 , h Ω i ]. If h Ω i i − 1 < 0 the evolution go es b etw een the source C 3 and the sink C 6 (past and future attractors). If h Ω i i − 1 > 0, it starts at C 3 and go es upw ards to h Ω i → ∞ to wards maximal expansion. In the collaps- ing stage, the p oint C 3 acts as a sink. If h Ω i = 1, we ha ve the Einstein–de Sitter universe and the evolution is constrained to the p oin t C 3 with co ordinates [0 , 0 , 1]. Also, we can identify the sink C 6 , whose co ordinates are [0 , 0 , 0], as Minko wski space–time in the Milne represen- tation. C. RSC’s The solution curv e corresponding to a RSC is also c har- acterized b y ∆ ( m ) = S = 0 with v arying h Ω i , thus the phase space ev olution of that particular curv e will also b e along the line [0 , 0 , h Ω i ], just like the FLR W sub– case. This is not surprising because a RSC is a privileged isotropic observer. The fact that the RSC’s hav e a separate evolution from the rest of the solution curves preven ts the criti- cal p oin t C 1 (a past and future attractor asso ciated with the initial and collapsing singularities) to b e global in all cases. How ev er, C 1 is a global attractor for all “off cen- ter” curves in configurations with RSC and for all curves in configurations lac king RSC (wormhole topology and v acuum limit). D. V acuum case It is known that the L TB metric (1) can describ e the Sc hw arzchild-Krusk al space–time in terms of the w orld lines of its radial geo desic test observ ers with 4-velocity u a = δ a 0 . The equation of motion of these geo desics is iden tical to (3) with M = M 0 , so that M 0 = 0 = ρ and the Ricci tensor v anishes, while the constant M 0 can b e identified with the “ Sc hw arzschild mass” (see [5, 6, 8]). The function 3 R is the scalar 3-curv ature of the hypersurfaces of simultaneit y of these geodesic con- gruences and the av erage scalar curv ature is related to their binding energy p er unit rest mass of the test parti- cle E i ( r ) = − K/ 2 = ( − 1 / 12) h 3 R i i Y 2 i = ( − 1 / 2) H 2 0 h k i i . This particular case can also b e constructed by assum- ing a Dirac delta distribution for ρ i so that a nonzero h ρ i i follo ws from the in tegral definition (10). Also, without resorting to a Dirac delta, this v acuum case follows from the 3+1 formalism in equations (34) b y setting ρ = 0, so a definition of a nonzero av erage densit y follows from the W eyl tensor comp onent in (7) and (27) as: E = − κc 2 6 h ρ i = − h m i H 2 0 2 = − M 0 / Y 3 i ` 3 (74) Since m = ρ = 0 but h m i > 0, as long as Γ is finite (see section X) we hav e the follo wing sufficient condition to c haracterize the Sch warzc hild-Krusk al v acuum: ∆ ( m ) = − 1 for all ( ξ , r ) (75) whic h defines an in v ariant set, since all solution curv es with initial condition ∆ ( m ) i = − 1 will necessarily be constrained to the plane ∆ ( m ) = − 1, parametrized b y [ S, h Ω i ]. The phase space evolution and critical p oints dep ends on the sign of h Ω i i , which dep ends in turn on the binding energy of the radial geo desics. F or zero binding energy h k i i = 0 = h Ω i i − 1 (the so–called Lema ˆ ıtre co ordinates for Sch warzsc hild space- time [8]) the critical points are the global future and past attractors: C 1 : source , C 2 : sink , (76) and phase space evolution is the line [ S, − 1 , 1], from the source C 1 to the sink C 2 . F or geodesics with negativ e binding energy: h k i i > 0 or h Ω i i > 1 (the so–called Novik o v coordinates [9]) the critical points are the same as those of the zero energy case, but C 2 is now a saddle. The phase space evolution of the curves b egin at the global past attractor (source) C 1 , approach the saddle C 2 and go up wards tow ards di- v erging h Ω i at the maximal expansion, terminating again at the global future attractor (sink) C 1 . F or geo desics with p ositive binding energy ( h k i i < 0 or 0 < h Ω i i < 1) the critical p oints are: C 1 : source , C 5 : sink , C 2 : saddle , C 4 : saddle . (77) In this case the curves begin at the global past attrac- tor (source) C 1 , terminate at the global future attractor (sink) C 5 , approaching the saddles C 2 and C 4 . This case is examined in subsection F of next section (see fig- ure 18). W e examine the phase space evolution of solution curv es of sev eral general case dust L TB configurations in the follo wing section. IX. L TB DUST CONFIGURA TIONS A. Expanding low density universe A low densit y p erp etually expanding dust configura- tion follo ws by selecting h m i i (0) < 1 and negativ e h k i i 10 for all r defined in an initial hypersurface 3 T i with an op en topology having a RSC. This leads to h Ω i i < 1 or to “hyperb olic” dynamics. Such a configuration can b e realized with the following initial v alue functions m i = 0 . 9 1 + tan 2 r , k i = − tan 2 r 1 + tan 2 r , Y i = H − 1 0 tan r, (78) so that we hav e an initial density clump ( − 1 < ∆ ( m ) i ≤ 0) but a curv ature void (0 ≤ ∆ ( k ) i < 2 / 3). It is easy to sho w that these initial conditions yield h Ω i i < 1 for all r , so that 0 < h Ω i < 1 for all ξ . These conditions also comply with the no–shell–crossing conditions (93), so dust lay ers monotonously expand p erp etually from an initial big bang. As shown by figure 2 the solution curves start at the past attractor or source C 1 (big bang), evolv e to wards the saddle C 3 and terminate at the future attractor given b y the line of sinks C 5 in the plane h Ω i = 0. The ev o- lution of the RSC is constrained to C 3 for all time. The function ∆ ( m ) remains negativ e for all ξ , thus we ha ve densit y clumps at all 3 T . FIG. 2: Lo w density configuration. Dust lay ers expand from the past attractor (source) C 1 , denoting the big bang singularit y , approac hes the saddle C 3 and drops tow ards the future attractor giv en b y the line of sinks C 5 . Notice that (78) implies that the evolution of the RSC is constrained to C 3 for all time. B. Lo w density void formation An interesting v ariation of a p erp etually expanding configuration, similar to the one presented in subsection A and figure 2, is that where an initial density clump ( − 1 < ∆ ( m ) i < 0) evolv es for all curv es into density v oids with ∆ ( m ) > 0. In this subsection we provide a numeric example of dust configurations discussed in reference [10] (see also [11, 12]). An initial dust clump transforming in to a v oid emerges b y choosing the follo wing initial v alue functions m i = m 01 + m 00 1 + α 2 0 tan 2 r , m 00 = 0 . 05 , m 01 = 0 . 9 , α 0 = 2 . 0 , k i = k 01 + k 00 1 + β 2 0 tan 2 r , k 00 = − 5 . 0 , k 01 = − 1 . 0 , β 0 = 0 . 5 , Y i = H − 1 0 tan r, (79) whic h comply with the no–shell–crossing conditions (93). The ev olution of this configuration is similar to the pre- vious one: dust lay ers also expand perp etually (“hy- p erb olic” dynamics) from an initial big bang and the h yp ersurfaces 3 T ha v e an op en top ology with a RSC. Ho wev er, the contrast function, ∆ ( m ) , now passes from negativ e to p ositive as ξ increases, indicating that ini- tial clumps transform in to v oids. This can b e seen very clearly if w e plot (see figure (3)) the normalized den- sit y profiles along the h yp ersurfaces 3 T from the function ρ/ρ c = m/m c = h m i (1 + ∆ ( m ) ) /m c , where the subindex c denotes ev aluation along the RSC marked by r = 0. The densit y as function of τ and r can be readily found b y solving numerically the system (60) for initial conditions (79). This transition from clumps to voids is depicted by figure 3. FIG. 3: Low density voids: densit y profiles The picture depicts the normalized radial density profile ρ/ρ c for different h yp ersurfaces 3 T marked by constant τ . Notice ho w an initial clump at τ = 0 evolv es into a voids profile for latter τ . This transition from clumps to voids can also b e seen from the form of the con trast densit y function ∆ ( m ) giv en 11 in terms of ξ and r , display ed in figure 4 which sho ws ho w ∆ ( m ) is initially negativ e but b ecomes positive for all solution curv es with r constan t (except the RSC at r = 0) as ξ increases. FIG. 4: Density con trast function The picture depicts the function ∆ ( m ) passing for all 0 < r < π / 2 from negative to p ositiv e indicating that an initial clump evolv es into a void profile. Notice ho w ∆ ( m ) = 0 at the RSC r = 0, and asymp- totically at r → π / 2 (or Y i → ∞ ). F or solution curves with r close to the RSC, ∆ ( m ) remains close to zero but is neverthe- less p ositive. The evolution of this configuration in phase space is sho wn by figure 5. Solution curves start at the past at- tractor (source) C 1 (big bang), approach the saddle C 3 and ev olv e to w ards the future attractor (line of sinks) C 5 in the plane h Ω i = 0. The RSC evolv es from the source C 3 to wards the sink C 6 . How ever, no w the curv es ev olv e in suc h a aw ay that all curv es (sav e the RSC) hit the line of sinks C 5 at p ositive v alues of ∆ ( m ) (passage from clumps to v oids). C. Re–collapsing high density univ erses Re–collapsing high density configurations follow b y se- lecting h m i i (0) > 1 and h k i i > 0 for all r , leading to h Ω i i ≥ 1 for all r , thus h Ω i ≥ 1 for all ξ . Dust lay- ers expand from an initial big bang, reach a maximal expansion and then collapse (the so–called “elliptic” dy- namics [2, 5, 6]). These configurations are compatible with 3 T i ha ving either a closed top ology with tw o RSC, an op en top ology with one RSC or a wormhole top ology without RSC (see App endix A). W e examine the closed and op en case separately b elow, and the wormhole case in section E. The case with spherical “closed” topology can b e ob- FIG. 5: Lo w density voids The evolution is similar to that of Fig 1, except that ∆ ( m ) b ecomes p ositive as all solution curv es drop to the line of sinks C 5 . tained with the following initial v alue functions: m i = m 01 + m 00 − m 01 1 + sin 2 r , m 00 = 3 . 9 , m 01 = 1 . 1 k i = k 01 + k 00 − k 01 1 + sin 2 r , k 00 = 2 . 9 , k 01 = 0 . 22 Y i = H − 1 0 sin r, (80) whic h comply with (91) at the turning p oint r = π / 2. The phase space evolution is illustrated in figure 6. The solution curves start at the past attractor (source) C 1 (big bang) and evolv e upw ards as h Ω i → ∞ , which cor- resp onds to maximal expansion giv en b y H → 0 and ξ → ln[ h Ω i i / ( h Ω i i − 1)]. The RSC also evolv es upw ards from the source C 3 . F or off center solution curv es this critical p oint is a saddle. Since the collapsing stage cannot b e describ ed with solution curves of system (51), w e can examine this stage with system (60) by plotting the curves [ S ( τ , r ) , ∆ ( m ) ( τ , r ) , h Ω i ( τ , r )], where S and h Ω i are de- fined by (44) and (43). As shown by figure 7, the curves come down wards from infinite h Ω i and S to the critical p oin t C 1 , which is now a future attractor or sink, asso ci- ated with a second curv ature singularity (“big crunch”). As in the expanding state, curves near the RSC approac h the saddle C 3 , while the RSC evolv es along the line [0 , 0 , h Ω i ] from infinite h Ω i to h Ω i = 1 at the sink C 3 . The case with open top ology with one RSC can b e 12 FIG. 6: Phase space ev olution of a high density re– collapsing univ erse with closed top ology . Solution curv es evolv e from the past attractor or source C 1 , approach the saddle C 3 and diverge as ξ → ln[ h Ω i i / ( h Ω i i − 1)], mark- ing the maximal expansion of dust lay ers. Only the range 0 ≤ r ≤ π / 2 is shown, with the RSC at r = 0 evolving from the source C 3 up wards. The off–center la yers approaching the saddle C 3 are those with r ≈ 0 while those further aw ay from C 3 are close to the “turning v alue” r = π / 2. Lay ers mark ed by π / 2 ≤ r ≤ π , including the second RSC at r = π , w ould ha ve iden tical ev olution as those sho wn in the figure, with those marked by r ≈ π near the saddle C 3 . FIG. 7: Phase space ev olution of the collapsing stage in the re–collapsing universe with closed top ology . Only the end stage of the collapse is describ ed by plotting the phase space v ariables by means of system (60). Solution curv es ev olve from infinite v alues of h Ω i and S to the p oint C 1 , which is now a future attractor or sink. The curves also approac h the saddle C 3 . Only the range 0 ≤ r ≤ π / 2 is sho wn, with the RSC at r = 0 evolving from infinite v alues of h Ω i to the C 3 , whic h is now a sink. La yers marked b y π / 2 ≤ r ≤ π , including the second RSC at r = π , would hav e iden tical ev olution as those shown in the figure, with those mark ed b y r ≈ π near the saddle C 3 . obtained with the following initial v alue functions: m i = 5 . 0 1 + tan 2 r , k i = 2 . 0 1 + tan 2 r , Y i = H − 1 0 tan r, (81) As in the case with closed top ology , dust lay ers emerge from a big bang singularity reach a maximal expansion and then re-collapse in a big crunc h, but no w dust lay ers far from the RSC at r = 0 re-collapse in very long time whic h b ecomes infinite in the limit r → π / 2. The ev olution of the solution curves of system (51) in phase space is sho wn b y figure 8. This ev olution is similar to the closed spherical case depicted in figure 6. Since initial conditions yield h Ω i i ≥ 1 w e hav e h Ω i ≥ 1 for all ξ . The curves (except the RSC) start at the past attractor or source C 1 (big bang), approac h the saddle C 3 and div erge upw ards as ξ → ln[ h Ω i i / ( h Ω i i− 1)], corresponding to maximal expansion ( H → 0). The RSC evolv es from the source C 3 to wards h Ω i → ∞ . FIG. 8: Phase space evolution of high density re– collapsing univ erse with open top ology . Since h Ω i → ∞ as ξ → h Ω i i / ( h Ω i i − 1), we plot arctan h Ω i , The curves start at the past attractor or source C 1 (big bang) approac h the saddle C 3 and div erge upw ards (maximal expansion). The RSC evolv es from the source C 3 to h Ω i → ∞ . The collapsing stage is not describ ed. As in the case with spherical top ology , the collapsing stage cannot describ ed with (51). W e use then the solu- tion curv es of (60) to plot the phase space v ariables in the collapsing stage in figure 9, whic h is v ery similar to figure 7: the critical p oint C 1 is now a future attractor or sink asso ciated with the collapsing singularit y (“big crunc h”), while curves near the RSC approac h the saddle C 3 and the RSC evolv es from infinite h Ω i down wards to h Ω i = 1 along the line [0 , 0 , h Ω i ]. 13 FIG. 9: Phase space evolution of the collapsing stage in the re–collapsing universe with op en top ology . Phase space v ariables are plotted by solving system (60). Solution curv es ev olve from infinite v alues of h Ω i and S to the p oint C 1 , which is no w a future attractor or sink. The RSC at r = 0 evolv es from infinite v alues of h Ω i to the sink C 3 . D. Structure formation scenario A very imp ortant type of configuration, asso ciated with structure formation scenarios, follows by ha ving a lo cal region around a RSC evolving from a big bang sin- gularit y and collapsing in to a black hole (second singular- it y), while the region “outside” expands p erp etually as a sort of “cosmic bac kground”. This type of configurations ha ve also b een examined with the L TB original v ariables (see [1, 5, 13, 14, 15, 16] and references quoted therein). The structure formation scenario requires that h Ω i i − 1 b e p ositive at the RSC and b ecomes negativ e for large r . It can b e constructed by choosing m i > 1 near the RSC that becomes m i < 1 asymptotically , while k i m ust b e p ositiv e at the RSC and become asymptotically negativ e. Suc h a choice will yield for re–collapsing curv es h Ω i i > 1 near the RSC, so that h Ω i ≤ 1 and diverges for ξ → ln[ h Ω i i / ( h Ω i i − 1)]. F or expanding curv es w e ha ve 0 < h Ω i i ≤ 1 and 0 < h Ω i ≤ 1 for all ξ . The following initial v alue functions fulfill these conditions: m i = m 01 + m 00 − m 01 1 + r 6 , m 00 = 1 . 5 , m 01 = 0 . 2 k i = k 01 + k 00 − k 01 1 + r 6 , k 00 = 0 . 5 , k 01 = − 0 . 5 Y i = H − 1 0 r , (82) W e hav e depicted in all the graphs the collapsing solu- tion curves in red and the p erp etually expanding ones in blue. The collapsing curves evolv e as the solution curves of figures 6 and 8, except that no w these curv es only co ver a bounded range of r around the RSC at r = 0, with expanding curv es filling the asymptotic range of r . F or the collapsing curves the function H obtained by solving the system (60) for initial conditions (82) c hanges sign (passes from p ositive to negative). This function is plotted in figure 10 showing collapsing and perp etually expanding lay ers resp ectively in red and blue colors. FIG. 10: The function H in a structure formation scenario. This function is obtained by solving the system (60) for initial conditions (82). F or re–collapsing curves with h Ω i i > 1 (red color) this function passes from p ositiv e (ex- panding) to negative (collapsing), while for p erp etually ex- panding curv es with h Ω i i ≤ 1 (blue) it remains alwa ys p osi- tiv e. It is also interesting to visualize the the dimensionless densit y m (see figure 11). F or collapsing curv es (red) around the RSC at r = 0 this function decreases from in- finite v alues as dust lay ers expand from the initial singu- larit y , reaches a minimal v alue (maximal expansion) and then increases again as curv es go into a collapsing singu- larit y . How ev er, for p erp etually expanding curves (blue) this function decreases monotonically from the initial big bang. The ev olution of the solution curves of (51) in phase space is more complicated than for previously presented configurations. Collapsing curv es (red) with h Ω i i > 1 ev olve as curves in re–collapsing high densit y configura- tions of figures 6 and 8, starting at the past attractor or source C 1 (big bang), approaching the saddle C 3 and div erging upw ards ( h Ω i → ∞ ) as ξ → ln[ h Ω i i / ( h Ω i i − 1)]. P erp etually expanding curv es (blue) with 0 < h Ω i i ≤ 1 resem ble the curv es of the low density universe of fig- ure 2, starting at the same past attractor or source C 1 (big bang) approac hing the saddle C 3 and descending to wards the future attractor given b y line of sinks C 5 . This is sho wn in figures 12 and 13. As shown by figure 12, all curves start at the past attractor C 1 (big bang), those near the RSC ( r = 0) ap- proac h the saddle C 3 (see figure 12) while those near the split from collapsing to expanding approach the saddle 14 FIG. 11: Densit y for a structure formation scenario. The initial big bang singularity is no w asso ciated to infinite densit y that decays as dust lay ers expand. Curves (red) with h Ω i i > 1 reach a minimal density (at maximal expansion) and then density reb ounds to infinity at the collapse. F or the remaining curves with h Ω i i ≤ 1 (blue) density just dilutes monotonously as τ gro ws. C 2 (see figure 13). All collapsing curv es end up div erging up wards h Ω i → ∞ as ξ → ln[ h Ω i i / ( h Ω i i − 1)] (maximal expansion). Blue curves with with 0 < h Ω i i ≤ 1 are p er- p etually expanding with 0 < h Ω i i ≤ 1 and so 0 < h Ω i ≤ 1 for all ξ . Expanding curves also start at the source C 1 (big bang), but do not approach the saddle C 3 , instead they approach the saddle C 2 and end up in the line of sinks C 5 . The RSC evolv es from the source C 3 up wards. Both the collapsing and expanding curves approach the saddle C 2 , which splits the former curves going upw ards to h Ω i → ∞ and the latter down wards tow ards the line of sinks C 5 . The curve with h Ω i i = 1 remains in the plane h Ω i = 1, ev olving from C 1 to wards C 2 . See figure 13. The collapsing curves b ehav e near the collapse as the curv es in figures 7 and 9. W e sho w in figure 14 the col- lapsing stage of the curv es plotted in figure 13, for a range of r around v alues where h Ω i i − 1 passes from p ositive (collapsing curv es) to negative (expanding curves). No- tice how collapsing curves end up falling into the critical p oin t C 1 , whic h is a future attractor or sink. As the curv es approac h the splitting saddle C 2 , the collapsing ones go to infinite h Ω i and S , while expanding ones fall in to the sinks C 5 . Because of the restricted domain of r that w e used to illustrate the b ehavior near C 2 , the fig- ure do esn’t show curv es near the RSC approaching the saddle C 3 and the RSC evolving from infinite h Ω i into C 3 , which is a sink for this worldline. ! arctan( ) arctan( ) " (m) arctan( ) S ! = 1 C 1 C 2 C 3 C 4 C 5 FIG. 12: Phase space picture of the structure forma- tion scenario. Red color corresp onds to collapsing curves with h Ω i i > 1, while p erp etually expanding curves with 0 < h Ω i i < 1 are depicted in blue. All “off center” curves start at the past attractor C 1 , p erp etually expanding curves descend into the future attractor giv en b y the line of sinks C 3 , collapsing curves approach the saddle C 3 and go upw ards to div erging h Ω i (maximal expansion b efore collapse). The sad- dle C 2 splits expanding from collapsing curv es. The curv e with h Ω i i = 1 (not shown) would go from C 1 to wards this saddle. ! = 1 ! " (m) S C 1 C 2 C 3 C 4 C 5 FIG. 13: Close up of the solution curves near the sad- dle C 2 .Red color corresponds to la yers with h Ω i i > 1 and blue ones to h Ω i i < 1. The plot only depicts solution curves that corresp ond to r v alues near the change from re–collapsing to p erpetually expanding b ehavior (where H b ecomes negativ e). Notice how this saddle splits collapsing curves going up wards to h Ω i → ∞ while expanding ones go down wards to line of sinks C 5 . Notice how the curv e corresp onding to h Ω i i = 1 remains in the plane h Ω i = 1, starting from the source C 1 and ending in the saddle C 2 . 15 FIG. 14: Collapsing stage for curves of figure 13. The curv es come from solving (60) for initial conditions (82). The end state of collapsing curves (red) is the point C 1 , which is now a sink. Collapsing curves go upw ards taking infinite v alues of h Ω i and S . Notice ho w the saddle C 2 splits col- lapsing curves going upw ards to h Ω i → ∞ while expanding ones (blue) go do wn wards to the future attractor given b y the sinks C 5 . E. Re–collapsing wormhole In dust configurations without a RSC (“wormhole” top ology), the 3 T can b e homeomorphic to either S 2 × R or to S 2 × S 1 . These configurations were examined in reference [17] (see also [5, 6]) and can be constructed with Y i ( r ) ha ving no zero es in all its domain. Since any c hoice of suc h function will necessary fulfill Y 0 i ( r ∗ ) = 0 for at least a “turning point” r = r ∗ , these configurations can only b e regular if they comply with the re gularit y condition (91), which implies a re–collapsing “elliptic” dynamics with h Ω i i − 1 > 0 (see App endix A). In the integral definitions of h ρ i and h 3 Ri in (10) and (11) w e to ok the RSC as the lo wer bound of the integrals. This is a boundary condition on these definitions ensuring that regularity at that RSC is fulfilled [2, 5, 6]. If there are no RSC then this boundary condition can be specified b y demanding that m i and k i v anish at an asymptotic v alue of r along the 3 T i . Such v alue can b e then tak en as the low er b ound for the integrals. How ever, while this solv es the mathematical problem, without a RSC the in terpretation of h ρ i and h 3 Ri as av erage functions b ecomes less clear. A choice of initial v alue functions for the re–collapsing w ormhole is given by: m i = m 00 (1 + α 0 ) 1 + α 0 sec 2 r , m 00 = 1 . 1 , α 0 = 5 . 0 k i = k 00 (1 + β 0 ) 1 + β 0 sec 2 r , k 00 = 0 . 48 , β 0 = 3 . 0 Y i = H − 1 0 sec r, (83) Notice that the “turning p oint” (or “throat”) in this ex- ample is lo cated at r = r ∗ = 0 with (91) taking the form h k i i (0) = 1. Also, both m i → 0 and k i → 0 as r → ± π / 2. Therefore, we can define h m i i and h k i i by integrals like (10) and (11) b y taking the low er in tegration limit at r = − π / 2 and integrating forward in r up to π / 2. The metric function Y = Y i ` obtained from (60) for this case is depicted in figure 15. FIG. 15: Curv ature radius of dust lay ers with w orm- hole top ology . Notice how for each h yp ersurface 3 T with τ constan t, Y = Y i ` has no zeros in its regular range and takes its minimal v alue at the “throat” at r = 0. Dust lay ers near the “throat” re–collapse in very a short interv al of τ , while la yers further aw ay re–collapse in m uch longer times. The evolution of the w ormhole configuration in phase space is depicted in figure 16. Only curves in the range 0 ≤ r < π / 2 are displa yed, since (b ecause of the sym- metric construction of the example) curv es in the range − π / 2 < r ≤ 0 b ehav e identically . This evolution is sim- ilar to that of high densit y re–collapsing universes, with all solution curv es starting at the past attractor or source C 1 asso ciated with an initial big bang singularity and ris- ing tow ards maximal expansion h Ω i → ∞ . The lack of a RSC is evident and this means that C 1 is a global past attractor. It is in teresting to notice how the “throat” and curves near it rapidly leap in to diverging h Ω i because they reac h maximal expansion and collapse very fast (see figure 15), but curves close to r = π / 2 sp end a long time very near the plane h Ω i = 1 and approach the saddle C 3 where 16 FIG. 16: Dust wormhole There is no RSC. Notice ho w curv es near the “throat” at r = 0 rapidly reach maximal ex- pansion with h Ω i diverging, while curves close to the asymp- totic v alue r = π / 2 take a long time to reach this stage and approac h the saddle p oint C 3 with ∆ ( m ) = S = 0. The lack of a RSC implies that C 1 is now a global past attractor. S = ∆ ( m ) = 0. W e to ok as example a v ery simple sym- metric form of 3 T i with only one turning p oint, but it is p ossible to choose any num b er of suc h p oints [17]. It is also p ossible to construct a configuration whose hy- p ersurfaces 3 T are homeomorphic to torii S 2 × S 1 (see [5, 6]). The collapsing stage for w ormhole configurations is v ery s imilar to the stages depicted in figures 7 and 9. W e show in figure 17 the expanding and collapsing stage for curves near the “throat” r = 0 of this configuration. It is easy to recognize some of the curv es plotted in figure 16, which contin ue to wards infinite v alues of h Ω i and S , in order to plunge down wards to end in the critical p oint C 1 , which is a future attractor or sink in this stage. F. V acuum limit The evolution in phase space takes place in the plane ∆ ( m ) = − 1. As we commen ted in the previous section, for zero binding energy h Ω i i = 1 (Lema ˆ ıtre co ordinates [8]) it go es from source C 1 to sink C 2 , while for negative binding energy h Ω i i > 1 (No viko v co ordinates [9]) it go es from source C 1 , approaches saddle C 2 and go es upw ards to h Ω i → ∞ . The phase space evolution in this case is qualitativ ely analogous to that of a wormhole top ology , whic h is not surprising since the Sch warzsc hild–Krusk al space–time has no RSC, hence C 1 will be a global past and future attractor. W e provide here an example with p ositive binding en- ergy h Ω i i < 1. As shown in figure 18 the solution curves emerge from the global past attractor source C 1 , some FIG. 17: Dust wormhole: collapsing stage The figure describ es the expanding and collapsing stages. The curves of figure 16 near the “throat” r = 0 can b e recognized and seen going upw ards into infinite v alues of h Ω i and S , all in order to plunge down wards into C 1 . Notice that the lack of a RSC mak es this critical p oint a global past and future attractor. curv es approach the saddle C 2 and some approac h the saddle C 4 and all terminate in the global future attractor or sink C 5 . FIG. 18: V acuum case The solution curves corresp ond to the Sc hw arzsc hild–Krusk al space–time in coordinates giv en b y radial geodesic congruence that is perp etually expanding with p ositiv e binding energy (thus h Ω i i < 1). Solution curves are confined to the plane ∆ ( m ) = − 1. The curves evolv e from the global past attractor or source C 1 to wards the global future attractor or sink in C 5 with S = 0, approaching the saddles C 2 and C 4 . 17 X. DISCUSSION AND CONCLUSION W e hav e examined the class of dust L TB mo dels in terms of the techniques generically known as “dynamical systems” as they evolv e in a 3–dimensional phase space parametrized b y the v ariables [ S, ∆ ( m ) , h Ω i ], whose ph ys- ical and geometric in terpretation is in tuitive and straight- forw ard. W e can say in general that if initial condi- tions are selected so that “shell crossing” singularities are a voided (condition (68)), then the phase space v ariables [ S, ∆ ( m ) ] will b e b ounded. Also, the sign of h Ω i i ( r 0 ) − 1 for an y given solution curve r = r 0 determines if it will ev olve for all ξ in the domain of h Ω i ( ξ , r 0 ) − 1 that has the same s ign. The v arious inv arian t subspaces, partic- ular cases and phase space evolution for representativ e dust configurations has b een presented in detail in the previous sections. Except for the worldlines asso ciated with RSC’s, all “off center” solution curves b egin at a past attractor giv en b y a rep elling source C 1 , c haracterized b y S = 1 / 2, h Ω i = 1 and ∆ ( m ) = − 1 and asso ciated with a big bang initial singularity . It is easy to show that C 1 is indeed the only past attractor for all ‘off cen ter” solution curv es: the asymptotic past regime for all curves is giv en by the limit ξ → ∞ (or ` → 0), so from (63a) we ha ve h Ω i → 1 in this limit. On the other hand, we can write ∆ ( m ) and ∆ ( k ) as ∆ ( m ) = m i h m i i Γ − 1 , ∆ ( k ) = k i h k i i Γ − 1 . (84) But, from the analytic forms giv en in App endix C and for finite and regular m i and k i , we ha ve in the limit ` → 0 the following asymptotic b ehavior: Γ ∼ ` − 3 / 2 for h Ω i i = 1 and Γ ∼ ` − 1 for h Ω i i 6 = 1. Thus, from (84), w e ha ve ∆ ( m ) → − 1 and ∆ ( k ) → − 1 in this limit, but then, from (66) and bearing in mind that h Ω i → 1, w e ha ve S → 1 / 2. As shown in figures 7, 9, 14 and 17, all re–collapsing curv es also terminate in C 1 , whic h is in this case a fu- ture attractor for these curves and is asso ciated with the collapsing (“big crunch”) singularity . It is only the fact that RSC’s do not evolv e nor terminate (collapse) in C 1 what prev ents this attractor to be a global one, though in configurations lacking a RSC (w ormhole top ology and v acuum limits) this attractor is indeed global. In order to prov e that near C 1 w e hav e a self–similar regime we use the definitions (6), (17), (36) and (44). The condition S = 1 / 2 can b e reduced to Y C 1 = r (1 + ϑ ) 2 / 3 , Y 0 C 1 = 1 + ϑ/ 3 (1 + ϑ ) 1 / 3 , with ϑ ≡ τ r , (85) where τ = H 0 ct and we hav e eliminated t wo functions of the form a ( ct ) and b ( r ) by trivially re–lab eling the time and radial coordinates. Since h Ω i = 1 implies h k i i = 0 = h 3 R i i , the metric (1) near C 1 tak es the form of the self–similar sub–case of L TB metrics [18] ds 2 = − dτ 2 + A 2 ( ϑ ) dr 2 + r 2 B 2 ( ϑ ) ( dθ 2 + sin 2 θ dφ 2 ) , with A ( ϑ ) = Y 0 C 1 , r B ( ϑ ) = Y C 1 . (86) The fact that ∆ ( m ) = − 1 do esn’t imply that C 1 is a Sc hw arzschild v acuum, but that Γ → ∞ as ` → 0. Ev al- uating Γ for the self–similar metric forms (85) and (86) with Y i = Y C 1 ( t i , r ), we get Γ = 3 + ϑ 3 + ϑ i 1 + ϑ i 1 + ϑ (87) with ϑ i = τ i /r , which shows how Γ → ∞ at the lo cus of the big bang singularit y ϑ = − 1 where Y C 1 = r B ( ϑ ) → 0. Notice that the self–similar case has no shell crossing singularit y since Γ v anishes at ϑ = − 3 but the evolution range is ϑ > − 1 (for collapsing stages we would hav e Y and Y 0 giv en by the same expressions as (85), but with a min us sign, so that the “big crunch” happ ens b efore the lo cus of Γ = 0). Therefore, we can consider the source C 1 as marking a self–similar limit and since all “off center” solution curv es start at this critical point (associated with a big bang singularit y), w e can sa y that dust L TB mo dels hav e a self–similar b ehavior near this singularity . Since collaps- ing curv es in re–collapsing configurations also end in C 1 (a sink), w e ha v e self–similar b ehavior near the collapsing singularit y . Another in teresting feature worth remarking is the role of the saddle C 2 in the “structure formation” configura- tions discussed in section IX-D. As shown by figures 13 and 14, this saddle splits collapsing curves h Ω i i > 1 from the expanding ones 0 < h Ω i i < 1. The former evolv e bac k into C 1 (no w a sink) while the latter fall into the line of sinks C 5 . This b ehavior denotes a basic instabil- it y of the inv ariant subspace of the “parab olic” evolution h Ω i = 1, since for an y solution curv e in this space any arbitrarily small perturbation δ  1 on initial condi- tions h Ω i i = 1 + δ will trigger a radically different “repul- siv e” ev olution aw ay from h Ω i = 1, either to a collapsing regime ( δ > 0) or to p erp etual expansion ( δ < 0). The splitting saddle C 2 is mark ed by ∆ ( m ) = − 1, S = − 1 / 3 and h Ω i = 1, it can b e shown from the ana- lytic solutions in App endix C and from numerical tests that Γ is finite as the curves approach these v alues (Γ div erges for later times as the collapsing curves actu- ally collapse). Therefore, the fact that ∆ ( m ) → − 1 for curves approaching this saddle means that these curv es exp erience near C 2 unstable conditions similar to those of a Sch w arzschild–Krusk al v acuum parametrized b y Lema ˆ ıtre coordinates made with radial geo desics with zero binding energy ( h k i = 0 or h Ω i = 1). In fact, the p oin t C 2 represen ts the global future attractor (sink) for this marginally b ound v acuum configuration. The asso ciation of the saddle C 2 with a Sc hw arzschild– Krusk al v acuum roughly conforms with an intuitiv e pic- ture. A v ery rough, hand w a ving, description of condi- tions near C 2 could b e: collapsing dust lay ers “retreat” 18 in wards tow ards a smaller collapsing region, while ex- panding lay ers “retreat” outw ards in a fast expansion, th us w e ha ve a sort of unstable “ev acuation” effect that is roughly similar to a Sch warzsc hild–Krusk al v acuum with the Sch w arzschild mass being the accum ulated effective mass M 0 = (1 / 2) h m i i Y 3 i of the collapsing lay ers. It is kno wn [5, 6] that in this type of configuration the more “external” lay ers in the collapsing region ( h Ω i i > 1) take infinite time to collapse, so that the resulting blac k hole p erp etually tak es accretion from these “border” dust la y- ers. How ev er, we can assume that the mass con tribution from this accretion is small and lengthy to come, so that once most “inner” la yers with h Ω i i > 1 hav e collapsed the expanding observers close to h Ω i i = 1 do p erceiv e a sort of appro ximate Sch w arzschild–Krusk al v acuum with zero binding energy . W e feel that the qualitative numerical treatment pre- sen ted in this article can serve for studying applications of astroph ysical interest of dust L TB solutions. P erhaps for this purp ose the system (60) can b e more useful than the dynamical system (51). It is also interesting to see this metho dology generically applied to generalizations of these solution, such as L TB–de Sitter (dust with cosmo- logical constant) or the most general source compatible with the L TB metric: a mixture of dust and an inho- mogeneous and anisotropic fluid, which can resp ectively mo del cold dark matter and dark energy . Extensions of this work along these lines are presently under consider- ation. Ac knowledgmen ts The author ac knowledges support from Instituto de F ´ ısica, Universidad de Guana juato. APPENDIX A: Regularity a nd singularities in dust L TB solutions An extensiv e discussion of regularity conditions and imp ortan t geometric features of dust L TB solutions can b e found in [5, 6]. These features were also discussed with the v ariables of this article in [2]. W e present here a brief review. Regular Symmetry Centers (RSC) and top ology of the 3 T i L TB dust solutions admit up to tw o RSC’s, which are the regular w ordlines made up by the space–time evolu- tion of the fixed p oin t of the rotation group SO(3) where the orbits ha ve zero area. In terms of the initial v alue functions, the RSC is the comoving worldline marked by a zero of the function Y i ( r ). Since the freedom to choose a radial co ordinate in the L TB metric means that this function can b e arbitrarily chosen, it is useful to select it as Y i = H − 1 0 f ( r ), where f ( r ) is a (at least) C 2 real func- tion that describes in a simple manner the top ological configuration or homeomorphic class for an initial 3 T i . W e ha ve the following choices: • “Op en top ology”: 3 T i homeomorphic to R 3 , f ( r ) = tan r 0 ≤ r < π / 2 One RSC at r = 0 . (88) . • “Closed top ology”: 3 T i homeomorphic to S 3 , f ( r ) = sin r 0 ≤ r ≤ π Two RSCs at r = 0 , π . (89) . • “W ormhole”: 3 T i homeomorphic to S 2 × R , f ( r ) = sec r − π / 2 < r < π / 2 Zero RSCs . (90) . F or the op en top ology f ( r ) is a monotonously increasing function. In fact, we could ha ve just selected f = r for simplicit y , but the c hoice f = tan r allows one to examine the asymptotic behavior asso ciated with Y i → ∞ in solution curves of (51). In the w ormhole topology the function Y i m ust b e p os- itiv e and without zero es. Numerical tests show that monotonous functions without (at least) one “turning v alue” (zero of Y 0 i ) yield finite prop er radial distances along the 3 T i as r → ±∞ (dep ending on the c hoice of Y i ). This indicates that suc h functions are inadequate c hoices of the radial coordinate, as they do not pro vide a full analytic extension of the manifold. Therefore, we select Y i > 0 with at least a turning v alue r = r ∗ within the open range of r . In the choice (90) this v alue is r ∗ = 0 (the “throat” of the wormhole). In the “closed” topology w e m ust ha v e a turning v alue, whic h with the choice (89) is r ∗ = π / 2. Since in b oth cases (closed and w ormhole top ologies) w e ha ve Y 0 i ( r ∗ ) = 0, regularity of the metric function g rr Y 0 (1 − K ) 1 / 2 = ` Γ Y 0 i [1 − H 2 0 h k i i Y 2 i ] 1 / 2 requires (see [2, 5, 6] and [19]) that 3 R i b e selected so that K ( r ∗ ) = 1, or equiv alently: 1 6 h 3 R i i ( r ∗ ) = H 2 0 h k i i ( r ∗ ) = 1 Y 2 i ( r ∗ ) . (91) Giv en a choice of homeomorphic class for the 3 T i all reg- ular 3 T hav e the same class (see [2]). All configurations examined numerically corresp onding to closed or w orm- hole top ologies comply with (91). 19 Singularities L TB dust solutions present tw o t yp es of curv ature sin- gularities marked b y ( ct, r ) v alues such that: ` ( ct, r ) = 0 , “bang” or “crush” , (92a) Γ( ct, r ) = 0 , “shell crossing” (92b) The “bang” or “crush” singularities are an inheren t fea- ture and cannot b e av oided, though it is alwa ys p ossible to select initial conditions so that “shell crossing” do not arise (either for all the ev olution time or for a given range t ≥ t i ). Initial conditions to a void unph ysical shell crossing sin- gularities w ere given in terms of the original L TB v ari- ables by Hellaby and Lake [7] (see [5, 6]), and in terms of the initial v alue functions defined in this article by [2]. F or broadly generic configurations that exclude very sp ecial (and unstable) features, we hav e: for 0 < h Ω i i ≤ 1 : − 1 < ∆ ( m ) i ≤ 0 , − 2 / 3 < ∆ ( k ) i < 2 / 3 , (93a) for h Ω i i > 1 : − 1 < ∆ ( m ) i ≤ 0 , − 2 / 3 < ∆ ( k ) i ≤ 0 , 2 π P i  ∆ ( m ) i − 3 2 ∆ ( k ) i  ≥ [ P i Q i − 1] ∆ ( m ) i +  3 2 P i Q i − 1  ∆ ( k ) i ≥ 0 , , (93b) where P i = √ 2 − x x 3 / 2 , Q i = arccos(1 − x ) − p x (2 − x ) , x ≡ h 3 R i i κc 2 ρ i = 2 h k i i h m i i = 2[ h Ω i i − 1] h Ω i i . Equations (93) clearly show how initial conditions c haracterized by densit y and 3–curv ature clumps (see (15)) lead generically to av oidance of a shell crossing sin- gularit y if 0 < h Ω i i ≤ 1, and as a corollary , to finite and b ounded v alues of ∆ ( m ) and S along solution curves. This regularity is more stringent for the case h Ω i i > 1. Ho wev er, even if (93) hold at the 3 T i it is still p ossible to ha ve densit y or 3–curv ature voids forming along the evo- lution of the curv es, all without violating (68). This was discussed by [10], but their metho dology is too compli- cated. It is more clearly seen here from the fact that the scaling laws (23) and (24) do not, necessarily , forbid ∆ ( m ) and/or ∆ ( k ) from rev ersing the sign of the initial ∆ ( m ) i and ∆ ( k ) i . The configuration discussed in section IX–B is a numerical example of this type of “clump turning in to v oid” configuration. APPENDIX B: An equiv alent dynamical system W e ha ve examined the evolution of solution curves of the autonomous system (51) in the 3–dimensional phase space parametrized b y [ S, ∆ ( m ) , h Ω i ]. How ever, even if (51) has only deriv atives with resp ect to the ev olution parameter ξ , the differen tial equations are still PDE’s. W e pro vide b elo w formal proof that suc h a system of autonomous evolution PDE’s is equiv alen t to a system of ODE’s with initial conditions restricted b y the fulfillment of space-like constrain ts. Let F ( r ) ⊂ R and F ( ξ ; r ) ⊂ R b e, resp ectively , the reg- ularit y range of r and the maximal range of extendibility of ξ for a giv en r , the system (51) can b e asso ciated with the flow Φ ξ : R 3 → R 3 of the v ector field ∂ /∂ ξ given by Φ ξ ( ~ X i ) = ~ X , (95) where ~ X i : F ( r ) → R 3 is the initial state represented as a curv e in R 3 giv en parametrically as ~ X i ( r ) = [ S i ( r ) , ∆ ( m ) i ( r ) , h Ω i ( r ) i ] , (96) and ~ X : F ( ξ ; r ) × F ( r ) → R 3 is a solution surface as- so ciated with ~ X i b y means of (95). This surface can b e represen ted parametrically as the following surface in R 3 : ~ X ( ξ , r ) = [ S ( ξ , r ) , ∆ ( m ) ( ξ , r ) , h Ω( ξ , r ) i ] , (97) so that ~ X i ( r ) = ~ X (0 , r ). Given a point ~ X i ( r 0 ) ∈ R 3 there is a unique orbit or solution curve C r 0 : F ( ξ ; r 0 ) → R 3 of (51) which can b e represented parametrically as: C r 0 ( ξ ) = ~ X ( ξ , r 0 ) , (98) so that C r 0 (0) = ~ X i ( r 0 ). The set of solution curv es C r ( ξ ) for all r ∈ F ( r ) generates the solution surface ~ X ( ξ , r ). In order to relate (51) with a system of ordinary dif- feren tial equations, w e notice that ev ery solution curve (98) associated with a fixed r = r 0 is a solution of the follo wing system: dS 0 dξ = S 0 (3 S 0 − 1) + 1 2  ∆ ( m ) 0 + S 0  h Ω i 0 , (99a) d ∆ ( m ) 0 dξ = 3 S 0 h 1 + ∆ ( m ) 0 i , (99b) d h Ω i 0 dξ = h Ω i 0 [ h Ω i 0 − 1] , (99c) where the functions ( S 0 , ∆ ( m ) 0 , h Ω i 0 ) : F ( ξ ; r 0 ) → R are giv en by S 0 = S ( ξ , r 0 ) , (100a) ∆ ( m ) 0 = ∆ ( m ) ( ξ , r 0 ) , (100b) h Ω i 0 = h Ω( ξ , r 0 ) i , (100c) 20 and comply with the initial conditions: [ S 0 (0) , ∆ ( m ) 0 (0) , h Ω i 0 (0)] = ~ X i ( r 0 ) . (101) The fact that every solution curve of (51) will b e a solu- tion of a system lik e (99) for initial conditions giv en as in (101), suggests considering a system of ordinary dif- feren tial equations asso ciated with the same flo w as (95) ds dξ = s (3 s − 1) + 1 2 ( δ + s ) ω , (102a) dδ dξ = 3 s [ 1 + δ ] , (102b) dω dξ = ω [ ω − 1] , (102c) for an initial state: ~ x 0 = [ s (0) , δ (0) , ω (0)] = [ s 0 , δ 0 , ω 0 ] , (103) and the solution surface given by a parametric form sim- ilar to (97): ~ x = [ s ( ξ , s 0 ) , δ ( ξ , δ 0 ) , ω ( ξ , ω 0 )] , (104) so that ~ x 0 = ~ x (0). Since the solution surface of this system can dep end smo othly on initial conditions ~ x 0 , eac h one of the functions ( s, δ, ω ) solving (102) is a one–parameter family of functions. In principle, there are many wa ys in which the initial state (103) can b e parametrized. In particular, this state can b e given as any curve [ s 0 ( α ) , δ 0 ( α ) , ω 0 ( α )] in R 3 that intersects the solution curves ~ x ( ξ , ~ x 0 ) only in one p oint, hence w e can alwa ys find suitable one–parameter functions s 0 ( α ) , δ 0 ( α ) , ω 0 ( α ) such that ~ x is expressible as: ~ x ( ξ , α ) = [ s ( ξ , s 0 ( α )) , δ ( ξ, δ 0 ( α )) , ω ( ξ, ω 0 ( α ))] . (105) If, in particular, the initial state (103) is parametrized as: ~ x 0 ( α ) = [ s 0 ( α ) , δ 0 ( α ) , ω 0 ( α )] = ~ X i ( α ) , (106) with α ∈ F ( r ) , then for every α there is a C ∞ inclusion map ψ : R 3 → R 3 mapping every solution curve of (51) giv en by (98) for a fixed r ∈ F ( r ) to a solution curve of the system (102) giv en by (105) with a fixed α = r . That is ψ ( ~ X ( ξ , r )) = ~ x ( ξ , α ) , (107) and so, the solution surfaces of (51) are subsets of the so- lution surfaces of (102) obtained by restricting the initial states by the constraint (106). Notice that (102) admits many solution curv es that do not comply with (106), and so are incompatible with the constraints (53)–(56c) and so they are incompatible with the evolution of an inhomogeneous dust source (L TB solution) complying with Einstein’s field equations ex- pressed in terms of the 3+1 decomp osition [4] under the regularit y assumptions that w e hav e made (see Appendix A). Ho wev er, the solution curv es of (102) satisfying (106) are equiv alent (under (107)) to solution curves of (51) asso ciated with the flow (95). Since (102) is a system of autonomous ODE’s, geometric features such as criti- cal p oin ts and inv ariant spaces are w ell defined, so the qualitativ e analysis of (51) can b e p erformed on (102) as long as we only consider its solution surfaces satisfying (106). In other words, w e can say that the dynamical sys- tem asso ciated to inhomogeneous dust L TB solutions is a sort of “reduced” version of (102) with its initial states restricted by (106). Since, as we hav e shown in sections IV and V, the space–lik e constrain ts are satisfied at all space–lik e slices once they hold at the initial 3 T i , the in- corp oration of the radial dep endence of the solutions only requires this restriction on the initial states. In order to simplify the notation we found it useful to k eep the same names for the v ariables S, ∆ ( m ) , h Ω i and to keep referring to (51), with the understanding that the system we are really considering is (102), thus ev ery men tion of solution curves and surfaces of (51) refers to solution curves and surfaces of (102) mapp ed b y (107) and complying with (106). App endix C: Analytic solutions Analytic solutions for L TB dust solutions follo w from solving the ev olution equation (3): ˙ Y 2 = 2 M Y − K. Using the v ariables defined in sections I I I and IV, this equation b ecomes equation (59):  ∂ ` ∂ τ  2 = H 2 i  h Ω i i ` − ( h Ω i i − 1)  . where: ` = Y Y i , τ = H 0 ct 2 M = h m i i H 2 0 Y 3 i , K = h k i i H 2 0 Y 2 i , h Ω i i = h m i i H 2 i = h m i i h m i i − h k i i , h Ω i i − 1 = h k i i H 2 i = h k i i h m i i − h k i i , W e summarize the solutions of this equation b elow P arab olic solutions W e hav e: h Ω i i = 1, so that h k i i = 0 and H 2 i = h m i i . These are the only ones that can b e giv en explicitly as ` = ` ( τ , r ) from: τ = τ i ± 2 ( ` 3 / 2 − 1) 3 H i , (108) 21 where the ± sign describ es expanding (+) and collaps- ing ( − ) dust lay ers, resp ectively characterized b y τ > τ bb , ` ,τ > 0 and τ < τ bb , ` ,τ < 0. Notice that the “big bang time” τ bb = H 0 ct bb is no longer an indep endent function but m ust b e found b y setting ` = 0 in (108): τ bb = τ i ∓ 2 3 H i , (109) where now expanding and collapsing la y ers resp ectively corresp ond to “ − ” and “+”. Hyp erb olic solutions W e hav e: 0 < h Ω i i < 1, so that h k i i < 0 and H 2 i = h m i i + |h k i i| . In this case we cannot obtain ` = ` ( τ , r ) in closed explicit form, but in parametric form: [ ` ( η , r ) , τ ( η , r )] or implicitly as τ = τ ( `, r ). The solutions are given in parametric form as: ` ( η , r ) = h Ω i i [cosh η − 1] 2 (1 − h Ω i i ) , (110a) τ ( η , r ) = τ i ± h Ω i i [(sinh η − η ) − (sinh η i − η i )] 2 H i (1 − h Ω i i ) 3 / 2 , (110b) where η i = arccosh  2 h Ω i i − 1  , sinh η i = 2 p 1 − h Ω i i h Ω i i . (111) and ± denote expanding (+) and collapsing (–) lay ers. The co ordinate lo cus of the central singularity (the “big bang time” function τ bb ) can be obtained b y setting η = 0 in (110b), ( so that ` = 0) leading to τ bb = τ i ∓ h Ω i i (sinh η i − η i ) 2 H i (1 − h Ω i i ) 3 / 2 , (112) where no w “–” and “+” respectively correspond to ex- panding and collapsing lay ers. Elliptic solutions W e hav e: h Ω i i > 1, so that h k i i > 0 and H 2 i = h m i i − h k i i . The parametric solutions are ` ( η , r ) = h Ω i i [1 − cos η ] 2( h Ω i i − 1) , (113a) τ ( η , r ) = τ i + h Ω i i [( η − sin η ) − ( η i − sin η i )] 2 H i ( h Ω i i − 1) 3 / 2 , (113b) where η i = arccos  2 h Ω i i − 1  , sin η i = 2 p h Ω i i − 1 h Ω i i . (114) As with the hyperb olic case, the lo cus of the central sin- gularit y (the “big bang time”, τ bb ), follows by setting η = 0 in (113b), we get τ bb = τ i − h Ω i i ( η i − sin η i ) 2 H i ( h Ω i i − 1) 3 / 2 , (115) Dust lay ers in elliptic solutions reach a maximal expan- sion as ` ,τ = 0 in (59), corresp onding to ` = ` max = h Ω i i h Ω i i − 1 , η = π , (116a) = τ bb + π h Ω i i 2 H i ( h Ω i i − 1) 3 / 2 , (116b) F or ` ,τ < 0 dust lay ers collapse and so a collapsing phase is given b y parametric solution (113b) for π < η < 2 π ). The co ordinate lo cus of the collapsing singularity (“big crunc h”), τ bc , follo ws from setting η = 2 π in (113b): τ bc = τ bb + π h Ω i i H i ( h Ω i i − 1) 3 / 2 . (117) Since the expanding and collapsing singularities are resp ectiv ely mark ed b y η = 0 and η = 2 π , then the la yers’ ev olution is symmetric with respect to η = π marking the maximal expansion (116). How ever, as shown by (115), (116) and (117), the maximal expansion time dep ends on the IVF’s h Ω i i and H i , so it does not coincide (in general) with any 3 T (hypersurface of constant τ ) and the ev olution for general IVF’s h Ω i i , H i is not “time- symmetric” with resp ect to the maximal expansion. While for the h yp erb olic and elliptic cases we do not ha ve an explicit closed form solution ` = ` ( τ , r ), it is p ossible to use the parametric solutions presented so far to visualize graphically the function ` ( τ , r ). This can b e done by plotting the 3–dimensional parametric surfaces: [ τ ( η , r ) , r, ` ( η , r )] , (118) where the radial dep endence is completely specified by the IVF’s that enter in (110)–(111) and (113)–(114). [1] A Krasinski, Physics in an Inhomo gene ous Universe , (Cam bridge Univ ersity Press, Cambridge, 1998). [2] R.A. Sussman and L. Garc ´ ıa–T rujillo, Class. Quantum 22 Gr av , 19 , 2897 – 2925, (2002). [3] J W ainwrigh t and G F R Ellis (Eds.), Dynamic al Systems in Cosmolo gy (Cambridge Universit y Press, Cambridge, 1997). [4] G. F. R. Ellis and H. v an Elst, Cosmolo gic al Mo dels , Carg ` ese Lectures 1995. LANL e–print arXiv: gr-qc/9812046 [5] D. R. Matrav ers and N. P . Humphreys, Gen. R el. Gr av. , 33 , 531-552, (2001). [6] N.P . Humphreys, R. Maartens and D. R. Ma- tra vers, R egular spheric al dust sp ac etimes , LANL e–print arXiv: gr-qc/9804023v1 [7] C. Hellaby and K. Lak e, Ap. J. 290 , 381, (1985). See Ap. J. 300 , 461, (1986) for erratum of this pap er. [8] L.D. Landau and E.M. Lifshitz, The Classic al The ory of Fields, Course of The or etic al Physics V olume 2 . F ourth Revised English Edition. P ergamon Press, 1975. See § 102 to § 104. [9] C Misner, K. Thorne and J.A. Wheeler, Gravitat ion , F reeman Press, 1972. [10] N. Mustapha and C. Hellaby , Gen. R el. Gr av. , 33 , 455- 477, (2001). [11] W. Bonnor and A. Chamorro, Ap. J. 361 , 21, (1990). [12] K. Bolejk o, A. Krasinski and C. Hellaby , Mon. Not. R oy. Astr on. So c. , 362 , 213-228, (2005). [13] A. Krasinski and C. Hellaby , Ph ys.Rev. D65 (2002) 023501 [14] A. Krasinski and C. Hellaby , Ph ys.Rev. D69 (2004) 023502 [15] A. Krasinski and C. Hellaby , Ph ys.Rev. D69 (2004) 043502 [16] A. Krasinski and C. Hellaby , Ph ys.Rev. D73 (2006) 023518 [17] C. Hellaby , Class. Quantum Gr av. , 4 , 635, (1987). [18] B.J. Carr and A.A. Coley , Phys.R ev. D62 , 044023, (2000) [19] W. Bonnor, Class. Quantum Gr av. , 2 , 781, (1985).