Interacting dark energy in the early 2020s: a promising solution to the $H_0$ and cosmic shear tensions

TL;DR

This paper investigates a dark matter-dark energy interaction model that, when tested against recent cosmological data, alleviates the H_0 and cosmic shear tensions, suggesting a promising unified solution.

Contribution

It introduces and tests a specific dark sector interaction model with current observational data, demonstrating its potential to resolve major cosmological tensions.

Findings

Reduces H_0 tension to statistically plausible levels.

Significantly alleviates the sigma_8 tension.

Supports dark sector interactions as a viable solution.

Abstract

We examine interactions between dark matter and dark energy in light of the latest cosmological observations, focusing on a specific model with coupling proportional to the dark energy density. Our data includes Cosmic Microwave Background (CMB) measurements from the Planck 2018 legacy data release, late-time measurements of the expansion history from Baryon Acoustic Oscillations (BAO) and Supernovae Type Ia (SNeIa), galaxy clustering and cosmic shear measurements from the Dark Energy Survey Year 1 results, and the 2019 local distance ladder measurement of the Hubble constant $H_0$ from the Hubble Space Telescope. Considering Planck data both in combination with BAO or SNeIa data reduces the $H_0$ tension to a level which could possibly be compatible with a statistical fluctuation. The very same model also significantly reduces the $\Omega_{\rm m}-\sigma_8$ tension between CMB and cosmic shear measurements. Interactions between the dark sectors of our Universe remain therefore a promising joint solution to these persisting cosmological tensions.

Tables

Table 1. Table 1: Constraints on selected cosmological parameters within the coupled vacuum ξ ​ Λ 𝜉 Λ \xi\Lambda CDM model, considering either the Planck 2018 legacy dataset alone, or the same dataset in combination with: the CMB lensing power spectrum reconstructed from the CMB temperature four-point function; a combination of Baryon Acoustic Oscillation distance measurements; and distance moduli measurements from the Pantheon Supernovae Type Ia catalogue. Constraints are reported as 68% C.L. intervals. In the case of Ω c ​ h 2 subscript Ω c superscript ℎ 2 \Omega_{\rm c}h^{2} for the Planck and Planck + lensing dataset combinations the quantity quoted is the 95% C.L. upper limit.

Parameter	Planck	Planck+lensing	Planck+BAO	Planck+Pantheon
${\Omega }_b{h}^2$	$0.0224\pm 0.0002$	$0.0224\pm 0.0002$	$0.0224\pm 0.0002$	$0.0224\pm 0.0002$
${\Omega }_c{h}^2$			${0.095}_{-0.008}^{+0.022}$	${0.103}_{-0.007}^{+0.013}$
$n_s$	$0.966\pm 0.004$	$0.966\pm 0.004$	$0.965\pm 0.004$	$0.964\pm 0.004$
$\tau$	$0.054\pm 0.008$	$0.053\pm 0.007$	$0.054\pm 0.008$	$0.054\pm 0.008$
$\xi$	$-{0.54}_{-0.28}^{+0.12}$	$-{0.51}_{-0.29}^{+0.12}$	$-{0.22}_{-0.05}^{+0.21}$	$-{0.15}_{-0.07}^{+0.12}$
$H_0\left[\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}\right]$	${72.8}_{-1.5}^{+3.0}$	${72.8}_{-1.6}^{+3.3}$	${69.4}_{-1.5}^{+0.9}$	${68.6}_{-1.0}^{+0.8}$

Table 2. Table 2: Constraints on selected cosmological parameters within the coupled vacuum ξ ​ Λ 𝜉 Λ \xi\Lambda CDM model, considering three dataset combinations all involving the R19 Gaussian prior on H 0 subscript 𝐻 0 H_{0} based on the latest local distance measurement from HST , in combination with: the Planck 2018 legacy dataset (without CMB lensing); the Planck 2018 legacy dataset including the CMB lensing power spectrum reconstructed from the CMB temperature four-point function; and the Planck 2018 legacy dataset (without CMB lensing) in addition to a combination of Baryon Acoustic Oscillation distance measurements. Constraints are reported as 68% C.L. intervals. In the case of Ω c ​ h 2 subscript Ω c superscript ℎ 2 \Omega_{\rm c}h^{2} for the Planck + R19 and Planck + lensing + R19 dataset combinations the quantity quoted is the 95% C.L. upper limit.

