The vanilla concordance model of cosmology, or standard $\Lambda$CDM model (the current standard model of cosmology with flat three-dimensional geometry), is based on the Friedmann-Lemaître-Robertson-Walker (FLRW) metric and has been a rather successful paradigm for the phenomenological description of the universe for more than three decades [1,2]. Its consolidation after a solid observational underpinning, however, was only possible in the late nineties [3]. The vanilla model has remained robust and unbeaten for a long time, as it is essentially consistent with a large body of observations. These have indeed provided strong support for a spatially flat and accelerating universe in the present time. The ultimate cause of such an acceleration is currently unknown, but it is attributed to an energy component in the universe popularly called “dark energy" (DE), which may adopt a large number of picturesque forms depending on the favorite theoretical preference of different cosmologists  see e.g. [4] for a large variety of options. The DE constitutes $\sim 70\%$ of the total energy density of the universe and presumably possesses enough negative pressure as to produce the observed cosmic acceleration. Nevertheless, the nature of the DE remains still a complete mystery. The simplest candidate is the cosmological term in Einstein’s equations, $\Lambda$, usually assumed to be constant, which is why it is usually called the cosmological constant (CC) [5,6]. Consistent observational measurements of $\Lambda$ (treating it as a mere fit parameter) made independently in the last quarter of a century using distant type Ia supernovae (SNIa), the baryonic acoustic oscillations (BAO) and the anisotropies of the cosmic microwave background (CMB), have put the very foundations of the concordance $\Lambda$CDM model of cosmology [7,8,9,10,11,12,13,14,15].

Despite the vanilla $(\Lambda$CDM) model fares relatively well with the current observational data, it traditionally suffers from a variety of problems of different kinds which seriously challenge its credibility. For a long time people somehow decided to turn a blind eye on the deepest questions and also on different spots and wrinkles which perturb that flawless and immaculate condition. The profound theoretical problems (and the practical wrinkles as well) are nonetheless still there alive and kicking, whether we wish to look at them or not. First and foremost, the hypothetical existence of dark matter (DM) still lacks of direct observational evidence. On a deeper level of mystery, the nature and origin of the DE (the dominant component of the cosmic energy budget) still lies in the limbo of the most unfathomable cosmological riddles. Basically because if we admit the simplest proposal for the DE, that is to say, the cosmological constant $\Lambda$, one has to cope with the ‘cosmological constant problem’[16], perhaps the most inscrutable of all problems in theoretical physics and cosmology ever [17]. It manifests itself in a dual manner, to wit: the fine-tuning problem associated with the large value of $\Lambda$ predicted by most theoretical approaches (“the old CC problem” [16]); and also what has become customary to call the ‘cosmic coincidence problem’ [18], see also [17,19,20] for a discussion of these enigmas, which lie in the interface between cosmology and quantum field theory (QFT).

