Can rotation solve the Hubble Puzzle?

ABSTRACT

The discrepancy between low and high redshift Hubble constant |$H_0$| measurements is the highest significance tension within the concordance Lambda cold dark matter paradigm. If not due to unknown systematics, the Hubble Puzzle suggests a lack of understanding of the universe’s expansion history despite the otherwise spectacular success of the theory. We show that a Gödel inspired slowly rotating dark-fluid variant of the concordance model resolves this tension with an angular velocity today |$\omega _0 \simeq 2\times 10^{-3}$| Gyr⁻¹. Curiously, this is close to the maximal rotation, avoiding closed time-like loops with a tangential velocity less than the speed of light at the horizon.

distance scale, large-scale structure of Universe, cosmology: observations

1 INTRODUCTION

The Hubble tension, the inconsistency of the late and early time measurements of the universe’s expansion rate, emerges as the most significant chink in the otherwise shiny armour of the concordance Lambda cold dark matter (⁠|$\Lambda$|CDM) model (see e.g. for reviews by Verde, Treu & Riess 2019; Di Valentino et al. 2021a; Kamionkowski & Riess 2023). The discrepancy has been established in a wide range of data sets and reached a |$5\sigma$| significance between cepheid-calibrated local supernovae and cosmic microwave background (CMB) measurements (for counterpoint and calibration uncertainties see Freedman et al. 2024).

The CMB constraints at recombination are indirect: they assume an expansion history governed by the |$\Lambda$|CDM model. The latest analyses of Planck CMB maps imply a Hubble constant |$H_{\mathrm{CMB}} = 67.4 \pm 0.5 \: \mathrm{km} \mathrm{s}^{-1} \mathrm{Mpc}^{-1}$| (Aghanim et al. 2020).

Type Ia supernovae directly constrain the late-time (local) expansion rate. In a definitive study of Riess et al. (2022) used the Hubble Space Telescope (HST) to observe Cepheid variables in the host galaxies of 42 Type Ia supernovae (SNe Ia) crucial for calibrating the local Hubble constant (⁠|$H_{\mathrm{SNe}}$|⁠). They utilized all suitable SNe Ia discovered at redshift |$z \le 0.01$| over the past four decades, significantly expanding the sample size with observations from over 1000 HST orbits. They performed geometric calibration of Cepheids using Gaia EDR3 parallaxes, masers in NGC 4258, and detached eclipsing binaries in the Large Magellanic Cloud. Their baseline result is |$H_{\mathrm{SNe}} = 73.04 \pm 1.07 \:\mathrm{km} \mathrm{s}^{-1} \mathrm{Mpc}^{-1}$|⁠, with systematic uncertainties, closely aligned with the median of various analysis variants. Notably, they found a significant 5|$\sigma$| discrepancy with the Planck CMB analysis.

A burst of activity resulted, including alternative (tip of the red giant branch) calibrations by Freedman et al. (2019), and extensions or modifications of |$\Lambda$|CDM by Di Valentino, Melchiorri & Silk (2020a), such as massive neutrino or weakly interacting massive particles (WIMP) models of Pan & Knox (2015), dark photon of Aboubrahim, Klasen & Nath (2022), and an extended dark sector by Di Valentino et al. (2020b). Next, we propose rotating space–time as a novel solution. Gödel Gödel (1947) introduced a rotating universe followed by Heckmann & Schücking (1955, 1956a, b) and Heckmann (1961) later Silk (1966) and Hawking (1969). Visualization of the Gödel’s universe is made by Buser, Kajari & Schleich (2013). While anisotropies in a variety of Bianchi models with large vector perturbations corresponding to rotation are tightly constrained from Planck CMB data by Saadeh et al. (2016), generalizations of the Gödel model by Obukhov (2000) with global rotation are still viable and free of the pathologies of the original. This paper considers a Newtonian approximation of these models in the context of the Hubble anomalies.

