{'fraction_non_alphanumeric': 0.040231468724166435, 'fraction_numerical': 0.014151871826162263, 'mean_word_length': 4.2232959802611285, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact spherically-symmetric inhomogeneous model with n perfect fluids

We present the exact equations governing the dynamics of a spherically-symmetric inhomogeneous model with n decoupled and non-comoving perfect fluids. Thanks to the use of physically meaningful quantities we write the set of 3+2n equations in a concise and transparent way. The n perfect fluids can have general equations of state, thus making the model extremely flexible to study a large variety of cosmological and astrophysical problems. As applications we consider a model sourced by two non-comoving dust components and a cosmological constant, and a model featuring dust and a dark energy component with negligible speed of sound.


Introduction
Spherically-symmetric models are of interest in the study of the nonlinear inhomogeneities of the universe for two reasons. First, spherical symmetry is often a reasonable working approximation for the inhomogeneities in the real universe. Second, spherical symmetry allows to exactly solve the Einstein's equations, a highly nontrivial task. Recent research has focused mainly on the Lemaître model [1] (see [2][3][4][5] for recent contributions) which describes the dynamics of a spherically-symmetric perfect fluid, and great attention has received 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 been used to describe the large-scale inhomogeneities of the late universe as, for example, in Swisscheese [8][9][10][11][12][13][14][15], void  and inhomogeneous [37][38][39] models for apparent acceleration (see [40,41] for recent reviews). Besides being useful in understanding the role of voids in the process of structure formation [42,43], the LTB model is also suitable to study the dynamics of spherical collapse in an expanding universe, black-hole formation included [44][45][46][47][48].
In the present paper we extend the Lemaître model to the case of n decoupled and non-comoving perfect fluids with general equations of state. Our main result is the set of 3 + 2n exact equations governing the dynamics of the model, which we write in terms of physically meaningful quantities and express in a concise and transparent way. Many are the possible applications, both at early and late times. At late time, one can study the evolution of overdensities and underdensities in a universe where dark energy is not the cosmological constant and, in particular, can be inhomogeneous (see, for example, [49,50] and references therein). This possibility may also induce an inhomogeneous variation of fundamental constants, as for instance the fine structure constant if a coupling between the dark energy and the electromagnetic field is allowed (see, for example, [51][52][53][54][55]). Moreover, dark matter and baryons can be described as two separate fluids, the latter possibly featuring pressure. As stressed in [5], the introduction of non-dust equations of state can have a nonnegligible effect on the cosmological models and this should be taken into consideration while interpreting cosmological datasets. This is, however, a non-trivial task due to the increased number of free parameters, and such difficulties increase when one considers n fluid components with the consequence that more astronomical and astrophysical input is needed to constrain the available parameter space. Finally, at early times the contribution of radiation can be included, which may be relevant both for the understanding of the evolution of the standard post-inflation inhomogeneities and for the correct modeling of the very large underdensities typical of void models [34].
This paper is organized as follows. In Section 2 we will go through all the details of our formalism, and in Sections 3 and 4 we will present the numerical results for the case of two non-comoving dust components in a flat ΛCDM universe and for spherical collapse in the presence of dark energy with negligible speed of sound, respectively. Then in Section 5 we will compare our findings to previous work dealing with exact solutions. Finally we will give our conclusions in Section 6. In Appendix A we discuss the expansion tensor and in Appendix B the case of the Lemaître metric, possibly in a non-comoving frame. The definitions of the functions used to model the numerical example of Section 3 are given in Appendix C.

Metric and Einstein tensor
A spherically-symmetric metric can be written as where the lapse function e λ , the scale function Y and the curvature function E depend on the coordinate time t and coordinate radius r, dΩ 2 = dθ 2 + sin 2 θ dφ 2 and we have set c = 1. A prime denotes partial derivation with respect to the coordinate radius r, whereas a dot denotes partial derivation with respect to the coordinate time t. Comma and semicolon signs will not denote derivatives. We will use a reference frame comoving with the arbitrary four-velocity field u α rf , which has components u α rf = (e −λ , 0, 0, 0). Moreover, we assume that Y = 0 and E > −1/2 in order to always have a regular g rr > 0. We will discuss more general evolutions and shell crossing in forthcoming work.
The nontrivial components of the Einstein tensor G αβ ≡ R αβ − 1 2 R g αβ for the metric (2.1) are then: where we have defined radial and angular expansion rates:

2)
and also radial and angular acceleration rates:

