# Relativistic approach to the kinematics of large-scale peculiar motions

Eleni Tsaprazi<fnref n="1" />,<fnref n="2" />,<fnref n="4" />, Christos G. Tsagas<fnref n="1" />,<fnref n="3" />,ᵃ

Received: 28 May 2020 / Accepted: 31 July 2020 / Published online: 19 August 2020
© The Author(s) 2020

**Abstract** We consider the linear kinematics of large-scale peculiar motions in a perturbed Friedmann universe. In so doing, we take the viewpoint of the “real” observers that move along with the peculiar flow, relative to the smooth Hubble expansion. Using relativistic cosmological perturbation theory, we study the linear evolution of the peculiar velocity field, as well as the expansion/contraction, the shear and the rotation of the bulk motion. Our solutions show growth rates considerably stronger than those of the earlier treatments, which were mostly Newtonian. On scales near and beyond the Hubble radius, namely at the long-wavelength limit, peculiar velocities are found to grow as $a^2$, in terms of the scale factor, instead of the Newtonian $a^{1/2}$-law. We attribute this to the fact that, in general relativity, the energy flux, triggered here by the peculiar motion of the matter, also contributes to the local gravitational field. In a sense, the bulk flow gravitates, an effect that has been bypassed in related relativistic studies. These stronger growth-rates imply faster peculiar velocities at horizon crossing and higher residual values for the peculiar-velocity field. Alternatively, one could say that our study favours bulk peculiar flows larger and faster than anticipated.

## 1 Introduction

Large-scale peculiar motions, also referred to as “bulk flows”, are an established observational fact, confirmed by many surveys extending out to scales of several hundred Mpc (e.g. see [1] and references therein). Peculiar velocities are predicted by all structure formation scenarios, as the direct and unavoidable outcome of the ever increasing inhomogeneity and anisotropy of the post-recombination universe.

The subject has a fairly long research history that goes back several decades. Nevertheless, although there are many structure-formation studies that incorporate peculiar velocities, essentially all the work that focuses on the evolution of the peculiar-velocity field is Newtonian, or quasi-Newtonian, in nature (see [2–4] and [5,6] respectively). In addition, to the best of our knowledge, all these studies are conducted in the rest-frame of the smooth Hubble expansion. The latter is defined as the coordinate system where the dipole in the Cosmic Microwave Background (CMB) spectrum vanishes. However, even at the linear level, there are subtle differences between the Newtonian and the relativistic treatments, given the very different way the two theories approach issues as fundamental as the time-vs-space relation and the nature of gravity itself. Moreover, no real observer in the universe follows the CMB frame, but we all move relative to it. Our Local Group of galaxies, for example, “drifts” with respect to the smooth Hubble flow at approximately 600 km/sec. For these reasons, the present work looks into the question of large-scale peculiar velocities and of their evolution by employing relativistic cosmological perturbation theory and by adopting the view point of observers living in typical galaxies (like our Milky Way) and moving relative to the smooth Hubble expansion. Put another way, all our calculations are done in a coordinate system moving with respect to the CMB frame. Our aim is to identify possible differences between the two approaches and thus provide a better theoretical understanding of the peculiar-motion kinematics, especially in view of the anticipated data from upcoming bulk-flow surveys. Note that the linear nature of our analysis implies that our scales of interest are large enough to neglect all nonlinear effects, which in practice means scales in excess of 100 Mpc.

We begin with a brief introduction to the analysis of peculiar velocities within the framework of relativistic cosmology in Sect. 2. In so doing, we assume a perturbed Friedmann-

<footnote n="1">¹ Section of Astrophysics, Astronomy and Mechanics, Department of Physics, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece</footnote>

<footnote n="2">² Sorbonne Universités, UPMC Univ Paris 06, CNRS, Institut d'Astrophysique de Paris (IAP), 98bis Boulevard Arago, Paris 75014, France</footnote>

<footnote n="3">³ DAMTP, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote n="4">⁴ Present address: Department of Physics, The Oscar Klein Centre, Albanova University Center, Stockholm University, 106 91 Stockholm, Sweden</footnote>

<footnote>ᵃ e-mail: tsagas@astro.auth.gr (corresponding author)</footnote>
Robertson-Walker (FRW) universe filled with a pressureless fluid that can be baryonic, or low-energy cold dark matter (CDM), or both. We also adopt the 1 + 3 covariant and gauge invariant approach to cosmological perturbations (see [8,9] for recent reviews) and allow for two families of observers. These are the “fictitious” ideal observers that follow the CMB-frame and their “real” counterparts moving relative to it. The latter represent typical observers in our universe, living in galaxies like our Milky Way, which makes their coordinate system the natural frame to use in cosmological studies.

In addition to the peculiar velocity field itself, we investigate the full spectrum of the peculiar kinematics, namely the volume expansion/contraction, the shear and the rotation of the bulk motion. We establish the linear sources of cosmological peculiar motions in Sect. 3. As expected, we find that such flows are triggered by the growing inhomogeneity of the post-recombination universe. Our analysis also provides the full set of differential equations monitoring the peculiar velocity field, its expansion/contraction-rate its shear distortion and its rotation after decoupling. The homogeneous parts of these formulae accept analytic (power-law) solutions. These show increase, both in absolute terms and relative to the background Hubble expansion (see Sect. 4). More specifically, the peculiar velocity ($\tilde{v}$) is found to grow as $\tilde{v} \propto a^2$, where $a = a(t)$ is the cosmological scale-factor. This implies that the strength of the peculiar flow relative to the universal expansion, described by the ratio $\tilde{v}/v_H$ – with $v_H$ representing the Hubble velocity on the corresponding length, grows as $\tilde{v}/v_H \propto a^{5/2}$. At the same time, the expansion/contraction-rate ($\tilde{\vartheta}$) of the bulk flow shows a scale-dependent growth, which reaches its maximum ($\tilde{\vartheta} \propto a$) on super-Hubble scales and tends to a constant (i.e. $\tilde{\vartheta} \to \text{constant}$) inside the horizon. Finally, the peculiar shear and the peculiar vorticity ($\tilde{\zeta}$ and $\tilde{\omega}$ respectively) are found to grow as $\tilde{\zeta}, \tilde{\omega} \propto a$.

To the best of our knowledge, the growth rates given above are considerably stronger than those previously reported in the theoretical literature. The Newtonian/quasi-Newtonian treatments, in particular, lead to $v \propto a^{1/2}$ in the CMB frame [2–6], which implies $v/v_H \propto a$ after decoupling. The same result was also recently obtained by applying the Newtonian version of the 1 + 3 formalism to the bulk-flow frame [7]. The considerable difference between the relativistic and the Newtonian results stems from the fundamentally different way the two theories treat the gravitational field. More specifically, in relativity, the energy flux contributes to the stress-energy tensor of the matter and therefore to the local gravitational field (e.g. see [8,9]). In our case, this purely relativistic effect ensures a flux-contribution to the stress-energy tensor due to the (peculiar) motion of the matter. Put another way, the bulk flow itself gravitates. This, in turn, feeds into the conservation laws and eventually appears in the equations monitoring cosmological perturbations and, more specifically, the evolution

