© The Authors 2024

1 Introduction

The comprehension of the recessional velocity of extragalactic objects hinges upon the redshift-distance relation derived from homogeneous solutions of Einstein’s field equations within the framework of the Friedmann–Lemaître–Robertson-Walker (FLRW) metric. Within this coordinate system, the spatial coordinates are comoving with the matter that fills the Universe, attributing the cosmological redshift to the expansion of ‘space’ along the path of photons. However, observations in the nearby Universe have revealed a notable scatter around the linear redshift-distance relation, primarily driven by the presence of peculiar velocities of sources inducing a Doppler shift upon photon emission (Peebles 1993; Peacock 1998). The debate in the literature regarding whether cosmological and Doppler shifts are equivalent phenomena is ongoing (Bunn & Hogg 2009; see also Peacock 1998, 2008; Whiting 2004; Chodorowski 2007). As noted in Bunn & Hogg (2009), the interpretation of space expansion is not entirely satisfactory from a relativistic standpoint, as it relies on a specific choice of coordinates and thus lacks covariance. Conversely, Räsänen (2009) introduced a covariant approach to redshift calculations for an inhomogeneous dust universe. While this approach represents a significant advancement in our understanding of the nature of redshifts in a manner consistent with the covariant framework of general relativity, it has been relatively overlooked. This is largely because, on average, the results do not seem to deviate from the background FLRW expansion in a statistically homogeneous and isotropic universe (Räsänen 2009; Rasanen 2010).

This paper aims to re-evaluate the cosmological redshift by addressing the influence of inhomogeneities and irreversible processes. In Sec. 2, we begin by providing a detailed analysis of the standard interpretation of the redshift within the FLRW metric and discuss the limitations of this approach by using non-comoving coordinates. In the subsequent section (Sec. 3), building on Räsänen (2009), we introduce a covariant approach to the redshift calculation, emphasising the role of peculiar motions and local volume expansion factors. We propose a new approach to the redshift calculation by introducing a tracer of time-like geodesics: the cosmological redshift is associated with the local expansion rate along these time-like geodesics, while the Doppler shifts are associated with the potential non-geodesic motions of the source and observer. Following this, in Sec. 4 we present our findings on the impact of irreversible processes on the cosmological redshift. We argue that in scenarios characterised by irreversible evolution, such as the dissipation of large-scale compression modes through matter virialisation and associated entropy production, the averaged expansion factor is dominated by expansion in voids. This leads to significant violations of the second law of thermodynamics if the Universe’s evolution adheres to the background evolution on average. Lastly, in Sec. 5 we discuss the broader implications of our findings for cosmological models and observations. We conclude in Sec. 6 on the necessity of re-evaluating the influence of inhomogeneities and irreversible processes to better understand the complex dynamics of the Universe.

2 Redshift in a homogeneous expanding universe in non-comoving coordinates

We began with a standard flat FLRW metric with a scale factor, a(t), such that $$d{s}^2=-c^{2}d{t}^2+a^2(d{r}^2+r^{2}d{\mathrm{\Omega }}^2).$$(1)

The proper distance, R, is $$R=a(t)r,$$(2)

and the total velocity is then $$\frac{dR}{dt}=H(t)R+a(t)\frac{dr}{dt},$$(3)

with $H\equiv \dot{a}/a$. Within the total velocity, one can discern a contribution originating from the recession velocity due to the expansion of space, and another contribution arising from peculiar velocities v_p = a(t)dr/dt. In the late evolution of the Universe, peculiar motions are non-relativistic in the weak field limit. Perturbation theory then leads to the following formulation of the redshift with a cosmological contribution from space expansion and a contribution from the Doppler shift associated with peculiar motions: $$1+z\approx \exp \left({\int }_{t_{\text{src }}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}Hdt\right)(1+\frac{v_p}{c}),$$(4)

assuming for simplicity that the observer at time t_dst has no peculiar motion, and the source at time t_src a peculiar motion v_p at emission. With this approach, the inhomogeneities have no significant impact on the redshift, provided that the peculiar velocities remain small v_p ≪ c.

We next assumed a perfectly homogeneous FLRW universe with no peculiar motion v_p = 0. We wished to apply a local coordinate transformation in a non-comoving frame, characterised by an expansion factor a_s(t) different from a(t) such that $$\begin{array}{rl}{r}_d& =r{a}_d(t),\\ t_d& =t+\alpha (r,t),\end{array}$$(5)

with a_s = a/a_d and α(r, t) a function that will be specified later, such that t_d ≈ t at the leading order in $O\left({\dot{a}}_{d}r/c\right)$. We used such a coordinate transformation, locally, such that ${\dot{a}}_{d}r/c$ is indeed much smaller than unity. With this transformation, the proper distance R is $$R=a(t)r=a_s(t)r_d,$$(6)

with a_s = a/a_d. The total velocity is then $$\frac{dR}{dt}=H_s(t)R+a_s(t)\frac{d{r}_d}{d{t}_d},$$(7)

with $H_s\equiv {\dot{a}}_s/a_s$. Within the total velocity, one can discern now a contribution originating from the recession velocity due to the expansion of ‘space’, characterised by the scale factor a_s, and another contribution arising from non-comoving velocities v_p,r = a_sH_dr_d.

We chose the function α(r, t) such that the metric in (r_d, t_d) coordinates has only perturbative deviations from a homogeneous metric with a scale factor a_s. The coordinate transformation verifies $$\begin{array}{rl}d{r}_d& =a_{d}dr+r{\dot{a}}_{d}dt,\\ d{t}_d& =(1+\dot{\alpha })dt+{\alpha }^{\mathrm{\prime }}dr,\end{array}$$(8)

