Running vacuum in quantum field theory in curved spacetime: renormalizing $\rho_{vac}$ without $\sim m^4$ terms

The $\Lambda$-term in Einstein's equations is a fundamental building block of the `concordance' $\Lambda$CDM model of cosmology. Even though the model is not free of fundamental problems, they have not been circumvented by any alternative dark energy proposal either. Here we stick to the $\Lambda$-term, but we contend that it can be a `running quantity' in quantum field theory (QFT) in curved spacetime. A plethora of phenomenological works have shown that this option can be highly competitive with the $\Lambda$CDM with a rigid cosmological term. The, so-called, `running vacuum models' (RVM's) are characterized by the vacuum energy density, $\rho_{vac}$, being a series of (even) powers of the Hubble parameter and its time derivatives. Such theoretical form has been motivated by general renormalization group arguments, which look plausible. Here we dwell further upon the origin of the RVM structure within QFT in FLRW spacetime. We compute the renormalized energy-momentum tensor with the help of the adiabatic regularization procedure and find that it leads essentially to the RVM form. This means that $\rho_{vac}(H)$ evolves as a constant term plus dynamical components ${\cal O}(H^2)$ and ${\cal O}(H^4)$, the latter being relevant for the early universe only. However, the renormalized $\rho_{vac}(H)$ does not carry dangerous terms proportional to the quartic power of the masses ($\sim m^4$) of the fields, these terms being a well-known source of exceedingly large contributions. At present, $\rho_{vac}(H)$ is dominated by the additive constant term accompanied by a mild dynamical component $\sim \nu H^2$ ($|\nu|\ll1$), which mimics quintessence.


Introduction
The cosmological constant (CC) term, Λ, in Einstein's equations has been for some three decades a fundamental building block of the 'concordance' or standard ΛCDM model of cosmology [1]. The model, however, was phenomenologically favored only as of the time Λ became a physically measured quantity some twenty years ago [2]. Nowadays Λ, or more precisely the associated current cosmological parameter Ω 0 Λ = ρ 0 vac /ρ 0 c became a precision quantity [3]. Here ρ 0 vac = Λ/(8πG N ) is the (vacuum) energy density induced by Λ, G N is Newton's constant and ρ 0 c = 3H 2 0 /(8πG N ) is the current critical density. The accurate knowledge of Ω 0 Λ around 0.7 is an important observational achievement, but it does not mean that we fully understand its nature and origin at a fundamental level. The cosmological constant problem [4,5] is a preeminent example of a fundamental theoretical conundrum, which actually affects all forms of dark energy (DE) [6][7][8][9][10]. The abstruse theoretical problems, though, are not the only nagging ones afflicting the concordance model. In practice the ΛCDM appears to be currently in tension with some important measurements, most significantly the discordant values of the current Hubble parameter H 0 obtained independently from measurements of the local and the early universe [11]. Whether these tensions are the result of as yet unknown systematic errors is not known, but there remains perfectly upright the possibility that a deviation from the ΛCDM model could provide an explanation for such discrepancies [12]. As it has been shown in the literature, models mimicking a time-evolving Λ (and hence a dynamical vacuum energy density ρ Λ ) could help in alleviating these problems, see e.g. [13][14][15][16][17][18] and [19][20][21][22][23][24][25].
In this paper, we would like to further dwell upon the theoretical possibility of having a dynamical vacuum energy density (VED), ρ vac , in the framework of quantum field theory (QFT) in curved spacetime [26][27][28][29]. Above all we wish to focus on the dynamics associated to the running vacuum model (RVM) [30][31][32]; for a review, see [33][34][35] and references therein. For related studies, see e.g. [36,37] and [38][39][40], some of them extending the subject to the context of supersymmetric theories [41,42] and also to supergravity [43]. More recently the matter has also been addressed successfully in the framework of the effective action of string theories [44,45]. Here, however, we aim at the computation of the VED in QFT in a curved background, specifically in the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric. We proceed by renormalizing the energy-momentum tensor using the adiabatic regularization prescription (ARP) [26,27]. This renormalization method is based on the WKB approximation of the field modes in the expanding universe. We perform the calculation in two ways, one through a modified form of the ARP [40] and the second (presented in one appendix) involving dimensional regularization (DR). The common result is that the properly renormalized VED, obtained upon inclusion of the renormalized value of ρ Λ at a given scale, does not contain the unwanted contributions proportional to the fourth power of the particle masses (∼ m 4 ) and hence it is free from large induced corrections to the VED. This is tantamount to subtracting the Minkowskian contribution from the curved spacetime result, as we show. In addition, we find that the final expression for the VED adopts the RVM form for the current universe, namely it contains not only the usual constant term but also one that evolves with the square of the Hubble rate (∼ νH 2 , with |ν| ≪ 1). The latter represents only a mild (dynamical) correction to the constant contribution and it can mimic quintessence or phantom DE depending on the sign of ν.
The structure of the article is as follows. In Sec. 2 we define our framework, which consists of a neutral scalar field non-minimally coupled to gravity, and compute the classical energy-momentum tensor (EMT). In Sections 3 and 4 we address the quantum fluctuations in the adiabatic vacuum through the WKB expansion of the field modes in the FLRW background. We discuss the adiabatic regularization of the EMT. In Sec. 5 we proceed to renormalize the EMT in the FLRW context using the adiabatic prescription, which is then needed in Sec. 6 to extract the precise form of ρ vac from the renormalized zero-point energy (ZPE) up to terms of adiabatic order 4, which in our case means up to O(H 4 ). We show that the relation between values of the VED at different scales is free from quartic powers of the masses. We also demonstrate that our renormalization procedure gives the same result as subtracting the Minkowskian contribution from the curved spacetime result. In Sec. 7, we provide the connection of the computed VED in this work with the running vacuum model (RVM), which had been derived before from the general point of view of the renormalization group in curved spacetime. The final discussion and a summary of the conclusions is presented in Sec. 8. Three appendices at the end furnish complementary material. Specifically, Appendix A defines our conventions and collects some useful formulas. Appendix B reconsiders the main parts of the renormalization of the EMT using dimensional regularization and the standard counterterm procedure, starting of course from the same WKB expansion of the field modes. Finally, Appendix C discusses alternative identifications of the VED leading to generalized forms of the RVM which had already been anticipated from the renormalization group approach in previous works. where T matter µν is the EMT of matter. They can conveniently be rewritten as where ρ Λ ≡ Λ/(8πG N ) is the VED associated to Λ. The latter contributes a term T Λ µν ≡ −ρ Λ g µν to the total EMT. However, in general, there will be also other contributions to the total VED, in particular those associated to the quantum fluctuations of the fields, and also to their classical ground state energy (if it is nonvanishing). For simplicity we will suppose that there is only one (matter) field contribution to the EMT on the right hand side of (2.2) in the form of a real scalar field, φ, and such contribution will be denoted T φ µν . Hence the total EMT reads T tot µν = T Λ µν + T φ µν . We neglect the incoherent matter contributions (e.g. from dust and radiation) for the kind of QFT considerations made in this study, as they can be added without altering the QFT aspects. 1 A list of geometric quantities of interest here are shown in the Appendix A, where we also define our conventions.
Suppose that the scalar field is non-minimally coupled to gravity and that it does not couple to itself 2 . The part of the action involving φ, then, reads where ξ is the non-minimal coupling of φ to gravity. In the special case ξ = 1/6, the massless (m = 0) action is conformally invariant, i.e. symmetric under simultaneous local Weyl rescalings of the metric and the scalar field: g µν → e 2α(x) g µν and φ → e −α(x) φ, for any local spacetime function α(x). However, we will keep ξ general since our scalar field will be massive. The field φ obeys the Klein-Gordon (KG) equation In the case of general non-minimal coupling ξ, the EMT can be computed upon straightforward calculation: (2.5) In the following, we are going to consider the spatially flat FLRW metric in the conformal frame. Introducing the conformal time, η, we have ds 2 = a 2 (η)η µν dx µ dx ν , with η µν = diag(−1, +1, +1, +1) the Minkowski metric. We will denote the derivative with respect to the conformal time by ′ ≡ d/dη. The corresponding Hubble rate is H(η) ≡ a ′ /a. Since dt = adη, the usual Hubble rate with respect to the cosmic time, H(t) =ȧ/a (with˙≡ d/dt), is related to the former through H(η) = aH(t). We will present most of our calculations in terms of the conformal time, but at the end it will be useful to express the VED in terms of the usual Hubble rate H(t), as this will ease the comparison with the RVM results in the literature. Because our metric is conformally flat, g µν = a 2 (η)η µν , we have the inverse g µν = a −2 (η)η µν and √ −g = a 4 (η), and as a result the action (2.3) can be rewritten If we perform the field redefinition φ = ϕ/a and disregard total derivatives, the previous action becomes the following functional of ϕ: where we have used (cf. Appendix A) R = 6a ′′ /a 3 . The above field redefinition enables us to have a simpler field equation for ϕ as if we were in Minkowski space (with conformal time) and an effective time-dependent mass different from that in (2.4). Computing δS[ϕ]/δϕ = 0 from (2.7) we find: 8) 2 We will not consider a possible contribution from a classical potential for φ in our analysis, which in general should also involve quantum corrections and hence leading to an effective potential. Here we wish to concentrate mainly on the zero-point energy of the quantum fields, which in itself is already rather cumbersome.
where˜ ϕ ≡ η µν ∂ µ ∂ ν ϕ = −ϕ ′′ + ∇ 2 ϕ is the Minkowskian box operator acting on ϕ in conformal coordinates x µ = (η, x). The above equation for ϕ is, of course, equivalent to (2.4) for the original field φ, as one can check by computing the curved spacetime box operator in the conformal metric.

