First evidence of running cosmic vacuum: challenging the concordance model

Despite the fact that a rigid $\Lambda$-term is a fundamental building block of the concordance $\Lambda$CDM model, we show that a large class of cosmological scenarios with dynamical vacuum energy density $\rho_{\Lambda}$ and/or gravitational coupling $G$, together with a possible non-conservation of matter, are capable of seriously challenging the traditional phenomenological success of the $\Lambda$CDM. In this paper, we discuss these"running vacuum models"(RVM's), in which $\rho_{\Lambda}=\rho_{\Lambda}(H)$ consists of a nonvanishing constant term and a series of powers of the Hubble rate. Such generic structure is potentially linked to the quantum field theoretical description of the expanding Universe. By performing an overall fit to the cosmological observables $SNIa+BAO+H(z)+LSS+BBN+CMB$ (in which the WMAP9, Planck 2013 and Planck 2015 data are taken into account), we find that the class of RVM's appears significantly more favored than the $\Lambda$CDM, namely at an unprecedented level of $\gtrsim4.2\sigma$. Furthermore, the Akaike and Bayesian information criteria confirm that the dynamical RVM's are strongly preferred as compared to the conventional rigid $\Lambda$-picture of the cosmic evolution.


INTRODUCTION
As of about twenty years ago dark energy (DE) has become an observational fact of the first magnitude in physics (Riess et al. 1998;Perlmutter et al. 1999) and the most recent observations do not cease to corroborate its existence as the prime cause for the acceleration of the Universe (Planck XIII 2015;Planck XIV 2015).
Next year it will be the centenary of the cosmological constant (CC) term, Λ, in Einstein's equations. The Λ-term is usually considered the simplest possible explanation for the DE and is an essential ingredient of the so-called concordance or ΛCDM model. On the theoretical side, however, the explanation is not so easy. In quantum field theory (QFT) the large value predicted for Λ, or equivalently for the corresponding vacuum energy density associated to it, ρ Λ = Λ/(8πG) (G being Newton's gravitational coupling), as compared to the measured one generates the old CC problem, for reviews see (Weinberg 1989;Sahni & Starobinsky 2000;Padmanabhan 2003;Peebles & Ratra 2003;Solà 2013). It is probably one of the most fundamental and unsolved conundrums of theoretical physics. The acclaimed finding of the Higgs boson of the standard model of particle physics at the LHC actually bolsters even more the problem since it adopts a more experimental basis. In fact, the associated electroweak (EW) vacuum energy density reads |ρ EW | ∼ M 2 H /G F , where M H ≃ 125 GeV is the measured Higgs boson mass and G F is Fermi's constant. Thus, in this most realistic situation, the CC problem appears on comparing the two quantitites ρ EW ∼ 10 9 GeV 4 to ρ Λ ∼ 10 −47 GeV 4 , which differ by an appalling amount of 56 orders of magnitude.
With such a state of affairs, cosmologists have felt sola@fqa.ub.edu adriagova@fqa.ub.edu decruz@fqa.ub.edu 1 Departament de Física Quàntica i Astrofísica, and Institute of Cosmos Sciences, Univ. de Barcelona, Av. Diagonal 647, E-08028 Barcelona, Catalonia, Spain motivated to look for many other sources of DE beyond Λ. For example, scalar fields in cosmology, φ, have been used since long ago, most conspicuously in the context of Brans-Dicke theories (Brans & Dicke 1961), where G ∝ 1/φ(t), and subsequently in general scalar-tensor theories. But soon also played a role as a strategy to endow the vacuum and the CC of some time dependence in a QFT context, Λ = Λ(φ(t)), and in some cases with the purpose to adjust dynamically its value. Some of the old approaches to the CC problem from the scalar field perspective are the works by Endo & Fukui (1977,1982; Fujii (1982); Dolgov (1983); Abbott (1985); Zee (1985); Barr (1987); Ford (1987); Weiss (1987). Among the proposed dynamical mechanisms, let us mention the cosmon model (Peccei, Solà & Wetterich 1987), which was subsequently discussed in detail by Weinberg (1989). In all cases, a more or less obvious form of fine tuning underlies the adjusting mechanisms. For this reason scalar fields were later used mostly to ascribe a possible evolution to the vacuum energy with the hope to explain the cosmic coincidence problem, giving rise to the notion of quintessence and the like, cf. (Peebles & Ratra 1988a,b;Wetterich 1988Wetterich , 1995Caldwell, Dave & Steinhardt 1998;Zlatev, Wang & Steinhardt 1999;Amendola 2000), among many other alternatives.
The old CC problem is a problem of fundamental nature that shows the profound interconnection among different branches of modern physics. Some of the above old works aimed at solving the problem at a time when it was thought that Λ = 0, so it was expected that some symmetry or some dynamical mechanism could help. But the task became much harder when it was realized that the CC value is nonvanishing and actually very small in particle physics units (ρ Λ ∼ 10 −47 GeV 4 ).
In this work we will not face the CC problem as such, not even the cosmic coincidence problem. Our main aim is much more modest. Taking into account the current amount and quality of the cosmological data on SNIa+BAO+H(z)+LSS+BBN+CMB, we wish to put to the test the possibility that the Λ-term and its associated vacuum energy density, ρ Λ = Λ/(8πG), could actually be dynamical ("running") quantities whose rhythms of variation might be linked to the Universe's expansion rate, H. The idea is to check if this possibility helps to improve the description of the overall cosmological data as compared to the rigid assumption Λ =const. inherent to the concordance ΛCDM model. For the class of models being considered we do not make any direct association of the Λ and G running with the dynamical evolution of scalar fields. The proposal being investigated here can be motivated in QFT in curved spacetime (cf. Solà (2013);  and references therein) and we want to show that it can be currently tested. Although a simple Lagrangian description of these models at the level of standard scalar fields is not available, attempts have been made in the literature (Solà 2008(Solà , 2013 and in any case this is of course something that one would eventually hope to find. There is, however, no guarantee that such description is possible in terms of a simple local action (Solà 2008).
Our main aim here is phenomenological. We will argue upon carefully confronting theory and observations that the idea of running vacuum models (RVM's) can be highly competitive, if not superior, to the traditional ΛCDM framework. The first serious indications of dynamical vacuum energy (at the ∼ 3σ c.l.) were reported in Solà, Gómez-Valent & de Cruz Pérez (2015). Earlier comprehensive studies hinted also at this possibility but remained at a lower level of significance, see e.g. Basilakos, Plionis & Solà (2009);Grande et al. (2011); Gómez-Valent, Solà & Basilakos (2015) 2 . Remarkably, in the present work the reported level of evidence is significantly higher than in any previous work in the literature (to the best of our knowledge). While Occam's razor says that "Among equally competing models describing the same observations, choose the simplest one", the point we wish to stress here is that the RVM's are able to describe the current observations better than the ΛCDM, not just alike. For this reason we wish to make a case for the RVM's, in the hope that they could shed also some new light on the CC problem, e.g. by motivating further theoretical studies on these models or related ones.

2
Recent claims that the ΛCDM may not be the best description of our Universe can also be found in e.g. The plan of the paper is as follows. In section 2 we describe the different types of running vacuum models (RVM's) that will be considered in this study. In section 3 we fit these models to a large set of cosmological data on distant type Ia supernovae (SNIa), baryonic acoustic oscillations (BAO's), the known values of the Hubble parameter at different redshift points, the large scale structure (LSS) formation data, the BBN bound on the Hubble rate, and, finally, the CMB distance priors from WMAP and Planck. We include also a fit of the data with the standard XCDM parametrization, which serves as a baseline for comparison. In section 4 we present a detailed discussion of our results, and finally in section 5 we deliver our conclusions.

TWO BASIC TYPES OF RVM'S
In an expanding Universe we may expect that the vacuum energy density and the gravitational coupling are functions of the cosmic time through the Hubble rate, thence ρ Λ = ρ Λ (H(t)) and G = G(H(t)). Adopting the canonical equation of state p Λ = −ρ Λ (H) also for the dynamical vacuum, the corresponding field equations in the Friedmann-Lemaître-Robertson-Walker (FLRW) metric in flat space become formally identical to those with strictly constant G and Λ: The equations of state for the densities of relativistic (ρ r ) and dust matter (ρ m ) read p r = (1/3)ρ r and p m = 0, respectively. Consider now the characteristic RVM structure of the dynamical vacuum energy: where G can be constant or a function G = G(H; ν, α) depending on the particular model. The above expression is the form that has been suggested in the literature from the quantum corrections of QFT in curved spacetime (cf. Solà (2013);  and references therein). The terms with higher powers of the Hubble rate have recently been used to describe inflation, see e.g. Lima, Basilakos & Solà (2013 and Solà (2015), but these terms play no role at present and will be hereafter omitted. The coefficients ν and α have been defined dimensionless. They are responsible for the running of ρ Λ (H) and G(H), and so for ν = α = 0 we recover the ΛCDM, with ρ Λ and G constants. The values of ν and α are naturally small in this context since they can be related to the β-functions of the running. An estimate in QFT indicates that they are of order 10 −3 at most (Solà 2008), but here we will treat them as free parameters of the RVM and hence we shall determine them phenomenologically by fitting the model to observations. As previously indicated, a simple Lagrangian language for these models that is comparable to the scalar field DE description may not be possible, as suggested by attempts involving the anomaly-induced action (Solà 2008(Solà , 2013. Two types of RVM will be considered here: i) type-G models, when matter is conserved and the running of ρ Λ (H) is compatible with the Bianchi identity at the expense of a (calculable) running of G; ii) type-A models, Note. -The best-fit values for the ΛCDM, XCDM and the RVM's, including their statistical significance (χ 2 -test and Akaike and Bayesian information criteria, AIC and BIC, see the text). The large and positive values of ∆AIC and ∆BIC strongly favor the dynamical DE options (RVM's and XCDM) against the ΛCDM (see text). We use 90 data points in our fit, to wit: 31 points from the JLA sample of SNIa, 11 from BAO, 30 from H(z), 13 from linear growth, 1 from BBN, and 4 from CMB (see S1-S7 in the text for references). In the XCDM model the EoS parameter ω is left free, whereas for the RVM's and ΛCDM is fixed at −1. The specific RVM fitting parameter is ν eff , see Eq. (6) and the text. For G1 and A1 models, ν eff = ν. The remaining parameters are the standard ones (h, ω b , ns, Ωm). The quoted number of degrees of freedom (dof ) is equal to the number of data points minus the number of independent fitting parameters (5 for the ΛCDM, 6 for the RVM's and the XCDM. The normalization parameter M introduced in the SNIa sector of the analysis is also left free in the fit, cf. Betoule et al. (2014), but it is not listed in the table). For the CMB data we have used the marginalized mean values and standard deviation for the parameters of the compressed likelihood for Planck 2015 TT,TE,EE + lowP data from Huang, Wang & Wang (2015), which provide tighter constraints to the CMB distance priors than those presented in Planck XIV (2015). Note. -Same as in Table 1, but excluding from our analysis the BAO and LSS data from WiggleZ, see point S5) in the text.
in contrast, denote those with G =const. in which the running of ρ Λ must be accompanied with a (calculable) anomalous conservation law of matter. Both situations are described by the generalized local conservation equation ∇ µ GT µν = 0, whereT µν = T µν +ρ Λ g µν is the total energy-momentum tensor involving both matter and vacuum energy. In the FLRW metric, and summing over all energy components, we find If G and ρ Λ are both constants, we recover the canonical conservation lawρ m +ρ r + 3Hρ m + 4Hρ r = 0 for the combined system of matter and radiation. For type-G models Eq. (4) boils down toĠ(ρ m + ρ r + ρ Λ ) + Gρ Λ = 0 sinceρ m + 3Hρ m = 0 andρ r + 4Hρ r = 0 for separated conservation of matter and radiation, as usually assumed. Mixed type of RVM scenarios are possible, but will not be considered here.
We can solve analytically the type-G and type-A models by inserting equation (3) into (1) and (2), or using one of the latter two and the corresponding conservation law (4). It is convenient to perform the integration using the scale factor a(t) rather than the cosmic time. For type-G models the full expression for the Hubble function normalized to its current value, E(a) = H(a)/H 0 , can be found to be where Ω i = ρ i0 /ρ c0 are the current cosmological parameters for matter and radiation, and we have defined Note that E(1) = 1, as it should. Moreover, for ξ, ξ ′ → 1 (i.e. |ν, α| ≪ 1) ν eff ≃ ν − α and ν ′ eff ≃ ν − (4/3)α. In the radiation-dominated epoch, the leading behavior of Eq. (5) is ∼ Ω r a −4ξ ′ , while in the matter-dominated epoch is ∼ Ω m a −3ξ . Furthermore, for ν, α → 0, E 2 (a) → 1 + Ω m (a −3 − 1) + Ω r (a −4 − 1). This is the ΛCDM form, as expected in that limit. Note that the following constraint applies among the parameters: c 0 = H 2 0 Ω Λ − ν + α Ω m + 4 3 Ω r , as the vacuum energy density ρ Λ (H) must reproduce the current value ρ Λ0 for H = H 0 , using Ω m + Ω r + Ω Λ = 1. The explicit scale factor dependence of the vacuum energy density, i.e. ρ Λ = ρ Λ (a), ensues upon inserting (5) into (3). In addition, since the matter is conserved for type-G models, we can use the obtained expression for ρ Λ (a) to also infer the explicit form for G = G(a) from (1). We refrain  30, 6.18, 11.81, 19.33, 27.65 (corresponding to 1σ, 2σ, 3σ, 4σ and 5σ c.l.) after marginalizing over the rest of the fitting parameters indicated in Table 1. We display the progression of the contour plots obtained for model G2 using the 90 data points on SNIa+BAO+H(z)+LSS+BBN+CMB, as we evolve from the high precision CMB data from WMAP9, Planck 2013 and Planck 2015 -see text, point S7). In the sequence, the prediction of the concordance model (ν eff = 0) appears increasingly more disfavored, at an exclusion c.l. that ranges from ∼ 2σ (for WMAP9), ∼ 3.5σ (for Planck 2013) and up to 4σ (for . Subsequent marginalization over Ωm increases slightly the c.l. and renders the fitting values indicated in Table 1, which reach a statistical significance of 4.2σ for all the RVM's. Using numerical integration we can estimate that ∼ 99.81% of the area of the 4σ contour for Planck 2015 satisfies ν eff > 0. We also estimate that ∼ 95.47% of the 5σ region also satisfies ν eff > 0. The corresponding AIC and BIC criteria (cf. Table 1) consistently imply a very strong support to the RVM's against the ΛCDM. from writing out these cumbersome expressions and we limit ourselves to quote some simplified forms. For instance, the expression for ρ Λ (a) when we can neglect the radiation contribution is simple enough: where ρ c0 = 3H 2 0 /8π G 0 is the current critical density and G 0 ≡ G(a = 1) is the current value of the gravitational coupling. Quite obviously for ξ = 1 we recover the ΛCDM form: ρ Λ = ρ c0 (1 − Ω m ) = ρ c0 Ω Λ =const. As for the gravitational coupling, it evolves logarithmically with the scale factor and hence changes very slowly 3 . It suffices to say that it behaves as is a smooth function of the scale factor. We can dispense with the full expression here, but let us mention that f (a) tends to one at present irrespective of the values of the various parameters Ω m , Ω r , ν, α involved in it; and f (a) → 1 in the remote past (a → 0) for ν, α → 0 (i.e. ξ, ξ ′ → 1). As expected, G(a) → G 0 for a → 1, and G(a) has a logarithmic evolution for ν ′ eff = 0. Notice that the limit a → 0 is relevant for the BBN (Big Bang Nucleosynthesis) epoch and therefore G(a) should not depart too much from G 0 according to the usual bounds on BBN. We shall carefully incorporate this restriction in our analysis of the RVM models, see later on.
Next we quote the solution for type-A models. As indicated, in this case we have an anomalous matter conservation law. Integrating (4) for G =const. and using (3) in it one finds ρ t (a) ≡ ρ m (a) + ρ r (a) = ρ m0 a −3ξ + ρ r0 a −4ξ ′ . We have assumed, as usual, that there is no exchange of energy between the relativistic and non-relativistic components. The standard expressions for matter and radiation energy densities are recovered for ξ, ξ ′ → 1. The normalized Hubble function for type-A models is simpler than for type-G ones. The full expression including both matter and radiation reads: (9) From it and the found expression for ρ t (a) we can immediately derive the corresponding ρ Λ (a): (10) Once more for ν, α → 0 (i.e. ξ, ξ ′ → 1) we recover the ΛCDM case, as easily checked. In particular one finds ρ Λ → ρ Λ0 =const. in this limit.

FITTING THE VACUUM MODELS TO THE DATA
In order to better handle the possibilities offered by the type-G and type-A models as to their dependence on the two specific vacuum parameters ν, α, we shall refer to model G1 (resp. A1) when we address type-G (resp. type-A) models with α = 0 in Eq. (3). In these cases ν eff = ν. When, instead, α = 0 we shall indicate them by G2 and A2, respectively. This classification scheme is used in Tables 1-2 and 5-7, and in Figs. 1-6. In the tables we are including also the XCDM (cf. Section 4) and the ΛCDM.
To this end, we fit the various models to the wealth of cosmological data compiled from distant type Ia supernovae (SNIa), baryonic acoustic oscillations (BAO's), the known values of the Hubble parameter at different redshift points, H(z i ), the large scale structure (LSS) formation data encoded in f (z i )σ 8 (z i ), the BBN bound on the Hubble rate, and, finally, the CMB distance priors from WMAP and Planck, with the corresponding correlation matrices in all the indicated cases. Specifically, we have used 90 data points (in some cases involving compressed data) from 7 different sources S1-S7, to wit: S1) The SNIa data points from the SDSS-II/SNLS3    Table 5 of that reference), including the correlations among these data encoded in the provided covariance matrices. We also use 2 data points based on D A (z i )/r s (z d ) and D H (z i )/r s (z d ) at z = 2.34, from the combined LyaF analysis Delubac et al. (2015). The correlation coefficient among these 2 points are taken from Aubourg et al. (2015) Chen, Kumar & Ratra (2016), where the authors make only use of Hubble parameter data in their analyses. We find, however, indispensable to take into account the remaining data sets to derive our conclusions on dynamical vacuum, specially the BAO, LSS and CMB observations. This fact can also be verified quite evidently in Figures 5-6, to which we shall turn our attention in Section 4. S5) f (z)σ 8 (z): 13 points. These are referred to in the text as LSS (large scale structure formation). The actual fitting results shown in Table 1 make use of the LSS data listed in Table 4, in which we have carefully avoided possible correlations among them (see below). Let us mention that although we are aware of the existence of other LSS data points in the literature concerning some of the used redshift values in our Table 4 -cf. e.g. Percival et al. (2004) 2014) -we have explicitly checked that their inclusion or not in our numerical fits has no significant impact on the main result of our paper, that is to say, it does not affect the attained 4σ level of evidence in favor of the RVM's. This result is definitely secured in both cases, but we have naturally presented our final results sticking to the most updated data.
The following observation is also in order. We have included both the WiggleZ and the CMASS data sets in our analysis. We are aware that there exists some overlap region between the CMASS and WiggleZ galaxy samples. But the two surveys have been produced independently and the studies on the existing correlations among these observational results (Beutler et al. 2016;Marín et al. 2016) show that the correlation is small. The overlap region of the CMASS and WiggleZ galaxy samples is actually not among the galaxies that the two surveys pick up, but between the region of the sky they explore. Moreover, despite almost all the WiggleZ region (5/6 parts of it) is inside the CMASS one, it only takes a very small fraction of the whole sky region covered by CMASS, since the latter is much larger than the WiggleZ one (see, e.g. Figure 1 in Beutler et al. (2016)). In this paper, the authors are able to quantify the correlation degree among the BAO constraints in CMASS and WiggleZ, and they conclude that it is less than 4%. Therefore, we find it justified to include the WiggleZ data in the main table of results of our analysis (Table 1), but we provide also the fitting results that are obtained when we remove the WiggleZ data points from the BAO and f (z)σ 8 (z) data sets (see Table 2). The difference is small and the central values of the fitting parameters and their uncertainties remain intact. Thus the statistical significance of Tables 1 and 2 is the same. S6) BBN: we have imposed the average bound on the possible variation of the BBN speed-up factor, defined as the ratio of the expansion rate predicted in a given model versus that of the ΛCDM model at the BBN epoch (z ∼ 10 9 ). This amounts to the limit |∆H 2 /H 2 Λ | < 10% (Uzan 2011). S7) CMB distance priors: R (shift parameter) and ℓ a (acoustic length) and their correlations with (ω b , n s ). For WMAP9 and Planck 2013 data we used the covariance matrix from the analysis of , while for Planck 2015 data those of Huang, Wang & Wang (2015). Our fitting results for the last case are recorded in all our tables (except in Table 5 where we test our fit in the absence of CMB distance priors R and ℓ a ). We display the final contour plots for all the cases, see Figs Notice that G1 and A1 have one single vacuum parameter (ν) whereas G2 and A2 have two (ν, α). There is nonetheless a natural alignment between ν and α for general type-G and A models, namely α = 3ν/4, as this entails ξ ′ = 1 (i.e. ν ′ eff = 0) in Eq. (6). Recall that for G2 models we have G(a) ∼ G 0 a 4(1−ξ ′ ) deep in the radiation epoch, cf. Eq. (8), and therefore the condition ξ ′ = 1 warrants G to take the same value as the current one, G = G 0 , at BBN. For model G1 this is not possible (for ν = 0) and we adopt the aforementioned |∆H 2 /H 2 Λ | < 10% bound. We apply the same BBN restrictions to the A1 and A2 models, which have constant G. With this setting all the vacuum models contribute only with one single additional parameter as compared to the ΛCDM: ν, for G1 and A1; and ν eff = ν − α = ν/4, for G2 and A2.
For the statistical analysis, we define the joint likelihood function as the product of the likelihoods for all the data sets. Correspondingly, for Gaussian errors the total χ 2 to be minimized reads: (11) Each one of these terms is defined in the standard way, for some more details see e.g. Gómez-Valent, Solà & Basilakos (2015), although we should emphasize that here the correlation matrices have been included. The BAO part was split as indicated in S2) and S3) above. Also, in contrast to the previous analysis of Solà, Gómez-Valent & de Cruz Pérez (2015), we did not use here the correlated Omh 2 (z i , z j ) diagnostic for H(z i ) data. Instead, we use As for the linear structure formation data we have computed the density contrast δ m = δρ m /ρ m for each vacuum model by adapting the cosmic perturbations formalism for type-G and type-A vacuum models. The matter perturbation, δ m , obeys a generalized equation which depends on the RVM type. For type-A models it reads (as a differential equation with respect to the cosmic time) where Ψ ≡ −ρ Λ ρm . For ρ Λ =const. we have Ψ = 0 and Eq. (13) reduces to the ΛCDM form 4 . For type-G models the matter perturbation equation is explicitly given in Solà, Gómez-Valent & de Cruz Pérez (2015). From here we can derive the weighted linear growth f (z)σ 8 (z), where f (z) = d ln δ m /d ln a is the growth factor and σ 8 (z) is the rms mass fluctuation amplitude on scales of R 8 = 8 h −1 Mpc at redshift z. It is computed from with W a top-hat smoothing function (see e.g. Gómez-Valent, Solà & Basilakos (2015) for details).
The linear matter power spectrum reads P (k, p) = P 0 k ns T 2 (k, p), where p = (h, ω b , n s , Ω m , ν eff ) is the fitting vector for the vacuum models we are analyzing (including the ΛCDM, for which ν eff = 0 of course), and T ( p, k) is the transfer function, which we take from Bardeen et al. (1986), upon introducing the baryon density effects through the modified shape parameter Γ (Peacock & Dodds 1994; Sugiyama 1995). We have also explicitly checked that the use of the effective shape of the transfer function provided in Eisenstein & Hu (1998) does not produce any change in our results. The expression (14) at z = 0 allows us to write σ 8 (0) in terms of the power spectrum normalization factor P 0 and the primary parameters that enter our fit for each model (cf . Table 1). We fix P 0 from in which we have introduced the vector of fiducial parameters p Λ = (h Λ , ω b,Λ , n s,Λ , Ω m,Λ , 0). This vector is defined in analogy with the fitting vector introduced before, but all its parameters are fixed and taken to be equal to those from the Planck 2015 TT,TE,EE+lowP+lensing analysis (Planck XIII 2015) with ν eff = 0. The fiducial parameter σ 8,Λ is also taken from the aforementioned Planck 2015 data. However, δ m,Λ (0) in (15) is computable: it is the value of δ m (z = 0) obtained from solving the perturbation equation of the ΛCDM using the mentioned fiducial values of the other parameters. Finally, from σ 8 (z) = σ 8 (0)δ m (z)/δ m (0) and plugging (15) in (14) one finds:   (Table 4) and the predicted curves by the RVM's, XCDM and the ΛCDM, using the best-fit values in Table 1. Shown are also the values of σ8(0) that we obtain for all the models. The theoretical prediction of all the RVM's are visually indistinguishable and they have been plotted using the same (blue) dashed curve.
the ΛCDM, all models become normalized to the same fiducial model defined above. The results for f (z)σ 8 (z) in the various cases are displayed in Fig. 4 together with the LSS data measurements (cf. Table 4). We will further comment on these results in the next section. Table 1 and Figures 1-2 present in a nutshell our main results. We observe that the effective vacuum parameter, ν eff , is neatly projected non null and positive for all the RVM's. The presence of this effect can be traced already in the old WMAP9 data (at ∼ 2σ), but as we can see it becomes strengthened at ∼ 3. It is also interesting to gauge the dynamical character of the DE by performing a fit to the overall data in terms of the well-known XCDM parametrization, in which the DE is mimicked through the density ρ X (a) = ρ X0 a −3(1+ω) associated to some generic entity X, which acts as an ersatz for the Λ term; ρ X0 being the current energy density value of X and therefore equivalent to ρ Λ0 , and ω is the (constant) equation of state (EoS) parameter for X. The XCDM trivially boils down to the rigid Λ-term for ω = −1, but by leaving ω free it proves a useful approach to roughly mimic a (non-interactive) DE scalar field with constant EoS. The corresponding fitting results are included in all our tables along with those for the RVM's and the ΛCDM. In Table 1 (our main table) and in Fig. 3, we can see that the best fit value for ω in the XCDM is ω = −0.916 ± 0.021. Remarkably, it departs from −1 by precisely 4σ.

