# Coupled and uncoupled early dark energy, massive neutrinos and the
cosmological tensions

Adrià Gómez-Valent1,2<EMAIL_ADDRESS>Ziyang Zheng3<EMAIL_ADDRESS>heidelberg.de Luca Amendola3<EMAIL_ADDRESS>Christof
Wetterich3<EMAIL_ADDRESS>Valeria Pettorino4
<EMAIL_ADDRESS>1 Dipartimento di Fisica, Università di Roma Tor
Vergata, via della Ricerca Scientifica 1, I-00133 Roma, Italy  2 INFN, Sezione
di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica 1,
I-00133 Roma, Italy  3 Institut für Theoretische Physik, Ruprecht-Karls-
Universität Heidelberg, Philosophenweg 16, D-69120 Heidelberg, Germany  4 AIM,
CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris
Cité, F-91191 Gif-sur-Yvette, France

###### Abstract

Some cosmological models with non-negligible dark energy fractions in
particular windows of the pre-recombination epoch are capable of alleviating
the Hubble tension quite efficiently, while keeping the good description of
the data that are used to build the cosmic inverse distance ladder. There has
been an intensive discussion in the community on whether these models enhance
the power of matter fluctuations, leading de facto to a worsening of the
tension with the large-scale structure measurements. We address this pivotal
question in the context of several early dark energy (EDE) models, considering
also in some cases a coupling between dark energy and dark matter, and the
effect of massive neutrinos. We fit them using the Planck 2018 likelihoods,
the supernovae of Type Ia from the Pantheon compilation and data on baryon
acoustic oscillations. We find that ultra-light axion-like (ULA) EDE can
actually alleviate the $H_{0}$ tension without increasing the values of
$\sigma_{12}$ with respect to those found in the $\Lambda$CDM, whereas EDE
with an exponential potential does not have any impact on the tensions. A
coupling in the dark sector tends to enhance the clustering of matter, and the
data limit a lot the influence of massive neutrinos, since the upper bounds on
the sum of their masses are too close to those obtained in the standard model.
We find that in the best case, namely ULA, the Hubble tension is reduced to
$\sim 2\sigma$.

Cosmology: observations – Cosmology: theory – cosmological parameters – dark
energy – dark matter

##  I Introduction

The standard model of cosmology, also known as $\Lambda$CDM, relies, among
other assumptions, on the existence of cold dark matter (CDM) and a constant
(immutable) vacuum energy density that pervades space and is in charge of the
late-time acceleration of the universe Peebles and Ratra (2003). Despite its
apparent simplicity, this model has shown an undeniable ability to explain a
wide variety of cosmological observations with an incredible accuracy, ranging
from the anisotropies of the cosmic microwave background (CMB) Hinshaw  _et
al._ (2013); Aghanim  _et al._ (2020), to the baryon acoustic oscillations
(BAO) imprinted in the distribution of matter in the universe Eisenstein  _et
al._ (2005); Cole  _et al._ (2005), and to the Hubble diagram of supernovae of
Type Ia (SNIa) Perlmutter  _et al._ (1999); Riess  _et al._ (1998). This is
why the $\Lambda$CDM is sometimes also referred to as the concordance model of
cosmology. However, apart from very important and long-standing theoretical
problems associated to the value of the vacuum energy density (see e.g.
Weinberg (1989) and the review Solà Peracaula (2022)), we have also witnessed
in the last years the stubborn persistence of significant tensions between the
model and some data, which make the $\Lambda$CDM less concordant than
previously thought. Two of the most prominent ones are the Hubble (or $H_{0}$)
Verde  _et al._ (2019) and $\sigma_{8}$ (or $S_{8}$) tensions Di Valentino
_et al._ (2021a). These are the tensions we will focus our attention on in
this work, although there are others, as those linked to the CMB anomalies,
see e.g. Perivolaropoulos and Skara (2022); Aluri  _et al._ (2022). The
$H_{0}$ tension is between the Hubble constant measured by the SH0ES team
model-independently with the cosmic distance ladder method Riess  _et al._
(2022), and the CMB-inferred value by the Planck collaboration Aghanim  _et
al._ (2020), which assumes $\Lambda$CDM. These values read, respectively,
$H_{0}=(73.04\pm 1.04)$ km/s/Mpc and $H_{0}=(67.36\pm 0.54)$ km/s/Mpc111Here,
as well as in the evaluation of the $S_{8}$ tension, we take the value
obtained from the TT,TE,EE+lowE+lensing analysis by Planck Aghanim  _et al._
(2020)., and are in $4.8\sigma$ tension. On the other hand, the $S_{8}$
parameter is defined as $S_{8}=\sigma_{8}(\Omega_{m}^{(0)}/0.3)^{0.5}$, with
$\Omega_{m}^{(0)}$ the current matter density fraction, $\sigma_{8}$ the rms
of mass fluctuations at scales of $R_{8}=8h^{-1}$ Mpc, and $h$ the reduced
Hubble parameter. The $S_{8}$ tension is found between, again, the
$\Lambda$CDM-based inference by Planck Aghanim  _et al._ (2020),
$S_{8}=0.832\pm 0.013$, and the value measured by some weak lensing surveys.
For instance, the combined tomographic weak gravitational lensing analysis of
the Kilo Degree Survey (KiDS +VIKING-450) and the Dark Energy Survey (DES-Y1)
Joudaki  _et al._ (2020) led to the measurement
$S_{8}=0.762^{+0.025}_{-0.024}$, under the assumption of the standard model.
The latter is in $2.5\sigma$ tension with Planck. This tension is actually
still compatible with a statistical fluctuation Nunes and Vagnozzi (2021),
although the first hints for a tension between the $\Lambda$CDM and the large-
scale structure (LSS) data appeared already almost one decade ago Macaulay
_et al._ (2013) and have persisted over the years Gil-Marín  _et al._ (2017);
Gómez-Valent and Solà (2017); Nesseris  _et al._ (2017); Gómez-Valent and Solà
Peracaula (2018).

These tensions, if not explained by unaccounted systematic uncertainties in
the data, might be calling for new physics to reestablish the cosmic
concordance. Many models have been already invoked to try to alleviate them,
with more or less success, cf. the dedicated reviews Perivolaropoulos and
Skara (2022); Di Valentino  _et al._ (2021b). Some of the most promising
proposals modify the physics in the pre-recomination epoch either by
considering changes in the gravitational theory Solà Peracaula  _et al._
(2019a); Ballesteros  _et al._ (2020); Braglia  _et al._ (2020); Solà
Peracaula  _et al._ (2020); Braglia  _et al._ (2021); Benevento  _et al._
(2022), an early dark energy (EDE) component Poulin  _et al._ (2019);
Niedermann and Sloth (2021); Ye and Piao (2020); Gómez-Valent  _et al._
(2021); Wang and Piao (2022), the effect of primordial magnetic fields
Jedamzik and Pogosian (2020), varying atomic constants Liu  _et al._ (2020);
Sekiguchi and Takahashi (2021) or running vacuum models Solà Peracaula  _et
al._ (2021). These changes are usually accompanied also by some features in
the post-recombination era, and the fitting performance of the various
proposals can vary significantly depending on the particular data set under
consideration.

In this paper we study some cosmological models with a dark energy (DE) scalar
field in interaction with dark matter, in which the DE does not only play a
role in the late-time universe, but can also leave some imprints in the
radiation- and matter-dominated epochs (RDE and MDE, respectively). These
models fall into the category of early dark energy models. While EDE is
already present in standard coupled DE scenarios Wetterich (1995); Amendola
(2000a), here we consider the effect of scalar field potentials that have a
direct impact on the universe’s dynamics already in the RDE and/or MDE. Recent
investigations of uncoupled EDE scenarios with a non-negligible peak of the
EDE fraction before recombination have proved to be efficient in loosening the
Hubble tension Poulin  _et al._ (2019, 2021); Hill  _et al._ (2022); Smith
_et al._ (2022), while respecting the tight upper bounds on the EDE fraction
at the decoupling time Gómez-Valent  _et al._ (2021). This is because they are
able to increase the energy budget of the universe in that era and, hence,
decrease the sound horizon at the baryon-drag epoch, which forces a larger
Hubble rate at low redshifts to explain the correct location of the CMB peaks
and the BAO data. Some authors have pointed out, though, that these models
might worsen the tension with the data on LSS, see e.g. Hill  _et al._ (2020);
Gómez-Valent  _et al._ (2021). There is still no clear consensus on these
matters (see Smith  _et al._ (2021); Murgia  _et al._ (2021)), and further
research is needed to shed some light on the problem at hand. In particular,
in this work we address the question whether the LSS estimators
$\sigma_{8}/S_{8}$ are appropriate for quantifying the tension with the LSS
data. We will see that, as firstly argued in Sánchez (2020), they could lead
to biased results when the posterior distribution for $H_{0}$ is significantly
different from the $\Lambda$CDM one. If so, this can have of course a direct
impact on our understanding of the cosmological tensions, as we already showed
in Gómez-Valent  _et al._ (2021). The LSS tension could be directly
intertwined with the Hubble tension. This issue has been discussed in the
context of the ultra-light axion-like (ULA) EDE in Gómez-Valent (2022a), where
the author concluded that actually this model is able to alleviate the Hubble
tension while keeping the mass fluctuations at scales of $12$ Mpc very close
to the preferred value in the $\Lambda$CDM and, hence, without any clear
enhancement of the clustering at these scales. Significant differences at
other scales might exist, though. ULA typically enhances the power spectrum at
small scales, and suppresses it at large scales. In this paper we will
elaborate more on this idea.

We concluded in our previous work Gómez-Valent  _et al._ (2021) that a
coupling between DE and DM could help us, in principle, to alleviate both
tensions at a time, since it could allow the dark matter energy density to be
larger than in the standard model at the last-scattering surface (which is
needed to counteract the enhancement of the early integrated Sachs-Wolfe
effect Vagnozzi (2021)) while retrieving standard values of $\omega_{\rm
cdm}$. This possibility has been already explored in the literature Karwal
_et al._ (2022); McDonough  _et al._ (2022), where the authors obtained
constraints using data on CMB and other data sets, incorporating in all cases
the prior on $H_{0}$ from SH0ES Riess  _et al._ (2022). Here we perform our
fitting analyses sticking to the combination CMB+SNIa+BAO, which contains the
main ingredients to build the inverse distance ladder Aubourg  _et al._
(2015); Cuesta  _et al._ (2015); Feeney  _et al._ (2019); Camarena and Marra
(2020a). We do not include any prior on the Hubble constant. We deem this
important in order to assess the real ability of these models to alleviate the
tensions by checking what is the actual room for shifts in the parameters of
interest according to this robust cosmological data set. We provide
constraints not only on $\sigma_{8}/S_{8}$, but also on $\sigma_{12}$ and the
related quantity $S_{12}$ (see Sánchez (2020)), and discuss the Hubble tension
also in terms of the absolute magnitude of SNIa. Constraints on these
parameters have not been previously reported for the coupled and uncoupled EDE
models explored in this paper.

We study ULA with and without coupling in the dark sector, but not only. In
Gómez-Valent  _et al._ (2021) we showed that EDE models with scaling solutions
during the matter- and radiation-dominated eras, as those found when the
scalar field dynamics is governed by an exponential potential in the MDE and
RDE, are not efficient concerning the tensions. Here we consider the case of
DE with an exponential potential during the MDE and RDE with and without
coupling with DM, and with a cosmological constant triggering acceleration at
late times. We present constraints on it, and compare our results with the
$\Lambda$CDM, the uncoupled model, and the coupled model with a trivial flat
potential, together with the coupled and uncoupled ULA.

On top of that, we study the effect of massive neutrinos in our analyses of
uncoupled ULA and coupled quintessence with a cosmological constant. It is
very well known that massive neutrinos suppress the matter power spectrum at
low scales, see e.g. Lesgourgues and Pastor (2006), so they could compensate
some of the enhancement effects found in ULA when a minimal setup with only
one massive neutrino of 0.06 eV is considered, and also those introduced by
the fifth force in coupled dark energy models. Understanding up to what extent
this is possible is one of the main goals of this work.

The paper is organized as follows. In Sec. II we present the main cosmological
equations at the background and linear perturbations level for the coupled
dark energy models studied in this work. In Sec. III we discuss the most
relevant phenomenological aspects of these models in order to understand their
main signatures in the cosmological observables. In Sec. IV we describe the
data sets and the methodology employed to constrain the models, and in Sec. VI
we present and discuss our results. Our conclusions are provided in Sec. VII.

##  II Cosmological equations in coupled dark energy models

We consider a perturbed flat Friedmann-Lemaître-Robertson-Walker universe with
the spacetime line element written in the synchronous gauge,

$ds^{2}=a^{2}(\tau)[-d\tau^{2}+\\{\delta_{ij}+h_{ij}(\tau,\vec{x})\\}dx^{i}dx^{j}]\,,$
(1)

with $\tau$ the conformal time and $\vec{x}$ the spatial comoving vector. The
two scalar degrees of freedom contained in the matrix $h_{ij}$,

$h_{ij}=\int
d^{3}k\,e^{-i\overrightarrow{k}\cdot\overrightarrow{x}}\left[h\,\hat{k}_{i}\hat{k}_{j}+6\eta\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)\right]\,,$
(2)

are its trace $h(\tau,\vec{k})$ and $\eta(\tau,\vec{k})$, with $\vec{k}$ the
comoving wavenumber Ma and Bertschinger (1995). We do not discuss here the
vector and tensor perturbations.

Regarding the matter content in the universe, we consider the fields described
in the standard model of particle physics and some extension accounting for
the neutrino masses, cold dark matter and dark energy. DM is treated as a
pressureless perfect fluid, whereas DE is described in terms of a scalar field
$\phi$ with an associated potential $V(\phi)$. We consider these two
components to be in interaction. The cross-talk between DM and DE modifies
their conservation equations, which take the following form,

$\nabla^{\mu}T^{\phi}_{\mu\nu}=+Q_{\nu}\qquad;\qquad\nabla^{\mu}T^{\rm
dm}_{\mu\nu}=-Q_{\nu}\,.$ (3)

The source vector $Q_{\nu}$ is determined by the concrete nature of the
interaction. In this work we consider the conformal coupling, namely

$Q_{\nu}=\kappa\beta T^{\rm dm}\nabla_{\nu}\phi\,,$ (4)

with $\kappa=\sqrt{8\pi G}=m_{P}^{-1}$ the inverse of the reduced Planck mass
and $T^{\rm dm}$ the trace of the DM energy-momentum tensor. We use, for
simplicity, a constant dimensionless coupling $\beta$. The 0-component of (3)
leads to a modified equation describing the anomalous dilution law for DM,

$\bar{\rho}^{\prime}_{\rm dm}+3\mathcal{H}\bar{\rho}_{\rm
dm}=-\beta\kappa\bar{\rho}_{\rm dm}\phi^{\prime}\,,$ (5)

and the modified Klein-Gordon (KG) equation for the scalar field,

$\phi^{\prime\prime}+2\mathcal{H}\phi^{\prime}+a^{2}\frac{dV}{d\phi}=\beta\kappa
a^{2}\bar{\rho}_{\rm dm}\,.$ (6)

The primes denote derivatives with respect to the conformal time and
$\mathcal{H}=a^{\prime}/a$. This two equations are valid at the background
level, with $\phi$ the mean (background) value of the scalar field and
$\bar{\rho}_{\rm dm}$ the background DM energy density. Eq. (5) can be
trivially solved in terms of $\phi$. If we assume that the number of DM
particles is conserved throughout the cosmic history the corresponding number
density reads $n_{\rm dm}(a)=n_{\rm dm}^{(0)}a^{-3}$ and we find that their
mass evolves according to

$m_{\rm dm}(\phi)=m_{\rm dm}^{(0)}e^{\beta\kappa(\phi^{(0)}-\phi)}\,,$ (7)

with the superscripts $(0)$ denoting current quantities. The dynamics of
$\phi$ is governed by eq. (6), and depends, of course, on the intensity of the
interaction and the particular shape of the scalar field potential. These
features characterize the model and determine the evolution of the DM mass
with the expansion as well. The Friedmann and pressure equations take the same
form as in the $\Lambda$CDM, but substituting the constant dark energy density
$\rho_{\Lambda}=\Lambda/\kappa^{2}$ and pressure $p_{\Lambda}=-\rho_{\Lambda}$
by

$\bar{\rho}_{\phi}=\frac{(\phi^{\prime})^{2}}{2a^{2}}+V(\phi)\qquad{\rm
and}\qquad\bar{p}_{\phi}=\frac{(\phi^{\prime})^{2}}{2a^{2}}-V(\phi)\,.$ (8)

At the linear perturbations level, and following the notation of Ma and
Bertschinger (1995), we obtain in momentum space the following set of coupled
differential equations Amendola (2000b):

$\mathcal{H}h^{\prime}-2\eta
k^{2}=\kappa^{2}a^{2}\left[\delta\rho+\frac{\phi^{\prime}}{a^{2}}\delta\phi^{\prime}+\frac{dV}{d\phi}\delta\phi\right]\,,$
(9)

$2\eta^{\prime}k^{2}=\kappa^{2}[a^{2}(\bar{\rho}+\bar{p})\theta+k^{2}\phi^{\prime}\delta\phi]\,,$
(10)

$-h^{\prime\prime}-2\mathcal{H}h^{\prime}+2\eta
k^{2}=3\kappa^{2}a^{2}\left[\delta
p+\frac{\phi^{\prime}}{a^{2}}\delta\phi^{\prime}-\frac{dV}{d\phi}\delta\phi\right]\,,$
(11)

