© The Authors 2026

1. Introduction

Large-scale peculiar motions, commonly known as bulk flows, are an observational certainty confirmed by numerous surveys over several decades. These surveys largely agree on the direction of the reported bulk flows, but they seem to disagree on their size and speed. Moreover, while many of the reported bulk flows have sizes and velocities within the theoretical boundaries of the current cosmological model (e.g. see Nusser & Davis 2011; Ma & Scott 2013; Scrimgeour et al. 2016; Qin et al. 2019), there has been an increasing number of recent surveys reporting bulk peculiar flows that are considerably faster and deeper than theoretically expected (e.g. see Watkins et al. 2009; Colin et al. 2011; Salehi et al. 2021; Watkins et al. 2023; Whitford et al. 2023). To these claims, one could also add the controversial ‘dark flows’, with sizes and velocities well in excess of the Lambda cold dark matter (ΛCDM) limits (Kashlinsky et al. 2008, 2011). Typically, the surveys that agree with the ΛCDM expectations cover scales smaller than 100 Mpc, whereas those that disagree extend beyond the aforementioned threshold. Currently, the reason for this discrepancy remains unknown and a subject of debate within the cosmological community since, if confirmed, the bulk-flow question could prove a serious problem for the concordance ΛCDM model. The reader is referred to Tsagas et al. (2026) for a recent review on cosmological peculiar motions, with more details and references.

When it comes to peculiar motions, the ΛCDM predictions are built on Newtonian theory, which argues for a relatively mild linear growth rate of v ∝ t^1/3 for the peculiar-velocity field (v) (Peebles 1976, 1980). This was also the result of the ‘quasi-Newtonian’ treatments (Maartens 1998; Ellis et al. 2001), though they reduce to mere Newtonian despite their initially relativistic profile. One of the reasons for remaining within the Newtonian limits was probably the anticipated technical difficulty of addressing the bulk-flow question relativistically. Another reason concerns the widespread expectation that when dealing with non-relativistic peculiar motions, Einstein’s theory was unlikely to drastically change the existing Newtonian picture. It is therefore not surprising that fully relativistic linear studies of peculiar velocities have only recently appeared in the literature (Tsaprazi & Tsagas 2020; Miliou & Tsagas 2024; Tsagas 2026). These first treatments argue for a considerably stronger peculiar growth rate, namely for v ∝ t^4/3, which has the potential of addressing both the higher speeds and the greater depths of the bulk flows reported in Watkins et al. (2009), Colin et al. (2011), Salehi et al. (2021), Watkins et al. (2023), Whitford et al. (2023), for example. The reason for the profound difference between the results of the Newtonian and quasi-Newtonian studies on the one hand and the relativistic studies on the other, can be tracked down to the fundamentally different way the two theories treat the gravitational field and its sources and more specifically to the entirely different role played by the ‘peculiar flux’ in these two theoretical approaches.

Peculiar flows are nothing but matter in motion, and moving matter means a non-zero energy flux (e.g. see Tsagas et al. 2008; Ellis et al. 2012). In general relativity, as opposed to Newtonian gravity, energy fluxes contribute to the stress-energy tensor and therefore gravitate (Tsaprazi & Tsagas 2020). The gravitational input of the peculiar flux, which survives at the linear perturbative level, is not accounted for in the Newtonian treatments nor in their quasi-Newtonian counterparts, though for different reasons. In the studies using Newtonian treatments, this happens by default, while in the quasi-Newtonian studies, it is due to the strict constraints the quasi-Newtonian approximation imposes on the perturbed spacetime. These compromise the relativistic nature of the approach and reduce it to mere Newtonian for all practical purposes. Although the problematic nature of the quasi-Newtonian approach has been known (e.g. see related ‘warning’ comments in Sect. 6.8.2 of Ellis et al. 2012 – pages: 150 and 151), this has not been widely appreciated. The likely reason is the absence of a direct comparison with a fully relativistic analysis (until very recently; see Tsagas 2026 for such a comparison).

In a fully relativistic linear study, the gravitational input of the peculiar flux feeds into the associated conservation laws and subsequently into the linear evolution formulae of the peculiar-velocity field. In so doing, the relativistic treatment brings into play the effects of density perturbations, namely of structure formation, in addition to those of the peculiar flux. All this drastically modifies the agent driving the linear peculiar velocities. As a result, unlike its Newtonian and quasi-Newtonian counterparts, the relativistic driving agent no longer decays with time. This in turn ensures that the relativistic differential equations and their associated solutions differ considerably. More specifically, the general relativistic analysis supports a substantially stronger linear growth for peculiar velocities during the Einstein-de Sitter phase of the universe. We demonstrate in this work that even in the ‘minimalist’ scenario where only the gravitational input of the peculiar flux is taken into account, the velocity growth rate is v ∝ t. The latter could potentially increase to the v ∝ t^4/3 rate of Tsaprazi & Tsagas (2020), Miliou & Tsagas (2024) when structure-formation effects are also accounted for – a possibility worthy of further investigation. In either case, the resulting velocity field has the potential to explain the speed and the depth of the fast bulk flows reported in Watkins et al. (2009), Colin et al. (2011), Salehi et al. (2021), Watkins et al. (2023), Whitford et al. (2023). Put another way, general relativity could naturally relax the typical (Newtonian) ΛCDM limits to accommodate the fast bulk peculiar flows.

This work provides a direct comparison between the Newtonian a,d quasi-Newtonian studies and the relativistic analysis of linear peculiar velocities on a Friedmann-Robertson-Walker (FRW) background with flat spatial sections. The linear nature of all the studies implies that their results formally apply to relatively large scales, typically near and beyond the 100 Mpc threshold, where one feels confident to ignore non-linear effects.

We begin our comparison by setting the theoretical framework of the approaches. We identify the point where they start to diverge from each other and explain why this happens. We then derive the linear evolution laws of the associated peculiar-velocity fields, highlight the differences between the corresponding solutions, and discuss their implications for addressing the ongoing bulk-flow question. In the process, we also highlight the role of the Newtonian acceleration and the relativistic four-acceleration as the driving forces of linear peculiar velocities. More specifically, we show that the stronger the four-acceleration, the faster the driven peculiar-velocity field. Finally, we demonstrate a way of recovering the relativistic results from the Newtonian setup by selectively accounting for certain source-free terms in Poisson’s formula that are typically sidestepped. Although, this recovery is rather ad hoc and not uniquely determined, it provides useful analogies between the two types of study. For instance, one could argue that ignoring the above mentioned source-free terms in the Newtonian studies is analogous to bypassing the purely general relativistic gravitational input of the peculiar flux.

