Non-comoving baryons and cold dark matter in cosmic voids

We examine the fully relativistic evolution of cosmic voids constituted by baryons and cold dark matter (CDM), represented by two non-comoving dust sources in a $\Lambda$CDM background. For this purpose, we consider numerical solutions of Einstein's field equations in a fluid-flow representation adapted to spherical symmetry and multiple components. We present a simple example that explores the frame-dependence of the local expansion and the Hubble flow for this mixture of two dusts, revealing that the relative velocity between the sources yields a significantly different evolution in comparison with that of the two sources in a common 4-velocity (which reduces to a Lemaitre-Tolman-Bondi model). In particular, significant modifications arise for the density contrast depth and void size, as well as in the amplitude of the surrounding over-densities. We show that an adequate model of a frame-dependent evolution that incorporates initial conditions from peculiar velocities and large-scale density contrast observations may contribute to understand the discrepancy between the local value of $H_0$ and that inferred from the CMB.

In the present article we examine the fully general relativistic evolution of a spherically symmetric cosmic void, assuming as matter source a mixture of non-interactive baryons and CDM species, each evolving along a different 4-velocity. Specifically, we consider the evolution of a generic cosmic void suitably embedded in a Λ CDM background. Since CDM is the dominant clustering source, we assume a frame in which its 4-velocity is comoving, whereas the baryons evolve along a non-comoving 4-velocity that defines a non-trivial field of spacelike relative velocities with respect to the CDM frame. Consequently the 4-velocities of the two dust sources are related by a boost, and the energy-momentum tensor of the mixture (as described in the comoving frame) no longer has the form of a perfect fluid, but that of a complicated "fluid" energy-momentum tensor that contains effective pressures and energy flux terms associated with the relative velocity field.
In order to deal with the general relativistic dynamics for this energy-momentum tensor, we do not resort to the traditional metric based methods of [72][73][74][75], but consider the system of first order partial differential equations (PDEs) provided by the "fluid-flow" (or "1+3") representation of Einstein's equations [76][77][78] in terms of evolution equations and constraints for the covariant quantities associated with the CDM comoving 4-velocity. These covariant variables are (i) kinematic: expansion scalar, 4-acceleration, shear tensor and relative velocity, all computed from the CDM frame; (ii) source terms: the total energy density, pressure and energy flux that arise from the relative velocity and (iii) the electric Weyl tensor. These equations (evolutions and constraints) must be supplemented by spacelike constraints.
The plan of the paper is as follows. In Section 2 , we introduce a model for the evolution of a generic mixture of fluids in a spherically symmetric spacetime. This model is specialized to the case of two non-interacting dust-like fluids, namely CDM and baryons, in Section 3. We examine the numerical solutions of the resulting system of partial differential equations (PDEs) in a void formation scenario. Through representative numerical examples, we look at the influence of the relative baryon-CDM velocity on the evolution and present-day final structure. Our results are summarized and discussed in Section 4. Finally, we have included three appendices that complement the main text. Appendix A provides the general tensorial evolution equations of the 1 + 3 description; while in Appendix B and Appendix C, we present the dimensionless baryon-CDM system of PDEs and show its (single-fluid) LTB limit.