$\kappa\bar{\rho}_{\rm dm}a^{2}\beta\delta_{\rm
dm}=\delta\phi^{\prime\prime}+2\mathcal{H}\delta\phi^{\prime}+\left(k^{2}+a^{2}\frac{d^{2}V}{d\phi^{2}}\right)\delta\phi+\frac{h^{\prime}}{2}\phi^{\prime}\,,$
(12)

$\delta_{\rm dm}^{\prime}=-\left(\theta_{\rm
dm}+\frac{h^{\prime}}{2}\right)-\kappa\beta\delta\phi^{\prime}\,,$ (13)

$0=\theta_{\rm dm}^{\prime}+(\mathcal{H}-\beta\kappa\phi^{\prime})\theta_{\rm
dm}+k^{2}\beta\kappa\delta\phi\,.$ (14)

The first three correspond to the 00, 0i and ii perturbed Einstein equations,
respectively. The fourth one is the perturbed KG equation, and the last two
are obtained from the perturbed 0 and i components of the conservation
equation for dark matter. The conservation equations for the other species
remain the same as in the standard model at all orders in perturbation theory,
since there is no direct coupling between them and $\phi$. The quantities
$\delta\rho\equiv\sum_{j}\delta\rho_{j}$ and $\delta p\equiv\sum_{j}\delta
p_{j}$ are the sum of the perturbed densities and pressures of the various
cosmological components, without including the scalar field contribution;
$\delta_{\rm dm}$ and $\theta_{\rm dm}$ are the DM density contrast and
velocity gradient, respectively; $\delta\phi$ is the perturbation of the
scalar field; and
$(\bar{p}+\bar{\rho})\theta\equiv\sum_{j}(\bar{p}_{j}+\bar{\rho}_{j})\theta_{j}$.
We use adiabatic initial conditions.

In order to see how the coupling impacts the large-scale structure formation
processes it is illustrative to obtain the equations for the baryon and dark
matter density contrasts at deep subhorizon scales in the matter- and DE-
dominated universe, when radiation can be safely neglected. They read

$\delta_{\rm b}^{\prime\prime}+\mathcal{H}\delta^{\prime}_{\rm
b}-\frac{\kappa^{2}a^{2}}{2}[\bar{\rho_{b}}\delta_{b}+\bar{\rho}_{\rm
dm}\delta_{\rm dm}]=0\,,$ (15) $\delta_{\rm
dm}^{\prime\prime}+\left[\mathcal{H}-\frac{\beta\kappa\phi^{\prime}k^{2}}{k^{2}+a^{2}m_{\phi}^{2}}\right]\delta^{\prime}_{\rm
dm}+(\delta_{b}^{\prime}-\delta^{\prime}_{\rm
dm})\frac{\beta\kappa\phi^{\prime}a^{2}m_{\phi}^{2}}{k^{2}+a^{2}m_{\phi}^{2}}-\frac{\kappa^{2}a^{2}}{2}\left[\bar{\rho}_{b}\delta_{b}+\bar{\rho}_{\rm
dm}\delta_{\rm
dm}\left(1+\frac{2\beta^{2}k^{2}}{k^{2}+a^{2}m^{2}_{\phi}}\right)\right]=0\,.$
(16)

These equations encode the non-trivial growth of the matter perturbations,
which depends a lot on the details of the model under study. Notice that we
have also included the effect of the DE mass $m^{2}_{\phi}(a)\equiv
d^{2}V/d\phi^{2}$. A larger mass confines the fifth force effects to smaller
scales, i.e. to larger $k$’s. Nevertheless, the DE potentials that are able to
trigger the late-time acceleration of the universe must be extremely flat at
late times to explain the current data, see e.g. Aghanim  _et al._ (2020);
Solà Peracaula  _et al._ (2019b). Hence, $m_{\phi}$ can be neglected at low
redshifts. However, if the shape of $V(\phi)$ allows for a non-negligible
effective DE mass in previous stages of the cosmic expansion, competing with
the physical wave modes $k/a$ inside the horizon, this could leave an imprint
on the large-scale structure of the universe, so it is useful to keep the mass
terms in these equations and discuss their effects on a case-by-case basis.

It is easy to check that in the limit $\beta\to 0$ of Eqs. (15)-(16) we
recover the results of the $\Lambda$CDM and uncoupled dark energy models, and
in the limit $m_{\phi}\to 0$ we retrieve the well-known result of coupled dark
energy with a nearly flat potential, see e.g. Barros  _et al._ (2019); Gómez-
Valent  _et al._ (2020),

$\delta_{\rm
dm}^{\prime\prime}+\left[\mathcal{H}-\beta\kappa\phi^{\prime}\right]\delta^{\prime}_{\rm
dm}-\frac{\kappa^{2}a^{2}}{2}\left[\bar{\rho}_{b}\delta_{b}+\bar{\rho}_{\rm
dm}\delta_{\rm dm}(1+2\beta^{2})\right]=0\,.$ (17)

By combining (9) and (12) one can see that in these coupled DE scenarios the
dark energy component does not cluster at deep subhorizon scales, since when
$k^{2}\gg\mathcal{H}^{2}$ we have

$\delta\rho_{\phi}=\frac{\phi^{\prime}}{a^{2}}\delta\phi^{\prime}+\frac{dV}{d\phi}\delta\phi\sim\left(\frac{H}{k}\right)^{2}\times\left[\beta^{2},\beta\alpha\right]\times\mathcal{O}(\delta\rho_{m})\,,$
(18)

with $\alpha$ the slope of the potential. This means that galaxies are tracers
of the underlying distribution of dark matter, as in the $\Lambda$CDM.

In this work we study models with a potential of the following form,

$V(\phi)=V_{0}+V_{\rm ede}(\phi)\,.$ (19)

By construction, $V_{\rm ede}$ is not important at $z\lesssim 1$, but can have
a sizable impact at $z\gg 1$, i.e. during the matter and/or radiation-
dominated epochs. $V_{0}$, instead, is a constant term that is responsible of
the current accelerated phase of the universe. Therefore,
$V_{0}=\mathcal{O}(m_{P}^{2}H_{0}^{2})$ and dominates the expansion at late
times. It is basically the energy density associated to the cosmological
constant. For simplicity we do not consider more complicated shapes of the
potential at low redshifts, just not to introduce possible degeneracies with
the features that characterize $V_{\rm ede}$. This would hinder the
interpretation of our results.

We focus in this paper on three alternative forms of $V_{\rm ede}$ with and
without considering the effect of the coupling in the dark sector, i.e.
studying both the cases with $\beta=0$ and $\beta\neq 0$. We dedicate the next
section to list them and explain their main phenomenological aspects.

##  III Phenomenology of several coupled and uncoupled EDE models in a
nutshell

###  III.1 Constant potential

Standard model: If $V(\phi)=V_{0}$ and $\beta=0$ the model reduces to the
$\Lambda$CDM, since in this case the DE has no dynamics and there is no
interaction between DM and DE. This is the only model studied in this paper in
which the dark energy density is a pure rigid constant at all times.

Coupled dark energy with a flat potential: If $V(\phi)=V_{0}$ but $\beta\neq
0$ the interaction in the dark sector is active, and this gives rise to a very
rich phenomenology. The first coupled DE models of this sort where presented
and studied in Wetterich (1995); Amendola (2000a). See also Amendola and
Tsujikawa (2015). In this paper we refer to the model of coupled DE with flat
potential as CDE_const. For values of the coupling of the order
$\beta\lesssim\mathcal{O}(10^{-1})$, as the ones preferred by the data
Amendola (2001); Pettorino  _et al._ (2012); Pettorino (2013); Xia (2013); Ade
_et al._ (2016); van de Bruck  _et al._ (2017); van de Bruck and Mifsud
(2018); Barros  _et al._ (2019); Agrawal  _et al._ (2021); Gómez-Valent  _et
al._ (2020); Gómez-Valent (2022a); Goh  _et al._ (2022), the scalar field has
a negligible energy density before the matter-radiation equality time, $t_{\rm
eq}$. In particular, its kinetic energy is derisory, and we can take as
initial condition $\phi^{\prime}_{\rm ini}\equiv\phi^{\prime}(z_{\rm ini})=0$,
with $z_{\rm ini}=10^{14}$. In this model only the derivatives of the scalar
field enter the equations, so the cosmological evolution of the relevant
quantities is not affected by the initial value of $\phi$. We opt to set
$\phi_{\rm ini}=0$. The interaction starts to play a role around $t_{\rm eq}$,
when the value of $\bar{\rho}_{\rm dm}$ becomes comparable to the one of the
radiation energy density and the source term in the rhs of the KG equation (6)
starts to accelerate the scalar field. In this model the product
$\beta\phi^{\prime}>0$ regardless of the sign of $\beta$, i.e. the rhs of eq.
(5) is negative, so there is a decay of the DM mass that happens mostly during
the matter-dominated era. The model has an almost exact scaling solution in
that epoch, with an EDE fraction $\Omega^{\rm MD}_{\rm ede}\approx
2\beta^{2}/3$ Amendola (2000a)222It would be exact if also baryons were
coupled to the DE. This is a possibility that we do not contemplate in this
paper, just to automatically pass the stringent local constraints on fifth
forces.. Thus, there is a non-zero amount of kinetic EDE during the post-
recombination epoch, even when the potential is a simple constant, as in the
case we are studying here. An increase of the absolute value of the coupling
enhances the fifth force and, for fixed initial conditions, it makes the
universe to be more decelerated and DM to cluster more efficiently in the
matter-dominated epoch. If we consider that baryons are also coupled to DE as
a first approximation, the deceleration parameter and density contrast read,
respectively, $q=-\ddot{a}a/\dot{a}^{2}=\frac{1}{2}+\beta^{2}$ and
$\delta_{\rm m}\sim a^{1+2\beta^{2}}$.

Shorthand | Model
---|---
CDE_const | Coupled DE with constant potential
EXP | Uncoupled DE with exponential potential
CDE_EXP | Coupled DE with exponential potential
ULA | Uncoupled DE with ULA potential
CDE_ULA | Coupled DE with ULA potential
Table 1: Shorthands for the different EDE models studied in this paper.

A larger mass of the DM particles at the last scattering surface could help,
in principle, to alleviate the Hubble tension, since this would decrease the
sound horizon at the baryon-drag epoch, $r_{d}$, which would require in turn a
larger Hubble function at late times to keep fixed the location of the first
peak of the CMB temperature power spectrum measured by Planck. In practice,
though, current data put very tight constraints on $\beta$ and limit strongly
the capability of the model of loosening the $H_{0}$ tension Gómez-Valent  _et
al._ (2020)333Weaker constraints on $\beta$ are obtained when the low
multiples from Planck are combined with ACT and SPT-3G CMB data, which allows
for larger values of $H_{0}$. However, the alleviation of the Hubble tension
is in this case induced by a reduction of the statistical power of the overall
data set Goh  _et al._ (2022).. Nevertheless, past studies found a persistent
peak in the posterior distribution of the coupling under various data sets and
with different statistical significance, which is certainly an interesting
feature of the model, see e.g. Pettorino (2013); Ade  _et al._ (2016); Gómez-
Valent  _et al._ (2020). This peak has recently been shown not to be induced
by volume effects introduced in the marginalization process Gómez-Valent
(2022a).

We deem useful to present in this work the fitting results for the
$\Lambda$CDM and coupled DE with a constant potential and compare them with
those obtained with the coupled and uncoupled EDE models that we describe in
the next two subsections, III.2 and III.3.

###  III.2 Exponential potential

We also study the case of EDE with the following exponential scalar field
potential

$V(\phi)=V_{0}+\frac{\Omega^{\rm RD}_{\rm ede}\rho_{r}(a_{\rm
ini})}{3(1-\Omega^{\rm RD}_{\rm
ede})}\exp\left[\frac{-2\kappa}{\sqrt{\Omega^{\rm RD}_{\rm
ede}}}(\phi-\phi_{\rm ini})\right]\,.$ (20)

During the RDE the model has a scaling solution with an EDE fraction equal to
$\Omega^{\rm RD}_{\rm ede}$ Wetterich (1988); Copeland  _et al._ (1998). This
parameter is left free in our Monte Carlo analyses, with the flat prior
$\Omega^{\rm RD}_{\rm ede}\in[0,1]$. It is the only additional parameter with
respect to the $\Lambda$CDM. The constant factor that is multiplying the
exponential term in formula (20) has been chosen to make the model to be
already in the scaling regime at $z_{\rm ini}$, and $\phi_{\rm ini}$ can be
safely set to $0$, since the dynamics of the scalar field depends only on the
difference $\phi-\phi_{\rm ini}$ 444 Notice that the scaling regime in the MDE
and RDE would be never reached if the constant factor $\lambda$ in the
exponential $e^{-\lambda\kappa\phi}$ was $\lambda<\sqrt{2}$ Copeland  _et al._
(1998); Amendola (2000a). In the latter case the exponential potential would
only become important at low redshifts. Values of $\lambda<\sqrt{2}$ are
required to produce the late-time accelerated phase of the universe if
$V_{0}=0$, of course. Constraints on the model with $V_{0}=0$ and an
exponential potential with $\lambda<\sqrt{2}$ (and other similar potentials)
have been obtained in many previous works in the literature, see again
Amendola (2001); Pettorino  _et al._ (2012); Pettorino (2013); Xia (2013); Ade
_et al._ (2016); van de Bruck  _et al._ (2017); van de Bruck and Mifsud
(2018); Barros  _et al._ (2019); Agrawal  _et al._ (2021); Gómez-Valent  _et
al._ (2020); Gómez-Valent (2022a). We remark that the model CDE_EXP is very
different from the latter, since we automatically have $\lambda>\sqrt{2}$ due
to the physical range of values allowed for $\Omega_{\rm ede}^{\rm RD}$. The
constant term in the potential (20), $V_{0}$, is crucial to ensure the
phenomenological viability of the model..

Figure 1: Diagram in the plane $(\beta,\Omega_{\rm ede}^{\rm RD})$ of CDE_EXP.
The black thick line is the border between regions RGI and RGII. Below it and
in the line $\Omega_{\rm ede}^{\rm RD}=0$ we have region RGI (21), whereas
above it we have RGII (22). The blue part corresponds to the region of
parameter space that leads to a lower growth of matter fluctuations than the
$\Lambda$CDM during the scaling regime in the MDE, $\delta\sim a^{n}$ with
$n<1$, whereas the pink one indicates the region at which $n>1$. The border
between the latter is set by the curves $\Omega_{\rm ede}^{\rm
RD}=4\beta^{2}/(-1+\sqrt{5/2})^{2}+\mathcal{O}(4)$ if $\beta\geq 0$ and
$\Omega_{\rm ede}^{\rm RD}=4\beta^{2}/(1+\sqrt{5/2})^{2}+\mathcal{O}(4)$ if
$\beta<0$. See Sec. III.2 for details.

In the MDE the model has two fixed points. The model is attracted by one or
the other depending on the values of $\Omega^{\rm RD}_{\rm ede}$ and the
coupling $\beta$ Amendola (2000a). These fixed points are not stable due to
the non-zero constant term $V_{0}$ in the potential, which forces the model to
abandon the scaling regime when matter is sufficiently diluted. The universe
enters then the phase of late-time acceleration. Let us analyze first the
phenomenology in the case when there is no coupling in the dark sector.

Exponential potential with no coupling: We have $\Omega^{\rm RD}_{\rm ede}\neq
0$, but $\beta=0$. We call this model “EXP”. It is a particular case of the
double exponential potential firstly studied in Barreiro  _et al._ (2000). The
EDE fraction in the MDE reads in this case $\Omega_{\rm ede}^{\rm
MD}=3\Omega_{\rm ede}^{\rm RD}/4$, i.e. it decreases by $25\%$ from the RDE to
the MDE Wetterich (1988); Copeland  _et al._ (1998). EDE delays the matter-
radiation equality time and slows the growth of matter perturbations down with
respect to the $\Lambda$CDM. During the MDE the matter density contrast at
subhorizon scales grows as $\delta_{m}\sim a^{1-\frac{9}{20}\Omega^{\rm
RD}_{\rm ede}}$. The presence of DE near last scattering affects the relation
between the position of CMB peaks Doran  _et al._ (2001); Doran and Lilley
(2002). The maximum value of $\Omega^{\rm MD}_{\rm ede}$ is actually extremely
constrained by the CMB data, basically due to the stringent bounds imposed on
the EDE fraction at the last scattering surface. We found in Gómez-Valent  _et
al._ (2021) an upper bound $\Omega^{\rm MD}_{\rm ede}\lesssim 0.3\%$ at
$1\sigma$ c.l. using the Planck CMB data and the supernovae of Type Ia (SNIa)
from the Pantheon compilation. This tight constraints render the model
extremely close to the $\Lambda$CDM, in practice, so its impact on the
cosmological tensions is small. Here we also provide constraints on this model
to ease the comparison with the results obtained when we take into account the
effect of the coupling.

Coupled dark energy with an exponential potential: We call this model
“CDE_EXP” along this paper. Here we consider the most general scenario, with
both $\Omega^{\rm RD}_{\rm ede}\neq 0$ and $\beta\neq 0$. The concrete scaling
solution during the MDE is usually reached at $z\sim\mathcal{O}(100)$ and
depends, as mentioned before, on the values of these two parameters Amendola
(2000a).

Figure 2: Upper plot: EDE fraction as a function of the redshift for the
CDE_EXP models listed in the legend; Lower plot: EoS parameter of the
matter+DE fluid. During the MDE $\Omega_{\rm ede}(z)$ and $w_{\rm m\phi}(z)$
tend to the scaling solutions discussed in the main text.

