Growth rate of spherical voids with non-comoving dark matter and baryons

We present numerical solutions to Einstein’s equations describing large spherical cosmic voids constituted by two components: dark matter and baryons, with a non-vanishing initial relative velocity, in an asymptotically homogeneous background compatible with the ΛCDM concordance model. We compute numerically the evolution of such configurations in the dark matter frame, with a hypothetical homogeneous distribution of baryons, but respecting the values dictated by the concordance model for the average baryon-to-dark matter density ratio. We reproduce the well-known formation of overdensities at the edge of the void and recover the Lemaître–Tolman–Bondi solutions in the comoving limit of our simulations. We compute the average growth factor of matter fluctuations and find that it departs significantly from the linear perturbative prescription even in the comoving case, where the non-linearity of inhomogeneities has an impact.


Introduction
Since the first evidence of the current accelerated expansion, many models have been proposed to introduce a dark energy component [1][2][3], and subsequently observables have been sought for probing such alternatives to a cosmological constant both at the background level [4,5], as well as in the statistics of structures [6][7][8][9].In particular, the growth function f = d log D/d log a is an observable widely studied in the galaxy distribution [10,11] in the search for deviations from the standard theory [12].In theory, this is customarily calculated from the evolution of perturbations within the model and parametrised with the growth index γ, where f = Ω γ m [13].The observational counterpart is constrained through the Redshift Space Distortions (RSD) manifest in the dipole and quadrupole of the galaxy power spectrum [14].The different prescriptions for the growth function in several dark energy models has motivated the advancement of surveys of great precision in the determination of this observable [15][16][17][18].For example, several models propose a scale-dependence of this function, aside from the redshift dependence allowed by the above parametrisation.The experimental evidence allows for such dependence and deviations from the standard ΛCDM prescription are still open.
In the precise characterization of this observable it is important to account for possible effects from velocity bias (see e.g.[19]).Specifically, if the baryon field is not exactly comoving with the dark matter structures, then the relevance of a bias in the determination of the growth function demands special attention.Since the magnitude of such bias is largely unknown, its effect could well lie beyond a perturbative one, and nonlinear models of structure formation are required to account for its influence on the growth function.
While it is standard to describe the evolution of cosmological structures using linear perturbations on a FLRW background, it is insufficient when quantities such as the density contrast δ reach non-linear values.Alternative approaches for the description of non-linear phenomena are needed such as "spherical collapse" [20,21] and inhomogeneous cosmological models.Inhomogeneous cosmological models have the advantage over linear perturbations and spherical collapse that they are exact solutions to the field equations that reduce to FLRW models at the appropriate limit [22].These models can be used in a variety of scenarios (see [23] for a comprehensive discussion of examples) including the modeling of cosmological structures.Important examples of the usefulness of these types of solutions include the quasispherical Szekeres Class I models, where it is possible to model arrays of multiple structures arranged in a spheroidal manner [24][25][26].
Relevant to this study, a particular advantage of inhomogeneus cosmological solutions is the ability to describe several matter components, such as non-comoving fluids with relative velocities between them.This description can range from simple models containing a single structure [27] to more complex arrangements of structures [28].Relative velocities are important to model since peculiar velocities are measured in supernovae data which can lie within a range of 328 -620 km/s [19].Characterizing the effect of (non-linear) velocity bias on the growth function is imperative to determine the possible values of this observable within the standard model of cosmology.Such results may help interpret observations within the ΛCDM model before considering more exotic theories.
In this paper we look at a toy model realization of the effect of non-linear bias [29] in cosmic voids, which are a suitable scenario to test gravity [30][31][32][33].With the technique deveolped in a previous paper [27] we study the late-time evolution of spherically symmetric voids with non-comoving components of dark matter and baryons.We take the dynamical system of the Einstein Field Equations (EFEs) and solve for the evolution of two fluid components out of the proper frame.We show how these components are evolved numerically from high redshifts to the present cosmic time.Our formalism is suited for inhomogeneities with arbitrary amplitude (particularly beyond the perturbative regime).We look at the growth factor f (z) for departures of the averaged void and emerging overdensity, and compare it to the perturbative prescription.We find that even a small (perturbative) initial curvature inhomogeneity yields sizable differences with the linear prescription, mostly for non-comoving fluids, and to a lesser extent for the comoving case.
The paper is organized as follows.Sec. 2 reviews the 1+3 formalism of general relativity to describe general fluids.In Sec. 3 we derive the system of equations which describe the evolution of a dark matter plus baryons configuration with a non-trivial profile of relative velocity.In Sec. 4, we present the evolution of a large void with two non-comoving components, in an otherwise homogeneous ΛCDM universe.Finally in Sec. 5 we discuss the prospects of detection of the effect that a relative velocity brings to the discrepancy of the growth factor from its perturbative prescription.