4)
respectively, and also the following auxiliary quantities: Note that for λ = 0 andĖ = 0 the usual LTB expressions are recovered.
In the previous equations d/dτ rf = u α rf ∂ α = e −λ ∂/∂t is the derivative with respect to the proper time of the comoving observer. As we show in Appendix A, the components of the acceleration of u α rf are a r, rf = λ and a t, rf = 0. We see therefore that, for a geodesic referenceframe velocity field, the lapse function λ depends only upon time and can be rescaled such that λ = 0 and d/dτ rf = ∂/∂t. From Eqs. (2.2-2.3) it follows then that in a non-geodesic reference frame the comoving observers measure with different proper times. Later we will identify 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 = λ is sourced by pressure gradients which push the observers out of the freely-falling geodesics.

Einstein's equations and energy-momentum tensor
From the Einstein's equations G αβ = κ T αβ we can obtain four independent dynamical equations. We choose to form two of them with the following combinations: where κ = 8πG and T αβ is the energy-momentum tensor for an ideal fluid source comprised of n decoupled components: where the energy-momentum tensor of the i:th perfect fluid component is: 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 . The isotropic pressure is related to the energy density by p = w ρ, where the equation of state parameter 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 where v i, c is the coordinate comoving peculiar velocity of the i:th component relative to the reference frame. The proper peculiar velocities and the gamma factors are instead given by respectively. The covariant components of the total energy-momentum tensor are then: 14) As the third dynamical equation we will consider the combination coming from −G t t + G r r + G θ θ + G φ φ , which gives the generalization of the acceleration equation: where we used the fact that the acceleration scalar of the reference-frame velocity is The term proportional to a rf gives a "spurious" contribution to the acceleration and vanishes if we use a geodesic reference frame, in which the total acceleration is A R + 2A A , similarly to LTB models.
Finally, the last independent equation will be simply the G tr = κ T tr component, which reads: Eq. (2.17) shows that the evolution of the curvature is due to two distinct causes. The first is the effect of having an inhomogeneous multicomponent fluid and goes to zero in the FLRW limit where the peculiar velocities vanish or if only one fluid is present and its reference frame is adopted. 1 It is sourced by the energy flux in the radial direction: we remind indeed that the curvature function E may be interpreted as the total energy of a given shell (see Eq. (2.28)). The second contribution toĖ is due to the fact that, generally, we are using a non-geodesic reference frame (a rf = 0) in which the total energy of a shell r is not conserved.

Conservation equations
As we are considering an ideal fluid source comprised of n decoupled components, each i:th component will be conserved individually, i.e., ∇ α T αβ i = 0, where ∇ α denotes the covariant derivative. For each i:th component we obtain the two following equations 2 where we omit the index i and multiply for clarity's sake by e λ and v c Y 2 1+2E , respectively: where the expansion scalar Θ and the acceleration scalar a have been calculated in Appendix A. By taking the combinations 1+2E ∇ α T αr and e λ ∇ α T αt − 1 vc ∇ α T αr it is possible to obtain the relativistic energy-conservation and Euler equations, respectively, for non-comoving fluids: where, analogously to the total derivative with respect to the proper time τ , we defined the convective derivative with respect to σ along the four-acceleration a µ : Note that in the rest frame of the i:th fluid component the derivative with respect to τ reduces to (−g tt ) −1/2 ∂/∂t and the derivative with respect to σ to g respectively, where a r, rf = λ and a t, rf = 0 (see Appendix A). The last equation shows that the pressure gradients are the cause for the non-geodesity of the reference frame relative to theī:th fluid component Finally, because of the (twice contracted) Bianchi identity it holds ∇ µ G µν = 0 and so: Consequently, two conservation equations will not be independent and can be discarded. Equivalently, we may replace, if convenient, any other two dynamical equations with the latter two extra conservation equations.

Gravitating mass
The effective gravitating total mass F is given by from which follows that the curvature function E can be interpreted as the total energy per unit of mass of a shell, and that to it contribute the kinetic energy per unit of mass and the potential energy per unit of mass due to the total gravitating mass up to that shell. Note that thanks to spherical symmetry one is able to define a potential energy also in cases far away from 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 the separate masses of the i:th fluid components: where the angular "spatial" expansion rate S A = v c Y Y has been introduced in Appendix A. 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 mass F i is sourced by the gravitational influence of pressures and velocities of the other fluid components, differently from the case of the conservation equations (2.20-2.21) which apply to 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 discussed in Appendix B.