In region

${\rm RGI:}\quad\beta>0\qquad{\rm and}\quad\Omega_{\rm ede}^{\rm
RD}<\frac{16\beta^{2}}{(3+2\beta^{2})^{2}}$ (21)

the EDE fraction and deceleration parameter in the scaling regime read
$\Omega_{\rm ede}^{\rm MD}\approx 2\beta^{2}/3$ and $q\approx 1/2+\beta^{2}$,
respectively, and the equation of state (EoS) of the composite matter+DE fluid
$w_{\rm m\phi}\equiv{\bar{p}}_{\phi}/(\bar{\rho}_{\rm
m}+\bar{\rho}_{\phi})\approx 2\beta^{2}/3$. The DM density contrast
$\delta_{m}\sim a^{1+2\beta^{2}}$ at subhorizon scales, as it happens also in
the model CDE_const described in Sec. III.1. It is important to remark,
though, that in this model the growth of perturbations can be very different
to the one found in CDE_const in previous stages of the cosmic history due to
the effect of the non-trivial EDE potential (20). We will discuss this
explicitly later on. The exponential potential in this model does not play any
role during the scaling regime, since DE only has kinetic energy.

Instead, in region

${\rm RGII:}\quad\beta\leq 0\qquad{\rm or}\quad\Omega_{\rm ede}^{\rm
RD}>\frac{16\beta^{2}}{(3+2\beta^{2})^{2}}\,,$ (22)

we have the following scaling formulas

$\Omega_{\rm ede}^{\rm MD}\approx\frac{3}{4}\Omega_{\rm ede}^{\rm
RD}-\frac{\beta}{2}\sqrt{\Omega_{\rm ede}^{\rm RD}}+\mathcal{O}(4)\,,$ (23)
$q\approx\frac{1}{2}+\frac{3}{4}\beta\sqrt{\Omega_{\rm ede}^{\rm
RD}}+\mathcal{O}(4)\,,$ (24) $w_{\rm
m\phi}\approx\frac{\beta}{2}\sqrt{\Omega_{\rm ede}^{\rm
RD}}+\mathcal{O}(4)\,,$ (25)

where $\mathcal{O}(4)\equiv\mathcal{O}(\beta^{2}\Omega_{\rm ede}^{\rm
RD},\beta(\Omega_{\rm ede}^{\rm RD})^{3/2})$, i.e. we consider
$\beta\sim\mathcal{O}(1)$ and $\sqrt{\Omega_{\rm ede}^{\rm
RD}}\sim\mathcal{O}(1)$. Negative values of $\beta$ produce larger EDE
fractions and a less decelerated universe during the scaling period, contrary
to what we find when $\Omega_{\rm ede}^{\rm RD}=0$, i.e. in CDE_const. In this
epoch the growing mode of the matter density contrast at subhorizon scales
evolves as $\delta_{m}\sim a^{n}$, with

$n=1+\frac{6}{5}\beta^{2}+\frac{6}{5}\beta\sqrt{\Omega_{\rm ede}^{\rm
RD}}-\frac{9}{20}\Omega_{\rm ede}^{\rm RD}+\mathcal{O}(4)\,.$ (26)

Notice that in this case there exists a region in the plane
$(\beta,\Omega_{\rm ede}^{\rm RD})$ that leads to values of $n<1$, i.e. that
lets matter perturbations to grow slower than in the standard model. We call
this region RGIIb to distinguish it from the region with an enhanced growth,
RGIIa. Of course, all these expressions reduce to the ones found under the
condition (21) when $\beta>0$ and $\Omega_{\rm ede}^{\rm RD}\to
16\beta^{2}/9+\mathcal{O}(4)$. Moreover, as expected, we also retrieve the
expressions for the model EXP in the limit $\beta\to 0$. The limit
$\Omega_{\rm ede}^{\rm RD}\to 0$ is more subtle and counterintuitive. One
would expect to recover the scaling solutions of the model CDE_const, but this
only happens if we perform the limit in the region RGI. In RGIIa this limit
leads to different results. In the presence of $\Omega_{\rm ede}^{\rm RD}\neq
0$, regardless how small it is, there is a clear asymmetry between positive
and negative values of the coupling, which is not present when $\Omega_{\rm
ede}^{\rm RD}$ is exactly equal to zero. We find, for instance, that during
the scaling period, for extremely small values of the EDE fraction in the RDE,
$\Omega_{\rm ede}^{\rm MD}$ is negligible if $\beta<0$, whereas it is
proportional to $\beta^{2}$ if $\beta>0$.

Figure 3: Evolution of the mass ratio $g(z)$ (27) in the CDE_EXP model for
four different sets of parameters $(\Omega_{\rm ede}^{\rm RD},\beta)$. The
second one (red dotted line) is located in the region RGI; the others are in
RGII. Notice that for $\Omega_{\rm ede}^{\rm RD}\neq 0$ if $\beta<0$ the DM
mass grows with the expansion, whereas it decreases if $\beta>0$. See the main
text for details. Figure 4: Left plot: Relative differences between the
density contrast of DM computed in the CDE_EXP models specified in the table
and the $\Lambda$CDM (our reference model), as a function of the redshift and
for $k=0.1$ Mpc-1. For the standard model we use the values of the parameters
obtained from the Planck 2018 TTTEEEE analysis Aghanim  _et al._ (2020). In
the inner plot we show the evolution of $\delta_{\rm dm}(z)$ in the
$\Lambda$CDM, to understand when these relative differences lead to important
absolute differences; Right plot: The same, but for the linear matter power
spectrum at $z=0$, in the relevant range of $k$’s.

In Fig. 1 we show a diagram of the $(\beta,\Omega_{\rm ede}^{\rm RD})$-plane,
with the two regions defined by the conditions (21) and (22), and also
indicating with different colors the regions that enhance or suppress the
growth of matter perturbations during the scaling time in the MDE with respect
to the $\Lambda$CDM when the initial conditions and the other parameters are
kept fixed in the analysis. It is clear that the presence of EDE in the RDE
introduces an asymmetry between the solutions found for negative and positive
values of the coupling constant $\beta$, which was not present in CDE_const.
This asymmetry, though, does not generate big changes in the power $n$ (26),
which satisfies $|n-1|<2\%$ for all the points in the region of parameter
space shown in Fig. 1. We will see later that, given some fixed initial
conditions, the largest deviations from the $\Lambda$CDM regarding the growth
of matter perturbations happen during the MDE, but before the scaling period,
i.e. in the redshift range $50\lesssim z\lesssim z_{\rm eq}$, with $z_{\rm
eq}=z(t_{\rm eq})$.

We show in Fig. 2 the evolution of the EDE fraction and the EoS of the
matter+DE fluid for several values of the CDE_EXP parameters.

It is also important to understand how the DM mass evolves in this model,
since it can make the radiation-equality time to happen earlier or later in
the cosmic history, and has an obvious impact on the LSS as well. In contrast
to model CDE_const, in CDE_EXP the mass of the dark matter particles does not
remain constant during the RDE due to the velocity of the scalar field
triggered by the potential, and can also increase throughout the history of
the universe depending on the sign of $\beta$. Using Eq. (7) it is
straightforward to show that

$g(a)\equiv\frac{d}{da}\left(\frac{m_{\rm dm}(a)}{m_{\rm
dm}^{(0)}}\right)=-\frac{m_{\rm dm}}{m^{(0)}_{\rm
dm}}\left[\beta\frac{d(\kappa\phi)}{da}\right]\,.$ (27)

During the RDE $g_{\rm RD}\propto-\beta\sqrt{\Omega_{\rm ede}^{\rm RD}}$, and
during the scaling period in the MDE

$g_{\rm MD}\propto\left\\{\begin{array}[]{ll}-\beta^{2}&\quad{\rm
in\,\,RGI}\\\ -\beta\sqrt{\Omega_{\rm ede}^{\rm RD}}&\quad{\rm
in\,\,RGII}\end{array}\right.$ (28)

Thus, the DM mass in the region RGI of parameter space decreases with the
scale factor, since $\beta>0$. In region RGII, instead, it decreases if
$\beta>0$ and increases if $\beta<0$. We have verified numerically that it is
also the case between the RDE and the scaling period in the MDE. Hence, we
find that if $\Omega_{\rm ede}^{\rm RD}\neq 0$ the DM mass grows with the
expansion if $\beta<0$ and decreases if $\beta>0$. If $\Omega_{\rm ede}^{\rm
RD}=0$ (i.e. in CDE_const) $m_{\rm dm}$ always decreases. In Fig. 3 we plot
the evolution of the DM mass for different sets of the parameters $\Omega_{\rm
ede}^{\rm RD}$ and $\beta$. It is clear that in order to keep the DM mass
around the recombination time $|\Delta m_{\rm dm}/m^{(0)}|\lesssim 5\%$, with
$|\beta|\sim 0.04$, we have to demand $\Omega_{\rm ede}^{\rm
RD}\lesssim\mathcal{O}(1\%)$, which is a similar upper bound to the one
obtained in the uncoupled scenario. Thus, we do not expect the non-null
coupling to relax a lot the existing tight constraints on $\Omega_{\rm
ede}^{\rm RD}$. We will actually see in Sec. VI that this is precisely the
order of magnitude that we get from our fitting analyses.

Now we discuss in more detail the evolution of the DM density contrast in the
CDE_EXP model, fixing the initial conditions and exploring different values of
$\beta$ and $\Omega_{\rm ede}^{\rm RD}$. In the left plot of Fig. 4 we show
the relative differences between the density contrast of DM computed in the
$\Lambda$CDM and the other models, as a function of the redshift. In all cases
the effect of the DE mass can be neglected in good approximation, since the
comoving perturbation scale under study ($k=0.1$ Mpc-1) is larger than the
product $am_{\phi}$ after horizon crossing, see Appendix A. The enhancement or
suppression of the DM density contrast affects the evolution of the baryon
density perturbations through equation (15), and leads also to the increase or
decrease of the total matter power spectrum, see the right plot in Fig. 4.

First of all let us analyze the simplest scenario, in which the coupling is
switched off and there is a non-null EDE fraction (yellow curve). Before the
scale reenters the horizon at $z\sim 10^{6}$ the DM density contrast remains
constant, $1\%$ below the $\Lambda$CDM value due to the $1\%$ additional
radiation content in form of DE. Apart from this, nothing differs from the
standard model, except that during the MDE the matter fraction is lower and DE
is non-zero. This suppresses the perturbations growth.

Figure 5: CMB temperature (TT) and polarization (EE) auto-correlation spectra,
and their cross-correlation (TE), obtained for four different CDE_EXP
scenarios. In the top-right plot we zoom in the region $l\gtrsim 500$ to
better distinguish the differences between the TT spectra in that multipole
range. The black curves correspond to the reference $\Lambda$CDM model, the
same employed in Fig. 4.

What does it happen, instead, when both, $\beta$ and $\Omega_{\rm ede}^{\rm
RD}$, are sizable? Once the scale reenters the horizon, before the matter-
radiation equality time the matter fraction is negligible. The scalar field
has been accelerated during the RDE, though, so it has some velocity and the
friction term in the equation for $\delta_{\rm dm}$ (17) modifies the
evolution of the density contrast, since it is proportional to
$-\beta\phi^{\prime}$. If $\beta<0$ ($\beta>0$) the friction is enhanced
(suppressed). These departures from the $\Lambda$CDM do not have a huge impact
in terms of absolute differences, since they happen when the density
perturbations are still quite small (see the inner plot in Fig. 4). After
$z_{\rm eq}$, the matter fraction in the universe starts to be important and
so also the source term in Eq. (17). The typical variation of $\kappa\phi$
triggered by the potential during the RDE is of order one and positive for the
values employed to build the orange, red, blue and cyan curves ($\Omega_{\rm
ede}^{\rm RD}=10^{-3},10^{-2}$). This statement is easy to prove. During the
RDE the scalar field does not feel the coupling and evolves as

$\kappa\phi(a)=\kappa\phi_{\rm ini}+2\sqrt{\Omega_{\rm ede}^{\rm
RD}}\ln\left(\frac{a}{a_{\rm ini}}\right)\,,$ (29)

Thus,

$\kappa\phi(z_{\rm eq})-\kappa\phi_{\rm ini}\sim 50\sqrt{\Omega_{\rm ede}^{\rm
RD}}\sim\mathcal{O}(1)$ (30)

for the aforementioned range of values of $\Omega_{\rm ede}^{\rm RD}$. This
basically means that at $z_{\rm eq}$ the mass of the DM particles (7) takes
the value

$m_{\rm dm}(z_{\rm eq})=m_{\rm ini}[1-\beta\times\mathcal{O}(1)]\,.$ (31)

Contrary to what we find when there is no important EDE fraction in the RDE,
the source term in Eq. (17) contains a contribution that is proportional to
$\beta$, instead of $\beta^{2}$. For $\beta>0$ the mass is smaller than
$m_{\rm ini}$ and the DM density contrast decreases with respect to the one in
the $\Lambda$CDM. If the coupling is negative the opposite happens. This
behavior can be clearly appreciated in Fig. 4. When the system reaches the
scaling regime, at $z\sim 50$, we can apply the formulas of $\delta_{\rm dm}$
shown before. For the values of $(\beta,\Omega_{\rm ede}^{\rm RD})$ contained
in the region RGIIb (blue curve) of the parameter space there is a suppression
of the DM growth during the scaling period. Nevertheless, as we have already
pointed out before, the largest fraction of the relative differences are
generated between the matter-radiation equality time and the beginning of the
scaling in the MDE.

The brown and green curves in fig. 4 deserve also some comments. Due to the
very small value of $\Omega_{\rm ede}^{\rm RD}$ in these cases, the friction
term plays no role before $t_{\rm eq}$, and this is why we do not observe any
bump at $z>z_{\rm eq}$ in the plot of $\Delta\delta_{\rm dm}/\delta_{\rm dm}$.
The mass of the DM particles at that time is very close to the initial one,
$m_{\rm ini}$, so the source term in Eq. (17) receives only a correction of
order $\beta^{2}$. Hence, regardless of the sign of $\beta$, $\delta_{\rm dm}$
starts to grow after $t_{\rm eq}$. The growth is faster if $\beta>0$, since
the increase of the scalar field and the DM mass is more efficient in this
case.

We illustrate in Fig. 5 the behavior of the CMB temperature and polarization
spectra and their cross-correlation in the the CDE_EXP models. For fixed
initial conditions, a negative coupling makes the mass of the DM particles to
increase with the expansion (see again Fig. 3), which leads to a larger DM
energy density compared to the $\Lambda$CDM. The relative difference is larger
in the late-time universe and this produces an increase of the angle subtended
by the sound horizon at the decoupling time. These phenomena, which are
enhanced by larger values of $\Omega_{\rm ede}^{\rm RD}$, in turn, produce a
decrease of the amplitude of the CMB peaks and their shift towards lower
multipoles. This explains why $\beta$ is positively correlated with $H_{0}$.

###  III.3 ULA potential

We also want to test the performance of the coupled and uncoupled dark energy
models ruled by the ultra-light axion-like potential Poulin  _et al._ (2018,
2019),

$V(\phi)=V_{0}+m^{2}f^{2}[1-\cos(\phi/f)]^{3}\,.$ (32)

This potential introduces three additional parameters with respect to the
$\Lambda$CDM, to wit: $m$, $f$ and the initial value of the scalar field, all
of them with dimensions of energy in natural units. The latter is usually
written in terms of the dimensionless quantity $\theta_{\rm ini}=\phi_{\rm
ini}/f$, with $\theta_{\rm ini}\in[0,\pi]$. Deep in the RDE, and regardless of
the coupling, the expansion rate is much larger than the mass of $\phi$, i.e.
$H\gg m$, so the scalar field has practically no dynamics and $V\simeq const.$
This allows us to take as initial condition $\phi^{\prime}_{\rm ini}=0$.

Figure 6: Upper plots: EDE fractions obtained with different ULA and CDE_ULA
models. We use in all cases the same values of $\omega_{b}$, $V_{0}$,
$\bar{\rho}_{\rm dm}(z_{\rm ini})$, $\log_{10}(f/{\rm eV})$ and
$\log_{10}(m/{\rm eV})$. In the left plot we set $\theta_{\rm ini}=2.9$, and
in the right plot $\theta_{\rm ini}=2.3$. We show the results for
$\beta=0,+0.05,-0.05$; Lower plots: The value of the mass ratio $m_{\rm
dm}(z)/m_{\rm dm}^{(0)}$ for the models studied in the upper plots. Positive
values of $\beta$ lead to a larger $f_{\rm ede}$, but do not necessarily make
the DM mass to decrease with the expansion. For sufficiently low values of
$\theta_{\rm ini}$ a $\beta>0$ produce $m_{\rm dm}^{(0)}>m_{\rm dm}(z_{\rm
ini})$. Notice also that a very significant fraction of the total variation of
the DM mass happens between $z_{\rm max}$ and the decoupling time. See the
main text for more comments.

We discuss the models with and without coupling separately in the subsequent
paragraphs:

ULA with no coupling: In this case $\beta=0$. We refer to this model simply as
“ULA”. The scalar field gets accelerated when $H\sim m$. At that moment it
starts to roll down the potential and eventually oscillates around its
minimum. During the oscillatory phase its energy density decays faster than
radiation, and after the decay the model reduces essentially to the standard
model, with $V(\phi)\approx V_{0}$. If the constant value of the potential
during the RDE is large enough, and if $m\sim H(z_{eq})$, it is possible to
generate a peak in the EDE fraction $f_{\rm ede}$ close to $z_{eq}$, at
$z_{\rm max}$, that can lead to a decrease of the sound horizon and,
therefore, force larger values of the Hubble rate at low redshifts. This is
required to keep the position of the first peak of the CMB temperature angular
power spectrum as it is measured by Planck, and also to keep the good fit to
the BAO data. The decay of the EDE density is sufficiently fast to respect the
tight constraints on the DE fraction around the recombination time Pettorino
_et al._ (2013); Gómez-Valent  _et al._ (2021). ULA was firstly proposed as a
viable option to alleviate the Hubble tension in Poulin  _et al._ (2019).
After that, some studies pointed out that this model tends to enhance the
large-scale structure in the universe, see e.g. Hill  _et al._ (2020); D’Amico
_et al._ (2021). In order to mitigate the increase of the early integrated
Sachs-Wolfe (iSW) effect introduced by EDE the model needs to increase
$\omega_{\rm cdm}$ and $n_{s}$ as well, see e.g. Vagnozzi (2021). This, in
turn, enhances the LSS if the other parameters of the theory (including
$V_{0}$) do not allow for any compensation of this effect. There have been
intensive discussions on this issue in the literature, see e.g. Hill  _et al._
(2020); D’Amico  _et al._ (2021); Murgia  _et al._ (2021); Smith  _et al._
(2021); Gómez-Valent  _et al._ (2021), and Hill  _et al._ (2022); Poulin  _et
al._ (2021) for the results obtained incorporating also the data from the
Atacama Cosmology Telescope Aiola  _et al._ (2020) and SPT-3G 2018 Smith  _et
al._ (2022). These discussions have been mainly based on the interpretation of
the posterior distributions obtained in several fitting analyses for the
quantities $\sigma_{8}$ and $S_{8}=\sigma_{8}(\Omega^{(0)}_{m}/0.3)^{0.5}$,
which tend to peak at larger values than those found for the $\Lambda$CDM and,
therefore, to worsen the $\sigma_{8}/S_{8}$ tension. As firstly argued in
Sánchez (2020), these quantities are sometimes difficult to interpret as clean
LSS estimators, though. The reason is simple. The rms of mass fluctuations at
the scale $R$ reads,

$\sigma_{R}^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}dk\,k^{2}P(k)[W(kR)]^{2}\,,$
(33)