We start this work with the traditional study of cosmological linear peculiar velocities, which employs Newtonian gravity and proceeds along the lines of Peebles (1976, 1980) and Nusser et al. (1991) in Sect. 2. We provide the covariant analogue of the Newtonian analysis in Sect. 3.3.1 to familiarise the reader with the covariant formalism and connect the standard Newtonian analysis of Sect. 2 with the relativistic approach of Sect. 3.3.3. In the latter section, we also discuss the physical reasons the Newtonian and quasi-Newtonian studies of linear peculiar velocities arrive at considerably weaker results compared to those of the relativistic treatment. In so doing, we highlight the key role played by the four-acceleration, keeping in mind that it is the driving force of the linear peculiar-velocity field (see Sect. 3.3.4). We also outline a way of bridging the aforementioned differences between the Newtonian and the relativistic studies by (selectively) including extra contributions to the Newtonian gravitational potential (see Sect. 4). Finally, we close with a discussion of our results and of their role in the quest for an answer to the bulk-flow puzzle.

2. Newtonian treatment of peculiar velocities

Here, we employ Newtonian theory and adopt co-moving coordinates {x}. These coordinates are related to their corresponding physical counterparts {r} by means of x = r/a, where a = a(t) is the cosmological scale factor. As a result, the velocity (u) of a typical observer in the Newtonian version of a perturbed Einstein-de Sitter universe decomposes as

$$\begin{array}{c}\hfill \mathbf{u}=H\mathbf{r}+\mathbf{v},\end{array}$$(1)

with H = ȧ/a representing the Hubble parameter. Also, $\mathbf{v}=a\dot{\mathbf{x}}$, which is the observer’s peculiar velocity relative to the universal rest frame, which has typically been identified with the rest-frame of the cosmic microwave background (CMB). Within the adopted cosmological framework, the linear evolution of the peculiar-velocity field is monitored by the associated Euler equation, namely by Peebles (1976, 1980)

$$\begin{array}{c}\hfill \dot{\mathbf{v}}+H\mathbf{v}=-\frac{1}{a}\mathrm{\nabla }\varphi ,\end{array}$$(2)

where ϕ is the peculiar gravitational potential and ∇ is the gradient in co-moving coordinates. The latter acts as a linear source of peculiar-velocity perturbations and satisfies Poisson’s formula, that is

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\varphi =4\pi G{a}^2{\rho }_0\delta ,\end{array}$$(3)

with δ = (ρ − ρ₀)/ρ₀ being the familiar density contrast, ρ the matter density, and ρ₀ its background value.

Differentiating Eq. (2) in time, recalling that H = 2/3t and $\dot{H}=2/3t^2$ to zero order, while keeping up to linear order terms leads to

$$\begin{array}{c}\hfill \ddot{\mathbf{v}}+2H\dot{\mathbf{v}}-\frac{1}{2}{H}^2\mathbf{v}=\ddot{\mathbf{v}}+\frac{4}{3t}\dot{\mathbf{v}}-\frac{2}{9t^2}\mathbf{v}=-\frac{1}{a}{\partial }_t\left(\mathrm{\nabla }\varphi \right).\end{array}$$(4)

One may ignore the right-hand side of the above equation by assuming that the gradient of the peculiar gravitational potential varies slowly in time, for example. Then, (4) reduces to

$$\begin{array}{c}\hfill 9t^2\frac{{\mathrm{d}}^2\mathbf{v}}{\mathrm{d}{t}^2}+12t\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}-2\mathbf{v}=0,\end{array}$$(5)

which accepts the power-law solution

$$\begin{array}{c}\hfill \mathbf{v}={\mathcal{C}}_1{t}^{1/3}+{\mathcal{C}}_2{t}^{-2/3},\end{array}$$(6)

on all scales. Therefore, according to Newtonian physics, linear peculiar velocities are expected to grow as v ∝ t^1/3 after recombination and throughout the Einstein-de Sitter epoch.

One can also refine solution (6) by incorporating the source term on the right-hand side of (4). Following Peebles (1976, 1980), we started by recalling that the peculiar gravitational acceleration is

$$\begin{array}{c}\hfill \mathbf{g}=-\frac{1}{a}\mathrm{\nabla }\varphi =Ga{\rho }_0{\displaystyle \int \delta }\frac{{\mathbf{x}}^{\prime }-\mathbf{x}}{||{\mathbf{x}}^{\prime }-\mathbf{x}{||}^{\mathbf{3}}}{\mathrm{d}}^3{\mathbf{x}}^{\prime }.\end{array}$$(7)

Given that ϕ and δ also depend on time, after a relatively straightforward calculation, the temporal derivative of the above leads to the linear time-evolution equation for ∇ϕ, namely

$$\begin{array}{c}\hfill {\left(\mathrm{\nabla }\varphi \right)}^{·}=-H\mathrm{\nabla }\varphi -\frac{3}{2}a{H}^2\mathbf{v}.\end{array}$$(8)

Note that in deriving this relation we have used the background expressions ρ₀ = 3H²/8πG and ${\dot{\rho }}_0=-3H{\rho }_0$ together with the linear result $\dot{\delta }=-\mathrm{\nabla }·\mathbf{v}$. The latter follows from the continuity equation, which reads $\dot{\rho }=-\mathcal{H}\rho$ at the linear level, with ρ and ℋ being the perturbed matter density and Hubble parameter, respectively. Setting ρ = ρ₀(1 + δ) and ℋ = H + ∇v, substituting into the linear continuity equation, and keeping up to first-order terms, we obtained the desired relation, namely, at $\dot{\delta }=-\mathrm{\nabla }·\mathbf{v}$. Finally, by combining Eqs. (2), (4) and (8), we arrived at the differential equation