with $\dot{\alpha }$ the derivative with respect to t and α′ the derivative with respect to r. This system can be inverted to give $$\begin{array}{rl}{a}_{d}dr& =((1+\dot{\alpha })d{r}_d-r{\dot{a}}_{d}d{t}_d)/(1+\dot{\alpha }-r{H}_d{\alpha }^{\mathrm{\prime }}),\\ dt& =(d{t}_d-{\alpha }^{\mathrm{\prime }}/a_{d}d{r}_d)/(1+\dot{\alpha }-r{H}_d{\alpha }^{\mathrm{\prime }}),\end{array}$$(9)

with $H_d={\dot{a}}_d/a_d$. By injecting these equations into ds², the nondiagonal metric coefficient is $$g_{r_d{t}_d}=2(c^2{\alpha }^{\mathrm{\prime }}/a_d-a_s^2(1+\dot{\alpha })r{\dot{a}}_d)/{(1+\dot{\alpha }-r{H}_d{\alpha }^{\mathrm{\prime }})}^2.$$(10)

We imposed $$\alpha (r,t)=\frac{1}{2}\frac{a^2{H}_d{r}^2}{c^2},$$(11)

such that the metric in (r_d, t_d) coordinates is given by $$d{s}^2=-c^{2}Ad{t}_d^2+a_s^2(Bd{r}_d^2+r_d^{2}d{\mathrm{\Omega }}^2)-2cd{t}_d{a}_{s}d{r}_{d}C,$$(12)

with $$\begin{array}{rl}& A=(1-a_s^2{H}_d^2{r}_d^2/c^2)/{(1+\dot{\alpha }-a_s^2{H}_d^2{r}_d^2/c^2)}^2,\\ & B=((1+\dot{\alpha }{)}^2-a_s^2{H}_d^2{r}_d^2/c^2)/{(1+\dot{\alpha }-a_s^2{H}_d^2{r}_d^2/c^2)}^2,\\ & C=\dot{\alpha }{a}_s{H}_d{r}_d/c/{(1+\dot{\alpha }-a_s^2{H}_d^2{r}_d^2/c^2)}^2.\end{array}$$(13)

At the leading order in O(a_sH_dr_d/c), we get $$\begin{array}{rl}\alpha & =O\left(a_s^2{H}_d{r}_d^2/c^2\right)\\ \dot{\alpha }& =O\left(a_s^2{H}_d^2{r}_d^2/c^2\right)\\ {\alpha }^{\mathrm{\prime }}& =O\left(a_s^2{H}_d{r}_d/c^2\right)\\ A& =1+O\left(a_s^2{H}_d^2{r}^2/c^2\right)\\ B& =1+O\left(a_s^2{H}_d^2{r}^2/c^2\right)\\ C& =O\left(a_s^3{H}_d^3{r}_d^3/c^3\right).\end{array}$$(14)

The metric in (r_d, t_d) coordinates has therefore only perturbative deviations from a homogeneous metric with a scale factor a_s. In this coordinate system and using the same computation of the redshift as usually done for weak perturbations around a homogeneous metric, at the leading order in O(a_sH_dr_d/c), the redshift of light emitted by a source at r_d = a_dr_src, t_d = t_src towards an observer at r_d = 0, t_d = t_dst can be expressed as follows: $$1+z\approx \exp \left({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}{H}_{s}dt\right)(1+\frac{H_d{a}_s{a}_d{r}_{\mathrm{s}\mathrm{r}\mathrm{c}}}{c}),$$(15)

with two contributions: one stemming from the ‘space’ expansion along the photon path of the homogeneous metric with scale factor a_s(t) and another originating from a Doppler term from non-comoving velocities only at emission. At the leading order in O(H_sr/c), this is equivalent to z ≈ (H_s + H_d)a_sa_dr_src/c = Har_src/c only when H_S is independent of time. This method of calculating redshift is not covariant, as it essentially converts part of the cosmological redshift along the photon’s path into a Doppler shift that occurs solely at the point of emission.

3 Covariant computation of the redshift in an inhomogeneous universe

Following Räsänen (2009), we began by defining a set of observers that are tracing time-like geodesics, characterised by their four-velocity (denoted as u^α). Consequently, these observers satisfy the conditions u^β∇βu^α = 0 and u^αu_α = −1. It is important to note that these observers may not necessarily move in tandem with the matter filling the Universe, as this matter can experience acceleration due to non-gravitational forces and may not strictly adhere to time-like geodesics. However, especially on large scales and within the context of dark matter, we can consider the dark matter fluid to be a useful tracer of these time-like geodesics.

To further elucidate on this, we introduced the projection tensor, denoted as h_αβ, which operates within the tangent space orthogonal to u^α. It is defined as h_αβ = g_αβ + u_α u_β. With this tensor, we decomposed the covariant derivative of u^α, $${\mathrm{\nabla }}_{\beta }{u}_{\alpha }=\frac{1}{3}{h}_{\alpha \beta }\theta +{\sigma }_{\alpha \beta }+{\omega }_{\alpha \beta },$$(16)

with θ = ∇_αu^α, the volume expansion rate, ${\sigma }_{\alpha \beta }={\mathrm{\nabla }}_{(\alpha }{u}_{\beta )}-h_{\alpha \beta }\theta /3$ the shear tensor, and ${\omega }_{\alpha \beta }={\mathrm{\nabla }}_{[\beta }{u}_{\alpha ]}$ the vorticity tensor. For simplicity, we assume in the rest of this paper that the shear and vorticity tensors can be neglected.