with $R$ the radius of the spherical top-hat window function $W(kR)$ used to
smooth the matter density field. The quantity $\sigma_{8}$ is computed using
$R_{8}=8h^{-1}$ Mpc, which depends on the reduced Hubble parameter
$h=H_{0}/(100\,{\rm km/s/Mpc})$. This means that in each step of a Monte Carlo
analysis the value of $\sigma_{8}$ is not obtained using the same scale, just
because the value of $h$ changes and so does also $R_{8}$. Thus, the fact that
the posterior distribution of $\sigma_{8}$ is shifted towards larger values in
a concrete model with respect to another model does not necessarily imply that
matter fluctuations are larger in the former. If the posterior distribution
for $H_{0}$ also peaks at higher values we systematically evaluate the rms of
mass fluctuations at lower scales and this can produce an enhancement of
$\sigma_{8}$ even if the shape of the matter power spectrum is the same or
very similar in both models.

The author of Sánchez (2020) suggested to use the rms of mass fluctuations in
spheres of radius $R_{12}=12$ Mpc, $\sigma_{12}$, and the related parameter
$S_{12}=\sigma_{12}(\omega_{m}/0.14)^{0.4}$, instead of $\sigma_{8}$ and
$S_{8}$. $R_{12}$ does not depend on $h$ and, hence, $\sigma_{12}$ is free
from the aforementioned bias that might affect some models with posteriors for
$h$ that differ from those found in the $\Lambda$CDM. Fitting values of these
parameters were already obtained for the first time in Gómez-Valent  _et al._
(2021) from the reconstruction of the EDE fraction with various data sets,
showing that the values of $\sigma_{12}/S_{12}$ are in all cases lower than
$\sigma_{8}/S_{8}$ and that the former lie closer to the $\Lambda$CDM results,
which means that these matters can have a direct impact on the quantification
of the tension between the models and the LSS data.

Figure 7: Upper plots: Relative differences between the density contrast of DM
computed in the CDE_ULA models and the $\Lambda$CDM (our reference model), as
a function of the redshift and for $k=0.1$ Mpc-1. For the standard model we
use the values of the parameters obtained from the Planck 2018 TTTEEEE
analysis Aghanim  _et al._ (2020), and for the CDE_ULA models we use the same
values of the parameters as in Fig. 6; Lower plots: Relative differences of
the linear matter power spectrum at $z=0$, using the same models.

In Gómez-Valent (2022a) it was shown that ULA actually alleviates the $H_{0}$
tension while keeping the value of $\sigma_{12}$ extremely close to the ones
obtained in the standard model. This might be pointing out that ULA does not
worsen the tension with the LSS data that much, if at all. In this paper we
will elaborate more on this observation, and will show explicitly that the
best-fit ULA model does not suffer from any enhancement of the parameter
$\sigma_{12}$. The equation that governs the evolution of matter fluctuations
at $k\gg\mathcal{H}$ from the MDE onwards is the same as in the $\Lambda$CDM,
i.e.

$\frac{d^{2}\delta_{m}}{da^{2}}+\frac{3}{2a}\frac{d\delta_{m}}{da}[2-\Omega_{m}(a)]-\frac{3}{2a^{2}}\Omega_{m}(a)\delta_{m}(a)=0\,.$
(34)

During the MDE $\Omega_{m}\approx 1$, so $\delta_{m}\sim a$ both in ULA and
the concordance model. Larger values of $n_{s}$ make matter fluctuations to
cluster more efficiently at large $k^{\prime}s$, though (at $k\gtrsim 0.1$
Mpc-1). Nevertheless, ULA is able to accommodate values of the cosmological
constant much higher than in the $\Lambda$CDM, making $\Omega_{m}$ to remain
close to 0.3 at low redshifts. This helps to decrease the amplitude of $P(k)$
at larger scales (lower $k$’s). $P(k)$ is therefore larger than the one in the
$\Lambda$CDM at small scales, but lower at large scales for typical best-fit
values (cf. Sec. VI and Fig. 12). We will see later that there is a
compensation in the computation of $\sigma_{12}$ for ULA with formula (33)
(which involves an integral over a wide range of $k$’s) that renders its value
stable. This might be pointing out that the assessment of the tension with the
LSS data is more subtle than previously thought, and that it cannot be
certainly carried out by just looking at the posteriors for
$\sigma_{8}/S_{8}$. However, as we have mentioned before, despite the
similarity in the values of $\sigma_{12}$ in ULA and the $\Lambda$CDM, there
exist important differences in the shape of the matter power spectra. Probing
observationally the clustering at different scales and measuring unambiguous
LSS estimators (see e.g. Semenaite  _et al._ (2022)) can be a good way of
arbitrating this issue. We will come back to this discussion in Sec. VI.

Coupled dark energy with an ULA potential: We call the ULA model with a non-
null coupling (i.e. with $\beta\neq 0$) “CDE_ULA”. We suggested in the
conclusions of our previous work Gómez-Valent  _et al._ (2021) that a coupling
in the dark sector could help, in principle, some EDE models with a non-
negligible fraction of DE in some windows of the pre-recombination epoch to
suppress the growth of matter fluctuations, while keeping large values of
$H_{0}$. Some authors have already worked in this direction Karwal  _et al._
(2022); McDonough  _et al._ (2022). In this paper we also explore this
possibility for coupled ULA, trying to understand more in detail the
phenomenology of the model and taking special care in the interpretation of
the fitting results by providing constraints on some parameters that can help
us to better assess the tensions, and also by studying the impact of volume
effects in the marginalization process, which have been already found to be
relevant in the uncoupled ULA model Herold  _et al._ (2022); Gómez-Valent
(2022a).

Let us begin analyzing how the coupling affects the evolution of the EDE
fraction and the dark matter mass in CDE_ULA. If $z_{\max}\sim z_{\rm eq}$,
the scalar field feels first the coupling, rather than the curvature of the
potential. This is why positive values of $\beta$ allow to climb the latter up
before rolling it down. This has two effects: (i) the moment at which $m\sim
H$ happens at later times (i.e. lower redshifts) due to the increase of
$\bar{\rho}_{\rm\phi}$; and (ii) $f_{\rm ede}$ is larger. This can be seen in
the two upper plots of Fig. 6, for two different values of $\theta_{\rm ini}$.
If $z_{\max}\gtrsim z_{\rm eq}$ these effects are also observed, but are less
prominent and eventually become negligible when $z_{\max}\gg z_{\rm eq}$, just
because the coupling does not have time to act alone, without the influence of
the potential.

Deep in the RDE the scalar field is in good approximation frozen and the DM
mass remains almost constant, but the coupling injects some dynamics before
$z_{\rm max}$ if $z_{\rm max}\sim z_{\rm eq}$. In this phase of the cosmic
expansion if $\beta>0$ the scalar field grows mildly, whereas if the coupling
is negative $\phi$ decreases. When the scalar field dynamics starts to be
dominated by the potential, close to $z_{\rm max}$, we find that negative
$\beta$’s produce a decrease of $m_{\rm dm}$ due to Eq. (27) and the fact that
$d\phi/da<0$. If $\beta>0$, instead, the DM mass can increase or decrease
depending on the value of the scalar field at $z_{\rm max}$. If $\phi(z_{\rm
max})>\pi f$ the DM mass decreases because $d\phi/da>0$, see again Eq. (27).
If $\phi(z_{\rm max})<\pi f$ the DM mass increases because the opposite
happens. This is why the initial condition of the scalar field plays also a
role, as it is clear from the lower plots in Fig. 6. If $\theta_{\rm ini}$ is
close enough to $\pi$ and $\beta>0$ the scalar field can oscillate around
$\phi=2\pi f$ after the decay instead of around $\phi=0$. If $z_{\rm max}\gg
z_{\rm eq}$ this possibility is much less probable and requires either a fine
tuning of $\theta_{\rm ini}$ to $\pi$ or larger (positive) values of the
coupling to happen. Notice that the largest fraction of the change of $m_{\rm
dm}$ happens right after $z_{\max}$, during the rolling-down period. The mean
value of the DM mass is kept constant in the oscillatory phase.

| Planck18 | Planck18+SNIa+BAO
---|---|---
Parameter |  $\Lambda$CDM | CDE_const |  $\Lambda$CDM | CDE_const
$10^{2}\omega_{b}$ |  $2.239\pm 0.015$ (2.240) |  $2.235^{+0.014}_{-0.015}$ (2.229) |  $2.247^{+0.013}_{-0.014}$ (2.252) |  $2.239^{+0.014}_{-0.015}$ (2.242)
$\omega_{\rm cdm}$ |  $0.1203^{+0.0011}_{-0.0013}$ (0.1209) |  $0.1187^{+0.0028}_{-0.0012}$ (0.1194) |  $0.1190\pm 0.0008$ (0.1189) |  $0.1187\pm 0.0008$ (0.1180)
$n_{s}$ |  $0.966\pm 0.004$ (0.966) |  $0.965^{+0.004}_{-0.005}$ (0.961) |  $0.968\pm 0.004$ (0.966) |  $0.966\pm 0.004$ (0.964)
$\tau_{\rm reio}$ |  $0.056^{+0.007}_{=0.006}$ (0.054) |  $0.055^{+0.007}_{-0.008}$ (0.055) |  $0.058^{+0.007}_{-0.008}$ (0.061) |  $0.057^{+0.007}_{-0.008}$ (0.054)
$\sigma_{12}$ |  $0.807\pm 0.008$ (0.809) |  $0.811^{+0.009}_{-0.010}$ (0.819) |  $0.799\pm 0.007$ (0.798) |  $0.808^{+0.009}_{-0.010}$ (0.808)
$H_{0}$ |  $67.47^{+0.56}_{-0.48}$ (67.29) |  $68.55^{+0.56}_{-1.73}$ (68.31) |  $68.03^{+0.35}_{-0.36}$ (68.15) |  $68.43^{+0.43}_{-0.53}$ (68.95)
$\beta$ | $-$ |  $<0.040$ (0.040) | $-$ |  $0.028^{+0.017}_{-0.013}$ (0.041)
$r_{d}$ |  $146.93^{+0.25}_{-0.27}$ (146.75) |  $146.88^{+0.26}_{-0.28}$ (146.69) |  $147.18^{+0.20}_{-0.22}$ (147.15) |  $146.99^{+0.27}_{-0.24}$ (146.89)
$M$ | $-$ | $-$ |  $-19.408\pm 0.010$ (-19.404) |  $-19.397^{+0.012}_{-0.015}$ (-19.381)
$S_{8}$ |  $0.832^{+0.012}_{-0.013}$ (0.836) |  $0.827^{+0.015}_{-0.013}$ (0.838) |  $0.818^{+0.010}_{-0.009}$ (0.816) |  $0.824^{+0.010}_{-0.011}$ (0.821)
$\sigma_{8}$ |  $0.814\pm 0.005$ (0.814) |  $0.827^{+0.007}_{-0.021}$ (0.833) |  $0.811\pm 0.006$ (0.810) |  $0.822^{+0.009}_{-0.013}$ (0.827)
$S_{12}$ |  $0.813^{+0.010}_{-0.009}$ (0.817) |  $0.814\pm 0.010$ (0.823) |  $0.803\pm 0.008$ (0.801) |  $0.810^{+0.009}_{-0.010}$ (0.809)
$\Delta\chi^{2}_{\rm min}$ | $-$ | 1.93 | $-$ | 2.39
$\Delta{\rm AIC}$ | $-$ | -0.07 | $-$ | 0.39
Table 2: Mean and $1\sigma$ uncertainties of the individual parameters of the
$\Lambda$CDM and CDE_const models, together with the corresponding best-fit
values (between parenthesis). $H_{0}$ and $r_{d}$ are expressed in km/s/Mpc
and Mpc, respectively. In the last two rows we report the differences
$\Delta\chi^{2}_{\rm min}\equiv\chi^{2}_{\Lambda,{\rm min}}-\chi^{2}_{{\rm
CDE\\_const},{\rm min}}$ and $\Delta{\rm AIC}\equiv{\rm AIC}_{\Lambda}-{\rm
AIC}_{\rm CDE\\_const}$ to better assess the fitting performance of the
aforementioned models. We discuss these results in Sec. VI.

In the upper plots of Fig. 7 we show how the wave mode $k=0.1$ Mpc-1 of the DM
density contrast evolves in the ULA and CDE_ULA models with respect to the
$\Lambda$CDM, for fixed initial conditions and two alternative values of
$\theta_{\rm ini}$. As expected, the larger is $f_{\rm ede}$ the larger is the
relative decrease of $\delta_{\rm dm}$ at $\sim z_{\rm max}$, since dark
energy fights against the aggregation of matter. After the decay of the scalar
field, in ULA the relative difference with respect to the $\Lambda$CDM result
remains constant, since the equation for the density contrast (34) is the same
to the one in the standard model. In CDE_ULA, though, the fifth force makes
$\delta_{\rm dm}$ to grow faster in the MDE. This growth can compensate in
most of the cases the negative bump generated around $z_{\rm max}$, giving
rise to an increase of the DM density contrast at low redshifts, but there are
also cases in which this compensation is not complete (see e.g. the green
curve in the left upper plot of Fig. 7). The oscillations of the scalar field
deep in the MDE (see the inner plots in Fig. 6) make the friction term of Eq.
(16) to be inefficient, so in CDE_ULA the equation of the DM density constrast
reads in good approximation,

$\delta_{\rm dm}^{\prime\prime}+\mathcal{H}\delta^{\prime}_{\rm
dm}-\frac{\kappa^{2}a^{2}}{2}\left[\bar{\rho}_{b}\delta_{b}+\bar{\rho}_{\rm
dm}\delta_{\rm dm}(1+2\beta^{2})\right]=0\,.$ (35)

During the scaling regime in the MDE, when the matter fraction oscillates
around $\Omega_{m}\sim 1-2\beta^{2}/3$, $\delta_{\rm dm}\sim
a^{1+6\beta^{2}/5}$ McDonough  _et al._ (2022), so matter fluctuations grow in
this epoch slower than in CDE_const, but faster than in the concordance model.
The relative differences between the matter power spectrum of the CDE_ULA
models and the $\Lambda$CDM are shown in the lower plots of Fig. 7. We remark
that these plots have been obtained fixing the initial conditions. By changing
them, and the parameters of the EDE scalar field potential we can make the
density perturbations to evolve differently, of course. Here we only wanted to
show what is the impact of the coupling. For a precise assessment of the
ability of these models to describe the current cosmological data and loosen
the cosmological tensions we are forced to perform a detailed fitting
analysis. We devote the next sections to explain the fitting strategy (data
sets and methodology) and report our results.

Figure 8: One-dimensional posterior distributions for some parameters of the
$\Lambda$CDM and CDE_const, and the corresponding contour plots at $1\sigma$
and $2\sigma$ c.l. The constraints from neutrino experiments (cf. Sec. IV) do
not allow $\sum m_{\nu}$ to be lower than $\sim 0.06$ eV. The green vertical
bands in the first column correspond to the measurement of $H_{0}$ by the
SH0ES team Riess  _et al._ (2022), whereas the gray horizontal bands are the
weak lensing constrain on $S_{8}$ from KiDS+VIKING-450+DES-Y1, which is
obtained under the assumption of the $\Lambda$CDM Joudaki  _et al._ (2020).
Both are provided at $1\sigma$ and $2\sigma$ c.l. See the comments in Sec. VI.

