Non-comoving Cold Dark Matter in a $\Lambda$CDM background

We examine the evolution of peculiar velocities of arrays of inhomogeneous cosmic structures in a $\Lambda$CDM background identified as a comoving CMB frame. These arrays are constructed by smoothly matching to this cosmological background regions of Szekeres-II models with their imperfect fluid energy--momentum tensor reinterpreted as non-comoving dust, keeping only first order terms in $v/c$. We find a peculiar velocity field with magnitudes compatible with values reported in the literature. We obtain a a contrast of the Hubble scalar of the same order of magnitude of the current tension on $H_0$. While the models cannot describe the virialization process, we show through a representative example that structures of galactic cluster mass reach the onset of this process at redshifts around $z\sim 3$.


Introduction
It is a well known fact that cosmic structures at different scales are not comoving with a frame of reference associated with the CMB frame identified as the frame of a ΛCDM background. This fact follows from measured and inferred peculiar velocities between Cold Dark Matter (CDM) structures and the CMB [9]. These peculiar velocities are clearly non-relativistic (up to 3000 km/s), which justifies studying their dynamical evolution by means of Newtonian gravity as a good approximation [27]. However, relativistic effects might not Sebastián Nájera Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México (ICN-UNAM), A. P. 70 -543, 04510 México D. F., México. E-mail: sebastian.najera@correo.nucleares.unam.com be negligible when considering the superposition of nonrelativistic peculiar velocities on scales comparable with the Hubble horizon. Studies of large-scale peculiar velocities [21,13] have shown a connection between them and the anisotropy and inhomogeneity on large but still subhorizon scales.
Since the scale dependence of these velocities is still an open topic, it is worth studying this dependence by means of a general relativistic approach that complies with their non-relativistic magnitudes. This approach involves considering congruences of observers with different 4-velocities, leading to distinct energy momentum tensors so that peculiar velocities result from momentum and energy fluxes between the congruences [8]. Different congruences of observers define different Hubble flows characterized by the kinematic quantities associated with their 4-velocities. In particular, Tsagas has examined the effect of considering peculiar velocities under this approach on the interpretation of cosmological observations (see [25,26] and references therein). Another example of non-comoving matter can be found in [10], where evolution of cosmic voids constituted by baryons and CDM is analyzed.
The simplest example of an observed effect that can be attributed to non-comoving observers is the peculiar velocity of our local Hubble flow with respect to the CMB frame, which manifests itself in the large observed CMB temperature dipole. In particular, we can define the CMB frame as the one in which this dipole vanishes [3]. Besides relativistic models of peculiar nonrelativistic velocities in the context of non-comoving observers [27], there are non-perturbative models based on exact and numerical solutions of Einstein's equations that consider non-comoving observers in complete generality (tilted models), [7,10,14].