We next defined k^α, the tangent vector of the null geodesics that satisfies $k^{\beta }{\mathrm{\nabla }}_{\beta }{k}^{\alpha }=0$ and k^αk_α = 0. The tangent vectors u^α and k^α are parallel propagated with respect to the time-like and null geodesics, respectively, but not with respect to each other. Consequently, the photon momentum changes along the timelike geodesics and the redshift is defined by $$1+z=\frac{E_{\mathrm{s}\mathrm{r}\mathrm{c}}}{E_{\mathrm{d}\mathrm{s}\mathrm{t}}},$$(17)

with E_src the photon energy at emission by the source and E_dst the photon energy at the location of the observer. The energy can be computed from E = −u_αk^α, and following Räsänen (2009) we decomposed k^α into a component parallel and a component orthogonal to the time-like geodesics k^α = E(u^α + e^α) with u_αe^α = 0 and e_αe^α = 1 (and consequently h_αβe^α e^β = 1). The evolution of the energy along the null geodesic can then be followed with the affine parameter λ, $$\begin{array}{rl}{\mathrm{\partial }}_{\lambda }E& \equiv k^{\alpha }{\mathrm{\nabla }}_{\alpha }E\\ & =-k^{\alpha }{k}^{\beta }{\mathrm{\nabla }}_{\alpha }{u}_{\beta }\\ & =-E^2{e}^{\alpha }{e}^{\beta }{\mathrm{\nabla }}_{\alpha }{u}_{\beta }\\ & =-\frac{E^2}{3}\theta .\end{array}$$(18)

Without vorticity, the hypersurfaces of constant proper time are orthogonal to u^α, and t(λ) is monotonic. We could then invert the relation between λ and t to obtain dλ = dt/E; hence, $$1+z=\exp \left({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}\frac{\theta (t,\mathit{x}(t))}{3}dt\right).$$(19)

A more detailed demonstration including vorticity and shear can be found in Räsänen (2009). The main difference here is that we take a set of observers following time-like geodesics that are not necessarily co-moving with matter filling the universe. For a homogeneous metric with scale factor $\tilde{a}(t),\tilde{H}(t)=\dot{\tilde{a}}/\tilde{a}$, and for small three-velocity ${\stackrel{\mathbf{~}}{\mathit{u}}}_{\mathit{g}}$, the volume expansion factor is given by $\theta \approx 3\tilde{H}+\mathbf{\nabla }\cdot {\tilde{\mathit{u}}}_{\mathit{g}}$. Applying Eq. (19) to the non-comoving coordinates defined above, we get $$\begin{array}{rl}1+z& \approx \exp \left({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}\frac{3H_s+\mathbf{\nabla }\cdot \left(d{r}_d/d{t}_d{\mathit{e}}_{r_d}\right)}{3}dt\right)\\ & \approx \exp \left({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}(H_s+H_d)dt\right).\end{array}$$(20)

With this covariant approach, ‘space-time’ geometry characterised by the expansion rate H_s and non-comoving velocities of observers following time-like geodesics characterised by H_dr_d contribute both as an expansion factor along the trajectory of the photons. For a homogeneous FLRW solution with a Hubble expansion rate H, the comoving coordinate system would lead to $\tilde{H}=H$ and ${\stackrel{\mathbf{~}}{\mathit{u}}}_{\mathit{g}}=0$ and the non-comoving coordinate system defined in Sec. 2 to $\tilde{H}=H_s$ and ${\tilde{\mathit{u}}}_{\mathit{g}}=d{r}_d/d{t}_d{\mathit{e}}_{r_d}$. In both coordinate systems, a consistent calculation of the cosmological redshift is obtained with θ = 3H.

As a simple example, we considered a universe with G = 0 such that H = 0 and a(t) = 1. In such a universe the cosmological redshift should be exactly zero, since there is no expansion along the photon path and there should be only a Doppler shift from the velocity difference between the source and the observer. We assumed for the sake of simplicity that they are both at rest such that the total redshift is zero. We also considered an arbitrary coordinate system with a scale factor a_s(t) and an expansion rate H_s(t), we thus have a_d(t) = 1/a_s(t) and H_d(t) = −H_s(t). In such a coordinate system, the perturbative redshift calculation gives $$\begin{array}{rl}& 1+z\approx \exp ({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}{H}_s(t)dt)(1-\frac{H_s\left(t_{\mathrm{s}\mathrm{r}\mathrm{c}}\right)r_{\mathrm{s}\mathrm{r}\mathrm{c}}}{c}),\\ & 1+z\approx (1+{\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}(H_s(t)-H_s\left(t_{\mathrm{s}\mathrm{r}\mathrm{c}}\right))dt),\end{array}$$(21)

which demonstrates that the Doppler shift at emission cannot compensate for the expansion of the coordinate system along the photon path as soon as this expansion is time dependent. However, with the covariant calculation, the redshift is exactly zero since θ = 3H_s(t) − ∇(H_s(t)r_d) = 0. This further demonstrates that the perturbative formulation of the redshift with the Hubble expansion rate H instead of the local expansion rate θ along time-like geodesics is not covariant.