of peculiar-velocity perturbations (see Sects. 2.4 and 3.1 below). To our best knowledge, the relativistic effect reflecting the bulk-flow flux contribution to the local gravitational field, has been typically bypassed in related relativistic studies and the reader is referred to Sect. 3.2 for further discussion and some characteristic examples.

Isolating and solving only the homogenous parts of the relativistic differential equations means that the resulting solutions are expected to hold on large enough scales, where the inhomogeneous components (consisting of spatial gradients) are negligible. With this in mind, the cautious approach should be to apply our power-law solutions on lengths near and beyond the Hubble radius, but not well inside the horizon. It may turn out that our long-wavelength results reduce to their Newtonian counterparts inside the Hubble scale, in which case the relativistic analysis could provide the “initial conditions” for the subsequent Newtonian treatments. In any case, stronger growth rates for the peculiar-velocity field on super-Hubble lengths, implies faster peculiar velocities at horizon-crossing, which in turn suggests higher residual values today. Therefore, at this stage, it is fair to say that the relativistic analysis seems to favour a number of surveys reporting bulk peculiar flows larger and faster than it is generally anticipated [10–16].

## 2 The peculiar velocity field

In relativity neither time nor space retain their absolute Newtonian notion. Instead, observers moving with respect to each other have their own measure of time and space and they generally experience different versions of what we call “reality”.

### 2.1 The 4-velocity “tilt”

Let us consider two families of observers with (timelike) 4-velocities $u_a$ and $\tilde{u}_a$ respectively. Let us also assume that the latter family has peculiar velocity $\tilde{v}_a$ relative to the former. The three aforementioned velocity fields are then related by the familiar Lorentz boost

$$ \tilde{u}_a = \tilde{\gamma}(u_a + \tilde{v}_a), \qquad (1) $$

where $u_a u^a = -1 = \tilde{u}_a \tilde{u}^a$, $u_a \tilde{v}^a = 0$, $\tilde{\gamma} = 1/\sqrt{1 - \tilde{v}^2}$ and $\tilde{v}^2 = \tilde{v}_a \tilde{v}^a < 1$ (see Fig. 1).<fnref n="1" /> The “tilt” between the two 4-velocity vectors is determined by the (hyperbolic) angle $\beta$, defined by $\cosh \beta = -\tilde{u}_a u^a = \tilde{\gamma}$. The latter ensures that $\beta = \ln(\tilde{\gamma} + \sqrt{\tilde{\gamma}^2 - 1})$, with $\beta > 0$ since $\tilde{\gamma} > 1$ (see [17] for further discussion and additional technical details).

<footnote n="1">¹ Throughout this manuscript Latin indices run from 0 to 3, while Greek ones take values from 1 to 3. Our metric tensor ($g_{ab}$) has signature $[-1, 1, 1, 1]$ and we have also set $c = 1$.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1 Observers ($O_1$ and $O_2$) living inside the bulk-flow domain (D) and moving with (local) peculiar velocity $\tilde{v}_a$ relative to the Hubble expansion. The 4-velocities $u_a$ and $\tilde{u}_a$ define the CMB frame and the coordinate system of the peculiar motion respectively (see Eq. (1))</figcaption></figure>

Note that, when dealing with large-scale bulk peculiar motions, the local and the mean velocities of the flow ($\tilde{v}$ and $\tilde{V}$ respectively) are typically related via the integral $\tilde{V} = (3/4\pi r^3) \int_{x<r} \tilde{v}dx^3$, where $\tilde{v}^2 = \tilde{v}_a \tilde{v}^a$ and $r$ is the (effective) radius of the bulk flow.

## 2.2 The 1 + 3 threading

Introducing two 4-velocities means that there are two temporal directions (along $u_a$ and $\tilde{u}_a$ respectively) and two 3-dimensional spatial sections (orthogonal to $u_a$ and $\tilde{u}_a$). Projecting onto these 3-spaces is achieved by using the symmetric projection tensors

$$h_{ab} = g_{ab} + u_a u_b \quad \text{and} \quad \tilde{h}_{ab} = g_{ab} + \tilde{u}_a \tilde{u}_b, \qquad (2)$$

with $h_{ab}u^a = 0 = \tilde{h}_{ab}\tilde{u}^b$ and $h_a^a = 3 = \tilde{h}_a^a$ [8,9]. Note that, when there is no vorticity, these two projectors also act as the metric tensors of their associated 3-dimensional hypersurfaces.

We may now proceed to define temporal and spatial differentiation relative to the two timelike coordinate systems introduced in Sect. 2.1 before. In particular, the time-derivatives in the $u_a$ and the $\tilde{u}_a$ frames are denoted by

$$\cdot = u^a \nabla_a \quad \text{and} \quad \cdot' = \tilde{u}^a \nabla_a, \qquad (3)$$

respectively. The corresponding spatial gradients, on the other hand, are

$$D_a = h_a^b \nabla_b \quad \text{and} \quad \tilde{D}_a = \tilde{h}_a^b \nabla_b. \qquad (4)$$

Using the above, one can decompose any spacetime variable, operator and equation into their timelike and spacelike parts (relative to the $u_a$ or the $\tilde{u}_a$ field), thus achieving an 1 + 3 threading of the host spacetime into time and space [8,9].

## 2.3 Kinematic decomposition

The irreducible kinematic variables of the two 4-velocity fields emerge by decomposing their gradients. More specifically, we have [18]

$$\nabla_b u_a = \frac{1}{3} \Theta h_{ab} + \sigma_{ab} + \omega_{ab} - A_a u_b. \qquad (5)$$

In the above $\Theta = D^a u_a$, $\sigma_{ab} = D_{[b}u_a)$, $\omega_{ab} = D_{[b}u_a]$ and $A_a = \dot{u}_a$ are respectively the expansion/contraction scalar, the shear tensor, the vorticity tensor and the 4-acceleration vector of the $u_a$-field.<fnref n="2" /> Positive values for $\Theta$ mean expansion, while in the opposite case we have contraction. Nonzero shear implies changes in the shape (under constant volume) of the associated fluid element, non-vanishing vorticity ensures rotation, while nonzero 4-acceleration indicates the presence of non-gravitational forces. In an exactly analogous way we may write