Parameter	Planck+R19	Planck+lensing+R19	Planck+BAO+R19
${\Omega }_b{h}^2$	$0.0224\pm 0.0002$	$0.0224\pm 0.0002$	$0.0223\pm 0.0002$
${\Omega }_c{h}^2$			${0.062}_{-0.017}^{+0.020}$
$n_s$	$0.966\pm 0.004$	$0.967\pm 0.004$	$0.964\pm 0.004$
$\tau$	$0.053\pm 0.008$	$0.052\pm 0.007$	$0.053\pm 0.008$
$\xi$	$-{0.66}_{-0.13}^{+0.09}$	$-{0.62}_{-0.14}^{+0.09}$	$-0.47\pm 0.14$
$H_0\left[\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}\right]$	${74.0}_{-1.0}^{+1.2}$	${74.0}_{-1.1}^{+1.4}$	$71.7\pm 1.1$

Equations

$${\dot{\rho }}_c+3\mathcal{H}{\rho }_c$$

$${\dot{\rho }}_x+3\mathcal{H}(1+w){\rho }_x$$

$$Q=\xi \mathcal{H}{\rho }_x,$$

$${\dot{\delta }}_c$$

$${\dot{\theta }}_c$$

$${\dot{\delta }}_x$$

$${\dot{\theta }}_x$$

$${\delta }_x^{\mathrm{in}}(\eta )$$

$${\theta }_x^{\mathrm{in}}(x)$$

Full text

††journal: Phys. Dark Univ.

1 Introduction

The concordance $\Lambda$CDM cosmological model has been incredibly successful at describing cosmological observations at high and low redshift [1, 2, 3, 4, 5]. Yet, as uncertainties on cosmological parameters keep shrinking, a number of weaknesses have emerged: one of the most intriguing ones is the “$H_{0}$ tension”, referring to the mismatch between the value of the Hubble constant $H_{0}$ inferred from Planck Cosmic Microwave Background (CMB) data and direct local distance ladder measurements [6, 7]. In the past decade we have witnessed the tension between these two values grow in significance level from $2\sigma$ to $4.4\sigma$: the latest determinations from the Planck 2018 results and from the observations of Large Magellanic Cloud Cepheids by the Hubble Space Telescope (HST; measurement denoted as R19 hereafter) give $h=(0.6737\pm 0.0054)$ [8] and $h=(0.7403\pm 0.0142)$ [9] respectively, with $h\equiv H_{0}/(100\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1})$ the reduced Hubble constant. A very appealing possibility is that the $H_{0}$ discrepancy might be a hint of physics beyond the canonical $\Lambda$CDM model. The most economic possibilities in this direction involve phantom dark energy or some form of dark radiation, but a number of more complex scenarios have been studied, e.g. [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109].

On the other hand, tensions between cosmic shear surveys (such as [110, 111, 112]) and CMB measurements have also emerged [111, 113, 114]. For instance, the quantity $S_{8}\equiv\sigma_{8}(\Omega_{\rm m}/0.3)^{0.5}$ as measured by the KiDS weak lensing survey was shown to be in $2.6\sigma$ tension with the same quantity as measured by Planck [113, 111] (see also [115, 116] for previous analyses of CFHTLenS data). Focusing on the joint galaxy clustering and lensing likelihoods from the Dark Energy Survey (DES) [5, 117, 118], the Planck collaboration found modest tension with the DES results when galaxy clustering measurements are included, as the latter prefer an $\approx 2.5\sigma$ lower value of $S_{8}$ [8].

It is worth clarifying that the situation concerning the $S_{8}$/cosmic shear tension is less well defined with respect to that of the $H_{0}$ tension. In fact, various weak lensing analyses report varying degrees of tension. Earlier we have discussed the analyses which show strongest tension with Planck. However, other analyses show no tension, or only a mild amount of tension: these include for instance a combination of KiDS and GAMA galaxy clustering [119], a re-analysis of KiDS after accounting for a more careful treatment of intrinsic galaxy shape noise, multiplicative shear calibration uncertainties, and bin angular scale [120], and DES results when only considering cosmic shear but not galaxy clustering and shear-galaxy cross-correlations [5, 117, 118]. Despite these uncertainties, a number of exotic scenarios have been advocated in the past to alleviate the $S_{8}$ tension, see for instance [112, 14, 16, 121, 30, 33, 36, 122, 123, 24, 124, 41, 51, 59, 125, 126, 127].

Within the $\Lambda$CDM model, dark matter (DM) and dark energy (DE) behave as separate fluids not sharing interactions beyond gravitational ones. Historically, interactions between DM and DE were originally introduced to alleviate the coincidence problem [128, 129], although it is now understood that the energy exchange required for this problem to be addressed by these means is too large and observationally excluded. Nonetheless, interactions between DM and DE cannot be excluded on general grounds, and it is thus worth constraining them against the available wealth of precision cosmological data. This has motivated a large number of studies based on models where DM and DE share interactions, usually referred to as interacting dark energy (IDE) models (see e.g. [130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239], for a recent comprehensive review see [200]). Several studies in the literature have been devoted to exploring whether DM-DE interactions may help resolve the enduring $H_{0}$ tension, see e.g. [240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259].