A mixture of multiple fluids in spherical symmetry
A spherically symmetric spacetime is characterized by the line element, where the metric coefficients N, B and Y are functions of the radial and time coordinates. Notice that this metric contains as a particular case the Robertson-Walker line element, a solution recovered in our examples at scales larger than ∼ 100 Mpc at z = 0.
Regarding the evolution of a mixture of non-comoving fluids, the fluid elements of each species will evidently present its own 4-velocity. In absence of specific criteria, a useful choice of "reference frame" is the one of the dominant species 2 . Hence the family of fundamental observers will evolve with comoving 4-velocity, which defines the convective fluid-flow (or time derivative) and space-like gradients (orthogonal to u µ ): The kinematic parameters associated with a spherically symmetric fluid as measured by the fundamental observers are the expansion scalar, 4-acceleration and shear tensor (the vorticity vanishes identically): The Einstein field equations can be recast [76][77][78] as a first order system of "1 + 3" evolution and constraint equations involving these kinematic parameters, together with the energy density ρ, isotropic p and anisotropic pressure π µν and energy flux q µ of the source given by the energy-momentum tensor, as well as the electric Weyl tensor E µν = C µανβ u α u β (the magnetic Weyl tensor identically vanishes). This system (displayed in Appendix A) is ideal for a numerical treatment in which all variables have a clear physical and geometric meaning. Since it is based on a 4-velocity flow it is fully covariant, and thus it is readily applicable to spacetimes that are not spherically symmetric. For the spherically symmetric metric (1) we have and the space-like symmetric trace-free tensors σ µ ν and E µ ν can be written as Here W and Σ are scalar functions, with Ψ 2 the conformal invariant of Petrov type D spacetimes, and e µ ν = h µ ν − 3n µ n ν = Diag [0, −2, 1, 1] is the tensor basis that serves as eigenframe for spacelike symmetric trace free tensors in Petrov type D spacetimes, with h µν = g µν + u µ u ν and n µ = √ g rr δ r µ the projection tensor and a spacelike normal vector tangent to the orbits of SO(3) (note thatė µ ν = 0). The 4-velocity of the other non-comoving components are related to u µ via the relative velocity measured by the fundamental observers v µ (i) , defined such that v µ (i) u µ = 0. Then, the 4-velocity of the i-th fluid reads, where "i" labels the components and v 2 (i) = g µν v µ (i) v ν (i) . The total energy-momentum tensor is made up of all the contributions from the different species, and in general it will no longer be the energy-momentum tensor of a perfect fluid, but where ρ, p, π µν and q µ are the energy density, isotropic and anisotropic pressures 3 , and the energy flow measured by the fundamental observers along u µ . These components are determined by projecting the energy-momentum tensor parallel and orthogonal to u µ [78,79]: Although the total energy-momentum tensor is always conserved, the energy-momentum tensors of the individual components are not necessarily conserved. If there are non-gravitational interactions between them, they satisfy ∇ ν T µν (i) = J µ (i) , where J (i) is the rate of energy and momentum densities transfer between species i as measured in the u µ -frame. In absence of non-gravitational interaction these energymomentum tensors are separately conserved: J µ (i) = 0 for all i.

A mixture of non-interacting perfect fluids
We now focus on the case of a mixture of non-interacting fluids, each one a perfect fluid with a suitable equation of state in its intrinsic frame (denoted with * ): In this way, the total energy-momentum tensor follows from adding up the corresponding tensors of the dynamically significant species as seen from the u µ frame (eqs. (9) and (10)). Explicitly, if we choose the fundamental observers those along u µ (0) , then the contribution to the total energy-momentum tensor of the "0" component reads On the other hand, the energy-momentum tensor of the i-th component comoving with velocity u µ (i) (see eq. (8)) takes the form with the dynamical quantities given by [78,79]: where for spherical symmetry spacetimes the anisotropy tensor can be written as π The dynamics of this fluid mixture can be determined from the first order "1 + 3" fluid flow representation of Einstein's field equations given in Appendix A, by direct substitution of ρ, p, π µν , q µ by (14a)-(14d), with andu µ , σ µν , E µν obtained from (5) and (7). For spherical symmetry this system can be further simplified and complemented by evolution equations for the metric functions Y, B and χ ≡ Y that can be obtained from (5) and (7). Since each one of the proper tensorsu µ , σ µν , E µν is fully determined by a single scalar function A, Σ , W , the tensorial system in Appendix A becomes a dynamical system involving only scalar functions: This system must be complemented by the following constraints: Here K is the spatial curvature: and M is the Misner-Sharp function that is characteristic of spherically symmetric spacetimes. Notice that this function furnishes an expression for W through the constraint (17a) that allows us to eliminate the evolution equation (16f) forẆ (and W in (16e)), though, since M is fully expressible through (19) in terms of the variables of (16), we do not need to use this function explicitly to integrate this system (we just eliminate M with (19)) 4 . Besides these constraints, we also need to supplement the system with the conservation equation for each fluid: Finally, the radial component velocity V (i) of the i-th fluid in (15) can be determined algebraically from (14a) and (14c).

Void evolution from a mixture of two decoupled dusts
The main characteristic of cosmic voids is the underdensity profile that depends on the (roughly) radial distance on Mpc scales. However, the usual single fluid (dust) approach generally focuses on the void dimensions (size and the depth of the density contrast) and the value of the local expansion [80], while the relative velocity between the dynamically significant species is usually ignored. In this section we show that the multiple components scenario brings important modifications to the evolution of cosmic voids.