arXiv:2011.11192v2 [gr-qc] 24 Nov 2020
Since spherical symmetry is too idealized and limited, it is useful to consider the class of exact solutions derived by Szekeres [19,15] that (in general) do not admit isometries and thus enhance the available degrees of freedom for applications in cosmology. These solutions are classified in two classes (I and II, Szekeres-I and Szekeres-II hereafter), each one of which subdividing in three subclasses: quasi-spherical, quasi-flat and quasi-hyperbolic, depending on their symmetrical limits.
The quasi-spherical Szekeres-I models are regarded as the most suitable for cosmological applications and thus have been widely used to address the limitations of the spherically symmetric Lemaître-Tolman-Bondi (LTB) models (see comprehensive discussion in [19,15,2]) In these models mass-energy and all physical and geometric objects appear as the superposition of a dipole on top of the LTB monopole of spherical symmetry [19], thus allowing for the description of two or more structures: typically a central monopole (over-density or void) evolving together with an elongated wall type structure ("pancakes") corresponding to the dipole. This extra degrees of freedom provided a significant enhancement to the "Big Void" models that were proposed ten years ago to account for observations without resorting to dark energy or a cosmological constant [4]. More recently, these models were used to describe multipole configuration involving an arbitrary number of structures in a ΛCDM background [22,24], providing an appealing coarse grained rendering of CDM structures at 100 Mpc scales that can be made consistent with the observed cosmography [23].
While most applications to cosmology involve the quasi-spherical Szekeres-I models (see review in [2]), Szekeres-II models have also been used for studying inflationary scenarios [1] and probing the structure growth factor [12,18]. However, we believe that the potential of Szekeres-II models for cosmological applications has remained largely unexplored.
In a recent paper [17] we proved that quasi-plane Szekeres-II models admit a smooth matching, along an arbitrary number of matching interfaces, with spatially flat FLRW models. This fact leads to appealing toy models of evolving arrays of multiple inhomogeneous and anisotropic "pancake-like" cosmic structures (regions of Szekeres-II models) embedded in a homogenous and isotropic background (we briefly review these models in section 6).
In the present paper we extend the results of [17] by considering Szekeres-II regions whose energy-momentum tensor has an imperfect fluid form with nonzero energy flux, thus generalizing the exact "heat conducting dust" solutions found by Goode [11], whose source is no longer interpreted as a dissipative fluid (difficult to justify for CDM sources), but as dust in a non-comoving frame (neglecting the subdominant baryon contribution) with non-relativistic peculiar velocities. Since the resulting Szekeres-II regions can be matched to a ΛCDM background with peculiar velocities vanishing at the matching interfaces, the comoving frame associated with this background can be regarded as the CMB frame in which the dynamical contribution of the photon gas is neglected. A similar model was derived in [6] but only considering Szekeres-II models with a comoving dust source.
Having set up the models, we find that their free parameters allow for the description of evolving CDM structures falling into the comoving CMB frame with peculiar velocities consistent with observed values. Specific numerical examples are provided in full. By computing the Hubble scalar we show that its contrast with respect to its ΛCDM background value H 0 is entirely determined by the shear tensor, producing fluctuations of H 0 of the same order of magnitude ∼ 10 % difference that has emerged in the "H 0 tension", though peculiar velocities have a negligible ∼ 0.01 % effect on present day values of H 0 between comoving and non-comoving frames.
While we show that the parameters allow for an evolution of Szekeres-II structures that is free from shell crossing singularities, we argue that these shell crossings mark the limit of validity of the dust description of CDM which necessarily breaks down at the onset of the virialization process. Hence, we present a numerical example of a CDM structure with mass M ∼ 10 15 M (roughly the mass of a galactic cluster) arriving to the onset of virialization at cosmic times z ∼ 3 that are compatible with structure formation scenarios derived from numerical n-body simulations [5].
The section by section description of the paper is as follows. We discuss in section 2 the definition and handling of peculiar velocities between two general imperfect fluid congruences with distinct 4-velocities. In section 3 we introduce Szekeres-II models whose energy momentum tensor is that of an imperfect fluid with nonzero energy flux (generalizing Goode's solutions [11]), showing in section 4 how to re-interpret this energy flux in terms of a field of non-relativistic peculiar velocities as a linear limit of the formalism presented in section 2. In section 5 we present the physical, kinematical and geometric variables of the system as covariant perturbations in a FLRW background. We provide in Section 6 a brief review of [17] and specify the parameters that will be used in section 7 to examine in full detail the peculiar velocity field and the relevant physical, kinematic and geometric variables of the models, including the evaluation of a difference in H 0 between the CDM and CMB frames that is of the order of magnitude to the H 0 tension. The range of validity of the models (i.e. of the dust description for CDM) is examined in section 8 providing a numerical example of galactic cluster structures starting to virialize at the expected redshifts. Finally, we provide two appendices: Appendix A discusses the consistency of assuming peculiar velocities defined up to first order in v/c and in Appendix B how the FLRW limit for Szekeres-II models can also arise without the smooth matchings with FLRW regions that we have considered.

Imperfect fluids in terms of peculiar velocities
Given a 4-velocity field the most general form of the energy-momentum tensor is given by where h ab = u a u b + g ab is the projection operator, ρ, p, π ab , q a are the mass-energy density, the isotropic pressure, the spacelike tracefee anisotropic pressure tensor and the spacelike energy flux vector. The energy-momentum tensor (1) is also referred to as describing "imperfect fluids" because of the terms π ab , q a that are usually identified with dissipative stresses, shear viscosity and heat conduction, in thermal and hydrodynamical systems (see examples in [15]). However, this interpretation is not suitable for gravity dominated long range interacting cosmic sources: CDM is best described at large scales as dust, a description also applicable to baryons whose internal energy and thermal dissipative effects are also negligible in these scales [20].
A more useful interpretation for the "imperfect fluid" terms π ab , q a in a cosmological context follows as 4momentum fluxes identified with peculiar velocities associated with a non-comoving 4-velocity with respect to a comoving frame. In particular, we can assume a 4-velocity comoving with a ΛCDM background for the CMB frame, with CDM and baryons evolving along a 4-velocity field that is not comoving with respect to this frame.
To examine the connection between energy-momentum tensors associated with different frames, we consider two general congruences of spacetime observers with different 4-velocities u a andû a . Choosing u a as a comoving 4-velocity, the non-comoving 4-velocityû a is given by the generalization of the Special Relativity boost where 1 v a is the spacelike peculiar velocity measured by the observer u a and v a u a = 0. Following [16] and assuming the energy-momentum for the non-comoving frame to have the general form (1), the relations between dynamical quantities of the non-comoving and comoving energy-momentum tensors are given by where (as in [16]) we have written terms linear in v a outside the curly brackets. Notice that the isotropic and anisotropic pressure are mostly associated with therms that are non-linear (at least quadratic) in v a .
3 Szekeres-II model with energy flux

The general case: geometric and kinematic properties
Szekeres-II models are exact solutions of Einstein's equations characterized by the line element where X = X(t, x i ) with x i = w, x, y. The canonical orthonormal tetrad e a (α) associated to this solution is: where g ab e a (α) e b (β) = η (α)(β) . Szekeres-II models in a comoving frame are compatible with the most general energy-momentum tensor (1) with u a = e a (0) = δ a 0 and we will assume henceforth a nonzero cosmological constant Λ > 0.
In general Szekeres-II models do not admit isometries, hence ρ, p, and the components of π ab and q a depend in general on the four coordinates x a = t, x i . The rest frames (hypersurfaces orthogonal to u a ) are conformally flat, while the sign of k determines the sign of the constant curvature of the 2-surfaces marked by t and w constant, leading to three general sub-classes: quasi-spherical (k = 1), quasi-plane (k = 0) and quasihyperbolic (k = −1) models. The latter sub-classes contain axially symmetric and higher symmetry particular cases with spherical (k = 1), plane (k = 0) and hyperbolical (k = −1) symmetries.
The spherical, plane and hyperbolic Kantowski-Sachs spacetimes (X = X(t)) are the natural homogeneous limits of the models in general. All models admit smooth matchings with their Kantowski-Sachs homogeneous subcases, though as proven in [17], the general quasi-plane models admit also a smooth matching with spatially flat FLRW models along 3-dimensional matching hypersurfaces marked by w = const.
In a comoving frame the only nonzero kinematic parameters are the expansion scalar Θ = u a ;a and shear tensor σ ab = u (a;b) − (Θ/3)h ab given by The expansion tensor Θ a b and its three eigenvalues Θ a b = λ (i) δ a b : provide a covariant description of the kinematic anisotropy by identifying two equivalent principal directions λ (2) = λ (3) along e a (x) and e a (y) , which are clearly different from λ (1) along e a (w) . The inhomogeneity and anisotropy of the models can also be appreciated from the local rate of change of redshift z along null geodesics [8] an expression that must be integrated along light rays, parametrized by the affine parameter ϑ, with tangent null vector k a . Note that redshift distribution of light received by local observations along a comoving worldline (x i fixed) is isotropic only if h b a Θ ,b = 0 and σ ab = 0 along the wordline.

The generalized Goode solution
In what follows we will consider the quasi-plane (k = 0) particular case of (1) given by p = Π ab = 0 but ρ, Λ, q a = 0 that generalizes the quasi-plane "heat conducting dust" exact solution found by Goode [11]. Expressed in cylindrical coordinates, x = r cos φ, y = r sin φ, the metric takes the following form: where where we notice that X has the form of a monopole with two independent superposed dipoles A(w, r, φ) and B(w, r, φ), becoming only a single dipole A(w, r, φ) in the perfect fluid subcase (Q = 0). The main advantage of cylindrical coordinates is dealing with a bounded coordinate 0 ≤ φ ≤ 2π and the fact that they mark in a simple way the class of privileged observers along the curve r = 0 parametrized by w (or x = y = 0) that is a geodesic and an integral curve of a Killing vector of the rest frames. Axial symmetry follows by the restriction (14)- (15) so that X becomes independent of φ.
The functions S(t) and F (t, w) are found by solving the differential equations with (16) formally identical to the Friedman equation for the pressure in FLRW models, thus suggesting an identification of S(t) with an FLRW factor, an identification that becomes rigorous by characterizing this function as identical to the scale factor of the FLRW spacetime that can be smoothly matched to the Szekeres-II model along surfaces of constant w (see [17]).

Non-comoving CDM
The energy-momentum tensor associated with the Szekeres-II models under consideration is whose interpretation in Goode's original article [11] was as "heat conducting" dust with Λ = 0. We propose a wholly different interpretation as non-comoving CDM modeled as dust in which q a is proportional to its peculiar velocity with respect to a comoving CMB frame associated with a ΛCDM background. This novel interpretation of (18) follows from (3)-(6) withû a associated with a non-comoving dust energy-momentum tensor:T ab = (ρ + Λ)û aûb , hencep =π ab =q a = 0. With respect to the comoving CMB frame and using (3)-(6) this energy-momentum tensor takes the form (1) with with ρ and q a linear with respect to v a but isotropic and anisotropic pressures quadratic in v a . However, observed and inferred peculiar velocities of large scale structures with respect to the CMB frame are clearly non-relativistic, thus it is well justified to keep only terms that first order in v a in (19)- (20). This leads to the following energy-momentum tensor in the comoving frame which is (up to first order in v a ) exactly like (18) with q a ≈ρ v a and ρ ≈ρ (since γ ≈ 1). This energymomentum tensor in the comoving frame can be regarded as an approximation to a more realistic one given by where the comoving CMB radiation and non-comoving baryon densities: ρ cmb and ρ b can be neglected in comparison with the non-comoving CDM density and Λ (i.e. ρ cmb + ρ ≈ ρ CDM ). We believe that this interpretation of their energy-momentum tensor furnishes a solid physical connection to the models under consideration.
To complement the interpretation of energy flux as non-comoving CDM, it is useful to compute the Hubble scalar for the non-comoving 4-velocityΘ =ĥ a bû b ;a withû a defined by (2). In the linear regime of peculiar velocities v a v a /c 2 1 we obtain the following relation whereĥ ab = g ab +û aûb ≈ h ab and v a = q a /ρ (at linear order).

Dynamical variables
In order to work with dimensionless variables we normalize the dynamical variables with respect to the present day critical density κ/(3H 2 0 ) where κ = 8πG/c 4 and H 0 the present day Hubble length. The energy density can be expressed as the sum of a purely time dependent ΛCDM density (a solution of (16)) plus a term depending on all coordinates that can be conceived as an exact fluctuation over this homogeneous background (see [17] for details) where S is given by the analytic solution of (16) with τ = H 0 t the dimensionless time and Ω m 0 is the Omega factor associated with CDM (neglecting the baryon contribution).
The peculiar velocity field follows from q a = ρ v a (first order on v a /c). From the field equations we have q a = q r δ r a + q φ δ φ a with q w = 0, hence: Ω qr ≡ κqr To examine the inhomogeneity and anisotropy of the Szekeres-II region it is useful to consider the contrasts of the normalized density Ω ρ and Hubble scalar with respect to their ΛCDM values where (from (8)) It is important to remark that (see [17]) the quantities δΩ ρ and δH are covariant fluctuations respectively related to the electric Weyl and shear tensors (E ab and σ ab ) where ξ ab = e a (w) e b (w) − 3h ab . The peculiar velocities are connected to the magnetic Weyl tensor H ab where η abc = − √ −gε abcd u d is the Levi-Civitta antisymmetric volume form. The variables δ Ωρ , δ H and v a determine the inhomogeneity and anisotropy of the Szekeres-II regions through (33)-(34) in a coordinate independent manner. In fact, these quantities satisfy evolution equations that reduce in the linear limit to covariant dust perturbations in the comoving gauge [17].
It is important to remark that first order in v a does not (necessarily) imply that gradients of peculiar velocities are also small. The general conditions for self consistency of the linear approximation to peculiar velocities are presented in Appendix A.
6 Multiple ΛCDM matchings: Pancake models Szekeres-II models in general do not satisfy a strict Copernican principle at any scale. However, we can achieve an approximation to a Copernican principle by considering an array of Szekeres-II regions embedded and evolving jointly with a ΛCDM background as in the "pancake models" derived and discussed in [17]. In fact, in that paper we gave a quick illustrative example of Szekeres-II regions characterized by an energymomentum tensor like (18), with q a associated with CDM peculiar velocities as in (21), but with Λ = 0 (hence the FLRW background was an Einstein de Sitter model). We will consider in the following sections the case Λ > 0 under the framework of the pancake models of [17] with more detail and depth.
The junction conditions for these matching (evaluated at arbitrary fixed values of w) are where a(τ ) is the FLRW scale factor. Note that these matchings can be performed with a single (but arbitrary) FLRW background, but the Szekeres-II patches can be different and need not correspond to the same source (as long as conditions (35) hold at the matching hypersurfaces). In particular, we consider a pancake configuration in which the FLRW spacetime is a ΛCDM model matched with a Szekeres-II region at two hypersurfaces marked by constant w.
To comply with the junction conditions (35) we choose the following forms for the parameters in (12)-(15): -For the α parameters in (14) α 0 (w) = cos 2 (ν 0 w), with the remaining free functions given as C i sin 2 (ν i w), hence the matching interfaces are at w = 0, 2π/ν 0 , while the arbitrary constants C i must comply with 0 < C i < 1 in order to fulfill the compatibility conditions (see Appendix A). -For the β parameters in (15) the compatibility conditions and v a v a 1 require β 2 i 1 and β i β j 1 with i, j = 0, 1, 2, 3, hence these free functions must have the form β i (w) = i ϕ i (w), with ϕ i (w) bounded functions and the constants i complying with i j 1. -In the analytic form (27) of S(τ ) for the ΛCDM regions we define present cosmic time from S(τ 0 ) = 1 leading to τ 0 = 0.9662. We select the values Ω m 0 = 0.3 and Ω Λ 0 = 0.7 so that present dayΩ 0 = 1. The present day peculiar velocity in (28)-(29) and the density and Hubble scalar contrasts in (30)-(31) become -The function F follows from the numerical solution of (17), with initial conditions given by F ,τ (τ 0 , w) = 0 and F (τ 0 , w) = 0 ϕ 0 (w), where 0 is a constant that satisfies 0 < 0 < 1, while ϕ 0 (w) is a sinusoidal function complying with the form stated in the specification of the β parameters in (15) discussed above. We selected initial conditions at present cosmic time only as a matter of convenience, as they can be chosen at any fixed τ . Since (17) is a second order linear ODE, its general solution must be of the form F = ϕ + (w)F + (τ, w)+ϕ − (w)F − (τ, w), guaranteeing that initial conditions F (τ, w i ) = F ,τ (τ, w i ) = 0 can always be fulfilled by fixed values w = w i at any value of τ .
It is worth commenting that X has a very weak dependence on the angular coordinate φ, thus it is possible to get a robust notion of all quantities in terms of a single representative angle. After various trials, we choose the following numerical values for constant parameters: 1 = 0.7, 2 , 3 = C c1 = C c0 = 0.001, while the rest of the constants C i were taken as random numbers such that 0 < C i < 1, with numerical trials showing very weak dependence on the choice of these numbers. In general, it is necessary to test numerically values for the constants i , C ci and forms or the functions c i to avoid shell crossings and to obtain peculiar velocities whose magnitudes comply with the range of values for peculiar velocities of large scale structures found in the literature.

Results
All plots in figures 1-4 correspond to τ 0 = H 0 t 0 fixed at present time S(τ 0 ) = 1. Plotted quantities are displayed as functions of w and r normalized with the Hubble length. As mentioned before, w = 0 and w = π mark the matching hypersurfaces between the Szekeres-II region and the ΛCDM background, with the the Szekeres-II region encompassing the range 0 ≤ w ≤ π with its edges separated a comoving distance of ∼ 1.5 times the Hubble radius at τ 0 . The matching hypersurfaces are depicted in 1b, 2b, 3b, 4b by bold vertical lines, with the ΛCDM background extending for w < 0 and w > π. The horizontal red line depicts the zero contrast level corresponding to the ΛCDM background. Figure 1 displays the profile of the present day density contrast [δ Ω ] 0 given by (36) as a function of r (panel (a)), w (panel (b)), both normalized by the Hubble length, and φ (panel (c)). The panels reveal nearly the same value density contrast [δ Ω ] 0 ∼ 0.4 with respect to the ΛCDM background in the directions of r and w, with a very weak dependence on φ. While the Szekeres-II regions are clearly inhomogeneous and anisotropic for every observer, these graphs show relatively small local variations of [δ Ω ] 0 among observers. In fact, the free parameters allow to adjust the scale variation of the density contrast depending on a desired set of limits. Figure 2 depicts the profiles of the present day contrast of the Hubble scalar [δ H ] 0 given by (36), as a function of r (constant w, φ), w (constant r, φ) and φ (with constant r, w). As with the density contrast, the contrast of the Hubble scalar shows small local variation in different directions, as well as weak angular depen-dence, thus allowing for a a to controlled description of a desired level of inhomogeneity and anisotropy. Figure 3 displays from (23) the present day values of whereΘ 0 /3 and Θ 0 /3 are the present day Hubble scalars in the comoving and non-comoving frames, as a function of r (constant w, φ) and w (constant r, φ) (we omit plotting [ϑ/3] 0 as a function of φ because angular dependence is similar to that of figures 1 and 2). Both panels show that ϑ 0 /H 0 ∼ 10 −4 , a ratio three orders of magnitude below the 10 % difference observed in the H 0 tension. These values are consistent with nonrelativistic peculiar velocities and with the compatibility conditions in A. Figure 4 displays the present day radial velocity [v r ] 0 given by (37) as a function ofr (with w constant), w (with constant r) and φ (with constant r, w) is depicted in figure 4. Numerical values of the velocities are fractions of c, with their magnitude in the expected range |v r | < 800 km/s. As with the contrasts [δ Ωρ ] 0 and [δ H ] 0 , the radial velocities have similar values in different directions and very weak angular dependence.
We examine the time variation of the radial profile (for fixed w and φ) of the density and Hubble scalar contrasts (30)-(31) in figure 5. Both of these contrasts take near constant shapes that steadily decrease from z = 2 to their present values.

Structure formation
Following a careful parameter selection it is possible to obtain configurations free from shell crossings at least up to scales within the Hubble horizon, though some parameter combinations lead to divergent peculiar velocities even without shell crossings. This divergent behavior and the shell crossings signal the limit of validity of the description of CDM as dust. Thus, we restrict the parameters of the models to |v r | < 0.01c, a reasonable range of validity for structure formation involving nonrelativistic conditions, so that spacetime points where this bound is violated can be associated with the onset of virialization whose proper description is beyond the scope of these models.
After several numerical trials we found how to set up the free in order to control the placing the locus marking the beginning of shell crossings (and divergent peculiar velocities) at specific spatial positions and cosmic times (measured by redshifts of the ΛCDM region). In particular, there is sufficient parameter freedom to describe structure formation scenarios in which the onset   of virialization takes place at redshift values compatible with observations [28], for example, with |v r | → 0.1c at z ≈ 3. Considering redshifts in the ΛCDM region given by S(t) = 1/(1 + z) we plot v r in figure 6 as function of r and z. Notice that the velocities tend to increase their magnitude with increasing z.
We illustrate how peculiar velocities can become larger than the bound |v r | < 0.01c for z < 3 by plotting the time evolution of v r in figure 6. Panel (a) depicts the profile of v r as function of r (w constant) for z = 0, 1, 2, 3. The limit velocity 0.01c is reached at z = 3, marking the onset of virialization. Panel (b) depicts the profile of v r as function of z with w constantes for various values of r. Again, the onset of virialization occurs at z = 3.
Finally, we examine in figure 7 the time evolution of ϑ/3H 0 = (Θ − Θ)/3H 0 defined in (23) normalized by H 0 . This difference between the Hubble scalar in the comoving and non-comoving frames remains small in consistency with the compatibility conditions. However, this difference starts becoming large from z = 3 onwards which sets the limits of validity of the nonrelativistic approximation of peculiar velocities relating the two frames. However, we can argue that at this point the models cease to be valid marking the onset of virialization.
From the locus of the shell crossing in the example displayed in figure 6 we estimated the approximate conserved mass of the structure undergoing virialization from where the integral was kept as z i = 3 and v r ≤ 0.1c was satisfied along the integration domain. The estimated mass roughly corresponds to a galactic cluster whose onset of virialization at z = 3 is plausible.

Final discussion and conclusions
We have found for the Szekeres-II models under consideration an appealing physical interpretation as models that describe CDM and dark energy modeled as a Λ term with the novelty of incorporating peculiar velocities v a = q a /ρ for a non-comoving CDM source, all this in the context of appealing "pancake models" of cosmological inhomogeneities described by regions of Szekeres-II solutions embedded by smooth matchings to a ΛCDM background, introduced in previous work [17].    We have also provided a complementary view to previous work looking at the effects of cosmological sources (for example baryons and CDM) evolving along different 4-velocity frames [10]). In order to illustrate the effects of local inhomogeneity and anisotropy brought by the models we compared their dynamical variables with their values in the ΛCDM background. For this purpose, we considered a configuration made of a single Szekeres-II region extending 1.5 times the Hubble radius in the w direction, smoothly matched to a ΛCDM background on both extremes. We obtained (see figures 1, 2) by numerical integration of the field equations the present cosmic time contrasts respect to this background of the density and Hubble scalar ([δ Ωρ ] 0 and [δ H ] 0 from (36)), also at different cosmic times (see figure 5). In all quantities the variation with respect to the angular coordinate was very weak, thus identifying an anisotropy based on differences along (essentially) two directions: r and w (as suggested by looking at the metric (7) in rectangularlike coordinates x, y instead of r, φ).
Figures figures 1, 2 reveal present cosmic time values of density and Hubble scalar contrasts within the Szekeres-II region respectively varying from zero to maximal values of 0.4 and 0.1 roughly in the same pattern along both directions r and w. We tested various combinations of initial conditions and found roughly the same variation patterns with different maximal values, thus indicating a relatively mild deviation of from local isotropy that can be controlled by suitable choices of free parameters at least in scales up to the Hubble radius.
The behavior of present day radial peculiar velocity [v r ] 0 is displayed in figure 4. Notice how [v r ] 0 vanishes at r = 0 and approximately increases linearly with r for all w within the Szekeres-II region, while in the w direction it takes larger values as r increases (vanishing as required by junction conditions at the w values marking the matching interface). This pattern illustrates how r = 0 for varying w denotes the coordinate locus of privileged observers, analogous to observers along the symmetry center of spherical symmetry. However, this is not a mere coordinate effect, as the curve along r = 0 parametrized by w is a spacelike geodesic and Killing vector of the hypersurfaces orthogonal to the comoving 4-velocity (see [17]).
We found free parameter choices that lead to an evolution free from shell crossings with peculiar velocities remaining in the non-relativistic regime |v r | < 0.01c at least in scales of the order of the Hubble radius in the main directions r and w. However, we also found free parameter combinations that lead to shell crossings around z = 3, with peculiar velocities growing and even diverging, thus identifying these spacetime points as marking the onset of virialization when a model based on a dust description of CDM is no longer valid. We estimated a CDM mass of ∼ 10 15 contained in a region associated with these shell crossings that can be identified with a large galactic cluser.
We fully acknowledge the limitations of the models we have studied in this paper: they are basically toy models of inhomogeneities in a ΛCDM background that are valid only in the scales and cosmic times in which CDM can be modeled as dust. The novelty of our approach (with respect to previous usage of Szekeres models in this context) is that we consider the class Szekeres-II and that CDM is not comoving with the frame associated with the CMB and the ΛCDM background. Evidently, these toy models cannot describe a highly complex process like virialization, but we can assume that shell crossings (which are a generic feature) can mark the limit of validity of the models due to the onset of this process.
Nevertheless, we believe that these toy models have a valuable potential for cosmological applications: first, they allow to study the observed non-relativistic peculiar velocities in the framework of an exact solution of General Relativity, thus potentially contributing to improve our understanding of the role of peculiar velocities in cosmic dynamics, not only at local deep subhorizon scales, but even at scales comparable to the Hubble horizon. Also, a better understanding of the dynamics of peculiar velocities and Hubble flows from different congruences of observers can contribute to address the H 0 tension. Finally, the models can serve as an exact solution to probe numerical codes in the emerging field of numerical relativistic cosmology.