In this work we (re)assess whether IDE cosmologies still provide a viable solution to the $H_{0}$ tension in light of the latest Planck and HST measurements. We find that IDE provides an interesting solution to the $H_{0}$ tension, which is brought below the $1\sigma$ level when considering Planck data alone. Intriguingly, when combining the latest Planck and HST measurements we find very strong indications for an interaction between the two dark components. These findings are, however, softened when including late-time measurements of the expansion history from Baryon Acoustic Oscillations and Supernovae Type Ia. We find that IDE also provides a promising solution to the $S_{8}$ tension between Planck and DES.

The rest of this paper is organized as follows. In Sec. 2 we revisit the basic background and perturbation equations of the IDE model we consider, discussing also stability considerations. In Sec. 3 we discuss the analysis method adopted, as well as the cosmological datasets considered. Our results are presented in Sec. 4, where we argue that the IDE model considered can partially address the $H_{0}$ and cosmic shear tensions, whereas possible model-dependence issues pertaining to the use of BAO and SNeIa data when constraining IDE models and assessing their ability to address the $H_{0}$ tension are discussed in Sec. 4.1. Finally, we provide closing remarks in Sec. 5.

2 Interacting dark energy revisited

We consider model featuring interactions between DM and DE with energy exchange proportional to the DM four-velocity, studied in earlier works such as [144, 260, 147, 154]. We assume a pressureless cold DM component and a DE component with equation of state (EoS) $w$, and denote the DM and DE energy densities by $\rho_{c}$ and $\rho_{x}$ respectively. At the background level, the DM-DE coupling modifies the continuity equations for the two dark fluids as follows [147]:

[TABLE]

where the dot denotes derivative with respect to conformal time $\tau$, and ${\cal H}\equiv\dot{a}/a$ is the conformal Hubble rate. In the notation of Eqs. (1,2), $Q>0$ and $Q<0$ indicate energy transfer from DE to DM and viceversa. We choose to focus on an interacting dark energy model, studied in various earlier works, wherein the coupling $Q$ takes the following form [144, 147]:

[TABLE]

where $\xi$ is a dimensionless coupling governing the strength of the DM-DE interaction. We shall refer to the resulting model as $\xi\Lambda$CDM model, or coupled vacuum model.

The presence of the DM-DE coupling also modifies the evolution of perturbations. In synchronous gauge, the linear perturbation equations for the evolution of the DM and DE density perturbations $\delta$ and velocity divergences $\theta$ are given by [144, 147, 261]:

[TABLE]

It is worth pointing out that the specific coupling chosen in Eq. (3) is purely phenomenological, introduced at the level of the continuity equations as shown in Eq. (1,2). In other words, it is not derived from a specific action (e.g. as in specific coupled quintessence models), and its phenomenological nature should be kept in mind. 111Nonetheless, given a certain field model of DM and DE, it is in principle possible to work backwards and reconstruct the specific type of Lagrangian interaction term between the DM and DE fields which could give rise to a chosen coupling function $Q$, as was done recently for instance in [262]. However, there is no guarantee that this procedure would return a Lagrangian term which is well-motivated from first-principle field-theoretical considerations. This explains for instance the difference between our results and those of the Planck collaboration, who in their 2015 dark energy and modified gravity paper [188] also considered constraints on coupled DE models. In [188] the Planck collaboration studied coupled quintessence models featuring an exponential coupling between a DM field and a scalar DE field (such a coupling could be motivated by Weyl scaling scalar-tensor theories), with an inverse power-law potential for the latter. The model and corresponding coupling functions were based on a specific well-motivated Lagrangian, whereas our choice of coupling function is phenomenological.

When introducing DM-DE interactions, we also need to modify the initial conditions for our Boltzmann system appropriately. To do so, we follow the earlier work of [261] and set adiabatic initial conditions for all species, which in turn requires appropriately modifying the initial conditions for $\delta_{x}$ and $\theta_{x}$. As shown in [261], following the gauge invariant formulation of [263], one can identify the appropriate initial conditions for the density contrast and velocity divergence of each species by writing the coupled Boltzmann equations in matrix form, solving the appropriate eigenvalue/eigenvector problem, and studying the modes which dominate the time evolution (which in this case will be those corresponding to the largest eigenvalue). These initial conditions turn out to be (see e.g. [261, 172, 244, 257, 264]):