Lastly, it is imperative to highlight the precise physical entity that should be traced by the local volume expansion rate and ${\stackrel{\mathbf{~}}{\mathit{u}}}_{\mathit{g}}$. It is crucial to note that the local volume expansion rate should not mirror the fluid velocity but rather correspond to the velocity field of a hypothetical set of freefalling observers following time-like geodesics (Räsänen 2009). The matter filling the universe can only serve as a tracer of these time-like geodesics if the acceleration due to non-gravitational forces can be deemed negligible. Hence, we advocate for the use of expressions that avoid the term ‘expansion of space’, as it lacks covariance, and instead, we refer to the expansion along time-like geodesics when discussing cosmological redshift. Then, the Doppler shifts can be defined as originating from the potential non-geodesic motion of the source and observer at emission and reception.

In an inhomogeneous universe with an arbitrary expanding coordinate system defined by $\tilde{a}(t),\tilde{H}(t)$, the general expression of the redshift can therefore be defined as follows: $$\begin{array}{rl}1+z=& \exp \left({\int }_{t_{\mathrm{s}\mathrm{r}\mathrm{c}}}^{t_{\mathrm{d}\mathrm{s}\mathrm{t}}}\frac{\theta (t,\mathit{x}(t))}{3}dt\right)\\ & \times (1+\frac{v_{\mathrm{D}\mathrm{o}\mathrm{p}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{r}}(t_{\mathrm{s}\mathrm{r}\mathrm{c}},\mathit{x}\left(t_{\mathrm{s}\mathrm{r}\mathrm{c}}\right))}{c}-\frac{v_{\text{Doppler }}(t_{\mathrm{d}\mathrm{s}\mathrm{t}},\mathit{x}\left(t_{\mathrm{d}\mathrm{s}\mathrm{t}}\right))}{c}),\end{array}$$(22)

with $$\begin{array}{c}\theta (t,\mathit{x})=3\tilde{H}(t)+\mathbf{\nabla }\cdot {\tilde{\mathit{u}}}_{\mathrm{g}},\\ v_{\text{Doppler }}(t,\mathit{x})=\tilde{a}(t)({\tilde{\mathit{u}}}_{\text{fluid }}-{\tilde{\mathit{u}}}_{\mathrm{g}})\cdot \mathit{n}.\end{array}$$(23)

Here, we define the cosmological redshift as the integral along the line of sight of the local expansion rate along time-like geodesics, and the Doppler shifts as the Doppler terms from non-geodesic motions of the source and the observer in the non-relativistic limit. These definitions are now fully covariant as θ and the velocity difference in the Doppler shift are invariant under a change of the coordinate system. Fundamentally, one can see that the cosmological contribution should not transform into a Doppler one or vice versa under a change of the coordinate system. We considered only the expansion rate in Eq. (22) since vorticity does not contribute to the redshift and shear is negligible on large scales if structures have no preferred orientation (Räsänen 2009).

At large scales and especially with dark matter, the matter within the Universe is freefalling and can be taken as a tracer of the dynamics along time-like geodesics. However, in the cosmological standard model, the coordinate system is comoving with the FLRW background and is not comoving with these time-like geodesics. Consequently, the conventional redshift calculation overlooks the inhomogeneous dynamics along these geodesics, attributing a Doppler shift solely to peculiar motions at the point of emission, mirroring our example using non-comoving coordinates within a homogeneous FLRW universe. In contrast, the covariant calculation of the redshift highlights that peculiar motions responsible for the formation of cosmic structures do indeed contribute to the redshift along the photon path through the local expansion rate.

This insight allowed us to reconsider the influence of inhomogeneities on the cosmological redshift. A backreaction is not a question of small peculiar velocities relative to the speed of light. It requires a deeper inquiry: whether, in a statistically homogeneous and isotropic universe, the averaged local volume expansion rate along time-like geodesics provides a close approximation to the background expansion 3H, which might not be the case even if the Doppler shifts of the source and observer relative to time-like geodesics are negligible.

4 Impact of irreversible processes on the cosmological redshift

The condition ⟨θ⟩ ≈ 3H can be expressed equivalently as ⟨∇ · u_g⟩ ≈ 0. When we start from a homogeneous universe where u_g = 0 and consider a reversible evolution, one might naturally expect, for example based on hydrodynamics equations, that the condition ⟨∇ · u_g⟩ ≈ 0 would hold. Essentially this implies that inhomogeneous expansion and compression modes evolve smoothly and compensate for one another on average, resulting in ⟨θ⟩ ≈ 3H (see the top panel of Fig. 1). It is a common belief that the cosmological principle (i.e. homogeneity and isotropy of the Universe at large scales) is sufficient to ensure such a condition. In previous investigations addressing the back-reaction problem (Buchert & Ehlers 1997; Buchert 2008), which examines the impact of inhomogeneities on the cosmological redshift, it has been generally assumed that this condition holds (Buchert & Ehlers 1997; Räsänen 2009). Consequently, it would appear that back-reaction is negligible in the non-relativistic limit (Ishibashi & Wald 2006; but see also Buchert et al. 2015). However, we aimed to re-evaluate this conclusion since structure formation does not have inhomogeneous expansion and compression modes that compensate for each other on average: compression modes are dissipated by virialisation and the associated entropy production and cannot compensate for the expansion of the voids. Hence, in the presence of irreversible processes, ⟨θ⟩ ≠ 3H. The link between expansion and entropy production can be made explicit by using the second principle of thermodynamics assuming that a quasi-equilibrium defined by a temperature T and a stress tensor with pressure σ has been reached by virialisation, $$\begin{array}{rl}TdS& =dE+\sigma dV,\\ Tds& =de+Z{k}_{b}Td(\log V),\end{array}$$(24)