##  IV Data sets

These are the data sets employed in our fitting analyses:

CMB data: We employ the full Planck 2018 TTTEEE+lowE+lensing likelihood
Aghanim  _et al._ (2020), i.e. the data on the temperature (TT) and
polarization (EE) anisotropies of the cosmic microwave background (CMB), their
cross-correlations (TE), and the data on the CMB lensing reconstruction. We
vary in our Monte Carlo (MC) runs the 21 Planck nuisance parameters together
with the cosmological ones. We denote this data set as Planck18, in short.

Supernovae and BAO: We also fit our models to a richer data set that
incorporates the observational ingredients usually used to build the cosmic
inverse distance ladder Aubourg  _et al._ (2015); Cuesta  _et al._ (2015);
Feeney  _et al._ (2019); Camarena and Marra (2020a), which is relevant for the
discussion on the $H_{0}$ tension. We combine Planck18 with the data on SNIa
from the Pantheon compilation Scolnic  _et al._ (2018) and the BAO information
from several galaxy surveys. We call this data set Planck18+SNIa+BAO. The
absolute magnitude of the SNIa, $M$, is left free in the MC runs. We display
its value in our tables in all cases in which we employ SNIa. It is a relevant
quantity. The SH0ES collaboration has measured $M$ using the calibration of
the SNIa in the first steps of the (direct) cosmic distance ladder, $M_{\rm
SH0ES}=-19.253\pm 0.027$ Riess  _et al._ (2022), which is independent from the
underlying cosmology. This value is fully compatible with the one inferred in
model-independent studies of low-redshift cosmological data sets Gómez-Valent
(2022b); Benisty  _et al._ (2022), but much higher than the one obtained in
cosmological analyses of the $\Lambda$CDM using the CMB data from Planck and
high-redshift SNIa. The latter usually lies close to $M\sim-19.40$. The
$H_{0}$ tension can be thought of as a tension between the locally measured
value of $M$ and the one inferred from cosmological studies in the context of
the standard model, see e.g. Camarena and Marra (2020b); Efstathiou (2021). It
can be therefore enlightening to discuss the $H_{0}$ tension directly in terms
of $M$. We do so in Sec. VI.

| Planck18 | Planck18+SNIa+BAO
---|---|---
Parameter | EXP | CDE_EXP | EXP | CDE_EXP
$10^{2}\omega_{b}$ |  $2.236\pm 0.016$ (2.235) |  $2.256^{+0.020}_{-0.019}$ (2.247) |  $2.248^{+0.014}_{-0.013}$ (2.250) |  $2.247^{+0.016}_{-0.026}$ (2.256)
$\omega_{\rm cdm}$ |  $0.1210^{+0.0015}_{-0.0015}$ (0.1201) |  $0.1229^{+0.0023}_{-0.0030}$ (0.1211) |  $0.1191^{+0.0008}_{-0.0009}$ (0.1192) |  $0.1189^{+0.0009}_{-0.0009}$ (0.1195)
$n_{s}$ |  $0.964\pm 0.004$ (0.965) |  $0.964\pm 0.005$ (0.965) |  $0.968\pm 0.004$ (0.969) |  $0.967\pm 0.004$ (0.968)
$\tau_{\rm reio}$ |  $0.056^{+0.007}_{-0.008}$ (0.056) |  $0.051\pm 0.008$ (0.042) |  $0.060^{+0.007}_{-0.008}$ (0.058) |  $0.055^{+0.010}_{-0.007}$ (0.051)
$\sigma_{12}$ |  $0.805^{+0.009}_{-0.008}$ (0.803) |  $0.814^{+0.011}_{-0.015}$ (0.821) |  $0.797\pm 0.007$ (0.802) |  $0.811^{+0.010}_{-0.015}$ (0.824)
$H_{0}$ |  $67.20^{+0.64}_{-0.63}$ (67.45) |  $65.43^{+1.95}_{-1.48}$ (66.54) |  $68.01\pm 0.37$ (67.97) |  $68.28^{+0.58}_{-0.70}$ (67.72)
$\Omega_{\rm ede}^{\rm RD}$ [%] |  $<0.22$ (0.06) |  $<0.29$ (0.03) |  $<0.13$ (0.00) |  $<0.13$ (0.00)
$\beta$ | $-$ |  $-0.061^{+0.018}_{-0.038}$ (-0.079) | $-$ |  $-0.008^{+0.090}_{-0.115}$ (-0.055)
$r_{d}$ |  $146.65^{+0.42}_{-0.28}$ (146.97) |  $146.83^{+0.41}_{-0.35}$ (147.09) |  $147.07^{+0.24}_{-0.22}$ (147.10) |  $146.99^{+0.36}_{-0.48}$ (147.12)
$M$ | $-$ | $-$ |  $-19.409\pm 0.010$ (-19.405) |  $-19.401^{+0.016}_{-0.019}$ (-19.417)
$S_{8}$ |  $0.833\pm 0.014$ (0.827) |  $0.855\pm 0.021$ (0.852) |  $0.817\pm 0.010$ (0.821) |  $0.829^{+0.011}_{-0.017}$ (0.846)
$\sigma_{8}$ |  $0.810^{+0.008}_{-0.006}$ (0.809) |  $0.803^{+0.019}_{-0.021}$ (0.820) |  $0.808\pm 0.006$ (0.812) |  $0.825^{+0.009}_{-0.015}$ (0.833)
$S_{12}$ |  $0.813\pm 0.011$ (0.809) |  $0.826^{+0.014}_{-0.016}$ (0.829) |  $0.801\pm 0.008$ (0.805) |  $0.814^{+0.010}_{-0.014}$ (0.829)
$\Delta\chi^{2}_{\rm min}$ | 2.75 | 7.51 | -0.15 | 3.75
$\Delta{\rm AIC}$ | 0.75 | 3.51 | -2.15 | -0.25
Table 3: As in Table 2, but for the EXP and CDE_EXP models.

Our BAO data set includes the following points:

  * •

$D_{V}/r_{d}$ at $z=0.122$ provided in Carter  _et al._ (2018), with $D_{V}$
being the dilation scale,

$D_{V}(z)=\left[D_{M}^{2}(z)\frac{cz}{H(z)}\right]^{1/3}\,,$ (36)

and $D_{M}=(1+z)D_{A}(z)$ the comoving angular diameter distance. This data
point combines the dilation scales previously reported by the 6dF Galaxy
Survey (6dFGS) Beutler  _et al._ (2011) at $z=0.106$ and the Sloan Digital Sky
Survey (SDSS) Main Galaxy Sample at $z=0.15$ Ross  _et al._ (2015).

  * •

The anisotropic BAO data ($D_{A}(z)/r_{d}$, $H(z)r_{d}$) measured by BOSS
using the LOWZ ($z=0.32$) and CMASS ($z=0.57$) galaxy samples Gil-Marín  _et
al._ (2017).

  * •

The dilation scale measurements by WiggleZ at $z=0.44,0.60,0.73$ Kazin  _et
al._ (2014).

  * •

$D_{A}(z)/r_{d}$ at $z=0.81$ measured by the Dark Energy Survey (DES) Abbott
_et al._ (2019).

  * •

The anisotropic BAO data from the extended BOSS Data Release 16 (DR16) quasar
sample at $z=1.48$ Neveux  _et al._ (2020).

  * •

The anisotropic BAO information obtained from the Ly$\alpha$ absorption and
quasars of the final data release (SDSS DR16) of eBOSS, at $z=2.334$ du Mas
des Bourboux  _et al._ (2020).

A comment on BBN: The only models among those described in Sec. III that could
modify significantly the big bang nucleosynthesis (BBN) processes are those
with the exponential potential discussed in Sec. III.2, i.e. EXP and CDE_EXP,
since we assume that the scaling solution in these models is already reached
at $z_{\rm ini}=10^{14}\gg z_{\rm BBN}$. Nevertheless, we will see in Sec. VI
that we get very tight constraints on $\Omega_{\rm ede}^{\rm RD}$ from the CMB
data, i.e. using CMB data alone we obtain a way lower upper bounds than those
imposed by the BBN data Uzan (2011), so even in these models EDE does not have
a sizable impact on the physics at the BBN epoch. On the other hand, the upper
value of the flat prior employed for $z_{\rm max}$ in ULA and CDE_ULA (cf.
Sec. III.3), is five orders of magnitude lower than $z_{\rm BBN}$. Hence, the
EDE fraction is completely negligible at the BBN time and therefore does not
affect the nucleosynthesis physics either. Some effects can be introduced by
shifts of the baryon energy density in this model, but they are in any case
very small Seto and Toda (2021). In summary, we do not need to consider any
BBN constraint in our study.

Massive neutrinos: In most of our fitting analyses we use two massless and one
massive neutrino of $0.06$ eV. However, we also constrain some models using
the Planck18+SNIa+BAO data set allowing the three neutrino masses to vary
simultaneously in the Monte Carlo runs, assuming a normal ordering, i.e.
$m_{1}<m_{2}<m_{3}$. The latter is a sufficient assumption because the
constraints on the sum of the neutrino masses $\sum m_{\nu}$ are not very
sensitive to the neutrino hierarchy Loureiro  _et al._ (2019). We consider the
constraints from solar and atmospheric neutrinos, which read, respectively
Esteban  _et al._ (2017): $\Delta m_{21}^{2}\equiv
m_{2}^{2}-m_{1}^{2}=(7.50^{+0.19}_{-0.17})\times 10^{-5}\,{\rm eV^{2}}$ and
$\Delta m_{31}^{2}\equiv|m_{3}^{2}-m_{1}^{2}|=(2.524^{+0.039}_{-0.040})\times
10^{-3}\,{\rm eV^{2}}$. We actually vary $m_{1}$, $\Delta m^{2}_{21}$ and
$\Delta m^{2}_{31}$ in the Monte Carlo, and compute the masses of the two
heavier neutrinos using

Figure 9: Triangle plot for the EXP model, with the parameters that are
relevant for the discussion of the cosmological tensions. As in Fig. 8 and the
other triangle plots in this paper, we show the constraint from SH0ES Riess
_et al._ (2022) in green and the one from KiDS+VIKING-450+DES-Y1 Joudaki  _et
al._ (2020) in gray.

$m_{2}=\sqrt{m_{1}^{2}+|\Delta m^{2}_{21}|}\quad;\quad
m_{3}=\sqrt{m_{1}^{2}+|\Delta m^{2}_{31}|}$ (37)

Massive neutrinos reduce the clustering of matter below the free streaming
scale (see e.g. Lesgourgues and Pastor (2006)) and, therefore, could have an
impact on our analyses if the data give room to larger values of $\sum
m_{\nu}$.

We display the fitting results for the models discussed in Sec. III under the
Planck18 and Planck18+SNIa+BAO data sets in Sec. VI.

Figure 10: Contour plots at $1\sigma$ and $2\sigma$ c.l. and one-dimensional
posterior distributions for the CDE_EXP model, obtained with the Planck18 and
Planck18+SNIa+BAO data sets.

##  V Method

We have implemented the background and linear perturbations equations of Sec.
II and the various scalar field potentials presented and discussed in Sec. III
in our own modified versions of the Einstein-Boltzmann system solver CLASS
Blas  _et al._ (2011). All of them have gone, of course, through the
corresponding validation processes in order to check their correct
functioning. For instance, we have checked that we recover the expected
scaling solutions in CDE_EXP and their nested models CDE_const and EXP; and
for ULA we have checked that we are able to reproduce other results found in
the literature.

We have used MontePython Audren  _et al._ (2013) to sample the likelihood with
the Metropolis-Hastings algorithm Metropolis  _et al._ (1953); Hastings
(1970). It is constructed from the data sets described in Sec. IV. The Markov
chains obtained in the output have been processed also with MontePython to
obtain the marginalized constraints for the individual parameters and the
absolute best-fit values, and with GetDist Lewis (2019) to generate the
confidence contour lines in all the relevant planes of parameter space and the
marginalized one-dimensional posteriors. The minimum values of the $\chi^{2}$
in each run, i.e. $\chi^{2}_{\rm min}$, and the corresponding best-fit values
of the parameters have been found applying the procedure described in Appendix
D.1 of Schöneberg  _et al._ (2022). We report our results in the tables and
figures of the next section.

For some parameters of CDE_EXP, ULA and CDE_ULA we have also applied the
profile distribution method to get rid of the volume effects that could affect
them, see e.g. Trotta (2017). We comment on these results in Sec. VI. For
previous analyses using this statistical technique see Ade  _et al._ (2014);
Herold  _et al._ (2022); Gómez-Valent (2022a); Campeti and Komatsu (2022). We
use the same approach explained in Gómez-Valent (2022a), i.e. we compute the
individual profile distributions directly from the Markov chains, instead of
carrying out a series of minimization runs. This allows us to save a
considerable amount of computing time.

In order to minimize the computational time spent in our Monte Carlo runs, we
have tried to modify CLASS avoiding the use of the shooting method. This is
impossible, though, in the case of ULA with and without coupling, since in
this case we need to sample the distribution using the pair of parameters
$(f_{\rm ede},z_{\rm max})$ instead of $(m,f)$, cf. Sec. III.3. This is more
reasonable from the cosmological perspective, and not doing so could hide the
evidence for a non-null EDE fraction $f_{\rm ede}$ due to prior issues Hill
_et al._ (2020). We explain the procedure that we have followed to implement
the shooting method for these models in Appendix B.

##  VI Results and discussion

We report our main fitting results in Tables 2-5 and Figs. 8-15. In Appendix C
we present some additional tables, where we list the contribution of the
$\chi^{2}_{i}$’s of the individual data sets to the minimum value
$\chi^{2}_{\rm min}$ in each model.

Figure 11: Confidence regions at $1\sigma$ and $2\sigma$ c.l. and one-
dimensional posterior distributions obtained with the Planck18+SNIa+BAO data
set with the $\Lambda$CDM and ULA. The parameters $f_{\rm ede}$, $z_{\rm max}$
and $\theta_{\rm ini}=\phi_{\rm ini}/f$ are the maximum EDE fraction, the
redshift at which there is the peak and the initial value of the scalar field
at $z_{\rm ini}=10^{14}$, respectively (cf. Sec. III). For ULA+$\sum m_{\nu}$
we allow the sum of the neutrino masses to vary in the Monte Carlo analysis,
using the approach described in Sec. IV.

We discuss first the output of the analysis of the $\Lambda$CDM and CDE_const.
These two models are nested. The former is described in terms of six
cosmological parameters, whereas the latter has an additional parameter, the
coupling $\beta$. It controls the departure of CDE_const from the
$\Lambda$CDM. We find that the CMB data from Planck does not exclude (when
used alone) values of $\beta\sim 0.07$ at $2\sigma$ c.l., but the value of the
coupling is still fully compatible with 0, being $\beta<0.04$ at $1\sigma$
c.l. (see the second column of Table 2). In Fig. 8 we observe a positive
correlation between the coupling and the Hubble parameter, as expected
Pettorino  _et al._ (2013); Ade  _et al._ (2016); Gómez-Valent  _et al._
(2020). Also between $\beta$ and the LSS estimator $\sigma_{12}$. The reasons
have been already discussed in Sec. III. A larger coupling increases the
growth of matter perturbations during the MDE and requires a larger expansion
rate at low redshifts to not alter the location of the CMB peaks and the
description of the BAO data. The contours in the $(H_{0},\beta)$-plane show
that it is possible to obtain values of $H_{0}$ compatible with the distance
ladder measurement from SH0ES Riess  _et al._ (2022) at $2\sigma$ c.l., due
the long right tale of the distribution. There is, though, a shift towards
larger values of $\sigma_{12}$, which nevertheless is not drastically big.
When we add the SNIa+BAO data, we can see that the values of $\beta\gtrsim
0.06$ get strongly excluded. This, in turn, removes the tail of the one-
dimensional posterior of $H_{0}$, limiting its value to
$H_{0}=(68.43^{+0.43}_{-0.53})$ km/s/Mpc at $1\sigma$ c.l. and keeping the
$H_{0}$ tension at more than $3\sigma$. Although the difference between
$\Delta\chi^{2}=\chi^{2}_{\Lambda,{\rm min}}-\chi^{2}_{\rm
CDE\\_const,min}=2.39$ is small and certainly there is no significant
statistical preference for an interaction in the dark sector, it is worth to
point out the existence of a peak in the posterior of $\beta$ when we include
the full data set. The maximum is located $\sim 2\sigma$ away from 0, and the
mean reads $\beta=0.028^{+0.017}_{-0.013}$ at 68.2% c.l. This peak has also
been found in past analyses using similar data configurations, see Pettorino
_et al._ (2013); Ade  _et al._ (2016); Gómez-Valent  _et al._ (2020), and it
was recently shown in Gómez-Valent (2022a) that it is not sourced by spurious
marginalization effects. It allows a better description of the CMB data (cf.
Table 7).

By allowing the mass of the neutrino masses to vary in the analysis we find
very similar results, see Table 5. The peak is still present, and slightly
displaced to the right ($\beta=0.033\pm 0.015$) due to the positive
correlation between the coupling and $\sum m_{\nu}$, which allows to
compensate the enhancement of the matter power spectrum caused by the fifth-
force. This positive correlation, which is more evident in the range
$\beta>0.03$, also explains the weaker constraints on $\sum m_{\nu}$, which in
CDE_const is $<0.18$ eV at $95\%$ c.l., whereas in the $\Lambda$CDM $<0.12$
eV. Lower values of the parameter $\sigma_{12}$ are allowed due to the
suppression of power at low scales introduced by the massive neutrinos, but
the differences are small. The contours in the $(H_{0},\beta)$-plane of Fig. 8
show, as expected, that the large values of $\beta$ are less able to lead to
larger values of the Hubble parameter than in the case in which we consider a
minimal neutrino setup with only one massive neutrino of 0.06 eV.