[TABLE]

where $\eta=k\tau$ and $\delta_{\gamma}^{\rm in}(\eta)$ denotes the initial conditions for the photon density perturbations. We therefore set the initial conditions for the DE density contrast and velocity divergence following Eq. (9).

In the presence of DM-DE interactions, care must be given to the stability of the interacting system. For $w=-1$ (i.e. interacting vacuum), IDE models can suffer from gravitational instabilities [144, 265]. However, even when $w\neq-1$, one has to worry about early-time instabilities, leading to curvature perturbations blowing up on superhorizon scales. For the $\xi\Lambda$CDM model in which $Q\propto\rho_{x}$, the instabilities present when $w\neq-1$ are absent if the signs of $\xi$ and $(1+w)$ are opposite [144, 265, 266, 261, 161] (see also [267, 268, 269, 270, 271, 272, 273] for alternative approaches to avoiding these instabilities).

3 Methodology and cosmological observations

We consider an IDE model characterized by the coupling given by Eq. (3). The model is described by the usual six cosmological parameters of $\Lambda$CDM ($\Omega_{\rm b}h^{2}$, $\Omega_{\rm c}h^{2}$, $\theta_{\rm s}$, $A_{\rm s}$, $n_{\rm s}$, and $\tau$), in addition to the DM-DE coupling $\xi$. To circumvent the instability problem discussed above, we fix the DE EoS to $w=-0.999$. The rationale behind this approach (already followed in [172, 244]) is that for $w$ sufficiently close to $-1$ the effect of DE perturbations in Eqs. (6,7) is basically unnoticeable: consequently, these equations are essentially only capturing the effect of the DM-DE coupling $\xi$, while at the same time ensuring the absence of gravitational instabilities present when $w$ is strictly equal to $-1$. Such a model provides therefore a rather accurate surrogate for a $\Lambda$CDM+$\xi$ cosmology, and we shall refer to this model as $\xi\Lambda$CDM model (or coupled vacuum model). In order to avoid early-time instabilities discussed above, we then need to impose that $(1+w)$ and $\xi$ have opposite signs. Since we fixed $w=-0.999$, we impose $\xi<0$, and are therefore considering a model where energy flows from DM to DE. In a more realistic analysis, one should also vary the DE EoS (paying attention to the instability issue, which requires $\xi$ to change sign depending on whether the DE EoS lies in the quintessence or phantom regime), as we have done in our later companion paper [257].

Data-wise, we first consider measurements of CMB temperature and polarization anisotropies, as well as their cross-correlations, from the Planck 2018 legacy data release [8, 274]. This dataset is referred to as Planck TT,TE,EE+lowE in [8], whereas we refer to it simply as Planck. We then include measurements of the CMB lensing power spectrum reconstructed from the CMB temperature four-point function [274]. We refer to this dataset as lensing.

In addition to CMB data, we then consider Baryon Acoustic Oscillation (BAO) measurements from the 6dFGS [275], SDSS-MGS [276], and BOSS DR12 [4] surveys, a combination which we refer to as BAO. We further include Type Ia Supernovae (SNeIa) distance moduli measurements from the Pantheon sample [277], referring to this dataset as Pantheon. As pointed out in a number of earlier works [13, 278, 279, 49, 75], it is important to consider BAO and SNeIa distance measurements, from which one can construct an inverse distance ladder, when constraining late-time deviations from $\Lambda$CDM. In fact, studies based on the inverse distance ladder approach suggest that finding late-time solutions to the $H_{0}$ tension which can fit BAO and SNe data is challenging (albeit not impossible), see e.g. Fig. 1 in [75].

We also include galaxy clustering and cosmic shear measurements, as well as their cross-correlations, from the Dark Energy Survey combined-probe Year 1 results [5, 117, 118]. In particular, we consider the shear-shear, galaxy-galaxy, and galaxy-shear correlation functions, with the combination often referred to as the “$3\times 2$pt likelihood”. We refer to this dataset as DES.

Finally, we also consider a Gaussian prior on the Hubble constant $H_{0}=74.03\pm 1.42\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, consistent with the latest measurement by HST in [9]. We refer to this prior as R19. We focus on the ability of our model to reduce the tension with the R19 measurement (and subsequently consider the corresponding prior) as this is the measurement which is most discrepant with the Planck+BAO dataset combination assuming $\Lambda$CDM, and hence sets a harder task for our model. However, we note that there are alternative local measurements of $H_{0}$ which are also discrepant with the same dataset combination with various degrees of tension, including: local SNeIa calibrated with the Tip of the Red Giant Branch [280], strong lensing time delays from distant quasars as measured e.g. by the H0LiCOW team [281], and water megamasers [282]. See e.g. Fig. 1 in [283] for a visual summary of some of the main local measurements of $H_{0}$ beyond R19.