$$\begin{array}{c}\hfill \ddot{\mathbf{v}}+3H\dot{\mathbf{v}}-H^2\mathbf{v}=\ddot{\mathbf{v}}+\frac{2}{t}\dot{\mathbf{v}}-\frac{4}{9t^2}\mathbf{v}=0\end{array}$$(9)

and the associated power-law solution Peebles (1976, 1980)

$$\begin{array}{c}\hfill \mathbf{v}={\mathcal{C}}_1{t}^{1/3}+{\mathcal{C}}_2{t}^{-4/3},\end{array}$$(10)

on all scales. Comparing solutions (6) and (10) made it clear that the inclusion of the non-zero right-hand side of (4) has no physical and/or practical significance. The rate of the growing mode in the former of the two solutions has remained the same and only the decaying mode has changed.

Our last note is on the peculiar acceleration g = −∇ϕ/a (see Eq. (7)), which is the driving force of the linear v field in Newtonian theory. By combining (2) with solution (10), one immediately realises that the Newtonian force decays with time as

$$\begin{array}{c}\hfill \mathbf{g}\propto t^{-2/3}\propto a^{-1},\end{array}$$(11)

since a ∝ t^2/3, after matter–radiation equality. In Sect. 3.3.4 we demonstrate that the same time-decay law also applies to the four-acceleration that drives linear peculiar velocities in the quasi-Newtonian studies. In contrast, the four-acceleration driving the peculiar-velocity field in the relativistic treatments is either constant (at a minimum) or increases with time (see Sect. 3.3.4). This difference offers another way of explaining why the Newtonian/quasi-Newtonian studies lead to the mediocre v ∝ t^1/3 growth and why general relativity supports faster and deeper bulk flows than generally anticipated.

3. Covariant treatment of peculiar velocities

The reported large-scale peculiar motions are believed to have started as weak peculiar-velocity perturbations that were amplified to the observed bulk flows by structure formation. In this section we provide a ‘unified’ covariant study looking into the linear evolution of peculiar velocities in cosmology. In so doing, we also facilitate direct comparison between the Newtonian, the quasi-Newtonian, and the relativistic treatments of the subject.

3.1. CMB and bulk-flow frames

Studies of peculiar velocities require two relatively moving coordinate systems, with one of them following the peculiar motion of the matter and the other the reference frame of the universe. The latter has been typically identified with the coordinate system of the CMB, namely the frame relative to which we define and measure peculiar motions.

In Newtonian theory, the velocities of the aforementioned observers are related via the familiar Galilean transformation

$$\begin{array}{c}\hfill {\stackrel{\sim }{u}}_{\alpha }=u_{\alpha }+v_{\alpha },\end{array}$$(12)

where ${\stackrel{\sim }{u}}_{\alpha }$ and u_α are the velocities of the bulk-flow observers and of their CMB partners, respectively. The term v_α is the peculiar velocity of the former with respect to the latter (with α = 1, 2, 3). In relativity, the Galilean transformation is replaced by the (also familiar) Lorentz boost, which in the case of non-relativistic peculiar motions reads

$$\begin{array}{c}\hfill {\stackrel{\sim }{u}}_a=u_a+v_a.\end{array}$$(13)

Here, in contrast to Eq. (12), ${\stackrel{\sim }{u}}_a$ and u_a are the (timelike) 4-velocities of the aforementioned observers, while v_a is the (spacelike) velocity of the peculiar motion (with a = 0, 1, 2, 3). Put another way, $u_a{u}^a=-1={\stackrel{\sim }{u}}_a{\stackrel{\sim }{u}}^a$ and u_av^a = 0 by construction. Also, v² = v_av^a ≪ 1 when dealing with non-relativistic peculiar flows. We note that although the transformations (12) and (13) may seem formally identical, they are different because the former involves spacelike vectors only, whereas in the latter the two 4-velocity vectors are timelike.

3.2. The key role of the peculiar flux

The key difference between the Newtonian and the relativistic study of peculiar motions follows from the fact that in relativity, the nature of the matter fields changes between the relatively moving frames and the associated observers. In particular, a cosmic medium that appears perfect to one reference frame will appear imperfect to any other relatively moving coordinate system. The imperfection takes the form of an effective energy flux vector, which is solely triggered by relative-motion effects. If we consider a perturbed Friedmann universe with zero spatial curvature, then when setting the pressure (both the isotropic and the viscous) to zero, we arrive at the linear relations

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }=\rho ,\stackrel{\sim }{p}=0,{\stackrel{\sim }{\pi }}_{\mathit{ab}}=0,\end{array}$$(14)

and

$$\begin{array}{c}\hfill {\stackrel{\sim }{q}}_a=q_a-\rho v_a,\end{array}$$(15)

between the CMB and the bulk-flow (i.e. the tilded) frames. In the above ρ, p, π_ab, and q_a are respectively the density, the isotropic pressure, the viscosity, and the energy flux of the matter relative to the CMB frame, with their tilded analogues measured in the matter frame. Following Eq. (14), the density and the pressure (both the isotropic and the viscous) of the matter are the same in the two coordinate systems at the linear level. On the other hand, according to Eq. (15), there is a difference between the energy fluxes entirely due to relative-motion effects. In fact, Eq. (15) ensures that even if the flux vanishes in one frame, there is a non-zero peculiar flux in all the other relatively moving coordinate systems, entirely due to relative-motion effects. Clearly, one cannot set the flux to zero in both frames simultaneously without switching the peculiar velocity off as well. Put another way, in the presence of peculiar motions, there is always a non-zero energy flux, which in turn means that the cosmic medium is no longer perfect.

When assuming a vanishing flux in the matter frame, we may set ${\stackrel{\sim }{q}}_a=0$. Then, Eq. (15) guarantees a non-zero peculiar flux equal to q_a = ρv_a in the CMB frame. The latter survives at the linear perturbative level, and in relativity, it contributes to the perturbed energy-momentum tensor and therefore to the relativistic gravitational field. As a result, the linear conservation laws of the energy and the momentum densities acquire flux-related terms and take the forms

$$\begin{array}{c}\hfill \dot{\rho }=-\mathrm{\Theta }\rho -{\mathrm{D}}^a{q}_a\mathrm{and}\rho A_a=-{\dot{q}}_a-4H{q}_a,\end{array}$$(16)

