Relativistic approach to the kinematics of large-scale peculiar motions

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 a2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^2$$\end{document}, in terms of the scale factor, instead of the Newtonian a1/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^{1/2}$$\end{document}-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.


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. a e-mail: tsagas@astro.auth.gr (corresponding author) 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][3][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 largescale 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-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 (ṽ) is found to grow asṽ ∝ 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ṽ/v H -with v H representing the Hubble velocity on the corresponding length, grows as v/v H ∝ a 5/2 . At the same time, the expansion/contarctionrate (θ) of the bulk flow shows a scale-dependent growth, which reaches its maximum (θ ∝ a) on super-Hubble scales and tends to a constant (i.e.θ → constant) inside the horizon. Finally, the peculiar shear and the peculiar vorticity (ς and˜ respectively) are found to grow asς,˜ ∝ 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 ∝ a 1/2 in the CMB frame [2][3][4][5][6], which implies v/v H ∝ 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 conservations laws and eventually appears in the equations monitoring cosmological perturbations and, more specifically, the evolu-tion 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][11][12][13][14][15][16].

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".

The 4-velocity "tilt"
Let us consider two families of observers with (timelike) 4velocities u a andũ a respectively. Let us also assume that the latter family has peculiar velocityṽ a relative to the former. The three aforementioned velocity fields are then related by the familiar Lorentz boost where u a u a = −1 =ũ aũ a , u aṽ a = 0,γ = 1/ √ 1 −ṽ 2 and v 2 =ṽ aṽ a < 1 (see Fig. 1). 1 The "tilt" between the two 4-velocity vectors is determined by the (hyperbolic) angle β, defined by cosh β = −ũ a u a =γ . The latter ensures that β = ln(γ + γ 2 − 1), with β > 0 sinceγ > 1 (see [17] for further discussion and additional technical details). Note that, when dealing with large-scale bulk peculiar motions, the local and the mean velocities of the flow (ṽ andṼ respectively) are typically related via the integral V = (3/4πr 3 ) x<rṽ dx 3 , whereṽ 2 =ṽ aṽ a and r is the (effective) radius of the bulk flow.

The 1 + 3 threading
Introducing two 4-velocities means that there are two temporal directions (along u a andũ a respectively) and two 3dimensional spatial sections (orthogonal to u a andũ a ). Projecting onto these 3-spaces is achieved by using the symmetric projection tensors h ab = g ab + u a u b andh ab = g ab +ũ aũb , with h ab u a = 0 =h abũ b and h a a = 3 =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 timederivatives in the u a and theũ a frames are denoted by respectively. The corresponding spatial gradients, on the other hand, are 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ũ a field), thus achieving an 1 + 3 threading of the host spacetime into time and space [8,9].