for a closed system of N particles (dN=0) with an entropy per particle s = S/N, an internal energy per particle e = E/N, and a compressibility factor Z = σV/Nk_bT. Assuming the evolution of an irreversible system from an initial state and the evolution of a reversible system that reaches the same internal energy (e), temperature (T), and compressibility factor (Z) and has the same entropy as the initial state, $$\begin{array}{rl}T\mathrm{\Delta }s& =\mathrm{\Delta }e+Z{k}_{b}T\log \left(V_{\text{irrev }}/V_{\text{ini }}\right),\\ 0& =\mathrm{\Delta }e+Z{k}_{b}T\log \left(V_{\text{rev }}/V_{\text{ini }}\right),\end{array}$$(25)

implies that $$\frac{V_{\text{irrev }}}{V_{\text{ini }}}=\frac{V_{\text{rev }}}{V_{\text{ini }}}\exp \left(\frac{\mathrm{\Delta }s}{k_{b}Z}\right).$$(26)

We recall that the expansion rate can be directly linked to the variation of volume by using the equation of mass conservation with a Lagrangian time derivative: $$\frac{d\log V}{dt}=\theta ,$$(27)

which implies the following expression of the redshift for motions that follow time-like geodesics in an irreversible system (assuming spatial ergodicity along the photon path), $$\begin{array}{rl}1+z& ={\left(\frac{V_{\text{irrev }}\left(t_{\text{dst }}\right)}{V_{\text{ini }}\left(t_{\text{src }}\right)}\right)}^{1/3},\\ & ={\left(\frac{V_{\text{rev }}}{V_{\text{ini }}}\right)}^{1/3}\exp \left(\frac{\mathrm{\Delta }s}{3k_{b}Z}\right),\\ & =\frac{a\left(t_{\text{dst }}\right)}{a\left(t_{\text{src }}\right)}\exp \left(\frac{\mathrm{\Delta }s}{3k_{b}Z}\right),\end{array}$$(28)

where a(t) is the scale factor of the reversible background evolution. Equation (28) demonstrates the direct link between entropy production and the cosmological redshift in an irreversible system. The second principle of thermodynamics states that the entropy production is always positive with Δs > 0, the cosmological principle is therefore not sufficient to ensure that the expansion of an irreversible system follows the background expansion without violation of the second principle.

In an adiabatic system, the occurrence of an irreversible process and the subsequent creation of entropy is possible through the dissipation of large-scale motions into microscopic motions. This process, driven by virialisation, is referred to as ‘violent relaxation’ within the context of structure formation in the Universe, in the presence of non-collisional dark matter (White 1996). When a significant portion of the large-scale compression modes dissipates due to virialisation, they can no longer compensate for the expansion occurring in cosmic voids, inevitably leading to the conclusion that ⟨θ⟩ ≠ 3H (see the bottom panel of Fig. 1). From numerical simulations (Haider et al. 2016), we estimated that 50% of the dark matter mass is completely virialised in haloes and 45% virialised in two directions in filaments (Eisenstein et al. 1997). A realistic universe has therefore virialised most of its mass and almost completely the initial compression modes that appeared during the formation phase of the large-scale structures. Furthermore, it is crucial to emphasise that an inhomogeneous system undergoing irreversible processes cannot, on average, precisely mimic the behaviour of the FLRW background. The background evolves akin to a reversible system (constant entropy), while an inhomogeneous system subject to entropy generation would essentially need to violate the second law of thermodynamics with processes generating negative local entropy to align its averaged entropy with the background entropy. Additionally, the process of dark matter virialisation entails the acquisition of microscopic stress through microscopic virialised motions. In this context, the dark matter fluid ceases to be truly pressureless after virialisation, resulting in a nonzero stress tensor within the Jeans equations. This can be seen as a phase transition to a bound state similar to the formation of a liquid but in the context of statistical mechanics for non-ideal self-gravitating fluids (Tremblin et al. 2022). A central point to note is that the equilibrium observed in large-scale structures is not a balance between gravitational attraction and stresses induced by non-gravitational interactions. Dark matter, by definition, is a tracer for dynamics along time-like geodesics. Therefore, its virialisation marks the dissipation of large-scale dynamics along these geodesics into microscopic motions that continue to follow geodesic trajectories in freefall. This fundamental aspect underscores why the virialisation of large-scale structures inevitably influences the cosmological redshift.

The presence of entropy creation, characterised by ⟨∇ · u_g⟩ ≠ 0, necessarily implies the existence of discontinuities in the three-velocity field. Remarkably, dark matter virialisation occurs subsequent to shell crossing, resulting in the formation of caustics, which are essentially discontinuities in the three-velocity field. In the context of entropy creation through violent relaxation, these discontinuities can be likened to non-collisional accretion shocks (see Parks et al. 2012, 2017 for non-collisional shocks in plasmas) and are associated with a microscopic stress tensor that is also discontinuous across the shocks (i.e. caustics). The nature of such discontinuities within the framework of general relativity may initially appear unclear. We propose characterising these discontinuities by utilising exact solutions of Einstein’s field equations in spherical symmetry. These solutions are locally comoving with the time-like geodesics and are known as generalised Lemaître-Tolman-Bondi solutions, when accounting for a microscopic stress tensor and virialisation (Lasky & Lun 2006).

The metric for such solutions with a fluid characterised by an energy density ρ and pressure σ is given by $$d{s}^2=-N(r,t{)}^{2}d{t}^2+\frac{1}{1+2E(r,t)}{R}^{\mathrm{\prime }2}d{r}^2+R(r,t{)}^{2}d{\mathrm{\Omega }}^2,$$(29)