We use the Akaike information criterion (AIC) Akaike (1974) to penalize the
additional degrees of freedom and perform a fairer comparison of the
statistical performance of the models. For a sufficiently large number of data
points, the AIC is defined as

${\rm AIC}=\chi^{2}_{\rm min}+2{\rm k},$ (38)

with k the number of parameters of the model. The difference between the AIC
values obtained in two models (1 and 2) can be used as an approximation of the
logarithm of the Bayes ratio $B_{12}$, ${\rm AIC}_{12}={\rm AIC}_{2}-{\rm
AIC}_{1}\approx 2\ln{B_{12}}$, with $B_{12}=\mathcal{E}_{2}/\mathcal{E}_{1}$
and $\mathcal{E}_{i}$ the Bayesian evidence of model $i$. The small
differences $\Delta{\rm AIC}\lesssim\mathcal{O}(1)$ between the $\Lambda$CDM
and CDE_const confirm that from the Bayesian point of view the two models
perform similarly well in the description of the Planck18 and
Planck18+SNIa+BAO data sets. There is a slight decrease of the $\chi^{2}_{\rm
min}$ in CDE_const, but the addition of the coupling $\beta$ is not especially
favored in the light of Occam’s razor.

EDE with an exponential potential (20) is already very constrained by the CMB
data, which limits $\Omega_{\rm ede}^{\rm RD}<0.22\%$ at 1$\sigma$ c.l.,
basically due to the tight upper bound that the EDE fraction has to satisfy at
the decoupling time Gómez-Valent  _et al._ (2021). Our fitting results for
this model are reported in Table 3. There is no room in this model for an
alleviation of the cosmological tensions. The addition of SNIa+BAO gives even
less margin, since $\Omega_{\rm ede}^{\rm RD}<0.13\%$ at 1$\sigma$ c.l. In
Fig. 9 we can see that the constraints on $\sigma_{12}$ and $H_{0}$ are
actually very similar to those found in the $\Lambda$CDM, with a negative
correlation between these two parameters and with a peak of the posterior
lying very far away from the SH0ES bands. We have checked that volume effects
are unimportant in this model.

CDE_EXP allows us to explain much better than the $\Lambda$CDM the CMB data.
The corresponding $\chi^{2}_{\rm min}$ is 7.51 units below the values obtained
in the standard model, and this is thanks to a $2\sigma$ departure of the
coupling from zero, with the posterior of $\beta$ peaking in the negative
region. Notice that due to the non-trivial potential (20), the coupling can
acquire larger values than in CDE_const. The regions RGIIa and RGI of
parameter space are preferred (see Figs. 1 and 10), being the former more
probable when we allow for larger values of the EDE fraction in the RDE. This
is supported also by the AIC. We obtain $\Delta{\rm AIC}=+3.51$, which
according to Jeffreys’ scale corresponds to a moderate evidence for CDE_EXP.
However, this hint for new physics is washed out by the SNIa and BAO data, cf.
Table 3. In Fig. 10 we cannot appreciate any decrease of $\sigma_{12}$
compared to the $\Lambda$CDM, although the model allows for larger values of
$H_{0}$. The tension is only mildly loosened, though, at a similar level at
which it is in CDE_const. Under the Planck18+SNIa+BAO data set we find a peak
at $\beta>0$. This data set seems to prefer a decrease of the DM mass with the
expansion. The constraints on $\Omega_{\rm ede}^{\rm RD}$ are similar to those
obtained in the absence of coupling.

| Planck18+SNIa+BAO
---|---
Parameter | ULA | CDE_ULA
$10^{2}\omega_{b}$ |  $2.265^{+0.016}_{-0.022}$ (2.269) |  $2.271^{+0.016}_{-0.022}$ (2.290)
$\omega_{\rm cdm}$ |  $0.1225^{+0.0016}_{-0.0041}$ (0.1237) |  $0.1211^{+0.0011}_{-0.0029}$ (0.1201)
$n_{s}$ |  $0.974^{+0.005}_{-0.009}$ (0.978) |  $0.974^{+0.005}_{-0.007}$ (0.978)
$\tau_{\rm reio}$ |  $0.059\pm 0.007$ (0.060) |  $0.054^{+0.009}_{-0.010}$ (0.049)
$\sigma_{12}$ |  $0.799\pm 0.007$ (0.796) |  $0.822^{+0.018}_{-0.021}$ (0.839)
$H_{0}$ |  $69.08^{+0.56}_{-1.20}$ (69.75) |  $68.99^{+0.46}_{-0.81}$ (69.24)
$f_{\rm ede}$ [%] |  $<4.80<9.46$ (5.17) |  $<2.24<7.70$ (1.13)
$z_{max}$ |  $7806^{+1000}_{-6000}$ (3874) |  $5317^{+1600}_{-2500}$ (4724)
$\theta_{ini}$ |  $>2.01$ (2.63) |  $>1.39$ (1.83)
$\beta$ | $-$ |  $-0.017^{+0.110}_{-0.083}$ (-0.089)
${\rm log_{10}(}m{\rm/eV)}$ |  $-26.66^{+0.34}_{-0.79}$ (-27.30) |  $-26.86^{+0.34}_{-0.60}$ (-27.10)
${\rm log_{10}(}f{\rm/eV)}$ |  $26.51^{+0.21}_{-0.39}$ (26.51) |  $26.52^{+0.42}_{-0.51}$ (26.63)
$r_{d}$ |  $145.26^{+2.04}_{-0.78}$ (144.37) |  $145.72^{+1.51}_{-0.40}$ (145.89)
$M$ |  $-19.376^{+0.017}_{-0.037}$ (-19.355) |  $-19.380^{+0.012}_{-0.025}$ (-19.373)
$S_{8}$ |  $0.825^{+0.011}_{-0.012}$ (0.823) |  $0.845^{+0.019}_{-0.021}$ (0.859)
$\sigma_{8}$ |  $0.819^{+0.008}_{-0.012}$ (0.822) |  $0.842^{+0.018}_{-0.021}$ (0.862)
$S_{12}$ |  $0.811^{+0.010}_{-0.013}$ (0.811) |  $0.831^{+0.018}_{-0.021}$ (0.846)
$\Delta\chi^{2}_{\rm min}$ | 5.41 | 7.87
$\Delta{\rm AIC}$ | -0.59 | -0.13
Table 4: As in Tables 2 and 3, but for ULA and CDE_ULA.

ULA offers the interesting possibility of generating a peak in the EDE
fraction around $z_{\rm eq}$ while respecting the very tight constraints on
$\Omega_{\rm ede}$ at recombination Gómez-Valent  _et al._ (2021) (see Fig.
6). We show in Table 4 and the contours plots of Fig. 11 the constraints on
the main parameters of the model obtained with the data set Planck18+SNIa+BAO.
There is a large positive correlation between $f_{\rm ede}$ and the Hubble
parameter. Values of $H_{0}\sim 71$ km/s/Mpc accompanied by values of $f_{\rm
ede}\sim 10\%$, fall inside the $2\sigma$ region according to the marginalized
posterior. This seems to render the $H_{0}$ tension at the $\sim 2\sigma$
level. However, it has been recently demonstrated in Herold  _et al._ (2022);
Gómez-Valent (2022a) that volume effects play an important role in the
marginalization process in this model. The constraints obtained with the
profile distributions, which are not subject to these issues, are quite
different, and favor larger values of $f_{\rm ede}$ and $H_{0}$ (cf. Fig. 5
and Table II in Gómez-Valent (2022a)). The constraints derived with the
marginalized posteriors read $f_{\rm ede}<0.048$ and
$H_{0}=(68.4^{+1.3}_{-0.5})$ km/s/Mpc555Here we provide the location of the
peak of the posterior for $H_{0}$ and the corresponding uncertainties, instead
of the mean., whereas with the profile distribution we obtain $f_{\rm
ede}=0.052^{+0.022}_{-0.021}$ and $H_{0}=(69.8^{+0.9}_{-1.0})$ km/s/Mpc. A
peak in $f_{\rm ede}$ is found with the profile distribution approach,
$2\sigma$ away from 0 Herold  _et al._ (2022); Gómez-Valent (2022a).

Figure 12: Left plots: Comparison of the linear matter power spectrum of the
$\Lambda$CDM and ULA using the best-fit values obtained in the fitting
analysis with the Planck18+SNIa+BAO data set (cf. Sec. IV); Upper right plot:
This is the term with the window function that enters the equation of
$\sigma_{12}$ (33). Plotting this function is useful to see the range of $k$’s
that is more important in the computation of the rms of mass fluctuations at
$R_{12}=12$ Mpc; Lower right plot: Relative difference between the power
spectra of the two models, $\Delta_{P}\equiv(P_{\rm ULA}-P_{\Lambda{\rm
CDM}})/P_{\Lambda{\rm CDM}}$.

Regarding the LSS estimators, it is true that the model prefers values of
$S_{8}$ and $\sigma_{8}$ larger than in the $\Lambda$CDM, but due to the
points firstly raised in Sánchez (2020) and subsequently discussed also in
Gómez-Valent (2022a) and Sec. III of this paper, we think it is safer to
discuss the impact of the model on the linear perturbations in terms of the
parameter $\sigma_{12}$, although we also provide the fitting values of
$\sigma_{8}$, $S_{8}$ and $S_{12}$ in our tables. The constraint on
$\sigma_{12}$ obtained in ULA is almost identical to the one obtained in the
standard model, which means that ULA, when analyzed under the
Planck18+SNIa+BAO data set, does not prefer an enhancement of the rms of mass
fluctuations at scales of 12 Mpc. Actually one can see in the contours in the
$(H_{0},\sigma_{12})$-plane of Fig. 11 that it is possible to reach the
$2\sigma$ band of SH0ES keeping a low $\sigma_{12}\sim 0.79-0.8$. In Fig. 12
we show the shape of the matter power spectrum obtained with the best-fit
parameters in the $\Lambda$CDM and ULA. There is a significant relative
decrease of power at large scales in ULA ($k\lesssim 0.1$ Mpc-1), and an
increase at low scales. This two opposite behaviors compensate each other in
the computation of $\sigma_{12}$ (cf. formula (33)), yielding $\sigma_{12}\sim
0.8$. Although the model requires a larger value of $\omega_{\rm cdm}$ to
counteract the early iSW effect caused by $f_{\rm ede}$, it is also able to
accommodate a larger value of the cosmological constant, which is positively
correlated with the current matter energy density. Despite being larger than
in the $\Lambda$CDM, $\Omega^{(0)}_{m}$ still takes reasonable values, around
$\sim 0.3$, see Fig. 13 ($\Omega_{m}^{(0)}\sim 0.26$ in the standard model).
This larger matter fraction and the increase of $n_{s}$ cause the enhancement
of $P(k)$ at large $k$’s, which is compensated by the decrease at low $k$’s
that is induced by the bluer tilt of the primordial power spectrum.
Constraints on the power spectrum at low scales will be important to further
assess the viability of this model. According to the value of $\Delta{\rm
AIC}$ (see again Table 4) the model is not statistically preferred under the
Planck18+SNIa+BAO data set with respect to the $\Lambda$CDM, although it
allows an alleviation of the $H_{0}$ tension.

If we leave the sum of the neutrino masses $\sum m_{\nu}$ free in the Monte
Carlo analysis with ULA, we find very similar results to the case in which we
use only one massive neutrino of $0.06$ eV (cf. Table 5). We obtain $\sum
m_{\nu}<0.13$ eV at $95\%$ c.l., which is very close to the upper bound
derived with the $\Lambda$CDM, $\sum m_{\nu}<0.12$ eV. A larger mass of the
neutrino masses increases the uncertainties of $\sigma_{12}$, but only
slightly. This does not have a significant impact on the cosmological
tensions, as it has been also recently reported in Reeves  _et al._ (2022),
where the authors obtained $\sum m_{\nu}<0.15$ eV at $95\%$ c.l. using Planck
and BOSS galaxy clustering data, and showed that the constraints on $\sum
m_{\nu}$ are not affected by marginalization issues. We obtain a somewhat
tighter constraint, closer to the lower bound of the inverted mass hierarchy.

Figure 13: Scatter and contour plots at $1\sigma$ and $2\sigma$ c.l. in the
$(\omega_{\rm cdm},V_{0})$ and $(\omega_{\rm cdm},\Omega^{(0)}_{m})$ planes of
the ULA model, obtained from the fitting analysis with the Planck18+SNIa+BAO
dataset. We indicate the value of $\sigma_{12}$ at each point of the scatter
plot. There is a significant positive correlation between the current CDM
energy density and the value of the cosmological constant, which allows to
have $\Omega_{m}^{(0)}\sim 0.3$ in the $1\sigma$ region even for large values
of $\omega_{\rm cdm}\sim 0.13$. See the comments in the main text. Figure 14:
Confidence regions at $1\sigma$ and $2\sigma$ c.l. and one-dimensional
posterior distributions obtained with the Planck18+SNIa+BAO data set with the
$\Lambda$CDM and CDE_ULA.

If we extend the model to incorporate a coupling in the dark sector we do not
gain that much. The value of $\Delta{\rm AIC}$ is very close to zero. The
results for CDE_ULA are shown in Table 4 and Fig. 14. As expected, when
$\theta_{\rm ini}\lesssim\pi$ there is a symmetry between positive and
negative values of the coupling, which is broken when $\theta_{\rm ini}$ is
not sufficiently close to $\pi$. In all cases there is a preference for a
decay of the dark energy mass with the cosmic expansion (as found also with
CDE_EXP), and this is why for $\theta_{\rm ini}\ll\pi$ the fit prefers a
negative $\beta$. This explains why the peak at $\beta<0$ is larger. The
latter is, though, still compatible with 0 at only $2\sigma$ c.l.

| Planck18+SNIa+BAO [3 massive neutrinos]
---|---
Parameter |  $\Lambda$CDM | CDE_const | ULA
$10^{2}\omega_{b}$ |  $2.248^{+0.014}_{-0.013}$ (2.246) |  $2.238\pm 0.015$ (2.243) |  $2.267^{+0.017}_{-0.022}$ (2.271)
$\omega_{\rm cdm}$ |  $0.1187\pm 0.0008$ (0.1186) |  $0.1181^{+0.0010}_{-0.0008}$ (0.1177) |  $0.1222^{+0.0016}_{-0.0040}$ (0.1257)
$n_{s}$ |  $0.968\pm 0.004$ (0.970) |  $0.966\pm 0.004$ (0.966) |  $0.974^{+0.005}_{-0.010}$ (0.982)
$\tau_{\rm reio}$ |  $0.059^{+0.007}_{-0.008}$ (0.064) |  $0.058^{+0.007}_{-0.009}$ (0.056) |  $0.060^{+0.008}_{-0.009}$ (0.066)
$\sigma_{12}$ |  $0.798^{+0.007}_{-0.008}$ (0.803) |  $0.805^{+0.009}_{-0.011}$ (0.801) |  $0.797^{+0.009}_{-0.008}$ (0.801)
$H_{0}$ |  $67.89\pm 0.36$ (67.98) |  $68.34^{+0.45}_{-0.55}$ (68.64) |  $68.91^{+0.62}_{-1.20}$ (69.99)
$\sum m_{\nu}$ |  $0.06<x<0.12$ (0.07) |  $0.06<x<0.18$ (0.06) |  $0.06<x<0.13$ (0.07)
$f_{\rm ede}$ [%] | $-$ | $-$ |  $<5.05<9.47$ (7.59)
$z_{max}$ | $-$ | $-$ |  $8200^{+1700}_{-6000}$ (6709)
$\theta_{ini}$ | $-$ | $-$ | $>2.11(2.89)$
$\beta$ | $-$ |  $0.033\pm 0.015$ (0.031) | $-$
${\rm log_{10}(}m{\rm/eV)}$ | $-$ | $-$ |  $-26.64^{+0.39}_{-0.82}$ (-26.78)
${\rm log_{10}(}f{\rm/eV)}$ | $-$ | $-$ |  $26.49^{+0.27}_{-0.40}$ (26.43)
$r_{d}$ |  $147.32^{+0.22}_{-0.21}$ (147.37) |  $147.10^{+0.32}_{-0.23}$ (147.29) |  $145.31^{+2.10}_{-0.80}$ (143.50)
$M$ |  $-19.412^{+0.011}_{-0.010}$ (-19.410) |  $-19.400^{+0.012}_{-0.015}$ (-19.393) |  $-19.381^{+0.019}_{-0.037}$ (-19.346)
$S_{8}$ |  $0.816^{+0.009}_{-0.010}$ (0.821) |  $0.820\pm 0.011$ (0.813) |  $0.822^{+0.012}_{-0.013}$ (0.833)
$\sigma_{8}$ |  $0.808\pm 0.007$ (0.814) |  $0.819^{+0.010}_{-0.013}$ (0.817) |  $0.815^{+0.010}_{-0.013}$ (0.829)
$S_{12}$ |  $0.800\pm 0.008$ (0.806) |  $0.806^{+0.010}_{-0.011}$ (0.801) |  $0.808^{+0.011}_{-0.014}$ (0.820)
$\Delta\chi^{2}_{\rm min}$ | $-$ | 1.25 | 3.94
$\Delta{\rm AIC}$ | $-$ | -2.75 | -4.06
Table 5: As in the previous tables, but for the three models: $\Lambda$CDM,
CDE_const and ULA, with the Planck18+SNIa+BAO data set, and allowing the sum
of the neutrino masses $\sum m_{\nu}$ to vary in the Monte Carlo, as explained
in Sec. IV. The constraints on this parameter are given in eV at $95\%$ c.l.
See the discussion of these results in Sec. VI. Figure 15: Left plots:
Marginalized constraints on $M$ and $H_{0}$ (in km/s/Mpc) at 1$\sigma$ and
2$\sigma$ c.l. for all the models studied in this paper (cf. Sec. III),
obtained with the Planck18+SNIa+BAO data set. The green bands correspond to
the baseline values of $M$ and $H_{0}$ measured by SH0ES Riess  _et al._
(2022). ULA is the model that alleviates the Hubble tension the most. No
significant improvement is achieved by considering a coupling in the dark
sector (CDE_ULA) or more massive neutrinos (ULA+ $\sum m_{\nu}$); Upper right
plot: Same as the left plots, but for $\sigma_{12}$. The bands in this case
are taken from the marginalized Planck18+SNIa+BAO constraint on $\sigma_{12}$
in $\Lambda$CDM, see Table 2. The constraint on $\sigma_{12}$ from ULA is
almost identical to the one from the standard model; Lower right plot:
Marginalized constraint on $S_{8}$ at 1$\sigma$ and 2$\sigma$ c.l. for all the
models. The green band indicates the combined measurement of $S_{8}$ by
KiDS+VIKING-450 + DES-Y1, again under the assumption of the standard model
Joudaki  _et al._ (2020).