We modify the Boltzmann solver CAMB [284] to incorporate the effect of the DM-DE coupling as in Eqs. (4,7). We sample the posterior distribution of the cosmological parameters by making use of Markov Chain Monte Carlo (MCMC) methods, through a modified version of the publicly available MCMC sampler CosmoMC [285]. We monitor the convergence of the generated MCMC chains through the Gelman-Rubin parameter $R-1$ [286], requiring $R-1<0.02$ for our MCMC chains to be considered as converged. We impose flat priors on all cosmological parameters unless otherwise stated. In particular, as required by stability considerations, we impose $\xi<0$ at the prior level.

Finally, we use our MCMC chains to compute the Bayesian evidence for the $\xi\Lambda$CDM model (for different choices of datasets) using the MCEvidence code [287]. We then compute the (logarithm of the) Bayes factor with respect to $\Lambda$CDM, $\ln B$, with a value $\ln B>0$ indicating that the $\xi\Lambda$CDM model is preferred. We qualify the strength of the obtained values of $\ln B$ using the modified version of the Jeffreys scale provided in [288].

4 Results

Our main results are shown in Tab. 1, Tab. 2, and Fig. 1. In particular, in Tab. 1 we show results obtained either with the Planck dataset alone, or combining the latter with the lensing, BAO, and Pantheon datasets, one at a time. In Tab. 2 we instead show results involving the R19 prior, considering in particular the Planck+R19, Planck+lensing+R19, and Planck+BAO+R19 dataset combinations.

As shown in Tab. 1, from the Planck dataset alone, the value of the Hubble constant $H_{0}$ inferred within the $\xi\Lambda$CDM model is $H_{0}=72.8^{+3.0}_{-1.5}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. While the uncertainty is larger than that reported in [8] within the standard $\Lambda$CDM scenario, the central value has significantly shifted upwards. Indeed, this value is perfectly consistent with the HST measurement of [9], showing an agreement well below the $1\sigma$ level. Therefore, within the $\xi\Lambda$CDM model, the $H_{0}$ tension is solved when considering Planck data alone. In addition, we find $\xi=-0.54^{+0.12}_{-0.28}$, an apparent $>4\sigma$ detection of DM-DE interactions.

As we see from Tab. 1, including the lensing dataset does not alter our conclusions. In particular, from the Planck+lensing dataset combination we still find a strong detection of DM-DE interactions, with $\xi=-0.51^{+0.12}_{-0.29}$, and a high value of $H_{0}=72.8^{+3.3}_{-1.6}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. Including the lensing dataset has therefore mostly enlarged the uncertainties and very slightly driven the results towards $\Lambda$CDM, but to an extent which is unable to alter our conclusions neither qualitatively nor quantitatively.

The reason for such a high value of $H_{0}$ from CMB measurements alone can be found in the strong degeneracy between $H_{0}$ and $\xi$, as depicted in the left panel of Fig. 1. The origin of this degeneracy resides in the fact that for the $\xi\Lambda$CDM model, the background evolution of the DM energy density has an extra contribution proportional to the absolute value of $\xi$ and growing with $(1+z)^{3}$. Due to the presence of this extra term, the amount of DM today, $\Omega_{\rm c}$, must be smaller. However, the acoustic peak structure of the CMB (and in particular the relative height of odd and even peaks, as well as the overall height of all peaks) accurately fixes the value of $\Omega_{\rm c}h^{2}$: in order to accommodate a lower value of $\Omega_{c}$, a higher value of $H_{0}$ is required. An inverse correlation between $\xi$ and $H_{0}$ is therefore expected, and is reflected in the contours in the left panel of Fig. 1.

Note that even if the Planck dataset alone shows a preference for a non-zero negative $\xi$ at $>95\%$ C.L., this is likely due to a volume effect, i.e. more models with $\xi<0$ are compatible with Planck than models with $\xi=0$. This explanation is supported by the fact that the best-fit $\chi^{2}$ for $\xi\neq 0$ is almost the same as the best-fit $\chi^{2}$ for $\Lambda$CDM. Computing the Bayes factor for the $\xi\Lambda$CDM model with respect to $\Lambda$CDM for the Planck dataset we find $\ln B=1.3$. According to the modified Jeffreys scale of [288], this indicates a positive preference for the $\xi\Lambda$CDM model. When we also include the lensing dataset, we find that the Bayes factor is reduced to $\ln B=0.9$, which still indicates a preference for the $\xi\Lambda$CDM model, but at a weak level.