with R′ the derivative with respect to r and with a unit choice such that c = 1. N(r, t) is the lapse function and E(r, t) the local curvature. The following Hamiltonian constraint equation can be derived from the Arnowitt–Deser–Misner (ADM) formalism (Lasky & Lun 2006) $$\frac{1}{2}{u}^2=(\frac{GM}{R}+E)\text{, }$$(30)

with $u\equiv \dot{R}/N,\dot{R}$ the derivative with respect to t and $M(r,t)=4\pi {\int }_0^{R(r,t)}\rho {\tilde{r}}^{2}d\tilde{r}$, with the following evolution equations: $$\begin{array}{rl}& \frac{\dot{E}}{N}=-\frac{1+2E}{\rho +\sigma }\frac{1}{R^{\mathrm{\prime }}}\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }r}u,\\ & \frac{\dot{\rho }}{N}=-N(\rho +\sigma )\frac{1}{R^2{R}^{\mathrm{\prime }}}\frac{\mathrm{\partial }\left(R^{2}u\right)}{\mathrm{\partial }r},\\ & \frac{\dot{u}}{N}=-\frac{GM}{R^2}-4\pi G\sigma R-\frac{1}{R^{\mathrm{\prime }}}\frac{1+2E}{\rho +\sigma }\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }r}.\end{array}$$(31)

The Euler equation provides the following relation between the lapse N and the pressure N′/N = −σ′/(ρ + σ). It is evident at this stage that the local curvature E(r, t) undergoes an evolution reminiscent of the behaviour of kinetic energy in classical Newtonian hydrodynamics. This similarity suggests that the curvature can be linked to microscopic dissipation through the presence of discontinuities. To establish this connection, we proceeded under the assumption that the fluid can be characterised by an internal microscopic energy, denoted as e, such that the total energy density ρ can be expressed as ρ = ρ_m + ρ_me, with ρ_m representing the mass density. It is important to emphasise that our analysis holds validity for both collisional and non-collisional systems. The microscopic internal energy can be defined for a non-collisional system as the microscopic stress tensor in the context of the Jeans equations. In the non-collisional limit, we do not have simple closure relations that establish a connection between the internal energy and the stress tensor with other variables. However, it is worth noting that such a relation is not a requisite element within our analysis. In the non-relativistic limit E, e ≪ 1, σ ≪ ρ. In this limit, and defining τ = 1/ρ_m, $$\begin{array}{rl}\dot{E}& =-\frac{\tau }{R^{\mathrm{\prime }}}\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }r}u,\\ \dot{\tau }& =\frac{\tau }{R^{\mathrm{\prime }}}\frac{1}{R^2}\frac{\mathrm{\partial }\left(R^{2}u\right)}{\mathrm{\partial }r},\\ \dot{e}& =-\frac{1}{{\rho }_m{R}^{\mathrm{\prime }}}\frac{1}{R^2}\frac{\mathrm{\partial }\left(R^{2}u\right)}{\mathrm{\partial }r}\sigma ,\\ \dot{u}& =-\frac{GM}{R^2}-\frac{\tau }{R^{\mathrm{\prime }}}\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }r}.\end{array}$$(32)

By employing the transformation dm = ρ_mR²dr/R′, where m represents the so-called ‘mass’ variable, one can readily identify the hydrodynamics equations in Lagrangian coordinates (see p. 8 of Godlewski & Raviart 1996 or p. 16 of Després 2017) with an additional equation governing the evolution of curvature. From this system, we derived the following conservative equation that establishes a coupling between curvature and internal energy: $$\dot{E}+\dot{e}=-\frac{\mathrm{\partial }}{\mathrm{\partial }m}\left(\sigma u{R}^2\right)\text{. }$$(33)

We assumed here that a discontinuity is located at r_i with R_i(t) = R(r_i, t), with the possibility for u and E to jump. By using standard Rankine-Hugoniot relations, $${\dot{m}}_i(E_l+e_l-E_r-e_r)=({\sigma }_l{u}_l-{\sigma }_r{u}_r)R_i^2,$$(34)

with ${\dot{m}}_i(t)$, the mass flow through the discontinuity. It is apparent from this relation that the behaviour of curvature bears a resemblance to that of kinetic energy. The second law of thermodynamics, encompassing the dissipation of large-scale motions and the concomitant generation of entropy, thus dictates that the evolution of curvature follows the condition ΔE < 0. By applying the Hamiltonian constraint and using the continuity of GM/R across the discontinuity, we derived the following relationship: $$\frac{1}{2}({\dot{R}}_l^2-{\dot{R}}_r^2)=E_l-E_r.$$(35)

Under the assumption that the solution tends to become homogeneous on each side of the interface, we arrived at the following jump condition concerning the local volume expansion rate: $$\frac{1}{18}({\theta }_l^2-{\theta }_r^2)=\frac{1}{R_i^2}(E_l-E_r).$$(36)