$$\nabla_b \tilde{u}_a = \frac{1}{3} \tilde{\Theta} \tilde{h}_{ab} + \tilde{\sigma}_{ab} + \tilde{\omega}_{ab} - \tilde{A}_a \tilde{u}_b, \qquad (6)$$

for the gradient of the $\tilde{u}_a$-field. Note that, in cosmological studies, the volume scalar defines the scale factor of the universe by means of $\dot{a}/a = \Theta/3$. Also, the shear, the vorticity and the 4-acceleration are all spacelike, namely $\sigma_{ab}u^b = 0 = \omega_{ab}u^b = A_a u^a$ by construction, with exactly analogous constraints applying to their “tilded” counterparts (i.e. $\tilde{\sigma}_{ab}\tilde{u}^b = 0 = \tilde{\omega}_{ab}\tilde{u}^b = \tilde{A}_a\tilde{u}^a$).

No realistic fluid-flow is absolutely rigid, but instead it is expected to expand or contract, to change shape and to rotate, even by small amounts. It is therefore plausible to argue that large-scale bulk peculiar motions should behave in a similar manner. The expansion/contraction-rate, the shear distortion and the rotation of the $\tilde{v}_a$-field, technically speaking the irreducible variables of the peculiar kinematics, are obtained by splitting its spatial gradient. More specifically, taking the view point of the tilded observer, we have [19,20]

$$\tilde{D}_b \tilde{v}_a = \frac{1}{3} \tilde{\Theta} \tilde{h}_{ab} + \tilde{\zeta}_{ab} + \tilde{\omega}_{ab}, \qquad (7)$$

with $\tilde{\Theta} = \tilde{D}^a \tilde{v}_a$, $\tilde{\zeta}_{ab} = \tilde{D}_{[b} \tilde{v}_a)$ and $\tilde{\omega}_{ab} = \tilde{D}_{[b} \tilde{v}_a]$ representing the peculiar expansion/contraction, the peculiar shear and the peculiar vorticity respectively.<fnref n="3" />

<footnote n="2">² We use round brackets for symmetrisation, square for antisymmetrisation and angled brackets to denote the symmetric and trace-free part of second-rank tensors. Consequently, $\sigma_{ab} = D_{[b}u_a) = D_{(b}u_a) - (D^c u_c/3)h_{ab}$ by default.</footnote>

<footnote n="3">³ The scalar $\tilde{\Theta}$ monitors the volume expansion/contraction of the bulk motion. When $\tilde{\Theta} > 0$ the peculiar flow expands and in the opposite case it contracts. The peculiar shear follows changes in the shape of the bulk flow, say from spherical to ellipsoidal, while the vorticity contains information about its rotation. At the linear level $|\tilde{\Theta}|/\Theta \ll 1$, $\tilde{\zeta}/\Theta \ll 1$ and $\tilde{\omega}/\Theta \ll 1$, with $2\tilde{\zeta}^2 = \tilde{\zeta}_{ab}\tilde{\zeta}^{ab}$ and $2\tilde{\omega}^2 = \tilde{\omega}_{ab}\tilde{\omega}^{ab}$.</footnote>
## 2.4 Linear relations between the two frames

So far we have not imposed any constraints on the peculiar velocity, which means that our definitions and our formulae hold for arbitrarily fast relative motions in a general spacetime. Hereafter, we will only consider non-relativistic drift velocities with $\tilde{v}^2 \ll 1$ and

$$ \tilde{u}_a \approx u_a + \tilde{v}_a, \qquad (8) $$

since $\tilde{\gamma} \simeq 1$.<fnref n="4" /> Treating the peculiar velocity field as a perturbation, we assume that the host spacetime is an almost-FRW universe. Finally, we identify the $u_a$-field with the coordinate system of the Hubble flow and place the tilded observers in a typical galaxy like our Milky Way.<fnref n="5" /> Then, the kinematic variables defined in Sect. 2.3 are related by the linear expressions [5,6]

$$ \tilde{\Theta} = \Theta + \tilde{\vartheta}, \quad \tilde{\sigma}_{ab} = \sigma_{ab} + \tilde{\zeta}_{ab}, \quad \tilde{\omega}_{ab} = \omega_{ab} + \tilde{\varpi}_{ab} \quad (9) $$

and

$$ \tilde{A}_a = A_a + \tilde{v}'_a + \frac{1}{3} \Theta \tilde{v}_a, \qquad (10) $$

where $\tilde{v}'_a = \tilde{u}_a^b \nabla_b \tilde{v}_a$ is the time derivative of the peculiar velocity in the tilded frame (not to be confused with the “peculiar gravitational acceleration” of the Newtonian treatments). Of particular interest for our purposes is Eq. (10). According to this relation, we cannot set both 4-acceleration vectors to zero simultaneously. When $A_a$ vanishes in the CMB frame, for example, the tilded observers measure a nonzero 4-acceleration solely because of their peculiar motion (i.e. $\tilde{A}_a = \tilde{v}'_a + (\Theta/3)\tilde{v}_a$ when $A_a = 0$).

In a similar manner, one can show that the dynamical variables measured in the two frames are related by

$$ \tilde{\rho} = \rho, \quad \tilde{p} = p, \quad \tilde{q}_a = q_a - (\rho + p)\tilde{v}_a \qquad (11) $$

and

$$ \tilde{\pi}_{ab} = \pi_{ab}, \qquad (12) $$

to first approximation [5,6]. Here, $\rho$ is the energy density, $p$ is the isotropic pressure, $q_a$ is the energy flux and $\pi_{ab}$ is the viscosity of the matter as measured in the $u_a$-frame, while

their tilded counterparts are associated with the $\tilde{u}_a$-field. Expression (11c) has special significance in this study, since it ensures an additional energy-flux contribution as a result of relative motion alone. In general relativity, the energy flux gravitates as well, since it contributes to the stress-energy tensor of the matter. Therefore, the peculiar motion of the matter also contributes to the local gravitational field. In a sense, the bulk flow itself gravitates.

In what follows we will use the above linear relations to study the kinematic evolution of large-scale peculiar motions in a perturbed Friedmann universe filled with pressureless matter (baryonic or/and CDM).

# 3 Linear sources of peculiar flows

Large-scale peculiar motions are treated as a result of the increasing inhomogeneity and anisotropy of our universe, due to the ongoing structure-formation process. The latter starts in earnest after recombination, once the baryons have decoupled from the background radiation field.

## 3.1 Peculiar velocities

Let us assume an almost-FRW cosmology filled with pressureless matter. This can be baryonic, or low-energy CDM, or a mixture of both. The assumption of an unperturbed FRW background, ensures that all perturbations (including those in the peculiar-velocity field) vanish there. This makes our analysis gauge invariant and therefore free of any gauge-related issues, in line with the Stewart & Walker lemma [21]. Also, the absence of pressure means that we can set $A_a = 0$ in the CMB frame.<fnref n="6" /> Then, relation (10), which also serves as the linear propagation equation of the peculiar velocity in the bulk-flow frame, recasts as