As the Planck and R19 datasets are now consistent, it is possible to combine them. When considering the Planck+R19 combination, we find an even stronger indication for non-zero $\xi$, inferring $\xi=-0.66^{+0.09}_{-0.13}$, an apparent $>5\sigma$ detection of DM-DE interactions. The Hubble constant is instead constrained to be $H_{0}=74.0^{+1.2}_{-1.0}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. Computing the Bayes factor, we find the extremely high value $\ln B=10.0$, indicating a very strong preference for the $\xi\Lambda$CDM model. Similarly, we can also consider the Planck+lensing+R19 dataset combination, given the mutual compatibility between Planck+lensing and R19. Doing so, we find, as already seen earlier, that the lensing dataset has mostly enlarged the uncertainties and very slightly drives the results towards $\Lambda$CDM. This is however insufficient to qualitatively or quantitatively alter our earlier conclusions, as we find $\xi=-0.62^{+0.09}_{-0.14}$ and $H_{0}=74.0^{+1.4}_{-1.1}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$.

As we discussed earlier, it is important to also consider BAO and SNeIa distance measurements when assessing the possibility of addressing the $H_{0}$ tension through late-time deviations from $\Lambda$CDM [13, 278, 279, 49, 75]. We therefore further investigate the impact of including late-time measurements of the expansion history from the BAO and Pantheon datasets. This results in all our previous findings being mildened. In fact, for the both the Planck+BAO and Planck+Pantheon dataset combinations, we infer a slightly lower value of $H_{0}$ within the $\xi\Lambda$CDM model than the value previously inferred from Planck alone. In particular, from the Planck+BAO dataset combination we find $H_{0}=69.4^{+0.9}_{-1.5}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, while from the Planck+Pantheon dataset combination we find $H_{0}=68.6^{+0.8}_{-1.0}\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. While both figures are clearly not high enough to indicate a strong resolution to the $H_{0}$ tension, they do nonetheless bring the tension down to the $\approx 3\sigma$ level. For the Planck+BAO dataset combination, in particular, the tension is brought down to the $2.6\sigma$ level, at which point one might start to view the $H_{0}$ tension as simply a statistical fluctuation within the $\xi\Lambda$CDM model.

For both the Planck+BAO and Planck+Pantheon dataset combinations, the hint for a non-zero coupling is still present, albeit at a much lower statistical significance, just above $1\sigma$. This is further confirmed when we compute the Bayes factor for the $\xi\Lambda$CDM model with respect to the $\Lambda$CDM model. We find $\ln B=-0.7$ for the Planck+BAO dataset combination, corresponding to a weak preference for $\Lambda$CDM, and $\ln B=-1.5$ for the Planck+Pantheon dataset combination, corresponding to a positive preference for $\Lambda$CDM. These results indicate that the increased model complexity of the $\xi\Lambda$CDM model is not warranted by a sufficient overall improvement in fit when including late-time measurements of the expansion history, unlike when the Planck dataset alone or the Planck+R19 dataset combination are considered.

As we see from the left panel of Fig. 1, the Planck+R19 and Planck+BAO dataset combinations are mildly consistent, since their $2\sigma$ credible regions overlap to a non-negligible extent. Furthermore, as we saw earlier, the value of $H_{0}$ obtained for the Planck+BAO dataset combination only shows a $2.6\sigma$ tension with R19. This means that we can consider the Planck+BAO+R19 dataset combination, although with caution. From this combination, we find $\xi=-0.47\pm 0.14$ and $H_{0}=71.7\pm 1.1\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. These figures are certainly intriguing and further highlight that, when the BAO dataset is included, the $\xi\Lambda$CDM model remains an interesting contender to address the $H_{0}$ tension (this is no longer true if the Pantheon dataset is included).

Overall, we find that when late-time measurements of the expansion history are included, our findings are mildened. While the $\xi\Lambda$CDM model clearly can no longer be considered an extremely promising solution to the $H_{0}$ tension, it is still able to bring the tension down to an interesting level, where it could be attributable (or at least almost attributable) to a statistical fluctuation. However, our findings also reinforce previous results indicating that finding late-time solutions to the $H_{0}$ tension is challenging (but not impossible) [13, 278, 279, 49, 75]. We also remind the reader that, within the $\xi\Lambda$CDM model, the Planck and BAO/Pantheon datasets are in mild tension, and it is therefore not completely clear whether one is allowed to combine them in first place. In Sec. 4.1, we discuss possible model-dependence issues pertaining to the use of BAO and SNeIa data when constraining IDE models and assessing their ability to address the $H_{0}$ tension.