The 1+3 Formalism
In General Relativity there are two main formalisms for the splitting of spacetime in space and time.The 3+1 formalism deals with the foliation of spacetime through a 3D spacelike hypersurface Σ t 0 at t = t 0 and a timelike vector field (orthogonal to these surfaces) which defines the time evolution to find hypersurfaces Σ t at t ̸ = t 0 .In this way, the spacetime is foliated by the different Σ t hypersurfaces.In contrast the 1+3 formalism threads spacetime by defining a timlike vector field, usually the 4-velocity u a of a set of special observers, and through it defining orthogonal spacelike hypersurfaces.We will give a brief discussion on this formalism, an in depth discussion can be found in Ref. [34].
The core of the 1+3 formalism is the 4-velocity u a (usually comoving) and the projector operator h ab defined as It is worth mentioning that this tensor works both as a projector tensor from the 4-dimensional spacetime to the 3D hypersurfaces as well as a metric on those hypersurfaces.Having these quantities one can split the covariant derivative of u a in terms of its irreducible kinematical quantities where ω ab is the vorticity tensor, σ ab the shear tensor, Θ the expansion scalar, and ua = u b ∇ b u a the acceleration 1-form.These quantities are used to split the Einstein equations into their 1+3 form together with the electric and magnetic parts of the Weyl tensor E ac = C abcd u b u d and H ab = η acd C cd be u e /2 respectiveley, where η acd = η acde u e = − |g|ϵ acde u e is the volume form of the spacelike hypersurfaces.Now, for a general stress-energy tensor, with elements and we have the covariant 1+3 evolution equations, derived from the Einstein field equations: ) ) ) ) ) -3 -the above system is complemented by the constraint equations derived from the Ricci identities ) ) ) (2.17)

Spherically Symmetric Equations for Non-Comoving fluids
We will now move to the problem in question: The evolution of a spherically symmetric void constituted by two non-comoving, self-gravitating pressureless fluids.We first present in detail the equations that describe the evolution of our configuration.We first show how a relative velocity between fluid components turns the matter fields oberver-dependent and secondly, we use a spherically symmetric metric to obtain the explicit equations for the two non-comoving dust components.

Relative velocity between fluids
Between two or more non-comoving fluids, there is a transformation between frames depending on the different relative velocities.For relative velocities v a I between a component with 4velocity u a and the components with 4-velocity u a I , we use the inverse transformation given by where , and the subindex I denotes the I-th fluid component with a 4-velocity u a I .This results in the following components of a transformed stress-energy tensor for multiple fluids with a barotropic equation of state The * index in ρ * I denotes that it is measured in its respective proper frame, and w I = p * I /ρ * I .