respectively. In the above equation, Θ = D^au_a is the three-divergence of the u_a field (with Θ = 3H in the background), $A_a={\dot{u}}_a$ is the associated four-acceleration, and D_a is the 3D covariant derivative operator. We note that according to Eq. (16b), the four-acceleration is non-zero, despite the absence of pressure. This is a direct consequence of the fact that in the presence of relative motions, there is always a non-zero energy flux in the system and the cosmic medium can no longer be treated as perfect¹. As it turns out (see Sect. 3.3 next), it is the peculiar-flux input to the gravitational field that separates the relativistic studies of peculiar motions from the rest.

3.3. Newtonian and quasi-Newtonian versus relativistic treatments

As mentioned above, the observed bulk peculiar flows began as peculiar-velocity perturbations, and they were subsequently amplified by structure formation. The driving agent of the amplification process depends on the theory that one employs to study the evolution of the peculiar-velocity field. In Newtonian physics, the driving force is gravity and, more specifically, the gravitational acceleration. Mathematically, this is reflected in the linear evolution equation of the v_α field, which reads

$$\begin{array}{c}\hfill {\dot{v}}_{\alpha }+H{v}_{\alpha }=-{\partial }_{\alpha }\varphi ,\end{array}$$(17)

with ϕ being the Newtonian gravitational potential, and the gradient is in physical coordinates². In relativity, gravity is not a force but the manifestation of spacetime curvature. As a result, in a relativistic treatment, linear peculiar velocities are driven by non-gravitational forces, and the evolution of linear peculiar velocities is monitored by

$$\begin{array}{c}\hfill {\dot{v}}_a+H{v}_a=-A_a,\end{array}$$(18)

where $A_a={\dot{u}}_a$ is the four-acceleration measured in the reference CMB frame³. It is therefore clear that the linear evolution of the v_a field is determined by the form and the impact of the associated driving agent.

3.3.1. Newtonian treatment

Newtonian peculiar-velocity studies start from the differential equation Eq. (17). Recalling that $\dot{H}=-3H^2/2$ to zero order, the linearised time-derivative of the former leads to

$$\begin{array}{c}\hfill {\ddot{v}}_{\alpha }+2H{\dot{v}}_{\alpha }-\frac{1}{2}{H}^2{v}_{\alpha }=-{\partial }_{\alpha }\dot{\varphi },\end{array}$$(19)

since ${\partial }_t({\partial }_{\alpha }\varphi )={\partial }_{\alpha }\dot{\varphi }-H{\partial }_{\alpha }\varphi$ to first approximation. As noted in Sect. 2, ignoring the right-hand side of the above equation leads to the power-law solution v = 𝒞₁t^1/3 + 𝒞₂t^−2/3. However, by appealing to the Jeans’ swindle (e.g. see Binney & Tremaine 1988) and writing Poisson’s equation as ∂²ϕ = ρδ/2⁴, a relatively lengthy but fairly straightforward calculation yields

$$\begin{array}{c}\hfill {\partial }_{\alpha }\dot{\varphi }=\frac{1}{2}H{\partial }_{\alpha }\varphi +\frac{3}{2}H{\dot{v}}_{\alpha },\end{array}$$(20)

which when substituted into the right-hand side of Eq. (19) gives

$$\begin{array}{c}\hfill {\ddot{v}}_{\alpha }+3H{\dot{v}}_{\alpha }-H^2{v}_{\alpha }=0,\end{array}$$(21)

with H = 2/3t in the background. The above equation can be solved to give

$$\begin{array}{c}\hfill v={\mathcal{C}}_1{t}^{1/3}+{\mathcal{C}}_2{t}^{-4/3}={\mathcal{C}}_3{a}^{1/2}+{\mathcal{C}}_4{a}^{-2}\end{array}$$(22)

on all scales, and it is in complete agreement with solution (10) in Sect. 2 and with Peebles (1976, 1980). Therefore, as expected, incorporating the right-hand side of Eq. (19) into the calculation only changes the decaying mode of the solution.

3.3.2. Quasi-Newtonian treatment

Expression (18) is also the starting point of the quasi-Newtonian studies of linear peculiar velocities in cosmology (Maartens 1998; Ellis et al. 2001). However, the agreement with the relativistic analysis stops at this point because of the strict constraints imposed on the perturbed spacetime by the quasi-Newtonian approximation. These constraints compromise the relativistic nature of the analysis and eventually lead to Newtonian equations and results (see Sect. 6.8.2 in Ellis et al. 2012 – pages 150 and 151 – for related ‘warning’ comments). More specifically, the quasi-Newtonian approach starts with setting the vorticity to zero at the linear level. This has the benefit of simplifying the mathematics because it allows one to express the four-acceleration as the spatial gradient of a scalar potential (φ) so that

$$\begin{array}{c}\hfill A_a={\mathrm{D}}_a\phi .\end{array}$$(23)

There is a downside, however, because φ is an ad hoc potential with no expression to describe it, while its linear evolution follows from the (non-uniquely determined) ansatz $\dot{\phi }=-\mathrm{\Theta }/3$ (Maartens 1998). In addition, the aforementioned ansatz requires setting the linear shear to zero as well. There is more to it, however. When it comes to the study of peculiar velocities, the most serious downside of the quasi-Newtonian treatment is that it inadvertently bypasses the gravitational contribution of the peculiar flux. This severely compromises the relativistic nature of the study and makes it identical to its purely Newtonian counterpart for all practical purposes. Indeed, combining Eqs. (18) and (23) the differential equation governing the linear evolution of the peculiar-velocity field takes the form (Maartens 1998; Ellis et al. 2001)

$$\begin{array}{c}\hfill {\ddot{v}}_a+3H{\dot{v}}_a-H^2{v}_a=0,\end{array}$$(24)

and accepts the power-law solution

$$\begin{array}{c}\hfill v={\mathcal{C}}_1{t}^{1/3}+{\mathcal{C}}_2{t}^{-4/3}={\mathcal{C}}_3{a}^{1/2}+{\mathcal{C}}_4{a}^{-2},\end{array}$$(25)