All objects within our universe rotate, including planets, stars, solar systems, galaxies, and galaxy clusters. Moreover, black holes, spherically symmetric objects with horizons, display near maximal rotation as presented by Daly (2019). The idea that everything revolves (⁠|$\pi \alpha \nu \tau \alpha$|  |$\kappa \upsilon \kappa \lambda o \upsilon \tau \alpha \iota$|⁠)¹ naturally extends to the whole universe, as hinted by recent claims of anisotropic Hubble expansion in X-ray observations by Migkas et al. (2021). Furthermore, a plausible syllogism is that the universe has near-maximal rotation, motivated by cosmologies where the universe is the interior of a black hole (Pathria 1972). There are many proposed solutions to the Hubble Puzzle (e.g. Di Valentino et al. 2021b) and any modification of the standard model expansion and growth history has to consider the entire concordance model (e.g. Knox & Millea 2020). Nevertheless, as we show later, the average rotation effect has a similar functional form to that of dark photons (Fabbrichesi, Gabrielli & Lanfranchi 2021; Aboubrahim et al. 2022), one of the promising contenders (Cyr-Racine, Ge & Knox 2022) for solving the Hubble Puzzle. Therefore, exploring how a rotating model would affect the Hubble constant is worthwhile. In the next section, we outline our methodology, Section 3 presents the results, while the last section contains our conclusions.

2 METHODOLOGY

The expansion history of Newtonian cosmological simulations is in precise agreement with Friedmann models (e.g. Rácz et al. 2018). Thus, we expect that the classical framework is sufficient for an initial estimate of rotational effects on the Hubble constant; we leave general relativistic considerations for future work.

Describing the evolution of the Hubble parameter in Newtonian non-rotating and rotating universe models is still challenging. The Sedov–von Neumann–Taylor blast wave models inspired us to construct a non-relativistic dark fluid model. We apply a non-linear partial differential equation system describing a non-viscous, non-relativistic, and self-gravitating fluid with zero thermal conductivity (Euler–Poisson system) and solve it with a time-dependent Sedov-type self-similar ansatz. This analytic approach incorporates various scaling mechanisms and describes different time decay scenarios of Taylor (1950). The resulting dynamical model is consistent with direct solutions from the Friedmann equations by Szigeti, Barna & Barnaföldi (2023). We generalize our method of intermediate asymptotic analysis of the hydrodynamical description for the rotating dark-fluid universe to investigate the effect of rotation on the Hubble-constant anomaly. Our partial differential equations read as follows:

$$\begin{eqnarray} \partial _t \rho + \mathrm{div} (\rho \boldsymbol{u}) &= 0 \ , \end{eqnarray}$$

(1a)

$$\begin{eqnarray} \partial _t (\rho \boldsymbol{u}) + \mathrm{div} (\rho \boldsymbol{u} \otimes \boldsymbol{u}) &= - \nabla P(\rho)-\rho \nabla \Phi + \rho {\mathit{ g}^{\ast }} \ , \end{eqnarray}$$

(1b)

$$\begin{eqnarray} \nabla ^2 \Phi &= 4 \pi G \rho , \end{eqnarray}$$

(1c)

where |$\rho ,{{\bf \boldsymbol{ u}}}, P, \Phi , {\mathit{ g}}$| are the fluid density, the fluid velocity vector, the pressure, the gravitation potential and the external force, respectively. This system has been investigated previously by Deng, Xiang & Yang (2003) and Wong, Yeung & Yuen (2020). Goldreich & Weber (1980) studied the homologously collapsing stellar cores with an adiabatic exponent. Later Yuen (2009) gave analytically periodic solutions to the 3-dimensional Euler–Poisson equations of gaseous stars with negative cosmological constants, commonly used to describe a dark-fluid system. In this model, dark matter and dark energy are two different aspects of the same substance, the ‘dark fluid’ Farnes (2018). As illustrated in Fig. 1, the self-similar solution using the specific dark fluid equation of state (EOS) yields results consistent with |$\Lambda$|CDM within the relevant time range.

Figure 1.

Time evolution (log–log scale) of the Hubble parameter for non-rotating (analytical) and rotating (numerical) solutions at different |$\omega _0$| rotation parameter values as of today. Small figures show the evolution (normal scale) at the decoupling period |$t_{\rm CMB}$| and today |$t_0 = t_{\rm today}$|⁠.

• Spherical symmetry: Initially, we assumed an ideal fluid with spherical symmetry. Therefore, the multidimensional partial differential equation system reduces to a one-dimensional, radius-dependent ordinary system. We assume a linear EoS by Horedt (2004), |$P (\rho) = w \rho$| with the effective w asymptotically tending to −1 at |$t_{\infty }$|⁠. The resulting dark-fluid models describe a mixture of dark matter and dark energy in a non-rotating, expanding universe (Szigeti et al. 2023), as long as we neglect dark matter fluctuations (see extensive studies by Guo & Zhang 2007). The Euler–Poisson equation in the spherical limit is the following:

  $$\begin{eqnarray} \partial _t \rho + (\partial _r \rho) \boldsymbol{ u} + (\partial _r \boldsymbol{ u}) \rho + \dfrac{2\boldsymbol{ u}\rho }{r} &= 0 \ , \end{eqnarray}$$

  (2)

   

  $$\begin{eqnarray} \partial _t \boldsymbol{ u} + (u \partial _r) \boldsymbol{ u} &= - \dfrac{1}{\rho } \partial _r P-\partial _r \Phi (r) + g^{\ast } \ , \end{eqnarray}$$

  (3)

   

  $$\begin{eqnarray} \dfrac{1}{r^2} \dfrac{\mathrm{ d}}{\mathrm{ d}r}\left(r^2 \partial _r \Phi \right) & = 4 \pi \rho \ , \end{eqnarray}$$

  (4)

  where the |$u(r,t)$| radial flow velocity. Equations (1a)–(1c) are reduced to equations (2)–(4) due to the high similarity of the system as demonstrated by Sedov (1959).

• Cylindrical symmetry: If the rotation becomes significant enough, assuming spherical symmetry is no longer adequate. We assumed that the system is fully symmetric in the z direction. Thus, we extend our previous analyses of spherical flows to the Euler–Poisson equation in a cylindrical coordinate system.

  $$\begin{eqnarray} \partial _t \rho + (\partial _r \rho) \boldsymbol{ u} + (\partial _r \boldsymbol{ u}) \rho + \dfrac{\boldsymbol{ u}\rho }{r} &= 0 \ , \end{eqnarray}$$

  (5)

   

  $$\begin{eqnarray} \partial _t \boldsymbol{ u} + (\boldsymbol{ u} \partial _r) \boldsymbol{ u} &= - \dfrac{1}{\rho } \partial _r P-\partial _r \Phi (r) + g^{\ast } \ , \end{eqnarray}$$

  (6)

   

  $$\begin{eqnarray} \dfrac{1}{r} \dfrac{\mathrm{ d}}{\mathrm{ d}r}\left(r \partial _r \Phi \right) & = 4 \pi \rho \ . \end{eqnarray}$$

  (7)

  Holden et al. (2009) studied self-similar solutions for the infinite cylindrical collapse. Self-similar exponents similar to the spherical case exist orthogonal to the z-axis. The |$\boldsymbol{ u}(r,t)$| again has a similar meaning as defined in the spherical symmetric case.

We solved both sets of the equations by using the Sedov–Taylor ansatz for the velocity field |$\boldsymbol{ u}(r,t)$|⁠, the density |$\rho (r,t)$|⁠, and the gravitational potential density field |$\Phi (r,t)$| in both system. For spherical symmetry we applied the |$\boldsymbol{ u}(r,t) = t^{-\alpha } f(\eta)$|⁠, |$\rho (r,t) = t^{-\gamma } g (\eta)$|⁠, and |$\Phi (r,t) = t^{-\delta } h(\eta)$| notation, where |$f(\eta)$|⁠, |$g(\eta)$|⁠, and |$h(\eta)$| are the shape functions of the reduced ordinary differential equation system with the reduced variable, |$\eta = r/t^{\beta }$|⁠. The analogous shape functions for cylindrical symmetry depend on the z coordinate: |$u(z,t) = t^{-\alpha } k(\xi)$|⁠, |$\rho (z,t) = t^{-\gamma } l (\xi)$|⁠, and |$\Phi (z,t) = t^{-\gamma } h(\xi)$|⁠, where |$\xi = z/t^{\beta }$|⁠.

The real parameters |$\alpha , \beta , \gamma$|⁠, and |$\delta$| are the self-similar exponents responsible for the solution’s temporal decay and spreading. The numerical values of the exponents: |$\alpha = 0$|⁠, |$\beta = 1$|⁠, |$\gamma = 2$|⁠, and |$\delta = 0$|⁠. We assumed similarly to the spherical case for the cylindrical case that the dynamical variables depend only on the |$\eta = r/t^{\beta }$| variables. We performed all calculations in physical (non-expanding) rotating coordinates. The scale factor is the solution of the differential equation of |$\dot{a}(t) = \boldsymbol{ u}(r(t),t)$| (for details, see technical Szigeti et al. 2023), i.e. the approximately uniform streaming of the particles in physical coordinates represents the expansion of the universe. We restricted to expanding solutions of equations (2)–(7) with |$r \propto a(t)$|⁠. Consequently, the Hubble parameter is expressed as usual |$H(t) = \dot{a}/a$|⁠.

We solved the two systems under equivalent initial conditions in |$\eta$|⁠. We have seen that the two systems exhibit analogous behaviour in regions where the intermediate-asymptotic behaviour yields a valid solution. Thus, we will use the spherically symmetric limit in our following analyses (see in Fig. 1). We assumed spherical symmetry for our self-similar solutions of the Euler–Poisson system. This is justifiable far from the boundary, and it is numerically more stable Yuen (2009) and Peng & Lien (2012).

The same mathematical framework describes the Hubble parameter for the non-rotating (⁠|${g^{\ast }}=0$|⁠) and the rotating model (⁠|${g^{\ast }} \ne 0$|⁠) by setting initial conditions for the density, the velocity and the gravitational field: |$\boldsymbol{ u}(\eta _{\textrm {IC}}) = 0.5$| and |$\rho (\eta _{\textrm {IC}}) = 0.01$|  |$\mathrm{l}_{ \textrm { }P}^{-1}/\mathrm{m}^{3}$| in geometrized units (⁠|$G = c = 1)$| at |$t_{\textrm {CMB}} = t_{\textrm {IC}} = 380$| kyr to determine the shape functions.² Similarly, in case of cylindrical symmetric system the |$k(\xi _{ \textrm {IC}}) = 0.5$| and |$l(\xi) = 0.01$|⁠. These initial conditions are consistent with al. (2020), and, as we have shown, they reproduce the no-rotation Friedman solutions. These initial conditions are consistent with the Planck initial conditions at the decoupling period. We can add an effective rotational term to equation (1b),

$$\begin{eqnarray} {g^{\ast }} = 2r \omega ^2(t) \sin (\theta), \end{eqnarray}$$

(8)

where the angular velocity is |$\omega (t) = (\omega _0/t_0)/t$| and |$\theta$| is the polar angle. The effective centrifugal force is given in a non-inertial rotating frame. The Coriolis force vanishes since the velocity from expansion is always perpendicular to the rotation axis. Slow rotation can still be consistent with present observations. A slight global rotation still preserves a uniform CMB (Obukhov 2000; McEwen et al. 2013; Saadeh et al. 2016). Soon, such a rotation might be constrained by comparing the local inertial frame with that of quasars in Szapudi (2021).