$$ \tilde{v}'_a = -H\tilde{v}_a + \tilde{A}_a, \qquad (13) $$

with $H = \Theta/3$ being the Hubble parameter in the CMB frame. Therefore, at the linear level, the sole source of peculiar velocities is the 4-acceleration. This makes $\tilde{A}_a$ the key to the subsequent evolution of the $\tilde{v}_a$-field. It also means that the way the 4-acceleration is treated is crucial. In the Newtonian analogue of Eq. (13), the role of the 4-acceleration is played by the gradient of the gravitational potential (e.g. see [2–4,7]). Similarly, the quasi-Newtonian approach (which

<footnote n="4">⁴ Typical Newtonian studies, define peculiar velocities by introducing physical and comoving coordinates ($r^\alpha$ and $x^\alpha$ respectively, with $\alpha = 1, 2, 3$), related by $r^\alpha = ax^\alpha$. The time derivative of the latter leads to $v^\alpha = v_H^\alpha + v_\rho^\alpha$, where $v^\alpha = \dot{r}^\alpha$ is the physical (total) velocity, $v_H^\alpha = Hr^\alpha$ is the Hubble velocity and $v_\rho^\alpha = a\dot{x}^\alpha$ is the peculiar velocity. The above given relation between the three velocity fields is the Newtonian analogue of Eqs. (1) and (8). Having said that, we remind the reader that the 4-velocities $\tilde{u}_a$ and $u_a$ are both timelike vectors, whereas their Newtonian counterparts are all purely spatial.</footnote>

<footnote n="5">⁵ Although relativity postulates the absence of preferred coordinate systems, the universal expansion naturally selects the CMB frame as the reference system relative to which peculiar velocities should be defined and measured.</footnote>

<footnote n="6">⁶ Setting $A_a = 0$ means that the worldlines of the CMB frame are timelike geodesics, whereas those of the tilted observers are not. In principle, one could also set $\tilde{A}_a = 0$, in which case $A \neq 0$ (e.g. see [5,6]). Then, it is the tilted observers that move along timelike geodesics, while their Hubble-flow counterparts follow non-geodesic lines. Here, we are assuming that the universe is an FRW model relative to the CMB frame, which sets $A_a = 0$ by default. That aside, one can perform the analysis in either coordinate system and reach the same results.</footnote>
was the first – to the best of our knowledge – 1+3 covariant study of the issue) introduces an effective gravitational potential to account for the effects of the 4-acceleration [5,6]. The latter, however, pre-assumes that the perturbed spacetime is both irrotational and shear-free. Moreover, the necessary propagation formula of the 4-acceleration was obtained after introducing an ansatz for the time evolution of the aforementioned potential [5,6]. Here, we will take an alternative route. More specifically, by applying relativistic linear cosmological perturbation theory to the tiled frame, we will obtain analytical expressions for both the 4-acceleration vector and its time-derivative.

## 3.2 The 4-acceleration

Our starting point is the linear transformation law (11c), which guarantees that, even when the cosmic medium appears as a perfect fluid in the CMB frame, there is an energy-flux vector in the tiled frame solely due to the latter’s relative motion. More specifically, following (11c), we deduce that $\tilde{q}_a = -\rho\tilde{v}_a$ when $q_a = 0$. Consequently, there is a flux-contribution to the stress-energy tensor triggered by the peculiar motion of the matter, which then feeds into the energy and the momentum conservation laws and eventually reaches the evolution formulae of cosmological perturbations. In particular, when linearised in the tiled frame, the evolution formula of density inhomogeneities (e.g. see Eqs. (2.3.1) and (10.101) in [8,9]), reads [22]

$$ \tilde{\Delta}'_a = -\tilde{Z}_a + \frac{3aH}{\rho} (\tilde{q}'_a + 4H\tilde{q}_a) - \frac{a}{\rho} \tilde{D}_a \tilde{D}^b \tilde{q}_b. \quad (14) $$

Here, $\tilde{\Delta}_a = (a/\rho)\tilde{D}_a\rho$ and $\tilde{Z}_a = a\tilde{D}_a\tilde{\Theta}$ to first approximation, representing inhomogeneities in the matter distribution and in the universal expansion respectively [8,9]. Employing relations (10) and (11c), with $A_a = 0 = q_a$ and $p = 0$, it is straightforward to show that