Spherically Symmetric Inhomogeneous Equations
We study a spherical void density configuration of two distinct dust fluid components, Cold Dark Matter (CDM) and baryons, in a ΛCDM background, so we consider the following inhomogeneus metric For this metric, the 4-velocity is given by u a = δ a t /N (t, r).The covariant objects (Hubble scalar, shear tensor and the electric part of the Weyl tensor) are given, in terms of scalar functions and the basis tensor for spacelike, symmetric and trace free Petrov D spacetimes where n a = √ g rr δ r a , Σ = Σ(t, r), and W = W (t, r).Note that we are choosing non-rotating fluids so ω ab = 0 = H ab .Aditionally, we define the relative velocity, the energy flux and the anisotropic stress tensor, also in terms of scalar functions, as The variables can be complemented by the Misner-Sharp mass [36], which yields an alternative expression for the electric Weyl scalar Finally, by defining χ = ∂ r Y and K = 1 − χ 2 /B 2 , the 3-dimensional Ricci scalar (3) R and the extrinsic curvature K take the form As stated before, we are considering two dust components representing CDM and baryons.We make the choice of CDM as the comoving frame with 4-velocity u a in a ΛCDM background.The thermodynamical variables for the fluids are, according to the velocity transformations (3.2), given by It is important to note that since we are considering dust components, then N = 1 which implies that ua = 0. Finally, we take the initial Hubble scalar H * and the initial characteristic void size l * to obtain a set of dimensionless variables (3.10) -5 - We note that in order to obtain the standard redshift z we need to define a scale factor a.
For this we take the asymptotic radial value of the metric function B(t, r) since where a(t) is the FLRW scale factor.This also means that When we mention the redshift, it is obtained through these asymptotic values.The resulting evolution equations are then given by ) ) where Ψ = ∂ t Ψ/H * for any variable Ψ.The system is complemented by the Hamiltonian constraint, the Weyl constraint and the Misner-Sharp mass (where equation (3.13) is plugged into (3.6)) The previous system of equations is equivalent to that presented in (Ref.[27]).

A two-component spherical void
We study two particular cases of a spherically symmetric central void: one with an initial extrinsic curvature profile of higher amplitude than the other.For both cases we match an asymptotic background universe representing a flat ΛCDM cosmology with fiducial values for the density and Hubble parameter taken from the Planck 2018 results [37] 1 .
The initial profiles for the velocity, curvature and density (for the CDM) are given by the following Gaussian profiles  [38].Note that the comoving case of V c = 0 represents the Lemaitre-Tolman-Bondi (LTB) solution.We take the baryonic density profile as non-trivial only in the comoving case.That is: Since the LTB case imposes a null relative velocity between baryons and dark matter, the equations naturally yield similar radial profiles.On the other hand, that restriction is absent for a non-zero relative velocity, in which case we chose a homogeneous baryonic density profile as an initial condition (as measured in their proper frame of reference).The initial profiles for CDM, velocity and curvature are shown in Fig. 1.
To analize the growth of structure, we take the definition of the growth function f as given in (Ref.[39]).In the linear regime, this is defined in terms of the growing mode D + of the linear density contrast In a ΛCDM background, the analytical growth function in the linear regime obeys the powerlaw (see e.g.[40]) This prescribes the evolution, at different values of z, for the density contrast in the case of comoving baryons and dark matter as shown in Figs.5-6.We generalize the definition of the growth function, by first defining the mean density contrast in non-linear structures.We take the average given by The notation ⟨⟩ D denotes the quantities in brackets averaged over the domain D. Also ⟨⟩ .

D
denotes the time derivative of the averaged quantity.The motivation for adopting these definitions is presented in Appendix A. We divide the profiles in two averaged structures as shown in Fig. 4, representing a void region and the over-dense spherical shell.From the initial conditions described at the  3) and consider it as the void region while we take the average for the positive values representing the over-dense spherical ridge and consider it as the over-density region.Both of these regions are represented as the filled in regions in both graphs.
beginning of this section, the resulting growth function f after the numerical evolution for the k c = 0.05 case is shown in Fig. 5 within a range of z ∈ [0, 20] for multiple velocities.On the other hand, the case for k c = 0.009 is shown in Fig. 5 for the same range of z.They are shown compared to the analytic function obtained from linear perturbation theory.As we can see in Fig. 2, the density contrast of both components behave in a similar way .We can clearly see how the void expands as time passes and even achieves non-linear amplitudes for the density contrast (δ CDM ∼ −0.8 at z = 0.1).Additionally, we observe the formation of a spherical over-dense region surrounding the void region (as represented by the positive values of the density contrast) which also achieves non-linear values (δ ∼ 0.3 for the CDM and δ ∼ 0.4 for the baryons at z = 0.1).This is a generic result of the evolution of cosmic voids [21,41,42].We observe in all cases that the behaviour further deviates from the analytic function the higher the maximum velocity is.A thing to note is that for the case of no relative velocity between the matter components the CDM and baryonic density contrast and, as a result, their growth function behave exactly in the same manner.As the velocity increases, the discrepancy between the components' behaviour increases.For the total matter density contrast, the effect of the relative velocity is less pronounced as the separate components one.As one can see in both figures, for both curvature cases the numerical growth function gets closer to the analytic one for lower z values for the case of the over-dense ridge excepting the total growth function where the lower z values are where the function differs mostly from the analytic function.For the void region we see the opposite effect, where the lower ranges of redshift are the ones where the numerical and analytic growth functions differ the most with the cases of higher relative velocity exceeding a difference of 40%.For the case of the lower curvature, the deviation from the analytic function is still present but to a lesser extent than for the case of the higher curvature.