A numerical example of the two-component mixture
To stress the above it is illustrative to look at the case of a mixture of two dust fluids identified as CDM and baryonic matter, including a cosmological constant characterized by the present-day parameters from Planck 2015 [81], in order to accommodate a Λ CDM asymptotic background model. We consider the fundamental observers comoving with dark matter u µ DM = δ µ t ; consequently, the baryonic matter will have 4-velocity, and the energy-momentum tensor, eq. (9), will be the sum of the CDM and baryonic components: such that one can identify the following quantities (as defined in (14) and (15)), Note that the pressure, heat flow and anisotropic stress terms are zero when V = 0 (CDM and baryons with common 4velocity: LTB limit). Consequently, the evolution will be governed by the system of equations that results from substituting 0 → DM and i → B in eq. (16) and considering the energy-momentum tensor variables (22). The resulting dimensionless system of equations is presented explicitly in Appendix B.
We examine the numerical solutions of this two-dust system in a grid simulating a cosmic void of present-day radius ∼ 60 Mpc. Starting from linear initial conditions at z = 23 we follow its evolution until z = 0 (see below for the justification of this choice of initial redshift). The initial CDM density, spatial curvature, and the relative velocity profiles are taken as Gaussian functions of linear amplitude with respect to the background parameters. In all our simulations the baryonic density is initially homogeneous and equal to its value in the background (as seen in its intrinsic frame), In the expressions above µ c ∼ 0.01, k c ∼ 0.05, σ K = σ µ = 0.03, σ v = 0.025, and r = l * ξ . From this, the spatial curvature (K ) is derived from Eqs. (18) and (23b). The characteristic length is l * ∼ 60 Mpc, while the characteristic speed constant V c (and the maximum of the velocity V peak ) will be specified further below. Figure 1 shows the typical initial profiles used for the numerical analysis.

Evolution of density profiles
To look at the effect of the relative velocity on voids and wall formation we develop a code capable of handling test cases with given initial densities for each species, a given curvature profile, and a series of profiles for the relative velocity (V peak ). In Fig. 2 we display the baryon and CDM density contrasts at z = 0 for different initial velocity profiles. As a reference, we have included the case in which both baryons and CDM are comoving (LTB solution). We find that even nonrelativistic relative velocity values exert non-trivial effects on present-day configurations as density contrasts become non-linear. On the other hand, the void size depends on the sign of V , so that smaller voids result from initially negative values for the relative velocity. We also illustrate the evolution of the density contrast profiles for the specific case of V peak ∼ 7 × 10 −3 (corresponding to the red curves in Fig. 2), snapshots for different values of z are displayed in Fig. 3.
Notice that as the evolution proceeds the density contrast at the surrounding wall increases, reaching probably a shell-crossing singularity. We interpret this as the onset of an intricate virialization process, a stage of structure formation that marks the limit of validity of the dust model, and that lies beyond the scope of this work (discussed elsewhere in the literature [83,84]). Since our purpose is to look at the simultaneous evolution of the matter-energy components (CDM and baryons) within the void before the onset of virialization, we have chosen z = 23 as the initial time slice, simply because it is easier to set the initial conditions at this time than at, say, the linear regime of the last scattering time z 1100, well before gravitational clustering becomes dynamically significant. However, these initial conditions are idealized but not fine-tuned or unrealistic, they simply correspond to a spherically symmetric realization of the generic spectrum of random CDM and baryon perturbations, characteristic of the linear regime at the last scattering surface z 1100, which evolve to produce a void of the desired size.