The toughest conundrum of all is probably that of explaining the relation between $\Lambda$ and the vacuum energy density (VED): ${\rho }_{\mathrm{vac}}=\Lambda /\left(8\pi G_N\right)$, where $G_N$ is Newton’s constant. The interrelationship between VED and $\Lambda$ in the general quantum theory context has been assessed by theoretical physicists since more than a century ago, as of the days of W. Nernst and W. Pauli. At that time the issue was already troublesome [17]. But the most severe implications in the cosmological arena took shape only with the development of the formal aspects of QFT. It is in this modern theoretical context where the notion of VED seems to cause a serious conflict with the cosmological measurements, the reason being that the typical contribution from the vacuum fluctuations of any quantum field of mass m is expected (on mere dimensional grounds) to be proportional to the quartic power of its mass: ${\rho }_{\mathrm{vac}}\propto m^4$ [21]. Such a prediction must be compared with the measured value of $\Lambda$ expressed in terms of the corresponding VED, which is ${\rho }_{\mathrm{vac}}\sim {10}^{-47}$ GeV${}^4$ in natural units. This is extremely small in comparison with the energy density that one may estimate using any particle physics’ mass m from, say, the electronvolt scale to the mass scale of the weak gauge bosons in electroweak theory, $W^{\pm }$ and Z ($\sim 80,90$ GeV), the Higgs mass ($\sim 125$ GeV) and the top quark mass ($\sim 170$ GeV). The exception would be, of course, a millielectronvolt neutrino [17], but for any typical standard model particle the value of ${\rho }_{\mathrm{vac}}\propto m^4$ is mind-bogglingly too large, being indeed dozens of orders of magnitude astrayed as compared to the measured value of $\Lambda$, not to speak of the situation in the Grand Unified Theories (GUT’s), where the characteristic energy scale can reach $\sim {10}^{16}$ GeV. It is because of the cosmological constant problem phrased on these grounds that the VED option became outcast as if were to be blamed of all evils. The aforesaid notwithstanding, the criticisms usually have nothing better to offer, except to defend tooth and nail particular forms of the DE without providing any explanation about the genuine subject involved in the original discussion of this problem, which is, of course, to understand the role played by the VED in QFT and its fundamental relation with $\Lambda$. Most, if not all proposed forms of DE, are actually plagued with the same (purported) fine tuning illness that is attributed (in a way by fiat) to the vacuum option exclusively. This is certainly the case e.g. with the popular family of quintessence models, phantom fields and generalizations thereof, see e.g. [5,6,22,23] and references therein.

In recent years, however, new approaches to the notion of vacuum energy in QFT and its relation with the $\Lambda$-term suggest that these problems can be smoothed out to a large extent. In fact, the VED can be properly renormalized in QFT in curved spacetime, thereby offering a better theoretical context for the traditional vacuum energy approach to fit in with the observations. In the light of these developments, the quantum vacuum energy could well be after all the most fundamental explanation for the DE in our universe. See e.g. [17,19] as well as the latest formal developments in [24,25,26,27], summarized in [20].