It is clear from the contours in the $(\beta,\sigma_{12})$-plane that when we
move away from the uncoupled scenario, i.e. from $\beta=0$, $\sigma_{12}$
grows. Large values of $\sigma_{12}\sim 0.85$ are not excluded by the
Planck18+SNIa+BAO data set. As already mentioned, the latter prefers values of
$\beta$ $\sim 2\sigma$ away from 0, and this is what generates a broader
posterior for $\sigma_{12}$, peaking also at larger values compared to the
uncoupled ULA model. The required values of $f_{\rm ede}$ to alleviate the
Hubble tension can only be obtained in the vicinity of $\beta=0$, so the
preference of the data for a non-null coupling, despite being still mild,
inevitably leads to a tighter constraint on the maximum EDE fraction, which is
now below $7.7\%$ at $95\%$ c.l. Unfortunately, in the absence of a prior from
SH0ES, we do not find an alleviation of the tensions caused by the coupling.
We have verified that these results are not strongly affected by volume
effects.

Before closing the discussion on our results, we would like to remark that the
quantification of the Hubble tension, when formulated as a tension between the
value of the absolute value of SNIa $M$ obtained from the distance ladder by
SH0ES and the value obtained from the fitting analysis involving high-$z$
supernovae, give similar but not exactly equal results Gómez-Valent  _et al._
(2021). The latter tends to lead to slightly larger estimates of the tension.
In the left plots of Fig. 15 we summarize our constraints on $H_{0}$ and $M$
for all the models studied in this paper, and compare them to the measurement
carried out by the SH0ES Team Riess  _et al._ (2022). We find that ULA is the
model that alleviates more the $H_{0}$ tension, and that neither a constant
coupling between the dark components nor a larger sum of the neutrino masses
can help to further mitigate it.

In the right plots of Fig. 15 we summarize the constraints on $\sigma_{12}$
and $S_{8}$ obtained for the various models. Some comments are in order. All
the values of these LSS estimators are compatible with the result in the
$\Lambda$CDM. There is no model leading to a significant decrease of these
quantities. The constraints on $\sigma_{12}$ in ULA and ULA+$\sum m_{\nu}$ are
actually almost identical to those found in the standard model. The estimation
of the level of tension between the models and the weak lensing measurement
from KiDS+VIKING-450+DES-Y1 is not straightforward and can be misleading for
the issues discussed in Sánchez (2020), and also in Sec. III. Moreover, it is
important to bear in mind that the observational constraint on $S_{8}$ from
KiDS+VIKING-450+DES-Y1 (see the green bands shown in the lower right plot of
Fig. 15) has been obtained under the assumption of the $\Lambda$CDM, so it
could lead to an overestimation of the level of tension in the context of non-
standard cosmologies.

##  VII Conclusions

In this paper we have studied some coupled and uncoupled early dark energy
(EDE) models with different shapes of the EDE fraction, confronting them with
the CMB data from Planck, the Pantheon compilation of supernovae of Type Ia
and data on baryon acoustic oscillations, without including any prior on
$H_{0}$ or $M$ from SH0ES. We have considered three different forms of the
potential energy density for the scalar field, to wit: a constant potential,
an exponential potential (20) that produces a scaling regime in the matter-
and radiation-dominated epochs, and an ultra-light axion-like (ULA) potential
(32) able to generate a peak in $\Omega_{\rm ede}$, typically around the
matter-radiation equality time. We have discussed in detail the phenomenology
of these models. Motivated by the fact that EDE has some difficulties in
relieving simultaneously the $H_{0}$ and LSS tensions, we have explored in
this paper the impact of a coupling between EDE and dark matter, and the
presence of massive neutrinos. The coupling can make the dark matter mass to
decay with the expansion, leading to a faster decrease of the dark matter
energy density. This fact could in principle slow down the clustering of
matter in the universe, and massive neutrinos could also help to suppress the
amount of large-scale structure at low scales. Nevertheless, the coupling can
also enhance the matter power spectrum due to the fifth force. In order to
elucidate their net effect, we have put these ideas to the test, by performing
dedicated fitting analyses.

When the EDE potential is constant we find under the Planck18+SNIa+BAO data
set a preference for a non-null coupling at $2\sigma$ c.l. This result
resonates well with previous studies in the literature Pettorino (2013); Ade
_et al._ (2016); Gómez-Valent  _et al._ (2020), and is not caused by volume
effects Gómez-Valent (2022a). Nevertheless, if we penalize the additional
complexity of the model through the calculation of the AIC we find that the
improvement in the description of the data is not sufficient to justify the
need of the interaction in the dark sector. The $H_{0}$ tension is alleviated,
yes, but persists at $\sim 3\sigma$. Massive neutrinos give some room for
lower values of $\sigma_{12}$, but the gain is not substantial, and is
accompanied by a slight decrease of $H_{0}$ as well. The upper bound of the
sum of the neutrino masses reads $\sum m_{\nu}<0.18$ eV at $2\sigma$ c.l., and
is considerably larger than in the $\Lambda$CDM, $\sum m_{\nu}<0.12$ eV. This
might indicate that probing the neutrino mass hierarchy with cosmological data
in a model-independent way might be quite difficult. Constraints obtained in
the context of other models can be even more loosen, see e.g. Solà Peracaula
_et al._ (2020). This complication will presumably remain with the advent of
future data.

Moreover our results show that EDE models with an exponential potential have
almost no impact on the cosmological tensions. We already reached this
conclusion in Gómez-Valent  _et al._ (2021), but here we have also checked the
effects introduced by the coupling. The latter allows us to describe
considerably better the CMB data and we even find a moderate $\sim 2\sigma$
hint for a non-null interaction in the dark sector when only the Planck
likelihood is considered. Nevertheless, this signal is diluted when we combine
CMB with the SNIa and BAO data sets, and the tensions remain large.

ULA is able to decrease the $H_{0}$ tension below the $2\sigma$ level,
basically because it allows to have a larger EDE fraction around the matter-
radiation equation time while respecting the stringent bounds at recombination
Gómez-Valent  _et al._ (2021). The tension is even lower when it is quantified
with the profile distribution method Gómez-Valent (2022a). The values of the
LSS estimator $\sigma_{12}$ are similar to those found in the standard model.
We have seen, though, that the shape of the power spectrum exhibits some
important differences at small and large scales, so it could be useful to
prove in the future the non-linear regime to study the viability of this
model. For instance, there could be significant differences between the number
counts of clusters in ULA and the $\Lambda$CDM obtained using the best-fit
values of the parameters from our analysis. By allowing the mass of the
neutrinos to take larger values in the fitting analysis we get an upper bound
which is extremely close to the one obtained in the $\Lambda$CDM, $\sum
m_{\nu}<0.13$ eV at $2\sigma$ c.l., and its posterior peaks very close to the
minimum value allowed by the experiments on atmospheric and solar neutrinos,
at (0.07 eV). Massive neutrinos do not induce any significant shift on the
other parameters and, hence, they are unable to lower the tensions. This has
been reported recently also in Reeves  _et al._ (2022), and this conclusion is
not affected by marginalization issues of any kind. Unfortunately, when we
activate the coupling in the context of ULA we do not appreciate any
improvement regarding the tensions either.

The cosmological tensions tackled in this paper cannot be fully alleviated by
an interaction with a constant coupling between the two dark components nor
neutrinos with larger masses. In the best case, namely ULA, the $H_{0}$
tension is reduced to $\sim 2\sigma$, and the amount of clustering at linear
scales is similar to the $\Lambda$CDM. A more detailed treatment of the non-
linear scales in these models and the inclusion of weak lensing and galaxy
clustering data in the analysis could shed more light on the status of the
tensions, of course. This could be done by performing scale cuts on CMB and
weak lensing data to explore only the effect of linear scales or by modeling
accurately the non-linear power spectrum (on the lines of e.g. Casas  _et al._
(2016)) to exploit the statistical information contained in the full range of
scales covered by the CMB and LSS data. This goes beyond the scope of this
paper, and is left for a future work.

Acknowledgements

AGV is funded by the Istituto Nazionale di Fisica Nucleare (INFN) through the
project of the InDark INFN Special Initiative: “Dark Energy and Modified
Gravity Models in the light of Low-Redshift Observations” (n. 22425/2020). ZZ
is supported by the DFG Research Training Group “Particle Physics beyond the
Standard Model” (GRK 1940). LA acknowledges support from DFG project
456622116. The authors would like to thank Dr. Elmar Bittner for his precious
technical help in the use of the cluster of the ITP Heidelberg, which has been
crucial for the completion of this project.

Figure 16: Comoving wave mode at horizon crossing, $k_{\rm hor}$, and comoving
dark energy mass, $am_{\phi}$, as a function of the redshift for several
CDE_EXP (upper plot) and CDE_ULA (lower plot) models. The parameters that are
not specified explicitly in the legends are set as in Fig. 2 and the left
plots of Fig. 6, respectively.

## Appendix A: Impact of the dark energy mass on the LSS

In this brief appendix we discuss the impact that the non-zero dark energy
mass $m^{2}_{\phi}=d^{2}V/d\phi^{2}$ has on the LSS both in CDE_EXP and
CDE_ULA. In the upper plot of Fig. 16 we can see that for typical values of
$\Omega_{\rm ede}^{\rm RD}$ and the coupling in the CDE_EXP models, the
comoving modes that are relevant for the LSS, i.e. $10^{-2}\,{\rm
Mpc}^{-1}\lesssim k\lesssim 1\,{\rm Mpc}^{-1}$, are already larger than the
comoving mass of the DE at the horizon crossing, i.e. when the modes of
interest are equal to $k_{\rm hor}$ in the plot. The corresponding comoving
scales rapidly become much smaller than the range of the fifth force. This
means that the mass terms $m_{\phi}$ in the equation of the density contrast
of DM (16) can be safely neglected in this model.

| Planck18
---|---
$\Delta\chi^{2}_{i}$ |  $\Lambda$CDM | CDE_const | EXP | CDE_EXP
CMB highl | 2351.83 | 2.00 | 3.19 | 7.54
CMB EE lowl | 396.03 | 0.32 | -0.13 | -0.59
CMB TT lowl | 23.27 | -0.08 | -0.21 | 0.72
CMB lens | 8.81 | -0.31 | -0.10 | -0.16
Total | 2779.94 | 1.93 | 2.75 | 7.51
Table 6: Differences
$\Delta\chi^{2}_{i,j}\equiv\chi^{2}_{i,\Lambda}-\chi^{2}_{i,j}$ for the
individual data sets $i$ and non-standard cosmological models $j$, using the
best-fit values obtained in the analyses with Planck18 data. In the last row
we show the total differences, $\Delta\chi^{2}_{{\rm
min},j}=\chi^{2}_{\Lambda,{\rm min}}-\chi^{2}_{j,{\rm min}}$. For the
$\Lambda$CDM, instead, we show in the first column the individual values of
$\chi^{2}_{i,\Lambda}$, and their sum.

In the CDE_ULA model the evolution of the DE mass is very different from the
one in CDE_EXP. The aforementioned modes do not feel the non-null DE mass when
these scales reenter the horizon, since for these modes $k\gg am_{\phi}$
already at the crossing time. This is shown in the lower plot of Fig. 16. The
scales that reenter the horizon at $30\lesssim z\lesssim 10^{3}$, though, are
sensitive to the oscillatory behavior of the DE mass. Nevertheless, these
modes are not relevant for the LSS, and the mean value of $am_{\phi}$ is still
much smaller than $k_{\rm hor}$.

## Appendix B: Shooting method for ULA and CDE_ULA

We implement a shooting method in ULA and CDE_ULA in order to map every
combination of the input parameters $(f_{\rm ede},z_{\rm max})$ to the
corresponding parameters $(m,f)$ of the ULA potential, see formula (32) in
Sec. III. We apply the following steps:

| Planck18+SNIa+BAO
---|---
$\Delta\chi^{2}_{i}$ |  $\Lambda$CDM | CDE_const | EXP | CDE_EXP | ULA | CDE_ULA
CMB highl | 2349.83 | 1.19 | -0.4 | 2.28 | 2.58 | 4.33
CMB EE lowl | 397.90 | 1.87 | 1.12 | 2.19 | 0.84 | 2.16
CMB TT lowl | 23.51 | -0.79 | 0.73 | 1.06 | 1.88 | 2.08
CMB lens | 9.14 | 0.20 | 0.22 | 0.05 | -0.39 | -1.29
SNIa | 1025.77 | -0.26 | -1.19 | -0.16 | 0.03 | 0.14
BAO | 9.45 | 0.19 | -0.65 | -1.71 | 0.43 | 0.42
Total | 3815.61 | 2.39 | -0.15 | 3.75 | 5.41 | 7.87
Table 7: As in TABLE 6, but for the analyses with the Planck18+SNIa+BAO data
set. Planck18+SNIa+BAO [3 massive neutrinos]
---
$\Delta\chi^{2}_{i}$ |  $\Lambda$CDM | CDE_const | ULA
CMB highl | 2348.37 | -1.72 | 3.88
CMB EE lowl | 399.14 | 2.71 | -0.24
CMB TT lowl | 22.76 | -0.68 | 1.44
CMB lens | 8.77 | -0.74 | -0.41
SNIa | 1025.75 | 0.12 | -0.13
BAO | 9.64 | 0.50 | 0.42
$\Delta m^{2}_{21},\Delta m^{2}_{31}$ | 0.23 | 0.10 | -1.18
Total | 3814.66 | 1.25 | 3.94
Table 8: As in TABLE 7, but allowing the neutrino masses to vary in the Monte
Carlo runs. $\Delta m^{2}_{21},\Delta m^{2}_{31}$ refer to the priors employed
for these quantities, see Sec. IV.

  * •

We start with an initial guess for $\Big{(}{\rm log_{10}(}m{\rm/eV)},{\rm
log_{10}(}f{\rm/eV)}\Big{)}$. We build a grid around it and find the location
of $x^{(1)}\equiv\Big{(}{\rm log_{10}(}m{\rm/eV)^{(1)}},{\rm
log_{10}(}f{\rm/eV)^{(1)}}\Big{)}$ that leads to the minimum $\chi^{2}$ value
on the grid. The associated values of the cosmological (input) parameters are
labeled as $y^{(1)}\equiv\Big{(}f_{\rm ede}^{(1)},\rm
log_{10}(z_{max})^{(1)}\Big{)}$. This step is important to start the main part
of the shooting sufficiently close to our final target $y^{(F)}$.

  * •

We build the Jacobian matrix $J_{ij}^{(1)}\equiv\frac{\partial
y_{i}^{(1)}}{\partial x_{j}^{(1)}}$, approximating it as follows,

$J^{(1)}\approx\begin{pmatrix}\frac{y_{1}(x_{1}+\delta
x_{1},x_{2})-y_{1}^{(1)}}{\delta x_{1}}&\frac{y_{1}(x_{1},x_{2}+\delta
x_{2})-y_{1}^{(1)}}{\delta x_{2}}\\\ \frac{y_{2}(x_{1}+\delta
x_{1},x_{2})-y_{2}^{(1)}}{\delta x_{1}}&\frac{y_{2}(x_{1},x_{2}+\delta
x_{2})-y_{2}^{(1)}}{\delta x_{2}}\\\ \end{pmatrix}\,,$

with $\delta x_{1}$ and $\delta x_{2}$ two sufficiently small increments.

  * •

A new trial $x^{(2)}$ is related to $y^{(\rm F)}$ and the Jacobian matrix
$J^{(1)}$ through the relation:

$x^{(2)}=x^{(1)}+\Big{(}J^{(1)}\Big{)}^{-1}(y^{\rm(F)}-y^{(1)}),$ (39)