both of which are identical to their purely Newtonian analogues (compare to Eqs. (21) and (22) in Sect. 3.3.1). Therefore, according to the Newtonian and the quasi-Newtonian studies, linear peculiar velocities grow no faster than v ∝ t^1/3 ∝ a^1/2 during the Einstein-de Sitter era. However, this growth rate is not strong enough to explain the fast and deep bulk flows reported by several surveys over roughly the past fifteen years.

3.3.3. Relativistic treatment

The relativistic analysis of linear peculiar velocities imposes no constraints on the perturbed spacetime, and crucially, it also accounts for the gravitational input of the peculiar flux (Tsaprazi & Tsagas 2020; Miliou & Tsagas 2024; Tsagas 2026). The latter feeds into Einstein’s equations and emerges into the equations monitoring a cosmological model endowed with a peculiar-velocity field. As a result, there is no longer any need to introduce an ad hoc scalar potential (φ) for the four-acceleration (see Eq. (23)) or adopt the ansatz $\dot{\phi }=-\mathrm{\Theta }/3$ for its time evolution. Instead, the four-acceleration emerges naturally in the linear formulae governing the evolution of the density and the expansion gradients. These are respectively given by (Tsagas 2026)

$$\begin{array}{c}\hfill {\dot{\mathrm{\Delta }}}_a=-{\mathcal{Z}}_a-3aH{A}_a-a{\mathrm{D}}_a\vartheta \end{array}$$(26)

and

$$\begin{array}{c}\hfill \dot{{\mathcal{Z}}_a}=-2H{\mathcal{Z}}_a-\frac{3}{2}{H}^2{\mathrm{\Delta }}_a-\frac{9}{2}a{H}^2{A}_a+a{\mathrm{D}}_a{\mathrm{D}}^b{A}_b,\end{array}$$(27)

with Δ_a = (a/ρ)D_aρ and 𝒵_a = aD_aΘ representing perturbations in the matter density and the universal expansion, respectively⁵. Also, ϑ = D^av_a is the spatial divergence of the peculiar velocity. Expressions (26) and (27) provide the relativistic set of differential equations that govern the linear evolution of peculiar velocities in a perturbed Einstein-de Sitter universe. It is therefore clear that the relativistic analysis naturally couples the driving force of the peculiar-velocity field, namely the four-acceleration, to both the density and the expansion gradients.

It is important to note at this point that, had the gradient of the energy-conservation law (16a) been involved in Maartens (1998), that study would also have arrived at the relativistic expression (26) for the 4-acceleration, instead of the Newtonian-like Eq. (23) – see also footnote 5 here. This omission reduced the quasi-Newtonian analysis to purely Newtonian for all practical purposes.

By differentiating Eq. (26) with respect to time, substituting the right-hand side of Eq. (27) in the resulting expression, and then using Eq. (18), we arrived at (Tsagas 2026)

$$\begin{array}{c}\hfill {\ddot{v}}_a+H{\dot{v}}_a-\frac{3}{2}{H}^2{v}_a=\frac{1}{3aH}({\ddot{\mathrm{\Delta }}}_a+2H{\dot{\mathrm{\Delta }}}_a-\frac{3}{2}{H}^2{\mathrm{\Delta }}_a),\end{array}$$(28)

which is the relativistic formula governing the evolution of linear peculiar velocities in an Einstein-de Sitter cosmology. We note that in deriving the above equation, we also used the background relation $\dot{H}=-3H^2/2$ together with the linear commutation laws ${({\mathrm{D}}^a{v}_a)}^{·}={\mathrm{D}}^a{\dot{v}}_a-H{\mathrm{D}}^a{v}_a$ and ${({\mathrm{D}}_a\vartheta )}^{·}={\mathrm{D}}_a\dot{\vartheta }-H{\mathrm{D}}_a\vartheta$.

The expression (28) is a non-homogeneous differential equation, the homogeneous component of which reads⁶

$$\begin{array}{c}\hfill {\ddot{v}}_{\alpha }+H{\dot{v}}_{\alpha }-\frac{3}{2}{H}^2{v}_{\alpha }=0\end{array}$$(29)

and accepts the solution

$$\begin{array}{c}\hfill v={\mathcal{C}}_{1}t+{\mathcal{C}}_2{t}^{-2/3}={\mathcal{C}}_3{a}^{3/2}+{\mathcal{C}}_4{a}^{-1}\end{array}$$(30)

on all scales Tsagas (2026). Therefore the relativistic analysis has led to the considerably stronger growth rate of v ∝ t for linear peculiar velocities. Moreover, according to the theory of differential equations, solution (30) provides the minimum growth rate of the v field. Indeed, a fairly well known theorem states that the full solution of a non-homogeneous differential equation forms from the general solution of its homogeneous part–in our case from Eq. (30) – and from a partial solution of the full equation. Therefore, solving (28) in full will make a physical difference only if its partial solution grows faster than the fastest growing mode of the homogeneous solution. Put another way, the power-law v ∝ t is the minimum growth rate of the linear peculiar-velocity field. Indeed, by following an approach that accounts for some (though not all) of the effects of the density gradients seen on the right-hand side of Eq. (28), the studies of Tsaprazi & Tsagas (2020) and Miliou & Tsagas (2024) led to the growth rate of v ∝ t^4/3 for the linear v field. Overall, the relativistic analysis leads to considerably stronger growth rates for peculiar velocities than the Newtonian and quasi-Newtonian treatments. This means that general relativity supports faster and deeper residual bulk flows than is generally anticipated, perhaps the bulk flows such as those reported in the surveys of Watkins et al. (2009), Colin et al. (2011), Salehi et al. (2021), Watkins et al. (2023), Whitford et al. (2023).

3.3.4. The role of the four-acceleration

A complementary, less technical, but perhaps more physically intuitive way of demonstrating and explaining the faster growth rates of the relativistic analysis is by looking closer at the role of the four-acceleration. According to linear theory, the latter is the driving force of the v field (see Eq. (18) in Sect. 3.3). Therefore, the stronger the four-acceleration, the faster the peculiar velocity.

Following Eqs. (26) and (27), the input of the peculiar flux to the relativistic gravitational field also couples the peculiar four-acceleration to the density and the expansion gradients. If we ignore for the moment the effects of the aforementioned gradients and keep only that of the peculiar flux, then on using Eq. (18), expressions Eq. (26) and (27) reduce to