Full set of equations
The total number of independent equations is generally 4 + 2n − 2, where the first term is the number of nontrivial Einstein's equations, the second comes from the conservation equations and the last from the Bianchi identity. By including the definition of mass in Eq. (2.27) we have a total of 3 + 2n equations. We choose to use all the 2n conservation equations, two of which will replace the acceleration equation (2.16) and the evolution equation (2.30) for the total F . We now list for clarity the full set of 3 + 2n exact equations governing the dynamics of a spherically symmetric inhomogeneous model with a multi-component ideal fluid: where the unknown functions are Y , E, the total F and n copies of ρ i and v i, c . The referenceframe velocity u α rf can be taken to be geodesic so that a rf = λ = 0. Alternatively, one can identify the reference-frame velocity with the velocity of theī:th fluid component, u α rf = u ᾱ i , with the consequence that vī , c = 0 and λ becomes an unknown function in its stead. The equations of state w i have to be given as an external input. From Eq. (2.32) it seems clear that the dynamics is similar to the FLRW one, the complication being an array of interconnected equations 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,Ẏ > 0) it is: 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.

Initial and boundary conditions
Initial conditions. We will give initial conditions at the timet < t 0 , where t 0 is the present time. 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 inhomogeneous 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 the universe t A att, which gives an additional constraint between the total F and E. From Eq. (2.37), in the approximation that F (r, t) F (r,t) and E(r, t) E(r,t) for t <t, we can obtain: which can be used without solving the dynamical equations to demand att the homogeneous big 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 boundary condition 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 to point out that while the initial conditions can be given individually for the different fluids by considering Eq. (2.34) as n independent equations, this will not be true at later times for which only the total F matters. In other words, one cannot obtain the individual F i by integrating 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, in the example of Section 3, by Eq. (3.3).
If we identified the reference-frame velocity with the velocity of theī:th fluid component, u α rf = u ᾱ i , we may have to give initial conditions for λ (and not vī , 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 obtain with Eq. (2.25) the lapse function att 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 the expansion and acceleration scalars. Finally, λ(r, t) and F (r, t) are found by integration for the needed constant t. The equations that we have discarded thanks to the Bianchi identity may be used to cross-check the accuracy of the numerical integration.
Boundary conditions. In discussing the necessary boundary conditions we will particularize the analysis to a metric which is matched to an external FLRW model at some comoving coordinate radius r b . We will give boundary conditions at the initial time: the evolution equations will automatically maintain them at later times. Similarly to the embedding of an LTB model, the curvature and the integrated gravitating masses have to match the FLRW values, , respectively (in the chosen gauge the scale function is already matched as Y (r,t) = a(t)r). The latter condition for the case of v i, c (r,t) = 0 explicitly demands that ρ out 0 ρ i (r,t)r 2 dr, as it is easy to see from Eq. (2.34). According to the Darmois-Israel boundary conditions [59] the local density may be discontinuous through a timelike surface (constant r b ), but the pressure must 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 model in a geodesic reference frame where proper and coordinate times coincide. Even if not necessary it may be desirable to match smoothly the background, i.e., ρ i (r b ,t) = ρ out i (t) and 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 3 and E ∝ r 2 , as one can see from Eq. (2.27) where Y ∝ r and a finite H A = 0 is assumed (see also [37]). The latter also impliesẎ ∝ r and sȯ Y (0,t) = Y (0,t) = 0. Also vanishing at the center must be the velocities, v i, c (0,t) = 0, and the pressure gradient, p i (0,t) = 0, which implies λ (0,t) = 0. It seems again desirable to have densities and velocities smooth at the origin, ρ i (0,t) = v i, c (0,t) = 0.

Light propagation
The metric (2.1) implies that for radial null worldlines it is where u is an affine parameter. Let us consider a photon emitted at t(u) in a coordinate time ν(u) (emitting source and observer are taken to be in the rest frame). Expanding the right-hand side of Eq. (2.39) around t(u) and keeping terms up to first order in ν(u) we get The redshift is defined as 3 where the subscript 0 refers to values evaluated at the affine-parameter value u 0 at the observer's position. The derivative of z with respect to u is then Using (2.43) and remembering the expression for the rest-frame acceleration scalar (see Appendix A), Eq. (2.39) gives a pair of equations which determines the coordinates t and r on the light cone as functions of the redshift, and can be put in the following gauge invariant form: where r p is the proper radial distance. In the ΛLTB model whereĖ = a rf = 0, the previous equations reduce to the familiar form 3 Two non-comoving dust components in a flat ΛCDM universe As a first application of the formalism introduced in this paper, we will consider an inhomogeneous sphere embedded into a flat ΛCDM universe. The inhomogeneities will be given by two non-comoving dust components, one of which may be interpreted as baryonic matter and the other as dark matter. The reader should keep in mind, however, that these results are merely 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 relative velocity is set to zero. We will use the reference frame comoving with the component labelled by "M ", so that the relative Euler equation gives λ = 0. The other conservation equation reads: and can be directly integrated to give: wheret is the initial time. We will label the other dust component with "N " and use the simpler 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.

Initial and boundary conditions
First we set the flat ΛCDM model by fixing the present-day expansion rate to H 0 = 100h km s −1 Mpc −1 with h = 0.7 and the matter density parameter to Ω matter = 0.3. We then solve the background equations to find the evolution of the scale factor a(t) and so the initial timet corresponding toz = a(t 0 ) a(t) − 1, which we will fix illustratively toz = 5. By redshift z we will always mean the redshift relative to an observer in the background model, which will differ, for example, from the redshift relative to the observer at the center of the spherical inhomogeneity. The inhomogeneous sphere will have a comoving radius of r b = 100h −1 Mpc.
Next we give initial and boundary conditions att. For the scale function we set Y (r,t) = a(t)r. Then we have to give initial conditions for We stress again that the division of the total F into separate components is meaningful only at the initial time as only the total gravitating mass F enters the evolution equations. Note, however, that it is possible to define individual invariant mass components M i . For instance, for the M component it is: from which it is easy to verify using Eq. (3.1) that it is correctlyṀ M = 0 (whileḞ = 0). The remaining background matter density is split between the M and N components according to a fraction q such that it is ρ out M (t) = q ρ out matter (t) and ρ out N (t) = (1 − q) ρ out matter (t). We choose the profile of the M component to be: where the δ n function is given in Appendix C and δ M = δ M 1 + δ M 2 is the contrast in matter density between the center of the spherical inhomogeneity and the background value of ρ out M . This inhomogeneous profile satisfies the smoothness conditions ρ M (0,t) = ρ M (r b ,t) = 0 and continuously matches the background density ρ M (r b ,t) = ρ out M (t). Moreover, the matching requires thatF M (r) has the corresponding background value and so the following condition has to be satisfied: from which follows that it has to be δ M 1 δ M 2 < 0. δ M 1 > 0 will give an overdensity, while δ M 1 < 0 an underdensity. For a given model we will fix the quantities δ M 1, M 2 and α M 1, M 2 and we will get analytically the radius r t M by demanding Eq. (3.5). The gravitating mass is then simplyF M (r) = 4πa 3 (t) r 0 ρ M (r,t)r 2 dr, which has also an analytic expression. Next we needF N (r), which we will take to be: We can now find the initial curvature function E by demanding that the universe has the same aget for any r, i.e., from Eq. (2.38): t A (r,t) =t. The approximation from neglecting the time dependence of F and E for t <t is very good as the cosmological constant is negligible 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. We choose the initial velocity profile to be: v c (r,t) = δ 1 (r, r tv , δ v1,v2 , α v1,v2 ) (3.7) with δ v2 = −δ v1 so that there are no peculiar velocities at the center and at the border. The contrast δ v2 will give the maximum peculiar velocity in c units. Finally we can compute the initial density profile for the N component, which from Eq. (2.34) is: whereH A is given by Eq. (2.32) evaluated at the initial time. This terminates the necessary initial conditions.

Evolution
We will now show the exact general relativistic evolution of the system presented in the previous Section for the case of a central underdensity. 4 The left panel of Fig. 1 shows the In particular, this plot shows the relation between the coordinate radius r and the background comoving radius r bkg = Y /a which will be used in all the other plots. Right panel: Evolution of the gravitating mass F for z =5 (red), 3, 1 and 0 (blue) with the contribution due to Λ subtracted and normalized to unity at hole radius r b . The time evolution is purely an effect of the peculiar velocities between the two dust components. See Section 3.2 for more details.
evolution of the scale function Y : in the interior region the scale function is expanding faster because of the less matter and consequently negative curvature present, as also shown by the spatial Ricci scalar R in the left panel of Fig. 2 (see Appendix B for its definition). All the 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 1 (see right panel in Fig. 2) and thus the proper radial distance is Y . The evolution of the scale function Y is governed by the gravitating mass F and the curvature function E, as shown by Eq. (2.32). The evolution of the gravitating mass F is plotted in the right panel of Fig. 1. We have subtracted the time-dependent component due to Λ, which is F Λ = 4π 3 Y 3 ρ Λ , in order to show the genuine time dependence of the gravitating mass due to the dust components. As there is no pressure, the evolution displayed in the plot is due to the peculiar velocities between the two components. In particular the gravitating mass at half the radius is decreasing in accordance with the outgoing flow shown in the right panel of Fig. 4.
The evolution of the curvature function E, is plotted in the right panel of Fig. 2. This is again a genuine effect of the peculiar velocities between the two dust components, as a cosmological constant never sources the evolution of the curvature function E. We remind that 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 with the peculiar velocity field shown in the right panel of Fig. 4. As said before, the left panel of Fig. 2 shows the Ricci scalar R, which seems to follow a profile which is an average of the density contrasts of the two dust components, displayed in Fig. 3. We also remind that the curvature 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   . This definition again relies on the fact that r p Y . As expected these peculiar velocities are growing with time as the underdensity becomes emptier and emptier (see also the left panel of Fig. 1). Opposite is, however, the case of the peculiar velocities of the N component with respect to the M component (right panel of Fig. 4), which decrease with time: the gravitational attraction between the two components dumps 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 = 0 and r = r b . In Fig. 5 we re-plot the same curves, but with the velocity of the peak in density 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 both directions.
Finally, this model maintains the scale invariance valid for the simpler LTB model [9] if the initial peculiar velocities are properly scaled:ṽ c (r,t)/r b = v c (r,t)/r b where tilde marks the quantities relative to an inhomogeneous patch of different radius.

Collapse in the presence of dark energy with negligible speed of sound
As a second application of our general model we will consider an inhomogeneous sphere embedded into a flat wCDM universe. The inhomogeneities will be given by a pressureless matter component and by a dark energy source with constant equation of state w out < −1/3. This solution becomes the usual ΛLTB model if w = w out = −1 and the Lemaître model (possibly in non comoving coordinates) if the density of the dust component is set to zero.
We will use the reference frame comoving with the matter component, which we label with M , so that λ = 0 and Eqs. (3.1-3.2) hold. We will label the dark energy component with "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 negative pressure we cannot consider adiabatic perturbations for which w(r, t) = w out : this would indeed give a negative speed of sound (c 2 s = w out ) which would cause an unstable growth of perturbations [60]. We will consider instead the following r-dependent non-adiabatic equation of state: which implies the following speed of sound: so that for α ≤ 0 the speed of sound is positive (the adiabatic speed of sound is recovered with α = 1, see e.g. Ref. [61]). Particularly interesting is the case α = 0, for which the pressure gradients are absent and the speed of sound vanishes, with the consequence that matter and dark energy collapse together following the geodesic flow. The possibility of a dark energy fluid with negligible speed of sound has been investigated in the literature under various assumptions, and generally requires a non canonical scalar field like k-essence, as the standard quintessence models with canonical scalar fields always have c s = 1 (see, for example, [49,50,[62][63][64] and references therein). Apart from theoretical considerations, the fact that the speed of sound of dark energy may vanish opens up new observational consequences. Indeed, the absence of dark energy pressure gradients allows instabilities to develop on all scales, also on the small scales where dark matter perturbations become non-linear. We expect, therefore, dark energy to modify in a detectable manner the growth history of dark matter halos not only through its different background evolution but also by actively participating to the structure formation process, in the linear and non-linear regime [65][66][67][68]. The observational consequences of a clustering dark energy have been indeed extensively studied. See, for example, [69][70][71][72][73][74][75][76] for the 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 and dark 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 of sound as the peculiar velocities developed by the dark energy fluid will be very small (see the right panel of Fig. 9).

Initial and boundary conditions
We fix the background model by setting h = 0.7 and Ω M 0 = 0.3 as in Section 3. About the equation of state we will study the following three cases: w out = −0.8 and α = 0 in Section 4.2, w out = −0.8 and α = −10 −10 in Section 4.3 and w out = −1.2 and α = −10 −10 in Section 4.4. The inhomogeneous sphere will have a comoving radius of r b = 10h −1 Mpc and we will give initial conditions at the timet corresponding to the background redshiftz = 1000. We set 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 initial matter inhomogeneity as: is the background gravitating mass, δ M is the density contrast, and the parameters r 0 and ∆r characterize respectively size and steepness of the density profile. The parameter r offset is used to translate the profile in order to have a smooth origin. We will adopt the following values: r 0 = r offset = 0.25 r b , ∆r = 0.3 r 0 and δ M (t) = 1.5 · 10 −3 , that is, we have a central overdensity. The initial matter density is then: Next, we can find the initial curvature function E by demanding that the universe has the same aget for any r, i.e., from Eq. (2.38): t A (r,t) =t. AsF w (r) F M (r) att, this condition will be independent of the dark energy initial conditions. In particular at the initial 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 1 [77]: These initial conditions are free from decaying modes [78]. Finally, we take the dark energy component to have no initial peculiar velocities, which means that it will have, with respect to the background, the peculiar velocities Y ∆H of the matter component. Therefore, in order to have initial conditions without decaying modes in the dark energy component, we have to set the initial profile as: whereF out w = 4π 3 a 3 (t)r 3 ρ out w (t) and initial contrast is: which is valid during matter domination for sub-Hubble perturbations and c 2 s 1 [79]. The initial dark energy density is then: With this we have all the necessary initial conditions.

Evolution for zero speed of sound
The left panel of Fig. 6 shows the evolution of the scale function Y for w out = −0.8 and α = 0: in the interior region the scale function is expanding slower because of the more matter and consequently positive curvature present. All the other figures will be plotted with respect to the background comoving radius defined as r bkg = Y (r, t)/a(t), which as explained in Section 3.2 is a good definition of distance. The right panel shows the peculiar velocities of the 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 central overdensity becomes denser and denser. In Fig. 7 the evolution of the density contrasts of the matter and dark energy components is shown. The evolution, as remarked above, is realistic. In particular the dark energy component closely follows the evolution of the matter component, as one should expect from a fluid with vanishing speed of sound. Because of the absence of pressure gradients peculiar velocities between the two components never develop (if initially zero) and the two fluids both follow the same geodesic flow. 5 As it is clear from the plot, however, the growth rates are different. This has to be expected because the negative pressure generally slows down the evolution of perturbations. To be quantitative, we have plotted in right panel of Fig. 8 the evolution of the ratio of the contrasts, δ w /δ M , which we normalized to the value valid during matter domination, 1+wout 1−3wout . As one can see, the agreement is perfect when Ω M ∼ 1 (see left panel of Fig. 8), showing that the chosen initial conditions are indeed free from decaying modes. During the nonlinear dark-energy-dominated phase the ratio departs as expected from the matter-domination prediction, and in particular becomes different in overdense and underdense regions, as illustrated in the right panel of Fig. 8 where the blue line is relative to r bkg = 1h −1 Mpc (overdensity) and the green line to r bkg = 6h −1 Mpc (underdensity). This could be due to the fact that even though the two components have similar profiles, in the overdense region the matter component reaches the nonlinear evolution earlier (the dark energy is actually almost linear at the present time).
It is also interesting to compare the evolution to the case of no dark energy. We find that the matter perturbation grows to a maximum contrast of about ∼ 90. Alternatively, we find that an initial perturbation of δ M (t) = 1.1 · 10 −3 grows to the same present-day contrast as the initial perturbation of δ M (t) = 1.5 · 10 −3 studied in this Section.
Finally, let us comment on the fact that because the speed of sound is zero, there are no characteristic length scales associated to the dark energy clustering and the spherical collapse remains independent of the size of the object [65]. The characteristic length scale associated to the dark energy clustering is indeed the sound horizon scale, L s = a c s dt/a, which vanishes for c s = 0 so that clustering takes place on all scales. In other words, this model maintains the scale invariance valid for the simpler LTB model [9].

Evolution for c s ∼ 10 −5
We will now consider the model of the previous Section (w out = −0.8) with α = −10 −10 , which corresponds to a speed of sound of order c s ∼ 10 −5 , as one can see from the left panel of Fig. 9. This latter value implies a sound horizon of order L s ∼ 80h −1 kpc, much less than the inhomogeneity scale considered in this example, so that the overall dynamical evolution shown in Figs. 6-8 is basically unchanged. The new features of this case are pressure gradients and peculiar velocities. In order to understand their behavior it is useful to work out the following approximations valid for |δ w |, |α| 1: From the previous relations it is easy to understand how speed of sound and pressure gradients are related to the density contrast. In particular it is δ p ∼ c s ∼ −δ w , as one can also see by comparing the left panel of Fig. 9 with the right panel of Fig. 7. As a consequence, the pressure gradients will push the dark energy component out of the free-falling geodesic away from the overdense regions, as shown in the right panel of Fig. 9. One can indeed obtain at linear order (see, e.g., [80]) the approximate relationv p ∼ − c 2 s 1+wout δ w .

Evolution for w < −1
Finally, we will now consider the case of w out = −1.2 and α = −10 −10 . We chose the same order of magnitude for the (positive) speed of sound of the previous Section for comparison sake, see the left panel of Fig. 11. Clearly, we can decrease the magnitude of α in order to prevent possible pathologies [62]. As one can see from the left panel of Fig. 10, the evolution of the matter component is similar to the previous case, with a moderately higher contrast at late times probably due to the different expansion history: for w out = −1.2 the dark-energy dominance happens indeed later with respect to the case for w out = −0.8. 6 Quite different is instead the evolution of the dark energy component, which still follows the matter flow but develops an underdensity rather than an overdensity (see the right panel of Fig. 10). This can be understood from the fact that the relation δ w /δ M ∼ 1+wout 1−3wout changes sign for w out < −1. It is interesting to note that a similar behavior was found in the analysis of Ref. [66] of perturbations in scalar-tensor cosmologies, where the non-minimal coupling of the field is indeed responsible for the crossing of the so-called phantom divide (w out = −1). As one can see from the right panel of Fig. 11, the evolution of the peculiar velocities is instead approximately the same one of the previous Section, as in the relationv p ∼ − c 2 s 1+wout δ w both δ w and 1 + w out change sign.

Comparison with previous work on exact solutions
In this Section we will compare our findings to previous work dealing with exact solutions. As explained in Section 1 the model here presented is a generalization to the case of n decoupled and non-comoving perfect fluids of the Lemaître model [1] (see [2][3][4][5] for recent contributions) which describes the dynamics of a spherically-symmetric perfect fluid. Therefore, the model of Eqs. (2.32-2.36) straightforwardly reduces to the Lemaître model, possibly with the addition of a cosmological constant.
In the case of one source only, our equations can be used to describe the fluid in noncomoving possibly-geodesic coordinates, which may be desirable if the fluid has pressure gradients and its four-velocity undergoes acceleration. Within our formalism this is accomplished by simply taking a geodesic reference-frame four-velocity u α rf , i.e., by setting λ = 0 (see Eq. (A.18)). This can be useful to disentangle the effect of pressure gradients on the evolution of the model as shown, for instance, by Eq. (2.17) for the time evolution of the curvature function. In the latter equation, indeed, only the second term on the right-hand side contributes if comoving coordinates are used and only the first term on the right-hand side contributes if geodesic coordinates are used. Spherical models have been discussed previously in the literature in non-comoving coordinates, albeit with a different purpose. The main idea was indeed that an exact solution with a simple appearance in non-comoving coordinates may become extremely complex when transformed to the appropriate comoving system, or that the integration of the partial differential equations necessary to obtain the comoving coordinates may even represent an intractable mathematical problem. We refer the interested reader to [81][82][83][84][85][86][87] and to [88] and references therein.
Another line of research focused on interpreting two or more fluid components as a single effective fluid, see for example [89] and references therein. The work of [90,91] showed indeed that a mixture of fluids whose four-velocities lay on a two-plane can always be reinterpreted as a single anisotropic fluid. The model of this paper satisfies this last requirement and it would be interesting to investigate in future work the amount of anisotropic pressure generated by the evolution of a cosmological multi-component fluid modeling, e.g., baryons, dark matter and dark energy. As a simple example we can apply here the work of [90] to the model of Section 3 featuring two non-comoving dust components in a flat ΛCDM universe. Following [90], the energy-momentum tensor of the effective fluid is: and a short calculation shows that the rest energy density is: and that the components of the diagonalized stress tensor are: where p A and p R are the angular and radial pressures, respectively. Clearly, for γ → 1 (v c → 0) the effective fluid becomes simply the sum of the two fluids and p A = p R = 0. For a nonzero radial velocity v c = 0, however, the radial pressure is positive p R > 0 and the effective fluid is anisotropic. See [90] for more details. Finally, we would like to mention that the authors of [91] claim that the opposite process is also feasible, i.e., that given a particular anisotropic fluid, it is always possible to find a multi-fluid model with an equivalent energymomentum tensor. Note, however, that this last claim was argued not to be valid for a general anisotropic fluid [92] such as the anisotropic fluid models of compact objects like neutron stars [93,94]. Finally, [95] studied anisotropic pressure in the context of the LTB metric in order to investigate the evolution of inhomogeneities in an interacting photon-baryon mixture. A physically-plausible hydrodynamical description of cosmological matter in the radiative era between nucleosynthesis and decoupling was given, emphasizing its thermodynamical consistency. It would be interesting to develop such considerations also in the presence of a dark matter component using the multi-fluid model presented in this paper.

Conclusions
We have extended the Lemaître model to the case of n decoupled and non-comoving perfect fluids with general equations of state. We have expressed the full set of 3 + 2n exact equations governing the dynamics of the model in a concise and transparent way thanks to the use of physically meaningful quantities like, for example, expansion and acceleration rates, expansion and acceleration scalars, total derivatives along four-velocity and acceleration. This general 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 universe where dark energy is not the cosmological constant and, in particular, can be inhomogeneous. This possibility may also induce an inhomogeneous variation of fundamental constants if a coupling between the dark energy and the matter-radiation Lagrangian is allowed. Moreover, dark matter and baryons can be described as two separate fluids, the latter possibly featuring pressure. Pressure in general can have a non-negligible effect on the cosmological models and this should be taken into consideration while interpreting cosmological datasets. At early times the contribution of radiation can be included, which may be relevant for the understanding of the evolution of the inhomogeneities.
In this paper we focused on the formal understanding of the general solution, which we applied to two different setups. For the case of two non-comoving dust components in a flat ΛCDM universe we found rich and interesting physics as, for instance, the gravitational interaction between the two fluids which makes one component follow the collapse of the other, or the evolution of the peculiar velocities which are increasing with respect to the background but decreasing between the two fluids.
Finally, for the case of clustering dark energy, we have been able to follow in an exact way the collapse of an overdensity to which also the dark energy contributes, and we confirmed that the evolution agrees with the theoretical growth rates during matter domination. We also considered the case of a small speed of sound (c s ∼ 10 −5 ) and of a phantom equation of state, which also showed the expected dynamics. These applications show some of the interesting features of this n-fluid component solution, which we will further investigate in forthcoming work, both from a theoretical and observational point of view.
h αβ is often referred to as the "spatial metric" because it is the induced metric on the space slices, h αβ = h µ α h ν β g µν . Moreover, h αβ is the first fundamental form of the hypersurface and the expansion tensor is (minus) the extrinsic curvature, K αβ = −Θ αβ , which is the second fundamental form of the hypersurface.
The expansion scalar can be obtained without calculating explicitly the components of the expansion tensor by remembering that it is: g is the determinant of the metric and we have defined the following quantities: 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 the temporal, radial and angular expansion scalars: respectively, where we defined the quantities: 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 the expansion scalar for a geodesic four-velocity field is Θ geod. = Θ R /γ 2 + Θ A . In order to calculate the shear we need the non-vanishing components of Θ αβ : 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: The shear tensor σ αβ is then defined in terms of Θ αβ as 14) and its non-vanishing components are: where the invariant σ of the shear tensor is given by: We may now particularize the above expressions for the case of the reference-frame velocity u α rf : Note that, as Θ T vanishes, λ appears in the expansion scalar and shear only through the derivatives with respect to the proper time d/dτ rf = u α rf ∂ α = e −λ ∂/∂t. This shows that a non-geodesic reference frame only changes the proper time with which the expansion rate is measured, but not its functional form (different is the case for the acceleration rates of Eqs. (2.4-2.5)). Moreover, these last results show that the definitions in Eqs. (2.2-2.3) are relative 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 the Lemaître metric [1] in a possibly non-comoving reference frame. We will discuss some properties of this solution, which may help us better understand the multi-fluid non-comoving case. We will not discuss the dynamical evolution of the Lemaître metric, which can be easily deduced from the set of Eqs. (2.32-2.36).
Let us start by considering the Eqs. (2.30-2.31) forḞ and F . By taking the combinations giving the derivative of F with respect to τ and σ it is possible to obtain the following simpler expressions: where the Euclidean and actual volumes of the space slices are, respectively: which has the usual thermodynamical interpretation for a fluid in equilibrium (at constant entropy): for an increase dV of volume at constant t the total mass-energy is increased by the mass-energy density of the volume dV , and for an increase dV at constant r the total massenergy 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 volume element and so the gravitating (or Euclidean) mass F does not coincide with the invariant mass M which is related to the local density by [7]: We can then evaluate the time derivative of M at constant r: Note, however, that in this last equation not only a different volume is used with respect to Eq. (B.6), but that Eq. (B.9) is not valid for p = 0, which implies a time-dependent curvature (that can be interpreted as total energy per unit of mass). In this latter case the usual thermodynamical expression, where only the pressure at the boundary matters, does not seem to hold and the value of the pressure at any r matters as shown by the last integral in Eq. (B.8).
Finally, it is interesting to rewrite Eq. (2.32) as: where the Euclidean average for the generic quantity X is defined as: and R is the spatial Ricci scalar of the rest frame, which is the trace of the Ricci tensor of the metric g ij on the hypersurface of constant t: The factor of 1/6 in Eq. (B.10) can be understood by the fact that for the FLRW model it is R/6 = k/a 2 . From Eq. (B.10) follows that the expansion rate is sourced by the Euclidean average of the local density and curvature, and not by the actual averages. This is a potential source for "strong" backreaction effects [96,97] of the inhomogeneities on the evolution of the 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, 1[. In this paper we used W 1 and W 3 , which have the following explicit form: From the W n function we then built the following general contrast δ n : where ∆ = ∆ 1 + ∆ 2 is the contrast between the center (r = 0) of the spherical inhomogeneity and the border (r = r b ), and ∆ 2 is the contrast at r t . We have used δ n to model inhomogeneities in ρ M , F N and v c in Section 3.