$$ \tilde{A}_a = -\frac{1}{\rho}(\tilde{q}'_a + 4H\tilde{q}_a). \quad (15) $$

Keeping in mind that $\tilde{q}_a = -\rho\tilde{v}_a$, the latter combines with Eq. (14) to give

$$ \tilde{A}_a = \frac{1}{3H} \tilde{D}_a \tilde{\vartheta} - \frac{1}{3aH} (\tilde{\Delta}'_a + \tilde{Z}_a). \quad (16) $$

Note that in a Newtonian perturbative study Eq. (14) reduces to $\tilde{\Delta}'_a = -\tilde{Z}_a$, with $\tilde{\Delta}_a = (a/\rho)\partial_\alpha\rho$ and $\tilde{Z}_a = a\partial_\alpha\tilde{\Theta}$ [23]. The absence of an acceleration term in the above, explains why expression (16) has no close Newtonian analogue. There, as well as in the quasi-Newtonian treatments, the acceleration and the 4-acceleration are given by the gradient of the gravitational potential. All these make (16) the key “relativistic correction” and the practical reason for the

differences between the relativistic and the Newtonian/quasi-Newtonian results.

The literature also contains a number of relativistic structure-formation studies. To the best of our knowledge, these investigations have not encountered the role of the bulk-flow flux in the evolution of peculiar-velocity perturbations, as reported here. The reasons vary and the following examples are indicative of that. Technically speaking, closer to our analysis are perhaps the multi-fluid scenarios, where the species involved have their individual peculiar velocities. These studies, however, are typically performed in the energy frame (otherwise known as Landau-Lifshitz frame), where the total flux vector is set to zero (i.e. the individual fluxes cancel each other out – e.g. see [8,9] and also [24]). Then, there is no flux contribution to the stress-energy tensor and the relativistic effects of the bulk-flow motion described here are bypassed.<fnref n="7" /> There are also non-covariant relativistic approaches in the literature. These are usually gauge-dependent, with some of them adopting the comoving gauge, where the peculiar velocities vanish by default (e.g. see [25]). Other treatments allow for peculiar-velocity perturbations, but their evolution is not the focal point of the study. As a result, no explicit solutions are provided (e.g. see [22,26]).<fnref n="8" /> Overall, there are relatively few analytical treatments focusing on the kinematics of cosmological peculiar motions. Moreover, essentially all the available studies are performed in the CMB and not in the tiled frame, which is by construction the natural coordinate system to analyse peculiar flows.

Before closing this section, we remind the reader that the 4-acceleration is the only source of peculiar-velocity perturbations (see Eq. (13) in Sect. 3.1). Given that $\tilde{\vartheta}$ vanishes in the absence of drift motions, expression (16) leaves $\tilde{\Delta}'_a$ and $\tilde{Z}_a$ as the sole sources of linear peculiar velocities. Therefore, as expected, the $\tilde{v}_a$-field is induced by the increasing inhomogeneity of the post-recombination universe and more specifically by temporal variations in the density-gradients and by spatial variations in the universal expansion (described by $\tilde{\Delta}'_a$ and $\tilde{Z}_a$ respectively).

<footnote n="7">⁷ Ours is a single-fluid study. Nevertheless, one can still recover the results of the multi-fluid analysis “phenomenologically”, by setting $\tilde{q}_a = 0$ in Eq. (15) and therefore $A_a = 0$ in (13). The latter then leads to $\tilde{v}_a \propto a^{-1}$ and to $\tilde{\vartheta} \propto a^{-2}$, as in [24] and [8,9] respectively. Taken at face value, these solutions imply that the peculiar flows decay quickly with the universal expansion on all scales. A result that does not seem to agree with the observations.</footnote>

<footnote n="8">⁸ Assuming that the direct comparison between Eq. (13) here and expression (5.6) in [26] is permitted, our 4-acceleration vector could be expressed as a combination of zero and first-order terms within the post-Friedmann approximation scheme.</footnote>
## 4 Linear evolution of peculiar flows

Expression (13) requires an evolution formula for the 4-acceleration to close the system of the differential equations. In what follows, we will do so and also obtain analytic solutions for the $\tilde{v}_a$-field.

### 4.1 The key differential equations

Taking the time derivative of Eq. (16), recalling that $\dot{H} = -H^2[1 + (\Omega/2)]$ in the FRW background (with $\Omega = \kappa\rho/3H^2$ being the associated density parameter) and using the linear commutation law $(\tilde{D}_a\tilde{\vartheta})' = \tilde{D}_a\tilde{\vartheta}' - H\tilde{D}_a\tilde{\vartheta}$ between the temporal and the spatial derivatives of first-order scalars [5,6], provides the propagation formula

$$ \tilde{A}'_a = \frac{1}{2} H \Omega (\tilde{v}'_a + H \tilde{v}_a) + \frac{1}{3H} \tilde{D}_a \tilde{\vartheta}' - \frac{1}{3aH} (\tilde{\Delta}''_a + \tilde{Z}'_a), \quad (17) $$

which monitors the linear evolution of the 4-acceleration relative to the bulk-flow frame. We therefore have analytic linear expressions for both the 4-acceleration and its time derivative, while avoiding the restrictions of the quasi-Newtonian analysis (see Sect. 3.1 previously).

Substituting the right-hand side of (17) into the time-derivative of Eq. (13) leads to