Local expansion of the components
Let us now focus on the effects of a relative velocity in the measure of kinematic quantities. In the case of a single fluid, the comoving observers define a natural threading of the spacetime by the future-directed unit timelike vector field u µ . In our case, as we stressed before, the choice of the fundamental observers is not unique, and observers comoving with each fluid will measure the kinematic quantities with different magnitudes. In fact, due to a change of frame, the local expansion of CDM (H DM ≡ H in this frame) will depart from the expansion of the baryonic matter H B , which is given by where h B µν = u µ B u ν B + g µν is the projection tensor and u µ B is the 4-velocity of the baryonic matter given in eq. (21). We find by computing (24) a relation between both estimations of the expansion,  Fig. 1) set initially at z = 23. The red lines stand for V peak ∼ 7 × 10 −3 , the blue lines for V peak ∼ 5 × 10 −3 and the green ones for V peak ∼ −2.6 × 10 −3 . As a reference, we have provided the case without a relative velocity (LTB model), denoted by a black line. Note that the solution displayed by the solid blue curve includes a baryonic matter shell of width ∼ 10 Mpc and density contrast of the order the unity and peculiar velocity of ∼ 500 km/sec with respect to the CDM comoving frame. This configuration is roughly comparable with the dynamics of our local group, which has similar size and density contrast and a dipole velocity of ∼ 600 km/s associated with its local motion with respect to the CMB frame [82].
where in order to derive eq. (25) we have used the fact that V 1 (but its derivatives need not be small) and χ ≡ Y ,r . Hence, the difference in the local expansion due to a change of frame can be expressed as follows, Figure 4 shows the difference between the two expansions H B and H DM at z = 0 for the solutions whose density contrast are depicted with red and blue lines in Fig. 2, corresponding to V peak equal to 7 × 10 −3 (red lines) and 5 × 10 −3 (blue lines). Note that this difference can be of the order of km/(s Mpc), around the maximum of the baryonic matter density (even larger differences are expected at times close to virialization). This estimation is roughly that of the discrepancy between the values of H 0 reported by CMB and SNe observations [81,85,86], thus suggesting that considering a relative velocity between baryons and CDM may provide interesting clues to understand this issue (though this task lies beyond the scope of the present work).

The baryon-CDM relative velocity
Since for the baryon-CDM mixture the radial component of the relative velocity, as defined in Eq. (21), can be determined from the algebraic relation with Q B and ρ B given by Eq. (22d)-(22e), its evolution equation iṡ where we used (20b), (20c) and A ≡ 0. In order to relate this equation to the well-studied perturbative case, we drop the term of order of V 3 to obtain, which shows the connection between the time evolution of the relative velocity and the radial gradients of the velocity and the metric function B that generalizes the background scale factor. In a quasi-homogeneous perturbative regime B ∼ a(t, r) r and thus B ,r > 0 should hold, while V ,r > 0 should also hold because relative velocities increase from the center onwards as r increases. Therefore the derivative of (V 2 )/(2B 2 ) should be positive and thus the right-hand side of the equation above negative. As a consequence,V < 0 holds and relative velocities dilute asymptotically during cosmic expansion. However, this is not applicable to a nonperturbative regime where large gradients of the involved variables may occur and/or change signs, so that the relative velocity can be amplified by a local inhomogeneity.
We can obtain further information on the evolution of V by looking at a definition of peculiar velocity often used in a perturbative approach: the difference between the local Hubble flow relative to the Hubble flow of the background, which can be estimated as v pec = (H local − H FLRW )Y [87]. Then, once again neglecting the highest power of V in Eq. (26), we get which shows that such spatial velocity field is intrinsically related with the local homogeneities. At large scales (in a perturbative regime) this field evolves by approximately diluting as the inverse of the background scale factor, since in a regime approaching FLRW-like conditions our variables can be written as B ∼ a(r,t) and Y ∼ ra(r,t) (see e.g., [88] for a formal equivalence of LTB models with Cosmological Perturbation Theory in the linear regime). In an inhomogeneity, however, where the spatial gradients are not restricted to small values, the relative velocity and peculiar velocity must be found by a non-trivial evolution equation. In particular, for the numerical solutions showed in this section we found that the relative velocity decreases, but without following a trivial scaling law in the spatial region identified with the walls. Note that in such regions the gradients can be large and the local expansion is slower than the background expansion, in fact, in part it is locally collapsing.

Discussion and final remarks
We have considered the fully relativistic evolution of spherically symmetric cosmic voids made up of a mixture of two non-comoving dust components, identified as CDM and baryonic matter. Specifically, we looked at the effects of the baryon-CDM relative velocity on the void properties. We found that for baryons converging to the center of the void, as seen from the CDM frame, the final density profile shows an effective reduction on the size of the void (see Fig. 2). On the other hand, if the baryon component is receding from the centre, the void presents a deeper (baryonic) underdensity, and the walls manifest a larger density contrast as illustrated in Figs. 2 and 3. The existence of a relative velocity between baryons and CDM leads to a difference in the expansion of each component. We find that small initial differences in velocities between two components (of order 7-5 ×10 −3 ) yield important differences in local expansion of the order of km/(s Mpc), similar to the gap between local and CMB measurements of  5 Peculiar velocity of the void components with respect to the asymptotic background as defined in [87]. Note that the larger differences occur at the radii corresponding to the wall structure.
the expansion parameter H 0 [64,81,85,86]. Indeed, this last result may be part of the effects missing in the usual single frame analysis of peculiar velocities and local expansion (e.g. curvature effects [67,89], among others). Related to this, we find significant differences in the peculiar velocities of each component, defined as deviations from the asymptotic background (common) expansion, reached at large radii. Such differences could be interpreted as the velocity bias field, here evolved to non-linear stages. Fig. 5 shows that such bias manifests most prominently at the peak of the density contrast (walls of the void).
The spherical void model we work with is qualitatively analogous to earlier models [72][73][74][75]. As in these models, we obtain qualitatively analogous results that depict the expected streaming of baryons determined by the rapid void expansion, a characteristic also present in Newtonian models. However, the dynamical equations employed in the past are based on a numerical scheme constructed from the Misner-Sharp mass that is completely tied to spherical symmetry. As a contrast, the system of evolution equations and constraints here considered is based on covariant fluid-flow scalars that can be computed for any spacetime, regardless of the symmetry considerations.
Our approach to void dynamics could also represent an important improvement over the "silent models" of [67] that try to address this issue through "emergent" spatial curvature. Silent models (characterized by a non-rotating dust source with purely electric Weyl tensor) are theoretically handicapped by the conjecture stating that Einstein's equations may not be integrable in general under the "silence" assumption [90,91] (the models in [67] also neglect matching conditions among the different silent cells). By assuming dust sources with different (non-comoving) 4-velocities, the resulting models are based on similar physical assumptions but are no longer silent because of the non-trivial energy and momentum flux among the dust sources.
In conclusion, the fully relativistic evolution of baryons and CDM along different 4-velocity frames can provide important clues in understanding the observational tension in the estimation of the value of H 0 from local observations and through interpretation of the Planck data. A concrete example is furnished by the study of the Hubble flow in the non-spherical models examined in [64], which tries to understand this tension, but did not consider different 4velocities for the baryon and CDM components. This work could be improved by allowing for a non-comoving baryon 4-velocity that would provide more degrees of freedom as we have done in this paper. Likewise the multiple fluid approach can provide important corrections to the usual study of the process of formation and growing of large-scale structure in the universe. Finally, we emphasize the fact that the system of evolution equations and constraints used in our numerical modeling has been constructed with covariant fluid flow variables, and thus it is readily applicable (under certain restrictions) to examine self-gravitating systems that are much less idealized that those under the assumption of spherical symmetry. system of PDEs governing the dynamics of a 2-dust-fluid mixture:Ŷ At each time the velocity and intrinsic density of the second (non-comoving) dust is determined by: The system is complemented by the following constraints where, On the other hand, the equation (16f) (redundant) results, We employed the Method of Lines to integrate this system of PDEs. Proceeding in this way the PDEs were discretized along the radial variable, setting 1000 grid points within the interval r/l * ∈ [0, 0.2]. The resulting set of ordinary differential equations was integrated using an adaptive step-size Runge-Kutta of 4(5)-th order.

Appendix C: LTB limit
The LTB model is a general inhomogeneous spherically symmetric solution of the Einstein's equations for a single irrotational dust fluid as source T µν = ρu µ u ν . The time-synchronous metric can be cast as follows [92], ds 2 = −dt 2 + R 2 ,r (r,t) 1 + 2E(r) dr 2 + R 2 (r,t) dθ 2 + sin 2 (θ )dφ 2 . (C.17) For a comoving 4-velocity u µ = δ µ t the field equations reduce to: Following a similar approach to that used in the main text, we rewrite the Einstein's equations in terms of covariant objects associated with the 4-velocity and the energy-momentum and projection tensors, which leads to, together with the constraint (among others not listed here) where the expansion scalar, the eigenvalues of the shear and magnetic Weyl tensors as well as the spatial curvature take the simple form: These previous equations are recovered by specializing the system (16)-(19) (or its particular case (B.15)) to the single comoving fluid case. This can be checked by making the following substitutions: (C. 22)