Quantum fluctuations and adiabatic vacuum
Let us now move from classical to quantum field theory. We can take into account the quantum fluctuations of the field φ by considering the expansion of the field around its background (or classical mean field) value φ b : We identify the vacuum expectation value (VEV) of the field with the background value, that is to say, 0|φ(η, x)|0 = φ b (η), whereas we assume zero VEV for the fluctuation: 0|δφ|0 = 0 (as it will become evident from the mode expansion in terms of creation and annihilation operators to be discussed below). We will define the vacuum state to which we are referring to with more precision below. Given the above field decomposition into a classical plus fluctuating part, the corresponding EMT decomposes itself as µν is the contribution from the classical or background part, whereas T δ φ µν is the contribution from the quantum fluctuations. The 00-component of the latter is connected with the zero-point energy (ZPE) density of the scalar field in the FLRW background. Thus, the total vacuum contribution to the EMT reads The above equation says that the total vacuum EMT is made out of the contributions from the cosmological term and of the quantum fluctuations of the field. We will use later on a renormalized version of this equation and extract a relation satisfied by the renormalized VED. The field (3.1) obeys the curved spacetime KG equation (2.4) independently by the classical and quantum parts. Similarly, ϕ and δϕ obey separately the Minkowskian KG equation (2.8). Let us concentrate on the fluctuation δϕ. We can decompose it in Fourier frequency modes h k (η): Since φ = ϕ/a, the expansion of δφ is, of course, the same as that of (3.3) but divided by the scale factor a. Here A k and A † k are the (time-independent) annihilation and creation operators. Their commutation relations are the usual ones Notice that A k and h k have mass dimensions −3/2 and −1/2 in natural units, respectively. Upon substituting the Fourier expansion (3.3) in (˜ − m 2 eff (η))δϕ = 0 we find that the frequency modes of the fluctuations satisfy the (linear) differential equation with ω 2 k (m) ≡ k 2 +a 2 m 2 . As we can see, h k depends only on the modulus k ≡ |k| of the momentum. Because Ω k (η) is a nontrivial function of the conformal time, the modes cannot be found in a simple form. However, one can generate an approximate solution from a recursive self-consistent iteration based on the phase integral ansatz The latter is normalized through the Wronskian condition h ′ k h * k − h k h * ′ k = i, which insures that the standard equal-time commutation relations between the field operator ϕ and its canonical momentum, π ϕ = ϕ ′ , are preserved. The function W k in the above ansatz is solution of the differential equation obtained from inserting (3.6) into (3.5): Although this equation is non-linear, it can be solved using the WKB approximation. Taking into account that the WKB solution is valid for large k (i.e. for short wave lengths, as e.g. in geometrical Optics) the function Ω k is slowly varying for weak fields. This motivates a notion of vacuum called the adiabatic vacuum [48]. Rather than formulating it as the state without particles, we can at least say it is a state essentially empty of high frequency modes. Indeed, particles with definite frequencies cannot be defined in a curved background, since Ω k (η) is a function of time. Nonetheless a Fock space interpretation is still possible, and the adiabatic vacuum can be formally defined as the quantum state which is annihilated by all the operators A k of the above Fourier expansion, see [26][27][28][29] for details. Our VEV's actually refer to that adiabatic vacuum. In such conditions, the minimal excited state is h k ≃ e ikη / √ 2k, with k ≃ Ω k , and hence one can maintain an approximate particle interpretation of the quantized fields in a curved background provided the geometry is slowly varying. However, in general, the physical interpretation of the modes (3.5) with time varying frequencies must be sought in terms of field observables rather than in particle language. Thus, for a more physical interpretation of the vacuum effects of the expanding universe, we must compute the renormalized EMT in the FLRW background, and the first thing to do in order to carry out this task is to (adiabatically) regularize it.

Adiabatic Regularization of the Energy-Momentum tensor
For an adiabatic (slowly varying) Ω k , we can use Eq. (3.7) as a recurrence relation to generate an (asymptotic) series solution. In the gravitational context, such WKB approximation is organized through adiabatic orders and constitutes the basis for the adiabatic regularization procedure (ARP) 3 . The quantities that are taken to be of adiabatic order 0 are: k 2 and a. Of adiabatic order 1 are: a ′ and H. Of adiabatic order 2: a ′′ , a ′2 , H ′ and H 2 . Each additional derivative increases the adiabatic order by one unit. Therefore, the solution of the "effective frequency" W k is found from a WKB-type asymptotic expansion in powers of the adiabatic order: k is an adiabatic correction of order j. In this way we obtain an adiabatic expansion of the mode functions h k in powers of even order adiabatic terms (0, 2, 4, ...), such as a, a ′′ ∝ R, ω ′2 , ω ′′ , ω ′′2 , R 2 etc. The non-appearance of odd adiabatic orders is justified by arguments of general covariance, which forbid tensors of odd adiabatic order in the field equations.