The solution to the $H_{0}$ tension due to a lower intrinsic value for $\Omega_{\rm c}$ at present within the $\xi\Lambda$CDM model implies a much larger degeneracy in the $\Omega_{\rm m}-\sigma_{8}$ plane, reflected in the right panel of Fig. 1: the allowed contours from the Planck dataset follow a band, rather than reproducing the small region usually singled out. The reason is that once a coupling $\xi$ is switched on, the required DM energy density $\Omega_{\rm c}$ must be smaller as we have seen, implying that the clustering parameter $\sigma_{8}$ must be larger to have a proper normalization of the (lensing and clustering) power spectra. This effect can be understood from the scatter plot in the $\Omega_{\rm m}-\sigma_{8}$ plane in Fig. 2: as the absolute value of $\xi$ is increased, the allowed region bends towards larger (smaller) values of $\sigma_{8}$ ($\Omega_{\rm m}$).

The DES contours follow the expected $S_{8}\equiv\sigma_{8}(\Omega_{\rm m}/0.3)^{0.5}\simeq 0.79$ behavior [8] . Notice that the DES and Planck contours overlap for a very large fraction of the parameter space in the $\Omega_{\rm m}-\sigma_{8}$ plane, implying that the tension between Planck and DES is alleviated. Notice that this is not merely an effect due to the larger uncertainties in the Planck contours, but rather is due to the strong overlap between the two contours. We have found that removing one at a time the shear-shear, galaxy-galaxy, or galaxy-shear correlation functions from DES data does not qualitatively impact our results, with the overlap between the Planck and DES contours in the $\Omega_{\rm m}-\sigma_{8}$ plane remaining substantial.

4.1 Note on BAO and SNeIa measurements

A comment on BAO and SNeIa measurements is in order at this point. It is worth noting that the BAO measurements have been extracted under the assumption of $\Lambda$CDM. In fact, assuming a fiducial cosmology is needed both for converting the galaxy catalogue from celestial to comoving cartesian coordinates, as well as for performing the reconstruction procedure which sharpens the BAO peak. Therefore, fiducial cosmology assumptions and possible model-dependence certainly enter at multiple stages of BAO analyses. Whether or not this makes a significant difference in models beyond $\Lambda$CDM, is a non-trivial question which is still open. Aspects of this problem have been addressed in recent works, see e.g. [289, 290, 291, 292].

In particular [291] has shown that, in the presence of models which deviate significantly from $\Lambda$CDM at late times (in particular if large metric gradients are present), the usual Alcock-Paczynski scaling underlying the standard BAO measurements breaks down. In other words, the two key assumptions in BAO analyses that 1) differences in the distance-redshift relationship between true and fiducial cosmology scale linearly, and 2) differences in comoving clustering between true and fiducial cosmology can be nulled by the same free parameters used to null the non-BAO signal in the 2-point correlation function or galaxy power spectrum (see e.g. [290]), break down. As a result, it was suggested in [291] that current BAO measurements should not be extrapolated to constrain cosmologies where late-time distance measures differ more than a few percent from $\Lambda$CDM. At face value, this issue would appear to apply to the $\xi\Lambda$CDM model, where late-time distance measures deviate at percent-level from $\Lambda$CDM if the $H_{0}$ tension is to be solved. Work to quantify the extent to which this possible residual BAO model-dependence affects our results, and in particular how much BAO measurements can be trusted to constrain the $\xi\Lambda$CDM model, is underway, and we hope to report on this in the near future. In the meantime, we advise caution with regards to the interpretation of results stemming from dataset combinations involving the BAO dataset in this work.

With regards to SNeIa distance measurements, no specific model assumption was made in obtaining these. However, the interpretation thereof, or in other words the Pantheon likelihood, assumes that intrisic SNeIa luminosities do not evolve with redshift. Whether or not this is actually valid has been the subject of recent debate (see e.g. [293, 294]). More importantly, even if astrophysics is not responsible for redshift evolution of intrinsic SNeIa luminosities, such an evolution can occur in models of DE beyond the cosmological constant, particularly if these lead to time-dependent fundamental constants or to coupling between DE and baryons (as for instance in [295, 296, 297]): see [298, 299] for examples of the impact on SNeIa luminosities.

