A dynamical system approach to inhomogeneous dust solutions
AA dynamical system approach to inhomogeneous dust solutions.
Roberto A. Sussman ‡∗ Instituto de F´ısica, Universidad de Guanajuato, Loma del Bosque 103, Leon,Guanajuato, 37150, M´exico. ‡ On sabbatical leave from Instituto de Ciencias Nucleares,Universidad Nacional Aut´onoma de M´exico (ICN-UNAM), A. P. 70–543, 04510 M´exico D. F., M´exico.
We examine numerically and qualitatively the Lemaˆıtre–Tolman–Bondi (LTB) inhomogeneousdust solutions as a 3–dimensional dynamical system characterized by six critical points. One of thecoordinates of the phase space is an average density parameter, (cid:104) Ω (cid:105) , which behaves as the ordinaryΩ in Friedman-Lemaˆıtre–Robertson–Walker (FLRW) dust spacetimes. The other two coordinates, ashear parameter and a density contrast function, convey the effects of inhomogeneity. As long as shellcrossing singularities are absent, this phase space is bounded or it can be trivially compactified. Thisspace contains several invariant subspaces which define relevant particular cases, such as: “parabolic”evolution, FLRW dust and the Schwarzschild–Kruskal vacuum limit. We examine in detail the phasespace evolution of several dust configurations: a low density void formation scenario, high densityre–collapsing universes with open, closed and wormhole topologies, a structure formation scenariowith a black hole surrounded by an expanding background, and the Schwarzschild–Kruskal vacuumcase. Solution curves (except regular centers) start expanding from a past attractor (source) in theplane (cid:104) Ω (cid:105) = 1, associated with self similar regime at an initial singularity. Depending on the initialconditions and specific configurations, the curves approach several saddle points as they evolvebetween this past attractor and other two possible future attractors: perpetually expanding curvesterminate at a line of sinks at (cid:104) Ω (cid:105) = 0, while collapsing curves reach maximal expansion as (cid:104) Ω (cid:105) diverges and end up in sink that coincides with the past attractor and is also associated with selfsimilar behavior. PACS numbers: 12.60.Jv, 14.80.Ly, 95.30.Cq, 95.30.Tg, 95.35.+d, 98.35.Gi
I. INTRODUCTION
Inhomogeneous dust solutions with spherical symme-try are among the oldest, simplest and most useful ex-act solutions of Einstein’s equations. These solutionswere initially derived independently by Lemaˆıtre (1933)and Tolman (1934) and then re-derived by Bondi (1947),hence they are known as the Lemaˆıtre–Tolman–Bondi(LTB) solutions (see [1] and references quoted thereinfor a comprehensive review on these solutions and theirapplications).In practically all applications of these solutions thestandard original variables in which they were derivedare used, though alternative variables amenable to an“initial value” treatment have been proposed (see [2] andreferences quoted therein). While not the only type ofparametrization for these solutions, we find these vari-ables particularly useful for a numerical treatment.We propose in this paper to examine these modelswithin the framework of numeric and qualitative tech-niques known generically as “dynamical systems” [3].This approach requires as a first necessary step express-ing the evolution equations as a proper system of first or-der autonomous ODE’s (ordinary differential equations).This necessary requirement would, apparently, excludeinhomogeneous LTB solutions because their evolutionequations are necessarily PDE’s (partial differential equa-tions). However, under a “3+1” covariant decomposition ∗ Electronic address: [email protected] based on the 4–velocity field [4], the evolution equationsfor LTB dust solutions become first order autonomousPDE’s containing only time derivatives [3]. We provein this case how the spacelike gradient equations (con-straints) are automatically satisfied for all times onceinitial conditions are selected so that these constraintshold in an arbitrary initial hypersurface T i of constanttime. Rigorously (see Appendix B), the solution curves ofthese evolution equations are then equivalent to a properdynamical system built with ODE’s, under the strong re-striction of complying with the special set of initial con-ditions that is compatible with the spacelike constraintsof the PDE system.In order to describe LTB dust solutions as a dynamicalsystem we re–write the first order evolution equationsof the 3+1 decomposition [4] in terms of initial valuevariables defined in reference [2]. These variables areintroduced and generalized in section III and in sectionsIV and V we show how they lead in a very natural way toa system of three autonomous evolution equations thatis similar to those obtained with “expansion normalized”variables for other known space–times (FLRW, Bianchimodels, Kantowsky–Sachs, see [3]).Initial conditions (section VI) are selected from two“primitive” initial value functions (density, scalar 3–curvature) defined along an initial and regular hyper-surface, T i marked by constant comoving time. Thesefunctions determine the initial values for the phase spacevariables in such a way that the two spacelike constraintscharacterizing the 3+1 evolution equations for the LTBsolutions are satisfied at the T i and at subsequent T . Athird initial function determines the number of Regular a r X i v : . [ g r- q c ] N ov Symmetry Centers (RSC) at the T i (see Appendix A),and thus it determines the topological class of the T ,leading to hypersurfaces with “open” topology, (home-omorphic to R , one RSC); “closed” topology (homeo-morphic to S , two RSC) “wormhole” (homeomorphic to S × R or to S × S , zero RSC).The resulting 3–dimensional phase space is discussedin section VII, with its critical points, invariant subsetsand particular solutions given in section VIII. This spaceis parametrized by three functions: a dimensionless shearscalar, S , a density average parameter, (cid:104) Ω (cid:105) , and a den-sity 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 variables are eitherbounded or can be trivially compactified for all configu-rations. These functions have a straightforward physicalinterpretation: the average (cid:104) Ω (cid:105) behaves exactly as thestandard Ω in a FLRW dust spacetime: its initial val-ues (cid:104) Ω i (cid:105) determine the dynamical evolution of each so-lution curve: re–collapse (if (cid:104) Ω i (cid:105) − >
0) or perpetualexpansion (if (cid:104) Ω i (cid:105) − ≤ (cid:104) Ω (cid:105) defines invari-ant subspaces (cid:104) Ω (cid:105) = 0, (cid:104) Ω (cid:105) = 1, 0 < (cid:104) Ω (cid:105) < (cid:104) Ω (cid:105) >
1. The density contrast function ∆ ( m ) provides anon–local characterization of the type of inhomogeneityalong the T (rest frames of fundamental observers): a“clump” if − < ∆ ( m ) ≤ ( m ) ≥ ( m ) = 0 at the RSC. This variable also defines aninvariant subset: ∆ ( m ) = −
1, corresponding to vacuumSchwarzschild–Kruskal solutions. Another invariant sub-space is given by the line S = ∆ ( m ) = 0 marking thehomogeneous and isotropic FLRW subcase, containingas subcases the Einstein–de Sitter and Minkowski–Milneuniverses.We examine in detail in section IX the phase spaceevolution of several representative dust configurations:low density perpetually expanding (“hyperbolic” dynam-ics) with (cid:104) Ω i (cid:105) − ≤
0, as well as a case of void forma-tion scenario in which an initial density clump becomesa void (a numeric realization of configurations proposedin [10]). We also examine configurations with high den-sity re–collapsing (“elliptic”) dynamics with open, closedand wormhole topologies ( (cid:104) Ω i (cid:105) − > (cid:104) Ω i (cid:105) − > (cid:104) Ω i (cid:105)− ≤
0. We also examine the vacuum subcase, whichis the Schwarzschild–Kruskal spacetime given in comov-ing non–static coordinates made up by radial timelikegeodesics.All solution curves (except the RSC) for all configura-tions start their evolution at a past attractor (source) at (cid:104) Ω (cid:105) = 1 , S = 1 / , ∆ ( m ) = − (cid:104) Ω i (cid:105)− (cid:104) Ω (cid:105) . Curves with (cid:104) Ω i (cid:105)− (cid:104) Ω (cid:105) = 1 andterminate in a future attractor (sink) associated with the Einstein–de Sitter zero spacial curvature FLRW universe.Curves with (cid:104) Ω i (cid:105) − < (cid:104) Ω (cid:105) = S = 0, while curves with (cid:104) Ω i (cid:105) − > (cid:104) Ω (cid:105) → ∞ and then terminate at the same criticalpoint from which they started, though this point is nowa sink or a future attractor.The RSC’s of the configurations evolve separately fromthe rest of the curves along a line with ∆ ( m ) = S = 0,starting always from a source at (cid:104) Ω (cid:105) = 1 and going to-wards a sink at (cid:104) Ω (cid:105) = 0 or towards (cid:104) Ω (cid:105) → ∞ , respec-tively, for expanding and re–collapsing configurations. Inall cases some of the curves pass near a saddle close tothe RSC associated with homogeneity ( S = ∆ ( m ) = 0),while in the structure formation scenario some curvesalso approach another saddle point that splits collapsingcurves (cid:104) Ω i (cid:105) > (cid:104) Ω i (cid:105) ≤
1. The solu-tion curves of the vacuum case are confined to the plane∆ ( m ) = −
1, evolving towards (cid:104) Ω (cid:105) = 0 or (cid:104) Ω (cid:105) → ∞ , de-pending on whether the radial geodesics are bound ornot.We summarize and discuss the results from the articlein section X. We also provide three appendices dealingwith important issues: RSC’s and geometric propertiesof hypersurfaces T (Appendix A), treatment of evolu-tion equations given as PDE as a restricted dynamicalsystem (Appendix B) and analytic solutions in terms ofthe variables used in the article (Appendix C). II. LTB DUST MODELS IN THEIR ORIGINALVARIABLES.
The Lemaˆıtre–Tolman–Bondi (LTB) metric [1, 2, 5, 6]is the spherically symmetric line element ds = − c dt + Y (cid:48) − K dr + Y ( dθ + sin θdφ ) , (1)where Y = Y ( t, r ), Y (cid:48) = ∂Y /∂r and K = K ( r ). Themomentum-energy tensor usually associated with (1) isthat of a dust source: T ab = ρc u a u b , (2)where ρ = ρ ( t, r ) is the rest–mass density (in gm cm − ).For the metric (1) and source (2) with a comoving 4–velocity u a = δ a , Einstein’s field equations G a b = κT a b ,with κ = 8 πG/c , yield:˙ Y = 2 MY − K, (3)2 M (cid:48) = κρc Y Y (cid:48) , (4)where M = M ( r ) has units of length and ˙ Y = ∂Y /∂x .The only nonzero kinematic parameters are the ex-pansion scalar Θ = u a ; a and the shear tensor: σ ab = ∇ ( a ; b ) u − (Θ / h ab , given byΘ = 2 ˙ YY + ˙ Y (cid:48) Y (cid:48) , (5) σ a b = diag [0 , − , Σ , Σ] , Σ = 13 (cid:34) ˙ YY − ˙ Y (cid:48) Y (cid:48) (cid:35) . (6)Other important quantities are the “Electric” Weyl ten-sor: E ab = C abcd u c u d and the 3–dimensional Ricci ten-sor, R , of the hypersurfaces orthogonal to u a (markedby constant values of x = ct ): E a b = diag [0 , − E , E , E ] , E = M (cid:48) Y Y (cid:48) − MY = κ ρc − MY . (7)2( KY ) (cid:48) = R Y Y (cid:48) . (8)The usual treatment of LTB dust solutions consists insolving analytically the evolution equation (3). The so-lutions are usually classified in terms of the sign of K ,which determines the type of evolution of dust layers:perpetually expanding/collapsing “parabolic” ( K = 0)and “hyperbolic” ( K <
K >
III. VOLUME AVERAGES AND SCALINGLAWS
The form of ρ and R in terms of 2 M (cid:48) and ( KY ) (cid:48) in (4) and (8) suggests considering the following volumeaverages along an arbitrary hypersurface T i orthogonalto u a and marked by constant x i = ct i (cid:104) A i (cid:105) ≡ (cid:82) A i d V i (cid:82) d V i , d V i = Y i Y (cid:48) i sin θ dr dθ dφ, (9)where the subindex i denotes evaluation at t = t i and wehave assumed that T i is fully regular in the integrationrange. Applying (9) to (4) and (8) we get κc (cid:104) ρ (cid:105) = κc (cid:82) ρ Y Y (cid:48) dr (cid:82) Y Y (cid:48) dr = 2 MY , (10)16 (cid:104) R(cid:105) = 16 (cid:82) R Y Y (cid:48) dr (cid:82) Y Y (cid:48) dr = KY , (11)where we have dropped the subindex i and have takenthe lower bound of the integration range of r to be deter-mined by a suitable boundary condition on ρ i and R i ,for example a RSC. See Appendix A and section IX–E. The averages (cid:104) ρ (cid:105) and (cid:104) R(cid:105) are non–local quantities de-pending on the form of ρ and R along the integrationrange. We define the following “contrast functions” com-paring these quantities with their local counterparts:∆ ( m ) ≡ ρ − (cid:104) ρ (cid:105)(cid:104) ρ (cid:105) , ⇒ ρ = (cid:104) ρ (cid:105) [1 + ∆ ( m ) ] , (12)∆ ( k ) ≡ R − (cid:104) R(cid:105)(cid:104) R(cid:105) , ⇒ R = (cid:104) R(cid:105) [1 + ∆ ( k ) ] , (13)The interpretation of these these contrast functions [2]follows by using the definitions (4), (10) and (12) andintegrating by parts along the T , leading to∆ ( m ) = κc M (cid:90) ρ (cid:48) Y dr, (14)where this integral is evaluated from a RSC. Since M ≥ Y ≥ ρ ≥
0, the only quantity thatcan change sign is ρ (cid:48) , therefore, looking at the radialvariation of a density profile from the symmetry centeralong an arbitrary T we have∆ ( m ) (cid:26) ≤ ⇔ ρ (cid:48) ≤ , ρ ≤ (cid:104) ρ (cid:105) , density clump ≥ ⇔ ρ (cid:48) ≥ , ρ ≥ (cid:104) ρ (cid:105) , density void (15)Likewise, for the contrast function ∆ ( k ) , using the def-initions (8), (11) and (13) and integrating by parts alongthe T we obtain an analogous expression to (14):∆ ( k ) = 16 K (cid:90) R (cid:48) Y dr, (16)However, while a negative ρ is not physically interesting,there is no physical objection for R being positive ornegative, or changing sign in the allowed range of r . If R > r , we have the sameresults as with ∆ ( m ) in (15), but if R <
0, then curva-ture clumps or voids are defined by the opposite signs [2].We define now the following “scale factors” related tothe metric functions: (cid:96) ≡ YY i , (17)Γ ≡ Y (cid:48) /YY (cid:48) i /Y i = 1 + (cid:96) (cid:48) /(cid:96)Y (cid:48) i /Y i , (18)relating Y and Y (cid:48) evaluated at an arbitrary and a fiducial(or “initial”) hypersurface T i . Since (10) and (11) arevalid at any T , we obtain the scaling laws (cid:104) ρ (cid:105) = (cid:104) ρ i (cid:105) (cid:96) , (19) (cid:104) R(cid:105) = (cid:104) R i (cid:105) (cid:96) . (20)while from (4), (8) and (17)–(18) we get ρ = ρ i (cid:96) Γ , (21) R = 1 (cid:96) Γ (cid:20) R i + 13 (cid:104) R i (cid:105) (1 − Γ) (cid:21) , (22)Comparing (12)–(13) with (19)–(20) and (21)–(22) weget the scaling laws for ∆ ( m ) and ∆ ( k ) ( m ) = 1 + ∆ ( m ) i Γ , (23)23 + ∆ ( k ) = 2 / ( k ) i Γ , (24)The quantities Θ, Σ and E given in (5), (6) and (7)take the formsΘ = 3 ˙ (cid:96)(cid:96) + ˙ΓΓ , (25)Σ = − ˙Γ3Γ , (26) E = κc ρ − (cid:104) ρ (cid:105) ] = κc (cid:104) ρ (cid:105) ∆ ( m ) , (27) the evolution equation (3) becomes˙ (cid:96) (cid:96) = κc (cid:104) ρ (cid:105) − (cid:104) R(cid:105) , (28)while the LTB metric (1) takes the Friedmanian form ds = − c dt + (cid:96) (cid:20) Γ Y (cid:48) i dr − (cid:104) R i (cid:105) Y i + Y i (cid:0) dθ + sin θ dφ (cid:1)(cid:21) , (29)It is also helpful to express the spacial ( i.e. radial ) gra-dients of (cid:104) ρ (cid:105) and (cid:104) R(cid:105) in terms of the contrast functionsand Γ. With the help of (4), (8), (10), (11), (12) and (13)we obtain the following useful relations: (cid:104) ρ (cid:105) (cid:48) (cid:104) ρ (cid:105) = 3 Y (cid:48) Y ∆ ( m ) = 3 Y (cid:48) i Y i ∆ ( m ) Γ , (30) (cid:104) R(cid:105) (cid:48) (cid:104) R(cid:105) = 3 Y (cid:48) Y ∆ ( k ) = 3 Y (cid:48) Y ∆ ( k ) Γ , (31)Notice that these equations are strictly valid along allhypersurfaces T . IV. 3+1 DECOMPOSITION.
From a covariant 3+1 decomposition of spacetimebased on the 4–velocity field u a presented in [3, 4], thefield and energy balance equations for a dust model char-acterized by (1) and (2) are equivalent to the followingset of first order evolution equations:˙Θ3 = − (cid:18) Θ3 (cid:19) − σ ab σ ab − κc ρ = 0 , (32a)˙ σ (cid:104) ab (cid:105) = −
23 Θ σ ab − σ (cid:104) ac σ b (cid:105) c − E ab = 0 , (32b)˙ E (cid:104) ab (cid:105) = − κc ρ σ ab − Θ E ab + 3 σ (cid:104) ac E b (cid:105) c , (32c)˙ ρ = − ρ Θ , (32d) together with the constraints:˜ ∇ b σ b a −
23 ˜ ∇ a Θ = 0 , (33a)˜ ∇ b E b a − κc ∇ a ρ = 0 , (33b)where ˙ A ab ≡ u c ∇ c A ab is the convective derivative alongthe 4–velocity, ˜ ∇ a A bc ≡ h ad ∇ d A bc is the spacelike gra-dient, tangent to the T hypersurfaces and orthogonalto u a , and A (cid:104) ab (cid:105) ≡ A ( ab ) − (1 / A cc h ab is the spacelike,symmetric, trace–free part of a tensor A ab .Equations (32a)–(32d) are 4 tensorial equations for 4tensorial quantities Θ , ρ, σ ab , E ab , which together withthe constraints (33a)–(33b) form a completely deter-mined system of partial differential equations (equiva-lent to Einstein’s field equations). Bearing in mind that u a = δ a and considering the forms of the trace–free ten-sors σ ab and E ab in (6) and (7), the system (32) reducesto the following set of scalar equations:˙Θ3 = − (cid:18) Θ3 (cid:19) − − κc ρ = 0 , , (34a)˙Σ = −
23 Θ Σ + Σ − E , (34b)˙ E = −E − Σ (cid:104) κ ρc + 3 E (cid:105) , (34c)˙ ρ = − ρ Θ , (34d)and (33a)–(33b) becomeΣ (cid:48) + Θ (cid:48) Y (cid:48) Y = 0 , (35a) E (cid:48) + κρ (cid:48) c E Y (cid:48) Y = 0 . (35b)It is straightforward to express the system (34) and theconstraints (35) in therms of the new variables introducedin the previous sections. Considering (12), (23) and (25)–(27) and and defining H ≡ ˙ (cid:96)(cid:96) = Θ3 + Σ , (36)the system (34) becomes the following set of equivalentevolution equations:˙ H = − H − κ (cid:104) ρ (cid:105) , (37a)˙Σ = 3Σ − H Σ + κ (cid:104) ρ (cid:105) ∆ ( m ) , (37b) (cid:104) ρ (cid:105) ˙ = − (cid:104) ρ (cid:105) H, (37c)˙∆ ( m ) = 3[1 + ∆ ( m ) ] Σ , (37d)while the constraints (35b) and (35a) become respectively (cid:104) ρ (cid:105) (cid:48) − (cid:104) ρ (cid:105) ∆ ( m ) Γ Y (cid:48) i Y i = 0 , (38a) H (cid:48) − Y (cid:48) Y = H (cid:48) − Y (cid:48) i Y i = 0 . (38b)The constraint (33b) leads to (38a) which is identicalto (30), so this constraint is automatically satisfied at all t (all T ) by a system of evolution equations based onthe variables (cid:104) ρ (cid:105) and ∆ ( m ) .Bearing in mind (17), (18) and (26), the constraint(33a) in its form (38b) can be written as the integrabilitycondition for (cid:96) :Γ (cid:32) ˙ (cid:96)(cid:96) (cid:33) (cid:48) − ˙Γ (cid:18) Y (cid:48) i Y i + (cid:96) (cid:48) (cid:96) (cid:19) = Γ (cid:34)(cid:32) ˙ (cid:96)(cid:96) (cid:33) (cid:48) − (cid:18) (cid:96) (cid:48) (cid:96) (cid:19) ˙ (cid:35) = 0 . (39)Further, since (39) holds and (36) implies H = ˙ (cid:96)(cid:96) = (cid:20) κc (cid:104) ρ (cid:105) − (cid:104) R(cid:105) (cid:21) / , (40)the constraint (38b) implies (with the help of (23), (24),(30), and (31)):Σ = − κc (cid:104) ρ (cid:105) ∆ ( m ) − (cid:104) R(cid:105) ∆ ( k ) (cid:2) κc (cid:104) ρ (cid:105) − (cid:104) R(cid:105) (cid:3) / . (41)If the system (37) is self–consistent, then (41) must holdfor every T and so it must comply with (37b). In-serting (41) into (37b) and eliminating time derivativesof (cid:104) ρ (cid:105) , ∆ ( m ) , H, (cid:104) R(cid:105) and ∆ ( k ) with the help of (37a),(37c), (37d) and the scaling laws (23)–(24), it is straight-forward to see that (37b) is identically satisfied at every T . V. DIMENSIONLESS VARIABLES: ADYNAMICAL SYSTEM
In order to transform system (37a)–(37d) into a dy-namical system, we need to recast H, Σ and (cid:104) ρ (cid:105) in termsof dimensionless “expansion normalized” variables. [3].However, instead of using the expansion parameter Θ in(25), it turns out to be easier to use H defined in (36).This leads to: H ≡ HH , (42) (cid:104) Ω (cid:105) ≡ κc (cid:104) ρ (cid:105) H , (43) S ≡ Σ H , (44)where H is a characteristic length scale. Notice that wecan also define a dimensionless Ω parameter for the localdensity, which would be related to that of (43) by:Ω = κc ρ H = (cid:104) Ω (cid:105) [1 + ∆ ( m ) ] . (45)We also need to make the “dot” derivative ∂/∂ct a di-mensionless operator, so we define: ∂∂τ ≡ H ∂c∂t . (46)Inserting (42)–(46) into (37a)–(37d) we obtain the fol-lowing dimensionless system: H ,τ = − (cid:2) (cid:104) Ω (cid:105) (cid:3) H , (47a) S ,τ = H (cid:2) S (3 S −
1) + (cid:0) ∆ ( m ) + S (cid:1) (cid:104) Ω (cid:105) (cid:3) , (47b) (cid:104) Ω (cid:105) ,τ = H (cid:104) Ω (cid:105) [ (cid:104) Ω (cid:105) − , (47c)∆ ( m ) ,τ = 3 H S (cid:104) ( m ) (cid:105) . (47d)where ,τ = ∂/∂τ .The form of these equations suggests that further sim-plification follows if the three equations (47b)–(47d) canbe decoupled from (47a). Defining the following coordi-nate transformation: τ = τ ( ξ, ¯ r ) , r = ¯ r (48)so that for any function A = A ( τ, r ) we have A ( τ, r ) = A ( ξ ( τ, ¯ r ) , ¯ r ) and all partial time drivatives in (47) can beexpressed as (cid:20) ∂A∂τ (cid:21) r = ∂A∂ξ (cid:20) ∂ξ∂τ (cid:21) r = H ∂A∂ξ , (49)where ξ is selected so that: (cid:20) ∂ξ∂τ (cid:21) r = H , ⇒ ξ = ln (cid:96), (50)The introduction of the new variable ξ defined by (49)and (50) removes the dependence on H of the evolutionequations (47b), (47c) and (47d): ∂S∂ξ = S (3 S −
1) + 12 (cid:16) ∆ ( m ) + S (cid:17) (cid:104) Ω (cid:105) , (51a) ∂ ∆ ( m ) ∂ξ = 3 S (cid:104) ( m ) (cid:105) , (51b) ∂ (cid:104) Ω (cid:105) ∂ξ = (cid:104) Ω (cid:105) [ (cid:104) Ω (cid:105) − , (51c)while (47a) becomes: ∂∂ξ (ln H ) = − (cid:20) (cid:104) Ω (cid:105) (cid:21) , (52)The system (51) is formally a dynamical system con-structed with “expansion normalized” variables [3]. How-ever, it is still a system of PDE’s, even if it only containsderivatives with respect to ξ and r is basically a param-eter. As we show in Appendix B, this system of PDE’sequivalent authonomous ODE system with restricted ini-tial conditions (restricted in order to fulfill the spacelikeconstraints). VI. INITIAL CONDITIONS
Bearing in mind (4), (11), (12)–(13) and (40)–(41),intial conditions for the system (51) can be constructedthrough the following steps:1. Choose the topological class of the initial hyper-surface T i in terms of the number of RSC (seeAppendix A). A convenient choice is: Y i = H − f ( r ) , (53)where f ( r ) is a (at least a C ) function whose zeroscorrespond to RSC’s and H is the same character-istic length scale as in (42) and (46). Strictly speak-ing, the choice of f is a choice of radial coordinate(see Appendix A).2. Construct dimensionless quantities out of ρ i and R i : m i ≡ κc ρ i H , k i = R i H , (54)3. Obtain the orbit volume averages and contrastfunctions (cid:104) m i (cid:105) = κc H (cid:104) ρ i (cid:105) , (55a) (cid:104) k i (cid:105) = 16 H (cid:104) R i (cid:105) , (55b)∆ ( m ) i = m i (cid:104) m i (cid:105) − , (55c)∆ ( k ) i = k i (cid:104) k i (cid:105) − , (55d) 4. The initial forms of the remaining variables arethen: H i = [ (cid:104) m i (cid:105) − (cid:104) k i (cid:105) ] / , (56a) (cid:104) Ω i (cid:105) = (cid:104) m i (cid:105)(cid:104) m i (cid:105) − (cid:104) k i (cid:105) , (56b) S i = − (cid:104) m i (cid:105) ∆ ( m ) i − (cid:104) k i (cid:105) ∆ ( k ) i (cid:104) m i (cid:105) − (cid:104) k i (cid:105) , (56c)Just as we proceeded with (37), the fulfillment of con-straints (38a) and (38b) imply that the functional formsof H , (cid:104) Ω (cid:105) and S are the generalization of (56) for all times,that is: H = [ (cid:104) m (cid:105) − (cid:104) k (cid:105) ] / , (57a) (cid:104) Ω (cid:105) = (cid:104) m (cid:105)(cid:104) m (cid:105) − (cid:104) k (cid:105) , (57b) S = − (cid:104) m (cid:105) ∆ ( m ) − (cid:104) k (cid:105) ∆ ( k ) (cid:104) m (cid:105) − (cid:104) k (cid:105) , (57c)where (cid:104) m (cid:105) = κc (cid:104) ρ (cid:105) H = (cid:104) m i (cid:105) (cid:96) , (58a) (cid:104) k (cid:105) = (cid:104) R(cid:105) H = (cid:104) k i (cid:105) (cid:96) . (58b)Inserting (57) and (58) into (51) and eliminating deriva-tives of (cid:104) m (cid:105) and (cid:104) k (cid:105) by means of (24), (26), (37a), (37c)and (37d), we can see that (57) and (58) fully solve (47).Thus (47) propagates the initial data given by (53)–(56c)with the constraints (38a) and (38b) satisfied for all T .The standard approach to LTB dust solutions is basedon solving the evolution equation (3). Therefore, it is il-lustrative to re–write such equation in terms of the initialvalue variables that we have defined here: (cid:96) ,τ = H i (cid:20) (cid:104) Ω i (cid:105) (cid:96) − ( (cid:104) Ω i (cid:105) − (cid:21) . (59)This form of (3) shows how the sign of (cid:104) Ω i (cid:105)− K does it in (3), leading to “parabolic” ( (cid:104) Ω i (cid:105) − (cid:104) Ω i (cid:105) − <
0) and “elliptic” ( (cid:104) Ω i (cid:105) − > K = 0 , K < , K > H ,τ = −H − (cid:104) m (cid:105) , (60a) s ,τ = 3 s − H s + (cid:104) m (cid:105) ( m ) , (60b) (cid:104) m (cid:105) ,τ = − (cid:104) m (cid:105) H , (60c)∆ ( m ) ,τ = 3(1 + ∆ ( m ) ) s, (60d)where τ , H and (cid:104) m (cid:105) have been defined in (42), (46), (58)and s = Σ /H . This system will be specially useful forlooking at the collapsing stage of re–collapsing configu-rations, a feature that cannot be studied by (51) becausemaximal expansion takes place as H →
0, and so (cid:104) Ω (cid:105) di-verges and the solution curves cannot be extended furtherto study their collapsing stage. Also, (60) yields directlyquantities like (cid:104) m (cid:105) and H , while the metric functions (cid:96) and Γ in (29) follow from the variables in (60) as: (cid:96) = (cid:20) (cid:104) m i (cid:105)(cid:104) m (cid:105) (cid:21) / , Γ = 1 + ∆ ( m ) i ( m ) , (61)and the so–called “curvature radius” Y = √ g θθ can beeasily found as Y = Y i (cid:96) . VII. SOME QUALITATIVE AND ANALYTICRESULTS.
The solution curves [ S ( ξ, r ) , ∆ ( m ) ( ξ, r ) , (cid:104) Ω (cid:105) ( ξ, r )] of(51) evolve in a phase space which is a region of R parametrized by the coordinates [ S, ∆ ( m ) , (cid:104) Ω (cid:105) ]. The evo-lution variable is ξ , while each curve represents a funda-mental comoving observer labelled by a constant value r . It is highly desirable that the phase space coordinatesremain bounded as these curves evolve along their maxi-mal range of ξ for all r . We discuss in this section severaluseful analytic results and the relation between the phasespace variables and their initial values.Considering that ξ = ln (cid:96) and (cid:104) Ω i (cid:105) − (cid:104) k i (cid:105)(cid:104) m i (cid:105) − (cid:104) k i (cid:105) , (62a) (cid:104) Ω (cid:105) − (cid:104) k (cid:105)(cid:104) m (cid:105) − (cid:104) k (cid:105) = (cid:104) k i (cid:105) e ξ (cid:104) m i (cid:105) − (cid:104) k i (cid:105) e ξ , (62b)the functions H , (cid:104) Ω (cid:105) and S in (57)–(58) can be given as (cid:104) Ω (cid:105) = (cid:104) Ω i (cid:105)(cid:104) Ω i (cid:105) − [ (cid:104) Ω i (cid:105) − ξ , (63a) H = H i (cid:2) (cid:104) Ω i (cid:105) e − ξ − ( (cid:104) Ω i (cid:105) −
1) e − ξ (cid:3) , (63b) S = − (cid:104) Ω i (cid:105) ∆ ( m ) − [ (cid:104) Ω i (cid:105) −
1] ∆ ( k ) e ξ (cid:104) Ω i (cid:105) − [ (cid:104) Ω i (cid:105) −
1] e ξ ] . (63c)where we are emphasizing that these are now functionsof ( ξ, r ). Notice that even if, in general, the surfaces ofconstant ξ do not coincide with the T (surfaces of con-stant t or τ ), the initial hypersurface ξ = 0 does coincidewith the initial T i . A. Constraints on (cid:104) Ω (cid:105) From(57), we have the following important constraint: (cid:104) Ω (cid:105) − (cid:104) k (cid:105)H = (cid:104) k i (cid:105) (cid:96) H = (cid:104) Ω i (cid:105) − (cid:96) H . (64) Since (cid:104) Ω (cid:105) ≥
0, this constraint (together with (63a)) hasimportant qualitative consequences: it defines invariantsubsets associated with (cid:104) Ω (cid:105) . If initial conditions for arange of r are selected so that (cid:104) Ω i (cid:105) = 0 , < (cid:104) Ω i (cid:105) < , (cid:104) Ω i (cid:105) = 1, or (cid:104) Ω i (cid:105) >
1, then all solution curves forsuch range respectively comply with (cid:104) Ω (cid:105) = 0 , < (cid:104) Ω (cid:105) < , (cid:104) Ω (cid:105) = 1, or (cid:104) Ω (cid:105) > ξ .Notice from (63) that, irrespectively of the sign of (cid:104) Ω i ( r ) (cid:105)−
1, we have (cid:104) Ω (cid:105) → ξ → −∞ (or, equivlently,as (cid:96) → (cid:104) Ω i (cid:105) −
1, all solution curves of(51) will approach the plane (cid:104) Ω (cid:105) = 1 near the initial bigbang singularity (in the asymptotic range: ξ → −∞ or (cid:96) → (cid:104) Ω (cid:105) ≥
0, (64) implies that all solution curveswith initial conditions 0 ≤ (cid:104) Ω i (cid:105) ≤ ≤ (cid:104) Ω (cid:105) ≤ ξ , while (63) implies that the such curveswill evolve from the plane (cid:104) Ω (cid:105) = 1 in ξ → −∞ ( (cid:96) → (cid:104) Ω (cid:105) = 0 in the asymptotic range ξ →∞ (or (cid:96) → ∞ ).Equations (63) place no restriction in the range of ξ of solution curves when (cid:104) Ω i (cid:105) ≤
1, but curves with initialconditions given by (cid:104) Ω i (cid:105) > ξ will have a bounded range of maximal extendibility givenby (cid:104) Ω (cid:105) → ∞ , as ξ → ln (cid:104) Ω i (cid:105)(cid:104) Ω i (cid:105) − . (65)Such solution curves start from the plane (cid:104) Ω (cid:105) = 1 as ξ → −∞ ( (cid:96) →
0) and diverge for the limiting value of ξ = ln (cid:96) given by (65). However, this is the same finitevalue of ξ that makes H in (63b) vanish, so the blowingup (cid:104) Ω (cid:105) → ∞ of all solution curves in this case simplymarks the values of (cid:96) associated with the “maximal” ex-pansion of dust layers (maximal value of (cid:96) ) in the “turnaround” before the collapsing stage. Since (cid:104) Ω (cid:105) blows upthe collapsing stage cannot be described by (51) (thoughit can be described by (60)). B. Constraints on S and ∆ ( m ) In order to appreciate the relation between the coordi-nates of the phase space, [ S, ∆ ( m ) , (cid:104) Ω (cid:105) ], it is very usefulto express S given by (57c) as: S = − (cid:110) (cid:104) Ω (cid:105) ∆ ( m ) − [ (cid:104) Ω (cid:105) −
1] ∆ ( k ) (cid:111) , (66)where, from the scaling laws (23) and (24), we have∆ ( m ) = 1 + ∆ ( m ) i Γ − , (67a)∆ ( k ) = 2 / ( k ) i Γ − . (67b)Thus, for finite ∆ ( m ) i and ∆ ( k ) i , it is evident that a suf-ficient condition for ∆ ( m ) and ∆ ( k ) to remain boundedfor all ξ is that initial conditions are selected so that ashell crossing singularity does not emerge (see AppendixA and references [2, 5, 6, 7]). That is, we must have forall solution curves marked by ( ξ, r )Γ( ξ, r ) > , no shell crossing singularity, (68)where Γ has been defined in (18). Notice, from (63c), (66)and (67), that even if (68) holds, S (like (cid:104) Ω (cid:105) ) divergesas H → (cid:104) Ω i (cid:105) >
1, though it remains bounded if (cid:104) Ω i (cid:105) ≤
1. In order toovercome the blowing up of S and (cid:104) Ω (cid:105) in the cases when (cid:104) Ω i (cid:105) > < (cid:104) Ω i (cid:105) ≤ S, ∆ ( m ) , (cid:104) Ω (cid:105) ], of all solutioncurves remain bounded and restricted to a finite regionof R .Initial conditions that guarantee the fulfillment of (68)were derived in terms of original variables in reference [7](see also [5, 6]), in terms of initial value variables dis-cussed here in [2] (see Appendix A). VIII. CRITICAL POINTS AND PARTICULARCASES
In the general case in which none of the variables[ S, ∆ ( m ) , (cid:104) Ω (cid:105) ] is restricted to take any special constantvalue associated with particular cases, the coordinatesand nature of critical points of (51) are: • (cid:104) Ω (cid:105) = 1: C S = 1 / , ∆ ( m ) = − , source , (69a) C S = − / , ∆ ( m ) = − , saddle , (69b) C S = 0 , ∆ ( m ) = 0 , saddle , (69c) • (cid:104) Ω (cid:105) = 0: C S = 1 / , ∆ ( m ) = − , saddle , (70a) C S = 0 , ∆ ( m ) arbitrary , sink , (70b) C S = 0 , ∆ ( m ) = 0 , sink , (70c)The nature of some of the critical points of (51) changeswhen [ S, ∆ ( m ) , (cid:104) Ω (cid:105) ] take certain specific particular valuesthat correspond to various space–times that are particu-lar cases of dust LTB solutions. We examine these par-ticular cases and their critical points in the remaining ofthis section. The phase space and the location of criticalpoints and location of all particular cases is depicted infigure 1, FIG. 1:
Phase space
The figure shows all critical points andthe invariant subspaces given by the surfaces (cid:104) Ω (cid:105) = 1 and∆ ( m ) = − (cid:104) Ω (cid:105) = 1 (thecritical point C ) and the Minkowski–Milne space–time (M) (cid:104) Ω (cid:105) = 0 (the critical point C ). A. Parabolic (or “marginally bound”) evolution
For initial conditions (cid:104) Ω i (cid:105) = 1 equations (63) imply (cid:104) Ω (cid:105) = 1 for all ξ . This is an invariant set associatedwith the so–called “parabolic” solutions of (59) and cor-responding to (cid:104) k (cid:105) = (cid:104) R(cid:105) = K = 0 (see [5, 6]). A fullclosed analytic solution (compatible with (57)) is givenby: ∆ ( m ) ( ξ, r ) = ∆ ( m ) i [1 + ∆ ( m ) i ] e ξ/ − ∆ ( m ) i , (71a) H ( ξ, r ) = H i e − ξ/ , (71b) S ( ξ, r ) = −
12 ∆ ( m ) ( ξ, r ) . (71c)so that solution curves remain for all ξ in the line S = − ∆ ( m ) / (cid:104) Ω (cid:105) = 1 (see figure 1). Wehave all variables fully determined. It is straightforwardto verify that (71) are fully compatible and equivalent toanalytic solutions given in terms of (cid:96) ( ct, r ) (see [2] andAppendix C).Critical points on this case are: C : source , C : sink , (72)The phase space evolution goes along the line [ S, − S, C to C (respectively, past and future global at-tractors). Since the critical point C is characterized byconditions of homogeneity ∆ ( m ) = S = 0, it correspondsto the Einstein–de Sitter universe (dust FLRW modelwith zero spacial curvature). B. FLRW dust models
If ∆ ( m ) = S = 0 with unrestricted (cid:104) Ω (cid:105) , we have aninvariant set associated with the particular homogeneousand isotropic sub–case of dust FLRW models. The criti-cal points are in this case: C : source , C : sink , (73)Depending on the sign of (cid:104) Ω i (cid:105) −
1, the phase space evolu-tion of these models will take place in the line [0 , , (cid:104) Ω (cid:105) ].If (cid:104) Ω i (cid:105) − < C and the sink C (past and future attractors). If (cid:104) Ω i (cid:105) − >
0, it starts at C and goes upwards to (cid:104) Ω (cid:105) → ∞ towards maximal expansion. In the collaps-ing stage, the point C acts as a sink. If (cid:104) Ω (cid:105) = 1, wehave the Einstein–de Sitter universe and the evolutionis constrained to the point C with coordinates [0 , , C , whose coordinates are[0 , , C. RSC’s
The solution curve corresponding to a RSC is also char-acterized by ∆ ( m ) = S = 0 with varying (cid:104) Ω (cid:105) , thus thephase space evolution of that particular curve will alsobe along the line [0 , , (cid:104) Ω (cid:105) ], just like the FLRW sub–case. This is not surprising because a RSC is a privilegedisotropic observer.The fact that the RSC’s have a separate evolutionfrom the rest of the solution curves prevents the criti-cal point C (a past and future attractor associated withthe initial and collapsing singularities) to be global in allcases. However, C is a global attractor for all “off cen-ter” curves in configurations with RSC and for all curvesin configurations lacking RSC (wormhole topology andvacuum limit). D. Vacuum case
It is known that the LTB metric (1) can describe theSchwarzchild-Kruskal space–time in terms of the worldlines of its radial geodesic test observers with 4-velocity u a = δ a . The equation of motion of these geodesics isidentical to (3) with M = M , so that M (cid:48) = 0 = ρ and the Ricci tensor vanishes, while the constant M can be identified with the “Schwarzschild mass” (see[5, 6, 8]). The function R is the scalar 3-curvature ofthe hypersurfaces of simultaneity of these geodesic con-gruences and the average scalar curvature is related totheir binding energy per unit rest mass of the test parti-cle E i ( r ) = − K/ − / (cid:104) R i (cid:105) Y i = ( − / H (cid:104) k i (cid:105) .This particular case can also be constructed by assum-ing a Dirac delta distribution for ρ i so that a nonzero (cid:104) ρ i (cid:105) follows from the integral definition (10). Also, without resorting to a Dirac delta, this vacuum case follows fromthe 3+1 formalism in equations (34) by setting ρ = 0, soa definition of a nonzero average density follows from theWeyl tensor component in (7) and (27) as: E = − κc (cid:104) ρ (cid:105) = − (cid:104) m (cid:105) H − M /Y i (cid:96) (74)Since m = ρ = 0 but (cid:104) m (cid:105) >
0, as long as Γ is finite (seesection X) we have the following sufficient condition tocharacterize the Schwarzchild-Kruskal vacuum:∆ ( m ) = − ξ, r ) (75)which defines an invariant set, since all solution curveswith initial condition ∆ ( m ) i = − ( m ) = −
1, parametrized by[ S, (cid:104) Ω (cid:105) ]. The phase space evolution and critical pointsdepends on the sign of (cid:104) Ω i (cid:105) , which depends in turn onthe binding energy of the radial geodesics.For zero binding energy (cid:104) k i (cid:105) = 0 = (cid:104) Ω i (cid:105) − C : source , C : sink , (76)and phase space evolution is the line [ S, − , C to the sink C .For geodesics with negative binding energy: (cid:104) k i (cid:105) > (cid:104) Ω i (cid:105) > C is now a saddle. The phase space evolutionof the curves begin at the global past attractor (source) C , approach the saddle C and go upwards towards di-verging (cid:104) Ω (cid:105) at the maximal expansion, terminating againat the global future attractor (sink) C .For geodesics with positive binding energy ( (cid:104) k i (cid:105) < < (cid:104) Ω i (cid:105) <
1) the critical points are: C : source , C : sink , C : saddle , C : saddle . (77)In this case the curves begin at the global past attrac-tor (source) C , terminate at the global future attractor(sink) C , approaching the saddles C and C . Thiscase is examined in subsection F of next section (see fig-ure 18).We examine the phase space evolution of solutioncurves of several general case dust LTB configurationsin the following section. IX. LTB DUST CONFIGURATIONSA. Expanding low density universe
A low density perpetually expanding dust configura-tion follows by selecting (cid:104) m i (cid:105) (0) < (cid:104) k i (cid:105) r defined in an initial hypersurface T i with anopen topology having a RSC. This leads to (cid:104) Ω i (cid:105) < m i = 0 .
91 + tan r ,k i = − tan r r ,Y i = H − tan r, (78)so that we have an initial density clump ( − < ∆ ( m ) i ≤
0) but a curvature void (0 ≤ ∆ ( k ) i < / (cid:104) Ω i (cid:105) < r , so that 0 < (cid:104) Ω (cid:105) < ξ . These conditionsalso comply with the no–shell–crossing conditions (93),so dust layers monotonously expand perpetually from aninitial big bang.As shown by figure 2 the solution curves start at thepast attractor or source C (big bang), evolve towardsthe saddle C and terminate at the future attractor givenby the line of sinks C in the plane (cid:104) Ω (cid:105) = 0. The evo-lution of the RSC is constrained to C for all time. Thefunction ∆ ( m ) remains negative for all ξ , thus we havedensity clumps at all T . FIG. 2:
Low density configuration.
Dust layers expandfrom the past attractor (source) C , denoting the big bangsingularity, approaches the saddle C and drops towards thefuture attractor given by the line of sinks C . Notice that(78) implies that the evolution of the RSC is constrained to C for all time. B. Low density void formation
An interesting variation of a perpetually expandingconfiguration, similar to the one presented in subsectionA and figure 2, is that where an initial density clump ( − < ∆ ( m ) i <
0) evolves for all curves into density voidswith ∆ ( m ) >
0. In this subsection we provide a numericexample of dust configurations discussed in reference [10](see also [11, 12]).An initial dust clump transforming into a void emergesby choosing the following initial value functions m i = m + m α tan r ,m = 0 . , m = 0 . , α = 2 . ,k i = k + k β tan r ,k = − . , k = − . , β = 0 . ,Y i = H − tan r, (79)which comply with the no–shell–crossing conditions (93).The evolution of this configuration is similar to the pre-vious one: dust layers also expand perpetually (“hy-perbolic” dynamics) from an initial big bang and thehypersurfaces T have an open topology with a RSC.However, the contrast function, ∆ ( m ) , now passes fromnegative to positive as ξ increases, indicating that ini-tial clumps transform into voids. This can be seen veryclearly if we plot (see figure (3)) the normalized den-sity profiles along the hypersurfaces T from the function ρ/ρ c = m/m c = (cid:104) m (cid:105) (1+∆ ( m ) ) /m c , where the subindex c denotes evaluation along the RSC marked by r = 0. Thedensity as function of τ and r can be readily found bysolving numerically the system (60) for initial conditions(79). This transition from clumps to voids is depicted byfigure 3. FIG. 3:
Low density voids: density profiles
The picturedepicts the normalized radial density profile ρ/ρ c for differenthypersurfaces T marked by constant τ . Notice how an initialclump at τ = 0 evolves into a voids profile for latter τ . This transition from clumps to voids can also be seenfrom the form of the contrast density function ∆ ( m ) given1in terms of ξ and r , displayed in figure 4 which showshow ∆ ( m ) is initially negative but becomes positive forall solution curves with r constant (except the RSC at r = 0) as ξ increases. FIG. 4:
Density contrast function
The picture depicts thefunction ∆ ( m ) passing for all 0 < r < π/ ( m ) = 0 at the RSC r = 0, and asymp-totically at r → π/ Y i → ∞ ). For solution curves with r close to the RSC, ∆ ( m ) remains close to zero but is neverthe-less positive. The evolution of this configuration in phase space isshown by figure 5. Solution curves start at the past at-tractor (source) C (big bang), approach the saddle C and evolve towards the future attractor (line of sinks) C in the plane (cid:104) Ω (cid:105) = 0. The RSC evolves from the source C towards the sink C . However, now the curves evolvein such a away that all curves (save the RSC) hit theline of sinks C at positive values of ∆ ( m ) (passage fromclumps to voids). C. Re–collapsing high density universes
Re–collapsing high density configurations follow by se-lecting (cid:104) m i (cid:105) (0) > (cid:104) k i (cid:105) > r , leading to (cid:104) Ω i (cid:105) ≥ r , thus (cid:104) Ω (cid:105) ≥ ξ . Dust lay-ers expand from an initial big bang, reach a maximalexpansion and then collapse (the so–called “elliptic” dy-namics [2, 5, 6]). These configurations are compatiblewith T i having either a closed topology with two RSC,an open topology with one RSC or a wormhole topologywithout RSC (see Appendix A). We examine the closedand open case separately below, and the wormhole casein section E.The case with spherical “closed” topology can be ob- FIG. 5:
Low density voids
The evolution is similar to thatof Fig 1, except that ∆ ( m ) becomes positive as all solutioncurves drop to the line of sinks C . tained with the following initial value functions: m i = m + m − m r , m = 3 . , m = 1 . k i = k + k − k r , k = 2 . , k = 0 . Y i = H − sin r, (80)which comply with (91) at the turning point r = π/ C (big bang) and evolve upwards as (cid:104) Ω (cid:105) → ∞ , which cor-responds to maximal expansion given by H → ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − C . For off center solution curves thiscritical point is a saddle.Since the collapsing stage cannot be described withsolution curves of system (51), we can examinethis stage with system (60) by plotting the curves[ S ( τ, r ) , ∆ ( m ) ( τ, r ) , (cid:104) Ω (cid:105) ( τ, r )], where S and (cid:104) Ω (cid:105) are de-fined by (44) and (43). As shown by figure 7, the curvescome downwards from infinite (cid:104) Ω (cid:105) and S to the criticalpoint C , which is now a future attractor or sink, associ-ated with a second curvature singularity (“big crunch”).As in the expanding state, curves near the RSC approachthe saddle C , while the RSC evolves along the line[0 , , (cid:104) Ω (cid:105) ] from infinite (cid:104) Ω (cid:105) to (cid:104) Ω (cid:105) = 1 at the sink C .The case with open topology with one RSC can be2 FIG. 6:
Phase space evolution of a high density re–collapsing universe with closed topology . Solutioncurves evolve from the past attractor or source C , approachthe saddle C and diverge as ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − ≤ r ≤ π/ r = 0 evolving fromthe source C upwards. The off–center layers approachingthe saddle C are those with r ≈ C are close to the “turning value” r = π/
2. Layersmarked by π/ ≤ r ≤ π , including the second RSC at r = π ,would have identical evolution as those shown in the figure,with those marked by r ≈ π near the saddle C .FIG. 7: Phase space evolution of the collapsing stagein the re–collapsing universe with closed topology .Only the end stage of the collapse is described by plottingthe phase space variables by means of system (60). Solutioncurves evolve from infinite values of (cid:104) Ω (cid:105) and S to the point C , which is now a future attractor or sink. The curves alsoapproach the saddle C . Only the range 0 ≤ r ≤ π/ r = 0 evolving from infinite valuesof (cid:104) Ω (cid:105) to the C , which is now a sink. Layers marked by π/ ≤ r ≤ π , including the second RSC at r = π , would haveidentical evolution as those shown in the figure, with thosemarked by r ≈ π near the saddle C . obtained with the following initial value functions: m i = 5 .
01 + tan r ,k i = 2 .
01 + tan r ,Y i = H − tan r, (81)As in the case with closed topology, dust layers emergefrom a big bang singularity reach a maximal expansionand then re-collapse in a big crunch, but now dust layersfar from the RSC at r = 0 re-collapse in very long timewhich becomes infinite in the limit r → π/ (cid:104) Ω i (cid:105) ≥ (cid:104) Ω (cid:105) ≥ ξ .The curves (except the RSC) start at the past attractoror source C (big bang), approach the saddle C anddiverge upwards as ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105)− H → C towards (cid:104) Ω (cid:105) → ∞ . FIG. 8:
Phase space evolution of high density re–collapsing universe with open topology.
Since (cid:104) Ω (cid:105) → ∞ as ξ → (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − (cid:104) Ω (cid:105) , The curves startat the past attractor or source C (big bang) approach thesaddle C and diverge upwards (maximal expansion). TheRSC evolves from the source C to (cid:104) Ω (cid:105) → ∞ . The collapsingstage is not described. As in the case with spherical topology, the collapsingstage cannot described with (51). We use then the solu-tion curves of (60) to plot the phase space variables in thecollapsing stage in figure 9, which is very similar to figure7: the critical point C is now a future attractor or sinkassociated with the collapsing singularity (“big crunch”),while curves near the RSC approach the saddle C andthe RSC evolves from infinite (cid:104) Ω (cid:105) downwards to (cid:104) Ω (cid:105) = 1along the line [0 , , (cid:104) Ω (cid:105) ].3 FIG. 9:
Phase space evolution of the collapsing stage inthe re–collapsing universe with open topology . Phasespace variables are plotted by solving system (60). Solutioncurves evolve from infinite values of (cid:104) Ω (cid:105) and S to the point C , which is now a future attractor or sink. The RSC at r = 0 evolves from infinite values of (cid:104) Ω (cid:105) to the sink C . D. Structure formation scenario
A very important type of configuration, associatedwith structure formation scenarios, follows by having alocal region around a RSC evolving from a big bang sin-gularity and collapsing into a black hole (second singular-ity), while the region “outside” expands perpetually as asort of “cosmic background”. This type of configurationshave also been examined with the LTB original variables(see [1, 5, 13, 14, 15, 16] and references quoted therein).The structure formation scenario requires that (cid:104) Ω i (cid:105)− r .It can be constructed by choosing m i > m i < k i must bepositive at the RSC and become asymptotically negative.Such a choice will yield for re–collapsing curves (cid:104) Ω i (cid:105) > (cid:104) Ω (cid:105) ≤ ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − < (cid:104) Ω i (cid:105) ≤ < (cid:104) Ω (cid:105) ≤ ξ . The following initialvalue functions fulfill these conditions: m i = m + m − m r , m = 1 . , m = 0 . k i = k + k − k r , k = 0 . , k = − . Y i = H − r, (82)We have depicted in all the graphs the collapsing solu-tion curves in red and the perpetually expanding ones inblue. The collapsing curves evolve as the solution curvesof figures 6 and 8, except that now these curves onlycover a bounded range of r around the RSC at r = 0, with expanding curves filling the asymptotic range of r .For the collapsing curves the function H obtained bysolving the system (60) for initial conditions (82) changessign (passes from positive to negative). This function isplotted in figure 10 showing collapsing and perpetuallyexpanding layers respectively in red and blue colors. FIG. 10:
The function H in a structure formationscenario. This function is obtained by solving the system(60) for initial conditions (82). For re–collapsing curves with (cid:104) Ω i (cid:105) > (cid:104) Ω i (cid:105) ≤ It is also interesting to visualize the the dimensionlessdensity m (see figure 11). For collapsing curves (red)around the RSC at r = 0 this function decreases from in-finite values as dust layers expand from the initial singu-larity, reaches a minimal value (maximal expansion) andthen increases again as curves go into a collapsing singu-larity. However, for perpetually expanding curves (blue)this function decreases monotonically from the initial bigbang.The evolution of the solution curves of (51) in phasespace is more complicated than for previously presentedconfigurations. Collapsing curves (red) with (cid:104) Ω i (cid:105) > C (big bang), approaching the saddle C anddiverging upwards ( (cid:104) Ω (cid:105) → ∞ ) as ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − < (cid:104) Ω i (cid:105) ≤ C (big bang) approaching the saddle C and descendingtowards the future attractor given by line of sinks C .This is shown in figures 12 and 13.As shown by figure 12, all curves start at the pastattractor C (big bang), those near the RSC ( r = 0) ap-proach the saddle C (see figure 12) while those near thesplit from collapsing to expanding approach the saddle4 FIG. 11:
Density for a structure formation scenario.
The initial big bang singularity is now associated to infinitedensity that decays as dust layers expand. Curves (red) with (cid:104) Ω i (cid:105) > (cid:104) Ω i (cid:105) ≤ τ grows. C (see figure 13). All collapsing curves end up divergingupwards (cid:104) Ω (cid:105) → ∞ as ξ → ln[ (cid:104) Ω i (cid:105) / ( (cid:104) Ω i (cid:105) − < (cid:104) Ω i (cid:105) ≤ < (cid:104) Ω i (cid:105) ≤ < (cid:104) Ω (cid:105) ≤ ξ . Expanding curves also start at the source C (big bang), but do not approach the saddle C , insteadthey approach the saddle C and end up in the line ofsinks C . The RSC evolves from the source C upwards.Both the collapsing and expanding curves approach thesaddle C , which splits the former curves going upwardsto (cid:104) Ω (cid:105) → ∞ and the latter downwards towards the lineof sinks C . The curve with (cid:104) Ω i (cid:105) = 1 remains in theplane (cid:104) Ω (cid:105) = 1, evolving from C towards C . See figure13.The collapsing curves behave near the collapse as thecurves in figures 7 and 9. We show in figure 14 the col-lapsing stage of the curves plotted in figure 13, for a rangeof r around values where (cid:104) Ω i (cid:105) − C , which is a future attractor or sink. As thecurves approach the splitting saddle C , the collapsingones go to infinite (cid:104) Ω (cid:105) and S , while expanding ones fallinto the sinks C . Because of the restricted domain of r that we used to illustrate the behavior near C , the fig-ure doesn’t show curves near the RSC approaching thesaddle C and the RSC evolving from infinite (cid:104) Ω (cid:105) into C , which is a sink for this worldline. ! a rc t a n ( ) a rc t a n ( ) " ( m ) a rc t a n ( ) S ! = C C C C C FIG. 12:
Phase space picture of the structure forma-tion scenario.
Red color corresponds to collapsing curveswith (cid:104) Ω i (cid:105) >
1, while perpetually expanding curves with0 < (cid:104) Ω i (cid:105) < C , perpetually expanding curvesdescend into the future attractor given by the line of sinks C ,collapsing curves approach the saddle C and go upwards todiverging (cid:104) Ω (cid:105) (maximal expansion before collapse). The sad-dle C splits expanding from collapsing curves. The curvewith (cid:104) Ω i (cid:105) = 1 (not shown) would go from C towards thissaddle. ! = ! " (m) S C C C C C FIG. 13:
Close up of the solution curves near the sad-dle C .Red color corresponds to layers with (cid:104) Ω i (cid:105) > (cid:104) Ω i (cid:105) <
1. The plot only depicts solution curves thatcorrespond to r values near the change from re–collapsing toperpetually expanding behavior (where H becomes negative).Notice how this saddle splits collapsing curves going upwardsto (cid:104) Ω (cid:105) → ∞ while expanding ones go downwards to line ofsinks C . Notice how the curve corresponding to (cid:104) Ω i (cid:105) = 1remains in the plane (cid:104) Ω (cid:105) = 1, starting from the source C andending in the saddle C . FIG. 14:
Collapsing stage for curves of figure 13.
Thecurves come from solving (60) for initial conditions (82). Theend state of collapsing curves (red) is the point C , whichis now a sink. Collapsing curves go upwards taking infinitevalues of (cid:104) Ω (cid:105) and S . Notice how the saddle C splits col-lapsing curves going upwards to (cid:104) Ω (cid:105) → ∞ while expandingones (blue) go downwards to the future attractor given by thesinks C . E. Re–collapsing wormhole
In dust configurations without a RSC (“wormhole”topology), the T can be homeomorphic to either S × R or to S × S . These configurations were examined inreference [17] (see also [5, 6]) and can be constructedwith Y i ( r ) having no zeroes in all its domain. Since anychoice of such function will necessary fulfill Y (cid:48) i ( r ∗ ) = 0for at least a “turning point” r = r ∗ , these configurationscan only be regular if they comply with the regularitycondition (91), which implies a re–collapsing “elliptic”dynamics with (cid:104) Ω i (cid:105) − > (cid:104) ρ (cid:105) and (cid:104) R(cid:105) in (10) and(11) we took the RSC as the lower bound of the integrals.This is a boundary condition on these definitions ensuringthat regularity at that RSC is fulfilled [2, 5, 6]. If thereare no RSC then this boundary condition can be specifiedby demanding that m i and k i vanish at an asymptoticvalue of r along the T i . Such value can be then takenas the lower bound for the integrals. However, whilethis solves the mathematical problem, without a RSCthe interpretation of (cid:104) ρ (cid:105) and (cid:104) R(cid:105) as average functionsbecomes less clear.A choice of initial value functions for the re–collapsing wormhole is given by: m i = m (1 + α )1 + α sec r , m = 1 . , α = 5 . k i = k (1 + β )1 + β sec r , k = 0 . , β = 3 . Y i = H − sec r, (83)Notice that the “turning point” (or “throat”) in this ex-ample is located at r = r ∗ = 0 with (91) taking the form (cid:104) k i (cid:105) (0) = 1. Also, both m i → k i → r → ± π/ (cid:104) m i (cid:105) and (cid:104) k i (cid:105) by integrals like(10) and (11) by taking the lower integration limit at r = − π/ r up to π/
2. Themetric function Y = Y i (cid:96) obtained from (60) for this caseis depicted in figure 15. FIG. 15:
Curvature radius of dust layers with worm-hole topology.
Notice how for each hypersurface T with τ constant, Y = Y i (cid:96) has no zeros in its regular range and takesits minimal value at the “throat” at r = 0. Dust layers nearthe “throat” re–collapse in very a short interval of τ , whilelayers further away re–collapse in much longer times. The evolution of the wormhole configuration in phasespace is depicted in figure 16. Only curves in the range0 ≤ r < π/ − π/ < r ≤ C associated with an initial big bang singularity and ris-ing towards maximal expansion (cid:104) Ω (cid:105) → ∞ . The lack of aRSC is evident and this means that C is a global pastattractor.It is interesting to notice how the “throat” and curvesnear it rapidly leap into diverging (cid:104) Ω (cid:105) because they reachmaximal expansion and collapse very fast (see figure 15),but curves close to r = π/ (cid:104) Ω (cid:105) = 1 and approach the saddle C where6 FIG. 16:
Dust wormhole
There is no RSC. Notice howcurves near the “throat” at r = 0 rapidly reach maximal ex-pansion with (cid:104) Ω (cid:105) diverging, while curves close to the asymp-totic value r = π/ C with ∆ ( m ) = S = 0. The lackof a RSC implies that C is now a global past attractor. S = ∆ ( m ) = 0. We took as example a very simple sym-metric form of T i with only one turning point, but itis possible to choose any number of such points [17]. Itis also possible to construct a configuration whose hy-persurfaces T are homeomorphic to torii S × S (see[5, 6]).The collapsing stage for wormhole configurations isvery similar to the stages depicted in figures 7 and 9.We show in figure 17 the expanding and collapsing stagefor curves near the “throat” r = 0 of this configuration.It is easy to recognize some of the curves plotted in figure16, which continue towards infinite values of (cid:104) Ω (cid:105) and S ,in order to plunge downwards to end in the critical point C , which is a future attractor or sink in this stage. F. Vacuum limit
The evolution in phase space takes place in the plane∆ ( m ) = −
1. As we commented in the previous section,for zero binding energy (cid:104) Ω i (cid:105) = 1 (Lemaˆıtre coordinates[8]) it goes from source C to sink C , while for negativebinding energy (cid:104) Ω i (cid:105) > C , approaches saddle C and goes upwardsto (cid:104) Ω (cid:105) → ∞ . The phase space evolution in this case isqualitatively analogous to that of a wormhole topology,which is not surprising since the Schwarzschild–Kruskalspace–time has no RSC, hence C will be a global pastand future attractor.We provide here an example with positive binding en-ergy (cid:104) Ω i (cid:105) <
1. As shown in figure 18 the solution curvesemerge from the global past attractor source C , some FIG. 17:
Dust wormhole: collapsing stage
The figuredescribes the expanding and collapsing stages. The curves offigure 16 near the “throat” r = 0 can be recognized and seengoing upwards into infinite values of (cid:104) Ω (cid:105) and S , all in orderto plunge downwards into C . Notice that the lack of a RSCmakes this critical point a global past and future attractor. curves approach the saddle C and some approach thesaddle C and all terminate in the global future attractoror sink C . FIG. 18:
Vacuum case
The solution curves correspond tothe Schwarzschild–Kruskal space–time in coordinates given byradial geodesic congruence that is perpetually expanding withpositive binding energy (thus (cid:104) Ω i (cid:105) < ( m ) = −
1. The curves evolve from theglobal past attractor or source C towards the global futureattractor or sink in C with S = 0, approaching the saddles C and C . X. DISCUSSION AND CONCLUSION
We have examined the class of dust LTB models interms of the techniques generically known as “dynamicalsystems” as they evolve in a 3–dimensional phase spaceparametrized by the variables [ S, ∆ ( m ) , (cid:104) Ω (cid:105) ], whose phys-ical and geometric interpretation is intuitive and straight-forward. We can say in general that if initial condi-tions are selected so that “shell crossing” singularities areavoided (condition (68)), then the phase space variables[ S, ∆ ( m ) ] will be bounded. Also, the sign of (cid:104) Ω i (cid:105) ( r ) − r = r determines if it willevolve for all ξ in the domain of (cid:104) Ω (cid:105) ( ξ, r ) − C , characterized by S = 1 / (cid:104) Ω (cid:105) = 1 and ∆ ( m ) = − C is indeedthe only past attractor for all ‘off center” solution curves:the asymptotic past regime for all curves is given by thelimit ξ → ∞ (or (cid:96) → (cid:104) Ω (cid:105) → ( m ) and∆ ( k ) as∆ ( m ) = m i (cid:104) m i (cid:105) Γ − , ∆ ( k ) = k i (cid:104) k i (cid:105) Γ − . (84)But, from the analytic forms given in Appendix C andfor finite and regular m i and k i , we have in the limit (cid:96) → ∼ (cid:96) − / for (cid:104) Ω i (cid:105) = 1 and Γ ∼ (cid:96) − for (cid:104) Ω i (cid:105) (cid:54) = 1. Thus, from (84),we have ∆ ( m ) → − ( k ) → − (cid:104) Ω (cid:105) →
1, wehave S → / C , which is in this case a fu-ture attractor for these curves and is associated with thecollapsing (“big crunch”) singularity. It is only the factthat RSC’s do not evolve nor terminate (collapse) in C what prevents this attractor to be a global one, thoughin configurations lacking a RSC (wormhole topology andvacuum limits) this attractor is indeed global.In order to prove that near C we have a self–similarregime we use the definitions (6), (17), (36) and (44).The condition S = 1 / Y C = r (1 + ϑ ) / , Y (cid:48) C = 1 + ϑ/ ϑ ) / , with ϑ ≡ τr , (85)where τ = H ct and we have eliminated two functions ofthe form a ( ct ) and b ( r ) by trivially re–labeling the timeand radial coordinates. Since (cid:104) Ω (cid:105) = 1 implies (cid:104) k i (cid:105) =0 = (cid:104) R i (cid:105) , the metric (1) near C takes the form of the self–similar sub–case of LTB metrics [18] ds = − dτ + A ( ϑ ) dr + r B ( ϑ ) ( dθ + sin θ dφ ) , with A ( ϑ ) = Y (cid:48) C , rB ( ϑ ) = Y C . (86)The fact that ∆ ( m ) = − C is aSchwarzschild vacuum, but that Γ → ∞ as (cid:96) →
0. Eval-uating Γ for the self–similar metric forms (85) and (86)with Y i = Y C ( t i , r ), we getΓ = 3 + ϑ ϑ i ϑ i ϑ (87)with ϑ i = τ i /r , which shows how Γ → ∞ at the locus ofthe big bang singularity ϑ = − Y C = rB ( ϑ ) → ϑ = − ϑ > − Y and Y (cid:48) given by the same expressions as (85), but with aminus sign, so that the “big crunch” happens before thelocus of Γ = 0).Therefore, we can consider the source C as marking aself–similar limit and since all “off center” solution curvesstart at this critical point (associated with a big bangsingularity), we can say that dust LTB models have aself–similar behavior near this singularity. Since collaps-ing curves in re–collapsing configurations also end in C (a sink), we have self–similar behavior near the collapsingsingularity.Another interesting feature worth remarking is the roleof the saddle C in the “structure formation” configura-tions discussed in section IX-D. As shown by figures 13and 14, this saddle splits collapsing curves (cid:104) Ω i (cid:105) > < (cid:104) Ω i (cid:105) <
1. The former evolveback into C (now a sink) while the latter fall into theline of sinks C . This behavior denotes a basic instabil-ity of the invariant subspace of the “parabolic” evolution (cid:104) Ω (cid:105) = 1, since for any solution curve in this space anyarbitrarily small perturbation δ (cid:28) (cid:104) Ω i (cid:105) = 1 + δ will trigger a radically different “repul-sive” evolution away from (cid:104) Ω (cid:105) = 1, either to a collapsingregime ( δ >
0) or to perpetual expansion ( δ < C is marked by ∆ ( m ) = − S = − / (cid:104) Ω (cid:105) = 1, it can be shown from the ana-lytic solutions in Appendix C and from numerical teststhat Γ is finite as the curves approach these values (Γdiverges for later times as the collapsing curves actu-ally collapse). Therefore, the fact that ∆ ( m ) → − C unstable conditions similar tothose of a Schwarzschild–Kruskal vacuum parametrizedby Lemaˆıtre coordinates made with radial geodesics withzero binding energy ( (cid:104) k (cid:105) = 0 or (cid:104) Ω (cid:105) = 1). In fact, thepoint C represents the global future attractor (sink) forthis marginally bound vacuum configuration.The association of the saddle C with a Schwarzschild–Kruskal vacuum roughly conforms with an intuitive pic-ture. A very rough, hand waving, description of condi-tions near C could be: collapsing dust layers “retreat”8inwards towards a smaller collapsing region, while ex-panding layers “retreat” outwards in a fast expansion,thus we have a sort of unstable “evacuation” effect that isroughly similar to a Schwarzschild–Kruskal vacuum withthe Schwarzschild mass being the accumulated effectivemass M = (1 / (cid:104) m i (cid:105) Y i of the collapsing layers. It isknown [5, 6] that in this type of configuration the more“external” layers in the collapsing region ( (cid:104) Ω i (cid:105) >
1) takeinfinite time to collapse, so that the resulting black holeperpetually takes accretion from these “border” dust lay-ers. However, we can assume that the mass contributionfrom this accretion is small and lengthy to come, so thatonce most “inner” layers with (cid:104) Ω i (cid:105) > (cid:104) Ω i (cid:105) = 1 do perceive asort of approximate Schwarzschild–Kruskal vacuum withzero binding energy.We feel that the qualitative numerical treatment pre-sented in this article can serve for studying applicationsof astrophysical interest of dust LTB solutions. Perhapsfor this purpose the system (60) can be more useful thanthe dynamical system (51). It is also interesting to seethis methodology generically applied to generalizations ofthese solution, such as LTB–de Sitter (dust with cosmo-logical constant) or the most general source compatiblewith the LTB metric: a mixture of dust and an inho-mogeneous and anisotropic fluid, which can respectivelymodel cold dark matter and dark energy. Extensions ofthis work along these lines are presently under consider-ation. Acknowledgments
The author acknowledges support from Instituto deF´ısica, Universidad de Guanajuato.
APPENDIX A: Regularity and singularities in dustLTB solutions
An extensive discussion of regularity conditions andimportant geometric features of dust LTB solutions canbe found in [5, 6]. These features were also discussedwith the variables of this article in [2]. We present herea brief review.
Regular Symmetry Centers (RSC) and topology ofthe T i LTB dust solutions admit up to two RSC’s, which arethe regular wordlines made up by the space–time evolu-tion of the fixed point of the rotation group SO(3) wherethe orbits have zero area. In terms of the initial valuefunctions, the RSC is the comoving worldline marked bya zero of the function Y i ( r ). Since the freedom to choosea radial coordinate in the LTB metric means that thisfunction can be arbitrarily chosen, it is useful to select it as Y i = H − f ( r ), where f ( r ) is a (at least) C real func-tion that describes in a simple manner the topologicalconfiguration or homeomorphic class for an initial T i .We have the following choices: • “Open topology”: T i homeomorphic to R , f ( r ) = tan r ≤ r < π/ r = 0 . (88). • “Closed topology”: T i homeomorphic to S , f ( r ) = sin r ≤ r ≤ π Two RSCs at r = 0 , π. (89). • “Wormhole”: T i homeomorphic to S × R , f ( r ) = sec r − π/ < r < π/ . (90).For the open topology f ( r ) is a monotonously increasingfunction. In fact, we could have just selected f = r for simplicity, but the choice f = tan r allows one toexamine the asymptotic behavior associated with Y i →∞ in solution curves of (51).In the wormhole topology the function Y i must be pos-itive and without zeroes. Numerical tests show thatmonotonous functions without (at least) one “turningvalue” (zero of Y (cid:48) i ) yield finite proper radial distancesalong the T i as r → ±∞ (depending on the choice of Y i ). This indicates that such functions are inadequatechoices of the radial coordinate, as they do not providea full analytic extension of the manifold. Therefore, weselect Y i > r = r ∗ withinthe open range of r . In the choice (90) this value is r ∗ = 0(the “throat” of the wormhole).In the “closed” topology we must have a turning value,which with the choice (89) is r ∗ = π/
2. Since in bothcases (closed and wormhole topologies) we have Y (cid:48) i ( r ∗ ) =0, regularity of the metric function g rr Y (cid:48) (1 − K ) / = (cid:96) Γ Y (cid:48) i [1 − H (cid:104) k i (cid:105) Y i ] / requires (see [2, 5, 6] and [19]) that R i be selected sothat K ( r ∗ ) = 1, or equivalently:16 (cid:104) R i (cid:105) ( r ∗ ) = H (cid:104) k i (cid:105) ( r ∗ ) = 1 Y i ( r ∗ ) . (91)Given a choice of homeomorphic class for the T i all reg-ular T have the same class (see [2]). All configurationsexamined numerically corresponding to closed or worm-hole topologies comply with (91).9 Singularities
LTB dust solutions present two types of curvature sin-gularities marked by ( ct, r ) values such that: (cid:96) ( ct, r ) = 0 , “bang” or “crush” , (92a)Γ( ct, r ) = 0 , “shell crossing” (92b)The “bang” or “crush” singularities are an inherent fea-ture and cannot be avoided, though it is always possibleto select initial conditions so that “shell crossing” do notarise (either for all the evolution time or for a given range t ≥ t i ).Initial conditions to avoid unphysical shell crossing sin-gularities were given in terms of the original LTB vari-ables by Hellaby and Lake [7] (see [5, 6]), and in termsof the initial value functions defined in this article by[2]. For broadly generic configurations that exclude veryspecial (and unstable) features, we have:for 0 < (cid:104) Ω i (cid:105) ≤ − < ∆ ( m ) i ≤ , − / < ∆ ( k ) i < / , (93a)for (cid:104) Ω i (cid:105) > − < ∆ ( m ) i ≤ , − / < ∆ ( k ) i ≤ , πP i (cid:20) ∆ ( m ) i −
32 ∆ ( k ) i (cid:21) ≥ [ P i Q i −
1] ∆ ( m ) i + (cid:20) P i Q i − (cid:21) ∆ ( k ) i ≥ , , (93b)where P i = √ − xx / ,Q i = arccos(1 − x ) − (cid:112) x (2 − x ) ,x ≡ (cid:104) R i (cid:105) κc ρ i = 2 (cid:104) k i (cid:105)(cid:104) m i (cid:105) = 2[ (cid:104) Ω i (cid:105) − (cid:104) Ω i (cid:105) . Equations (93) clearly show how initial conditionscharacterized by density and 3–curvature clumps (see(15)) lead generically to avoidance of a shell crossing sin-gularity if 0 < (cid:104) Ω i (cid:105) ≤
1, and as a corollary, to finiteand bounded values of ∆ ( m ) and S along solution curves.This regularity is more stringent for the case (cid:104) Ω i (cid:105) > T i it is still possible tohave density or 3–curvature voids forming along the evo-lution of the curves, all without violating (68). This wasdiscussed by [10], but their methodology is too compli-cated. It is more clearly seen here from the fact that thescaling laws (23) and (24) do not, necessarily, forbid ∆ ( m ) and/or ∆ ( k ) from reversing the sign of the initial ∆ ( m ) i and ∆ ( k ) i . The configuration discussed in section IX–B isa numerical example of this type of “clump turning intovoid” configuration. APPENDIX B: An equivalent dynamical system
We have examined the evolution of solution curves ofthe autonomous system (51) in the 3–dimensional phasespace parametrized by [ S, ∆ ( m ) , (cid:104) Ω (cid:105) ]. However, even if(51) has only derivatives with respect to the evolutionparameter ξ , the differential equations are still PDE’s.We provide below formal proof that such a system ofautonomous evolution PDE’s is equivalent to a system ofODE’s with initial conditions restricted by the fulfillmentof space-like constraints.Let F ( r ) ⊂ R and F ( ξ ; r ) ⊂ R be, respectively, the reg-ularity range of r and the maximal range of extendibilityof ξ for a given r , the system (51) can be associated withthe flow Φ ξ : R → R of the vector field ∂/∂ξ given byΦ ξ ( (cid:126)X i ) = (cid:126)X, (95)where (cid:126)X i : F ( r ) → R is the initial state represented as acurve in R given parametrically as (cid:126)X i ( r ) = [ S i ( r ) , ∆ ( m ) i ( r ) , (cid:104) Ω i ( r ) (cid:105) ] , (96)and (cid:126)X : F ( ξ ; r ) × F ( r ) → R is a solution surface as-sociated with (cid:126)X i by means of (95). This surface can berepresented parametrically as the following surface in R : (cid:126)X ( ξ, r ) = [ S ( ξ, r ) , ∆ ( m ) ( ξ, r ) , (cid:104) Ω( ξ, r ) (cid:105) ] , (97)so that (cid:126)X i ( r ) = (cid:126)X (0 , r ). Given a point (cid:126)X i ( r ) ∈ R thereis a unique orbit or solution curve C r : F ( ξ ; r ) → R of(51) which can be represented parametrically as: C r ( ξ ) = (cid:126)X ( ξ, r ) , (98)so that C r (0) = (cid:126)X i ( r ). The set of solution curves C r ( ξ )for all r ∈ F ( r ) generates the solution surface (cid:126)X ( ξ, r ).In order to relate (51) with a system of ordinary dif-ferential equations, we notice that every solution curve(98) associated with a fixed r = r is a solution of thefollowing system: dS dξ = S (3 S −
1) + 12 (cid:16) ∆ ( m )0 + S (cid:17) (cid:104) Ω (cid:105) , (99a) d ∆ ( m )0 dξ = 3 S (cid:104) ( m )0 (cid:105) , (99b) d (cid:104) Ω (cid:105) dξ = (cid:104) Ω (cid:105) [ (cid:104) Ω (cid:105) − , (99c)where the functions ( S , ∆ ( m )0 , (cid:104) Ω (cid:105) ) : F ( ξ ; r ) → R aregiven by S = S ( ξ, r ) , (100a)∆ ( m )0 = ∆ ( m ) ( ξ, r ) , (100b) (cid:104) Ω (cid:105) = (cid:104) Ω( ξ, r ) (cid:105) , (100c)0and comply with the initial conditions:[ S (0) , ∆ ( m )0 (0) , (cid:104) Ω (cid:105) (0)] = (cid:126)X i ( r ) . (101)The fact that every solution curve of (51) will be a solu-tion of a system like (99) for initial conditions given asin (101), suggests considering a system of ordinary dif-ferential equations associated with the same flow as (95) dsdξ = s (3 s −
1) + 12 ( δ + s ) ω, (102a) dδdξ = 3 s [ 1 + δ ] , (102b) dωdξ = ω [ ω − , (102c)for an initial state: (cid:126)x = [ s (0) , δ (0) , ω (0)] = [ s , δ , ω ] , (103)and the solution surface given by a parametric form sim-ilar to (97): (cid:126)x = [ s ( ξ, s ) , δ ( ξ, δ ) , ω ( ξ, ω )] , (104)so that (cid:126)x = (cid:126)x (0). Since the solution surface of thissystem can depend smoothly on initial conditions (cid:126)x ,each one of the functions ( s, δ, ω ) solving (102) is aone–parameter family of functions. In principle, thereare many ways in which the initial state (103) can beparametrized. In particular, this state can be givenas any curve [ s ( α ) , δ ( α ) , ω ( α )] in R that intersectsthe solution curves (cid:126)x ( ξ, (cid:126)x ) only in one point, hencewe can always find suitable one–parameter functions s ( α ) , δ ( α ) , ω ( α ) such that (cid:126)x is expressible as: (cid:126)x ( ξ, α ) = [ s ( ξ, s ( α )) , δ ( ξ, δ ( α )) , ω ( ξ, ω ( α ))] . (105)If, in particular, the initial state (103) is parametrizedas: (cid:126)x ( α ) = [ s ( α ) , δ ( α ) , ω ( α )] = (cid:126)X i ( α ) , (106)with α ∈ F ( r ) , then for every α there is a C ∞ inclusionmap ψ : R → R mapping every solution curve of (51)given by (98) for a fixed r ∈ F ( r ) to a solution curve ofthe system (102) given by (105) with a fixed α = r . Thatis ψ ( (cid:126)X ( ξ, r )) = (cid:126)x ( ξ, α ) , (107)and so, the solution surfaces of (51) are subsets of the so-lution surfaces of (102) obtained by restricting the initialstates by the constraint (106).Notice that (102) admits many solution curves thatdo not comply with (106), and so are incompatible withthe constraints (53)–(56c) and so they are incompatiblewith the evolution of an inhomogeneous dust source (LTBsolution) complying with Einstein’s field equations ex-pressed in terms of the 3+1 decomposition [4] under the regularity assumptions that we have made (see AppendixA).However, the solution curves of (102) satisfying (106)are equivalent (under (107)) to solution curves of (51)associated with the flow (95). Since (102) is a systemof autonomous ODE’s, geometric features such as criti-cal points and invariant spaces are well defined, so thequalitative analysis of (51) can be performed on (102) aslong as we only consider its solution surfaces satisfying(106). In other words, we can say that the dynamical sys-tem associated to inhomogeneous dust LTB solutions isa sort of “reduced” version of (102) with its initial statesrestricted by (106). Since, as we have shown in sectionsIV and V, the space–like constraints are satisfied at allspace–like slices once they hold at the initial T i , the in-corporation of the radial dependence of the solutions onlyrequires this restriction on the initial states.In order to simplify the notation we found it useful tokeep the same names for the variables S, ∆ ( m ) , (cid:104) Ω (cid:105) andto keep referring to (51), with the understanding thatthe system we are really considering is (102), thus everymention of solution curves and surfaces of (51) refers tosolution curves and surfaces of (102) mapped by (107)and complying with (106). Appendix C: Analytic solutions
Analytic solutions for LTB dust solutions follow fromsolving the evolution equation (3):˙ Y = 2 MY − K. Using the variables defined in sections III and IV, thisequation becomes equation (59): (cid:20) ∂(cid:96)∂τ (cid:21) = H i (cid:20) (cid:104) Ω i (cid:105) (cid:96) − ( (cid:104) Ω i (cid:105) − (cid:21) . where: (cid:96) = YY i , τ = H ct M = (cid:104) m i (cid:105) H Y i , K = (cid:104) k i (cid:105) H Y i , (cid:104) Ω i (cid:105) = (cid:104) m i (cid:105)H i = (cid:104) m i (cid:105)(cid:104) m i (cid:105) − (cid:104) k i (cid:105) , (cid:104) Ω i (cid:105) − (cid:104) k i (cid:105)H i = (cid:104) k i (cid:105)(cid:104) m i (cid:105) − (cid:104) k i (cid:105) , We summarize the solutions of this equation below
Parabolic solutions
We have: (cid:104) Ω i (cid:105) = 1, so that (cid:104) k i (cid:105) = 0 and H i = (cid:104) m i (cid:105) .These are the only ones that can be given explicitly as (cid:96) = (cid:96) ( τ, r ) from: τ = τ i ± (cid:96) / − H i , (108)1where the ± sign describes expanding (+) and collaps-ing ( − ) dust layers, respectively characterized by τ >τ bb , (cid:96) ,τ > τ < τ bb , (cid:96) ,τ <
0. Notice that the “bigbang time” τ bb = H ct bb is no longer an independentfunction but must be found by setting (cid:96) = 0 in (108): τ bb = τ i ∓ H i , (109)where now expanding and collapsing layers respectivelycorrespond to “ − ” and “+”. Hyperbolic solutions
We have: 0 < (cid:104) Ω i (cid:105) <
1, so that (cid:104) k i (cid:105) < H i = (cid:104) m i (cid:105) + |(cid:104) k i (cid:105)| . In this case we cannot obtain (cid:96) = (cid:96) ( τ, r ) in closed explicit form, but in parametric form:[ (cid:96) ( η, r ) , τ ( η, r )] or implicitly as τ = τ ( (cid:96), r ). The solutionsare given in parametric form as: (cid:96) ( η, r ) = (cid:104) Ω i (cid:105) [cosh η − − (cid:104) Ω i (cid:105) ) , (110a) τ ( η, r ) = τ i ± (cid:104) Ω i (cid:105) [(sinh η − η ) − (sinh η i − η i )]2 H i (1 − (cid:104) Ω i (cid:105) ) / , (110b)where η i = arccosh (cid:18) (cid:104) Ω i (cid:105) − (cid:19) , sinh η i = 2 (cid:112) − (cid:104) Ω i (cid:105)(cid:104) Ω i (cid:105) . (111)and ± denote expanding (+) and collapsing (–) layers.The coordinate locus of the central singularity (the “bigbang time” function τ bb ) can be obtained by setting η = 0in (110b), ( so that (cid:96) = 0) leading to τ bb = τ i ∓ (cid:104) Ω i (cid:105) (sinh η i − η i )2 H i (1 − (cid:104) Ω i (cid:105) ) / , (112)where now “–” and “+” respectively correspond to ex-panding and collapsing layers. Elliptic solutions
We have: (cid:104) Ω i (cid:105) >
1, so that (cid:104) k i (cid:105) > H i = (cid:104) m i (cid:105) −(cid:104) k i (cid:105) . The parametric solutions are (cid:96) ( η, r ) = (cid:104) Ω i (cid:105) [1 − cos η ]2( (cid:104) Ω i (cid:105) − , (113a) τ ( η, r ) = τ i + (cid:104) Ω i (cid:105) [( η − sin η ) − ( η i − sin η i )]2 H i ( (cid:104) Ω i (cid:105) − / , (113b) where η i = arccos (cid:18) (cid:104) Ω i (cid:105) − (cid:19) , sin η i = 2 (cid:112) (cid:104) Ω i (cid:105) − (cid:104) Ω i (cid:105) . (114)As with the hyperbolic case, the locus of the central sin-gularity (the “big bang time”, τ bb ), follows by setting η = 0 in (113b), we get τ bb = τ i − (cid:104) Ω i (cid:105) ( η i − sin η i )2 H i ( (cid:104) Ω i (cid:105) − / , (115)Dust layers in elliptic solutions reach a maximal expan-sion as (cid:96) ,τ = 0 in (59), corresponding to (cid:96) = (cid:96) max = (cid:104) Ω i (cid:105)(cid:104) Ω i (cid:105) − , η = π, (116a)= τ bb + π (cid:104) Ω i (cid:105) H i ( (cid:104) Ω i (cid:105) − / , (116b)For (cid:96) ,τ < π < η < π ). The coordinate locus of the collapsing singularity(“big crunch”), τ bc , follows from setting η = 2 π in (113b): τ bc = τ bb + π (cid:104) Ω i (cid:105)H i ( (cid:104) Ω i (cid:105) − / . (117)Since the expanding and collapsing singularities arerespectively marked by η = 0 and η = 2 π , then the layers’evolution is symmetric with respect to η = π marking themaximal expansion (116). However, as shown by (115),(116) and (117), the maximal expansion time dependson the IVF’s (cid:104) Ω i (cid:105) and H i , so it does not coincide (ingeneral) with any T (hypersurface of constant τ ) andthe evolution for general IVF’s (cid:104) Ω i (cid:105) , H i is not “time-symmetric” with respect to the maximal expansion.While for the hyperbolic and elliptic cases we do nothave an explicit closed form solution (cid:96) = (cid:96) ( τ, r ), it ispossible to use the parametric solutions presented so farto visualize graphically the function (cid:96) ( τ, r ). This can bedone by plotting the 3–dimensional parametric surfaces:[ τ ( η, r ) , r, (cid:96) ( η, r )] , (118)where the radial dependence is completely specified bythe IVF’s that enter in (110)–(111) and (113)–(114). [1] A Krasinski, Physics in an Inhomogeneous Universe ,(Cambridge University Press, Cambridge, 1998). [2] R.A. Sussman and L. Garc´ıa–Trujillo,
Class. Quantum Grav , , 2897 – 2925, (2002).[3] J Wainwright and G F R Ellis (Eds.), Dynamical Systemsin Cosmology (Cambridge University Press, Cambridge,1997).[4] G. F. R. Ellis and H. van Elst,
CosmologicalModels , Carg`ese Lectures 1995. LANL e–printarXiv: gr-qc/9812046 [5] D. R. Matravers and N. P. Humphreys,
Gen. Rel. Grav. , , 531-552, (2001).[6] N.P. Humphreys, R. Maartens and D. R. Ma-travers, Regular spherical dust spacetimes , LANL e–printarXiv: gr-qc/9804023v1 [7] C. Hellaby and K. Lake,
Ap. J. , 381, (1985). See
Ap.J. , 461, (1986) for erratum of this paper.[8] L.D. Landau and E.M. Lifshitz,
The Classical Theory ofFields, Course of Theoretical Physics Volume 2 . FourthRevised English Edition. Pergamon Press, 1975. See § § Gravitation , Freeman Press, 1972.[10] N. Mustapha and C. Hellaby,
Gen. Rel. Grav. , , 455-477, (2001).[11] W. Bonnor and A. Chamorro, Ap. J. , 21, (1990).[12] K. Bolejko, A. Krasinski and C. Hellaby,
Mon. Not. Roy.Astron. Soc. , , 213-228, (2005).[13] A. Krasinski and C. Hellaby, Phys.Rev. D65 (2002)023501[14] A. Krasinski and C. Hellaby, Phys.Rev. D69 (2004)023502[15] A. Krasinski and C. Hellaby, Phys.Rev. D69 (2004)043502[16] A. Krasinski and C. Hellaby, Phys.Rev. D73 (2006)023518[17] C. Hellaby, Class. Quantum Grav. , , 635, (1987).[18] B.J. Carr and A.A. Coley, Phys.Rev.
D62 , 044023,(2000)[19] W. Bonnor,
Class. Quantum Grav. ,2