Relating different renormalization scales through the ARP
We start by defining the first term ω (0) k of the above WKB expansion and compute its first two derivatives: Notice that in this approach the WKB expansion is performed off-shell, i.e. we use the arbitrary mass scale M instead of the original mass m (both being parameters of adiabatic order zero).
In this fashion the ARP can be formulated in such a way that we can relate the adiabatically renormalized theory at two scales [40]. The mass scale M can play a role similar to the scale µ in DR, but it can be given a more physical meaning. When M is fixed at the physical mass of the quantized field (M = m) we expect to obtain the renormalized theory on-shell. By keeping the M -dependence we can subtract the EMT at such value, thus obtaining the renormalized theory at M . In the subtraction procedure, the divergences will be cancelled and the quadratic mass differences ∆ 2 ≡ m 2 − M 2 will appear in the correction terms relating the theory at the two renormalization scales. These differences must be reckoned as being of adiabatic order 2 since they appear in the WKB expansion together with other terms of the same adiabatic order [40]. For ∆ = 0 we recover M = m and corresponds to the usual ARP (where one renormalizes the theory only at the scale of the particle mass) [26,27]. We will use this procedure to explore the behavior of the VED throughout the cosmological evolution. The masses m could be associated to fields of the Standard Model of particle physics, but as we shall see it will be convenient to consider also the heavy fields of some Grand Unified Theory (GUT) and explore the behavior in the low energy domain M 2 ≪ m 2 . Needless to say, for the sake of simplicity, we model here all particles in terms of (real) scalar fields. We can see from Eq.(4.2) that the expansion in terms of an even number of derivatives of ω k (hence of even adiabatic order) is equivalent to an expansion in even powers of H and odd powers of H ′ (notice e.g. that H 2 and H ′ are homogeneous), as in both cases it involves an even number of derivatives of the scale factor. In this way the expansion is compatible with general covariance, as indicated above. For the current universe, the powers H 2 and H ′ are sufficient for the phenomenological description, as it is obvious from the fact that R = (6/a 2 )(H 2 + H ′ ), whereas the higher powers bring corrections which can be important in the early universe.