DISCUSSION
Obviously, given the significance of the above result it is highly convenient to compare it with previous analyses of the XCDM reported by the Planck and BOSS collaborations. The Planck 2015 value for the EoS parameter of the XCDM reads ω = −1.019 +0.075 −0.080 (Planck XIII 2015) and the BOSS one is ω = −0.97 ± 0.05 (Aubourg et al. 2015). These results are perfectly compatible with our own result for ω shown in Table 1 for the XCDM, but in stark contrast to our result their errors are big enough Note. -Same as in Table 1, but removing both the R-shift parameter and the acoustic length la from our fitting analysis. Note. -Same as in Table 1, but removing the points from the LSS data set from our analysis, i.e. all the 13 points on f σ8 .
as to be also fully compatible with the ΛCDM value ω = −1. This is, however, not too surprising if we take into account that none of these analyses included LSS data in their fits, as explicitly indicated in their papers 5 . In the absence of LSS data we would find a similar situation. In fact, as our Table 6 clearly shows, the removal of the LSS data set in our fit induces a significant increase in the magnitude of the central value of the EoS parameter, as well as the corresponding error. This happens because the higher is |ω| the higher is the structure formation power predicted by the XCDM, and therefore the closer is such prediction with that of the ΛCDM (which is seen to predict too much power as compared to the data, see Fig. 4). In these conditions our analysis renders ω = −0.991±0.040, which is definitely closer to (and therefore compatible with) the central values obtained by Planck and BOSS teams. In addition, this result is now fully compatible with the ΛCDM, as in the Planck 2015 and BOSS cases, and all of them are unfavored by the LSS observations. This is consistent with the fact that both information criteria, ∆AIC and ∆BIC, become now slightly negative in Table 6, which reconfirms that if the LSS data are not used the ΛCDM performance is comparable or even better than the other models. So in order to fit the observed values of f σ 8 , which are generally lower than the predicted ones by the ΛCDM, |ω| should decrease. This is exactly what happens for the XCDM, as well as for the RVM's, when the LSS data are included in our analysis (in combination with the other data, particularly with BAO and CMB data). It is apparent from Fig. 4 that the curves for these models are then shifted below and hence adapt significantly better to the data 5 Furthermore, at the time these analyses appeared they could not have used the important LSS and BAO results from (Gil-Marín et al. 2016b), i.e. those that we have incorporated as part of our current data set, not even the previous ones from ( Gil-Marín et al. 2016a). The latter also carry a significant part of the dynamical DE signature we have found here, as we have checked.
points. Correspondingly, the quality of the fits increases dramatically, and this is also borne out by the large and positive values of ∆AIC and ∆BIC, both above 10 (cf. Table 1).
The above discussion explains why our analysis of the observations through the XCDM is sensitive to the dynamical feature of the DE, whereas the previous results in the literature are not. It also shows that the size of the effect found with such a parametrization of the DE essentially concurs with the strength of the dynamical vacuum signature found for the RVM's using exactly the same data. This is remarkable, and it was not obvious a priori since for some of our RVM's (specifically for A1 and A2) there is an interaction between vacuum and matter that triggers an anomalous conservation law, whereas for others (G1 and G2) we do not have such interaction (meaning that matter is conserved in them, thereby following the standard decay laws for relativistic and non-relativistic components). The interaction, when occurs, is however proportional to ν eff and thus is small because the fitted value of ν eff is small. This probably explains why the XCDM can succeed in nailing down the dynamical nature of the DE with a comparable performance. However not all dynamical vacuum models describe the data with the same efficiency, see e.g Salvatelli et al. (2014), Murgia, Gariazzo & Fornengo (2016), Li, Zhang & Zhang (2016). A detailed comparison is made among models similar (but different) from those addressed here in Solà, de Cruz Pérez, Gómez-Valent & Nunes (2016). In the XCDM case the departure from the ΛCDM takes the fashion of "effective quintessence", whereas for the RVM's it appears as genuine vacuum dynamics. In all cases, however, we find unmistakable signs of DE physics beyond the ΛCDM (cf. Table 1), and this is a most important result of our work.
As we have discussed in Section 2, for models A1 and A2 there is an interaction between vacuum and matter. Such interaction is, of course, small because the fitted values of ν eff are small, see Table 1. The obtained values are in the ballpark of ν eff ∼ O(10 −3 ) and therefore this is also the order of magnitude associated to the anomalous conservation law of matter. For example, for the nonrelativistic component we have This behavior has been used in the works by Fritzsch & Solà (2012 as a possible explanation for the hints on the time variation of the fundamental constants, such as coupling constants and particle masses, frequently considered in the literature. The current observational values for such time variation are actually compatible with the fitted values we have found here. This is an intriguing subject that is currently of high interest in the field, see e.g. Uzan (2011) and Solà, ed. (2015). For models G1 and G2, instead, the role played by ν eff and ν ′ eff is different. It does not produce any anomaly in the traditional matter conservation law (since matter and radiation are conserved for type-G models), but now it impinges a small (logarithmic) time evolution on G in the fashion sketched in Eq. (8). Thus we find, once more, a possible description for the potential variation of the fundamental constants, in this case G, along the lines of the above cited works, see also Fritzsch, Nunes & Solà (2016). There are, therefore, different phenomenological possibilities to test the RVM's considered here from various points of view.
We may reassess the quality fits obtained in this work from a different point of view. While the χ 2 min value of the overall fit for any RVM and the XCDM is seen to be definitely smaller than the ΛCDM one, it proves very useful to reconfirm our conclusions with the help of the time-honored Akaike and Bayesian information criteria, AIC and BIC, see Akaike (1974); Sugiura (1978); Schwarz (1978); Burnham & Anderson (2002). They read as follows: , BIC = χ 2 min +n ln N . (18) In both cases, n is the number of independent fitting parameters and N the number of data points used in the analysis. To test the effectiveness of a dynamical DE model (versus the ΛCDM) for describing the overall data, we evaluate the pairwise differences ∆AIC (∆BIC) with respect to the model that carries smaller value of AIC (BIC) -in this case, the RVM's or the XCDM. The larger these differences the higher is the evidence against the model with larger value of AIC (BIC) -the ΛCDM, in this case. For ∆AIC and/or ∆BIC in the range 6−10 one may claim "strong evidence" against such model; and, above 10, one speaks of "very strong evidence" (Akaike 1974;Burnham & Anderson 2002). The evidence ratio associated to rejection of the unfavored model is given by the ratio of Akaike weights, e ∆AIC/2 . Similarly, e ∆BIC/2 estimates the so-called Bayes factor, which gives the ratio of marginal likelihoods between the two models (Amendola 2015;Amendola & Tsujikawa 2015). Table 1 reveals conspicuously that the ΛCDM appears very strongly disfavored (according to the above statistical standards) as compared to the running vacuum models. Specifically, ∆AIC is in the range 17−18 and ∆BIC around 15 for all the RVM's. These results are fully con-sistent and since both ∆AIC and ∆BIC are well above 10 the verdict of the information criteria is conclusive. But there is another remarkable feature to single out at this point, namely the fact that the simple XCDM parametrization is now left behind as compared to the RVM's. While the corresponding XCDM values of ∆AIC and ∆BIC are also above 10 (reconfirming the ability of the XCDM to improve the ΛCDM fit) they stay roughly 4 points below the corresponding values for the RVM's. This is considered a significant difference from the point of view of the information criteria. Therefore, we conclude that the RVM's are significantly better than the XCDM in their ability to fit the data. In other words, the vacuum dynamics inherent to the RVM's seems to describe better the overall cosmological data than the effective quintessence behavior suggested by the XCDM parametrization.
Being the ratio of Akaike weights and Bayes factor much bigger for the RVM's than for the ΛCDM, the former appear definitely much more successful than the latter.
The current analysis undoubtedly reinforces the conclusions of our previous study (Solà, Gómez-Valent & de Cruz Pérez 2015), with the advantage that the determination of the vacuum parameters is here much more precise and therefore at a higher significance level. Let us stand out some of the most important differences with respect to that work: 1) To start with, we have used now a larger and fully updated set of cosmological data; 2) The selected data set is uncorrelated and has been obtained from independent analysis in the literature, see points S1-S7) above and references therein; 3) We have taken into account all the known covariance matrices among the data; 4) In this work, h, ω b and n s are not fixed a priori (as we did in the previous one), we have now allowed them to vary in the fitting process. This is, of course, not only a more standard procedure, but also a most advisable one in order to obtain unbiased results. The lack of consensus on the experimental value of h is the main reason why we have preferred to use an uninformative flat prior -in the technical sense -for this parameter. This should be more objective in these circumstances, rather than being subjectively elicited -once more in the technical sense -by any of these more or less fashionable camps for h that one finds in the literature, Riess et al.  (2016), which is ∼ 3σ larger than the former); 5) But the most salient feature perhaps, as compared to our previous study, is that we have introduced here a much more precise treatment of the CMB, in which we used not only the shift parameter, R, (which was the only CMB ingredient in our previous study) but the full data set indicated in S7) above, namely R together with ℓ a (acoustic length) and their correlations with (ω b , n s ).
Altogether, this explains the substantially improved accuracy obtained in the current fitted values of the ν eff parameter as compared to (Solà, Gómez-Valent & de Cruz Pérez 2015).
In particular, in what concerns points 1-3) above we should stress that for the present analysis we are using a much  more complete and restrictive BAO data set. Thus, while in our previous work we only used 6 BAO data points based on the A(z) estimator (cf. Table 3 of Blake et al. (2011b)), here we are using a total of 11 BAO points (none of them based on A(z), see S2-S3). These include the recent results from Gil- Marín et al. (2016b), which narrow down the allowed parameter space in a more efficient way, not only because the BAO data set is larger but also owing to the fact that each of the data points is individually more precise and the known correlation matrices have been taken into account. Altogether, we are able to significantly reduce the error bars with respect to the ones we had obtained in our previous work. We have actually performed a practical test to verify what would be the impact on the fitting quality of our analysis if we would remove the acoustic length l a from the CMB part of our data and replace the current BAO data points by those used in (Solà, Gómez-Valent & de Cruz Pérez 2015). Notice that the CMB part is now left essentially with the R-shift parameter only, which was indeed the old situation. The result is that we recover the error bars' size shown in the previous paper, which are ∼ 4 − 5 times larger than the current ones, i.e. of order O(10 −3 ). We have also checked what would be the effect on our fit if we would remove both the data on the shift parameter and on the acoustic length; or if we would remove only the data points on LSS. The results are presented in Tables 5 and 6, respectively. We observe that the ∆AIC and ∆BIC values become 2 − 4 points negative. This means that the full CMB and LSS data are individually very important for the quality of the fit and that without any of them the evidence of dynamical DE would be lost. If we would restore part of the CMB effect on the fit in Table 5 by including the R-shift parameter in the fitting procedure we can recover, approximately, the situation of our previous analysis, but not quite since the remaining data sources used now are more powerful. It is also interesting to explore what would have been the result of our fits if we would not have used our rather complete SNIa+BAO+H(z)+LSS+BBN+CMB data set and had restricted ourselves to the much more limited one used by the Planck 2015 collaboration in the paper Planck XIV (2015). The outcome is presented in Table 7. In contrast to Planck XIII (2015), where no LSS (RSD) data were used, the former reference uses some BAO and LSS data, but their fit is rather limited in scope since they use only 4 BAO data points, 1 AP (Alcock-Paczynski parameter) data point, and one single LSS point, namely f σ 8 at z = 0.57, see details in that paper. In contradistinction to them, in our case we used 11 BAO and 13 LSS data points, some of them very recent and of high precision (Gil-Marín et al. 2016b). From Table  7 it is seen that with only the data used in Planck XIV (2015) the fitting results for the RVM's are poor enough and cannot still detect clear traces of the vacuum dynamics. In particular, the ∆AIC and ∆BIC values in that table are moderately negative, showing that the ΛCDM does better with only these data. As stated before, not even the XCDM parametrization is able to detect any trace of dynamical DE with that limited data set, as the effective EoS is compatible with ω = −1 at roughly 1σ (ω = −0.960 ± 0.033). This should explain why the features that we are reporting here have been missed till now.
We complete our analysis by displaying in a graphical way the contributions from the different data sets to our final contour plots in Figs. 1-3. We start an-alyzing the RVM's case. For definiteness we concentrate on the rightmost plot for model A2 in Fig. 2, but we could do similarly for any other one in Figs 1-2. The result for model A2 is depicted in Fig. 5, where we can assess the detailed reconstruction of the final contours in terms of the partial contours from the different SNIa+BAO+H(z)+LSS+BBN+CMB data sources. This reconstruction is presented through a series of three plots made at different magnifications. In the third plot of the sequence we can easily appraise that the BAO+LSS+CMB data subset plays a fundamental role in narrowing down the final physical region of the (Ω m , ν eff ) parameter space, in which all the remaining parameters have been marginalized over. This reconstruction also explains in very obvious visual terms why the conclusions that we are presenting here hinge to a large extent on considering the most sensitive components of the data. While CMB obviously is a high precision component in the fit, we demonstrate in our study (both numerically and graphically) that the maximum power of the fit is achieved when it is combined with the wealth of BAO and LSS data points currently available.
In Fig. 6 we show the corresponding decomposition of the data contours for the XCDM model as well. In the upper-left plot we display the two-dimensional contours at 1σ and 2σ c.l. in the (Ω m , ω) plane, found using only the LSS data set. The elliptical shapes are obtained upon applying the Fisher matrix formalism (Amendola 2015), i.e. assuming that the two-dimensional distribution is normal (Gaussian) not only in the closer neighborhood of the best-fit values, but in all the parameter space. In order to obtain the dotted contours we have sampled the exact distribution making use of the Metropolis-Hastings Markov chain Monte Carlo algorithm (Metropolis et al. 1953;Hastings 1970). We find a significant deviation from the ideal perfectly Gaussian case. In the upper-right plot we do the same for the combination BAO+LSS. The continuous and dotted contours are both elliptical, which remarkably demonstrates the Gaussian behavior of the combined BAO+LSS distribution. Needless to say the correlations among BAO and LSS data (whose covariance matrices are known) are responsible for that, i.e. they explain why the product of the non-normal distribution obtained from the LSS data and the Gaussian BAO one produces perfectly elliptical dotted contours for the exact BAO+LSS combination. Similarly, in the lower-left plot we compare the exact (dotted) and Fisher's generated (continuous) lines for the CMB data. Again, it is apparent that the distribution inferred from the CMB data in the (Ω m , ω) plane is a multivariate normal. Finally, in the lower-right plot we produce the contours at 1σ and 2σ c.l. for all the data sets in order to study the impact of each one of them. They have all been found using the Fisher approximation, just to sketch the basic properties of the various data sets, despite we know that the exact result deviates from this approximation and therefore their intersection is not the final answer. The final contours (up to 5σ) obtained from the exact distributions can be seen in the small colored area around the center of the lower-right plot. The reason to plot it small at that scale is to give sufficient perspective to appreciate the contour lines of all the participating data. The final plot coincides, of course, with the one in Fig.3, where it can be appraised in full detail.
As it is clear from Fig. 6, the data on the H(z) and SNIa observables are not crucial for distilling the final dynamical DE effect, as they have a very low constraining power. This was also so for the RVM case. Once more the final contours are basically the result of the combination of the crucial triplet of BAO+LSS+CMB data (upon taking due care in this case of the deviations from normality of the LSS-inferred distribution). The main conclusion is essentially the same as for the corresponding RVM analysis of combined contours in Fig. 5, except that in the latter there are no significant deviations from the normal distribution behavior, as we have checked, and therefore all the contours in Fig. 5 can be accurately computed using the Fisher's matrix method.
The net outcome is that using either the XCDM or the RVM's the signal in favor of the DE dynamics is clearly pinned down and in both cases it is the result of the combination of all the data sets used in our detailed analysis, although to a large extent it is generated from the crucial BAO+LSS+CMB combination of data sets. In the absence of any of them the signal would get weakened, but when the three data sets are taken together they have enough power to capture the signal of dynamical DE at the remarkable level of ∼ 4σ.