Kinematic decomposition
The irreducible kinematic variables of the two 4-velocity fields emerge by decomposing their gradients. More specifically, we have [18] In the above Θ = D a u a , σ ab = D b u a , ω ab = D [b u a] and A a =u a are respectively the expansion/contraction scalar, the shear tensor, the vorticity tensor and the 4-acceleration vector of the u a -field. 2 Positive values for Θ 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 for the gradient of theũ a -field. Note that, in cosmological studies, the volume scalar defines the scale factor of the universe by means ofȧ/a = Θ/3. Also, the shear, the vorticity and the 4-acceleration are all spacelike, namely σ ab u b = 0 = ω ab u b = A a u a by construction, with exactly analogous constraints applying to their "tilded" counterparts (i.e.σ abũ b = 0 =ω abũ b =Ã aũ 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ṽ 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] withθ =D aṽ a ,ς ab =D bṽa and˜ ab =D [bṽa] representing the peculiar expansion/contraction, the peculiar shear and the peculiar vorticity respectively. 3 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ṽ 2 1 and u a u a +ṽ a , sinceγ 1. 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. 5 Then, the kinematic variables defined in Sect. 2.3 are related by the linear expressions [5,6] Θ = Θ +θ,σ ab = σ ab +ς ab ,ω ab = ω ab +˜ ab (9) and whereṽ a =ũ b ∇ bṽ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.Ã a =ṽ a + (Θ/3)ṽ a when A a = 0). In a similar manner, one can show that the dynamical variables measured in the two frames are related bỹ ρ = ρ,p = p,q a = q a − (ρ + p)ṽ a (11) and to first approximation [5,6]. Here, ρ is the energy density, p is the isotropic pressure, q a is the energy flux and π ab is the viscosity of the matter as measured in the u a -frame, while 4 Typical Newtonian studies, define peculiar velocities by introducing physical and comoving coordinates (r α and x α respectively, with α = 1, 2, 3), related by r α = ax α . The time derivative of the latter leads to where v α =ṙ α is the physical (total) velocity, v α H = Hr α is the Hubble velocity and v α p = aẋ α 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ũ a and u a are both timelike vectors, whereas their Newtonian counterparts are all purely spatial. 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. their tilded counterparts are associated with theũ 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).

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.

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. 6 Then, relation (10), which also serves as the linear propagation equation of the peculiar velocity in the bulk-flow frame, recasts as v a = −Hṽ a +Ã a , with H = Θ/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Ã a the key to the subsequent evolution of theṽ 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][3][4]7]). Similarly, the quasi-Newtonian approach (which 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Ã a = 0, in which case A = 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.
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 tilded frame, we will obtain analytical expressions for both the 4-acceleration vector and its time-derivative.

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 tilded frame solely due to the latter's relative motion. More specifically, following (11c), we deduce thatq a = −ρṽ 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 tilded frame, the evolution formula of density inhomogeneities (e.g. see Eqs. (2.3.1) and (10.101) in [8,9]), reads [22] Here,Δ a = (a/ρ)D a ρ andZ a = aD aΘ 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 Keeping in mind thatq a = −ρṽ a , the latter combines with Eq. (14) to givẽ Note that in a Newtonian perturbative study Eq. (14) reduces toΔ α = −Z α , withΔ α = (a/ρ)∂ α ρ andZ α = a∂ αΘ [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. 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]). 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 tilded 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θ vanishes in the absence of drift motions, expression (16) leavesΔ a and Z a as the sole sources of linear peculiar velocities. Therefore, as expected, theṽ 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 bỹ Δ a andZ a respectively). 7 Ours is a single-fluid study. Nevertheless, one can still recover the results of the multi-fluid analysis "phenomenologically", by setting q a = 0 in Eq. (15) and thereforeÃ a = 0 in (13). The latter then leads toṽ a ∝ a −1 and toθ ∝ 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. 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.

Linear evolution of peculiar flows
Expression (13) requires an evolution formula for the 4acceleration to close the system of the differential equations.
In what follows, we will do so and also obtain analytic solutions for theṽ a -field.

The key differential equations
Taking the time derivative of Eq. (16), recalling thatḢ = −H 2 [1 + (Ω/2)] in the FRW background (with Ω = κρ/3H 2 being the associated density parameter) and using the linear commutation law (D aθ ) =D aθ − HD aθ between the temporal and the spatial derivatives of first-order scalars [5,6], provides the propagation formulã 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 timederivative of Eq. (13) leads tõ In addition, using the linear commutation lawsD bṽ a = (D bṽa ) + HD bṽa andD vṽ a = (D bṽa ) + 2H (D bṽa ) − (H 2 Ω/2)D bṽa , the linearised 3-gradient of the above reads whereΔ ab = aD bΔa andZ ab = aD bZa . 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.

Peculiar velocity
To this point, we have considered a perturbed Friedmann universe without imposing any constraints on its spatial curva-ture. Hereafter, we will confine to an Einstein-de Sitter background by setting Ω = 1 and H = 2/3t. However, Eqs. (18) and (19) do not accept analytic solutions even at the Ω = 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 (θ), in the matter density (Δ a ) and in the universal expansion (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 longwavelength solutions. 9 Isolating the homogeneous component of Eq. (18), while setting Ω = 1 and H = 2/3t at the same time, we have on all scales where its inhomogeneous part is subdominant. The above accepts the power-law solutioñ which shows increase asṽ ∝ t 4/3 ∝ a 2 (since a ∝ t 2/3 after equipartition) for linear peculiar-velocity perturbations. 10 Keeping only the growing mode, we deduce that v/v H ∝ t 5/3 ∝ a 5/2 after decoupling, with v H ∝ 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 ∝ t 1/3 ∝ a 1/2 (e.g. see [2][3][4][5][6] as well as [7]), which implies v/v H ∝ t 2/3 ∝ a (on an Ω = 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- 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). 10 Evaluating the integration constants and using the cosmological redshift parameter (z) solution (21) reads 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. 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][11][12][13][14][15][16].

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 whereΔ =Δ a a andZ =Z a a by default. Applying our longwavelength approximation, we retain only the first two terms on the right-hand side of the above. Then, after equipartition, Eq. (23) reads with a power-law solution of the form Therefore, throughout the dust epoch,θ 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θ on sub-horizon scales as well. Then, a simple Fourier decomposition, leads to the differential equatioñ for the n-th harmonic mode. 11 Here, λ H = 1/H is the Hubble radius, λ n = a/n is the physical scale of the perturbation (i.e. of the bulk flow) and n is the associated (comoving) wavenumber. The above has a scale-dependent solution. Indeed, introducing the scale parameter α = λ H /λ n -with α > 0, we may express the solution of Eq. (26) in terms of the α-parameter as with Well inside the horizon, where α 1, we recover solution (25) with β 1 2/3 and β 2 −4/3. At horizon crossing, namely at the α = 1 threshold, one obtains β 1 = 2( √ 22 − 2)/9 and β 2 = −2( √ 22 + 2)/9, which shows a slight decrease in the growth-rate ofθ, 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/α → 0), a simple (linear) Taylor expansion of (28) gives β 1,2 −[α 2 + 3 ∓ α 2 (1 + 3/α 2 )]/9, with β 1 0 and β 2 −2α 2 /9. In other words, on sufficiently small scales the valueθ tends to a constant.
The above results suggest that, after matter-radiation equality, the relative-strength ratioθ/H grows asθ/H ∝ t 5/3 ∝ a 5/2 on super-Hubble lengths. Well inside the horizon, on the other hand, we find thatθ/H ∝ t ∝ 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][3][4]). On these grounds, peculiar velocity perturbations that cross inside the horizon at decoupling, withθ/H 10 −5 at that time, could reach values as high asθ/H ∼ 10 −1/2 today. Finally, we should note that our growth rates are again significantly stronger than the Newtonian ones, according to which ϑ/H ∝ t 2/3 ∝ a (e.g. see [2][3][4] as well as [7]).

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 havẽ 11 We use the familiar harmonic splittingθ = nθ (n) Q (n) , wherẽ D aθ(n) = 0 and Q (n) are scalar harmonic functions withQ (n) = 0 and D 2 Q (n) = −(n/a) 2 Q (n) (e.g. see [8,9]). and ab = − respectively. 12 When pressureless (baryonic or/and CDM) matter dominates, we have a ∝ t 2/3 and H = 2/3t. In such a case, the homogeneous component of Eq. (29) recasts as 9t 2ς ab + 15tς ab − 8ς ab = 0.
The above accepts the power-law solutioñ 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˜ ∝ t 2/3 ∝ a. 13 These results imply that, after decoupling, the relative strength of the peculiar shear and that of the peculiar vorticity grows asς/H,˜ /H ∝ t 5/3 ∝ 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.

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 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Δ [ab] ,Z [ab] = 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. 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ς/H ∝ t 2/3 ∝ a and also that˜ /H ∝ t −1/3 ∝ a −1/2 after equipartition. 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][11][12][13][14][15][16] and [30][31][32][33][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 postrecombination 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 peculiarvelocity 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ṽ ∝ a 2 , a rate considerably stronger than the one reported in previous Newtonian studies (where v ∝ 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-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 powerlaw 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][11][12][13][14][15][16] for example.

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://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .