Introduction

Decades of progress in observational and theoretical cosmology have led to the consensus that our universe is well described by a flat Friedman–Robertson–Lemaitre metric and is currently comprised of around 5% baryons, 25% cold dark matter (CDM), and 70% dark energy in its simplest form—the cosmological constant Λ. Although this ΛCDM model fits many observations exquisitely well, its prediction for the present-day cosmic expansion rate, H0 = 67.36 ± 0.54 km/s/Mpc1, based on precise cosmic microwave background (CMB) radiation observations by the Planck satellite, do not compare well with direct measurements of the Hubble constant. In particular, the Supernovae H0 for the Equation of State (SH0ES) collaboration2, using Cepheid calibrated supernovae Type Ia, finds a much higher value of H0 = 73.5 ± 1.4 km/s/Mpc. This 4.2σ disagreement, known as the “Hubble tension”, has spurred much interest in modifications of the ΛCDM model capable of resolving it (cf.3 for a comprehensive list of references). Several other determinations of H0, using different methods, are also in some degree of tension with Planck, such as the Megamaser Cosmology Project4 finding 73.9 ± 3.0 km/s/Mpc or H0LiCOW5 finding \(73.{3}_{-1.8}^{+1.7}\) km/s/Mpc. It is worth noting that a somewhat lower value of 69.8 ± 2.5 km/s/Mpc was obtained using an alternative method for calibrating SNIa6.

Among the most precisely measured quantities in cosmology are the locations of the acoustic peaks in the CMB temperature and polarization anisotropy spectra. They determine the angular size of the sound horizon at recombination,

$${\theta }_{\star }\,\equiv\, \frac{{r}_{\star }}{D({z}_{\star })},$$
(1)

with an accuracy of 0.03%1. The sound horizon r is the comoving distance a sound wave could travel from the beginning of the universe to recombination, a standard ruler in any given model, and D(z) is the comoving distance from a present-day observer to the last scattering surface, i.e., to the epoch of recombination. D(z) is determined by the redshift-dependent expansion rate H(z) = h(z) × 100 km/s/Mpc which, in the flat ΛCDM model, depends only on two parameters (see Methods for details): Ωmh2 and h, where Ωm is the fractional matter energy density today and h = h(0) = H0/100 km/s/Mpc. Thus, given r and an estimate of Ωmh2, one can infer h from the measurement of θ. Using the Planck best fit values of Ωmh2 = 0.143 ± 0.001 and r = 144.44 ± 0.27 Mpc, obtained within the ΛCDM model1, yields a Hubble constant significantly lower than the more direct local measurements.

If the value of the Hubble constant was the one measured locally, i.e., h ≈ 0.735, it would yield a much larger value of θ unless something else in Eq. (1) was modified to preserve the observed CMB acoustic peak positions. There are two broad classes of models attempting to resolve this tension by introducing new physics. One introduces modifications at late times (i.e., lower redshifts), e.g., by introducing a dynamical dark energy or new interactions among the dark components that alter the Hubble expansion to make it approach a higher value today, while still preserving the integrated distance D in Eq. (1). In the second class of models, the new physics aims to reduce the numerator in Eq. (1), i.e., modify the sound horizon at recombination.

Late time modifications based on simple phenomenological parameterizations tend to fall short of fully resolving the tension7. This is largely because the baryon acoustic oscillation (BAO) and and supernovae (SN) data, probing the expansion in the 0 z 1 range, are generally consistent with a constant dark energy density. One can accommodate a higher value of H0 by making parameterizations more flexible, as e.g., in8,9, that allow for a non-monotonically evolving effective dark energy fluid. Such non-monotonicity tends to imply instabilities within the context of simple dark energy and modified gravity theories10 but can, in principle, be accommodated within the general Horndeski class of scalar-tensor theories11.

Early-time solutions aim to reduce r with essentially two possibilities: (i) a coincidental increase of the Hubble expansion around recombination or (ii) new physics that alters the rate of recombination. Proposals in class (i) include the presence of early dark energy12,13,14,15,16,17, extra radiation in either neutrinos18,19,20,21 or some other dark sector22,23,24,25,26,27, and dark energy–dark matter interactions28. Proposals in class (ii) include primordial magnetic fields29, non-standard recombination30, or varying fundamental constants31,32. In this work we show that any early-time solution which only changes r can never fully resolve the Hubble tension without being in significant tension with either the weak lensing (WL) surveys33,34 or BAO35 observations.

Results and discussion

The acoustic peaks, prominently seen in the CMB anisotropy spectra, are also seen as BAO peaks in the galaxy power spectra and carry the imprint of a slightly different, albeit intimately related, standard ruler—the sound horizon at the “cosmic drag” epoch (or the epoch of baryon decoupling), rd, when the photon drag on baryons becomes unimportant. As the latter takes place at a slightly lower redshift than recombination, we have rd ≈ 1.02r with the proportionality factor being essentially the same in all proposed modified recombination scenarios. More importantly for our discussion, the BAO feature corresponds to the angular size of the standard ruler at zz, i.e., in the range 0 z 2.5 accessible by galaxy redshift surveys. For the BAO feature measured using galaxy correlations in the transverse direction to the line of sight, the observable is

$${\theta }_{\perp }^{{\rm{BAO}}}({z}_{{\rm{obs}}})\,\equiv\, \frac{{r}_{{\rm{d}}}}{D({z}_{{\rm{obs}}})}\ ,$$
(2)

where zobs is the redshift at which a given BAO measurement is made. For simplicity, we do not discuss the line of sight and the “isotropic” BAO measurements36 here, but our arguments apply to them as well. It is well known that BAO measurements at multiple redshifts provide a constraint on rdh and Ωm.

In any particular model, r (and rd) is a derived quantity that depends on Ωmh2, the baryon density and other parameters. However, in this work, for the purpose of illustrating trends that are common to all models, we treat r as an independent parameter and assume that no new physics affects the evolution of the universe after recombination.

Without going into specific models, we now consider modifications of ΛCDM which decrease r, treating the latter as a free parameter and taking rd = 1.0184r. The relation between r and rd in different models that reduce the sound horizon is largely the same as the one in ΛCDM, hence we fix it at the Planck best fit ΛCDM value. For a given Ωmh2, Eq. (1) defines a line in the rdH0 plane, and since Eqs. (1) and (2) are the same in essence, a BAO measurement at each different redshift also defines a respective line in the rdH0 plane. However, the significant difference between z and zobs results in different slopes of the respective rd(h) lines (see Methods for details), as illustrated in Fig. 1. The latter shows the rd(h) lines from two different BAO observations, one at redshift z = 0.5 and another at z = 1.5, at Ωmh2 fixed to the Planck best fit ΛCDM value of 0.143, and the analogous lines defined by the CMB acoustic scale plotted for three values of Ωmh2: 0.143, 0.155, and 0.167. Both lines correspond to transverse BAO measurements. Slopes derived from the line of sight and isotropic BAO at the same redshift would be different, but the trend with increasing redshift is the same. The lines are derived from the central observational values and do not account for the uncertainties in \({\theta }_{\perp }^{{\rm{BAO}}}\) and θ (although the uncertainty in θ is so tiny that it would be difficult to see by eye on this plot). As anticipated, the slope of the rd(h) lines becomes steeper with increased redshift.

Fig. 1: A plot illustrating that achieving a full agreement between cosmic microwave background (CMB), baryonic acoustic oscillations (BAO) and SH0ES through a reduction of rd requires a higher value of Ωmh2.
figure 1

Shown are the lines of degeneracy between the sound horizon rd and the Hubble constant H0 defined by the CMB acoustic scale θ at three different values of Ωmh2: 0.143, 0.155, and 0.167. Also shown are the marginalized 68% and 95% CL bands derived from the combination of all current BAO data, and the ΛCDM based bounds from Planck. To demonstrate how the slope of the lines changes with redshift, we show two lines corresponding to the SDSS measurements of \({\theta }_{\perp }^{{\rm{BAO}}}\) at z = 0.51 and z = 1.554 at a fixed Ωmh2 = 0.143. The gray band shows the 68% and 95% CL determination of the Hubble constant by SH0ES.

Also shown in Fig. 1 are the marginalized 68% and 95% confidence levels (CL) derived from the combination of all presently available BAO observations in a recombination-model-independent way, namely, while treating rd as an independent parameter (see37 and Methods for details). The red contours show the ΛCDM based constraint from Planck, in good agreement with BAO at H0 ≈ 67 km/s/Mpc, but in tension with the SH0ES value shown with the gray band. In order to reconcile Planck with SH0ES solely by reducing rd, one would have to move along one of the CMB lines. Doing it along the line at Ωmh2 = 0.143 would quickly move the values of rd and H0 out of the purple band, creating a tension with BAO. Full consistency between the observed CMB peaks, BAO and the SH0ES Hubble constant could only be achieved at a higher value of Ωmh2 ≈ 0.167. However, unless one supplements the reduction in rd by yet another modification of the model, such high values of Ωmh2 would cause tension with galaxy WL surveys such as the Dark Energy Survey (DES)33 and the Kilo-Degree Survey (KiDS)34, which we illustrate next.

DES and KiDS derived strong constraints on the quantity \({S}_{8}\,\equiv\, {\sigma }_{8}{({{{\Omega }}}_{m}/0.3)}^{0.5}\), where σ8 is the matter clustering amplitude on the scale of 8 h−1 Mpc, as well as Ωm. The value of S8 depends on the amplitude and the spectral index of the spectrum of primordial fluctuations, which are well-determined by CMB and have similar best fit values in all modified recombination models. S8 also depends on the net growth of matter perturbations which increases with more matter, i.e., a larger Ωmh2.

The values of S8 and Ωm obtained by DES and KiDS are already in slight tension with the Planck best fit ΛCDM model, and the tension between KiDS and Planck is notably stronger than that between DES and Planck. Increasing the matter density aggravates this tension – a trend that can be seen in Fig. 2. The figure shows the 68% and 95% CL joint constraints on S8m by DES supplemented by the Pantheon SN sample38 (which helps by providing an independent constraint on Ωm), along with those by Planck within the ΛCDM model. The purple contours (Model 2) correspond to the model that can simultaneously fit BAO and CMB acoustic peaks at Ωmh2 = 0.155, i.e., the model defined by the overlap between the BAO band and the \({\theta }_{\star }^{(2)}\) (blue dashed) line in Fig. 1. The green contours (Model 3) are derived from the model with Ωmh2 = 0.167 corresponding to the overlap region between the \({\theta }_{\star }^{(3)}\) (green dotted) line and the BAO and SH0ES bands in Fig. 1 (see Methods for details). The figure shows that when attempting to find a full resolution of the Hubble tension, with CMB, BAO, and SH0ES in agreement with each other, one exacerbates the tension with DES and KiDS.

Fig. 2: The 68% and 95% confidence level bounds on S8 and Ωm.
figure 2

Shown are the constraints derived by fitting the ΛCDM model to a joint dataset of Dark Energy Survey (DES) and supernovae (SN) and to Planck, along with the contours for Model 2 and Model 3. Model 2 is defined by the simultaneous fit to baryonic acoustic oscillations (BAO) and cosmic microwave background (CMB) acoustic peaks at Ωmh2 = 0.155, i.e., the overlap between the BAO band and the \({\theta }_{\star }^{(2)}\) line in Fig. 1. Model 3 has Ωmh2 = 0.167 and corresponds to the overlap region between the \({\theta }_{\star }^{(3)}\) line and the BAO and SH0ES bands in Fig. 1.

We note that there is much more information in the CMB than just the positions of the acoustic peaks. It is generally not trivial to introduce new physics that reduces r and rd without also worsening the fit to other features of the temperature and polarization spectra39,40. Our argument is that, even if one managed to solve the Hubble tension by reducing r while maintaining a perfect fit to all CMB data, one would still necessarily run into problems with either the BAO or WL.

Surveying the abundant literature of the proposed early-time solutions to the Hubble tension, one finds that the above trends are always confirmed. Figure 3 shows the best fit values of rdh, H0, and S8 in models from Refs. 13,14,18,23,24,28,29,30,32. Note that there are other proposed early-time solutions to the Hubble tension. Figure 3 only shows the models for which explicit estimates of H0, Ωmh2, S8, and possibly rdh were provided. One can see that, except for the model represented by the red dot at the very right of the plot, corresponding to the strongly interacting neutrino model of18, solutions requiring low Ωmh2 are in tension with BAO, whereas solutions with higher Ωmh2 are in tension with DES and KiDS. This latter tension was previously observed and extensively discussed in the context of the early dark energy models41,42,43,44,45,46. As we have shown in this paper, it is part of a broader problem faced by all proposals aimed at reducing the Hubble tension in which the main change amounts to a reduction of rd.

Fig. 3: A compilation of values of Ωmh2, rdh, H0, and S8 predicted by some of the models aiming to relieve the Hubble tension by lowering the sound horizon.
figure 3

The best fit values of S8, H0, rdh (ac respectively), along with Ωmh2 (the horizontal axis), obtained within the models listed on the right. The horizontal bands show the 68% confidence level observational constraint on the corresponding parameter from different (types of) surveys. The sub-labels I and II in the list of models denote either different choices of model parameters within the same model, or constraints derived from different data combinations on the same model. The red square point with error bars represents the Planck best fit ΛCDM model1. With the exception of the red dot, corresponding to the model from18 with multiple modifications of ΛCDM fit to Planck temperature anisotropy data only, there is a consistent trend: models with low Ωmh2 either fail to achieve a sufficiently high H0 or are in tension with baryonic acoustic oscillations (BAO), and models with high values of Ωmh2 run into tension with the Dark Energy Survey (DES) or the Kilo-Degree Survey (KiDS).

In most of the models represented in Fig. 3, the effect of introducing new physics only amounts to a reduction in rd. We note that, in any specific model of a reduced rd, the best fit values of other cosmological parameters also change, which can affect the quality of the fit to various datasets. However, such changes, e.g., in the best fit value of the spectral index ns which affects S8, tend to be small for the models studied in the literature and have a minor impact compared to the effect of reducing rd, which is a pre-requisite for reconciling CMB with SH0ES. As we have argued, this will necessarily limit their ability to address the Hubble tension while staying consistent with the large scale structure data. Resolving the Hubble tension by new early-time physics without creating other observational tensions requires more than just a reduction of the sound horizon. This is exemplified by the interacting dark matter-dark radiation model25 and the neutrino model18 proposed as solutions. Here, extra tensions are avoided by supplementing the reduction in the sound horizon due to extra radiation by additional exotic physics: dark matter-dark radiation interactions in the first case and neutrino self-interactions and non-negligible neutrino masses in the second case. Consequently, with so many parameters, the posteriori probabilities for cosmological parameters are highly inflated over those for ΛCDM. It is not clear how theoretically appealing such scenarios are, and the model in18 seems to be disfavored by the CMB polarization data.

In conclusion, we have argued that any model which tries to reconcile the CMB inferred value of H0 with that measured by SH0ES by only reducing the sound horizon automatically runs into tension with either the BAO or the galaxy WL data. While we do not expect our findings to be surprising for the majority of the community, the novelty of our result is in isolating and clearly stating the essence of the problem—that the slopes of the rH0 degeneracy lines for BAO and CMB are vastly different, thus making it impossible to reconcile CMB with SH0ES by reducing r without violating BAO. We believe this very simple fact has not been stated before in this context in a model-independent way. With just a reduction of r, the highest value of the Hubble constants one can get, while remaining in a reasonable agreement with BAO and DES/KiDS, is around 70 km/s/Mpc. Thus, a full resolution of the Hubble tension will require either multiple modifications of the ΛCDM model or discovering systematic effects in one or more of the datasets.

Methods

The acoustic scale measurements from the CMB and BAO

The CMB temperature and polarization anisotropy spectra provide a very accurate measurement of the angular size of the sound horizon at recombination,

$${\theta }_{\star }\,=\,\frac{{r}_{\star }}{D({z}_{\star })}\ ,$$
(3)

where r is the sound horizon at recombination, or the comoving distance a sound wave could travel from the beginning of the universe to recombination, and D(z) is the comoving distance from a present-day observer to the last scattering surface, i.e., to the epoch of recombination. In a given model, r and D(z) can be determined from \({r}_{\star }\,=\,\mathop{\int}\nolimits_{{z}_{\star }}^{\infty }{c}_{s}(z){\rm{d}}z/H(z)\) and \(D({z}_{\star })\,=\,\mathop{\int}\nolimits_{0}^{{z}_{\star }}c\,{\rm{d}}z/H(z)\), where cs(z) is the sound speed of the photon–baryon fluid, H(z) is the redshift-dependent cosmological expansion rate and c is the speed of light. To complete the prescription, one also needs to determine z using a model of recombination.

The redshift dependence of the Hubble parameter in the ΛCDM model can be written as

$$h(z)\,=\,\sqrt{{{{\Omega }}}_{r}{h}^{2}{(1\,+\,z)}^{4}\,+\,{{{\Omega }}}_{m}{h}^{2}{(1\,+\,z)}^{3}\,+\,{{{\Omega }}}_{{{\Lambda }}}{h}^{2}}$$
(4)

where h(z) is simply H(z) in units of 100 km/s/Mpc, and h is the value at redshift z = 0. Here, Ωr, Ωm, and ΩΛ are the present-day density fractions of radiation, matter (baryons and CDM) and dark energy. From the precise measurement of the present-day CMB temperature T0 = 2.7255 K (however, also see47), and adopting the standard models of particle physics and cosmology, one knows the density of photons and neutrinos Ωrh2. Using the theoretically well motivated criticality condition on the sum of the fractional densities, i.e., Ωr + Ωm + ΩΛ = 1, one finds that h(z) is dependent only on two remaining quantities: Ωmh2 and h. The photon–baryon sound speed cs in Eq. (1) is determined by the ratio of the baryon and photon densities and is well-constrained by both Big Bang nucleosynthesis and the CMB. Fitting the ΛCDM model to CMB spectra also provides a tight constraint on Ωmh2, making it possible to measure h.

In alternative models, a smaller r is achieved by introducing new physics that reduces z through a modification of the recombination process or by modifying h(z) before and/or during recombination, or a combination of the two. In our analysis, we consider Eq. (3) while remaining agnostic about the particular model that determines the sound horizon. Namely, we treat r as an independent parameter. We assume, however, that after the recombination, the expansion of the universe is well described by Eq. (4), which is the case in many alternative models. Thus, our independent parameters are r, Ωmh2, and h, with the latter two determining D(z). The dependence of D(z) on the precise value of z is very weak, so that the differences in z in different models do not play a role.

The same acoustic scale is also imprinted in the distribution of baryons. There are three types of BAO observables corresponding to the three ways of extracting the acoustic scale from galaxy surveys36: using correlations in the direction perpendicular to the line of sight, using correlations in the direction parallel to the line of sight, and the angle-averaged or “isotropic” measurement. While our MCMC analysis includes all three types of the BAO data, for the purpose of our discussion it suffices to consider just the first type, which is the closest to CMB in its essence, but our conclusions apply to all three. Namely, we consider

$${\theta }_{\perp }^{{\rm{BAO}}}({z}_{{\rm{obs}}})\,\equiv\, \frac{{r}_{{\rm{d}}}}{D({z}_{{\rm{obs}}})}\ ,$$
(5)

where \({r}_{{\rm{d}}}\,=\,\mathop{\int}\nolimits_{{z}_{{\rm{d}}}}^{\infty }{c}_{s}(z){\rm{d}}z/H(z)\) is the sound horizon at the epoch of baryon decoupling, closely related to r, and zobs is the redshift at which a given BAO measurement is made. We adopt a fixed relation rd = 1.0184r that holds for the Planck best fit ΛCDM model and is largely unchanged in the alternative models.

As the distance integrals D(z) and D(zobs) in the denominators of Eqs. (3) and (5) are dominated by the matter density at low redshifts, one can safely neglect Ωrh2 and write

$${\theta }_{\star }\,=\,\frac{{r}_{\star }}{2998\ {\rm{Mpc}}}{\left(\mathop{\int}\nolimits_{0}^{{z}_{\star }}\frac{{\rm{d}}z}{{\omega }_{m}^{1/2}\sqrt{{(1\,+\,z)}^{3}\,+\,{h}^{2}/{\omega }_{m}\,-\,1}}\right)}^{-1},$$
(6)

where ωm = Ωmh2 and 2998 Mpc = c/100km/s/Mpc, and an analogous equation for BAO with the replacement \(({r}_{\star },{\theta }_{\star },{z}_{\star })\,\to\, ({r}_{{\rm{d}}},{\theta }_{\perp }^{{\rm{BAO}}},{z}_{{\rm{obs}}})\). For a given Ωmh2, Eq. (6) defines a line in the rdH0 plane. Similarly, a BAO measurement at each different redshift also defines a respective line in the rdH0 plane. Taking the derivative of r with respect to h one finds

$$\frac{\partial {r}_{\star }}{\partial h}\,=\,-\frac{h}{{\omega }_{m}}{\theta }_{\star }\mathop{\int}\nolimits_{0}^{{z}_{\star }}\frac{2998\ {\rm{Mpc}}\ {\rm{d}}z}{{\omega }_{m}^{1/2}{({(1\,+\,z)}^{3}\,+\,{h}^{2}/{\omega }_{m}\,-\,1)}^{3/2}}$$
(7)

and a completely analogous equation for BAO. It is important to realize that the derivative is very different for CMB and BAO due to the vast difference in redshifts at which the standard ruler is observed, z ≈ 1100 for CMB vs. zobs ~ 1 for BAO, resulting in different values of the integral in Eq. (7). This results in different slopes of the respective rd(h) lines. Note that the slopes of the rd(h) lines differ for the transverse, parallel and volume averaged BAO measured at the same redshift. While important for constraining cosmological parameters48, these differences are small compared to that caused by the big difference between the BAO and CMB redshifts.

Obtaining the contours and the r d(h) lines in Fig. 1

The marginalized joint rdH0 constraints from BAO were obtained using CosmoMC49 modified to work with rd as an independent parameter. The cosmological parameters we vary are rd, Ωmh2, and h, and the shown constraint is obtained after marginalizing over Ωmh2. The BAO data included the recently released Date Release (DR) 16 of the extended Baryon Oscillation Spectroscopic Survey (eBOSS)50 that includes BAO and redshift space distortions measurements at multiple redshifts from the samples of Luminous Red Galaxies (LRGs), Emission Line Galaxies (ELGs), clustering quasars (QSOs), and the Lyman-α forest. We use the BAO measurement from the full-shape auto- and cross-power spectrum of the eBOSS, LRGs, and ELGs51,52, the BAO measurement from the QSO sample53, and from the Lyman-α forest sample54. We combine these with the low-z BAO measurements by 6dF55 and the SDSS DR7 main Galaxy sample56.

The CMB and BAO lines shown in Fig. 1 were obtained by talking the measured value of θ or θ(zobs), fixing Ωmh2 at a certain value (provided for each line in the legend), varying h and deriving rd from Eqs. (3) and (5). We do not show the uncertainties around the individual lines because they are only meant to demonstrate the differences in slopes and the effect of different Ωmh2. The marginalized BAO and the Planck CMB contours provide a more accurate representation of the uncertainties involved.

The dependence of the CMB rd(h) lines on Ωmh2 may appear contradictory to the Ωmh2 dependence shown in Fig. 1 of a well-know paper by Knox and Millea40. There, increasing Ωmh2 moves the CMB best fit (rd, h) point in a direction orthogonal to where our CMB lines move. The reason for the difference is that their rd is a derived parameter obtained from the standard recombination model and, hence, depends on Ωmh2. In our derivation of the CMB lines, on the other hand, the Ωmh2 dependence only appears in D(z) and D(zobs).

Obtaining the S 8 constraints in Fig. 2

The joint DES + SN contours in Fig. 2 are obtained using the default version of CosmoMC and marginalizing over all relevant ΛCDM and nuisance parameters. To derive the Model 2 and Model 3 contours in Fig. 2, we fit the ΛCDM model to the BAO data using rd, Ωmh2, and h as a free parameters, supplemented by Gaussian priors on Ωmh2 and h, and with the primordial spectrum amplitude As and the spectral index ns fixed to their best fit ΛCDM values. The fit then generates constraints on S8 and Ωm as derived parameters. For Model 2, the Gaussian priors were Ωmh2 = 0.155 ± 0.0012, where we assumed the same relative uncertainty in Ωmh2 as for the Planck best fit ΛCDM model, and h = 0.71 ± 0.01, corresponding to the central value and the 1σ overlap between the CMB2 line and the BAO band. For Model 3, the priors were Ωmh2 = 0.167 ± 0.0013 and h = 0.735 ± 0.14.