3 RESULTS AND DISCUSSION

We numerically calculate the evolution equations following Szigeti et al. (2023), transforming the equations into the co-moving frame and applying the constraints detailed in Section (5). Fig. 1 shows the time dependence of the Hubble parameter for various angular frequencies from recombination until the present. Different initial rotations result in different |$H_0$| values today, but all solutions converge to zero at the asymptotic limit, |$t \rightarrow \infty$|⁠. The evolution of the Hubble parameter with initial rotation in |$\omega (t) \rightarrow 0$| limit approaches the non-rotating model. However, the limit is extrapolated due to numerical instabilities for extremely small |$\omega _0$| values today. As a test, the black solid line displays the standard |$\Lambda$|CDM result,³ in perfect agreement with our formalism (the orange solid line labelled non-rotating). The non-rotating self-similar solution is consistent between |$t_{\rm CMB}$| and |$t_{\rm today}$|⁠.⁴ Fig. 2 displays the Hubble constant (⁠|$H_0$|⁠) for the non-rotating (analytical) and various slowly rotating (numerical) cases as a function of the rotational parameter |$\omega _0$|⁠. The solid grey curve represents interpolation to the numerical calculations (markers), while the dashed curve is an extrapolation for |$\omega _0 \rightarrow 0$|⁠. Numerical extrapolation for the Hubble constant, |$H_0$| with |$\omega _0 = 0.002^{+0.001}_{-0.0009}$| Gyr|$^{-1}$| predicts a value today comparable to the measured by |$H_{\mathrm{SNe}}$|⁠. The present day |$\omega _0$| rotation corresponds to an initial, |$\omega (t_{\rm CMB})= 3.54^{+1.3}_{-1.2}$| Myr|$^{-1}$|⁠, where |$H_{\rm CMB}$| is measured at |$t_{\mathrm{CMB}}=380$| kyr. In Fig. 3, we illustrate the rotating universe. Its angular rotation parameter, |${\boldsymbol{\omega }}(t)$|⁠, is approximately

$$\begin{eqnarray} |{\boldsymbol{\omega }}(t)| = \omega _0 a^{-2}(t) \end{eqnarray}$$

(9)

from angular momentum conservation during matter domination, consistently with equation (8). Next, we estimate the maximal rotation of a dark-matter-filled universe and compare it with |$\omega _0 \simeq 0.002$| solving the Hubble tension. We require that the speeds remain below the speed of light within the observable horizon, hence |$\omega \lesssim H$|⁠, during the universe’s entire history. Taking |$H(a) \sim a^{-3/2}$|⁠, we limit |$\omega _0$| today as,

$$\begin{eqnarray} \omega _0 \lesssim H_0 a^{1/2}(t_{\mathrm{eq}}) \simeq 0.002\ \mathrm{Gyr}^{-1}, \end{eqnarray}$$

(10)

where |$t_{\mathrm{eq}}$| is the time of matter–radiation equality. Note that since at earlier times |$H(a)\simeq 1/a^2$| and |$\omega (a) \simeq 1/a$|⁠, the above condition is satisfactory for the entire evolution of the universe. Most remarkably, the allowed maximal rotation is approximately the same as the one required to solve the Hubble Puzzle. Our simplified argument neglected any late effects of Dark Energy on angular momentum, but there should be a reasonable estimate for our calculation. Our results are consistent with Heckmann & Schücking (1955, 1956a, b) and Heckmann (1961), despite the differences in techniques and their original motivation of removing the initial singularity at the big bang. The required a minimal rotation, |$\omega \simeq 0.03~\mathrm{Gyr}^{-1}$|⁠, is an order of magnitude larger than the maximal rotation avoiding closed time-like loops within the horizon.