$$\begin{array}{c}\hfill A_a=-\frac{1}{3H}{\mathrm{D}}_a\vartheta \mathrm{and}9H^2{A}_a=-2{\mathrm{D}}_a\dot{\vartheta }-4H{\mathrm{D}}_a\vartheta ,\end{array}$$(31)

both of which combine to the linear relation

$$\begin{array}{c}\hfill {\mathrm{D}}_a\dot{\vartheta }=-\frac{1}{2}H{\mathrm{D}}_a\vartheta .\end{array}$$(32)

By taking the time derivative of Eq. (31a), employing Eq. (32), and using the auxiliary background and linear relations given in Sect. 3.3.3 previously, one can easily show that Ȧ_a = 0. Accordingly, when only the flux-effects are accounted for, linear peculiar velocities are driven by a time-invariant four-acceleration. In such cases, the differential Equation (18) accepts the linearly growing solution

$$\begin{array}{c}\hfill v={\mathcal{C}}_{1}t+{\mathcal{C}}_2{t}^{-2/3}\end{array}$$(33)

on all scales and in full agreement with solution (30) obtained in Sect. 3.3.3.

Keeping the density and the expansion gradients on the right-hand side of Eq. (26) and taking its time derivative, we arrived at the following propagation formula for the linear peculiar four-acceleration (Tsaprazi & Tsagas 2020):

$$\begin{array}{c}\hfill {\dot{A}}_a-\frac{1}{2}H{A}_a=-\frac{1}{3H}{\mathrm{D}}_a\dot{\vartheta }-\frac{1}{3aH}({\ddot{\mathrm{\Delta }}}_a+{\dot{\mathcal{Z}}}_a).\end{array}$$(34)

Solving the homogeneous left-hand side of the above differential equation gives A_a ∝ t^1/3 (we recall that H = 2/3t in the background). Finally, after substituting A_a ∝ t^1/3 into the right-hand side of Eq. (18), we found that

$$\begin{array}{c}\hfill v={\mathcal{C}}_3{t}^{4/3}+{\mathcal{C}}_2{t}^{-2/3}\end{array}$$(35)

in agreement with Tsaprazi & Tsagas (2020) (see also Maglara & Tsagas 2022; Miliou & Tsagas 2024). Consequently, when some of the gradient effects are also accounted for, the driving force of the peculiar-velocity could grow in time (as A_a ∝ t^1/3). As a result, the associated v field could grow faster than linearly, namely as v ∝ t^4/3 – a possibility worthy of further investigation.

Contrary to its relativistic counterparts, the quasi-Newtonian four-acceleration decays in time as A_a ∝ t^−2/3 (Maartens 1998). Substituting this result into Eq. (18), the latter leads to a linear peculiar-velocity field that grows as v ∝ t^1/3 only. It goes without saying that the purely Newtonian studies lead to the same exact results (e.g. see Eq. (11) in Sect. 2). The fact that in the Newtonian and quasi-Newtonian studies, the linear peculiar velocities are driven by a decaying force, whereas the relativistic driving agent is constant or even increasing in time, explains both the rather mediocre v ∝ t^1/3 growth rate of the former treatments and the substantially stronger (v ∝ t, or even v ∝ t^4/3) rates of the latter.

3.3.5. The peculiar vorticity and shear

A generic peculiar-velocity field is expected to have non-zero vorticity and shear, and they are respectively defined as the antisymmetric and the symmetric traceless components of the v field. In other words, ϖ_ab = D_[bv_a] is the peculiar vorticity, and ς_ab = D_⟨bv_a⟩ is the peculiar shear. These kinematic variables do not necessarily evolve as the peculiar velocity itself. Indeed, by taking the spatial gradient of Eq. (29) and using the linear commutation laws ${\mathrm{D}}_b{\dot{v}}_a={({\mathrm{D}}_b{v}_a)}^{·}+H{\mathrm{D}}_b{v}_a$ and ${\mathrm{D}}_b{\ddot{v}}_a={({\mathrm{D}}_b{v}_a)}^{·}+2H{\partial }_t({\mathrm{D}}_b{v}_a)-(H^2/2){\mathrm{D}}_b{v}_a$, we arrive at the following differential equation

$$\begin{array}{c}\hfill {\ddot{\varpi }}_{\mathit{ab}}+3H{\dot{\varpi }}_{\mathit{ab}}-H^2{\varpi }_{\mathit{ab}}=0\end{array}$$(36)

for the linear evolution of the peculiar vorticity. Clearly, an exactly analogous differential formula also monitors the linear peculiar shear (ς_ab). Recalling that H = 2/3t during the Einstein-de Sitter epoch, when the above equation is solved, it gives

$$\begin{array}{c}\hfill \varpi ={\mathcal{C}}_1{t}^{1/3}+{\mathcal{C}}_2{t}^{-4/3}={\mathcal{C}}_3{a}^{1/2}+{\mathcal{C}}_4{a}^{-2},\end{array}$$(37)

with the same solution governing the evolution of the peculiar shear (ς) as well. Therefore, both ϖ and ς grow considerably slower than the peculiar-velocity field itself (we recall that v ∝ t – see solution (30) in Sect. 3.3.3 before). This result is in agreement with the earlier relativistic analysis of Tsaprazi & Tsagas (2020).

3.3.6. Impact of the late-time acceleration

So far our analysis has been confined to the Einstein-de Sitter epoch of the universe. Within the framework of the ΛCDM scenario, the earlier dust-dominated era is followed by a recent phase of accelerated expansion, where the energy density of the cosmos is typically dominated by a cosmological constant or by dark energy. Therefore, our results apply to the period between matter–radiation equality and the onset of cosmic acceleration. Once the latter epoch has started, the evolution of the peculiar-velocity field is expected to change drastically. More specifically, cosmic acceleration is expected to suppress the growth of peculiar velocities, in the same way it suppresses essentially all types of perturbations, such as those in the matter density.