The vanilla $\Lambda$CDM model, to which modern cosmological observations have converged in the last decades, is certainly an important triumph in our description of the main background features of the cosmic expansion and the large-scale structure formation processes in the universe. However, it is only a partial success. Its exceeding simplicity eventually turned into a perilous double-edge sword; in fact, the absence of any connection with fundamental physics is the literal expression of such a simplicity and is most likely at the root of its many shortages. In truth, the $\Lambda$CDM does not possess enough theoretical structure to explain the successfulness of the observations (e.g. the measured value of $\Lambda$) on a fundamental context, and at the same time it cannot even provide an explanation for other measurements that are threatening its viability. If we pay attention to the existing conflicts on several active fronts, the observational situation of the $\Lambda$CDM in the last decade or so does not seem to paint a fully rosy picture anymore. Beyond formal theoretical issues, a series of practical problems of more mundane nature than those mentioned above are piling up as well [28]. On a mere phenomenological perspective, it is particularly worrisome the situation with some “tensions” existing with the data. For example, it has long been known that there appear to exist potentially serious discrepancies between the CMB observations (based on the vanilla $\Lambda$CDM), and the local direct (distance ladder) measurements of the Hubble parameter today[29]. The persisting mismatch between these measurements is what has been called the “$H_0$-tension”. It is arguably the most puzzling open question within the current cosmological paradigm and it leads, if taken at face value, to a severe discrepancy of $\sim 5\sigma$ c.l. or more, between the mentioned observables. Many proposals have been put forward to shed some spark of light into that puzzling cosmological imbalance. Among the possibilities debated in the literature, it has been conjectured e.g. that it could stem from a possible intrinsic “running of $H_0\left(z\right)$ with the redshift” presumably connected with the differences that may appear in the (total) effective equation of state (EoS) of the universe between the vanilla cosmology and the actual FLRW model underlying the observations [30,31]. While these are an interesting possibilities, we are probably still far away from understanding the resolution of this conundrum on fundamental grounds. At the same time, there exists a smaller but appreciable ($\sim 2-3\sigma$) tension in the realm of the large-scale structure (LSS) growth data, called the “${\sigma }_8$-tension” [32]. It is concerned with the measurements of weak gravitational lensing at low redshifts ($z<1$). Such a tension is usually evaluated with the help of the parameter $S_8$ or, alternatively, by means of ${\sigma }_8$; recall that $S_8\equiv {\sigma }_8\sqrt{{\Omega }_m^0/0.3}$. It turns out that these measurements favor matter clustering weaker than that expected from the vanilla model using parameters determined by CMB measurements, see e.g. [33,34,35,36,37,38,39,40]. Recently, it has been claimed that $S_8$ values determined from $f{\sigma }_8$ increase with redshift in the $\Lambda$CDM[41], which, according to these authors, provides additional support to the fact that such a discrepancy may be physical in origin and with a value in the enhanced $2-4\sigma$ range. In the constant pursue for a possible late-time solution to these tensions, it has been argued that within the large class of models where the DE is treated as a fluid with EoS $w\left(z\right)$, solving the $H_0$ tension demands the phantom condition $w\left(z\right)<-1$ at some z, while solving both the $H_0$ and ${\sigma }_8$ tensions requires $w\left(z\right)$ to cross the phantom divide and/or other sorts of transitions, see e.g. [42,43,44,45,46,47]. Specific realizations of the noticed double condition for the DE fluid can be found in the literature, e.g. in the context of the $\Lambda$XCDM model [48,49], closely related to the idea of running vacuum to be discussed in the present work, see below. For detailed reviews on these tensions and other challenges afflicting the concordance $\Lambda$CDM model, see e.g.  [28,50,51] and the long list of references quoted there bearing relation to these matters.