Discussion
In the previous sections we have described the 1+3 splitting of space-time and the Einstein equations.We described how this splitting plus an inhomogeneous model helps us describe a two-component cosmological system that takes into account a relative velocity between components and uses fully non-linear equations.The effect of a relative velocity between two matter components on a spherical void has been previously studied in [27] showing the effects this velocity has on the evolution of the void.Similar to that previous work, we developed a numerical system and extended the work to analyze the growth factor f .We can see in figures 5 and 6 that even the inhomogeneus LTB case with zero relative velocity has a clear distinction from the perturbative case and the relative velocity between baryons and CDM further intensifies this difference.-13 -For the following discussion we focus on the low curvature, with k c = −0.009,case.In figure 7 we plotted the percentile difference between the analytic growth function of (4.4) and our numerical growth functions for the different velocities and the lower curvature compared to the expected margin of meassurement from the DESI collaboration [16] plotted in a redshift range of z ∈ [0, 2] for consistency with DESI's report.
We divert our attention to the over-dense ridge first: as we can see from the left handside plots of the figure for both the baryons and CDM separately the numerical plots lie mostly outside the ranges given by DESI.The one exception is the LTB plot corresponding to no relative velocity between components.All lines corresponding to all the relative velocity amplitudes tested lie within the margin of error of the closest redshift meassurement error available for both fluid components.The effect of having a higher value for |f /f (Ω) − 1| increases with an increase in the relative velocity distribution amplitude; this same effect of increase is seen more pronounced for the case of the baryonic component.Finally when the total density contrast is considered the effects of relative velocity in the growth function are much less pronounced.
By considering now the void region (right hand-side panels) we see a similar effect but more pronounced.Even for the case of CDM and no relative velocity the curve lies outside most of the values of DESI expected errors; while for the baryonic case the two curves corresponding to the higher velocity lie completely outside of all DESI expected errors.One difference between the behaviour of the void and over-dense regions is that the behaviour of the total density contrast for the void lies in between the CDM and baryonic components instead of being the most close to the DESI expected values.With all this in mind we conclude that such effects produced by relative velocity between field sources and, to a lesser extent, inhomogeneities is meassurable in future observational collaborations.
There are claims that the non-linearity of LTB models with large density gradients in the supercluster scale is merely an effect of using a comoving gauge that disappears when passing to the conformal Newtonian gauge and considering the values of realistic Newtonian peculiar velocities defined as v pec = Ṙ − HR, with R the areal radius of the LTB model (see [43,44]).According to these authors, this supports the claim that the Universe is quasi-Newtonian at these scales [45,46].These claims are sustained on a gauge transformation from the comoving to the conformal Newtonian gauge applied to the linear approximation of the LTB metric in comoving and non-comoving frames and neglecting higher order terms on v pec ≪ 1.However, the definition of these authors of peculiar velocities is completely artificial and ad hoc, while we have defined peculiar velocities properly in terms of the energy flux of an energy-momentum tensor in the context of two dust sources in comoving and non-comoving frames, even if we have also assumed that v pec ≪ 1.In fact, the reasoning of [43,44] might apply only to the idealized unrealistic case of a single dust source and an ad hoc artificial definition of peculiar velocities.As a contrast, we have considered in this paper a more realistic situation of two sources and a correct relativistic definition of peculiar velocities.Our results clearly shows that the dynamics of the non-comoving dust are relativistic and non-linear, distinct from the dynamics of dust in a conformal Newtonian gauge of linear perturbations, even if we also assumed peculiar velocities to be Newtonian and an almost LTB metric ((3.3) with N = 1).This issue is outside the scope of the present paper, so we will discussed in a separate article.