It is worth highlighting that these co-moving solutions can find utility in a broader context not necessarily tied to gravitational dynamics. While we have primarily considered particles following time-like geodesics, one can also investigate co-moving solutions with non-geodesic motions in the limit as G approaches zero. In this scenario, the solution becomes a spherical shock discontinuity akin to the Sedov solution in classical hydrodynamics. Within this context, θ_l = ∇ · u_g remains approximately constant within the volume undergoing expansion, while θ_r = 0 outside of it. These shocks dissipate large-scale motions into microscopic motions, accompanied by entropy creation and a distinctive jump in the local curvature, the evolution of which conforms to ΔE < 0. Commencing with a homogeneous system characterised by zero curvature and transitioning to a locally comoving spacetime with time-like geodesics, we find that virialisation leads to the generation of local negative curvature. Although the connection between virialisation and negative curvature has been discussed in existing literature (Roukema et al. 2013; Roukema 2018), our contribution lies in linking the curvature production by virialisation to the second law of thermodynamics. This implies that the averaged space-time curvature of a co-moving solution that traces the motion of time-like geodesics, cannot be zero in the presence of irreversible processes.

Fig. 1 Averaged expansion rate during the evolution of an inhomogeneous universe. Top: reversible evolution with expansion and compression modes. Down: irreversible evolution when compression modes are dissipated by virialisation through accretion shocks.

5 Link with observations and quantitative estimates

DESI Collaboration (2024c,a) present the release of the first science data from the Dark Energy Spectroscopic Instrument (DESI) project, which includes data from commissioning and Survey Validation phases conducted between December 2020 and June 2021. The authors highlight that the DESI collaboration has successfully validated its survey design and observing strategy, ensuring the data’s quality and completeness and that future releases will facilitate detailed studies of the Universe’s large-scale structure, providing critical insights into dark energy and the fundamental physics governing the Universe’s evolution. Following the initial data release, the DESI collaboration has unveiled its first-year analysis results of baryon acoustic oscillations based on extensive observations of galaxies, quasars, and the Lyman-α forest (DESI Collaboration 2024b). While the DESI data alone align with the Λ cold dark matter (CDM) model, the results deviate when generalised to the ω₀ω_aCDM model. The DESI data combined with cosmic microwave background (CMB) and supernova data suggest a rejection of the ΛCDM model in favour of the ω₀ω_aCDM model with significant confidence levels between 2 and 4σ (DESI Collaboration 2024b). The combined data indicate a preference for ω₀ > −1 and ω_a < 0, where ω(a) = ω₀ + ω_a(1 − a) represents the equation-of-state parameter of dark energy with a(t) the scale factor of the FLRW cosmology. If confirmed, these findings could have profound implications for our understanding of the Universe’s dynamics. In that context, Tada & Terada (2024) suggest that the nature of dark energy might be explained by a quintessential scalar field, although this field could potentially become a phantom in the past. They conclude that the DESI data challenge the ΛCDM model and supports a dynamic dark energy scenario, thereby contributing to our understanding of cosmic acceleration and the potential need for new physics beyond the standard cosmological model.

However, due to the covariance of general relativity, the current approach in the standard cosmological model can be reformulated as follows: the choice of an expanding background coordinate system should be arbitrary and introducing a cosmological constant to drive the dynamics of the coordinate system has no physical meaning other than selecting an expanding coordinate system with a_eff(t), H_eff(t), and three velocities u_eff,g for time-like geodesics such that ⟨θ⟩ ≈ 3H_eff and ⟨∇ · u_eff,g⟩ ≈ 0. This choice compensates for the averaged expansion rate along time-like geodesics in the coordinate system. If one prefers to use the FLRW solution of a fictive homogeneous universe without a cosmological constant, where the Hubble expansion rate is H(t) (with $H^2=8\pi G\overline{\rho }/3$), and the three velocities for time-like geodesics are u_g, then ⟨θ⟩ = 3H + ⟨∇ · u_g⟩ and ⟨∇ · u_g⟩ ≠ 0 as soon as virialisation and entropy production occur. We highlight that virialisation during structure formation is likely to replace dark energy because it precisely mimics the counter-intuitive anti-gravity effect needed for dark energy. The dissipation of large-scale compression modes by virialisation stabilises large-scale gravitational collapse (anti-gravity) while remaining a purely gravitational effect, as virialised dark matter continues to be in freefall.

In this effective coordinate system, we can use the second principle of thermodynamics to infer an effective evolution equation that takes into account entropy production, $$\begin{array}{rl}\frac{d(\log V_{\mathrm{e}\mathrm{f}\mathrm{f}})}{dt}& =\frac{d(\log V_{\mathrm{r}\mathrm{e}\mathrm{v}})}{dt}+\frac{1}{Z{k}_b}\frac{ds}{dt},\\ H_{\mathrm{e}\mathrm{f}\mathrm{f}}(t)& =H(t)+\frac{1}{3Z{k}_b}\frac{ds}{dt},\end{array}$$(37)

which can be rewritten as $$\begin{array}{rl}{H}_{\mathrm{e}\mathrm{f}\mathrm{f}}(t{)}^2& =\frac{8\pi G\overline{\rho }}{3}+\frac{{\mathrm{\Lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}}(t)c^2}{3},\\ {\mathrm{\Lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}}(t)& =\frac{2H}{c^{2}Z{k}_b}\frac{ds}{dt}+\frac{1}{3c^2{Z}^2{k}_b^2}{\left(\frac{ds}{dt}\right)}^2.\end{array}$$(38)