CONCLUSIONS
To conclude, the running vacuum models emerge as serious alternative candidates for the description of the current state of the Universe in accelerated expansion. These models have a close connection with the possible quantum effects on the effective action of QFT in curved spacetime, cf. Solà (2013) and references therein. There were previous phenomenological studies that hinted in different degrees at the possibility that the RVM's could fit the data similarly as the ΛCDM, see e.g. the earlier works by Basilakos, Plionis & Solà (2009), Grande et al. (2011), Basilakos, Polarski & Solà (2012, Basilakos & Solà (2014), as well as the more recent ones by Gómez-Valent, Solà & Basilakos (2015) and , including of course the study that precedes this work, Solà, Gómez-Valent & de Cruz Pérez (2015). However, to our knowledge there is no devoted work comparable in scope to the one presented here for the running vacuum models under consideration. The significantly enhanced level of dynamical DE evidence attained with them is unprecedented, to the best of our knowledge, all the more if we take into account the diversified amount of data used. Our study employed for the first time the largest updated SNIa+BAO+H(z)+LSS+BBN+CMB data set of cosmological observations available in the literature. Some of these data (specially the BAO+LSS+CMB part) play a crucial role in the overall fit and are substantially responsible for the main effects reported here. Furthermore, recently the BAO+LSS components have been enriched by more accurate contributions, which have helped to further enhance the signs of the vacuum dynamics. At the end of the day it has been possible to improve the significance of the dynamical hints from a confidence level of roughly 3σ, as reported in our previous study (Solà, Gómez-Valent & de Cruz Pérez 2015), up to the 4.2σ achieved here. Overall, the signature of dynamical vacuum energy density seems to be rather firmly supported by the current cosmological observations. Al- contours in blue and purple are the exact ones, whilst the red and green ellipses have been obtained using the Fisher's approximation. Upper-right plot: Same, but for the combination BAO+LSS. Lower-left plot: As in the upper plots, but for the CMB data. Lower-right plot: The Fisher's generated contours at 1σ and 2σ c.l. for all the data sets: SNIa (dotted lines), H(z) (solid lines), BAO (dot-dashed lines), LSS (dotted, very thin lines) and CMB (solid lines, tightly packed in a very small, segment-shaped, region at such scale of the plot). The exact, final, combined contours (from 1σ up to 5σ) can be glimpsed in the small colored area around the center. See the text for further explanations and Fig. 3 for a detailed view. ready in terms of the generic XCDM parametrization we are able to exclude, for the first time, the absence of vacuum dynamics (ΛCDM) at 4σ c.l., but such limit can be even surpassed at the level of the RVM's and other related dynamical vacuum models, see  and the review Solà (2016).
It may be quite appropriate to mention at this point of our analysis the very recent study of us , in which we have considered the well-known Peebles & Ratra scalar field model with an inverse power law potential V (φ) ∝ φ −α (Peebles & Ratra 1988a,b), where the power α here should, of course, not be confused with a previous use of α for model A2 in Sect. 2). In that study we consider the response of the Peebles & Ratra model when fitted with the same data sets as those used in the current work. Even though there are other recent tests of that model, see e.g. the works by Samushia (2009) This explains why the analysis of Solà, Gómez-Valent & de Cruz Pérez (2016) was able to show that a non-trivial scalar field model, such as the Peebles & Ratra model, is able to fit the observations at a level comparable to the models studied here. In fact, the central value of the α parameter of the potential is found to be nonzero at ∼ 4σ c.l., and the corresponding equation of state parameter ω deviates consistently from −1 also at the 4σ level. These remarkable features are only at reach when the crucial triplet of BAO+LSS+CMB data are at work in the fitting analysis of the various cosmological models. The net outcome of these investigations is that several models and parametrizations of the DE do resonate with the conclusion that there is a significant (∼ 4σ) effect sitting in the current wealth of cosmological data. The effect looks robust enough and can be unveiled using a variety of independent frameworks. Needless to say, compelling statistical evidence conventionally starts at 5σ c.l. and so we will have to wait for updated observations to see if such level of significance can eventually be attained. In the meanwhile the possible dynamical character of the cosmic vacuum, as suggested by the present study, is pretty high and gives hope for an eventual solution of the old cosmological constant problem, perhaps the toughest problem of fundamental physics.