Since we have not assumed our coupling function $Q$ to arise from any particular Lagrangian, it is hard (if not impossible) to assess what impact our IDE model should have on intrinsic SNeIa luminosities. In [300], a model-agnostic parametrization of a possible redshift-dependence in intrinsic SNeIa luminosities in the Pantheon likelihood showed that for the $\xi\Lambda$CDM model in question, the effect of this possible systematic is to slightly raise $H_{0}$ and enlarge the uncertainties thereof (see Tabs. 6 and 7 of [300]). This small but noticeable effect could help make the case for IDE models providing an interesting route towards addressing the $H_{0}$ tension even when SNeIa luminosity distance measurements are considered.

In summary, we certainly believe that the inclusion of BAO and SNeIa data when constraining late-time deviations from $\Lambda$CDM is extremely important, particularly in terms of assessing the feasibility of addressing the $H_{0}$ tension within these models. However, care must be taken with regards to the interpretation of these datasets and their use in constraining models beyond $\Lambda$CDM. This is the case for BAO measurements if late-time distance measures differ more than a few percent from $\Lambda$CDM in the target cosmology, or for SNeIa measurements if the DE component in the target cosmology couples to baryons and/or leads to time-dependent fundamental constants.

5 Conclusions

In this work, we have examined the persisting $H_{0}$ tension in light of the Planck 2018 legacy data release, late-time distance measurements of the expansion history from BAO and SNeIa, and the latest $1\%$ determination of $H_{0}$ from HST. We find that within an interacting dark energy model studied in previous works, with coupling function given by Eq. (3), the value of $H_{0}$ inferred by Planck is consistent with the latest local distance measurement well within $1\sigma$, representing a solution to the $H_{0}$ tension.

Our findings are, however, mildened when we include late-time BAO and SNeIa measurements. In this case, the value of $H_{0}$ we infer is not high enough that we can claim a definitive solution to the $H_{0}$ tension. Nonetheless, the tension is brought down to a level where it could be compatible with a statistical fluctuation within the interacting dark energy model we have considered, although the model is not favoured by Bayesian evidence considerations. At the same time, our interacting dark energy model appears promising in terms of alleviating the tensions between CMB and cosmic shear measurements. In particular, we observe a considerably improved overlap between the Planck and DES contours in the $\Omega_{\rm m}-\sigma_{8}$ plane.

Possible model-dependence issues pertaining to the use of BAO and SNeIa data when constraining interacting dark energy models and assessing their ability to address the $H_{0}$ tension are discussed in Sec. 4.1. It is certainly very important to include BAO and SNeIa measurements, from which an inverse distance ladder can be constructed, when constraining late-time deviations from $\Lambda$CDM. Nonetheless, we feel like advising slightly more caution when interpreting our results involving either or both the BAO or Pantheon datasets, due to these possible model-dependence issues we have discussed, and to their tension with the Planck dataset when interpreted within the framework of an interacting dark energy model.

It is worth keeping in mind that we have not varied the dark energy equation of state but fixed it to $w=-0.999$ (which in turn requires $\xi<0$), an approach already followed in the earlier [172, 244]. As explained at the start of Sec. 3, we did so to provide a surrogate for a coupled vacuum model, while bypassing the usual instabilities problem of the latter [144]. The case where $w$ is varied as well was considered in our later companion paper [257], and shown to decrease the Bayes factor in favour of the resulting IDE model. It is also worth recalling that the dark matter-dark energy coupling we have studied is purely phenomenological, introduced at the level of the continuity equations, and not coming from a specific Lagrangian.

To conclude, it is intriguing that the interacting dark energy model we have considered provides not only an interesting solution to the $H_{0}$ tension, but at the same time can alleviate the tension between CMB and cosmic shear measurements. Our findings suggest that the possibility of late-time new physics in the dark sector hiding behind these persisting cosmological tensions cannot yet be discarded. Future data (e.g. from Euclid or from gravitational wave standard sirens) will be vital towards further confirming the trends we have seen, or definitively discarding the model [228, 264].

Acknowledgements

E.D.V. acknowledges support from the European Research Council in the form of a Consolidator Grant with number 681431. A.M. is supported by TASP, iniziativa specifica INFN. O.M. is supported by PROMETEO II/2014/050, by the Spanish Grant FPA2017-85985-P of the MINECO, by the MINECO Grant SEV-2014-0398 and by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements 690575 and 674896. S.V. is supported by the Isaac Newton Trust and the Kavli Foundation through a Newton-Kavli fellowship, and acknowledges a College Research Associateship at Homerton College, University of Cambridge. O.M. would like to thank the hospitality of the Fermilab Theory Department. This work is based on observations obtained with Planck (www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We acknowledge the use of the Planck Legacy Archive.