In this context, the effective cosmological constant appears to vary with time, as suggested by the early DESI data release. It is important to note, however, that this does not imply the existence of new physics; rather, the effective cosmological constant serves as a proxy to account for the impact of entropy production on the averaged expansion rate along time-like geodesics. Entropy production is expected to begin during structure formation at the first shell crossing, peak during the virialisation of structures, and gradually decline to the present day. This behaviour could simulate the need for ω(a) < −1 in the past and ω(a) > −1 in the present. As noted in Tada & Terada (2024), The increase in ω(a) is correlated with the decrease in the Hubble constant, which is the opposite direction to solve the Hubble tension. The Hubble Tension refers to the significant discrepancy between two different methods of measuring the Hubble constant. The first method involves local measurements, such as observing the distances and redshifts of nearby galaxies. The second method involves the CMB and large-scale structure data, which rely on ΛCDM model and measurements from the early Universe (Kamionkowski & Riess 2023). The local measurements consistently yield a higher value for the Hubble constant (around 73 km/s/Mpc) compared to the CMB-based measurements (around 67 km/s/Mpc). Although the decline phase of entropy production after the virialisation of structures cannot explain the Hubble tension, the initial increase phase after the first shell crossing could play an important role in resolving the current tension. However, precise estimates of the evolution of Λ_eff(t) require dedicated numerical simulations to accurately constrain the evolution of entropy production during structure formation and subsequent evolution.

We can, however, use existing numerical simulations to provide a more quantitative order-of-magnitude estimate. We considered the illustrative example depicted in the top panel of Fig. 2. In this scenario, we assumed that within the bound and virialised large-scale structures with a characteristic size of d_lss, the local volume expansion rate, denoted as θ, is approximately negligible (θ ≈ 0). Conversely, in the expanding cosmic voids with a typical size of d_void, θ is substantially higher, with θ ≈ 6u_pre–shock/d_void ≫ 3H. The bottom panel of Fig. 2 shows an equivalent apparent velocity profile by integrating the gradients along the line of sight. Both profiles in Fig. 2 provide the same integrated redshift using Eq. (19) and illustrate how peculiar velocities resulting from structure formation can mimic a Hubble flow. Without discontinuities in Fig. 2, the velocity field divergence in the structures and in the voids would exactly compensates for one another such that ⟨∇ · u_g⟩ = 0. However, this is not the case with discontinuities: in the extreme case for which ∇ · u_g is strictly positive everywhere (i.e. with expanding motions in the lss regions in Fig. 2), ⟨∇ · u_g⟩ is clearly strictly positive on average. In a realistic configuration for the large-scale structures with θ ≈ 0, the dominant contributor to the redshift arises from the expansion occurring within the voids, accounting for a fraction of d_void/(d_lss + d_void). This expansion in the voids cannot be counterbalanced by compression within the structures, as it has dissipated due to violent relaxation and virialisation processes. Assuming typical structure sizes of approximately d_lss ≈ 3 Mpc and void sizes around d_void ≈ 30 Mpc, the average expansion rate, represented by ⟨θ⟩/3, attains values of the order of 70 km/s/Mpc with pre-shock velocities of roughly u_pre–shock ≈ 1200 km/s. This estimation falls within the same order of magnitude as the infall velocities observed in N-body numerical simulations (Zu & Weinberg 2013).

The forthcoming generation of high-resolution cosmological simulations on exascale supercomputers will be essential in advancing our understanding of entropy production through caustics (i.e. velocity discontinuities or non-collisional shocks). These simulations will enable more precise estimates of the averaged local expansion rate along time-like geodesics. Such estimates are crucial for interpreting dynamic dark energy scenarios, where dark energy could be replaced by entropy production during structure formation. These insights will be essential for analysing future data releases from the DESI collaboration and for interpreting observations from the Euclid satellite (Amendola et al. 2018). By combining Euclid’s weak lensing data, redshift surveys (including baryon acoustic oscillations, redshift distortions, and the full power spectrum shape), and CMB data from Planck, we can constrain dynamic dark energy scenarios more effectively.

Fig. 2 Impact of discontinuities on a velocity profile. Top: illustration of a discontinuous velocity profile that impacts the redshift-distance relation. Bottom: equivalent apparent velocity profile derived via the integration of the gradients along the line of sight.

6 Conclusions

This paper highlights the intricate relationship between the cosmological redshift and entropy production due to structure formation. We emphasise the need for the covariant calculation of the cosmological redshift, which accounts for peculiar motions with a local expansion rate along time-like geodesics.

Our analysis suggests that virialisation and the subsequent entropy production during structure formation play a crucial role in this dynamic. The time-varying nature of the effective cosmological constant, as indicated by recent DESI data, may be a manifestation of these processes rather than new physics. This perspective necessitates a reassessment of dark energy models, and we propose that what is currently attributed to dark energy might instead be the result of entropy production during structure formation.

Future high-resolution cosmological simulations on exascale supercomputers will be vital for refining our understanding of these phenomena. These simulations will help produce more accurate estimates of the averaged local expansion rate along time-like geodesics, providing essential insights for interpreting dynamic dark energy scenarios and forthcoming observational data from DESI and the Euclid satellite. This research opens new avenues for understanding the fundamental forces shaping our Universe, emphasising the need for continuous advancements in computational and observational cosmology.

Acknowledgements

We thank Quentin Vigneron, Thomas Guillet, and Pascal Wang for all the discussions that lead to this paper. We also thank David Elbaz and Matthias González for helpful comments on this manuscript.

References

All Figures

	Fig. 1 Averaged expansion rate during the evolution of an inhomogeneous universe. Top: reversible evolution with expansion and compression modes. Down: irreversible evolution when compression modes are dissipated by virialisation through accretion shocks.
In the text

	Fig. 2 Impact of discontinuities on a velocity profile. Top: illustration of a discontinuous velocity profile that impacts the redshift-distance relation. Bottom: equivalent apparent velocity profile derived via the integration of the gradients along the line of sight.
In the text