Computing the adiabatic orders and the regularized ZPE
To obtain the different orders, we start with the initial solution W k ≈ ω (0) k indicated in Eq.(4.2). For a = 1 this would yield the standard Minkowski space modes. But since a = a(η) we have to find a better approximation. Introducing that initial solution on the RHS of (3.7) and expanding it in powers of ω −1 k we may collect the new terms up to adiabatic order 2 to find ω (2) k . Next we iterate the procedure by introducing W k ≈ ω k on the RHS of the same equation, expand again in ω −1 k and collect the terms of adiabatic order 4, etc. Since this mathematical procedure implies an expansion in powers of ω −1 k ∼ 1/k ∼ λ (i.e. a short wavelength expansion) it is obvious that the UV divergent terms of the ARP are the ones containing the first lowest powers of 1/ω k , and hence are concentrated in the first adiabatic orders, whilst the higher adiabatic orders represent finite contributions [46][47][48]. The result is intuitive: for any given physical quantity, the UV divergences are concentrated in the first adiabatic orders whereas the higher orders must decay sufficiently quick at high momentum so as to make the corresponding integrals convergent and yielding a suppressed contribution. Although not involved in our calculation, if we take the electric current, for instance, the divergences are concentrated up to 3rd order, since here one has to include a vector potential with adiabatic order 1. In contrast, for the main quantity at stake in our case, the EMT, its regularization implies to work up to 4th adiabatic order, as we shall show in detail below -cf. Eq. (4.9). Upon renormalization, we will obtain a finite expression for the EMT.
Using (4.2) and working out the second and fourth order adiabatic terms of (4.1), one finds We are now ready to compute the energy density associated to the quantum vacuum fluctuations in curved spacetime with FLRW metric. We start from the EMT given in Eq. (2.5) with φ decomposed as in (3.1). However, we are interested just on the fluctuating part, and select the quadratic fluctuations in δφ only since, as previously indicated, we have zero VEV for the fluctuation itself. This follows from (3.3) and the definition of adiabatic vacuum, implying that the crossed terms ∝ δφ vanish. For the 00-component, related to the energy density of the vacuum fluctuations, we find To clarify the notation, notice that (δφ We may now substitute the Fourier expansion of δφ = δϕ/a, as given in (3.3), into Eq. (4.4) and apply the commutation relations (3.4). After symmetrizing the operator field product δφδφ ′ with respect to the creation and annihilation operators, we end up with the following expression in terms of the amplitudes of the Fourier modes of the scalar field: where we have integrated´d 3 k (2π) 3 (...) over solid angles and expressed the final integration in terms of k = |k|. The different terms of the above integral should be expanded up to 4th order in adiabatic expansion using the WKB approximations (4.3): Upon substituting the above WKB expansions in (4.5) and using the relations (4.3) and (4.2), the result can be phrased as follows after a significant amount of algebra: Let us note the presence of the ∆-dependent terms in the last two rows, which contribute at second (∆ 2 ) and fourth (∆ 4 ) adiabatic order.

Particular cases: ZPE with minimal coupling and in Minkowski spacetime
As a particular case of the cumbersome expression obtained above, let us consider what is left when the non-minimal coupling to gravity is absent (ξ = 0). Let us also fix the scale M at the physical mass of the particle (M = m), so that the ∆-terms vanish. Finally, let us project the UV-divergent terms of order H 2 and neglect those of higher adiabatic order. It is then easy to check that Eq. (20) boils down to the very simple expression where we recall that ω k (m) ≡ √ k 2 + a 2 m 2 . Formula (4.10) is in agreement with previous results found in the literature for ξ = 0, in the O(H 2 ) approximation [46] -see also [37,41]. Notice that k is the comoving momentum, whereas the physical momentum isk = k/a. Defining the physical energy modeω k (m) = k2 + m 2 , and keeping in mind that H = aH, we can re-express the above result as .
The dependence on the scale factor can be eliminated as soon as we write T 00 = −ρ vac g 00 = a 2 ρ vac and rephrase the above result in terms of ρ vac . The last quantity can be thought of as representing the VED associated to the quantum fluctuations, i.e. the aforementioned ZPE. We will, however, come back to this point later on. At the moment we note that with this interpretation we can retrieve also the very particular situation of Minkowskian spacetime, as indicated above. Setting a = 1 (hence H = 0) the previous expression maximally simplifies to where the last integral is just the well-known ZPE of the quantum field φ in flat spacetime [33,51], as it should be expected (in natural units, = 1). In what follows, unless stated otherwise, we will continue using comoving momenta, as it will ease the presentation. The previous formulas correspond to simpler situations, but as we have seen the ZPE in FLRW spacetime is much more complicated, and Eq. (4.9) constitutes a WKB approximation to it up to 4th adiabatic order. Of course, all the above forms of ZPE are UV divergent and require renormalization.

Renormalization of the ZPE in the FLRW background
Let us consider the ZPE part of the EMT, as given by Eq. (4.9). We can split it into two parts as follows: is the UV-divergent contribution, which involves ω k = √ k 2 + a 2 M 2 and the powers 1/ω n k up to n = 3.
The terms in (4.9) which are not in (5.2) are the ones which are finite (as they involve powers of 1/ω k higher than 3), and constitute the T δφ 00 N on−Div (M ) part of (5.1). Computing the (manifestly convergent) integrals with the help of Eq. (B.2) (for ǫ = 0) in Appendix B, the final result reads where the dots in the last expression correspond to higher adiabatic orders. Let us now take a closer look to the divergent part of the ZPE, Eq. (5.2). Since the complete adiabatic series is an asymptotic series representation of Eq. (4.4), there is some arbitrariness in the way of choosing the leading adiabatic order because, independently of our choice, such series does not really converge and only serves as an approximation, which is obtained after one cuts the series at some particular order. There is, however, a minimum number of steps to do in order to obtain a meaningful result.
To start with, let us set the arbitrary scale M to the on-shell mass value of the quantized scalar field, M = m, hence ∆ = 0 (cf. Sec.4.1). In such a case, the divergent part (5.2) reduces to Again, (5.4) is a bare integral, formally divergent and does not depend on any renormalization scale. The prescription we are going to follow in order to renormalize the ZPE (and, in general, the EMT) is somehow reminiscent of the momentum subtraction scheme, although is certainly different in many respects. In the latter the renormalized Green's functions and running couplings are obtained by subtracting their values at a renormalization point p 2 = M 2 (space-like in our metric, which becomes an Euclidean point after Wick rotation) or at the time-like one p 2 = −M 2 (depending on the kinematical region involved) [52,53]. Since for vacuum diagrams we do not have external momenta, here, instead, we renormalize the ZPE by subtracting the terms that appear up to 4th adiabatic order at the arbitrary mass scale M . This suffices to eliminate the divergent terms through the ARP, as it is amply discussed in the literature [26][27][28].

Renormalized ZPE off-shell
In view of the previous considerations, we will define the renormalized ZPE in curved spacetime at the scale M as follows: where we have used the fact that T δφ 00 N on−Div (m) − T δφ N on−Div (M ) yields precisely the last term of (5.5), as it follows immediately from Eq. (5.3). In these expressions, (0 − 4) indicates the expansion up to fourth adiabatic order and the dots in (5.5) denote finite terms of higher adiabatic order. Using now Eq. (5.4), we arrive at the result For better clarity, we will henceforth distinguish explicitly between the off-shell energy mode ω k (M ) = √ k 2 + a 2 M 2 (formerly denoted just as ω k ) and the on-shell one ω k (m) = √ k 2 + a 2 m 2 . On using simple manipulations, such as e.g.
etc. one can work out the renormalized result (5.6) into the following convenient form in which all the integrals are seen to be manifestly convergent since the power counting for all of them leads to ∼´dkk −3 in the UV region. The calculation of some of these convergent integrals can be a bit cumbersome, as not all of them can be dealt directly with Eq. (B.2). Owing to various cancellations, however, the final result can be cast in a rather compact form: We have checked this result with the help of Mathematica [54]. The obtained expression vanishes for M = m, which was already obvious from (5.6) or (5.8), since the integrand is proportional to various powers of ∆ and to expressions that cancel in that limit. This is also clear from the definition itself, Eq. (5.5). However, it should be emphasized that the vanishing result in the M = m limit occurs only because we have computed the on-shell value T δφ µν Ren (m) also up to adiabatic order 4 in Eq. (5.9). In general one can compute T δφ µν Ren (m) up to any desired adiabatic order, keeping however in mind the asymptotic character of the WKB series. But in all cases the subtracted term in Eq. (5.5) at the arbitrary scale M is always to be computed to adiabatic order 4, as this suffices to cancel all the existing divergences. Beyond 4th order one always obtains finite, subleading, corrections. These higher order finite effects satisfy the Appelquist-Carazzone decoupling theorem [55] since they become suppressed for large values of the physical mass m of the quantum field. In our study, however, we do not track these finite, subleading, contributions, but of course they are there and provide a nonvanishing on-shell value of the renormalized EMT as defined by (5.5) .
Noteworthy, the final renormalized ZPE in curved spacetime (5.9), although it is perfectly finite, still carries at this point quartic powers of the masses.
We have explicitly checked that the above direct subtraction procedure gives the same result as the conventional DR technique applied to the divergent integrals of (5.2), see Appendix B for a summary of that alternative calculation. Of course, the DR is used here only as an auxiliary tool to regularize the UV divergences by tracking the poles up to adiabatic order four, but we do not mean at all to renormalize the calculation through the minimal subtraction (MS) scheme [52,53]. In fact, as we have demonstrated, the above result can be fully obtained without any use of DR, if it is not desired. The truly guiding renormalization principle here is the one based on the ARP relating different scales, with or without the auxiliary use of DR in the intermediate steps.

Renormalized vacuum energy density
We remind the reader that in order to make possible the renormalization program in the context of QFT in curved spacetime, we need to count on the higher derivative (HD) terms in the classical effective action of vacuum [26], in addition to the usual Einstein-Hilbert (EH) term with a cosmological constant, Λ. In four dimensions, the HD part is composed of the O(R 2 ) terms, i.e. the squares of the curvature and Ricci tensors: R 2 and R µν R µν . No more HD terms are needed in our case since the one associated to the square of the Riemann tensor, R µνρσ R µνρσ , is not independent owing to the topological nature of the Euler's density in 4 dimensions, which involves all these HD terms together: E = R µνρσ R µνρσ − 4R µν R µν + R 2 . Moreover the square of the Weyl tensor, C 2 = R µνρσ R µνρσ − 2R µν R µν + (1/3) R 2 , exactly vanishes for conformally flat spacetimes such as FLRW. The full action, therefore, boils down to the relevant EH+HD terms mentioned above plus the matter part (consisting here of the scalar field φ only) with a non-minimal coupling to gravity, Eq. (2.3). Variation of the action with respect to the metric provides the modified Einstein's equations, which become extended as compared to their original form (2.2) as follows: where we use renormalized quantities and hence we have indicated explicitly the dependence of the parameters and of the EMT on the subtraction point M . The background part does not depend on it. The higher order tensor H µν is obtained by functionally differentiating R 2 with respect to the metric (see Appendix A). A further simplification is possible here since the corresponding term associated to the functional differentiation of the square of the Ricci tensor, called H µν , is not necessary since it is not independent of H (1) µν for FLRW spacetimes [26]. This follows from the aforementioned properties of the Euler density and the Weyl tensor for conformally flat spacetimes. The higher order tensor H (1) µν is indeed to be included in the extended field equations since it is anyway generated by the quantum fluctuations and is therefore indispensable for the renormalizability of the theory. The fact that Eq. (6.1) has been written with all couplings defined at some arbitrary mass scale M is because we have shown that the EMT used in our calculation is the renormalized one at that scale following the ARP. However, in the Appendix B we offer an alternative approach based on the more conventional counterterm procedure, starting from the bare parameters of the action.
Baring in mind that we wish to relate the theory at different renormalization points 4 , let us subtract the modified Einstein's equations (6.1) at the two scales M and M 0 . The classical (background) contribution T φ b µν cancels in the difference, since as noted it does not depend on the renormalization scale, and we find for the various couplings involved X = G −1 N , ρ Λ , a 1 . Using now the tensor pattern shown by the generalized field equations (6.1), and taking into account that we know the expression for the final renormalized form of the EMT within the ARP, namely Eq (5.9), we can derive by comparison the renormalization shift (or 'running') undergone by the couplings G −1 N , ρ Λ and a 1 in (6.2) between the two scales M and M 0 . Such identification is possible since we know the explicit expressions for G 00 and H 6.1 Vacuum energy density at different scales. Absence of ∼ m 4 terms.
Following our discussion in Sec. 4.2, let us provisionally define the vacuum state as that one satisfying p vac = −ρ vac and T vac µν = −ρ vac g µν . We will further discuss the significance of this identification in Appendix C. Equating the last expression to Eq. (3.2), and taking the 00-component of the equality (keeping also in mind that g 00 = −a 2 (η) in the conformal frame), we obtain Notice that we have included the dependence on the renormalization point since we are using the renormalized theory at that scale. The above equation says that the total VED at an arbitrary scale M is the sum of the renormalized contributions from the cosmological term plus that of the quantum fluctuations of the scalar field at that scale (i.e. the renormalized ZPE). Subtracting the renormalized result at two scales, M and M 0 , and using (6.2), we find: where the term f ρ Λ (m, M, M 0 ) has cancelled, and we have used the expressions for G 00 and H (1) 00 given in the Appendix A. From equations (6.4) and (6.6,) we finally obtain where, in addition, we have used Eq. (A.5) to re-express the final result in terms of the ordinary Hubble function in cosmic time (H = H/a) since it will be useful for further considerations. The result (6.9) is the value of the VED at the scale M once we know its value at another scale M 0 , i.e. it expresses the 'running' of the VED. Only in the case of conformally invariant fields (ξ = 1/6) the result would be the same at all scales, if the VED would receive only contributions from scalar fields. But in general, this is not the case since one has to add the contribution from fermions and vector boson fields, which we do not consider here, so in general the total VED appears a running quantity with the expansion. The running is slow for small H, as it depends on terms of the form O(H 2 ) times a mass scale squared, and on O(H 4 ) contributions, but not on quartic mass scales.

Equivalent approach: subtracting the Minkowskian contribution
It cannot be overstated that the above result (6.9) is free from quartic powers of the masses. These would still be present if we had subtracted just the ZPE part at different scales without including the renormalized ρ Λ . This is obvious from Eq.(5.9), where we can see that the problem actually stems from Minkowskian spacetime, see [33] for a discussion. The renormalized ZPE in flat spacetime is obtained from Eq.(5.9) in the limit a = 1 (which implies that H and all its derivatives are zero). Only the first term of it remains, although it is the one carrying the mentioned quartic powers. This term vanishes for M = m since the renormalized on-shell value was computed only up to fourth adiabatic order. As previously emphasized (cf. Sec. 5.1), this does not mean that the exact renormalized ZPE vanishes on-shell, of course. One still has to add the higher order adiabatic terms, but they are finite and subleading since they decouple for larger and larger values of the physical mass m (i.e. they satisfy the decoupling theorem [55]), and we have not tracked them explicitly. Our main aim here was to pick out just the leading contributions to the renormalized ZPE up to 4th adiabatic order. Because we compute the total VED, defined as the sum of the renormalized value of ρ Λ and the renormalized ZPE, the difference of VED values at two scales is free from the quartic powers of mass scales. Of course this is possible owing to the renormalized form for the ZPE that we have used, Eq. (5.5), which involves a subtraction of the on-shell value at another arbitrary mass scale. In the Appendix B, we provide an alternative calculation leading to the same result (6.9) and further comments on this fact.
The above observations suggests that we can recover the expression (6.9) for the VED by performing an analogy with the Casimir effect; that is to say, we may compute the expression for T δφ 00 in Minkowskian spacetime and substract it from its equivalent in curved spacetime. One should expect that the result appears only mildly evolving with the cosmic evolution through a function of the Hubble rate (which is the key term providing the departure of the FLRW background from Minkowskian spacetime) [33]. In fact, the subtraction of the Minkowskian spacetime result has been argued from different perspectives [37,41]. In the Minkowski limit, the subtraction of scales in Eq.
The result is indeed the same as in Eq. (6.8), and hence we end up once more with the formula (6.9) for the total VED after we cast H and its derivatives in terms of the ordinary Hubble rate, H.
In other words, we can reach again the same relation between the values of VED at two different scales, which does not involve ∼ m 4 contributions.

Running vacuum connection
As previously remarked in a footnote, the result we were aiming at and which is represented now by Eq. (6.9) does not provide the calculated value of the vacuum energy density at a given scale, e.g. it says nothing on the value of ρ vac (M 0 ) and hence it has no implication on the cosmological constant problem mentioned in the Introduction, that is to say, it has no bearing on it if such problem is meant to be the computation of the value itself of the VED at some point in the history of the universe. However, our result can be useful to explore the 'running' of the VED when we move from one scale to another. In other words, if ρ vac is known at some scale M 0 , we can use the obtained relation to compute the value of ρ vac at another scale M . Such connection of values from one renormalization point to another is what we have been calling "running" of the VED, and in fact it was suggested long ago from the point of view of the renormalization group in curved spacetime from different perspectives [30][31][32] -for a review of the running vacuum model (RVM), see [33][34][35] and references therein. Interestingly enough, it can provide also a framework for the possible time variation of the so-called fundamental constants of nature [56]. In order to illustrate a possible interpretation of Eq. (6.9) along the lines of the RVM, let us assume that we define the renormalized VED at some Grand Unified Theory (GUT) scale M 0 = M X , where typically M X ∼ 10 16 GeV is associated also with the inflationary scale. It is natural to assume that the fundamental parameters of cosmology, such as e.g. ρ vac , are primarily defined at that scale, which appeared from the very beginning in the history of the universe. By choosing a GUT scale we also insure that most matter fields can be active degrees of freedom to some extent.

RVM in the current universe
Let ρ vac (M X ) be the value of the VED at the GUT scale M 0 = M X . Despite the fact that ρ vac (M X ) is unknown, it can be related to the current value of the VED, ρ 0 vac , through the relation ρ vac (M = H 0 ) = ρ 0 vac , in which we choose the second scale M at today's numerical value of the Hubble parameter, H 0 . This quantity can be used as an estimate for the energy scale of the background gravitational field associated to the FLRW universe at present. Notice that this is precisely the kind of association originally made in the aforementioned references on the RVM [33]. Therefore, from (6.9) applied to the current universe, we find the connection between the vacuum densities at the two points: where we have neglected all terms of order O(H 4 ) (which include alsoḢ 2 , HḦ and H 2Ḣ ) for the present universe (H = H 0 ). This equation can be used to find out the unknown value of ρ vac (M X ), and can be conveniently written as where we have defined the 'running parameter' for the VED: For ξ = 1/6, it vanishes and there is no running, as expected from conformal invariance. For ξ = 0 and m 2 /M 2 Remarkably, for general ξ the structure obtained for ν eff is very close to that obtained within the RVM approach, see [33]. In such context, it defines the coefficient of the one-loop β-function for the renormalization group equation of ρ vac . The presence of the additional logarithmic piece ln H 2 0 /M 2 X appears in the direct QFT calculation employed here, but it does not make any difference in practice since it is constant and ν eff must be fitted directly to the observations as an effective coefficient. In our case we have simplified the theoretical calculation by considering just the contribution from one single scalar field to ν eff . We expect it to be small, i.e. |ν eff | ≪ 1, owing to the ratio M 2 X /M 2 P ∼ 10 −6 . However the final value could be much larger since ν eff depends on ξ and also on the multiplicity and nature (fermion/boson) of the fields involved, so we cannot predict ν eff with precision on mere theoretical grounds. It must be confronted against observations. Notice that the standard model particles make no significant contribution, since for all of them m 2 /M 2 X ≪ 1. Only particles near or of order of the GUT scale may contribute significantly. The accurate determination of ν eff can only be obtained by fitting the RVM to the overall cosmological data, as it has been done in detail e.g. in [15,16], where it has been found to be positive and of order 10 −3 . Substituting (7.2) into Eq. (6.9) and limiting ourselves once more to the late universe (where all terms of O(H 4 ) can be neglected), we can estimate the VED near our time by taking M of order of the energy scale defined by the numerical value of H around the current epoch: which matches the exact canonical form of the RVM formula [33]. Let us note that such approximation holds reasonably good even if we explore the CMB epoch since the departure of ν eff (H) from ν eff is less than 8% (for m ≃ M X ) or much less if m ≪ M X . The above equation is the one which has been used to fit the value of ν eff in a variety of works, such as e.g. [15,16]. As a matter of fact, such works have considered a more general form as well, in which a term proportional tȯ H is also present in the running equation for the VED. Such term can appear under conditions that are discussed in Appendix C.

Implications for the early universe: RVM-inflation
While the previous discussion obviously applies to the current universe only, since we have neglected the O(H 4 ) terms on the r.h.s. of Eq. (6.9), we should emphasize that these terms can play a significant role in the early universe. These terms are actually generated from the functional differentiation of the R 2 -term in the higher derivative part of the vacuum action (cf. Appendix A), and therefore they play a similar role as in the case of Starobinsky's inflation [57,58] In fact, Starobinsky's inflation is not actually triggered by an early epoch in which H =const. and is actually characterized by a final phase with rapid oscillations of the gravitational field, which is when the universe leaves the inflationary phase and enters the radiation epoch after a reheating period, see e.g. [34] for a summarized discussion. However, as indicated in the introduction, the effective generation of terms proportional to H 4 in the early universe is perfectly possible from string-inspired mechanisms, see [44,45], in which the ∼ H 4 power is generated in the early universe from the vacuum average of the (anomalous) gravitational Chern-Simons term ∼ M P α ′ b(x)R µνρσ (x) R µνρσ (x), which is characteristic of the bosonic part of the low-energy effective action of the string gravitational multiplet. Here b(x) is the Kalb-Ramond axion field and α ′ is the slope parameter (M s = 1/α ′ being the string scale). An effective metastable vacuum is conceivable in this context since such state can be sustained until the universe transits into the radiation phase, and this occurs only after the gravitational anomalies are cancelled. This must indeed happen because matter (relativistic and nonrelativistic particles) cannot coexist with gravitational anomalies. These can actually be cancelled by the chiral anomalies of matter itself, see [44,45] for details. Before such thing occurs, a metastable de Sitter period remains temporally active and can bring about inflation through the (anomalygenerated) H 4 term. The type of inflation produced by the H 4 -term -and, in general, by higher order (even) powers of H -is characteristic of RVM-inflation. The latter follows a different pattern as compared to Starobinsky's inflation, but graceful exit is still granted -see [59][60][61] for details and particularly [34] for a comparison with Starobinsky's inflation 5 .
It seems clear that the presence of the higher powers of the Hubble rate in the early universe can be very important from different perspectives. For example, as noted in [44], they could help eschewing the possible trouble of string theories with the 'swampland' criteria on the impossibility to construct metastable de Sitter vacua in the string framework [62], which if so it would forbid the existence of de Sitter solutions in a low energy effective theory of quantum gravity. The existence of the H 4 -terms does not depend on picking out a particular potential for the scalar field since, as we should recall here, no potential has been introduced at any time in our framework nor in that of [44,45]. Thus, the RVM string inflation approach could provide a loophole to the swampland no-go criterion applied to fundamental scalar fields. But, of course, to fully establish it requires of a detailed investigation in the context of string-induced RVM [44,45], which is certainly not the subject of the present paper. What is phenomenologically relevant, though, is that once these terms are available they can be used to build up a generalized form of VED, which reads as follows: in which c 0 is a constant of dimension +2 in natural units, closely related to Λ; H I is a dimension +1 scale related to inflation; and ν and α are dimensionless coefficients, the former being obviously related to ν eff from the previous section. Such extended expression for the VED involving both ∼ H 2 and ∼ H 4 terms can produce successful inflation with graceful exit in the early universe [34,[59][60][61] and leaves an effective form of dynamical VED for the present universe behaving as (7.6). Remarkably, that form has been positively confronted with the data [13][14][15][16][17][18]. From our direct QFT calculation, we have seen that the ∼ H 2 terms indeed apear (see also Appendix C for a more general case) whereas the higher order terms that we have obtained are more along the Starobinsky inflationary line. However, we cannot exclude the presence of the ∼ H 4 string-induced effective contributions, as discussed in [44,45]. Being these contributions nonvanishing for H = H I =const. and taking into account that the Starobinsky-like higher order terms just vanish in such regime, it is reasonable to expect that for large values of H I the ∼ H 4 terms (if available from string-induced origin) prove to be the dominant terms at the inflationary scale. If so, this could change dramatically our picture of inflation into a more RVM-like one.

Discussion and conclusions
We have devoted this paper to investigate the possible dynamics of vacuum in the context of quantum field theory in the expanding universe, and more specifically in FLRW spacetime. The quantum field theoretical context is well-known [26][27][28][29] but the difficulties are still of formidable magnitude. This is obviously so since we know that in this kind of business sooner or later we have to face a huge stumbling block on our way, which is the cosmological constant problem [4,5]. Such mystery is perhaps the greatest conceptual challenge faced by theoretical physics ever, owing to the mind boggling discrepancy existing between the measured value of the vacuum energy density (VED) and the typically predicted one by our most cherished QFT's, say quantum chromodynamics and specially the electroweak standard model, both being essential parts of what we call the standard model of particle physics, which in itself is a highly successful theory of the fundamental interactions. Even though tackling such problems may require the concepts and the sophisticated theoretical tools underlying quantum gravity and string theory [5], difficulties appear indeed in all fronts, and string theory might not be an exception. Indeed, as of some time we known that string theories somehow abhor de Sitter space, as 'swampland' conjectures point to the impossibility to construct metastable de Sitter vacua in such theories [62]. We remain simply agnostic about these problems but, if true, they add up more trouble to the list of conundrums that fundamental physics has to face when addressing the physics of vacuum in an expanding universe. In the meantime, we expect that some sort of provisional result should perhaps be possible within the -much more pedestrian -semiclassical QFT approach, in which quantum matter fields interact with an external gravitational field.
Specifically, in this work we have reconsidered the calculation of the renormalized energymomentum tensor (EMT) of a real quantum scalar field non-minimally coupled to the FLRW background. We have performed the calculation following two lines of approach based on adiabatic regularization and renormalization of the EMT. In both cases we started from the WKB expansion of the field modes in the FLRW spacetime. Then we defined an appropriately renormalized EMT by performing a substraction of its on-shell value (i.e. the value defined at the mass m of the quantized field) at an arbitrary renormalization point M . The resulting EMT becomes finite because we subtract the first four adiabatic orders (the only ones that can be divergent). Since the renormalized EMT becomes a function of the arbitrary scale M , we can compare the renormalized result at different epochs of the cosmic history characterized by different energy scales. In one of the approaches (presented in the main text) we have shown by direct calculation that the renormalized EMT defined in that way is finite. In another approach (left for Appendix B) we use dimensional regularization to subtract the poles of the low adiabatic orders. Here we use the more conventional method based on cancelling the poles using the counterterms associated to the fundamental parameters ρ Λ , G −1 N and a 1 (the coefficient of R 2 ). The two approaches concur to the same renormalized result. The next important point is the extraction of the VED from the renormalized EMT, which is composed not only of the zero-point energy part (involving the quantum fluctuations of the scalar field) but also of ρ Λ (M ), the renormalized value of ρ Λ at the scale M . Remarkably, the sum of these two quantities is free from quartic terms ∼ m 4 , which are usually responsible for the exceedingly large contributions to the VED.
We have also shown that the renormalized VED obtained from this QFT calculation takes on approximately the usual form of the running vacuum models (RVM's) [33], in which ρ vac = ρ vac (H) appears in the manner of an additive constant plus a series of powers of H (the Hubble rate) and its time derivatives. Originally, the RVM approach was motivated from general considerations involving the renormalization group in QFT in curved spacetime (cf. [33] and references therein). At the end of the day, we have been able to show that the RVM form of the VED for the current universe can be achieved from direct calculations of QFT in the FLRW spacetime. In it, all the terms made out of powers of H (and its time derivatives) are of even adiabatic order. This means that all these powers effectively carry an even number of time derivatives of the scale factor, which is essential to preserve the general covariance of the action. The lowest order dynamical component of the VED is just ∼ ν H 2 , where the dimensionless coefficient ν is naturally predicted to be small (|ν| ≪ 1), but must ultimately be determined experimentally by confronting the model to the cosmological data. That term is nevertheless sufficient to describe the dynamics of the vacuum in the current universe, while the higher order components can play a role in the early universe, and in particular for describing inflation. In fact, in previous works the model has been phenomenologically fitted to a large wealth of cosmological data and the running parameter ν has been found to be positive and in the ballpark of ∼ 10 −3 [13][14][15][16]. Let us finally mention that even though our QFT calculation has been simplified by the use of a single (real) quantum scalar field, further investigations will be needed to generalize these results for multiple fields, involving scalar as well as vector and fermionic components. Up to computational details, however, we expect that the main results of the renormalization program presented here should be maintained.

A Conventions and geometrical quantities
We use natural units, therefore = c = 1 and G N = 1/M P , where M P is the Planck mass. As for the conventions on geometrical quantities used throughout this work, they read as follows: signature of the metric g µν , (−, +, +, +); Riemann tensor, R λ µνσ = ∂ ν Γ λ µσ + Γ ρ µσ Γ λ ρν − (ν ↔ σ); Ricci tensor, R µν = R λ µλν ; and Ricci scalar, R = g µν R µν . Overall, these correspond to the (+, +, +) conventions in the classification by Misner-Thorn-Wheeler [63]. As usual, the Einstein tensor is defined through G µν = R µν − 1 2 Rg µν and the Einstein field equations read G µν + Λg µν = 8πG N T µν . The Christoffel symbols associated to the conformally flat metric ds 2 = a 2 (η)η µν dx µ dx ν , with η µν = diag(−1, +1, +1, +1), are the following: Recalling that the relation between the Hubble rate in conformal and cosmic times is H = aH, the Ricci scalar and the nonvanishing components of the curvature tensors are alternatively given by and We remind the reader that primes indicate differentiation with respect to conformal time and dots differentiation with respect to cosmic time. We also need the higher order curvature tensor (of adiabatic order 4) obtained by functionally differentiating the R 2 -term in the higher derivative vacuum action: Its 00-component in the conformally flat metric reads

B Combining adiabatic and dimensional regularization
In this appendix, we sketch the calculation of the regularized EMT by using dimensional regularization (DR). Let us nonetheless emphasize that while we will use minimal subtraction of poles as a regularization procedure, we do not intend to renormalize the theory with this prescription. If we would do that the renormalized vacuum energy would still exhibit the unwanted ∼ m 4 contributions. In the following, we show that after the ARP has been performed, the divergent integrals appearing in the intermediate calculations can be regularized through DR and then we can recover exactly the same result (6.9) for the renormalized VED.

B.1 Useful formulas
For our purposes it will suffice to focus on integrals of the form where k ≡ |k| and Q is an arbitrary energy scale. If we generalize it to N dimensions, where in the last step we have set N = 3 − 2ǫ. Notice that the scale µ has been introduced for dimensional purposes only and has no obvious physical meaning. It could have equally well been inserted in the original integral in the form d N k → µ 2ǫ d N k to insure that the integration measure has the same dimension as d 3 k. Poles appear for n ≤ 3 in the integral (B.1) and then we can use relations such as which parameterize the divergent behavior of Euler's Γ function near the origin and near −1, respectively, where γ E is Euler's constant. Similar expressions can be generated to parameterize the divergent behavior of Γ around other negative integers using the well-known functional relation Γ(x + 1) = x Γ(x).

B.2 Dimensionally regularized ZPE in FLRW spacetime
Next we summarize how to obtain the same expression for the renormalized VED as the one we have found in Sec. where the divergent and non-divergent contributions are the same ones as in equations (5.2) and (5.3), respectively. The order of adiabaticity of these expressions, therefore is the same as in the calculation presented in the main text, and we shall take this fact for granted hereafter. We should where in the second equality we have used once more ∆ 2 = m 2 − M 2 . At this stage, the DR procedure carries a dependence on the artificial mass scale µ. However, in our case µ will play no role since we are not just aiming at a conventional renormalization based on minimal subtraction, so µ serves only as an auxiliary variable which will eventually disappear from the renormalized result. We should emphasize that the relevant renormalization scale in our calculation is not µ but M . Dimensional regularization is used here only as a technique to display explicitly the divergences of the EMT and to enable their subtraction with the conventional counterterm procedure.

B.3 Counterterms
While the entire calculation can be carried out without any use of DR, provided one defines a properly subtracted EMT from the beginning with the ARP (cf. Sec. 5), we follow now the more conventional approach. Thus we remove the unphysical divergences of the EMT by generating counterterms from the coupling constants present in the extended gravitational action with the HD terms. The modified Einstein's equations read formally as in Eq. (6.1) but carrying the bare couplings, i.e. couplings which are formally UV-divergent and scale independent: We will focus on the 00-component of this equation since we are interested in the ZPE. Following the standard renormalization procedure, we split each of the bare couplings on the l.h.s of the above equation into the renormalized term (which depends on the renormalization point M ), and a counterterm (which does not depend on M ): We define the counterterms such that we can subtract the universal terms γ E and 4π of the DR procedure alongside with the poles, as it is conventional in the modified MS (or MS) [52,53]. That is why we have defined the quantity D ǫ in Eq. (B.6). As we can see, three 'primitive divergences' appear in the unrenormalized form of the EMT, which are proportional to ∼ m 4 , ∼ m 2 (ξ − 1/6) and (ξ − 1/6) 2 , respectively. These can be cancelled by the corresponding counterterms generated from the bare couplings in Eq. (B.11), i.e. the counterterms can now be precisely used to cancel the three divergent quantities proportional to D ǫ in Eq. (B.7). Using the 00-components of the geometric tensors given in Appendix A, they are readily found to be is the finite part left of the EMT (B.4) after removing the poles. Being finite we might be tempted to call it (provisionally) the renormalized ZPE, but in fact is not our final renormalized expression. Some further insight on it can be achieved by considering the term labelled FR, which is given by (B.8). If we apply the limit a = 1 (so that H and all its derivatives vanish) and project the result on-shell (M = m, hence ∆ = 0), the whole expression (B.8) or (B.9) shrinks to just one of the equivalent forms (B.15) This is nothing but the standard (one-loop) ZPE in flat spacetime, namely it is the renormalized form of the UV-divergent integral (4.12) within the MS. As we can see, Eq. (B.15) brings a explicit dependence on µ and above all it grows as the quartic power of the mass of the field. Because the total VED is the sum of (B.15) plus the renormalized ρ Λ -cf. Eq. (3.2) -we are led to face a huge contribution from the quartic term ∼ m 4 (for virtually every known particle, except a very light neutrino), which amounts to a large fine-tuning between these two quantities. This is odd, in fact unacceptable. As discussed in detail in [33], the flat space formula carries indeed the core of the cosmological constant problem [4] and the curved spacetime calculation just inherits it at this point, but it does not aggravate it further. Thus, not surprisingly the subtraction of this part leaves a well-behaved result (cf. Secc. 6.2). However, let us continue with our renormalization procedure and evade this conundrum within the present context.

B.4 Renormalized ZPE and absence of ∼ m 4 contributions
The problem stems from the tilded definition of the renormalized EMT given in (B.14), which is just a variant of the MS-renormalized one, although carrying off-shell ∆ 2 -corrections. However, a well-defined expression can be obtained if we call back anew our definition of renormalized EMT as in Eq. (5.5) of the main text. The prescription amounts to take the on-shell value (at the physical mass m) and subtract from it the terms up to 4th adiabatic order at some arbitrary mass scale M . This provides automatically an overall finite result, as we have proven in the main text without using DR. Taking into account that in this alternative procedure we have already removed the poles appearing in the intermediate steps with the help of DR, it suffices to perform the aforementioned subtraction directly with the finite expression (B.14): The µ-dependence has cancelled at this point, and as we can see this equation turns out to be exactly the same one as in Eq.(5.9). Therefore, from this point onwards we can reproduce the same renormalized VED (6.9), just starting from (6.7) and subtracting its value at the two scales M and M 0 . Once more the result is that the VED at the scale M can be related with its value at another scale M 0 without receiving any contribution from the quartic values of the mass scales or of the mass of the particle. Thus, on using this renormalization procedure we can get rid of the dependence on the quartic powers of the masses as well as on the spurious DR parameter µ. The lesson we can learn is the following. While the mere MS renormalization of the VED (based on using DR together with the subtraction of the poles by the counterterms) leaves a result which is explicitly dependent both on the artifical DR scale µ and on the quartic powers of the masses [38], the extended ARP technique [40] allows to relate the renormalized quantities at two different scales. With detailed calculations, which we have presented here through two different approaches (one of them not using DR at all), we have shown that we can avert the mentioned problems associated to a mere removal of the poles by the counterterms. The common final result is the running of the renormalized EMT in a FLRW background as a function of the Hubble rate, thus allowing to trace the VED evolution throughout the cosmic history. The result we have obtained is indeed much closer in spirit to the renormalization group approach of the RVM, cf. [33,34] and references therein -particularly [30][31][32] -in which such mild evolution of the vacuum energy density in terms of (even) powers of the Hubble rate was predicted on very general grounds. Here, for the first time, we have provided a detailed account from explicit QFT calculations under an appropriate renormalization scheme leading to a possible physical interpretation of the results. The outcome is Eq. (6.9).
We should perhaps repeat once more that such relation is not a prediction of the value of the CC and in general of the VED, as this is out of the scope of renormalization theory. Every renormalization calculation needs a set of renormalization conditions. Behind these renormalization conditions there is a set of physical (and hopefully known) inputs and from these observational inputs we can predict other physical results. In the present instance, this means that given the VED in one scale (entailing a physical input) we can predict its value at another scale. What is, however, distinctive in the kind of calculation we have presented here is the fact that the connection between the renormalized values of the VED at different points appears smooth enough, i.e. it does not involve ∼ m 4 terms, which are usually very large for ordinary particle masses in the standard model of particle physics (let alone in GUT's) and this suggests that no fine tuning is actually involved.
C Identification of the vacuum energy density.
In Minkowskian spacetime we know that the EMT of vacuum takes on the form T vac µν = −ρ vac η µν , being the Lorentz metric η µν the only geometric structure available in a flat background to construct a Lorentz-invariant quantity which can characterize the vacuum state. This allows us to identify the vacuum energy density (VED) from the general structure of the vacuum EMT. It is natural to generalize this identification by assuming that in curved spacetime the vacuum EMT should take the form T vac µν = −ρ vac g µν , with g µν the general background metric. We can formally motivate this result by starting from the EH action with a cosmological term, Varying the part involving the vacuum energy density (i.e. the second term on the r.h.s, which we call S Λ ) yields which provides the identification T vac µν = −ρ vac g µν . This is the line of approach that we have followed here. However, we should point out that such identification has some ambiguities, which as we shall argue below should not alter in a significant way the results that we have obtained and, remarkably enough, lead to a generalized form of the RVM which had actually been considered previously in the literature in different phenomenological formulations [13,14]. In this sense we believe this point deserves being mentioned here, see also [37,41,42].

C.1 More geometric structures for vacuum in curved spacetime
As we know, the vacuum effective action of QFT in the presence of gravity contains the higher derivative terms [26,27]. This is because in curved spacetime we have more geometric quantities to characterize the vacuum, and these structures are actually necessary for implementing the renormalizability of the semiclassical theory of quantized matter fields in an external gravitational background, as we have just seen in our discussion in Appendix B. By the same token one might expect a more general relation between the VED and the EMT, which we may schematize as follows: T vac µν = −ρ vac g µν − α 1 Rg µν − α 2 R µν + O(R 2 ) , (C. 3) in which O(R 2 ) stand for the higher derivative terms, and α i are parameters of dimension +2 in natural units. For g µν = η µν the previous ansatz just boils down to the flat spacetime form mentioned above. To illustrate the possible impact of the additional terms, let us first focus on the following specific form for the renormalized energy-momentum tensor: T vac µν = −ρ vac g µν + λ 16π 2 M 2 G µν + O(R 2 ) . (C.4) Here, G µν is Einstein's tensor (cf. Appendix A). We restrict hereafter our considerations to the late universe since we wish to assess what is the possible impact of the new terms on the parameters that can be directly fitted to observations. The appearance of the mass scale M is necessary for dimensional reasons. Furthermore, λ is an appropriately normalized (dimensionless) parameter. On equating this expression to Eq. (3.2) and considering the 00-component, we find a generalized form of (6.7): We may neglect as before the log evolution ofν eff (H) and approximateν eff (H 0 ) ≃ν eff . The coefficientν is dealt with as constant here, but in general the dimensionless couplings λ i may also be dependent on the scale M , although we expect that the renormalization effects should be logarithmic; and, therefore, following the same practice as withν eff (H), we have not considered these subleadig effects here for the sake of a simpler presentation. We can easily verify that for λ 1 = λ/2 and λ 2 = −λ the last two formulas reduce to (C.6) andν = 0, respectively, since in that case (C.3) boils down to (C.4). Finally, determining ρ vac (M X ) from the boundary condition ρ vac (H 0 ) = ρ 0 vac , we can write down (C.7) in a manifestly normalized way with respect to the current values: 10) in which H 0 andḢ 0 stand, of course, for the respective values of H andḢ at present. The above formula generalizes Eq. (7.6) by including the additional coefficientν, which accompanies the new dynamical term ∼Ḣ. We note that several generalized forms of the RVM containing dynamical terms beyond the canonical one H 2 were studied under different phenomenological conditions, and fitted as well to the data, in [13,14]. Here we have shown that these extended forms of the RVM, which were confronted to observations in the aforementioned references, can also appear as a result of the QFT calculations in the FLRW background.