The severity of some of these tensions and the huge number of proposals existing in the literature trying to explain them through a large disparity of ideas suggest that it is perhaps time to come to grips anew with the fundamentals of the theoretical formulations, such as quantum field theory and string theory. We have already pointed out to recent calculations claiming a more adequate renormalization prescription for the VED in quantum field theory in FLRW spacetime, and leading to the “running vacuum model” (RVM). It turns out that this QFT approach may have real impact not only on the more formal theoretical problems described in the beginning, but also on the practical issues concerning the aforementioned tensions. In fact, the resulting VED from the RVM leads to a time-varying vacuum energy density, and hence a time-varying (physical) $\Lambda$ as well, in which $\Lambda$ acquires a dynamical component through the quantum vacuum effects: $\Lambda \to \Lambda +\delta \Lambda$. The shift $\delta \Lambda$ is calculable in QFT since it depends on the contributions from the quantized matter fields (bosons and fermions). Upon appropriate renormalization one finds that $\delta \Lambda$ depends on a term of order of the Hubble rate $H\left(t\right)$ squared [24]: $\delta \Lambda \sim \nu \mathcal{O}\left(H^2\right)$ ($\left|\nu \right|\ll 1$). This is the typical form of the RVM. The connection of the latter with QFT can be motivated from semi-qualitative renormalization group arguments on scale-dependence, see the reviews [17,19]. In particular, we mention the old works [52,53] and also recent approaches along these lines, such as [54]. However, an explicit QFT calculation leading to that form of ‘running $\Lambda$’ (associated to the `running VED’) appeared only very recently [24,25,26,27].

A note of caution is in order here. Over time a large variety of cosmological models have been proposed to describe the DE and its possible dynamics. Apart from the aforesaid quintessence models and the like[5,6,22,23], there is a very populated habitat of models with generic time-dependent cosmological constant, the so-called “$\Lambda \left(t\right)$-cosmologies”. Many of these models, however, are of pure phenomenological nature since the time dependence of $\Lambda \left(t\right)$ is parameterized in an ad hoc manner. They might have a connection with fundamental theory, but it is not implemented in an explicit way in the corresponding papers. The list of models of this type is large and we will cite here only a few of them [55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74]; see also the old review [75]. In some cases, the parameterization is performed through a direct function of the cosmic time or of the scale factor, and sometimes as a function of the Hubble parameter, or even an hybrid combination of these various possibilities. Be as it may, the general and rather nonspecific class of the “$\Lambda \left(t\right)$-cosmologies” should not be confused with the “running vacuum models” (RVMs) discussed above, in which the running of $\Lambda$ stems from the quantum effects on the effective action of QFT in curved (FLRW) spacetime. In other words, the RVMs are to be understood in a much more restricted sense; in fact, one that is closer to fundamental aspects of QFT, and only this precise type of time-evolving VED cosmologies will be dealt with here. Let us finally note that, apart from the QFT formulation, a ‘stringy’ version of the RVMs is also available which can be very promising too [76,77,78,79]. The potential dynamics of the cosmic vacuum is, therefore, well motivated from different theoretical perspectives, and this fact further enhances the interest for the current study, whose main purpose is to focus exclusively on the phenomenological implications of the class of RVM models.

We should mention that the running vacuum framework has already been tested with considerable success in previous works over the years. It has been known for quite some time that the RVM-type of cosmological models can help in improving the overall fit to the cosmological observations and also in smoothing out the mentioned tensions as compared to the $\Lambda$CDM; see for instance [36,80,81,82,83,84,85,86,87,88], and [89,90] for a short summary. For this reason we believe it is worthwhile to keep on exploring the phenomenological consequences of the running vacuum in the light of the latest observations on all the main data sources: SNIa+BAO+$H\left(z\right)$+LSS+CMB. The state-of-the-art-phenomenological performance of the RVM has been reported not too long ago in [80]. In the current work, however, we definitely enhance the scope of the results presented in that paper by considering an updated cosmological data set in combination with an extended analysis of the CMB part. In point of fact, the main focus in this paper is to delve into the practical ability of the RVM to tackle the ${\sigma }_8$ and $H_0$ tensions versus the vanilla $\Lambda$CDM model. It is reassuring to find that the global fit to the cosmological observations can be improved within the running vacuum framework with respect to the $\Lambda$CDM. The optimal situation is when the VED presents a threshold in the recent past, where its dynamics becomes activated, and/or when the gravitational coupling is mildly running.

All in all, the dynamical DE models may offer a clue not only to relieve some high-brow aspects of the cosmological constant and coincidence problems, but also to straighten out some very practical ones, such as helping to modulate the processes of structure formation which may impinge positively on the ${\sigma }_8$-tension. Last but not least, they can help to explain the existing mismatch between the distinct values of $H_0$ derived from measurements of the local and the early universe.

The paper is organized as follows. In Section 2, we present the running vacuum model (RVM) from a phenomenological point of view and emphasize its connection with QFT in curved spacetime. For convenience, we introduce the model variant of the RVM which we call RRVM, as we did in [80]. In it, the VED can be expressed entirely in terms of the curvature scalar $\mathcal{R}$ (which is of order $\sim H^2$) at the background level and study two types of RRVMs: type I and type II, depending on whether the gravitational coupling G is fixed at its current local gravity value, $G=G_N$, or evolving mildly with the expansion, $G=G\left(H\right)$, a feature which in our case is respectively linked to the interaction or not of the evolving vacuum energy with cold dark matter (CDM). Type I is studied at length in Section 3, where we describe the background cosmological equations and its solution under appropriate conditions. At the same time we discuss the corresponding perturbations equations. Type II, on the other hand, is studied in detail in Section 4, where again we provide the background solution and analyze the perturbations. In Section 5 we enumerate and briefly describe the different sources of observational data employed in this paper and the methodology used to constrain the free parameters of the models under discussion. We also define the four characteristic datasets (Baseline and Baseline+SH0ES with and without CMB polarization data) that will be used to test the running vacuum models and their comparison with the vanilla $\Lambda$CDM model. The outcome of our analyses under the different datasets is presented and discussed in detail in Section 6. Finally, in Section 7 we summarize our findings and present the main conclusions of this study. In the Appendix A at the end of our work we include additional tables with a detailed breakdown of the different ${\chi }^2$ contributions from each observable.

As indicated, throughout our study we will assume that the background spacetime is FLRW with flat three-dimensional hypersurfaces. The general low-energy form of the vacuum energy density (VED) within the running vacuum model (RVM) has been explored phenomenologically on several previous occasions and with a remarkable degree of success, in the sense that in all cases it has proven to be rather competitive with the $\Lambda$CDM and even able to surpass the fitting performance of the latter, see e.g.[36,80,81,82,83,84,85,86,87,88]. Herein we shall test if it is still the case with the current wealth of observations and using the state-of-the-art methods of analysis of the cosmological data. The dynamical structure of the running VED adopts the perspicuous form[24,25]

Here the dot indicates derivative with respect to the cosmic time and $H=\dot{a}/a$ is the Hubble function. As we can see, the two leading dynamical terms of ${\rho }_{\mathrm{vac}}$ in Eq. (1) are both of $\mathcal{O}\left(H^2\right)$ since $\dot{H}\sim H^2$, this being true both in the matter and radiation dominated epochs. Despite that the higher order terms $\mathcal{O}\left(H^4\right)$ in the above expression are predicted as well in the QFT context [25], they are unimportant for the current universe and will be hereafter ignored. The additive parameter $c_0$ is constrained to satisfy ${\rho }_{\mathrm{vac}}\left(H_0\right)={\rho }_{\mathrm{vac}}^0$, where ${\rho }_{\mathrm{vac}}^0$ is the value of the VED today, and hence is connected with the physical value of the measured cosmological constant through ${\rho }_{\mathrm{vac}}^0={\Lambda }_{\mathrm{phys}}/\left(8\pi G_N\right)$. The bulk of this value is provided by $c_0$ as the two (dimensionless) coefficients $\nu$ and $\tilde{\nu }$ adjoined to the two dynamical terms in (1) are expected to be small ($|\nu ,\tilde{\nu }|\ll 1$)  [91]. They encode the running character of the vacuum at low energy and can be computed in QFT in curved spacetime from the contributions from the quantized bosons and fermion fields. The explicit calculation was first presented in [24,25] and was recently completed in [27]. In these references, it is shown that the above VED structure can be formally derived from quantum effects on the effective action of QFT in FLRW spacetime. In practice, however, the values of $\nu ,\tilde{\nu }$ must be fitted to the cosmological observations. Thus, in what follows we will focus exclusively on the phenomenological consequences of the RVM. What is important is that these coefficients are expected to be small and of order $\sim M_X^2/m_{\mathrm{Pl}}^2\ll 1$, where $m_{\mathrm{Pl}}\simeq 1.22\times {10}^{19}$ GeV is the Planck mass and $M_X\sim {10}^{16-17}$ GeV is of order of a typical GUT scale (or even a string scale slightly above it) times a multiplicity factor accounting for the number of heavy particles in the GUT [91]. For $\nu =\tilde{\nu }=0$ we recover the $\Lambda$CDM smoothly. This is a very welcome property of the RVM since DE models having no smooth $\Lambda$CDM limit, e.g. predicting a VED of the form ${\rho }_{\mathrm{vac}}\left(H\right)\propto H^2$ or a combination of $H^2$ and $\dot{H}$ (without any additive term), would be excluded owing to their absence of an inflexion point from deceleration into acceleration in the cosmic evolution, see [88,92,93]. The presence of the nonvanishing additive term $c_0$ is therefore crucial for the RVM to avoid this unwanted situation, something that other models (e.g. entropic and ghost models of the DE) cannot avoid and thereby get into trouble [94,95]. Holographic models with dynamical cutoff $L=H^{-1}$ also lack of an additive term in the DE and are also unfavored already at a pure cosmographic level [96]. In stark contrast, the condition $c_0\ne 0$ is always warranted within the class of the RVMs.

It is important to realize that the dynamics of the VED must preserve, of course, the Bianchi identity satisfied by the Einstein tensor. In practice this means that the total energy-momentum tensor (EMT), which receives the contributions from nonrelativistic matter, radiation and vacuum (assumed here to be ideal fluids), must be covariantly conserved, namely ${\nabla }^{\mu }{T}_{\mu \nu }^{\mathrm{tot}}=0$. The total EMT reads, where $U^{\mu }$ is the 4-velocity vector of the cosmic fluid. We have defined ${\rho }_t={\rho }_m+{\rho }_{\mathrm{ncdm}}+{\rho }_{\gamma }$, where ${\rho }_m={\rho }_{\mathrm{cdm}}+{\rho }_b$ denotes the contribution to the proper density of nonrelativistic matter from cold dark matter and baryons, ${\rho }_{\mathrm{ncdm}}$ (non-cold dark matter) corresponds to the energy density of neutrinos and ${\rho }_{\gamma }$ designates the energy density associated with photons. Analogous notations apply to the pressures. We shall, however, be more specific in our treatment of the various contributions to the EMT in the next section. Notice that in the above expression we have used $p_{\mathrm{vac}}=-{\rho }_{\mathrm{vac}}$ for the EoS of the vacuum fluid. Even though this condition may be violated slightly by quantum effects within a formal treatment of the subject in QFT [26], we shall nonetheless stick for now to the traditional EoS of the vacuum. We shall come back to this point later on. Upon expanding ${\nabla }^{\mu }{T}_{\mu \nu }^{\mathrm{tot}}=0$, it amounts to the local covariant conservation law in a FLRW universe where, in general, not only ${\rho }_{\mathrm{vac}}$ but also G may be functions of the cosmic time. This will depend on the particular implementation assumed for the matter sector. If we assume that there is an interaction of the VED with matter, then G can stay fixed at the usual value $G_N$ (the local gravity value), whereas if matter is locally conserved, then G must vary accordingly in order to preserve the covariant conservation law (3).

In order to ease the comparison with previous results, we shall adhere to the approach of [97] and assume $\tilde{\nu }=\nu /2$. In this way the RVM model is left with one single parameter and at the same time adopts the suggestive form in which $\mathcal{R}=12H^2+6\dot{H}$ is the curvature scalar1. That particular implementation is called, for obvious reasons, the RRVM, since it is a version of the RVM which involves the scalar of curvature [97]. One additional advantage is that it is automatically well-behaved in the radiation dominated epoch since in it $\mathcal{R}/H^2\ll 1$, and the standard BBN is not perturbed at all by the presence of vacuum energy. In the general case (1) such condition can also be fulfilled on assuming sufficiently small (absolute) values of $\nu ,\tilde{\nu }$ [85].

Finally, despite the general structure of the running VED is of the form (1), for convenience we define two types of RRVM scenarios. In type-I scenario the vacuum is in interaction with matter, whereas in type-II matter is conserved at the expense of an exchange between the vacuum and a slowly evolving gravitational coupling $G\left(H\right)$. The combined cosmological ‘running’ of these quantities insures the accomplishment of the Bianchi identity (and the associated local conservation law). In the following sections we study these two cases separately.

In this section, we consider the type-I RRVM scenario, in which the vacuum can be running at the expense of exchanging energy with matter. We will assume that only cold dark matter (CDM) is involved in such an exchange (therefore no baryons, neutrinos or photons are transferred to or from the vacuum). Whether it is the vacuum that generates new CDM or the CDM that disappears into the vacuum depends on the sign of the parameter $\nu$ in Eq. (4). For $\nu >0$, the vacuum decays into tiny amounts of CDM, whilst for $\nu <0$ some dark matter disappears into the vacuum. We do not presume which of these situations hold, we will fit the value (and sign) of $\nu$ to the cosmological data. This requires to solve the background and linear perturbation equations of the type-I running vacuum model, which we do in Secs. Section 3.1 and Section 3.2.

For type-II RRVM we have an entirely different sort of scenario, in which matter is strictly conserved, in particular dark matter, and hence no interaction of any sort is permitted between matter and vacuum. However, to switch off the energy exchange between matter and vacuum in a fully consistent way with the Bianchi identity, we must allow for the running of the gravitational coupling with the expansion, $G=G\left(H\right)$. Therefore, for type-II models we have both the running of the vacuum energy density ${\rho }_{\mathrm{vac}}$ and the running of G. Let us briefly see how it comes about. If matter is conserved we have ${\dot{\rho }}_t+3H({\rho }_t+p_t)=0$. The subscript t refers to the sum of all the species in the universe excluding the vacuum, as in previous sections. Whereupon the general Bianchi identity (3) boils down to4

Since the running of ${\rho }_{\mathrm{vac}}$ is still fixed by (5), the previous equation is essential to determine the running of G. Being ${\dot{\rho }}_{\mathrm{vac}}\propto \nu$, the sign of $\nu$ determines the sign of $\dot{G}$, i.e. if $\nu >0$ ($\nu <0$) G increases (decreases) with the cosmic expansion. Of course, we have to make sure that G evolves in a very mild way, which in fact turns out to be the case as we shall verify explicitly.

In what follows we bring forth the relevant background and linear perturbation equations for the type-II RRVM in Secs. Section 4.1 and Section 4.2, respectively. As it will become clear, the solution of this model type is more complicated.

We fit the $\Lambda$CDM model, the running vacuum models under consideration (the type-I RRVM, the type-I RRVM${}_{\mathrm{thr}.}$ and the type-II RRVM) and finally the XCDM [147] (also called wCDM), a generic parameterization of the dynamical DE whose dark energy EoS, $w_0$, is constant and is one parameter of the fit (expected to lie near $-1$). To test the response of the XCDM along with the relevant models under consideration can be useful, as it serves as a benchmark scenario for generic models of dynamical dark energy. We fit all these models to a large, robust and updated set of cosmological observations from all the main sources. Our data set involves observations from: i) distant type Ia Supernovae (SNIa); ii) baryonic acoustic oscillations (BAO); iii) a compilation of (differential age) measurements of the Hubble parameter at different redshifts ($H\left(z_i\right)$); iv) large-scale structure (LSS) formation data (specifically, an updated list of data points on the observable $f\left(z_i\right){\sigma }_8\left(z_i\right)$), and, finally, v) CMB Planck 2018 data of different sorts. A brief description now follows of each of these datasets along with the corresponding references:

SNIa: We consider the data from the so-called ‘Pantheon+’ compilation [148], which contains the apparent magnitudes and redshifts associated to 1701 light curves obtained from 1550 SNIa in the redshift range $0.001\le z\le 2.26$. See Sec. 2.2 of [149] for details of the theoretical formulae employed to take into account these data points. Interestingly, the new Pantheon+ compilation also includes the 77 light curves from the 42 SNIa in the host galaxies employed by the SH0ES team in their analysis [150,151]. The distance to the host galaxies has been measured with calibrated Cepheids. The inclusion of these luminosity distances in our dataset will be made clear by adding the label “+SH0ES”. They break the existing full degeneracy between $H_0$ and the absolute magnitude of SNIa, M, when only SNIa are considered in the analysis. The SH0ES calibration of the supernovae in conjunction with the cosmic distance ladder leads to larger preferred values of the Hubble parameter of $73.04\pm 1.04$ km/s/Mpc [150]. This large value as compared to Planck’s measurement ($67.36\pm 0.54$ km/s/Mpc, obtained from the TT,TE,EE+lowE+lensing data[152]), is at the root of the $\sim 5\sigma$$H_0$-tension.

BAO: We employ 13 data points on isotropic and anisotropic BAO estimators. See Table 1 to know the exact values and the corresponding references.

Cosmic chronometers: In our analyses we use 32 data points on the Hubble parameter $H\left(z_i\right)$ measured with the differential age technique [153]. They span the redshift range $0.07\le z\le 1.965$. We provide the complete list of data points and the corresponding references in Table 2. We have considered the effect of the known correlations between the various data points, as explained in [137]. See also Table 2 and its caption. The covariance matrix has been computed using the script provided in the following link 7.

LSS: 15 large-scale structure (LSS) data points between $0.01\lesssim z\lesssim 1.5$ embodied in the observable $f\left(z_i\right){\sigma }_8\left(z_i\right)$, which is known as the weighted linear growth rate, being $f\left(z\right)$ the so-called growth factor and ${\sigma }_8\left(z\right)$ the root mean square mass fluctuations on the $R_8=8h^{-1}$ Mpc scale. See Table 3 for the complete list of data points and the corresponding references. We can take advantage of the relation $f\left(z\right){\sigma }_8\left(z\right)=-(1+z)\frac{d{\sigma }_8\left(z\right)}{dz}$ to compute this quantity. The function ${\sigma }_8\left(z\right)$ involves the matter power spectrum, which is computed numerically by our modified version of the Einstein-Boltzmann code CLASS. It is important to note that this way of computing $f\left(z\right){\sigma }_8\left(z\right)$ can only be used provided that we are in the linear regime, since in this case and in our models the matter density contrast can be written as ${\delta }_m(a,k)=D\left(a\right)F\left(k\right)$, where the dependence on the scale factor and the comoving wave number k is factored out. The term $D\left(a\right)$ is known as the growth function and $F\left(k\right)$ encodes the initial conditions.

CMB: For the cosmic microwave background data, we utilize the full Planck 2018 TT,TE,EE+lowE likelihood [15]. It incorporates the information of the CMB temperature and polarization power spectra, and their cross-correlation. We refer to this dataset simply as “CMB”. We also test separately the Planck 2018 TT+lowE likelihood, which does not include the effect of the high-ℓ multipoles of the CMB polarization spectrum. This is done to check the impact of this particular dataset on our fitting results. It is also useful to compare with our previous analyses [80], in which only this type of CMB data were used. In our fitting scenarios, we indicate the removal of the high-ℓ CMB polarization data from the complete CMB likelihood with the label “CMB (No pol.)”.

As described in the preceding lines, for the SNIa we may or may not include the information provided by the SH0ES team whereas for the CMB we can consider the effect of the high-ℓ polarization data or not. An alternative calibration method of the absolute magnitude of SNIa based on the tip of the red giant branch [154,155] instead of Cepheids yields a measurement of $H_0$ somewhat in the middle of those provided by Planck [15] and SH0ES [150], $H_0=69.8\pm 1.7$ km/s/Mpc 8. In addition, it is also convenient to test the impact of the CMB polarization data from Planck, since previous works in the literature have found a moderate inconsistency between them and the Planck CMB temperature data, both in the 2015 [35] and 2018 [37] releases. This inconsistency could be due to a deficiency of the $\Lambda$CDM or the presence of unaccounted systematics in the data. Hence, these arguments motivate us to explore these four different datasets:

These are the four datasets that we employ to constrain our models. We present our fitting results in Table 4, Table 5, Table 6 and Table 7, Figure 1, Figure 2 and Figure 3, and also in the Table A1–Table A4 of Appendix A.