$$ \tilde{v}''_a = -H \left(1 - \frac{1}{2}\Omega\right) \tilde{v}'_a + H^2 (1 + \Omega) \tilde{v}_a + \frac{1}{3H} \tilde{D}_a \tilde{\vartheta}' - \frac{1}{3aH} (\tilde{\Delta}''_a + \tilde{Z}'_a). \quad (18) $$

In addition, using the linear commutation laws $\tilde{D}_b\tilde{v}'_a = (\tilde{D}_b\tilde{v}_a)' + H\tilde{D}_b\tilde{v}_a$ and $\tilde{D}_v\tilde{\vartheta}''_a = (\tilde{D}_b\tilde{v}_a)'' + 2H(\tilde{D}_b\tilde{v}_a)' - (H^2\Omega/2)\tilde{D}_b\tilde{v}_a$, the linearised 3-gradient of the above reads

$$ (\tilde{D}_b \tilde{v}_a)'' = -3H \left(1 - \frac{1}{6}\Omega\right) (\tilde{D}_b \tilde{v}_a)' + 2H^2 \Omega \tilde{D}_b \tilde{v}_a + \frac{1}{3H} \tilde{D}_b \tilde{D}_a \tilde{\vartheta}' - \frac{1}{3a^2 H} (\tilde{\Delta}''_{ab} + \tilde{Z}'_{ab}), \quad (19) $$

where $\tilde{\Delta}_{ab} = a\tilde{D}_b\tilde{\Delta}_a$ and $\tilde{Z}_{ab} = a\tilde{D}_b\tilde{Z}_a$. The last two differential equations govern the linear kinematics of peculiar motions, as seen by observers “living” inside these bulk peculiar flows. Both relations apply to all scales and hold in an almost-FRW universe with nonzero background curvature and a pressureless (baryonic or/and CDM) matter. In what follows, we will attempt to extract analytical solutions from these formulae.

### 4.2 Peculiar velocity

To this point, we have considered a perturbed Friedmann universe without imposing any constraints on its spatial curvature.

Hereafter, we will confine to an Einstein-de Sitter background by setting $\Omega = 1$ and $H = 2/3t$. However, Eqs. (18) and (19) do not accept analytic solutions even at the $\Omega = 1$ limit. This is not uncommon in analytical cosmological studies and the standard way around it is to focus on the so-called long-wavelength solutions. In our case, the inhomogeneous components of (18) and (19) are comprised of spatial gradients in the peculiar volume expansion/contraction ($\tilde{\vartheta}$), in the matter density ($\tilde{\Delta}_a$) and in the universal expansion ($\tilde{Z}_a$) – of their first and second time-derivatives in particular. The effect of 3-gradients weakens and becomes less prominent as one moves to progressively larger scales. On these grounds, spatial gradients are typically dropped on scales close and beyond the Hubble horizon. This may reduce the range of the solutions, but provides important information regarding their large-scale behaviour and until horizon-crossing at least. Therefore, hereafter, our study will focus on the long-wavelength solutions.<fnref n="9" />

Isolating the homogeneous component of Eq. (18), while setting $\Omega = 1$ and $H = 2/3t$ at the same time, we have

$$ 9t^2 \tilde{v}''_a + 3t \tilde{v}'_a - 8 \tilde{v}_a = 0, \quad (20) $$

on all scales where its inhomogeneous part is subdominant. The above accepts the power-law solution

$$ \tilde{v} = C_1 t^{4/3} + C_2 t^{-2/3} = C_3 a^2 + C_4 a^{-1} \quad (21) $$

which shows increase as $\tilde{v} \propto t^{4/3} \propto a^2$ (since $a \propto t^{2/3}$ after equipartition) for linear peculiar-velocity perturbations.<fnref n="10" /> Keeping only the growing mode, we deduce that $\tilde{v}/v_H \propto t^{5/3} \propto a^{5/2}$ after decoupling, with $v_H \propto a^{-1/2}$ giving the evolution of the Hubble velocity during the same period. These growth rates are significantly stronger than those reported in previous studies. Indeed, the Newtonian and quasi-Newtonian approaches give $v \propto t^{1/3} \propto a^{1/2}$ (e.g. see [2–6] as well as [7]), which implies $v/v_H \propto t^{2/3} \propto a$ (on an $\Omega = 1$ background). On these grounds, the relativistic analysis presented here seems to favour bulk motions larger and faster than it is generally expected. An alternative interpre-

$$ 3\tilde{v} = \left(\tilde{v}_0 + \frac{\tilde{v}'_0}{H_0}\right) \left(\frac{1+z_0}{1+z}\right)^2 + \left(2\tilde{v}_0 - \frac{\tilde{v}'_0}{H_0}\right) \left(\frac{1+z}{1+z_0}\right), \quad (22) $$

<footnote n="9">⁹ As required, we applied the long-wavelength approximation to the final formulae and not at an earlier stage. This ensures that all the linear effects have been consistently incorporated into Eqs. (18) and (19).</footnote>

<footnote n="10">¹⁰ Evaluating the integration constants and using the cosmological redshift parameter ($z$) solution (21) reads</footnote>

<footnote>where the zero suffix marks the initial time. The above provides the peculiar velocity at redshift $z < z_0$, when one knows the initial conditions. These are typically set at decoupling, namely at $z_0 = 10^3$. According to (22), essentially all the contribution to the residual peculiar velocity comes from higher redshifts, with $1 < z < z_0$. Therefore, allowing for a late-time accelerated epoch (typically starting at $z < 1$), should for all practical purposes leave our results unaffected.</footnote>
tation of our results is that the observed peculiar velocities could have started considerably weaker than anticipated.

The analytic results given above apply to scales where the inhomogeneous component of (18) is sub-dominant. This means that, although one could readily apply solution (21) beyond and near the Hubble radius, they should be cautious before doing so on scales well inside the horizon. There, the relativistic solution could reduce to its Newtonian counterpart, in which case it could also provide the initial conditions (near horizon-crossing) for a subsequent Newtonian treatment. Solving (18) on sub-Hubble lengths will probably require numerical treatment, with the analytical work providing the initial conditions. That aside, stronger growth-rates for the peculiar-velocity field on super-Hubble lengths, imply faster peculiar velocities at horizon-crossing. This in turn suggests higher residual values today and therefore provides theoretical support to a number of recent surveys reporting bulk peculiar flows larger and faster (sometimes considerably) than those typically expected [10–16].

## 4.3 Peculiar expansion/contraction

Let us go back to Eq. (19), which governs the linear evolution of the peculiar-velocity gradients. Assuming again a spatially flat FRW background, the trace of (19) gives

$$ \tilde{\vartheta}'' = -\frac{5H}{2}\tilde{\vartheta}' + 2H^2\tilde{\vartheta} + \frac{1}{3H}\tilde{D}^2\tilde{\vartheta}' - \frac{1}{3a^2H}(\tilde{\Delta}'' + \tilde{Z}') \quad (23) $$

where $\tilde{\Delta} = \tilde{\Delta}_a^a$ and $\tilde{Z} = \tilde{Z}_a^a$ by default. Applying our long-wavelength approximation, we retain only the first two terms on the right-hand side of the above. Then, after equipartition, Eq. (23) reads

$$ 9t^2\tilde{\vartheta}'' + 15t\tilde{\vartheta}' - 8\tilde{\vartheta} = 0 \quad (24) $$

with a power-law solution of the form

$$ \tilde{\vartheta} = C_1 t^{2/3} + C_2 t^{-4/3} = C_3 a + C_4 a^{-2} \quad (25) $$

Therefore, throughout the dust epoch, $\tilde{\vartheta}$ grows in tune with the dimensions of the universe, as long as the velocity perturbation remains outside the Hubble scale.

Including the Laplacian term does not prevent Eq. (23) from accepting an analytic solution. Although, this may seem at odds with our adopted approximation scheme, it will allow us to probe (at least to a certain extent) the evolution of $\tilde{\vartheta}$ on sub-horizon scales as well. Then, a simple Fourier decomposition, leads to the differential equation

$$ \tilde{\vartheta}''(n) = -\frac{5}{3t} \left[ 1 + \frac{2}{15} \left( \frac{\lambda_H}{\lambda_n} \right)^2 \right] \tilde{\vartheta}'(n) + \frac{8}{9t^2} \tilde{\vartheta}(n) \quad (26) $$

for the n-th harmonic mode.<fnref n="11" /> Here, $\lambda_H = 1/H$ is the Hubble radius, $\lambda_n = a/n$ is the physical scale of the perturbation (i.e. of the bulk flow) and $n$ is the associated (moving) wavenumber. The above has a scale-dependent solution. Indeed, introducing the scale parameter $\alpha = \lambda_H/\lambda_n$ – with $\alpha > 0$, we may express the solution of Eq. (26) in terms of the $\alpha$-parameter as

$$ \tilde{\vartheta} = C_1 t^{\beta_1} + C_2 t^{\beta_2} = C_3 a^{3\beta_1/2} + C_4 a^{3\beta_2/2} \quad (27) $$

with

$$ \beta_{1,2} = -\frac{\alpha^2 + 3 \mp \sqrt{a^4 + 6\alpha^2 + 81}}{9} \quad (28) $$

Well inside the horizon, where $\alpha \ll 1$, we recover solution (25) with $\beta_1 \simeq 2/3$ and $\beta_2 \simeq -4/3$. At horizon crossing, namely at the $\alpha = 1$ threshold, one obtains $\beta_1 = 2(\sqrt{22} - 2)/9$ and $\beta_2 = -2(\sqrt{22} + 2)/9$, which shows a slight decrease in the growth-rate of $\tilde{\vartheta}$, relative to the super-Hubble solution. The slowing-down effect continuous as we move to progressively smaller lengths. In fact, on scales deep inside the Hubble radius (where $1/\alpha \to 0$), a simple (linear) Taylor expansion of (28) gives $\beta_{1,2} \simeq -[\alpha^2 + 3 \mp \alpha^2(1 + 3/\alpha^2)]/9$, with $\beta_1 \simeq 0$ and $\beta_2 \simeq -2\alpha^2/9$. In other words, on sufficiently small scales the value $\tilde{\vartheta}$ tends to a constant.

The above results suggest that, after matter-radiation equality, the relative-strength ratio $\tilde{\vartheta}/H$ grows as $\tilde{\vartheta}/H \propto t^{5/3} \propto a^{5/2}$ on super-Hubble lengths. Well inside the horizon, on the other hand, we find that $\tilde{\vartheta}/H \propto t \propto a^{3/2}$. The divergence of the peculiar velocity field is typically related to the density contrast of the matter inhomogeneities (e.g. see [2–4]). On these grounds, peculiar velocity perturbations that cross inside the horizon at decoupling, with $\tilde{\vartheta}/H \simeq 10^{-5}$ at that time, could reach values as high as $\tilde{\vartheta}/H \sim 10^{-1/2}$ today. Finally, we should note that our growth rates are again significantly stronger than the Newtonian ones, according to which $\tilde{\vartheta}/H \propto t^{2/3} \propto a$ (e.g. see [2–4] as well as [7]).

## 4.4 Peculiar shear and vorticity

Taking the symmetric and traceless part of (19), we obtain the evolution formula of the peculiar shear. The skew component of the same relation, on the other hand, determines the propagation of the peculiar vorticity. More specifically, on an Einstein-de Sitter background, we have

$$ \tilde{\zeta}''_{ab} = -\frac{5}{2} H \tilde{\zeta}'_{ab} + 2H^2 \tilde{\zeta}_{ab} + \frac{1}{3H} \tilde{D}_{\langle b} \tilde{D}_{a\rangle} \tilde{\vartheta}' - \frac{1}{3a^2H} (\tilde{\Delta}''_{\langle ab\rangle} + \tilde{Z}'_{\langle ab\rangle}) \quad (29) $$

<footnote n="11">¹¹ We use the familiar harmonic splitting $\tilde{\vartheta} = \sum_n \tilde{\vartheta}_{(n)} Q^{(n)}$, where $\tilde{D}_a \tilde{\vartheta}_{(n)} = 0$ and $Q^{(n)}$ are scalar harmonic functions with $\dot{Q}^{(n)} = 0$ and $\tilde{D}^2 Q^{(n)} = -(n/a)^2 Q^{(n)}$ (e.g. see [8,9]).</footnote>
and

$$ \tilde{\omega}_{ab}'' = -\frac{5}{2} H \tilde{\omega}_{ab}' + 2H^2 \tilde{\omega}_{ab} - \frac{1}{3a^2 H} (\tilde{\Delta}_{[ab]}'' + \tilde{Z}_{[ab]}') \quad (30) $$

respectively.<fnref n="12" /> When pressureless (baryonic or/and CDM) matter dominates, we have $a \propto t^{2/3}$ and $H = 2/3t$. In such a case, the homogeneous component of Eq. (29) recasts as

$$ 9t^2 \tilde{\zeta}_{ab}'' + 15t \tilde{\zeta}_{ab}' - 8 \tilde{\zeta}_{ab} = 0. \quad (31) $$

The above accepts the power-law solution

$$ \tilde{\zeta} = C_1 t^{2/3} + C_2 t^{-4/3} = C_3 a + C_4 a^{-2}, \quad (32) $$

which holds on sufficiently large scales where the inhomogeneous part of (29) is subdominant (see also Sect. 4.2 earlier). On these long wavelengths, the peculiar shear grows proportionally to the dimensions of the post-recombination host universe. The same is also true for the peculiar vorticity, since (after equipartition) the homogeneous component of (30) accepts a power-law solution identical to (31). Therefore, on scales where the higher-order derivatives in the last term of (30) are negligible, the peculiar vorticity also grows as $\tilde{\omega} \propto t^{2/3} \propto a$.<fnref n="13" />

These results imply that, after decoupling, the relative strength of the peculiar shear and that of the peculiar vorticity grows as $\tilde{\zeta}/H$, $\tilde{\omega}/H \propto t^{5/3} \propto a^{5/2}$. Therefore, on large enough scales (near and outside the Hubble horizon), peculiar flows could start with negligibly small amounts of shear and vorticity and still reach cosmologically relevant magnitudes today.

## 5 Discussion

Studies of cosmological peculiar motions have a research history that goes back several decades. Nevertheless, the great majority of the available theoretical work is essentially Newtonian in nature, despite the fact that the scales involved are a good fraction of the Hubble horizon. Moreover, all the aforementioned studies are performed in the CMB frame, although

no real observer in the universe follows the smooth Hubble flow and despite the subtlety of the relative-motion effects. With these in mind, we have attempted a relativistic treatment of peculiar velocities in a “tilted” almost-FRW universe (with pressureless baryons or/and CDM) and conducted our analysis in the rest-frame of a typical galaxy (like our Milky Way), that moves relative to the universal expansion. Our main aim was to provide a better theoretical understanding of the bulk-motion kinematics. Further motivation came from an apparent disagreement between the current bulk-flow surveys, with some reporting larger magnitudes and scales for the peculiar velocity fields than others (see [10–16] and [30–34], respectively, for representative though incomplete lists).

Even at the linear level, the kinematics and the dynamics of the universe appear different in the coordinate system of the bulk peculiar flow than in the frame of the Hubble expansion, simply because of relative motion effects. For instance, although the cosmic medium may look like a perfect fluid in the CMB frame, it will appear imperfect to the real observers solely due to their motion with respect to the smooth universal expansion [5,6]. Taking these fairly well known relativistic effects into account and employing linear (relativistic) cosmological perturbation theory, enabled us to obtain analytical expressions for the linear sources of peculiar velocities. This allowed us to go a step further than the quasi-Newtonian approach, where an effective gravitational potential and an evolution ansatz were introduced to address the issue.

The linear analysis confirmed that peculiar motions are the result of the ongoing structure formation process, and more specifically of the increasing inhomogeneity of the post-recombination universe. Our study also provided the first (to the best of our knowledge) relativistic insight to the evolution of the full peculiar kinematics. Technically speaking, this was achieved through a set of four differential equations, monitoring the linear propagation of the peculiar velocity itself, as well as those of the associated irreducible (local) kinematic quantities. The latter are the expansion/contraction, the shear and the rotation of the bulk flow. Solving the aforementioned differential formulae analytically, we found substantial growth for all aspects of the peculiar motion on sufficiently large scales. Moreover, the strength of the peculiar-velocity field, relative to the background Hubble expansion of the post-recombination universe, was found to increase in time as well. In particular, peculiar velocities were found to grow as $\tilde{v} \propto a^2$, a rate considerably stronger than the one reported in previous Newtonian studies (where $v \propto a^{1/2}$). Stronger than the Newtonian rates were also obtained for the peculiar expansion/contraction, as well as for the peculiar shear and for the peculiar vorticity.

On theoretical grounds, the disagreement between the relativistic and the Newtonian results is due to the different way the two theories treat issues as fundamental as the grav-

<footnote n="12">¹² The last term on the right-hand side of Eq. (30) stems from Frobenius’ theorem (e.g. see [27,28]), which ensures that rotating spacetimes do not possess integrable 3-dimensional hypersurfaces. Alternatively, one could say that the 4-velocity field is not hypersurface orthogonal in the presence of rotation. As a result, the (covariant) spatial gradients of scalars do not commute in rotating spaces, which in turn guarantees that $\tilde{\Delta}_{[ab]}$, $\tilde{Z}_{[ab]} \neq 0$. The interested reader is referred to [29] for further discussion and for an application of the Frobenius theorem to the relativistic study of cosmological perturbations.</footnote>

<footnote n="13">¹³ To the best of our knowledge, the only Newtonian treatment of the peculiar shear and the peculiar vorticity is the one recently given in [7]. There, it was found that $\tilde{\zeta}/H \propto t^{2/3} \propto a$ and also that $\tilde{\omega}/H \propto t^{-1/3} \propto a^{-1/2}$ after equipartition.</footnote>
itational field itself. General relativity, in particular, advocates that, in addition to the energy density and the pressure (isotropic and/or anisotropic), the energy flux gravitates as well. Applied to peculiar motions, this principle ensures a flux-contribution to the energy-momentum tensor that is entirely due to the (peculiar) motion of the matter. This then feeds into the conservation laws and eventually leads to expression (16), which (together with Eqs. (18) and (19)) plays a central role in this study and it can be seen as the relativistic correction to the Newtonian analysis (see Sects. 2.4 and 3.2 for details). The quasi-Newtonian approach and other relativistic studies have also bypassed the aforementioned input of the bulk-flow flux to the local gravitational field and the typical reasons are discussed in Sects. 3.1 and 3.2 respectively.

We finally remind the reader that we obtained our power-law solutions by confining to the long-wavelength limit of the associated differential equations. This means that the aforementioned analytic results apply on large enough scales, with the typical threshold set by the Hubble radius. On smaller scales, the solutions reported here should be treated with caution, though they can still provide the initial conditions for future (analytical or numerical) studies. At this stage, our relativistic analysis does suggest that the linear growth of large-scale peculiar velocities can be considerably stronger than it is generally expected. Indeed, even if these stronger growth-rates are confined to super-Hubble lengths, the residual peculiar-velocity field should be larger than anticipated. On these grounds, one should not be surprised to measure bulk flows larger and faster than it is generally expected. Like those reported in [10–16] for example.

**Acknowledgements** ET was supported by a CNRS internship fellowship at the Institut d'Astrophysique de Paris (IAP), during the later stages of this work. CGT wishes to acknowledge support from a visiting fellowship by Clare Hall College and visitor support by DAMTP at the University of Cambridge, where part of this work was conducted. CGT also acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.), under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project Number: 789).

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited [Authors' comment: There is no data associated with our manuscript. It is a purely theoretical paper and no experimental data was used.].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-

right holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP³.

## References

1. P.A.R. Ade et al., Astron. Astrophys. 561, A97 (2014)

2. P.J.E. Peebles, *The Large-Scale Structure of the Universe* (Princeton University Press, Princeton, New Jersey, 1980)

3. A. Nusser, A. Dekel, E. Bertschinger, G.R. Blumenthal, Astrophys. J. 379, 6 (1991)

4. T. Padmanabhan, *Structure Formation in the Universe* (Cambridge University Press, Cambridge, 1993)

5. R. Maartens, Phys. Rev. D 5(8), 124006 (1998)

6. G.F.R. Ellis, H. van Elst, R. Maartens, Class. Quant. Grav. 18, 5115 (2001)

7. K. Filippou, C.G. Tsagas, arXiv:2003.01186

8. C.G. Tsagas, A. Challinor, R. Maartens, Phys. Rep. 465, 61 (2008)

9. G.F.R. Ellis, R. Maartens, M.A.H. MacCallum, *Relativistic Cosmology* (Cambridge University Press, Cambridge, 2012)

10. R. Watkins, H.A. Feldman, M.J. Hudson, Mon. Not. R. Astron. Soc. 392, 743 (2009)

11. A. Kashlinsky, F. Atrio-Barandela, D. Kocevski, H. Ebeling, Astrophys. J. 686, L49 (2009)

12. G. Lavaux, R.B. Tully, R. Mohayaee, S. Colombi, Astrophys. J. 709, 483 (2010)

13. J. Colin, R. Mohayaee, S. Sarkar, A. Shafieloo, Mon. Not. R. Astron. Soc. 414, 264 (2011)

14. C. Magoulas et al., IAU Symp. 308, 336 (2016)

15. A. Salehi, M. Yarahmadi, S. Fathi. arXiv:2001.01743

16. R. Mohayaee, M. Rameez, S. Sarkar. arXiv:2003.10420

17. A.R. King, G.F.R. Ellis, Commun. Math. Phys. 31, 209 (1973)

18. J. Ehlers, Akad. Wiss. Lit. Mainz Abh. Math.-Nat. Kl. 11, 792 (1961)

19. G.F.R. Ellis, C.G. Tsagas, Phys. Rev. D 66, 124015 (2002)

20. C.G. Tsagas, M.I. Kadiltzoglou, Phys. Rev. D 88, 083501 (2013)

21. J.M. Stewart, M. Walker, Proc. R. Soc. A 341, 49 (1974)

22. C.G. Tsagas, M.I. Kadiltzoglou, Phys. Rev. D 92, 043515 (2015)

23. G.F.R. Ellis, Mon. Not. R. Astron. Soc. 243, 509 (1990)

24. J.A.R. Cembranos, A.L. Maroto, H. Villarrubia-Rojo, JCAP 06, 041 (2019)

25. G. Rigopoulos, W. Valkenburg, Phys. Rev. D 86, 043523 (2012)

26. I. Milillo, D. Bertacca, M. Bruni, A. Maselli, Phys. Rev. D 92, 023519 (2015)

27. R.M. Wald, *General Relativity* (Chicago University Press, Chicago, 1984)

28. E. Poisson, *A Relativist's Toolkit* (Cambridge University Press, Cambridge, 2004)

29. G.F.R. Ellis, M. Bruni, J. Hwang, Phys. Rev. D 42, 1035 (1990)

30. A. Nusser, M. Davis, Astrophys. J. 736, 93 (2011)

31. S.J. Osborne, D.S.Y. Mak, S.E. Church, E. Pierpaoli, Astrophys. J. 737, 98 (2011)

32. E. Branchini, M. Davis, A. Nusser, Mon. Not. R. Astron. Soc. 424, 472 (2012)

33. S.J. Turnbull et al., Mon. Not. R. Astron. Soc. 420, 447 (2012)

34. Y.-Z. Ma, J. Pan, Mon. Not. R. Astron. Soc. 437, 1996 (2014)