Analytical relativistic studies that directly address and quantify the impact of a late ΛCDM phase on large-scale peculiar velocities would be very useful. At present, there is support for the aforementioned suppressing effect of the universal acceleration, but it is indirect. In particular, according to the Newtonian numerical work of Demchenko et al. (2016), dark energy inhibits the growth of peculiar velocities compared to the Einstein-de Sitter epoch. Studies looking into the linear evolution of peculiar velocities during a phase of de Sitter inflation have also shown exponential decay for the v field (Maglara & Tsagas 2022). Nevertheless, even if peculiar velocities are suppressed during the late Λ-dominated phase, at least some of the strong growth of the preceding Einstein-de Sitter period should survive through to the present.

On these grounds, one should expect to measure bulk velocities that are faster than anticipated by the Newtonian studies at higher redshifts but to also see a decrease in their mean value at lower redshifts (see Tsagas 2026; Tsagas et al. 2026 for further discussion). We note that a peculiar-velocity profile with increased values at higher redshifts (i.e. at earlier times) and a late-time suppression has been reported in the survey of Colin et al. (2011) and (more pronounced) in those of Watkins et al. (2023), Watkins & Feldman (2025). Put another way, our work appears to provide theoretical support to the aforementioned observational studies.

4. An effective gravitational potential

As we have already mentioned, it is possible to recover the relativistic solution through Newtonian analysis, although the mapping is neither unique nor physically equivalent in a strict sense. To do so, one needs to introduce a specific ansatz for the evolution of the gravitational potential (or, equivalently, for the integration constants that appear when Poisson’s equation is used in reverse). More specifically, in line with the covariant Newtonian study (see Sect. 3.3.1), the divergence of Eq. (20) gives

$$\begin{array}{c}\hfill {\partial }^2\dot{\varphi }=\frac{1}{2}H{\partial }^2\varphi +\frac{3}{2}H{\partial }^{\alpha }{\dot{v}}_{\alpha },\end{array}$$(38)

which upon integrating once in space yields

$$\begin{array}{c}\hfill {\partial }_{\alpha }\dot{\varphi }=\frac{1}{2}H{\partial }_{\alpha }\varphi +\frac{3}{2}H{\dot{v}}_{\alpha }+V_{\alpha }.\end{array}$$(39)

Here V_α could (in principle) have both a temporal and a spatial dependence, and it encodes the solenoidal (divergence-free) integration mode that is left undetermined by the divergence in Eq. (38) and is fixed only after specifying boundary conditions and/or a coordinate choice. When substituting the above equation into the right-hand side of Eq. (19), one arrives at the following (non-homogeneous) differential equation

$$\begin{array}{c}\hfill {\ddot{v}}_{\alpha }+3H{\dot{v}}_{\alpha }-H^2{v}_{\alpha }=-V_{\alpha }\end{array}$$(40)

for the linear evolution of the peculiar-velocity field. Expression (40) differs from Eq. (21) only by its non-zero right-hand side, the presence of which allows one to change the standard Newtonian evolution of the linear peculiar velocities. By allowing V_α ≠ 0, we effectively parameterise the freedom that the time evolution of ∇ϕ is not fixed by the Poisson equation alone. For instance, by choosing $V_{\alpha }=-(5/2)H{\dot{v}}_{\alpha }-H^2{v}_{\alpha }$, one can recover the relativistic growth rate of v ∝ t obtained previously in Sect. 3.3.3 (see Eqs. (29) and (30) there). It goes without saying, however, that our choice of V_α is not uniquely determined, and therefore this procedure should be viewed as illustrative mapping rather than as a physical equivalence between Newtonian gravity and general relativity.

It is worth pointing out that a vector function analogous to V_α was obtained in Peebles (1980) also after involving and integrating Poisson’s formula. There, the aforementioned vector was set to zero by an appropriate choice of coordinates or boundary conditions. This freedom is reminiscent of the fact that the relativistic peculiar flux is also a vector that can be set to zero (albeit only locally) by moving to the energy (or Landau-Lifshitz) frame. Moreover, even in the energy frame, the peculiar energy flux does not disappear but re-emerges ‘disguised’ as a ‘peculiar particle flux’. We remind the reader that in relativity and depending on the specifics of the problem in hand, one may adopt the energy frame where the energy flux vanishes but there is non-zero particle flux. Alternatively, one may choose to use the particle (or Eckart) frame when the situation is reversed (e.g. see Tsagas et al. 2008; Ellis et al. 2012).

5. Discussion

Large-scale peculiar motions, commonly referred to as bulk flows, are commonplace in our Universe and have been repeatedly verified by many surveys (e.g. see Nusser & Davis 2011; Ma & Scott 2013; Scrimgeour et al. 2016; Qin et al. 2019; Watkins et al. 2009; Colin et al. 2011; Salehi et al. 2021; Watkins et al. 2023; Whitford et al. 2023 for a representative, though incomplete, list). Although the aforementioned flows typically encompass domains of up to several hundred megaparsecs, we believe that they started as weak peculiar-velocity perturbations around the time of recombination that grew larger and faster by structure formation and by the increasing non-uniformity of the post-recombination Universe. Nevertheless, the evolution and the implications of these bulk peculiar motions are still largely unknown. Moreover, during the roughly past fifteen years the bulk-flow puzzle has become even more perplexing, after an increasing number of surveys have reported sizes and speeds well in excess of those allowed by the concordance cosmological model, namely, by the ΛCDM paradigm. In fact, if one is allowed to take into account the controversial dark flows, with sizes of the order of 1000 Mpc and speeds close to 1000 km/s (see Kashlinsky et al. 2008, 2011 for example), the situation becomes rather untenable.

Significantly, the ΛCDM constraints on the speed and the depth of the observed bulk peculiar flows are based entirely on Newtonian theory, and there is a serious difference between Newtonian gravity and general relativity when it comes to the study of peculiar motions, a difference which stems from the fundamentally different way the two theories treat the gravitational field and its sources. To begin with, we recall that bulk flows are nothing but matter in motion and that moving matter means a non-zero energy flux. In Newton’s theory only the density of the matter contributes to the gravitational field, via Poisson’s equation, while any energy fluxes that may exist have no gravitational input. In relativity, on the other hand, the pressure (both isotropic and viscous), as well as the fluxes also contribute to the energy momentum tensor and therefore to the (relativistic) gravitational field. It is the input of the peculiar flux that has been bypassed in essentially all the mainstream studies so far. Thus, accounting for the gravitational contribution (or not) is the key difference that separates the relativistic treatment from the Newtonian and quasi-Newtonian treatments of cosmological peculiar motions.

