Exact spherically-symmetric inhomogeneous model with n perfect fluids
PPrepared for submission to JCAP
Exact spherically-symmetricinhomogeneous model with n perfectfluids Valerio Marra a,b,c and Mikko P¨a¨akk¨onen b,c a Institut f¨ur Theoretische Physik, Universit¨at Heidelberg, Philosophenweg 16,69120 Heidelberg, Germany b Department of Physics, PL 35 (YFL), 40014 University of Jyv¨askyl¨a, Finland c Helsinki Institute of Physics, PL 64, 00014 University of Helsinki, FinlandE-mail: [email protected], mikko.u.paakkonen@jyu.fi
Abstract.
We present the exact equations governing the dynamics of a spherically-symmetricinhomogeneous model with n decoupled and non-comoving perfect fluids. Thanks to the useof physically meaningful quantities we write the set of 3 + 2 n equations in a concise andtransparent way. The n perfect fluids can have general equations of state, thus making themodel extremely flexible to study a large variety of cosmological and astrophysical prob-lems. As applications we consider a model sourced by two non-comoving dust componentsand a cosmological constant, and a model featuring dust and a dark energy component withnegligible speed of sound. Keywords: dark matter theory, dark energy theory, cosmological perturbation theory
ArXiv ePrint: a r X i v : . [ g r- q c ] D ec ontents c s ∼ − w < − W n and δ n functions 27 Spherically-symmetric models are of interest in the study of the nonlinear inhomogeneitiesof the universe for two reasons. First, spherical symmetry is often a reasonable workingapproximation for the inhomogeneities in the real universe. Second, spherical symmetryallows to exactly solve the Einstein’s equations, a highly nontrivial task. Recent researchhas focused mainly on the Lemaˆıtre model [1] (see [2–5] for recent contributions) whichdescribes the dynamics of a spherically-symmetric perfect fluid, and great attention hasreceived its pressureless limit, usually named the Lemaˆıtre-Tolman-Bondi (LTB) model [1,6, 7], which can also include a nonvanishing cosmological constant. LTB metrics have beenused to describe the large-scale inhomogeneities of the late universe as, for example, in Swiss-cheese [8–15], void [16–36] and inhomogeneous [37–39] models for apparent acceleration (see[40, 41] for recent reviews). Besides being useful in understanding the role of voids in the– 1 –rocess of structure formation [42, 43], the LTB model is also suitable to study the dynamicsof spherical collapse in an expanding universe, black-hole formation included [44–48].In the present paper we extend the Lemaˆıtre model to the case of n decoupled andnon-comoving perfect fluids with general equations of state. Our main result is the set of3 + 2 n exact equations governing the dynamics of the model, which we write in terms ofphysically meaningful quantities and express in a concise and transparent way. Many arethe possible applications, both at early and late times. At late time, one can study theevolution of overdensities and underdensities in a universe where dark energy is not thecosmological constant and, in particular, can be inhomogeneous (see, for example, [49, 50]and references therein). This possibility may also induce an inhomogeneous variation offundamental constants, as for instance the fine structure constant if a coupling between thedark energy and the electromagnetic field is allowed (see, for example, [51–55]). Moreover,dark matter and baryons can be described as two separate fluids, the latter possibly featuringpressure. As stressed in [5], the introduction of non-dust equations of state can have a non-negligible effect on the cosmological models and this should be taken into considerationwhile interpreting cosmological datasets. This is, however, a non-trivial task due to theincreased number of free parameters, and such difficulties increase when one considers n fluid components with the consequence that more astronomical and astrophysical input isneeded to constrain the available parameter space. Finally, at early times the contribution ofradiation can be included, which may be relevant both for the understanding of the evolutionof the standard post-inflation inhomogeneities and for the correct modeling of the very largeunderdensities typical of void models [34].This paper is organized as follows. In Section 2 we will go through all the details of ourformalism, and in Sections 3 and 4 we will present the numerical results for the case of twonon-comoving dust components in a flat ΛCDM universe and for spherical collapse in thepresence of dark energy with negligible speed of sound, respectively. Then in Section 5 we willcompare our findings to previous work dealing with exact solutions. Finally we will give ourconclusions in Section 6. In Appendix A we discuss the expansion tensor and in AppendixB the case of the Lemaˆıtre metric, possibly in a non-comoving frame. The definitions of thefunctions used to model the numerical example of Section 3 are given in Appendix C. A spherically-symmetric metric can be written asd s = − e λ d t + Y (cid:48) E d r + Y dΩ , (2.1)where the lapse function e λ , the scale function Y and the curvature function E depend onthe coordinate time t and coordinate radius r , dΩ = d θ + sin θ d φ and we have set c = 1.A prime denotes partial derivation with respect to the coordinate radius r , whereas a dotdenotes partial derivation with respect to the coordinate time t . Comma and semicolonsigns will not denote derivatives. We will use a reference frame comoving with the arbitraryfour-velocity field u α rf , which has components u α rf = ( e − λ , , , Y (cid:48) (cid:54) = 0 and E > − / g rr >
0. We will discuss moregeneral evolutions and shell crossing in forthcoming work.– 2 –he nontrivial components of the Einstein tensor G αβ ≡ R αβ − R g αβ for the met-ric (2.1) are then: G tt = e λ (cid:20) − EY ) (cid:48) Y Y (cid:48) + H A (cid:16) H A + 2 H R (cid:17)(cid:21) , G tr = 2 e λ (cid:18) λ (cid:48) H A − E d Y (cid:48) Y (cid:19) ,G rr = − Y (cid:48) E (cid:18) H A + 2 A A − EY (cid:48) Y λ (cid:48) − EY (cid:19) ,G θθ = G φφ / sin θ = − Y (cid:18) A R + A A + H R H A − E (cid:48) Y (cid:48) Y − EY (cid:48) λ (cid:48) F (cid:19) , where we have defined radial and angular expansion rates: H R = √ EY (cid:48) d Y (cid:48) √ E dτ rf = e − λ ˙ Y (cid:48) Y (cid:48) − E d , (2.2) H A = 1 Y dYdτ rf = e − λ ˙ YY , (2.3)and also radial and angular acceleration rates: A R = √ EY (cid:48) d Y (cid:48) √ E dτ = e − λ ¨ Y (cid:48) Y (cid:48) − E d H R + E d − E d − e − λ ˙ λH R , (2.4) A A = 1 Y d Ydτ = e − λ ¨ YY − e − λ ˙ λH A , (2.5)respectively, and also the following auxiliary quantities: E d = e − λ ˙ E E , E d = e − λ ¨ E E , F = Y (cid:48) Y − Y (cid:48)(cid:48) Y (cid:48) + E (cid:48) E + λ (cid:48) + λ (cid:48)(cid:48) λ (cid:48) . Note that for λ = 0 and ˙ E = 0 the usual LTB expressions are recovered.In the previous equations d/dτ rf = u α rf ∂ α = e − λ ∂/∂t is the derivative with respect tothe proper time of the comoving observer. As we show in Appendix A, the components of theacceleration of u α rf are a r, rf = λ (cid:48) and a t, rf = 0. We see therefore that, for a geodesic reference-frame velocity field, the lapse function λ depends only upon time and can be rescaled suchthat λ = 0 and d/dτ rf = ∂/∂t . From Eqs. (2.2-2.3) it follows then that in a non-geodesicreference frame the comoving observers measure with different proper times. Later we willidentify the reference-frame velocity u α rf with the velocity field of one of the fluid components,and it will turn out that the acceleration a r, rf = λ (cid:48) is sourced by pressure gradients whichpush the observers out of the freely-falling geodesics. From the Einstein’s equations G αβ = κ T αβ we can obtain four independent dynamical equa-tions. We choose to form two of them with the following combinations: Y Y (cid:48) G tt − Y ˙ Y G tr e λ ≡ (cid:0) e − λ Y ˙ Y − EY (cid:1) (cid:48) = κ Y Y (cid:48) T tt − Y ˙ Y T tr e λ , (2.6) Y ˙ Y G rr − Y Y (cid:48) G tr − Y (cid:48) / (1 + 2 E ) ≡ (cid:0) e − λ Y ˙ Y − EY (cid:1) ˙ = κ Y ˙ Y T rr − Y Y (cid:48) T tr − Y (cid:48) / (1 + 2 E ) , (2.7)– 3 –here κ = 8 πG and T αβ is the energy-momentum tensor for an ideal fluid source comprisedof n decoupled components: T αβ = n (cid:88) i =1 T αβi , (2.8)where the energy-momentum tensor of the i :th perfect fluid component is: T αβi = ρ i u αi u βi + p i h αβi , (2.9)where u αi , ρ i and p i are the four-velocity field, energy density and pressure, respectively,and h αβi = g αβ + u αi u βi is the projection tensor on the hypersurface orthogonal to u αi . Theisotropic pressure is related to the energy density by p = w ρ , where the equation of stateparameter w = w ( t, r ) is assumed to be a general function of t and r (see, for example, [56]for the case of a static fluid with anisotropic pressure). In the chosen coordinate system it is u αi = γ i ( e − λ , v i, c , , , (2.10)where v i, c is the coordinate comoving peculiar velocity of the i :th component relative to thereference frame. The proper peculiar velocities and the gamma factors are instead given by v i, p = Y (cid:48) E v i, c and γ i = 11 − v i, p , (2.11)respectively. The covariant components of the total energy-momentum tensor are then: T tt = e λ n (cid:88) i =1 ρ i (cid:2) (1 + w i ) γ i − w i (cid:3) , (2.12) T tr = − Y (cid:48) E e λ n (cid:88) i =1 ρ i (1 + w i ) v i, c γ i , (2.13) T rr = Y (cid:48) E n (cid:88) i =1 ρ i (cid:2) (1 + w i ) v i, p γ i + w i (cid:3) , (2.14) T θθ = T φφ / sin θ = Y n (cid:88) i =1 p i . (2.15)As the third dynamical equation we will consider the combination coming from − G tt + G rr + G θθ + G φφ , which gives the generalization of the acceleration equation: A ≡ A R + 2 A A − √ EY (cid:48) a rf (cid:18) Y (cid:48) Y + F (cid:19) = − κ n (cid:88) i =1 ρ i (1 + 3 w i ) , (2.16)where we used the fact that the acceleration scalar of the reference-frame velocity is a rf = √ EY (cid:48) λ (cid:48) (see Appendix A). The term proportional to a rf gives a “spurious” contributionto the acceleration and vanishes if we use a geodesic reference frame, in which the totalacceleration is A R + 2 A A , similarly to LTB models.Finally, the last independent equation will be simply the G tr = κ T tr component, whichreads: e − λ ˙ E E = κ Y √ E n (cid:88) i =1 ρ i (1 + w i ) v i, p γ i + e − λ ˙ Y √ E a rf . (2.17)– 4 –q. (2.17) shows that the evolution of the curvature is due to two distinct causes. The first isthe effect of having an inhomogeneous multicomponent fluid and goes to zero in the FLRWlimit where the peculiar velocities vanish or if only one fluid is present and its referenceframe is adopted. It is sourced by the energy flux in the radial direction: we remind indeedthat the curvature function E may be interpreted as the total energy of a given shell (seeEq. (2.28)). The second contribution to ˙ E is due to the fact that, generally, we are using anon-geodesic reference frame ( a rf (cid:54) = 0) in which the total energy of a shell r is not conserved. As we are considering an ideal fluid source comprised of n decoupled components, each i :thcomponent will be conserved individually, i.e., ∇ α T αβi = 0, where ∇ α denotes the covariantderivative. For each i :th component we obtain the two following equations where we omitthe index i and multiply for clarity’s sake by e λ and v c Y (cid:48) E , respectively: e λ ∇ α T αt = ( ρ + p ) γ (Θ + v p a ) + v c γ ( ρ + p ) (cid:48) + e − λ ˙ pv p γ + e − λ ˙ ργ = 0 , (2.18) v c Y (cid:48) E ∇ α T αr = ( ρ + p ) γ (cid:0) v p Θ + v p a (cid:1) + v p γ e − λ ( ρ + p )˙+ p (cid:48) v c γ + v c ρ (cid:48) γ v p = 0 , (2.19)where the expansion scalar Θ and the acceleration scalar a have been calculated in AppendixA. By taking the combinations e λ ∇ α T αt − v c Y (cid:48) E ∇ α T αr and e λ ∇ α T αt − v c ∇ α T αr it ispossible to obtain the relativistic energy-conservation and Euler equations, respectively, fornon-comoving fluids: dρdτ = − Θ ( ρ + p ) , (2.20) dpdσ = − a ( ρ + p ) , (2.21)where, analogously to the total derivative with respect to the proper time τ , we defined theconvective derivative with respect to σ along the four-acceleration a µ : ddτ = u µ ∂ µ = γe − λ ∂∂t + γv c ∂∂r , (2.22) ddσ ≡ a µ ∂ µ a = γv p e − λ ∂∂t + γ v c v p ∂∂r . (2.23)Note that in the rest frame of the i :th fluid component the derivative with respect to τ reduces to ( − g tt ) − / ∂/∂t and the derivative with respect to σ to g − / rr ∂/∂r . Eq. (2.21) isactually a contracted form of the standard Euler equation, which reads h µν ∂ ν p = − a µ ( ρ + p ).From Eq. (2.21) it is clear that the four-velocity of a dust ( p = 0) component is geodesic.By identifying the reference-frame velocity u α rf with the ¯ i :th fluid-component velocity u α ¯ i we can obtain the rest-frame expressions of Eqs. (2.20-2.21) by simply setting v c = 0: e − λ ˙ ρ ¯ i = − Θ rf ( ρ ¯ i + p ¯ i ) , (2.24) p (cid:48) ¯ i = − a r, rf ( ρ ¯ i + p ¯ i ) , (2.25) Note, however, that the cosmological constant never sources ˙ E . The other two equations relative to the angular components identically state the spherical symmetry ofthe energy-momentum tensor. – 5 –espectively, where a r, rf = λ (cid:48) and a t, rf = 0 (see Appendix A). The last equation shows thatthe pressure gradients are the cause for the non-geodesity of the reference frame relative tothe ¯ i :th fluid componentFinally, because of the (twice contracted) Bianchi identity it holds ∇ µ G µν = 0 and so: n (cid:88) i =1 ∇ µ T µνi = 0 . (2.26)Consequently, two conservation equations will not be independent and can be discarded.Equivalently, we may replace, if convenient, any other two dynamical equations with thelatter two extra conservation equations. The effective gravitating total mass F is given by2 GF ≡ H A Y − E Y , (2.27)where G is the gravitational constant. This quantity corresponds to the mass in the Schwarzschild’smetric if the metric of Eq. (2.1) was matched to the latter at a given radius r . See Refs. [4,57, 58] for more details. It is interesting to note that Eq. (2.27) may be rewritten as: E = 12 (cid:18) dYdτ rf (cid:19) − GFY , (2.28)from which follows that the curvature function E can be interpreted as the total energy perunit of mass of a shell, and that to it contribute the kinetic energy per unit of mass and thepotential energy per unit of mass due to the total gravitating mass up to that shell. Note thatthanks to spherical symmetry one is able to define a potential energy also in cases far awayfrom nearly Newtonian ones and that the potential energy is related to the curvature [7].In giving the initial conditions it will be useful to divide the total mass F into theseparate masses of the i :th fluid components: F = n (cid:88) i =1 F i . (2.29)Note, however, that this division is irrelevant for the dynamics as only the total F enters theevolution equations. By substituting Eqs. (2.12-2.14) and (2.27) into the (combinations ofthe) Einstein’s Eqs. (2.6-2.7) one obtains the equations for ˙ F and F (cid:48) , which using Eq. (2.29)read: e − λ n (cid:88) i =1 ˙ F i = − πY n (cid:88) i =1 ρ i γ i (cid:104) H A ( v i, p + w i ) + S A (1 + w i ) (cid:105) , (2.30) n (cid:88) i =1 F (cid:48) i = 4 πY n (cid:88) i =1 ρ i γ i v i, c (cid:104) H A (1 + w i ) v i, p + S A (1 + v i, p w i ) (cid:105) , (2.31)where the angular “spatial” expansion rate S A = v c Y (cid:48) Y has been introduced in AppendixA. Note that the Einstein’s equations give Eqs. (2.30-2.31); therefore to write Eq. (2.29)we have put the constants of integration to zero. Note also that the evolution of the i :th– 6 –ass F i is sourced by the gravitational influence of pressures and velocities of the other fluidcomponents, differently from the case of the conservation equations (2.20-2.21) which applyto the fluid components individually: decoupled fluids are indeed still gravitationally coupled.In the case of only one component Eqs. (2.30-2.31) have a simple interpretation as discussedin Appendix B. The total number of independent equations is generally 4 + 2 n −
2, where the first term is thenumber of nontrivial Einstein’s equations, the second comes from the conservation equationsand the last from the Bianchi identity. By including the definition of mass in Eq. (2.27) wehave a total of 3 + 2 n equations. We choose to use all the 2 n conservation equations, two ofwhich will replace the acceleration equation (2.16) and the evolution equation (2.30) for thetotal F . We now list for clarity the full set of 3 + 2 n exact equations governing the dynamicsof a spherically symmetric inhomogeneous model with a multi-component ideal fluid: H A = 2 Y n (cid:88) i =1 GF i + 2 EY , (2.32) e − λ ˙ E E = κ Y √ E n (cid:88) i =1 ρ i (1 + w i ) v i, p γ i + e − λ ˙ Y √ E a rf , (2.33) n (cid:88) i =1 F (cid:48) i = 4 πY n (cid:88) i =1 ρ i γ i v i, c (cid:104) H A (1 + w i ) v i, p + S A (1 + v i, p w i ) (cid:105) , (2.34) dρ i dτ i = − Θ i ( ρ i + p i ) i = 1 , . . . , n , (2.35) dp i dσ i = − a i ( ρ i + p i ) i = 1 , . . . , n , (2.36)where the unknown functions are Y , E , the total F and n copies of ρ i and v i, c . The reference-frame velocity u α rf can be taken to be geodesic so that a rf = λ = 0. Alternatively, one canidentify the reference-frame velocity with the velocity of the ¯ i :th fluid component, u α rf = u α ¯ i ,with the consequence that v ¯ i, c = 0 and λ becomes an unknown function in its stead. Theequations of state w i have to be given as an external input. From Eq. (2.32) it seems clear thatthe dynamics is similar to the FLRW one, the complication being an array of interconnectedequations defining the evolution of gravitating mass and curvature.From Eq. (2.32) one may obtain the age of the universe at a given coordinate radius r .For example, in the case of expansion (positive root, ˙ Y >
0) it is: t A ( r, t ) ≡ t − t B ( r ) = (cid:90) tt B ( r ) dt e − λ ˙ Y (cid:20) GF ( r, t ) Y + 2 E ( r, t ) (cid:21) − / , (2.37)which can be computed once the dynamical equations (2.32-2.36) are solved. The quantity t B is the so-called bang function also present in LTB models. We will give initial conditions at the time ¯ t < t , where t is the presenttime. First we fix the residual gauge freedom in the coordinate r by setting Y ( r, ¯ t ) = a (¯ t ) r ,where a ( t ) is the scale factor of the external FLRW model to which we will match the– 7 –nhomogeneous metric (see below). If one is not interested in embedding the metric, a (¯ t )may be taken as an arbitrary number. Then, E ( r, ¯ t ), F ( r, ¯ t ) and n copies of ρ i ( r, ¯ t ) and v i, c ( r, ¯ t ) are needed.The curvature function E may be specified by demanding a homogeneous age of theuniverse t A at ¯ t , which gives an additional constraint between the total F and E . FromEq. (2.37), in the approximation that F ( r, t ) (cid:39) F ( r, ¯ t ) and E ( r, t ) (cid:39) E ( r, ¯ t ) for t < ¯ t , we canobtain: t A ( r, ¯ t ) ≡ ¯ t − t B ( r ) (cid:39) (cid:90) Y ( r, ¯ t )0 dy e − λ ( r, ¯ t ) (cid:20) GF ( r, ¯ t ) y + 2 E ( r, ¯ t ) (cid:21) − / , (2.38)which can be used without solving the dynamical equations to demand at ¯ t the homogeneousbig bang, t B ( r ) = 0.Note then that in the set of Eqs. (2.32-2.36) there is not a dynamical equation for F ,which can always be found by integrating Eq. (2.34) for a constant t once the boundarycondition F (0 , t ) = 0 is given. In particular, Eq. (2.34) relates the initial conditions for F , ρ i and v i, c , which cannot be given independently, as it is intuitively clear. It is important topoint out that while the initial conditions can be given individually for the different fluidsby considering Eq. (2.34) as n independent equations, this will not be true at later timesfor which only the total F matters. In other words, one cannot obtain the individual F i byintegrating the corresponding i -piece of Eq. (2.34) and then combining them to obtain F .This also means that for t > ¯ t there is not a meaningful way to define the individual F i ,which are not physically relevant. We can instead define the invariant mass M as shown, inthe example of Section 3, by Eq. (3.3).If we identified the reference-frame velocity with the velocity of the ¯ i :th fluid component, u α rf = u α ¯ i , we may have to give initial conditions for λ (and not v ¯ i, c which is identically zero).Note, however, that ˙ λ only appears in the acceleration equation which we have discarded.As with F , it is therefore enough to give a boundary condition (see below) in order to obtainwith Eq. (2.25) the lapse function at ¯ t for any r .We have now all the functions and their r -derivatives on the hypersurface of constant¯ t . The next step is to evolve the model forwards in time along the worldlines of constant r ,which is done by simultaneously integrating the first-order in t differential equations listed in(2.32-2.36). While the t -derivatives for Y , E and ρ i are evident, ˙ v i, c enters only through theexpansion and acceleration scalars. Finally, λ ( r, t ) and F ( r, t ) are found by integration forthe needed constant t . The equations that we have discarded thanks to the Bianchi identitymay be used to cross-check the accuracy of the numerical integration. Boundary conditions.
In discussing the necessary boundary conditions we will partic-ularize the analysis to a metric which is matched to an external FLRW model at somecomoving coordinate radius r b . We will give boundary conditions at the initial time: theevolution equations will automatically maintain them at later times. Similarly to the embed-ding of an LTB model, the curvature and the integrated gravitating masses have to matchthe FLRW values, E ( r b , ¯ t ) = − k r b and F i ( r b , ¯ t ) = F out i ( r b , ¯ t ), respectively (in the chosengauge the scale function is already matched as Y ( r, ¯ t ) = a (¯ t ) r ). The latter condition forthe case of v i, c ( r, ¯ t ) = 0 explicitly demands that ρ out i (¯ t ) = r b (cid:82) r b ρ i (ˆ r, ¯ t ) ˆ r d ˆ r , as it is easyto see from Eq. (2.34). According to the Darmois-Israel boundary conditions [59] the lo-cal density may be discontinuous through a timelike surface (constant r b ), but the pressuremust be continuos p i ( r b , ¯ t ) = p out i (¯ t ). Furthermore the peculiar velocities have to vanish, v i, c ( r b , ¯ t ) = 0, as also the lapse function, λ ( r b , ¯ t ) = 0, as we want to describe the FLRW– 8 –odel in a geodesic reference frame where proper and coordinate times coincide. Even if notnecessary it may be desirable to match smoothly the background, i.e., ρ i ( r b , ¯ t ) = ρ out i (¯ t ) and ρ (cid:48) i ( r b , ¯ t ) = p (cid:48) i ( r b , ¯ t ) = λ (cid:48) ( r b , ¯ t ) = v (cid:48) i, c ( r b , ¯ t ) = 0.Finally, we discuss the regularity conditions at the center of the inhomogeneous metric.Because we do not want a singularity at the center it is F i (0 , ¯ t ) = E (0 , ¯ t ) = 0. In particular,near the origin it is F i ∝ r and E ∝ r , as one can see from Eq. (2.27) where Y ∝ r and a finite H A (cid:54) = 0 is assumed (see also [37]). The latter also implies ˙ Y ∝ r and so˙ Y (0 , ¯ t ) = Y (0 , ¯ t ) = 0. Also vanishing at the center must be the velocities, v i, c (0 , ¯ t ) = 0, andthe pressure gradient, p (cid:48) i (0 , ¯ t ) = 0, which implies λ (cid:48) (0 , ¯ t ) = 0. It seems again desirable tohave densities and velocities smooth at the origin, ρ (cid:48) i (0 , ¯ t ) = v (cid:48) i, c (0 , ¯ t ) = 0. The metric (2.1) implies that for radial null worldlines it is dtdu = ± drdu e − λ Y (cid:48) √ E , (2.39)where u is an affine parameter. Let us consider a photon emitted at t ( u ) in a coordinatetime ν ( u ) (emitting source and observer are taken to be in the rest frame). Expanding theright-hand side of Eq. (2.39) around t ( u ) and keeping terms up to first order in ν ( u ) we get dtdu + dνdu = − drdu e − λ Y (cid:48) √ E − drdu Y (cid:48) √ E (cid:16) H R − e − λ ˙ λ (cid:17) ν , (2.40)which, using Eq. (2.39) for dt/du , gives1 ν dνdu = − drdu Y (cid:48) √ E (cid:16) H R − e − λ ˙ λ (cid:17) . (2.41)The redshift is defined as z = e λ ν − e λ νe λ ν , (2.42)where the subscript 0 refers to values evaluated at the affine-parameter value u at theobserver’s position. The derivative of z with respect to u is then dzdu = − (1 + z ) drdu (cid:18) dλdu + 1 ν dνdu (cid:19) = − (1 + z ) drdu (cid:18) λ (cid:48) − Y (cid:48) √ E H R (cid:19) . (2.43)Using (2.43) and remembering the expression for the rest-frame acceleration scalar (see Ap-pendix A), Eq. (2.39) gives a pair of equations dtdz = 11 + z e − λ a rf − H R , drdz = −
11 + z a rf − H R √ EY (cid:48) , (2.44)which determines the coordinates t and r on the light cone as functions of the redshift, andcan be put in the following gauge invariant form: dτ rf dz = 11 + z a rf − H R , dr p dz = −
11 + z a rf − H R , (2.45) This definition is different from the one in [5]. The authors of [5] agree, however, that (2.42) is correct. – 9 –here r p is the proper radial distance. In the ΛLTB model where ˙ E = a rf = 0, the previousequations reduce to the familiar form dtdz = − Y (cid:48) (1 + z ) ˙ Y (cid:48) , drdz = √ E (1 + z ) ˙ Y (cid:48) . (2.46) As a first application of the formalism introduced in this paper, we will consider an inhomo-geneous sphere embedded into a flat ΛCDM universe. The inhomogeneities will be given bytwo non-comoving dust components, one of which may be interpreted as baryonic matter andthe other as dark matter. The reader should keep in mind, however, that these results aremerely illustrative of the multi-fluid model here presented. This solution becomes the usualΛLTB model if the density of the second dust component is set to zero or if their relativevelocity is set to zero.We will use the reference frame comoving with the component labelled by “ M ”, so thatthe relative Euler equation gives λ = 0. The other conservation equation reads:˙ ρ M ρ M = − ( H R + 2 H A ) = − ( J R J A )˙ J R J A , (3.1)and can be directly integrated to give: ρ M ( r, t ) ρ M ( r, ¯ t ) = J R ( r, ¯ t ) J A ( r, ¯ t ) J R ( r, t ) J A ( r, t ) = Y ( r, ¯ t ) Y (cid:48) ( r, ¯ t ) (cid:112) E ( r, ¯ t ) (cid:112) E ( r, t ) Y ( r, t ) Y (cid:48) ( r, t ) , (3.2)where ¯ t is the initial time. We will label the other dust component with “ N ” and use thesimpler notation γ N = γ , v N, p = v p , v N, c = v c as the corresponding M quantities are trivial.About Λ note that the corresponding conservation equations (2.35-2.36) trivially state that ρ Λ is uniform and constant. First we set the flat ΛCDM model by fixing the present-day expansion rate to H = 100 h km s − Mpc − with h = 0 . matter = 0 .
3. We thensolve the background equations to find the evolution of the scale factor a ( t ) and so the initialtime ¯ t corresponding to ¯ z = a ( t ) a (¯ t ) −
1, which we will fix illustratively to ¯ z = 5. By redshift z we will always mean the redshift relative to an observer in the background model, which willdiffer, for example, from the redshift relative to the observer at the center of the sphericalinhomogeneity. The inhomogeneous sphere will have a comoving radius of r b = 100 h − Mpc.Next we give initial and boundary conditions at ¯ t . For the scale function we set Y ( r, ¯ t ) = a (¯ t ) r . Then we have to give initial conditions for F ( r, ¯ t ) = ¯ F M ( r ) + ¯ F N ( r ) + ¯ F Λ ( r ). For Λ itis simply ¯ F Λ ( r ) = π a (¯ t ) r ρ Λ . We stress again that the division of the total F into separatecomponents is meaningful only at the initial time as only the total gravitating mass F entersthe evolution equations. Note, however, that it is possible to define individual invariant masscomponents M i . For instance, for the M component it is: M (cid:48) M = 4 πY Y (cid:48) √ E ρ M , (3.3)– 10 –rom which it is easy to verify using Eq. (3.1) that it is correctly ˙ M M = 0 (while ˙ F (cid:54) = 0).The remaining background matter density is split between the M and N componentsaccording to a fraction q such that it is ρ out M (¯ t ) = q ρ outmatter (¯ t ) and ρ out N (¯ t ) = (1 − q ) ρ outmatter (¯ t ).We choose the profile of the M component to be: ρ M ( r, ¯ t ) ρ out M (¯ t ) = 1 + δ ( r, r t M , δ M ,M , α M ,M ) (3.4)where the δ n function is given in Appendix C and δ M = δ M + δ M is the contrast in matterdensity between the center of the spherical inhomogeneity and the background value of ρ out M .This inhomogeneous profile satisfies the smoothness conditions ρ (cid:48) M (0 , ¯ t ) = ρ (cid:48) M ( r b , ¯ t ) = 0 andcontinuously matches the background density ρ M ( r b , ¯ t ) = ρ out M (¯ t ). Moreover, the matchingrequires that ¯ F M ( r ) has the corresponding background value and so the following conditionhas to be satisfied: ρ out M (¯ t ) = 3 r b (cid:90) r b ρ M (ˆ r, ¯ t ) ˆ r d ˆ r , (3.5)from which follows that it has to be δ M δ M < δ M > δ M < δ M , M and α M , M and we will get analytically the radius r t M by demanding Eq. (3.5). The gravitating mass isthen simply ¯ F M ( r ) = 4 πa (¯ t ) (cid:82) r ρ M (ˆ r, ¯ t )ˆ r d ˆ r , which has also an analytic expression.Next we need ¯ F N ( r ), which we will take to be:¯ F N ( r ) = ¯ F out N ( r ) + δ ( r, r t N , δ N ,N , α N ,N ) , (3.6)with δ N = − δ N so that ¯ F N ( r ) is correctly matched to the background value of ¯ F out N ( r ) = π a (¯ t ) r ρ out N (¯ t ).We can now find the initial curvature function E by demanding that the universe has thesame age ¯ t for any r , i.e., from Eq. (2.38): t A ( r, ¯ t ) = ¯ t . The approximation from neglectingthe time dependence of F and E for t < ¯ t is very good as the cosmological constant isnegligible at early times.Next, we have to set the peculiar velocities for the N component with respect to the M component. As there is not pressure, the case v c ( r, ¯ t ) = 0 gives trivially the dynamics of theΛLTB model with a density profile that is the average of the M and N density profiles. Wechoose the initial velocity profile to be: v c ( r, ¯ t ) = δ ( r, r t v , δ v ,v , α v ,v ) (3.7)with δ v = − δ v so that there are no peculiar velocities at the center and at the border. Thecontrast δ v will give the maximum peculiar velocity in c units. Finally we can compute theinitial density profile for the N component, which from Eq. (2.34) is: ρ N ( r, ¯ t ) = ¯ F (cid:48) N ( r )4 π ¯ Y γ ¯ v c (cid:16) ¯ H A ¯ v p + ¯ S A (cid:17) , (3.8)where ¯ H A is given by Eq. (2.32) evaluated at the initial time. This terminates the necessaryinitial conditions. We will now show the exact general relativistic evolution of the system presented in theprevious Section for the case of a central underdensity. The left panel of Fig. 1 shows the The actual parameters we use are q = 0 . α M = α M = 0 . δ M = − . δ M = 0 . r t N = 2 r b / α N = 0 . α N = 0, δ N = − δ N = 0 . F out N ( r b / r t v = 2 r b / α v = α v = 0, δ v = − δ v = 0 . · − /a (¯ t ). – 11 –
20 40 60 80 100020406080100 r Mpc h r bkg Y a Mpc h r bkg Mpc h F F Figure 1 . Left panel: evolution of the scale function Y for the times corresponding to the backgroundredshifts of z =5 (red), 3, 1 and 0 (blue). In particular, this plot shows the relation between thecoordinate radius r and the background comoving radius r bkg = Y /a which will be used in all theother plots. Right panel: Evolution of the gravitating mass F for z =5 (red), 3, 1 and 0 (blue) withthe contribution due to Λ subtracted and normalized to unity at hole radius r b . The time evolution ispurely an effect of the peculiar velocities between the two dust components. See Section 3.2 for moredetails. evolution of the scale function Y : in the interior region the scale function is expanding fasterbecause of the less matter and consequently negative curvature present, as also shown bythe spatial Ricci scalar R in the left panel of Fig. 2 (see Appendix B for its definition). Allthe other figures will be plotted with respect to the background comoving radius defined as r bkg = Y ( r, t ) /a ( t ). This is a good definition as the curvature function is E (cid:28) r p = (cid:82) Y (cid:48) d ˆ r √ E (cid:39) Y .The evolution of the scale function Y is governed by the gravitating mass F and thecurvature function E , as shown by Eq. (2.32). The evolution of the gravitating mass F isplotted in the right panel of Fig. 1. We have subtracted the time-dependent component dueto Λ, which is F Λ = π Y ρ Λ , in order to show the genuine time dependence of the gravitatingmass due to the dust components. As there is no pressure, the evolution displayed in the plotis due to the peculiar velocities between the two components. In particular the gravitatingmass at half the radius is decreasing in accordance with the outgoing flow shown in the rightpanel of Fig. 4.The evolution of the curvature function E , is plotted in the right panel of Fig. 2. Thisis again a genuine effect of the peculiar velocities between the two dust components, as acosmological constant never sources the evolution of the curvature function E . We remindthat E can be interpreted as the total energy per unit of mass of a given shell (see Eq. (2.28)),and one can see from the plot that energy seems to flow outwards, again in agreement withthe peculiar velocity field shown in the right panel of Fig. 4. As said before, the left panelof Fig. 2 shows the Ricci scalar R , which seems to follow a profile which is an average of thedensity contrasts of the two dust components, displayed in Fig. 3. We also remind that thecurvature function is related to the Euclidean average of R , as pointed out in Appendix B.Both the density contrasts of the two dust components show a realistic evolution, as one– 12 –
20 40 60 80 100 r bkg Mpc h a r b r bkg Mpc h E c Figure 2 . Left panel: evolution of the spatial Ricci scalar R plotted as the dimensionless combination R a r b for z =5 (red), 3, 1 and 0 (blue). Its profile is closely related to the density profiles of thetwo dust components displayed in Fig. 3. Right panel: evolution of the curvature/energy function E for z =5 (red), 3, 1 and 0 (blue). The plot shows an outgoing flow of energy in accordance with thepeculiar velocity field shown in the right panel of Fig. 4. See Section 3.2 for more details. r bkg Mpc h Ρ M Ρ M out r bkg Mpc h Ρ N Ρ N out Figure 3 . Evolution of the M component local density (left panel) and N component local density(right panel) for z =5 (red), 3, 1 and 0 (blue). Both the density contrasts show a realistic evolution,with overdense regions contracting and becoming thin shells (mimicking structures), and underdenseregions becoming larger and deeper (mimicking voids). Note in particular how the N componentinduces gravitationally a peak in the density of the M component. See Section 3.2 for more details. can see from Fig. 3. Overdense regions start contracting and they become thin shells (mim-icking structures), while underdense regions become larger and deeper (mimicking voids), andeventually they occupy most of the volume. An interesting feature of this multi-componentmodel is that the gravitational interaction between the two fluids is evident. The profile ofthe M component (left panel) develops indeed a peak exactly where the N component peaks:that is, an overdensity of one components drags in the collapse the other component. This– 13 –
20 40 60 80 1000.00.51.01.52.0 r bkg Mpc h Y H c r bkg Mpc h v p c Figure 4 . Left panel: evolution of the peculiar velocities v M ≡ Y ( H A − H out ) of the M componentwith respect to the background for z =5 (red), 3, 1 and 0 (blue). Velocities are increasing as theunderdensity becomes emptier. Right panel: evolution of the peculiar velocities of the N componentwith respect to the M component for z =5 (red), 3, 1 and 0 (blue). The gravitational attractionbetween the two components reduces their relative velocities. See Section 3.2 for more details. r bkg Mpc h v M v M ,peak c r bkg Mpc h v p v p ,peak c Figure 5 . Same as Fig. 4 but with the velocity of the peak subtracted. These plots show thatthe velocity flow near the peaks in density is natural: matter is falling towards the overdense regionswhich indeed grow with time as illustrated in Fig. 3. See Section 3.2 for more details. can be checked by decreasing the matter density allocated to the N component, with theeffect that the the peak at 90 h − Mpc in the M component disappears.Finally, in Fig. 4 we show the peculiar velocities of the two components. The left panelshows the peculiar velocities of the M component with respect to the background, whichare given by v M ≡ Y ( H A − H out ). This definition again relies on the fact that r p (cid:39) Y .As expected these peculiar velocities are growing with time as the underdensity becomesemptier and emptier (see also the left panel of Fig. 1). Opposite is, however, the case of– 14 –he peculiar velocities of the N component with respect to the M component (right panel ofFig. 4), which decrease with time: the gravitational attraction between the two componentsdumps their relative velocities. The plots of Fig. 4 are with respect to the “reference frame”of the background and of the center of the underdensity, as the velocities go to zero at r = 0and r = r b . In Fig. 5 we re-plot the same curves, but with the velocity of the peak indensity subtracted. Fig. 5 therefore better shows the velocity field close to the overdensity,which appears to be natural, as matter is falling towards the peak in the density from bothdirections.Finally, this model maintains the scale invariance valid for the simpler LTB model [9] ifthe initial peculiar velocities are properly scaled: ˜ v c ( r, ¯ t ) / ˜ r b = v c ( r, ¯ t ) /r b where tilde marksthe quantities relative to an inhomogeneous patch of different radius. As a second application of our general model we will consider an inhomogeneous sphereembedded into a flat w CDM universe. The inhomogeneities will be given by a pressurelessmatter component and by a dark energy source with constant equation of state w out < − / w = w out = − M , so that λ = 0 and Eqs. (3.1-3.2) hold. We will label the dark energy componentwith “ w ” and use the simpler notation w w = w , p w = p , γ w = γ , v w, p = v p , v w, c = v c as the corresponding M quantities are trivial. As we are discussing a fluid with negativepressure we cannot consider adiabatic perturbations for which w ( r, t ) = w out : this wouldindeed give a negative speed of sound ( c s = w out ) which would cause an unstable growth ofperturbations [60]. We will consider instead the following r -dependent non-adiabatic equationof state: w ( r, t ) = w out (cid:20) ρ w ( r, t ) ρ out w ( t ) (cid:21) α − , (4.1)which implies the following speed of sound: c s = ∂p w ∂ρ w = α w ( r, t ) , (4.2)so that for α ≤ α = 1, see e.g. Ref. [61]).Particularly interesting is the case α = 0, for which the pressure gradients are absentand the speed of sound vanishes, with the consequence that matter and dark energy collapsetogether following the geodesic flow. The possibility of a dark energy fluid with negligiblespeed of sound has been investigated in the literature under various assumptions, and gener-ally requires a non canonical scalar field like k -essence, as the standard quintessence modelswith canonical scalar fields always have c s = 1 (see, for example, [49, 50, 62–64] and refer-ences therein). Apart from theoretical considerations, the fact that the speed of sound of darkenergy may vanish opens up new observational consequences. Indeed, the absence of darkenergy pressure gradients allows instabilities to develop on all scales, also on the small scaleswhere dark matter perturbations become non-linear. We expect, therefore, dark energy tomodify in a detectable manner the growth history of dark matter halos not only through itsdifferent background evolution but also by actively participating to the structure formation– 15 –rocess, in the linear and non-linear regime [65–68]. The observational consequences of aclustering dark energy have been indeed extensively studied. See, for example, [69–76] forthe impact on large-scale structures and cosmic microwave background.Finally, we would like to point out that, as we will use the matter reference frame,the speed of sound given in Eq. (4.2) will not be in general in the dark energy rest frame,which is usually used in order to have a gauge independent expression for the speed of sound.However, as we will not introduce initial peculiar velocities, if c s = 0 then the matter anddark energy rest frames will always coincide because of the absence of pressure gradients.The latter will also be approximately true for the case of nonzero but negligible speed ofsound as the peculiar velocities developed by the dark energy fluid will be very small (seethe right panel of Fig. 9). We fix the background model by setting h = 0 . M = 0 . w out = − . α = 0 in Section4.2, w out = − . α = − − in Section 4.3 and w out = − . α = − − in Section4.4. The inhomogeneous sphere will have a comoving radius of r b = 10 h − Mpc and we willgive initial conditions at the time ¯ t corresponding to the background redshift ¯ z = 1000. Weset again the scale function at initial time to Y ( r, ¯ t ) = a (¯ t ) r .We have to give initial conditions for F ( r, ¯ t ) = ¯ F M ( r ) + ¯ F w ( r ). We model the initialmatter inhomogeneity as:¯ F M ( r )¯ F out M ( r ) = 1 + δ M (¯ t ) 1 − tanh(( r − r − r offset ) / r )1 + tanh( r / r ) , (4.3)where ¯ F out M = π a (¯ t ) r ρ out M (¯ t ) is the background gravitating mass, δ M is the density contrast,and the parameters r and ∆ r characterize respectively size and steepness of the densityprofile. The parameter r offset is used to translate the profile in order to have a smooth origin.We will adopt the following values: r = r offset = 0 . r b , ∆ r = 0 . r and δ M (¯ t ) = 1 . · − ,that is, we have a central overdensity. The initial matter density is then: ρ M ( r, ¯ t ) = ¯ F (cid:48) M ( r )4 πa (¯ t ) r . (4.4)Next, we can find the initial curvature function E by demanding that the universe hasthe same age ¯ t for any r , i.e., from Eq. (2.38): t A ( r, ¯ t ) = ¯ t . As ¯ F w ( r ) (cid:28) ¯ F M ( r ) at ¯ t , thiscondition will be independent of the dark energy initial conditions. In particular at theinitial time the model is very close to a (dust) LTB model, so that it is possible to find E analytically by means of the following expression valid for δ M (cid:28) E ( r, ¯ t ) (cid:39) − (cid:16) a (¯ t ) H out (¯ t ) r (cid:17) (cid:18) ¯ F M ( r )¯ F out M ( r ) − (cid:19) . (4.5)These initial conditions are free from decaying modes [78].Finally, we take the dark energy component to have no initial peculiar velocities, whichmeans that it will have, with respect to the background, the peculiar velocities Y ∆ H of thematter component. Therefore, in order to have initial conditions without decaying modes inthe dark energy component, we have to set the initial profile as:¯ F w ( r )¯ F out w ( r ) = 1 + δ w (¯ t ) 1 − tanh(( r − r − r offset ) / r )1 + tanh( r / r ) , (4.6)– 16 – Mpc h r bkg Y a Mpc h r bkg Mpc h Y H c Figure 6 . Left panel: evolution of the scale function Y for the times corresponding to the backgroundredshifts of z = 1000 (red), 100, 10, 2, 1, and 0 (blue). In particular, this plot shows the relationbetween the coordinate radius r and the background comoving radius r bkg = Y /a which will be usedin the other plots. Right panel: evolution of the peculiar velocities v M ≡ Y ( H A − H out ) of the M component with respect to the background. Velocities are directed inwards and increase in magnitudewith time as the central overdensity becomes denser. See Section 4.2 for more details. where ¯ F out w = π a (¯ t ) r ρ out w (¯ t ) and initial contrast is: δ w (¯ t ) = 1 + w out − w out δ M (¯ t ) , (4.7)which is valid during matter domination for sub-Hubble perturbations and c s (cid:28) ρ w ( r, ¯ t ) = ¯ F (cid:48) w ( r )4 πa (¯ t ) r . (4.8)With this we have all the necessary initial conditions. The left panel of Fig. 6 shows the evolution of the scale function Y for w out = − . α = 0:in the interior region the scale function is expanding slower because of the more matter andconsequently positive curvature present. All the other figures will be plotted with respectto the background comoving radius defined as r bkg = Y ( r, t ) /a ( t ), which as explained inSection 3.2 is a good definition of distance. The right panel shows the peculiar velocities ofthe M component with respect to the background. As there are no initial decaying modes,the peculiar velocities are natural, with the matter falling towards the peak in the density,as shown by the left panel of Fig. 6. Moreover, the velocities grow with time as the centraloverdensity becomes denser and denser.In Fig. 7 the evolution of the density contrasts of the matter and dark energy componentsis shown. The evolution, as remarked above, is realistic. In particular the dark energycomponent closely follows the evolution of the matter component, as one should expect froma fluid with vanishing speed of sound. Because of the absence of pressure gradients peculiar– 17 – r bkg Mpc h Ρ M Ρ M out r bkg Mpc h Ρ w Ρ w out Figure 7 . Evolution of the matter contrast (left panel) and dark-energy contrast (right panel) for z = 1000 (red), 100, 10, 2, 1, and 0 (blue). Both the density contrasts show a realistic evolution,with the overdense region becoming denser and denser. Note, in particular, how the dark energycomponent follows the evolution of the matter component, with the difference of a slower growth ratedue to the nonzero pressure. See Section 4.2 for more details. velocities between the two components never develop (if initially zero) and the two fluidsboth follow the same geodesic flow. As it is clear from the plot, however, the growth ratesare different. This has to be expected because the negative pressure generally slows down theevolution of perturbations. To be quantitative, we have plotted in right panel of Fig. 8 theevolution of the ratio of the contrasts, δ w /δ M , which we normalized to the value valid duringmatter domination, w out − w out . As one can see, the agreement is perfect when Ω M ∼ r bkg = 1 h − Mpc (overdensity) and the green line to r bkg = 6 h − Mpc (underdensity).This could be due to the fact that even though the two components have similar profiles, inthe overdense region the matter component reaches the nonlinear evolution earlier (the darkenergy is actually almost linear at the present time).It is also interesting to compare the evolution to the case of no dark energy. We findthat the matter perturbation grows to a maximum contrast of about ∼
90. Alternatively, wefind that an initial perturbation of δ M (¯ t ) = 1 . · − grows to the same present-day contrastas the initial perturbation of δ M (¯ t ) = 1 . · − studied in this Section.Finally, let us comment on the fact that because the speed of sound is zero, there areno characteristic length scales associated to the dark energy clustering and the sphericalcollapse remains independent of the size of the object [65]. The characteristic length scaleassociated to the dark energy clustering is indeed the sound horizon scale, L s = a (cid:82) c s dt/a ,which vanishes for c s = 0 so that clustering takes place on all scales. In other words, thismodel maintains the scale invariance valid for the simpler LTB model [9]. For the same reason the curvature is time independent, as it is easy to see from Eq. (2.33). – 18 – .01 0.1 1 10 100 10000.00.20.40.60.81.0 z M blue and w red Figure 8 . Left panel: redshift evolution of the matter and dark energy density parameters. Rightpanel: evolution of the ratio of the contrasts, δ w /δ M , normalized to the value valid during matterdomination, w out − w out . The blue (dot-dashed) line is relative to r bkg = 1 h − Mpc and the green (solid)line to r bkg = 6 h − Mpc, see Fig. 7. See Section 4.2 for more details. r bkg Mpc h Α c s r bkg Mpc h v p c Figure 9 . Left panel: evolution of the speed of sound c s for the times corresponding to thebackground redshifts of z = 1000 (red), 100, 10, 2, 1, and 0 (blue). Right panel: evolution of thepeculiar velocities of the dark-energy component with respect to the matter component. See Section4.3 for more details. c s ∼ − We will now consider the model of the previous Section ( w out = − .
8) with α = − − ,which corresponds to a speed of sound of order c s ∼ − , as one can see from the left panelof Fig. 9. This latter value implies a sound horizon of order L s ∼ h − kpc, much less thanthe inhomogeneity scale considered in this example, so that the overall dynamical evolutionshown in Figs. 6-8 is basically unchanged. The new features of this case are pressure gradientsand peculiar velocities. In order to understand their behavior it is useful to work out the– 19 – r bkg Mpc h Ρ M Ρ M out r bkg Mpc h Ρ w Ρ w out Figure 10 . Evolution of the matter contrast (left panel) and dark-energy contrast (right panel) for z = 1000 (red), 100, 10, 2, 1, and 0 (blue). Note that the evolution of the dark energy componentstill follows the matter flow but develops an underdensity instead of an overdensity. See Section 4.4for more details. following approximations valid for | δ w | , | α | (cid:28) c s (cid:39) α w out (1 − δ w ) , δ p ≡ p − p out | p out | (cid:39) − α δ w . (4.9)From the previous relations it is easy to understand how speed of sound and pressure gradientsare related to the density contrast. In particular it is δ (cid:48) p ∼ c (cid:48) s ∼ − δ (cid:48) w , as one can also seeby comparing the left panel of Fig. 9 with the right panel of Fig. 7. As a consequence, thepressure gradients will push the dark energy component out of the free-falling geodesic awayfrom the overdense regions, as shown in the right panel of Fig. 9. One can indeed obtain atlinear order (see, e.g., [80]) the approximate relation ˙ v p ∼ − c s w out δ (cid:48) w . w < − Finally, we will now consider the case of w out = − . α = − − . We chose the sameorder of magnitude for the (positive) speed of sound of the previous Section for comparisonsake, see the left panel of Fig. 11. Clearly, we can decrease the magnitude of α in order toprevent possible pathologies [62].As one can see from the left panel of Fig. 10, the evolution of the matter component issimilar to the previous case, with a moderately higher contrast at late times probably due tothe different expansion history: for w out = − . w out = − . Quite different is instead the evolution ofthe dark energy component, which still follows the matter flow but develops an underdensityrather than an overdensity (see the right panel of Fig. 10). This can be understood fromthe fact that the relation δ w /δ M ∼ w out − w out changes sign for w out < −
1. It is interestingto note that a similar behavior was found in the analysis of Ref. [66] of perturbations inscalar-tensor cosmologies, where the non-minimal coupling of the field is indeed responsible Similar considerations apply to the case of a cosmological constant which gives a present-day contrast ofabout ∼
8, roughly half way between the cases studied in Sections 4.2 and 4.4. – 20 – r bkg Mpc h Α c s r bkg Mpc h v p c Figure 11 . Left panel: evolution of the speed of sound c s for the times corresponding to thebackground redshifts of z = 1000 (red), 100, 10, 2, 1, and 0 (blue). Right panel: evolution of thepeculiar velocities of the dark-energy component with respect to the matter component. See Section4.4 for more details. for the crossing of the so-called phantom divide ( w out = − v p ∼ − c s w out δ (cid:48) w both δ (cid:48) w and 1 + w out changesign. In this Section we will compare our findings to previous work dealing with exact solutions. Asexplained in Section 1 the model here presented is a generalization to the case of n decoupledand non-comoving perfect fluids of the Lemaˆıtre model [1] (see [2–5] for recent contributions)which describes the dynamics of a spherically-symmetric perfect fluid. Therefore, the model ofEqs. (2.32-2.36) straightforwardly reduces to the Lemaˆıtre model, possibly with the additionof a cosmological constant.In the case of one source only, our equations can be used to describe the fluid in non-comoving possibly-geodesic coordinates, which may be desirable if the fluid has pressuregradients and its four-velocity undergoes acceleration. Within our formalism this is accom-plished by simply taking a geodesic reference-frame four-velocity u α rf , i.e., by setting λ (cid:48) = 0(see Eq. (A.18)). This can be useful to disentangle the effect of pressure gradients on theevolution of the model as shown, for instance, by Eq. (2.17) for the time evolution of thecurvature function. In the latter equation, indeed, only the second term on the right-handside contributes if comoving coordinates are used and only the first term on the right-handside contributes if geodesic coordinates are used. Spherical models have been discussed pre-viously in the literature in non-comoving coordinates, albeit with a different purpose. Themain idea was indeed that an exact solution with a simple appearance in non-comoving co-ordinates may become extremely complex when transformed to the appropriate comovingsystem, or that the integration of the partial differential equations necessary to obtain thecomoving coordinates may even represent an intractable mathematical problem. We referthe interested reader to [81–87] and to [88] and references therein.– 21 –nother line of research focused on interpreting two or more fluid components as a singleeffective fluid, see for example [89] and references therein. The work of [90, 91] showed indeedthat a mixture of fluids whose four-velocities lay on a two-plane can always be reinterpreted asa single anisotropic fluid. The model of this paper satisfies this last requirement and it wouldbe interesting to investigate in future work the amount of anisotropic pressure generated bythe evolution of a cosmological multi-component fluid modeling, e.g., baryons, dark matterand dark energy. As a simple example we can apply here the work of [90] to the model ofSection 3 featuring two non-comoving dust components in a flat ΛCDM universe. Following[90], the energy-momentum tensor of the effective fluid is: T αβ eff = ρ eff u α eff u β eff + S αβ eff , (5.1)and a short calculation shows that the rest energy density is: ρ eff = 12 ( ρ M + ρ N ) + 12 (cid:112) ( ρ M + ρ N ) + 4 ρ M ρ N ( γ − , (5.2)and that the components of the diagonalized stress tensor are: p A = 0 , (5.3) p R = −
12 ( ρ M + ρ N ) + 12 (cid:112) ( ρ M + ρ N ) + 4 ρ M ρ N ( γ − , (5.4)where p A and p R are the angular and radial pressures, respectively. Clearly, for γ → v c →
0) the effective fluid becomes simply the sum of the two fluids and p A = p R = 0. Fora nonzero radial velocity v c (cid:54) = 0, however, the radial pressure is positive p R > We have extended the Lemaˆıtre model to the case of n decoupled and non-comoving perfectfluids with general equations of state. We have expressed the full set of 3 + 2 n exact equa-tions governing the dynamics of the model in a concise and transparent way thanks to the useof physically meaningful quantities like, for example, expansion and acceleration rates, ex-pansion and acceleration scalars, total derivatives along four-velocity and acceleration. Thisgeneral solution can have many possible applications, both at early and late times.At late time, one can study the evolution of overdensities and underdensities in a uni-verse where dark energy is not the cosmological constant and, in particular, can be inho-mogeneous. This possibility may also induce an inhomogeneous variation of fundamental– 22 –onstants if a coupling between the dark energy and the matter-radiation Lagrangian is al-lowed. Moreover, dark matter and baryons can be described as two separate fluids, the latterpossibly featuring pressure. Pressure in general can have a non-negligible effect on the cos-mological models and this should be taken into consideration while interpreting cosmologicaldatasets. At early times the contribution of radiation can be included, which may be relevantfor the understanding of the evolution of the inhomogeneities.In this paper we focused on the formal understanding of the general solution, whichwe applied to two different setups. For the case of two non-comoving dust components in aflat ΛCDM universe we found rich and interesting physics as, for instance, the gravitationalinteraction between the two fluids which makes one component follow the collapse of theother, or the evolution of the peculiar velocities which are increasing with respect to thebackground but decreasing between the two fluids.Finally, for the case of clustering dark energy, we have been able to follow in an exact waythe collapse of an overdensity to which also the dark energy contributes, and we confirmedthat the evolution agrees with the theoretical growth rates during matter domination. Wealso considered the case of a small speed of sound ( c s ∼ − ) and of a phantom equationof state, which also showed the expected dynamics. These applications show some of theinteresting features of this n -fluid component solution, which we will further investigate inforthcoming work, both from a theoretical and observational point of view. Acknowledgments
It is a pleasure to thank Krzysztof Bolejko, Marie-No¨elle C´el´erier, Woei Chet Lim, LeandrosPerivolaropoulos, Ignacy Sawicki and Wessel Valkenburg for useful comments and discussions.M. P. acknowledges financial support from the Magnus Ehrnrooth Foundation.
A Expansion tensor
In this Appendix we will study the kinematical properties of a generic four-velocity field u α ,which may be identified with the i :th fluid-component velocity u αi or with the reference-framevelocity u α rf .The acceleration of u α is given by a α ≡ u µ ∇ µ u α , whose components are: a t = − e λ γ v p a and a r = γ v p v c a , (A.1)where the acceleration scalar a is a = (cid:113) g αβ a α a β = v − p (cid:16) Θ T + v p Θ R (cid:17) , (A.2)where Θ T and Θ R are defined below. The expansion tensor Θ αβ is then:Θ αβ = h µα h νβ ∇ ( µ u ν ) = ∇ ( α u β ) + a ( α u β ) , (A.3)where h αβ = g αβ + u α u β is the projection tensor on the hypersurface orthogonal to u α .In the second equality we have used the fact that u α ∇ β u α = 0. As we are considering aspherically-symmetric spacetime the velocity is irrotational, i.e., there is no vorticity: ω αβ = h µα h νβ ∇ [ µ u ν ] = 0 . (A.4)– 23 – αβ is often referred to as the “spatial metric” because it is the induced metric on the spaceslices, h αβ = h µα h νβ g µν . Moreover, h αβ is the first fundamental form of the hypersurface andthe expansion tensor is (minus) the extrinsic curvature, K αβ = − Θ αβ , which is the secondfundamental form of the hypersurface.The expansion scalar can be obtained without calculating explicitly the components ofthe expansion tensor by remembering that it is:Θ = ∇ α u α = ∂ α ( J u α ) J = ∂ α u α + d ln Jdτ , (A.5)where J = √− g = J T J R J A , g is the determinant of the metric and we have defined thefollowing quantities: J T = e λ , J R = Y (cid:48) √ E and J A = Y sin θ . (A.6)In the last expression of Eq. (A.5), the first term can be interpreted as the “newtonian”expansion scalar and the second as the contribution from the metric. We then define thetemporal, radial and angular expansion scalars:Θ T = ∂ t u t + d ln J T dτ = e − λ ˙ γ + γv c λ (cid:48) , (A.7)Θ R = ∂ r u r + d ln J R dτ = ( γv c ) (cid:48) + γ (cid:18) H R + S R (cid:19) , (A.8)Θ A = ∂ θ u θ + ∂ φ u φ + d ln J A dτ = γ (cid:18) H A + 2 S A (cid:19) , (A.9)respectively, where we defined the quantities: S R = v c Y (cid:48)(cid:48) Y (cid:48) − v c E (cid:48) E and S A = v c Y (cid:48) Y , (A.10)which are defined analogously to H R and H A of Eqs. (2.2-2.3) and can be understood as“spatial” expansion rates relevant for non-comoving fluids. From Eq. (A.2) follows that theexpansion scalar for a geodesic four-velocity field is Θ geod . = Θ R /γ + Θ A .In order to calculate the shear we need the non-vanishing components of Θ αβ :Θ tt = e λ γ v p (Θ T + Θ R ) , Θ rt = − e λ γ v p v c (Θ T + Θ R ) , (A.11)Θ rr = γ v p v c (Θ T + Θ R ) , Θ θθ = Θ φφ / sin θ = 12 Y Θ A . (A.12)Note that Θ tt and Θ rt are nonzero as we are using coordinates in general not comoving with u α . The components of the expansion tensor clearly satisfy:Θ = g αβ Θ αβ = Θ T + Θ R + Θ A . (A.13)The shear tensor σ αβ is then defined in terms of Θ αβ as σ αβ = Θ αβ −
13 Θ h αβ , (A.14)– 24 –nd its non-vanishing components are: σ tt = 2 √ e λ γ v p σ , σ rt = − √ e λ γ v p v c σ , (A.15) σ rr = 2 √ γ v p v c σ , σ θθ = σ φφ sin θ = − √ Y σ , (A.16)where the invariant σ of the shear tensor is given by: σ = 12 g αµ g νβ σ µν σ αβ = 13 (cid:18) Θ T + Θ R −
12 Θ A (cid:19) . (A.17)We may now particularize the above expressions for the case of the reference-framevelocity u α rf : a rf = v c v p λ (cid:48) , a r, rf = λ (cid:48) , a t, rf = 0 , (A.18)Θ rf = Θ R, rf + Θ A, rf = H R + 2 H A , √ σ rf = Θ R, rf −
12 Θ A, rf = H R − H A . (A.19)Note that, as Θ T vanishes, λ appears in the expansion scalar and shear only through thederivatives with respect to the proper time d/dτ rf = u α rf ∂ α = e − λ ∂/∂t . This shows that anon-geodesic reference frame only changes the proper time with which the expansion rateis measured, but not its functional form (different is the case for the acceleration rates ofEqs. (2.4-2.5)). Moreover, these last results show that the definitions in Eqs. (2.2-2.3) arerelative to the reference-frame velocity. B The Lemaˆıtre metric
It is useful to consider the case of only one fluid component, which could be seen as theLemaˆıtre metric [1] in a possibly non-comoving reference frame. We will discuss some prop-erties of this solution, which may help us better understand the multi-fluid non-comovingcase. We will not discuss the dynamical evolution of the Lemaˆıtre metric, which can be easilydeduced from the set of Eqs. (2.32-2.36).Let us start by considering the Eqs. (2.30-2.31) for ˙ F and F (cid:48) . By taking the combinationsgiving the derivative of F with respect to τ and σ it is possible to obtain the following simplerexpressions: dFdτ = − p dV e dτ = − πY p Θ A , (B.1) dFdσ = ρ dV e dσ , (B.2)where the Euclidean and actual volumes of the space slices are, respectively: V e = 4 π (cid:90) r Y Y (cid:48) d ˆ r = 4 π Y , V = 4 π (cid:90) r Y Y (cid:48) √ E d ˆ r . (B.3)Note that Eqs. (B.1-B.2) give the gravitating mass of a non-comoving fluid. In the fluid restframe (the actual Lemaˆıtre metric) the previous equations have the familiar form:˙ F = − πY ˙ Y p = − p ˙ V e , (B.4) F (cid:48) = 4 πY Y (cid:48) ρ = ρ V (cid:48) e , (B.5)– 25 –o that the first gives the time evolution of F , which is constant only for vanishing pressure,and the second links the local density to the integrated gravitating mass. Eqs. (B.4-B.5) canbe combined into the differential expression: dF = ρ dV e | t =const − p dV e | r =const , (B.6)which has the usual thermodynamical interpretation for a fluid in equilibrium (at constantentropy): for an increase dV of volume at constant t the total mass-energy is increased by themass-energy density of the volume dV , and for an increase dV at constant r the total mass-energy is decreased by the work of the pressure against dV . In this picture the coordinate r has to be seen as a label for the mass-energy density.It is interesting to point out that the previous equations use the Euclidean volumeelement and so the gravitating (or Euclidean) mass F does not coincide with the invariantmass M which is related to the local density by [7]: M (cid:48) = 4 π Y Y (cid:48) √ E ρ = ρ V (cid:48) . (B.7)We can then evaluate the time derivative of M at constant r :˙ M = 4 π (cid:90) r Y Y (cid:48) √ E ρ ( J R J A ρ )˙ J R J A ρ d ˆ r = 4 π (cid:90) r Y Y (cid:48) √ E ρ (cid:18) ˙ ρρ + e λ Θ rf (cid:19) d ˆ r = − π (cid:90) r Y Y (cid:48) √ E e λ Θ rf p d ˆ r = − p ˙ V , (B.8)where to go from the first to the second line we have used the conservation equation (2.35)and for the last equality we have assumed a homogeneous pressure p (cid:48) = 0. Note that if p (cid:48) = 0,then from Eqs. (2.36) and (2.33) follows that ˙ E = 0, even if M and F are evolving. Similarlyto Eq. (B.6), we can combine Eqs. (B.7-B.8) into: dM = ρ dV | t =const − p dV | r =const . (B.9)Note, however, that in this last equation not only a different volume is used with respectto Eq. (B.6), but that Eq. (B.9) is not valid for p (cid:48) (cid:54) = 0, which implies a time-dependentcurvature (that can be interpreted as total energy per unit of mass). In this latter case theusual thermodynamical expression, where only the pressure at the boundary matters, doesnot seem to hold and the value of the pressure at any r matters as shown by the last integralin Eq. (B.8).Finally, it is interesting to rewrite Eq. (2.32) as: H A = 2 GFY + 2 EY = κ (cid:104) ρ (cid:105) e − (cid:28) R (cid:29) e , (B.10)where the Euclidean average for the generic quantity X is defined as: (cid:104) X (cid:105) e = 4 π (cid:82) r Y Y (cid:48) X d ˆ rV e , (B.11)and R is the spatial Ricci scalar of the rest frame, which is the trace of the Ricci tensor ofthe metric g ij on the hypersurface of constant t : R = − EY ) (cid:48) Y Y (cid:48) . (B.12)– 26 –he factor of 1 / R / k/a . From Eq. (B.10) follows that the expansion rate is sourced by the Euclideanaverage of the local density and curvature, and not by the actual averages. This is a potentialsource for “strong” backreaction effects [96, 97] of the inhomogeneities on the evolution ofthe background (see, for example, Appendix B of [40] and also [98, 99]). C The W n and δ n functions The W n ( x, α ) step functions introduced in [48] are C n everywhere and go from 1 to 0 for x from 0 to 1 with a transition which is steeper for higher α ∈ [0 , W and W , which have the following explicit form: W ( x, α ) = ≤ x < α + sin (cid:104) π (cid:16) − x − α − α (cid:17)(cid:105) α ≤ x <
10 1 ≤ x , and W ( x, α ) = ≤ x < α π (cid:26) π (cid:20) − (cid:16) x − α − α (cid:17) (cid:21) − cos (cid:16) π x − α − α (cid:17)(cid:27) α ≤ x < α π (cid:20) − π (cid:16) x − α − α − (cid:17) + cos (cid:16) π x − α − α (cid:17)(cid:21) α ≤ x <
10 1 ≤ x From the W n function we then built the following general contrast δ n : δ n ( r, r t , ∆ , , α , ) = ∆ + ∆ W n (cid:18) rr t , α (cid:19) ≤ r < r t ∆ W n (cid:18) r − r t r b − r t , α (cid:19) r t ≤ r < r b r b ≤ r , (C.1)where ∆ = ∆ + ∆ is the contrast between the center ( r = 0) of the spherical inhomo-geneity and the border ( r = r b ), and ∆ is the contrast at r t . We have used δ n to modelinhomogeneities in ρ M , F N and v c in Section 3. References [1] G. Lemaitre, Annales Soc. Sci. Brux. Ser. I Sci. Math. Astron. Phys. A , 51 (1933).[2] K. Bolejko, Mon. Not. Roy. Astron. Soc. , 924 (2006)[3] P. D. Lasky and A. W. C. Lun, Phys. Rev. D , 084013 (2006).[4] A. A. H. Alfedeel, C. Hellaby, Gen. Rel. Grav. , 1935-1952 (2010).[5] P. D. Lasky and K. Bolejko, Class. Quant. Grav. , 035011 (2010).[6] R. C. Tolman, Proc. Nat. Acad. Sci. , 169 (1934).[7] H. Bondi, Mon. Not. Roy. Astron. Soc. , 410 (1947).[8] T. Biswas and A. Notari, JCAP , 021 (2008).[9] V. Marra, E. W. Kolb, S. Matarrese and A. Riotto, Phys. Rev. D , 123004 (2007). – 27 –
10] V. Marra, E. W. Kolb and S. Matarrese, Phys. Rev. D , 023003 (2008).[11] N. Brouzakis, N. Tetradis and E. Tzavara, JCAP , 008 (2008).[12] V. Marra, Padua@research ID588, arXiv:0803.3152 [astro-ph].[13] R. A. Vanderveld, E. E. Flanagan and I. Wasserman, Phys. Rev. D , 083511 (2008).[14] W. Valkenburg, JCAP , 010 (2009).[15] S. J. Szybka, Phys. Rev. D84 , 044011 (2011).[16] M. N. Celerier, Astron. Astrophys. , 63 (2000).[17] K. Tomita, Astrophys. J. , 38 (2000).[18] J. W. Moffat and D. C. Tatarski, Astrophys. J. , 17 (1995).[19] H. Alnes, M. Amarzguioui and O. Gron, Phys. Rev. D , 083519 (2006).[20] D. J. H. Chung and A. E. Romano, Phys. Rev. D , 103507 (2006).[21] K. Enqvist and T. Mattsson, JCAP , 019 (2007).[22] M. Tanimoto and Y. Nambu, Class. Quant. Grav. , 3843 (2007).[23] S. Alexander, T. Biswas, A. Notari and D. Vaid, JCAP , 025 (2009).[24] J. Garcia-Bellido, T. Haugboelle, JCAP , 003 (2008).[25] C. M. Yoo, T. Kai and K. i. Nakao, Prog. Theor. Phys. , 937 (2008).[26] J. P. Zibin, A. Moss and D. Scott, Phys. Rev. Lett. , 251303 (2008).[27] K. Kainulainen and V. Marra, Phys. Rev. D , 127301 (2009)[28] S. February, J. Larena, M. Smith and C. Clarkson, Mon. Not. Roy. Astron. Soc. , 2231(2010).[29] J. Sollerman et al. , Astrophys. J. , 1374 (2009).[30] E. W. Kolb and C. R. Lamb, arXiv:0911.3852 [astro-ph.CO].[31] P. Dunsby, N. Goheer, B. Osano and J. P. Uzan, JCAP , 017 (2010).[32] C. M. Yoo, K. i. Nakao and M. Sasaki, JCAP , 012 (2010).[33] T. Biswas, A. Notari, W. Valkenburg, JCAP , 030 (2010).[34] C. Clarkson, M. Regis, JCAP , 013 (2011).[35] A. Moss, J. P. Zibin, D. Scott, Phys. Rev. D83 , 103515 (2011).[36] V. Marra, M. Paakkonen, JCAP , 021 (2010).[37] N. Mustapha, C. Hellaby and G. F. R. Ellis, Mon. Not. Roy. Astron. Soc. , 817 (1997).[38] H. Iguchi, T. Nakamura and K. i. Nakao, Prog. Theor. Phys. , 809 (2002).[39] M. N. Celerier, K. Bolejko and A. Krasinski, Astron. Astrophys. , A21 (2010).[40] V. Marra, A. Notari, Class. Quant. Grav. , 164004 (2011).[41] K. Bolejko, M. N. Celerier and A. Krasinski, Class. Quant. Grav. , 164002 (2011).[42] R. K. Sheth and R. van de Weygaert, Mon. Not. Roy. Astron. Soc. , 517 (2004).[43] K. Bolejko, A. Krasinski and C. Hellaby, Mon. Not. Roy. Astron. Soc. , 213 (2005).[44] P. S. Joshi, I. H. Dwivedi, Phys. Rev. D47 , 5357-5369 (1993).[45] A. Krasinski and C. Hellaby, Phys. Rev. D , 043502 (2004).[46] J. T. Firouzjaee and R. Mansouri, Gen. Rel. Grav. , 2431 (2010). – 28 –
47] J. P. Mimoso, M. Le Delliou and F. C. Mena, Phys. Rev. D , 123514 (2010)[48] W. Valkenburg, arXiv:1104.1082 [gr-qc].[49] D. Bertacca, N. Bartolo, S. Matarrese, Adv. Astron. , 904379 (2010).[50] M. Li, X.-D. Li, S. Wang, Y. Wang, Commun. Theor. Phys. , 525 (2011).[51] J. D. Bekenstein, Phys. Rev. D (1982) 1527.[52] K. A. Olive and M. Pospelov, Phys. Rev. D (2002) 085044.[53] V. Marra and F. Rosati, JCAP , 011 (2005).[54] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D , 1753 (2006).[55] J. K. Webb, J. A. King, M. T. Murphy, V. V. Flambaum, R. F. Carswell andM. B. Bainbridge, Phys. Rev. Lett. , 191101 (2011).[56] J. Grande, L. Perivolaropoulos, Phys. Rev. D84 , 023514 (2011).[57] M. E. Cahill, G. C. McVittie, J. Math. Phys. , 1382 (1970).[58] M. E. Cahill, G. C. McVittie, J. Math. Phys. , 1392 (1970).[59] Junction Conditions in General Relativity, course notes by Charles Hellaby,[email protected][60] W. Hu, Phys. Rev. D71 , 047301 (2005).[61] E. W. Kolb, M. S. Turner, “The Early universe,” Front. Phys. , 1-547 (1990).[62] P. Creminelli, G. D’Amico, J. Norena, F. Vernizzi, JCAP , 018 (2009).[63] D. Bertacca, N. Bartolo, A. Diaferio, S. Matarrese, JCAP , 023 (2008).[64] E. A. Lim, I. Sawicki, A. Vikman, JCAP f 1005, 012 (2010).[65] P. Creminelli, G. D’Amico, J. Norena, L. Senatore, F. Vernizzi, JCAP f 1003, 027 (2010).[66] J. C. B. Sanchez, L. Perivolaropoulos, Phys. Rev. f D81, 103505 (2010).[67] E. Sefusatti, F. Vernizzi, JCAP f 1103, 047 (2011).[68] S. Anselmi, G. Ballesteros and M. Pietroni, JCAP , 014 (2011).[69] S. DeDeo, R. R. Caldwell and P. J. Steinhardt, Phys. Rev. D , 103509 (2003). [Erratum-ibid.D , 129902 (2004)].[70] L. Amendola, F. Finelli, C. Burigana, D. Carturan, JCAP , 005 (2003).[71] J. Weller and A. M. Lewis, Mon. Not. Roy. Astron. Soc. , 987 (2003).[72] R. Bean, O. Dore, Phys. Rev. D69 , 083503 (2004).[73] W. Hu and R. Scranton, Phys. Rev. D , 123002 (2004).[74] S. Hannestad, Phys. Rev. D , 103519 (2005).[75] P. S. Corasaniti, T. Giannantonio and A. Melchiorri, Phys. Rev. D , 123521 (2005).[76] D. Bertacca, A. Raccanelli, O. F. Piattella, D. Pietrobon, N. Bartolo, S. Matarrese,T. Giannantonio, JCAP , 039 (2011).[77] K. Van Acoleyen, JCAP , 028 (2008).[78] J. P. Zibin, Phys. Rev. D78 , 043504 (2008).[79] G. Ballesteros, J. Lesgourgues, JCAP , 014 (2010).[80] O. E. Bjaelde and Y. Y. Y. Wong, JCAP , 038 (2011).[81] McVittie, G. C., & Wiltshire, R. J. 1977, International Journal of Theoretical Physics, 16, 121. – 29 –
82] Bonnor, W. B., & Knutsen, H. 1993, International Journal of Theoretical Physics, 32, 1061.[83] H. Knutsen, Class. Quant. Grav. , 2817 (1995).[84] L. Herrera, W. Barreto, A. Di Prisco and N. O. Santos, Phys. Rev. D , 104004 (2002).[85] Davidson, W. 2003, Classical and Quantum Gravity, 20, 1835[86] M. Ishak, Phys. Rev. D , 124027 (2004).[87] L. Herrera, A. Di Prisco and J. Ibanez, Phys. Rev. D , 064036 (2011).[88] H. Stephani, D. Kramer, M. A. H. MacCallum, C. Hoenselaers, E. Herlt, “Exact solutions ofEinstein’s field equations,” Cambridge, UK: Univ. Pr. (2003) 701 P.[89] A. Krasinski, “Inhomogeneous cosmological models,” Cambridge University Press, Cambridge1997, 317 pp.[90] P. S. Letelier, Phys. Rev. D22 , 807 (1980).[91] P. S. Letelier, P. S. C. Alencar, Phys. Rev.
D34 , 343 (1986).[92] J. P. de Le´on, Phys. Rev.
D35 , 2060 (1987).[93] R. L. Bowers, E. P. T. Liang, Astrophys. J. , 657 (1974).[94] M. Cosenza, L. Herrera, M. Esculpi, L. Witten, Phys. Rev.
D25 , 2527 (1982).[95] R. A. Sussman and D. Pavon, Phys. Rev. D (1999) 104023.[96] E. W. Kolb, V. Marra, S. Matarrese, Phys. Rev. D78 , 103002 (2008).[97] E. W. Kolb, V. Marra and S. Matarrese, Gen. Rel. Grav. , 1399 (2010).[98] M. Mattsson, T. Mattsson, JCAP , 021 (2010).[99] R. A. Sussman, Class. Quant. Grav. , 235002 (2011)., 235002 (2011).