where the matrix $\Big{(}J^{(1)}\Big{)}^{-1}$ is the inverse of $J^{(1)}$. Eq.
(39) is nothing else than the Taylor expansion of $x$ around $y^{(1)}$
evaluated at $y^{(F)}$, truncated at linear order. We compute the new Jacobian
matrix $J^{(2)}$ following the previous steps to obtain $x^{(3)}$, and we
repeat this iterative process until reaching convergence (within the desired
precision).

## Appendix C: Tables with the individual $\Delta\chi^{2}_{i}$

This appendix contains Tables 6-8, where we list the contribution from the
individual data sets to the total $\chi^{2}_{\rm min}$ for each model and
fitting analysis carried out in this work. These tables complement the
information reported in the tables of Sec. VI.

## References

  * Peebles and Ratra (2003) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), arXiv:astro-ph/0207347 .
  * Hinshaw  _et al._ (2013) G. Hinshaw _et al._ (WMAP), Astrophys. J. Suppl. 208, 19 (2013), arXiv:1212.5226 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Eisenstein  _et al._ (2005) D. J. Eisenstein _et al._ (SDSS), Astrophys. J. 633, 560 (2005), arXiv:astro-ph/0501171 [astro-ph] .
  * Cole  _et al._ (2005) S. Cole _et al._ (2dFGRS), Mon. Not. Roy. Astron. Soc. 362, 505 (2005), arXiv:astro-ph/0501174 [astro-ph] .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 [astro-ph] .
  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 [astro-ph] .
  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * Solà Peracaula (2022) J. Solà Peracaula, Phil. Trans. Roy. Soc. Lond. A 380, 20210182 (2022), arXiv:2203.13757 [gr-qc] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, Nature Astron. 3, 891 (2019), arXiv:1907.10625 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021a) E. Di Valentino _et al._ , Astropart. Phys. 131, 102604 (2021a), arXiv:2008.11285 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Aluri  _et al._ (2022) P. K. Aluri _et al._ , (2022), arXiv:2207.05765 [astro-ph.CO] .
  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Joudaki  _et al._ (2020) S. Joudaki _et al._ , Astron. Astrophys. 638, L1 (2020), arXiv:1906.09262 [astro-ph.CO] .
  * Nunes and Vagnozzi (2021) R. C. Nunes and S. Vagnozzi, Mon. Not. Roy. Astron. Soc. 505, 5427 (2021), arXiv:2106.01208 [astro-ph.CO] .
  * Macaulay  _et al._ (2013) E. Macaulay, I. K. Wehus, and H. K. Eriksen, Phys. Rev. Lett. 111, 161301 (2013), arXiv:1303.6583 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2017) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, Mon. Not. Roy. Astron. Soc. 465, 1757 (2017), arXiv:1606.00439 [astro-ph.CO] .
  * Gómez-Valent and Solà (2017) A. Gómez-Valent and J. Solà, EPL 120, 39001 (2017), arXiv:1711.00692 [astro-ph.CO] .
  * Nesseris  _et al._ (2017) S. Nesseris, G. Pantazis, and L. Perivolaropoulos, Phys. Rev. D 96, 023542 (2017), arXiv:1703.10538 [astro-ph.CO] .
  * Gómez-Valent and Solà Peracaula (2018) A. Gómez-Valent and J. Solà Peracaula, Mon. Not. Roy. Astron. Soc. 478, 126 (2018), arXiv:1801.08501 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021b) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021b), arXiv:2103.01183 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2019a) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, Astrophys. J. Lett. 886, L6 (2019a), arXiv:1909.02554 [astro-ph.CO] .
  * Ballesteros  _et al._ (2020) G. Ballesteros, A. Notari, and F. Rompineve, JCAP 11, 024 (2020), arXiv:2004.05049 [astro-ph.CO] .
  * Braglia  _et al._ (2020) M. Braglia, M. Ballardini, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, K. Koyama, and D. Paoletti, Phys. Rev. D 102, 023529 (2020), arXiv:2004.11161 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2020) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, Class. Quant. Grav. 37, 245003 (2020), arXiv:2006.04273 [astro-ph.CO] .
  * Braglia  _et al._ (2021) M. Braglia, M. Ballardini, F. Finelli, and K. Koyama, Phys. Rev. D 103, 043528 (2021), arXiv:2011.12934 [astro-ph.CO] .
  * Benevento  _et al._ (2022) G. Benevento, J. A. Kable, G. E. Addison, and C. L. Bennett, Astrophys. J. 935, 156 (2022), arXiv:2202.09356 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Niedermann and Sloth (2021) F. Niedermann and M. S. Sloth, Phys. Rev. D 103, L041303 (2021), arXiv:1910.10739 [astro-ph.CO] .
  * Ye and Piao (2020) G. Ye and Y.-S. Piao, Phys. Rev. D 101, 083507 (2020), arXiv:2001.02451 [astro-ph.CO] .
  * Gómez-Valent  _et al._ (2021) A. Gómez-Valent, Z. Zheng, L. Amendola, V. Pettorino, and C. Wetterich, Phys. Rev. D 104, 083536 (2021), arXiv:2107.11065 [astro-ph.CO] .
  * Wang and Piao (2022) H. Wang and Y.-S. Piao, Phys. Lett. B 832, 137244 (2022), arXiv:2201.07079 [astro-ph.CO] .
  * Jedamzik and Pogosian (2020) K. Jedamzik and L. Pogosian, Phys. Rev. Lett. 125, 181302 (2020), arXiv:2004.09487 [astro-ph.CO] .
  * Liu  _et al._ (2020) M. Liu, Z. Huang, X. Luo, H. Miao, N. K. Singh, and L. Huang, Sci. China Phys. Mech. Astron. 63, 290405 (2020), arXiv:1912.00190 [astro-ph.CO] .
  * Sekiguchi and Takahashi (2021) T. Sekiguchi and T. Takahashi, Phys. Rev. D 103, 083507 (2021), arXiv:2007.03381 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2021) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, EPL 134, 19001 (2021), arXiv:2102.12758 [astro-ph.CO] .
  * Wetterich (1995) C. Wetterich, Astron. Astrophys. 301, 321 (1995), arXiv:hep-th/9408025 .
  * Amendola (2000a) L. Amendola, Phys. Rev. D 62, 043511 (2000a), arXiv:astro-ph/9908023 .
  * Poulin  _et al._ (2021) V. Poulin, T. L. Smith, and A. Bartlett, Phys. Rev. D 104, 123550 (2021), arXiv:2109.06229 [astro-ph.CO] .
  * Hill  _et al._ (2022) J. C. Hill _et al._ , Phys. Rev. D 105, 123536 (2022), arXiv:2109.04451 [astro-ph.CO] .
  * Smith  _et al._ (2022) T. L. Smith, M. Lucca, V. Poulin, G. F. Abellan, L. Balkenhol, K. Benabed, S. Galli, and R. Murgia, Phys. Rev. D 106, 043526 (2022), arXiv:2202.09379 [astro-ph.CO] .
  * Hill  _et al._ (2020) J. C. Hill, E. McDonough, M. W. Toomey, and S. Alexander, Phys. Rev. D 102, 043507 (2020), arXiv:2003.07355 [astro-ph.CO] .
  * Smith  _et al._ (2021) T. L. Smith, V. Poulin, J. L. Bernal, K. K. Boddy, M. Kamionkowski, and R. Murgia, Phys. Rev. D 103, 123542 (2021), arXiv:2009.10740 [astro-ph.CO] .
  * Murgia  _et al._ (2021) R. Murgia, G. F. Abellán, and V. Poulin, Phys. Rev. D 103, 063502 (2021), arXiv:2009.10733 [astro-ph.CO] .
  * Sánchez (2020) A. G. Sánchez, Phys. Rev. D 102, 123511 (2020), arXiv:2002.07829 [astro-ph.CO] .
  * Gómez-Valent (2022a) A. Gómez-Valent, Phys. Rev. D 106, 063506 (2022a), arXiv:2203.16285 [astro-ph.CO] .
  * Vagnozzi (2021) S. Vagnozzi, Phys. Rev. D 104, 063524 (2021), arXiv:2105.10425 [astro-ph.CO] .
  * Karwal  _et al._ (2022) T. Karwal, M. Raveri, B. Jain, J. Khoury, and M. Trodden, Phys. Rev. D 105, 063535 (2022), arXiv:2106.13290 [astro-ph.CO] .
  * McDonough  _et al._ (2022) E. McDonough, M.-X. Lin, J. C. Hill, W. Hu, and S. Zhou, Phys. Rev. D 106, 043525 (2022), arXiv:2112.09128 [astro-ph.CO] .
  * Aubourg  _et al._ (2015) E. Aubourg _et al._ , Phys. Rev. D 92, 123516 (2015), arXiv:1411.1074 [astro-ph.CO] .
  * Cuesta  _et al._ (2015) A. J. Cuesta, L. Verde, A. Riess, and R. Jimenez, Mon. Not. Roy. Astron. Soc. 448, 3463 (2015), arXiv:1411.1094 [astro-ph.CO] .
  * Feeney  _et al._ (2019) S. M. Feeney, H. V. Peiris, A. R. Williamson, S. M. Nissanke, D. J. Mortlock, J. Alsing, and D. Scolnic, Phys. Rev. Lett. 122, 061105 (2019), arXiv:1802.03404 [astro-ph.CO] .
  * Camarena and Marra (2020a) D. Camarena and V. Marra, Mon. Not. Roy. Astron. Soc. 495, 2630 (2020a), arXiv:1910.14125 [astro-ph.CO] .
  * Lesgourgues and Pastor (2006) J. Lesgourgues and S. Pastor, Phys. Rept. 429, 307 (2006), arXiv:astro-ph/0603494 .
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), arXiv:astro-ph/9506072 .
  * Amendola (2000b) L. Amendola, Mon. Not. Roy. Astron. Soc. 312, 521 (2000b), arXiv:astro-ph/9906073 .
  * Solà Peracaula  _et al._ (2019b) J. Solà Peracaula, A. Gómez-Valent, and J. de Cruz Pérez, Phys. Dark Univ. 25, 100311 (2019b), arXiv:1811.03505 [astro-ph.CO] .
  * Barros  _et al._ (2019) B. J. Barros, L. Amendola, T. Barreiro, and N. J. Nunes, JCAP 01, 007 (2019), arXiv:1802.09216 [astro-ph.CO] .
  * Gómez-Valent  _et al._ (2020) A. Gómez-Valent, V. Pettorino, and L. Amendola, Phys. Rev. D 101, 123513 (2020), arXiv:2004.00610 [astro-ph.CO] .
  * Amendola and Tsujikawa (2015) L. Amendola and S. Tsujikawa,  _Dark Energy: Theory and Observations_ (Cambridge Univ. Press, Cambridge, 2015).
  * Amendola (2001) L. Amendola, Phys. Rev. Lett. 86, 196 (2001), arXiv:astro-ph/0006300 .
  * Pettorino  _et al._ (2012) V. Pettorino, L. Amendola, C. Baccigalupi, and C. Quercellini, Phys. Rev. D 86, 103507 (2012), arXiv:1207.3293 [astro-ph.CO] .
  * Pettorino (2013) V. Pettorino, Phys. Rev. D 88, 063519 (2013), arXiv:1305.7457 [astro-ph.CO] .
  * Xia (2013) J.-Q. Xia, JCAP 1311, 022 (2013), arXiv:1311.2131 [astro-ph.CO] .
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A14 (2016), arXiv:1502.01590 [astro-ph.CO] .
  * van de Bruck  _et al._ (2017) C. van de Bruck, J. Mifsud, and J. Morrice, Phys. Rev. D95, 043513 (2017), arXiv:1609.09855 [astro-ph.CO] .
  * van de Bruck and Mifsud (2018) C. van de Bruck and J. Mifsud, Phys. Rev. D97, 023506 (2018), arXiv:1709.04882 [astro-ph.CO] .
  * Agrawal  _et al._ (2021) P. Agrawal, G. Obied, and C. Vafa, Phys. Rev. D 103, 043523 (2021), arXiv:1906.08261 [astro-ph.CO] .
  * Goh  _et al._ (2022) L. Goh, A. Gómez-Valent, and V. Pettorino, (2022), in preparation.
  * Wetterich (1988) C. Wetterich, Nucl. Phys. B 302, 668 (1988), arXiv:1711.03844 [hep-th] .
  * Copeland  _et al._ (1998) E. J. Copeland, A. R. Liddle, and D. Wands, Phys. Rev. D 57, 4686 (1998), arXiv:gr-qc/9711068 .
  * Barreiro  _et al._ (2000) T. Barreiro, E. J. Copeland, and N. J. Nunes, Phys. Rev. D 61, 127301 (2000), arXiv:astro-ph/9910214 .
  * Doran  _et al._ (2001) M. Doran, M. J. Lilley, J. Schwindt, and C. Wetterich, Astrophys. J. 559, 501 (2001), arXiv:astro-ph/0012139 .
  * Doran and Lilley (2002) M. Doran and M. Lilley, Mon. Not. Roy. Astron. Soc. 330, 965 (2002), arXiv:astro-ph/0104486 .
  * Poulin  _et al._ (2018) V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Phys. Rev. D 98, 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .
  * Pettorino  _et al._ (2013) V. Pettorino, L. Amendola, and C. Wetterich, Phys. Rev. D 87, 083009 (2013), arXiv:1301.5279 [astro-ph.CO] .
  * D’Amico  _et al._ (2021) G. D’Amico, L. Senatore, P. Zhang, and H. Zheng, JCAP 05, 072 (2021), arXiv:2006.12420 [astro-ph.CO] .
  * Aiola  _et al._ (2020) S. Aiola _et al._ (ACT), JCAP 12, 047 (2020), arXiv:2007.07288 [astro-ph.CO] .
  * Semenaite  _et al._ (2022) A. Semenaite _et al._ , Mon. Not. Roy. Astron. Soc. 512, 5657 (2022), arXiv:2111.03156 [astro-ph.CO] .
  * Herold  _et al._ (2022) L. Herold, E. G. M. Ferreira, and E. Komatsu, Astrophys. J. Lett. 929, L16 (2022), arXiv:2112.12140 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. M. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Gómez-Valent (2022b) A. Gómez-Valent, Phys. Rev. D 105, 043528 (2022b), arXiv:2111.15450 [astro-ph.CO] .
  * Benisty  _et al._ (2022) D. Benisty, J. Mifsud, J. L. Said, and D. Staicova, (2022), arXiv:2202.04677 [astro-ph.CO] .
  * Camarena and Marra (2020b) D. Camarena and V. Marra, Phys. Rev. Res. 2, 013028 (2020b), arXiv:1906.11814 [astro-ph.CO] .
  * Efstathiou (2021) G. Efstathiou, Mon. Not. Roy. Astron. Soc. 505, 3866 (2021), arXiv:2103.08723 [astro-ph.CO] .
  * Carter  _et al._ (2018) P. Carter, F. Beutler, W. J. Percival, C. Blake, J. Koda, and A. J. Ross, Mon. Not. Roy. Astron. Soc. 481, 2371 (2018), arXiv:1803.01746 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Kazin  _et al._ (2014) E. A. Kazin _et al._ , Mon. Not. Roy. Astron. Soc. 441, 3524 (2014), arXiv:1401.0358 [astro-ph.CO] .
  * Abbott  _et al._ (2019) T. M. C. Abbott _et al._ (DES), Mon. Not. Roy. Astron. Soc. 483, 4866 (2019), arXiv:1712.06209 [astro-ph.CO] .
  * Neveux  _et al._ (2020) R. Neveux _et al._ , Mon. Not. Roy. Astron. Soc. 499, 210 (2020), arXiv:2007.08999 [astro-ph.CO] .
  * du Mas des Bourboux  _et al._ (2020) H. du Mas des Bourboux _et al._ , Astrophys. J. 901, 153 (2020), arXiv:2007.08995 [astro-ph.CO] .
  * Uzan (2011) J.-P. Uzan, Living Rev. Rel. 14, 2 (2011), arXiv:1009.5514 [astro-ph.CO] .
  * Seto and Toda (2021) O. Seto and Y. Toda, Phys. Rev. D 103, 123501 (2021), arXiv:2101.03740 [astro-ph.CO] .
  * Loureiro  _et al._ (2019) A. Loureiro _et al._ , Phys. Rev. Lett. 123, 081301 (2019), arXiv:1811.02578 [astro-ph.CO] .
  * Esteban  _et al._ (2017) I. Esteban, M. C. González-García, M. Maltoni, I. Martínez-Soler, and T. Schwetz, JHEP 01, 087 (2017), arXiv:1611.01514 [hep-ph] .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 07, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 02, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Metropolis  _et al._ (1953) A. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Hastings (1970) W. Hastings, Biometrika 57, 97 (1970).
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Schöneberg  _et al._ (2022) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Trotta (2017) R. Trotta (2017) arXiv:1701.01467 [astro-ph.CO] .
  * Ade  _et al._ (2014) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 566, A54 (2014), arXiv:1311.1657 [astro-ph.CO] .
  * Campeti and Komatsu (2022) P. Campeti and E. Komatsu, (2022), arXiv:2205.05617 [astro-ph.CO] .
  * Akaike (1974) H. Akaike, IEEE Trans. Autom. Control 19, 716 (1974).
  * Reeves  _et al._ (2022) A. Reeves, L. Herold, S. Vagnozzi, B. D. Sherwin, and E. G. M. Ferreira, (2022), arXiv:2207.01501 [astro-ph.CO] .
  * Casas  _et al._ (2016) S. Casas, L. Amendola, M. Baldi, V. Pettorino, and A. Vollmer, JCAP 01, 045 (2016), arXiv:1508.07208 [astro-ph.CO] .