Figure 2.

The predicted |$H_0$| values from non-rotating (analytical) and rotating (numerical) models evaluated at different |$\omega _0$| values today. The solid curve interpolates the calculated values (markers). The continuing dashed curve extrapolates the |$\omega _0 \rightarrow 0$| case. The lower dashed (⁠|$H_{\mathrm{CMB}}$|⁠) (Aghanim et al. 2020) and the upper dashed (⁠|$H_{\mathrm{SNe}}$|⁠) (Riess et al. 2022) lines correspond to measurements with |$2\sigma$| uncertainty ranges. The shading approximates the prohibited region exceeding maximal rotation.

Figure 3.

Schematic view of the rotating (non-expanding) spherical and cylindrical physical coordinate systems, where the outflow of the particles, |$u(r(t),t)$| determines the expansion rate through |$\dot{a}(t)= u(r(t),t)$|⁠; |$\omega$| refers to the angular velocity. Even though |$u(r(t),t)$| has a formal dependence on r, our solutions produce an outflow uniform enough that it is well described by a single expansion rate, |$a(t)$|⁠.

4 CONCLUSION

We analyse the time evolution of the Hubble parameter within the Euler–Poisson model with a self-similar time-dependent Sedov-type scaling for a linearized dark-fluid EOS. This model is consistent with a Newton–Friedmann cosmology when the angular momentum is zero and facilitates the analysis of cosmologies with slow rotation.

We found that an angular speed near the maximal rotation |$\omega _0 \lesssim 0.002\ \mathrm{Gyr}^{-1}$| today predicts a Hubble constant consistent with local measurements when starting from an expansion rate consistent with the CMB. Extrapolation to the initial rotation of the early universe gives the values of |$\omega (t_{\rm CMB}) \approx 3.54^{+1.3}_{-1.2}$| Myr|$^{-1}$| for the time of the origin of the CMB. These tantalizing initial results have the caveat that we only focused on the Hubble constant. Further investigations contrasting the rotating model against the entire intertwined network of the concordance model observations, confirmation and development of numerical models using rotating cosmological N-body simulations⁵, and extension for a general relativistic treatment are left for future work.

ACKNOWLEDGEMENTS

We acknowledge the financial support by the Hungarian National Research, Development and Innovation Office (NKFIH) under contract no. OTKA K135515, 2021–4.1.2-NEMZ_KI-2024-00031, 2024–1.2.5-TÉT-2024-00022, and Wigner Scientific Computing Laboratory (WSCLAB, the former Wigner GPU Laboratory) IS acknowledges NASA grants N4-ADAP24-0021 and 24-ADAP24-0074 and thanks the hospitality of the MTA-CSFK Lendület ‘Momentum’ Large-Scale Structure (LSS) Research Group at Konkoly Observatory supported by a Lendület excellence grant by the Hungarian Academy of Sciences (MTA).

DATA AVAILABILITY

No new data were generated or analysed in support of this research. The software code developed for this article will be available upon reasonable request.

Footnotes

1

Panta kykloutai paraphrasing the ancient Greek philosopher Heraclitus.

2

In SI: |${\sim} 10^{-15} \mathrm{kg}/\mathrm{m}^{3}$| and |$h(\eta _{\textrm {IC}}) = 10^{-12}$|⁠.

3

The curve is evaluated by using the matter-dominated scale factor |$a(t) = \left(3 H_{\rm CBM;0} t/2 \right)^{2/3}\Omega _\mathrm{ m}^{1/3}$|⁠, with the value of |$\Omega _\mathrm{ m}=0.3089$|⁠.

4

An appropriate choice of initial conditions for the scale factor in equation (25) in the work of Szigeti et al. (2023) can be reduced, due to |$u_2 \sim \mathcal {O}(10^{-4})$|⁠.

5

Pál et al. 2025 in preparation.

Author notes

MTA Guest Professor 2023