By not incorporating the gravitational input of the peculiar flux, Newtonian and quasi-Newtonian analyses have both led to the same mediocre (v ∝ t^1/3) growth rate for the linear peculiar-velocity field (Peebles 1976, 1980; Maartens 1998; Ellis et al. 2001). Unfortunately, the v ∝ t^1/3 rate is too slow to explain the reported fast bulk flows. However, had the gradient of the energy-conservation law (16a) been included in Maartens (1998), Ellis et al. (2001), these quasi-Newtonian studies would also have arrived at the relativistic expression (26) for the 4-acceleration, instead of the Newtonian (23) – see Sect. 3.3.2 and 3.3.3. In that case, the quasi-Newtonian and the relativistic treatments would have arrived at the same results and any subsequent comparison between them would have been entirely unnecessary.

On the other hand, by including the peculiar-flux contribution to the gravitational field alone, the relativistic analysis has led to the considerably faster growth rate of v ∝ t (see Tsagas 2026 and also Sect. 3.3.3 here). Moreover, the standard theory of differential equations guarantees that v ∝ t is the minimum growth supported by the relativistic studies. To a large extent, the latter conclusion was already foreseen in the studies of Tsaprazi & Tsagas (2020)–Miliou & Tsagas (2024), which arrived at the faster linear growth rate of v ∝ t^4/3, after accounting for some of the inhomogeneity effects as well.

A complementary explanation for the aforementioned relativistic results involves the acceleration and the four-acceleration, which are the driving forces of the linear peculiar-velocity field in all three approaches. A closer inspection (see Sect. 3.3.4) revealed that the driving agent decays as ∝t^−2/3 in both the Newtonian and the quasi-Newtonian studies. In contrast, the relativistic four-acceleration is constant in time (at a minimum; see Tsagas 2026 and Sect. 3.3.3 here), or it could potentially increase as t^1/3 throughout the Einstein-de Sitter phase of the Universe (see Tsaprazi & Tsagas 2020 and Miliou & Tsagas 2024). This drastic difference explains why general relativity supports significantly faster peculiar-velocity fields. We also note that within the framework of the standard ΛCDM model, cosmic acceleration starts late into the Einstein-de Sitter phase of the Universe, and some (at least) of the increased velocity growth taking place during the Einstein-de Sitter epoch should survive through the subsequent accelerated phase and have consequently led to faster and deeper residual bulk flows today, such as those reported in Watkins et al. (2009), Colin et al. (2011), Salehi et al. (2021), Watkins et al. (2023), Whitford et al. (2023). Moreover, all of this happens within the standard ΛCDM model and without the need for any new free parameters.

The exact analytical results of the linear analysis presented here as well as in the earlier works of Tsaprazi & Tsagas (2020), Maglara & Tsagas (2022), Miliou & Tsagas (2024), and Tsagas (2026) can also provide a useful testing ground for the emerging numerical studies simulating the post-recombination universe. The fact that all the aforementioned theoretical treatments are gauge invariant is an additional advantage since their results are free of any gauge-related ambiguities⁷. We recall that (to the best of our knowledge) the numerical simulations rely on particular gauge choices when addressing the problems in hand.

Up to now, most of the simulations that employ relativistic techniques do not consider the evolution of peculiar velocities, with a couple of possible exceptions that looked at the effects of the peculiar vorticity (e.g. Barrera-Hinojosa et al. 2021). The vorticity was found to have a negligible impact, which was to be expected since it grows considerably slower that the peculiar velocity itself (even in the relativistic analysis; see Tsaprazi & Tsagas 2020 and also Sect. 3.3.5 here). Moreover, the simulations typically operate within the ΛCDM paradigm (e.g. see Aurrekoetxea et al. 2025 for a recent review), so the numerical codes are adapted to the recent accelerated phase of the Universe, where the peculiar-velocity field is suppressed (e.g. see Demchenko et al. 2016), instead of the preceding Einstein-de Sitter epoch, where peculiar velocities are strongly enhanced. As a result, so far there is no common ground to facilitate direct comparison between analytical and the numerical work on this topic. Nevertheless, the relativistic simulations are still under development, and there is room for testing and improvement, especially in view of repeated surveys reporting bulk flows in excess of the (Newtonian-based) ΛCDM limits. Hopefully, the results of the analytical relativistic studies will provide the motivation for extending and adapting the numerical codes to address the ongoing bulk-flow puzzle.

Observationally, the time evolution of the peculiar velocities can be mapped into a redshift evolution at low z, enabling a comparison between our theoretical predictions and the standard Newtonian expectation using peculiar-velocity surveys such as CosmicFlows-4. However, peculiar-velocity measurements are subject to a well-known degeneracy between the assumed background cosmology and local kinematic contributions. This degeneracy can be partially mitigated through independent distance indicators, such as Type Ia supernovae or Tully-Fisher relations, which provide redshift-independent distance estimates. Nevertheless, the use of such indicators introduces additional uncertainties that grow with redshift, thereby impacting the reconstruction of the peculiar-velocity field and limiting the robustness of the physical inferences. A detailed implementation of unbiased and statistically robust estimators – such as the minimum-variance reconstruction techniques applied to CosmicFlows-4 (Watkins et al. 2023; Watkins & Feldman 2025) – is beyond the scope of the present work and will be addressed in future studies.

In summary, general relativity can provide a simple and physically motivated theoretical solution to the ongoing bulk-flow puzzle by simply taking into account the (purely general relativistic) gravitational contribution of the peculiar flux. This can in turn relax the current (purely Newtonian) peculiar-velocity expectations of the ΛCDM model and, in so doing, allow the latter to naturally accommodate bulk flows faster and deeper than anticipated, such as those reported in Watkins et al. (2009), Colin et al. (2011), Salehi et al. (2021), Watkins et al. (2023), Whitford et al. (2023). The perceived discrepancy between the current cosmological model and the above bulk-flow reports thus may not reflect a generic problem in the ΛCDM paradigm, but instead indicate the use of an inappropriate gravitational theory in our studies of cosmological peculiar flows.

Acknowledgments

CGT acknowledges support from 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).

References