A Compatibility conditions
For the sake of completeness we include the compatibility conditions we presented in [17].
The energy-momentum tensor of a dust source in the frame of a non-comoving observer,û a iŝ Using the decomposition (1), we considered p = 0, Π ab = 0 and searched under what conditions, in the limit v a v a → 0, T ab =T ab . To order zero, lim v a v aT ab = ρ u a u b + 2u (a v b) , obtaining ρu (a v b) = u (a q b) . Even though we take v a v a 1, this does not imply the derivatives are small, so we must search conditions to first order. As both energy-momentum tensors are conserved, we obtain the second condition when we obtain an identity from this equation 0 = 0 in the limit v a v a → 0. From the previous conditions, and considering the derivatives of the γ factor: From (40)-(42) and the zero order relations it is straightforward to verify that Neglecting quadratic terms on v a we obtain the second condition: (v a v b ) ;b = 0.
This implies v a v b is constant, which we take as v a v a << 1 for consistency with our initial hypothesis v a v a << 1. Therefore our conditions for compatibility are With these considerations the energy conservation u a T ab ;b is proportional to (H 0 /c) 3 which justifies our approximation.

B The FLRW limit
The form of the line element presented in [11] is where X = e ν [A + BQ] + F , and e ν = (1 + k(x 2 + y 2 )) −1 . The perfect fluid case arises for B = 0, [11], this election can be made considering the form of B given by Goode, B = β 3 (w)r 2 + β 2 (w)r cos θ + β 1 (w)r sin θ + β 0 , by taking β 3 (w) = h 2 (w) = β 1 (w) = β 0 = 0. Krasinski presents the Szekeres-II with a perfect fluid source with the following line element ds 2 = −dt 2 + e 2α dw 2 + e 2β (dx 2 + dy 2 ), where e β = e ν S(t), e α = λ(t, w) + S(t)Σ(x, y, w), and λ a function determined by a differential equation. From the form of the functions stated in [11] and [15] it is easy to show that λ(t, w) = S(t)F (t, w), while Σ = Ae ν . Therefore, the condition stated by Krasinki at the end of section 2.1.1 holds. As mentioned in the introduction there is no natural FLRW limit for this model, quoting [15] "the FLRW limit results unnaturally in this subfamily: the additional symmetries of the FLRW models appear from nowhere." Further comments on how the FLRW models arise are found in B. It is worth noting that the FLRW limit is a limit in the space of parameters and not a limit in the manifold or an extension of the manifold itself.