B Removing σ 8 from the DESI report values
In figure 7 we present plots for the quantity |f /f (Ω) − 1| and place the DESI values for the expected errors in order to present a comparison.However the DESI collaboration presents results for the product f σ 8 , so we must remove the σ 8 from the DESI values for a proper comparison with our results.To that end, we start with equation (53) from [40]  Using the values given in tables 2.3 and 2.5 from [16] for the errors for ∆f σ 8 (z) obtaining for each z given (in the notation of that same reference) ∆f σ 8 (z) = 1 100 Finally, we obtain the desired value for ∆f (z) that we present in figure 7 ∆f (z) = ∆f σ 8 (z)σ 8 (z) − f (z)σ 8 (z)∆σ 8 (z) σ . (B.4) ab is the derivative projected along the flow of the 4-velocity and curl [ S ab ] = η cd(a ∇ c S d b) .

r 2 e 2 .( 4 . 1 )
−( r−0.01 0.025 ) 2 and K i = k c r 2 e −( r 0.03 ) From the equations (4.1) we take the two specific cases of differing curvature by using the values for the amplitude constant k c = −0.05 and k c = −0.009.The different velocities are modulated by the amplitude constant V c , and we present 6 different velocities for each of the curvature profiles: V c = 0, 0.2, 0.6, 1.0, 1.4, 1.8 corresponding to maximum relative velocities, expressed in km/s, of V max = 0, 28.472, 85.416, 142.361, 199.305, 256.249 respectively -All of these representing values compatible with the peculiar velocities according to statistics from recent surveys

Figure 1 :
Figure 1: Initial profiles of the void configuration: the upper left is the initial Gaussian dark matter density profile, the upper right is the initial relative velocity profile for V c = 1.8, and the lower panel is the graph of the initial intrinsic curvature normalized by the amplitude |k c |.

Figure 2 :
Figure 2: The left panel shows the evolution of Dark Matter's density contrast for different values of z.The right panel shows the evolution for the ratio between the two component's densities, ρ b /ρ CDM in the dark matter frame.Both panels correspond to the higher amplitude in curvature, k c = −0.05.Both figures are the result of a relative velocity profile with V c = 1.8.

Figure 3 :
Figure 3: Like both figures in Fig. 2 the left panel shows the evolution of Dark Matter's density contrast while the right panel corresponds to ρ b /ρ CDM for different values of z.Both figures are the result of a relative veocity profile with V c = 1.8 but with a curvature with an amplitude constant of k c = 0.009.Note the considerable amplitude difference in the profile is due to the curvature magnitude.

Figure 4 :
Figure4: As seen in both graphs, we take the region of negative δ values and average it over the volume as defined by the metric (3.3) and consider it as the void region while we take the average for the positive values representing the over-dense spherical ridge and consider it as the over-density region.Both of these regions are represented as the filled in regions in both graphs.

Figure 5 :
Figure 5: For the higher amplitude in curvature, k c = −0.05. the panels are aranged as follows: top left figure represents the growth function f of CDM for multiple velocities of the over-dense ridge, top right belongs to the void region for CDM, middle left and middle right are the growth functions for the baryonic component for the over-dense ridge and the void region respectively and, finally, the bottom right and bottom left are the growth function for the total matter density contrast for the over-dense and void regions respectively.All individual figures include the percentage difference between the different growth functions for a given velocity profile and the analytic linear function.

Figure 6 :
Figure6: The arrangement of these graphs is exactly the same as in Fig.5but for the evolution of the inhomogeneity with a curvature amplitude k c = 0.009.

Figure 7 :
Figure 7: The figures are arranged in the exact same way as the previous one.These graphs display the percentile difference that every graph in Fig. 6 possesses for a range of z ∈ [0, 2] with the added expected ranges to be measured by the DESI colaboration.